text:
function [yIntervalEst,PInvIntervalEst,yFwdPred,PInvFwdPred,yFwdEnd,PInvFwdEnd]=FPInfoIntervalSmoother(yFwdPred,PInvFwdPred,yFwdPrev,PInvFwdPrev,N,z,R,H,F,Q,u,hasLastPred)
%%FPINFOINTERVALSMOOTHER Update the Fraser-Potter information smoother
%                smoother over a sliding or growing interval given a new
%                measurement. The smoothed estimates over the entire
%                interval are available. If one already computed the
%                predicted state at the current time, it can be provided
%                rather than being recomputed. This can slide an interval
%                of estimates fixed at length N or grow an interval until
%                it is length N. See the example below for how one can grow
%                the interval until it is a desired length.
%
%INPUTS: yFwdPred The xDimXN set of predicted information states. If the
%                 interval is is increasing in size (has not yet reached
%                 the desired N), then the last dimension is NCur-1 in
%                 size, where NCur is the current length of z, unless
%                 hasLastPred=true, in which case it is NCur in length
%                 with the final prediction already computed. If
%                 hasLastPred=false, then the final column of xFwdPred is
%                 the state predicted to the time-step PRIOR to the current
%                 one. Otherwise, if hasLastPred=true, then the final
%                 column is the state predicted to the time of the latest
%                 measurement in z. This cannot be an empty matrix.
%     PInvFwdPred The xDimXxDimXN set of predicted inverse covariance
%                 matrices associated with yFwdPred. The size and meaning
%                 of the final dimension is the same as for yFwdPred.
%        yFwdPrev The xDimX1 posterior information state estimate at the
%                 time of the second to latest measurement in z. When
%                 calling this function repeatedly, this can be the yFwdEnd
%                 output from the last function call. If hasLastPred=true,
%                 then an empty matrix can be passed in place of this
%                 input.
%     PInvFwdPrev The xDimXxDim inverse covariance matrix associate with
%                 yFwdPrev. When calling this function repeatedly, this can
%                 be the PInvFwdEnd output from the last function call. If
%                 hasLastPred=true, then an empty matrix can be passed in
%                 place of this input.
%               N The positive scalar integer desired length of the
%                 smoothed interval; N>=2. If omitted or an empty matrix is
%                 passed, then N is set to the number of columns of z
%                 (NCur).
%               z The zDimXN set of measurements over the interval. z(:,N)
%                 is the latest measurement. If the interval is expanding,
%                 then z is zDimXNCur in length.
%               R The zDimXzDimXN set of measurement covariance matrices
%                 associated with the values in z. If the interval is
%                 expanding, then the last dimension is NCur in size. If
%                 all the matrices are the same, then a single zDimXzDim
%                 matrix can be passed.
%               H The xDimXxDimXN set of measurement matrices associated
%                 with the values in z. If the interval is expanding, then
%                 the last dimension is NCur in size. If all the matrices
%                 are the same, then a single xDimXzDim matrix can be
%                 passed.
%               F The xDimXxDimX(N-1) set of state transition matrices over
%                 the interval. F(:,:,N-1) is the state transition matrix
%                 to the time of the latest measurerment in z. If the batch
%                 is expanding, then the last dimension is NCur-1 in size.
%                 If all F are the same, then a single xDimXxDim matrix can
%                 be passed.
%               Q The xDimXxDimX(N-1) set of process noise covariance
%                 matrices. the meaning of the third dimension is the same
%                 as for F. If all Q are the same, then a single xDimXxDim
%                 matrix can be passed.
%               u The xDimX(N-1) set of state transition control inputs,
%                 where the second dimension has the same meaning as the
%                 third dimensions of F. . If an empty matrix is passed,
%                 then it is assumed that all control inputs are zero. If
%                 aall the u vectors are the same, then a single xDimX1
%                 vector can be passed.
%     hasLastPred This boolean value indicates whether yFwdPred and
%                 PInvFwdPred contain the prediction to the time of the
%                 latest measurement in z. The default if omitted or an
%                 empty matrix is passed is false.
%
%OUTPUTS: yIntervalEst The xDimXN set of smoothed information state
%                      estimates over the interval. The oldest estimate is
%                      yIntervalEst(:,1)
%                      and the newest yIntervalEst(:,N). If the interval is
%                      expanding, the second dimension will be size NCur.
%      PInvIntervalEst The xDimXxDimXN set of inverse covariance matrices
%                      associated with those in yIntervalEst. If the
%                      interval is expanding, the third dimension will be
%                      size NCur.
% yFwdPred,PInvFwdPred The xDimXN set of predicted information states and
%                      the associated xDimXxDimX(N-1) set of predicted
%                      inverse covariance matrices. If the interval is
%                      expanding, the final dimensions will be length
%                      NCur.
%   yFwdEnd,PInvFwdEnd The xDimX1 and xDimXxDim posterior information state
%                      and its inverse covariance matrix at the latest
%                      time.
%
%The Kalman smoothing technique utilizing forward and backwards filters is
%described in [1]. The use with information filters is related.
%
%EXAMPLE:
%Given an uninformative prior, the information smoother is initialized with
%the output of an information filter at the first step. Then, as more
%measurements are obtained, the batch length increases until it reached
%length N. After that, the interval slides. We also run an information
%filter and show that not only does the smoothed estimate at the end of the
%interval (the most smoothed) outperform the information filter, but that
%the estimate at the start of the interval is the same as the real-time
%information filter output. The measurement error reduction factor (MERF)
%and the NEES of the estimates (converted to target states) of the
%information filter and the information  smoother are compared. The MERF of
%the smoother is lower (better) and the NEES of both are close to 1,
%indicating estimator consistency.
% T=1;%Sample period.
% numRuns=500;
% N=4;%The length of the interval considered.
% numSteps=100;
%
% %Parameters for the noise process dynamics and measurement.
% H=[eye(2,2),zeros(2,2)];
% zDim=size(H,1);
% xDim=size(H,2);
% R=diag([40;40]);
% SR=chol(R,'lower');
%
% %Statistics for the initial state.
% x0Mean=[1e3;0;75;-50];
% x0Cov=diag([1e3^2;100^2;25^2;25^2]);
% x0S=chol(x0Cov,'lower');
%
% %Parameters for the state dynamics.
% q=processNoiseSuggest('PolyKal-ROT',9.8,1);
% F=FPolyKal(T,xDim,1);
% Q=QPolyKal(T,xDim,1,q);
% SQ=chol(Q,'lower');
%
% %The initial uninformative state.
% yInit=zeros(xDim,1);
% PInvInit=zeros(xDim,xDim,1);
%
% xTrue=zeros(xDim,numSteps+N-1,numRuns);
% z=zeros(zDim,numSteps+N-1,numRuns);
% ySmoothed=zeros(xDim,numSteps,numRuns);
% PInvSmoothed=zeros(xDim,xDim,numSteps,numRuns);
% yInfo=zeros(xDim,numSteps,numRuns);
% PInvInfo=zeros(xDim,xDim,numSteps,numRuns);
% yInfoAlt=zeros(xDim,numSteps,numRuns);
% PInvInfoAlt=zeros(xDim,xDim,numSteps,numRuns);
% for curRun=1:numRuns
%     %Draw the initial state.
%     xTrue(:,1,curRun)=x0Mean+x0S*randn(xDim,1);
%     %Get the first measurement.
%     z(:,1,curRun)=H*xTrue(:,1,curRun)+SR*randn(zDim,1);
%
%     %Get the states and measurements for the rest of the steps plus extra
%     %steps so that the Kalman filter and smoother cover the same range.
%     for curStep=2:(numSteps+N-1)
%         xTrue(:,curStep,curRun)=F*xTrue(:,curStep-1,curRun)+SQ*randn(xDim,1);
%         z(:,curStep,curRun)=H*xTrue(:,curStep,curRun)+SR*randn(zDim,1);
%     end
%
%     %Initialize the information filters and the information smoother with a
%     %single measurement.
%     [yUpdate,PInvUpdate]=infoFilterUpdate(yInit,PInvInit,z(:,1,curRun),R,H);
%     yFwdPred=yInit;
%     PInvFwdPred=PInvInit;
%     yFwdPrev=yUpdate;
%     PInvFwdPrev=PInvUpdate;
%
%     yInfo(:,1,curRun)=yUpdate;
%     PInvInfo(:,:,1,curRun)=PInvUpdate;
%
%     yInfoAlt(:,1,curRun)=yUpdate;
%     PInvInfoAlt(:,:,1,curRun)=PInvUpdate;
%
%     %Continue with the rest of the measurements.
%     for curStep=2:numSteps
%         zSpan=z(:,max(curStep-N+1,1):curStep,curRun);
%
%         [yIntervalEst,PInvIntervalEst,yFwdPred,PInvFwdPred,yFwdPrev,PInvFwdPrev]=FPInfoIntervalSmoother(yFwdPred,PInvFwdPred,yFwdPrev,PInvFwdPrev,N,zSpan,R,H,F,Q);
%         if(curStep>=N)
%             ySmoothed(:,curStep-N+1,curRun)=yIntervalEst(:,1);
%             PInvSmoothed(:,:,curStep-N+1,curRun)=PInvIntervalEst(:,:,1);
%         end
%         yInfoAlt(:,curStep,curRun)=yIntervalEst(:,end);
%         PInvInfoAlt(:,:,curStep,curRun)=PInvIntervalEst(:,:,end);
%
%         [yPred,PInvPred]=infoFilterDiscPred(yInfo(:,curStep-1,curRun),PInvInfo(:,:,curStep-1,curRun),F,Q);
%         [yUpdate,PInvUpdate]=infoFilterUpdate(yPred,PInvPred,z(:,curStep,curRun),R,H);
%         yInfo(:,curStep,curRun)=yUpdate;
%         PInvInfo(:,:,curStep,curRun)=PInvUpdate;
%     end
%
%     %Now, go N-1 steps more with the smoother, so that the information
%     %filter estimates and the information smoother estimates cover the same
%     %range of values.
%     for curStep=(numSteps+1):(numSteps+N-1)
%         zSpan=z(:,(curStep-N+1):curStep,curRun);
%         [yIntervalEst,PInvIntervalEst,yFwdPred,PInvFwdPred,yFwdPrev,PInvFwdPrev]=FPInfoIntervalSmoother(yFwdPred,PInvFwdPred,yFwdPrev,PInvFwdPrev,N,zSpan,R,H,F,Q);
%         ySmoothed(:,curStep-N+1,curRun)=yIntervalEst(:,1);
%         PInvSmoothed(:,:,curStep-N+1,curRun)=PInvIntervalEst(:,:,1);
%     end
% end
%
% %The maximum absolute difference between the standard information filter
% %estimate and the estimate at the front of the infomration smoother batch.
% %This should be 0.
% maxInfoDiff=max(abs(yInfo(:)-yInfoAlt(:)))
%
% %Get the target states. We shall skip the first one which is not completely
% %observable.
% xSmoothed=zeros(xDim,numSteps-1,numRuns);
% PSmoothed=zeros(xDim,xDim,numSteps-1,numRuns);
% xEst=zeros(xDim,numSteps-1,numRuns);
% PEst=zeros(xDim,xDim,numSteps-1,numRuns);
% for curRun=1:numRuns
%     for curStep=2:numSteps
%         xSmoothed(:,curStep-1,curRun)=PInvSmoothed(:,:,curStep,curRun)\ySmoothed(:,curStep,curRun);
%         PSmoothed(:,:,curStep-1,curRun)=inv(PInvSmoothed(:,:,curStep,curRun));
%
%         xEst(:,curStep-1,curRun)=PInvInfo(:,:,curStep,curRun)\yInfo(:,curStep,curRun);
%         PEst(:,:,curStep-1,curRun)=inv(PInvInfo(:,:,curStep,curRun));
%     end
% end
%
% %Limit the period plotted to the range where the target state is fully
% %observable and is common to the smoother and the filter.
% xTrue=xTrue(:,2:numSteps,:);
% z=z(:,2:numSteps,:);
% posTrue=xTrue(1:2,:,:);
%
% MERFSmoothed=calcMERF(posTrue,z,xSmoothed(1:2,:,:),0,true);
% MERFInfo=calcMERF(posTrue,z,xEst(1:2,:,:),0,true);
% NEESSmoothed=calcNEES(xTrue,xSmoothed,PSmoothed);
% NEESInfo=calcNEES(xTrue,xEst,PEst);
%
% figure(1)
% clf
% hold on
% plot(2:numSteps,MERFInfo,'-k','linewidth',4)
% plot(2:numSteps,MERFSmoothed,'--r','linewidth',2)
% h1=xlabel('Step');
% h2=ylabel('MERF');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
% legend('Information Filter','Information Smoother')
%
% figure(2)
% clf
% hold on
% plot(2:numSteps,NEESInfo,'-k','linewidth',4)
% plot(2:numSteps,NEESSmoothed,'--r','linewidth',2)
% h1=xlabel('Step');
% h2=ylabel('NEES');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
% legend('Information Filter','Information Smoother','location','southeast')
%
%REFERENCES:
%[1] D. C. Fraser and J. E. Potter, "The optimal linear smoother as a
%    combination of two optimum linear filters," IEEE Transactions on
%    Automatic Control, pp. 387-390, Aug. 1969.
%
%August 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(yFwdPred,1);
NCur=size(z,2);

if(nargin<12||isempty(hasLastPred))
    hasLastPred=false;
end

if(isempty(N))
    %Assume that the batch is already the desired size if the batch length
    %is not specified.
    N=NCur;
end

NFEnd=NCur-1;

if(size(F,3)==1)
    F=repmat(F,[1,1,NFEnd]);
end

if(size(Q,3)==1)
    Q=repmat(Q,[1,1,NFEnd]);
end

if(nargin<11||isempty(u))
    u=zeros(xDim,NFEnd);
elseif(size(u,2)==1)
    u=repmat(u,[1,NFEnd]);
end

if(size(H,3)==1)
    H=repmat(H,[1,1,NCur]);
end

if(size(R,3)==1)
    R=repmat(R,[1,1,NCur]);
end

if(hasLastPred==false)
    [yFwdPredEnd,PInvFwdPredEnd]=infoFilterDiscPred(yFwdPrev,PInvFwdPrev,F(:,:,NFEnd),Q(:,:,NFEnd),u(:,NFEnd));

    if(NCur==N&&size(yFwdPred,2)==N)
        %If the interval is sliding.
        yFwdPred=[yFwdPred(:,2:N),yFwdPredEnd];
        PInvFwdPred=cat(3,PInvFwdPred(:,:,2:N),PInvFwdPredEnd);
    else%The interval is growing; NCur is one more than the current length
        %of yFwdPred.
        yFwdPred(:,NCur)=yFwdPredEnd;
        PInvFwdPred(:,:,NCur)=PInvFwdPredEnd;
    end
end

%Run the information filter backwards.
yRev=zeros(xDim,NCur);
PInvRev=zeros(xDim,xDim,NCur);

[yRev(:,NCur),PInvRev(:,:,NCur)]=infoFilterUpdate(zeros(xDim,1),zeros(xDim,xDim),z(:,NCur),R(:,:,NCur),H(:,:,NCur));
for curStep=(NCur-1):-1:1
    [yRevPred,PInvRevPred]=infoFilterDiscPredRev(yRev(:,curStep+1),PInvRev(:,:,curStep+1),F(:,:,curStep),Q(:,:,curStep),u(:,curStep));
    [yRev(:,curStep),PInvRev(:,:,curStep)]=infoFilterUpdate(yRevPred,PInvRevPred,z(:,curStep),R(:,:,curStep),H(:,:,curStep));
end

PInvIntervalEst=PInvFwdPred+PInvRev;
yIntervalEst=yFwdPred+yRev;

if(nargout>4)
    yFwdEnd=yIntervalEst(:,end);
    PInvFwdEnd=PInvIntervalEst(:,:,end);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function varargout = spm_unlink(varargin)
% routine for silently deleting files on disk - a compiled routine
% FORMAT spm_unlink file1 file2 file3 file4...
%     OR spm_unlink('file1','file2','file3','file4',...)
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_unlink.m 1143 2008-02-07 19:33:33Z spm $

%-This is merely the help file for the compiled routine
error('spm_unlink.c not compiled - see Makefile')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function c=xxxc(a,b)
n=length(a);
c1=0;
c2=0;
c3=0;
for i=1:n
    c1=c1+a(i)*b(i);
    c2=c2+a(i)^2;
    c3=c3+b(i)^2;
end
c=c1/(sqrt(c2)*sqrt(c3));

repo: DeepFi


text:
%computes the correlation (matrix) sequence given a stable ss description
function [corrmat]=ss2corr_v000(A, Q, C, P0, dt, crind)

nout=size(C,1);
nout=nout*(nout+1)/2;
nst=size(A,1);
corrmat=zeros(nout,length(crind));

%compute the initial covariance
if (isempty(P0))
    %check whether system is stable or not
    eigenval=eig(A);
    if (~isempty(find(eigenval>=1,1)))
        disp('A is not stable');
        eigenval
        return;
    else
        %Use lyapunov eq.
        P0=dlyap(A,Q);
    end
end

%convert the discrete time into cont
crcind=crind*dt;
[Ac, Qc]=dc2dc_v000(A, Q, [], dt, 1,[]);

%compute the covariances
for in=1:length(crind)
    sr_a=crcind(in);
    Ad=expm(Ac*abs(sr_a));
    %[Ad, Qd]=dc2dc_v000(Ac, Qc, [], abs(crcind(in)), 0,[]);
    if (sr_a>0)
        mx_a=C*(Ad*P0)*C';
    elseif (sr_a<0)
        mx_a=C*(P0*Ad')*C';
    else
        mx_a=C*P0*C';
    end
    corrmat(:,in)=mat2vec(mx_a);
end

function vec=mat2vec(mat)
nrow=size(mat,1);
nout=nrow*(nrow+1)/2;
vec=zeros(nout,1);
lb=1;
ub=nrow;
for in=1:size(mat,1)
    vec(lb:ub)=diag(mat,in-1);
    lb=ub+1;
    ub=ub+nrow-in;
end

repo: ekfmonoslam


text:
function X = negatefactors(X,L)
if ~isempty(X.midfactors)
    for i = 1:length(X.midfactors)
        X.leftfactors{i} = -X.leftfactors{i};
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function results  = RS_DF_rstest(z);

rvec = (0:0.001:1);
P = size(z,1);

cumcumz = [];
for r = rvec
    cumcumz = [cumcumz (z<r)-r];
end

v = sum(cumcumz,1)/sqrt(P);

Qv = []; Qvabs = [];
for rs = 1:size(rvec,2);
    Qv = [Qv; v(1,rs)'*v(1,rs)];
    Qvabs = [Qvabs; abs(v(1,rs)')];
end

Kv = max(Qvabs);
CVMv = mean(Qv);

results(1,1) = Kv;
results(1,2) = CVMv;

author: european-central-bank

repo: BEAR-toolbox


text:
function [FutureCF AccruedInt] = czbondfuturecf(CouponRate, Issue, Maturity, Settle, Basis, Notional)
% =========================================================================
% CZBONDFUTURECF generates future (and current - settlement date) cash flows,
% cash flow dates and cash flow time fraction according to day count basis (30E/360, Actual/360)
% for Czech Government Bonds. Further, it calculates accrued interest.
%
% [Yield AccruedInt] = czgbfuturecf(CouponRate, Issue, Maturity,...
%       Settle, Basis, Notional)
%
% INPUTS:
%         CouponRate - Coupon rate in decimal form
%              Issue - Issue Date in serial date number
%                      (NOTE: convert date string using DATENUM)
%           Maturity - Maturity Date in serial date number
%             Settle - Settlment Date in serial date number
%
% OPTIONAL INPUT:
%           Notional - Notional of the bond (default is 100)
%              Basis - Day-count basis
%                      1 - 30E/360 (default)
%                      2 - Actual/360
%
% OUTPUT:
%            FutureCF - Future cash flow matrix
%                       - 1st row, nominal cash flows
%                       - 2nd row, cash flow dates
%                       - 3rd row, cash flow time fraction measured in years
%
%         AccruedInt - Accrued interest since last coupon
%
% NOTE:
%         It doesn't work take ex-coupon date into account
%         Normally, Czech Government Bonds are traded according to 30E/360 convention!
%
% USES: czbondcf
%
% Kamil Kladivko
% email: kladivko@gmail.com
% December 2010
% Cite as:
% Kladivko Kamil (2010). The Czech Treasury Yield Curve from 1999 to the Present,
% Czech Journal of Economics and Finance, 60(4): 307-335
% =========================================================================

    if nargin < 4, error('Need at least 4 inputs'); end
    if nargin == 4, Basis = 0; Notional = 100; end
    if nargin == 5, Notional = 100; end
    if all(Basis ~= [1, 2])
        warning('Unknown code for day count convention. Using 30E/360');
        Basis = 1;
    end
    % Generates CZGB cash flows
    CF = czbondcf(CouponRate, Issue, Maturity, Basis, Notional);
    TempCF = CF(:, 2:end); % Don't want first negative CF (notional of the bond)
    % Pick only furture and current (Settlment date) cash flows
    [SettYr, SettMo, SettDay] = datevec(Settle);
    FutureCF = TempCF(:, TempCF(2, :) >= datenum(Settle));
    % Find last cash flow date before settlment
    [junk FutureCFcount] = size(FutureCF);
    AccrualDate = CF(2, end-FutureCFcount);
    % Acrrued Interest
    switch Basis
        case 1 % 30E/360
            [AccrualYr, AccrualMo, AccrualDay] = datevec(AccrualDate);
            if SettDay == 31, SettDay = 30; end; if AccrualDay == 31, AccrualDay = 30; end
            DaysCount = 360*(SettYr - AccrualYr) + 30*(SettMo - AccrualMo) + (SettDay - AccrualDay);
            AccruedInt = CouponRate*Notional*DaysCount/360;
        case 2 % Actual/360;
            DaysCount = datenum(Settle) - AccrualDate;
            AccruedInt = CouponRate*Notional*DaysCount/360;
    end
    % Time Factor
    switch Basis
        case 1 % 30E/360
            DaysCount = zeros(1, FutureCFcount);
            for i = 1:FutureCFcount
                [NxtYr, NxtMo, NxtDay] = datevec(FutureCF(2,i));
                if NxtDay == 31, NxtDay = 30; end
                DaysCount(i) = 360*(NxtYr - SettYr) + 30*(NxtMo - SettMo) + (NxtDay - SettDay);
            end
            TimeFactor = DaysCount./360;
        case 2 % Actual/360;
              DaysCount =  FutureCF(2,:) - datenum(Settle);
              TimeFactor = DaysCount./360;
    end
    FutureCF = [FutureCF; TimeFactor];
end

author: Sable

repo: mcbench-benchmarks


text:
function z = colon( x, y )

if ~isa( x, 'cvxdual' ),
    z = cvx( x );
    d = y;
elseif ~isa( y, 'cvxdual' ),
    z = cvx( y );
    d = x;
else
    error( 'Usage: <dual var> : <expression> or <expression> : <dual var>' );
end

if inuse( d ),
    nm = cvx_subs2str( name( d ) );
    error( [ 'Dual variable "', nm(2:end), '" has already been used.' ] );
else
    z = setdual( z, name( d ) );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function res = hlp_aggregatestructs(structs,defaultop,varargin)
% Aggregate structs (recursively), using the given combiner operations.
% Result = hlp_aggregatestructs(Structs,Default-Op,Field-Ops...)
%
% This results in a single 1x1 struct which has aggregated values in its fields (e.g., arrays,
% averages, etc.). For a different use case, see hlp_superimposedata.
%
% In:
%   Structs    : cell array of structs to be aggregated (recursively) into a single struct
%
%   Default-Op : optional default combiner operation to execute for every field that is not itself a
%                struct; see notes for the format.
%
%   Field-Ops  : name-value pairs of field-specific ops; names can have dots to denote operations
%                that apply to subfields. field-specific ops that apply to fields that are
%                themselves structures become the default op for that sub-structure
%
% Out:
%   recursively merged structure.
%
% Notes:
%   If an operation cannot be applied, a sequence of fall-backs is silently applied. First,
%   concatenation is tried, then, replacement is tried (which never fails). Therefore,
%   function_handles are being concatenated up to 2008a, and replaced starting with 2008b.
%   Operations are specified in one of the following formats:
%   * 'cat': concatenate values horizontally using []
%   * 'replace': replace values by those of later structs (noncommutative)
%   * 'sum': sum up values
%   * 'mean': compute the mean value
%   * 'std': compute the standard deviation
%   * 'median': compute the median value
%   * 'random': pick a random value
%   * 'fillblanks': replace [] by values of later structs
%   * binary function: apply the function to aggregate pairs of values; applied in this order
%     f(f(f(first,second),third),fourth)...
%   * cell array of binary and unary function: apply the binary function to aggregate pairs of
%     values, then apply the unary function to finalize the result: functions {b,u} are applied in
%     the following order: u(b(b(b(first,second),third),fourth))
%
% Examples:
%   % calc the average of the respective field values, across structs
%   hlp_aggregatestructs({result1,result2,result3},'mean')
%
%   % calc the std deviation of the respective field values, across structs
%   hlp_aggregatestructs({result1,result2,result3},'std')
%
%   % concatenate the field values across structs
%   hlp_aggregatestructs({result1,result2,result3},'cat')
%
%   % as before, but use different operations for a few fields
%   hlp_aggregatestructs({result1,result2,result3},'cat','myfield1','mean','myfield2.subfield','median')
%
%   % use a custom combiner operation (here: product)
%   hlp_aggregatestructs({result1,result2,result3},@times)
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-05-04

% Copyright (C) Christian Kothe, SCCN, 2010, christian@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
% even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

warning off MATLAB:warn_r14_function_handle_transition

if nargin < 2
    defaultop = 'cat'; end
if ~iscell(structs)
    structs = {structs}; end
fieldops = hlp_varargin2struct(varargin);

% translate all ops (if they are specified as strings)
defaultop = translateop(defaultop);
fieldops = translateop(fieldops);

% aggregate, then finalize
res = finalize(aggregate(structs,defaultop,fieldops),defaultop,fieldops);

function res = aggregate(structs,defaultop,fieldops)
% we skip empty records
structs = structs(~cellfun('isempty',structs));

% for any struct array in the input, we first merge it recursively as if it were a cell array
for k=find(cellfun('length',structs)>1)
    structs{k} = hlp_aggregatestructs(structarray2cellarray(structs{k}),defaultop,fieldops); end

% at this point, we should have a cell array of structs
if ~isempty(structs)
    % we begin with the first struct
    res = structs{1};
    % and aggregate the remaining ones onto it
    for i=2:length(structs)
        si = structs{i};
        % proceeding field by field...
        for fn=fieldnames(si)'
            f = fn{1};
            % figure out which operation applies
            if isfield(fieldops,f)
                % a field-specific op applies
                if isstruct(fieldops.(f))
                    % ... which is itself a struct
                    fop = fieldops.(f);
                else
                    op = fieldops.(f);
                end
            else
                % the default op applies
                op = defaultop;
                fop = fieldops;
            end

            % now process the field
            if ~isfield(res,f)
                % field is not yet in the aggregate: just assign
                res.(f) = si.(f);
            else
                % need to aggregate it
                if isstruct(res.(f)) && isstruct(si.(f))
                    % both are a struct: recursively aggregate
                    res.(f) = aggregate({res.(f),si.(f)},op,fop);
                else
                    % they are not both structus
                    try
                        % try to apply the combiner op
                        res.(f) = op{1}(res.(f),si.(f));
                    catch
                        % didn't work: try to concatenate as fallback
                        try
                            res.(f) = [res.(f),si.(f)];
                        catch
                            % didn't work: try to assign as fallback (dropping previous field)
                            res.(f) = si.(f);
                        end
                    end
                end
            end
        end
    end
else
    % nothing to aggregate
    res = [];
end

function x = finalize(x,defaultop,fieldops)
% proceed field by field...
for fn=fieldnames(x)'
    f = fn{1};
    % figure out which operation applies
    if ~isempty(fieldops) && isfield(fieldops,f)
        % a field-specific op applies
        if isstruct(fieldops.(f))
            % ... which is itself a struct
            fop = fieldops.(f);
        else
            op = fieldops.(f);
        end
    else
        % the default op applies
        op = defaultop;
        fop = fieldops;
    end
    try
        % now apply the finalizer
        if isstruct(x.(f))
            % we have a sub-struct: recurse
            x.(f) = finalize(x.(f),op,fop);
        else
            % we have a regular element: apply finalizer
            x.(f) = op{2}(x.(f));
        end
    catch
        % for empty structs, x.(f) produces no output
    end
end

% translate string ops into actual ops, add the default finalizer if missing
function op = translateop(op)
if isstruct(op)
    % recurse
    op = structfun(@translateop,op,'UniformOutput',false);
else
    % remap strings
    if ischar(op)
        switch op
            case 'cat'
                op = @(a,b)[a b];
            case 'replace'
                op = @(a,b)b;
            case 'sum'
                op = @(a,b)a+b;
            case 'mean'
                op = {@(a,b)[a b], @(x)mean(x)};
            case 'median'
                op = {@(a,b)[a b], @(x)median(x)};
            case 'std'
                op = {@(a,b)[a b], @(x)std(x)};
            case 'random'
                op = {@(a,b)[a b], @(x) x(min(length(x),ceil(eps+rand(1)*length(x))))};
            case 'fillblanks'
                op = @(a,b)quickif(isempty(a),b,a);
            otherwise
                error('Unsupported combiner operation specified: %s',op);
        end
    elseif ~isa(op,'function_handle')
        error('The given combiner operation must be either a string or a function handle, but was: %s',hlp_tostring(op,1000));
    end
    % add finalizer if missing
    if ~iscell(op)
        op = {op,@(x)x}; end
end

% inefficiently turn a struct array into a cell array of structs
function res = structarray2cellarray(arg)
res = {};
for k=1:numel(arg)
    res = [res {arg(k)}]; end

% for the 'fillblanks' translate op
function val = quickif(cond,trueval,falseval)
if cond
    val = trueval;
else
    val = falseval;
end


text:
problemData.problemFileName = 'test2d_micro';
problemData.scale = 'MICRO';

designVarSettings.type        = 'LevelSet';
designVarSettings.initialCase = 'circleInclusion';

incrementalSchemeSettings.nSteps = 1;

costSettings.shapeFuncSettings = {
    struct('type','chomog_alphabeta','alpha',[1 0 0]','beta',[0 -1 0]')
    struct('type','perimeterConstraint','PerimeterTarget',5)
    };
costSettings.weights = [1 0.1];

constraintSettings.shapeFuncSettings = struct('type','volumeConstraint');

optimizerSettings.type       = 'AlternatingPrimalDual';
optimizerSettings.shallPrint = false;
optimizerSettings.settingsMonitor.showOptParams       = true;
optimizerSettings.settingsMonitor.refreshInterval     = 2;
optimizerSettings.settingsMonitor.shallDisplayDesignVar       = true;
optimizerSettings.settingsMonitor.shallShowBoundaryConditions = true;

optimizerSettings.uncOptimizerSettings.type       = 'SLERP';

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% This is a demonstration of deterministic neoclassical growth model in
% macro economic theory programed for phd course I'm taking in CMU
% for full reference of model See  "Recersive methods in Economi Dynamics", Nancy L stockey and robert E. Lucas. Jr.
% you can change model parameter for your need
%
% There is no population growth or technology progress
% Labor supply is inelastic
% Method=1----Deterministic, discrete state space model, No interpolation
% Method=2----Deterministic, continuous state space model, With interpolation
%
% Utility function is util.m
% Production function is cobb.m
% Plotting function vgc_glot.m
%
% Author Dan Li, Carnegie Mellon University
% Email :  danli@andrew.cmu.edu
% Homepage  www.danli.org
% April.2003

clc;
clear all;

% ------Model Parameter specification---------
global gamma alpha
gamma=1;        % risk aversion coefficient, log utility function
alpha=0.3;      % capital share
beta=0.96;      % discount factor
delta=.1;       %   depreciation rate
z=1  ;           % initial state for deteriministic case

method=input(['Method=?,\n'...
             '1 for Deterministic,No interpolation \n' ...
             '2 for Deterministic With Interpolation,\n']);

%------ Parameter of computation--------------------------
nk=input(['Number of grid?(20-200 for Method 1)\n'...
          '(<=30 for Method 2, <=200 for Method 3)>>']);   % number of grid
maxloop=40;         % number of maximum iteration
tol=10^(-3);        % tolerence level for value function iteration

tic; % start the clock

% ******Low and upper bound of state variable k ******
% Steady state value kstar,   f'(kstar)=1/beta-1+delta

kstar=(((1/(alpha*beta))-((1-delta)/alpha)))^(1/(alpha-1));
lowk=0.8*kstar;
upk=1.2*kstar;

k=linspace(lowk,upk,nk)';  % build up grid
[K,Kprime]=meshgrid(k);

v=zeros(size(k))*ones(size(z));     vn=v;
g=v;     gn=g;
loop=0;     dv=1;   dk=1;

if method==1
%--------------------------------------
%  No Interpolation Method, Method=1
%---------------------------------------
	C=cobb(K)+(1-delta)*K-Kprime;
	U=util(C);

	kpos=100*ones(size(k));% initial value of loop and distance of value

	while dv>tol & loop<=maxloop & dk
        loop=loop+1;
        [vn,kposn]=max(U+beta*v*ones(1,nk));  %new value vector
        vn=vn'; kposn=kposn';
        dv=max(abs((vn-v)./(v+eps)));
        dk=~all(kpos==kposn); % whether or not position exactly the same
        v=vn; kpos=kposn;
	end
	g=k(kpos);  %  policy  function of next period capital

elseif method==2;
%--------------------------------------
%  Interpolation Method, Method=2
%---------------------------------------
   	while dv>tol & dk>tol & loop<=maxloop
	loop=loop+1;
	pp = spline(k,v);
	vf=inline('-(util(-y+cobb(k)+(1-delta)*k)+beta*ppval(pp,y))','y','k','delta','beta','pp');
		for ik=1:nk
            [gn(ik),vn(ik)]=fminbnd(vf,lowk,upk,[],k(ik),delta,beta,pp) ;
            vn(ik)=-vn(ik);
		end;
	dk=max(abs((gn-g)./(g+eps)));
    dv=max(abs((vn-v)./(v+eps)));
	v=vn; g=gn;
    end;
end
%------------------------------------------------------------
% use calculated k g to get consumption policy and plot , save
%--------------------------------------------------------------
%Polc=cobb(k)*z+(1-delta)*k*ones(size(z))-g; %  policy function of consumption
Polc=cobb(k)+(1-delta)*k-g; %  policy function of consumption

vgc_plot(k,kstar,v,Polc,g, method);

%--------------so other output for observation of convergence--------------
t=toc  % read clock


text:
%This Matlab script can be used to reproduce Figure 7.2 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.
%
%Note: This script require additional software packages to be used, which
%need to be downloaded and installed separately. These packages are
%developed independently and are delivered with separate licenses.
%
%The downlink power allocation is optimized using CVX from CVX Research,
%Inc. (http://cvxr.com/cvx/). This script has been tested with CVX 2.1,
%using the solver Mosek, version 7.1.0.12. We discourage the use of the
%solvers SDPT3 and SeDuMi since these crashed during the test.

%Empty workspace and close figures
close all;
clear;

%Number of BSs
L = 16;

%Number of UEs per BS
K = 10;

%Number of BS antennas
M = 100;

%Define the pilot reuse factor
f = 2;

%Select the number of setups with random UE locations
nbrOfSetups = 50;

%Select the number of channel realizations per setup
nbrOfRealizations = 500;

%% Propagation parameters

%Communication bandwidth
B = 20e6;

%Total uplink transmit power per UE (mW)
p = 100;

%Total downlink transmit power per UE (mW)
rho = 100;

%Maximum downlink transmit power per BS (mW)
Pmax = K*rho;

%Compute downlink power per UE in case of equal power allocation
rhoEqual = (Pmax/K)*ones(K,L);

%Define noise figure at BS (in dB)
noiseFigure = 7;

%Compute noise power
noiseVariancedBm = -174 + 10*log10(B) + noiseFigure;

%Select length of coherence block
tau_c = 200;

%Use the approximation of the Gaussian local scattering model
accuracy = 2;

%Angular standard deviation in the local scattering model (in degrees)
ASDdeg = 10;

%Prepare to save simulation results
SE_MR_equal = zeros(K,L,nbrOfSetups);
SE_RZF_equal = zeros(K,L,nbrOfSetups);
SE_MMMSE_equal = zeros(K,L,nbrOfSetups);

SE_MR_maxmin = zeros(K,L,nbrOfSetups);
SE_RZF_maxmin = zeros(K,L,nbrOfSetups);
SE_MMMSE_maxmin = zeros(K,L,nbrOfSetups);

SE_MR_maxprod = zeros(K,L,nbrOfSetups);
SE_RZF_maxprod = zeros(K,L,nbrOfSetups);
SE_MMMSE_maxprod = zeros(K,L,nbrOfSetups);

%% Go through all setups
for n = 1:nbrOfSetups

    %Output simulation progress
    disp([num2str(n) ' setups out of ' num2str(nbrOfSetups)]);

    %Compute channel statistics for one setup
    [R,channelGaindB] = functionExampleSetup(L,K,M,accuracy,ASDdeg);

    %Compute the normalized average channel gain, where the normalization
    %is based on the noise power
    channelGainOverNoise = channelGaindB - noiseVariancedBm;

    %Generate channel realizations with estimates and estimation
    %error correlation matrices
    [Hhat,C,tau_p,~,H] = functionChannelEstimates(R,channelGainOverNoise,nbrOfRealizations,M,K,L,p,f);

    %Compute the signal and interference terms of the DL SEs using the
    %hardening bound in Theorem 4.6
    [signal_MR,interf_MR,signal_RZF,interf_RZF,signal_MMMSE,interf_MMMSE,prelogFactor] = functionComputeSINR_DL(H,Hhat,C,tau_c,tau_p,nbrOfRealizations,M,K,L,p);

    %Delete large matrices
    clear Hhat C H R;

    %Compute the SEs with equal power allocation
    SE_MR_equal(:,:,n) = functionComputeSE_DL_poweralloc(rhoEqual,signal_MR,interf_MR,prelogFactor);
    SE_RZF_equal(:,:,n) = functionComputeSE_DL_poweralloc(rhoEqual,signal_RZF,interf_RZF,prelogFactor);
    SE_MMMSE_equal(:,:,n) = functionComputeSE_DL_poweralloc(rhoEqual,signal_MMMSE,interf_MMMSE,prelogFactor);

    %Compute the SEs with max-min fairness power allocation
    disp('Solve max-min fairness problem'); %Output simulation progress
    SE_MR_maxmin(:,:,n) = functionPowerOptimization_maxmin(signal_MR,interf_MR,Pmax,prelogFactor);
    SE_RZF_maxmin(:,:,n) = functionPowerOptimization_maxmin(signal_RZF,interf_RZF,Pmax,prelogFactor);
    SE_MMMSE_maxmin(:,:,n) = functionPowerOptimization_maxmin(signal_MMMSE,interf_MMMSE,Pmax,prelogFactor);

    %Compute the SEs with max-min fairness power allocation
    disp('Solve max product SINR problem'); %Output simulation progress
    SE_MR_maxprod(:,:,n) = functionPowerOptimization_prodSINR(signal_MR,interf_MR,Pmax,prelogFactor);
    SE_RZF_maxprod(:,:,n) = functionPowerOptimization_prodSINR(signal_RZF,interf_RZF,Pmax,prelogFactor);
    SE_MMMSE_maxprod(:,:,n) = functionPowerOptimization_prodSINR(signal_MMMSE,interf_MMMSE,Pmax,prelogFactor);

end

%% Plot the simulation results

figure;
hold on; box on;

plot(sort(SE_MR_maxmin(:)),linspace(0,1,K*L*nbrOfSetups),'k--','LineWidth',1);
plot(sort(SE_MR_equal(:)),linspace(0,1,K*L*nbrOfSetups),'k-','LineWidth',1);
plot(sort(SE_MR_maxprod(:)),linspace(0,1,K*L*nbrOfSetups),'k-.','LineWidth',1);

plot(sort(SE_RZF_equal(:)),linspace(0,1,K*L*nbrOfSetups),'b-','LineWidth',1);
plot(sort(SE_MMMSE_equal(:)),linspace(0,1,K*L*nbrOfSetups),'r','LineWidth',1);

plot(sort(SE_RZF_maxprod(:)),linspace(0,1,K*L*nbrOfSetups),'b-.','LineWidth',1);
plot(sort(SE_MMMSE_maxprod(:)),linspace(0,1,K*L*nbrOfSetups),'r-.','LineWidth',1);

plot(sort(SE_RZF_maxmin(:)),linspace(0,1,K*L*nbrOfSetups),'b--','LineWidth',1);
plot(sort(SE_MMMSE_maxmin(:)),linspace(0,1,K*L*nbrOfSetups),'r--','LineWidth',1);

legend('Max-min fairness','Equal power','Max product SINR','Location','SouthEast');

xlabel('SE per UE [bit/s/Hz]');
ylabel('CDF');

author: emilbjornson

repo: massivemimobook


text:
function [ d_mid ] = searchForDistanceLoss( Loss_wanted, freq, land_sea, urban_suburban_rural, h_a, h_eff, h_2, tPct, lPct, mobile, indoor)
%searchForDistanceLoss finds the distance at which a certain loss value
%occurs given a certain set of parameters for the JTG5-6 model

d_upper = 1000;
d_lower = 0.001;
while 1
    d_mid = d_lower + (d_upper-d_lower)/2;
    Loss_mid = JTG5_6( freq , land_sea, urban_suburban_rural, h_a, h_eff, h_2, d_mid, tPct, lPct, mobile, indoor) ;%-tx_erp ;
    delta = abs(Loss_mid - Loss_wanted );
    if delta < 0.05
        break
    else
        if Loss_mid > Loss_wanted
            d_upper = d_mid;
        else
            d_lower = d_mid;
        end
    end
end
end


text:
function psnr = calc_psnr(V, W, H)

    % PSNR = 10 log10 (MAX^2/MSE)
    %
    %       MAX_VAL: Maximum value of pixels

    max_val = max(max(V));
    mse = calc_mse(V, W, H);
    psnr = 10 * log10 (max_val.^2/mse);

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [vMat,fMat] = spheretribydepth(depth)
% SPHERETRIBYDEPTH is a high-performance vectorized function for building
% a triangulation of a unit sphere based on recursive partitioning of each
% of Icosahedron faces into 4 triangles with vertices in the middles of
% original face
%
% Input:
%   depth: double[1,1] - depth of partitioning, use 1 for the first level of
%       Icosahedron partitioning, and greater value for a greater level
%       of partitioning
%
% Output:
%   vMat: double[nVerts,3] - (x,y,z) coordinates of triangulation
%       vertices
%   fMat: double[nFaces,3] - indices of face verties in vertMat
%
% Example:
%   [vMat, fMat] = spheretribydepth(4);
%   patch('Vertices',vMat,'Faces',fMat,'FaceColor','g','EdgeColor','k');
%
% $Author: Peter Gagarinov, PhD  <pgagarinov@gmail.com> $
% $Copyright: Peter Gagarinov, PhD,
%            Moscow State University,
%            Faculty of Computational Mathematics and Computer Science,
%            System Analysis Department 2011-2016 $
%
if ~(numel(depth)&&isnumeric(depth)&&depth>=0&&fix(depth)==depth)
    error('spheretri:wrongInput',...
        'depth is expected to be a not negative integer scalar');
end
[vMat,fMat]=icosahedron();
[vMat,fMat]=shrinkfacetri(vMat,fMat,0,depth,@normvert);
end
function x=normvert(x)
x=x./repmat(realsqrt(sum((x.*x),2)),1,3);
end

repo: OpenMAS


text:
function ind = squareNeighbors2(RowsCols,Nid)
% return the neighbor ids in a square grid
%
%

[col,row] = meshgrid(1:RowsCols(2),1:RowsCols(1));

dRow = [1 -1 0  0 1  1 -1 -1];
dCol = [0  0 1 -1 1 -1  1 -1];

if nargin == 2
  ind = calcInd(Nid);
else

  ind = zeros([RowsCols,length(dRow)]);
  for Nid = 1:length(dRow)
    ind(:,:,Nid) = calcInd(Nid);
  end

end

  function indLocal = calcInd(Nid)

  nrow = row + dRow(Nid);
  ncol = col + dCol(Nid);

  % ensure coordinates are within the range
  ncol = max(min(ncol,RowsCols(2)),1);
  nrow = max(min(nrow,RowsCols(1)),1);

  indLocal = sub2ind(RowsCols,nrow,ncol);

  end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function bytes = get_file_size(fname)
% gets file size in bytes, ensuring that symlinks are dereferenced
% MP: not sure who wrote this code, but they were careful to do it right on Linux
    bytes = NaN;
    if isunix
        cmd = sprintf('stat -Lc %%s %s', fname);
        [status, r] = system(cmd);
        if status == 0
            bytes = str2double(r);
        end
    end

    if isnan(bytes)
        o = dir(fname);
        bytes = o.bytes;
    end
end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
% Compute QCFC correlations and plot them (Parkes et al 2018, Ciric et al
% 2018, Power et al 2012). Computes edgewise correlations between
% connectivity estimates (stored in bs.connectivity.regions.r) and a QC metric (typically head
% motion) passed in by user
function qcfc(bs, qc_metric)

    flat_conn_mat = flatten_conn_matrices(bs);
    corrs = corr(flat_conn_mat, qc_metric);
    figure; histogram(corrs)
    fprintf('Mean (SD) corr = %3.2f (%3.2f)\n', mean(corrs), std(corrs));
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function gtdata = breakGts(gtdata)
% normalize IDs
[~, ~, ic] = unique(gtdata(:,2));
gtdata(:,2) = ic;
% break the groundtruth trajetory with multiple object categories
tracks = unique(ic);
newtrack_id = max(tracks) + 1;
for i = 1:numel(tracks)
    idx = gtdata(:, 2) == tracks(i);
    cls = gtdata(idx, 8);
    allcls = unique(cls);
    if(numel(allcls) > 1)
        for j = 2:numel(allcls)
            idxCls = idx & gtdata(:, 8) == allcls(j);
            gtdata(idxCls, 2) = newtrack_id;
            newtrack_id = newtrack_id + 1;
        end
    end
end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function [ deltaHD , deltaSEN , deltaPRO ] = CV2XMode4_common( lambda , Pt , distance, Psen , step_dB , noise , coding);

% CV2XMode4_common is a script that calculates the probability of packet
% loss due to half-duplex transmissions, the probability of packet loss
% due to a received signal power below the sensing power threshold and the
% probability of packet loss due to propagation effects for different
% Tx-Rx distances based on the models described in the following paper:
%
%    Manuel Gonzalez-Martín, Miguel Sepulcre, Rafael Molina-Masegosa, Javier Gozalvez,
%    "Analytical Models of the Performance of C-V2X Mode 4 Vehicular Communications",
%    IEEE Transactions on Vehicular Technology, Vol. 68, Issue 2, Feb. 2019. DOI: 10.1109/TVT.2018.2888704
%    Final version available at: https://ieeexplore.ieee.org/document/8581518
%    Post-print version available at: https://arxiv.org/abs/1807.06508
%
% CV2XMode4_common is called from the main script, CV2XMode4.
%
% Input parameters:
%   lambda: packet transmission frequency in Hz. .
%   Pt: transmission power in dBm.
%   distance: distance between tramsmitter and receiver in meters. It can be a vector with multiple distances.
%   Psen: sensing threshold in dBm.
%   step_dB: discrete steps to compute the PDF of the SNR and SINR in dB.
%   noise: noise corresponding to the DATA field of each message. Assumes a noise figure of 9dB and 10MHz channel (background noise of -95dBm). The total number of RBs in 10MHz is 50.
%   coding: ID of the coding used to identify the BLER curve
%
% Output metrics:
%    deltaHD: probability of packet loss due to half-duplex transmissions for different Tx-Rx distances
%    deltaSEN: probability of packet loss due to a received signal power below the sensing power threshold for different Tx-Rx distances
%    deltaPRO: probability of packet loss due to propagation effects for different Tx-Rx distances
%
% The equations that are identified with a number between brackets in this script are the ones
% that also appear in the paper so that they can be easily identified.

    D = length(distance);

    deltaHD(1:D) = lambda/1000;  % Equation (7)

    [PL_E_R, std_dev_E_R] = get_PL_SH(distance);   % Obtains pathloss and shadowing for different Tx-Rx distances.
    deltaSEN = 0.5 * ( 1 - erf( (Pt - PL_E_R - Psen)./(std_dev_E_R*sqrt(2)) ) );         % Equation (10)

    [SNR, PDF_SNR] = get_SINRdistribution( Pt - PL_E_R , -180 , std_dev_E_R , 3 , noise , Psen , step_dB); % Obtains the PDF of the SNR experienced by Rx (without interference).
    deltaPRO = get_BLER( SNR , PDF_SNR , coding , step_dB);   % Equation (13)

end

author: msepulcre


text:
function [xs,ys] = draw_Gaussian_1D(Mu,Var)
%DRAW_1D_GAUSSIAN Draws a Gaussian function onto axis.
%
%
%   input -----------------------------------------------------------------
%
%       o haxes, axis handle
%
%       o Mu & Sigma of Gaussian

Mu = Mu(:);

D = size(Mu,1);

if D == 1

    Max = Mu + 4.0*sqrt(Var);
    Min = Mu - 4.0*sqrt(Var);
    xs=linspace(Min,Max, 400);
    ys = normpdf(xs,Mu,sqrt(Var));
else

    disp('Not implemented yet');

end

end

author: epfl-lasa


text:
function [vs, output, v0] = C13ConfidenceInterval(v0, expdata, model, max_score, directions, majorIterationLimit)
%
%
% USAGE:
%
%    [vs, output, v0] = C13ConfidenceInterval(v0, expdata, model, max_score, directions, majorIterationLimit)
%
% INPUTS:
%    v0:                     set of flux vectors to be used as initial guesses.  They may be valid or not.
%    expdata:                experimental data
%    model:                  the standard model. Additional field .N (= null(S)) should also
%                            be provided. This is a basis of the flux space.
%    max_score:              maximum allowable data fit error
%
% OPTIONAL INPUTS:
%    directions:             ones and zeros of which reactions to compute (size = `n` x 1) or numbers of reactions to use  aka. [1; 5; 7; 8; 200] or
%                            reaction strings  aka.  {'GPK', 'PGL'}.  Ratios are possible with this input only.  Default = [] meaning  - do FVA with no ratios.
%    majorIterationLimit:    default = 10000
%
% OUTPUT:
%    vs:                     matrix
%    output:                 structure
%    v0:                     as in input

if nargin < 5
    directions = ones(size(v0,1),1);
end
if isempty(directions)
    directions = ones(size(v0,1),1);
end
t_start = clock;
printLevel = 3;
if nargin < 6
    majorIterationLimit = 1000;  %max number of iterations
end
diffInterval = 1e-5;         %gradient step size.
feasibilityTolerance = max_score/20; % how close you need to be to the max score.
logdirectory = strcat('temp', filesep);

% isratio(i):  >0 indicates not a ratio and points to reaction#
%              <0 indicates -numerator of ratio fraction.  In this case,
%              denom(i) stores the denominator
%
denom = zeros(size(directions));
if isnumeric(directions) %cannot be an actual ratio
    if max(directions) == 1
        isratio = find(directions);
    else
        isratio = directions;
    end
else % might be a ratio.  Gotta process strings
    isratio = zeros(size(directions));
    for i = 1:length(directions)
        if findstr(directions{i}, '/')
            [rxn1,rest] = strtok(directions{i}, '/');
            rxnID = findRxnIDs(model,rxn1);
            if rxnID == 0
                error('Unable to process rxn from list\n %s', directions{i});
            else
                isratio(i) = -rxnID;
            end
            rxn2 = rest(2:end);
            rxnID = findRxnIDs(model,rxn2);
            if rxnID == 0
                error('Unable to process rxn from list\n %s', rest(2:end))
            else
                denom(i) = rxnID;
            end
        else
            rxnID = findRxnIDs(model,directions{i});
            if rxnID == 0
                error('Unable to process rxn from list\n %s', directions{i});
            else
                isratio(i) = rxnID;
            end
        end
    end
end

numdirections = length(isratio);
numpoints = size(v0,2);

numiterations = numdirections*numpoints*2; % total number of iterations.

x0 = model.N\v0; % back substitute
scores = zeros(numpoints,1);

tProb.user.expdata = expdata;
tProb.user.model = model;
for i = 1:numpoints
    scores(i) = errorComputation2(x0(:,i),tProb);
end

% fit points if they are not currently feasible
v0(:,scores> max_score) = fitC13Data(v0(:,scores > max_score),expdata,model, majorIterationLimit);

if ~isfield(model, 'N')
   model.N = null(model.S);
end

x0 = model.N\v0; % back substitute

% safety check:
if (max(abs(model.S*v0))> 1e-6)
    display('v0 not quite in null space');
    pause;
end
if(max(abs(model.N*x0 - v0)) > 1e-6)
    display('null basis is weird');
    pause;
end

Name = 't2';
nalpha = size(model.N, 2);

x_L = -1000*ones(nalpha,1);
x_U = 1000*ones(nalpha,1);
[A, b_L, b_U] = defineLinearConstraints(model);

scores = zeros(numpoints,1);
% compute scores for all points.
for i = 1:numpoints
    scores(i) = errorComputation2(x0(:,i),tProb);
end
valid_index = scores < max_score + feasibilityTolerance;
fprintf('found %d valid points\n', sum(valid_index));

x0_valid = x0(:,valid_index);
x0_invalid = x0;
scores_valid = scores(valid_index);

% pre-compute unnecesary directions.
% if checkedbefore(i) ~= 0 then direction i is redundant
% if checkedbefore(i) = j then direction i and j are identical and do not
% need to be recomputed.
% checkedbefore(i) = j < 0 means that direction j is the same as i except
% for a sign switch.

checkedbefore = zeros(length(isratio),1);
for i = 2:length(isratio)
    if(isratio(i) < 0) % meaning it actually IS a ratio and no simplification possible
        continue
    end
    d = objective_coefficient(isratio(i), model);
    for j = 1:i-1
        if(isratio(j) < 0) % meaning it actually IS a ratio and no simplification possible
            continue
        end
        dj = objective_coefficient(isratio(j), model);
        if max(abs(dj - d))< 1e-4
            checkedbefore(i) = j;
            break;
        elseif max(abs(dj + d))< 1e-4
            checkedbefore(i) = -j;
            break;
        end
    end
end

% initialize variables;
outputv = 222*ones(numiterations,1);
outputexitflag = -222*ones(numiterations,1);
outputfinalscore = -222*ones(numiterations,1);
outputstruct = cell(numiterations,1);

csense = '';
for mm = 1:length(b_L),csense(mm,1) = 'L';end
for mm = 1:length(b_L),csense(mm+length(b_L),1) = 'G';end

fLowBnds = zeros(length(isratio), 2); %initialize but fill in later.
for rxn = 1:length(isratio)
    for direction = -1:2:1
        if isratio(rxn) < 0
            ration = objective_coefficient(-isratio(rxn),model);
            ratiod = objective_coefficient(denom(rxn),model);

            % in case RXN is a ratio
            Result = solveCobraNLP(...
                struct('lb', x_L, 'ub', x_U,...
                'name', Name,...
                'A', A,...
                'b_L', b_L, 'b_U', b_U,...
                'objFunction', 'ratioScore', 'g', 'ratioScore_grad',...
                'userParams', struct(...
                     'ration', direction*ration, 'ratiod', ratiod,... % set direction here too.
                     'diff_interval', diffInterval,'useparfor', true)...
                ),...
                'printLevel', 1, ...
                'iterationLimit', 1000);
        else
            d = objective_coefficient(isratio(rxn),model);
            Result = solveCobraLP(...
             struct('A', [A;A],'b',[b_U;b_L],'csense', csense, ...
                    'c', direction*d, ...
                    'lb', x_L,'ub', x_U, ...
                    'osense', 1),...
                    'feasTol',1e-7,'optTol',1e-7);
            if Result.stat ~= 1
                Result
                pause
            end
        end
        fLowBnds(rxn, (direction+3)/2) = direction*Result.obj; % fill in
    end
end

if ~exist (logdirectory, 'dir')
    if ~mkdir(logdirectory)
        error('unable to create logdirectory');
    end
end
clear d
%iterate through directions
parfor itnum = 1:numiterations
    if exist('ttt.txt', 'file') % abort w/o crashing if file 'ttt.txt' found in current directory
        fprintf('quitting due to file found\n');
        continue;
    end
    [rxn, direction, point] = getValues(itnum, numpoints); % translate itnum to rxn, direction and point
    % direction == 1 means minimize.  direction == -1 means maximize
    % (opposite of what you might think.

    if checkedbefore(rxn) ~= 0 %if this reaction maps to a previous reaction, we can skip
        continue;
    end

    fLowBnd = fLowBnds(rxn, (direction+3)/2); %get the absolute bound in the space w/o regards to C13 constraints.
    fprintf('reaction %d of %d, direction %d, lowerbound %f point %d of %d\n', rxn ,length(isratio), direction, fLowBnd, point, numpoints);
     % short circuit if x0 already close to a bound.
    if isratio(rxn) > 0
        di = objective_coefficient(isratio(rxn),model);
        obj1 = di'*x0_valid;
    else
        rationi = objective_coefficient(-isratio(rxn),model);
        ratiodi = objective_coefficient(denom(rxn),model);
        obj1 = (rationi'*x0_valid) ./ (ratiodi'*x0_valid);
    end

    if(any(abs(obj1-fLowBnd)<.0001))
        display('short circuiting');
        [nil, min_index] = min(obj1);
        outputv(itnum,1) = fLowBnd; % multiply by direction to correct sign.
        outputexitflag(itnum,1) = 111;
        outputfinalscore(itnum,1) = scores_valid(min_index);

    else % gotta actually do the computation.
        xinitial = x0_invalid(:,point);
        if isratio(rxn) > 0
            NLPsolution = solveCobraNLP(...
                struct('x0', xinitial, ...
                'lb', x_L, 'ub', x_U,...
                'name', Name,...
                'A', A, 'b_L', b_L, 'b_U', b_U,...
                'd', 'errorComputation2', 'dd', 'errorComputation2_grad',...
                'd_L', 0, 'd_U', max_score,...
                'c', di*direction, ... % direction of optimization
                'userParams', struct(...
                     'expdata', expdata,'model', model,'max_error', max_score,...
                     'diff_interval', diffInterval,'useparfor', true)...
                ),...
                'printLevel', printLevel, ...
                ...%'intTol', 1e-7, ...
                'iterationLimit', majorIterationLimit, ...
                'logFile', strcat(logdirectory, 'ci_', num2str(rxn),'x',num2str(direction),'x', point, '.txt'));
        else
            NLPsolution = solveCobraNLP(...
                struct('x0', xinitial, ...
                'lb', x_L, 'ub', x_U,...
                'name', Name,...
                'A', A, 'b_L', b_L, 'b_U', b_U,...
                'd', 'errorComputation2', 'dd', 'errorComputation2_grad',...
                'd_L', 0, 'd_U', max_score,...
                'objFunction', 'ratioScore', 'g', 'ratioScore_grad',...
                'userParams', struct(...
                     'expdata', expdata,'model', model,'max_error', max_score,...
                     'ration', direction*rationi,...
                     'ratiod', ratiodi,...
                     'diff_interval', diffInterval,'useparfor', false)...
                ),...
                'printLevel', printLevel, ...
                ...%'intTol', 1e-7, ...
                'iterationLimit', majorIterationLimit, ...
                'logFile', strcat(logdirectory, 'ci_', num2str(rxn),'x',num2str(direction),'x', point, '.txt'));
        end

        tscore = errorComputation2(NLPsolution.full, tProb);
        tbest = NLPsolution.obj;

        fprintf('reaction %d (%d), x %d; x=%f (%f); score=%f (%f)\n', rxn, length(isratio),direction, tbest,fLowBnd, tscore, max_score)

        outputv(itnum,1) = direction*tbest; % multiply by direction to correct sign.
        outputexitflag(itnum,1) = NLPsolution.origStat;
        outputfinalscore(itnum,1) = tscore;
        outputstruct{itnum,1} = NLPsolution;
    end
end

for itnum = 1:numiterations
    [rxn, direction, point] = getValues(itnum,numpoints);
    if direction == 1;
        output.minv(rxn,point) = outputv(itnum);
        output.minexitflag(rxn,point) = outputexitflag(itnum);
        output.minfinalscore(rxn,point) = outputfinalscore(itnum);
        output.minstruct(rxn,point) = outputstruct(itnum);
    else
        output.maxv(rxn,point) = outputv(itnum);
        output.maxexitflag(rxn,point) = outputexitflag(itnum);
        output.maxfinalscore(rxn,point) = outputfinalscore(itnum);
        output.maxstruct(rxn,point) = outputstruct(itnum);
    end
end

for i = 1:length(isratio)
    if checkedbefore(i) > 0 %short circuit if seen before.
        output.minv(i) = output.minv(checkedbefore(i));
        output.maxv(i) = output.maxv(checkedbefore(i));
        output.minexitflag(i) = output.minexitflag(checkedbefore(i));
        output.maxexitflag(i) = output.maxexitflag(checkedbefore(i));
        output.minfinalscore(i) = output.minfinalscore(checkedbefore(i));
        output.maxfinalscore(i) = output.maxfinalscore(checkedbefore(i));
        output.minstruct(i) = output.minstruct(checkedbefore(i));
        output.maxstruct(i) = output.maxstruct(checkedbefore(i));
    elseif checkedbefore(i) < 0
        output.minv(i) = -output.maxv(-checkedbefore(i));
        output.maxv(i) = -output.minv(-checkedbefore(i));
        output.minexitflag(i) = output.maxexitflag(-checkedbefore(i));
        output.maxexitflag(i) = output.minexitflag(-checkedbefore(i));
        output.minfinalscore(i) = output.maxfinalscore(-checkedbefore(i));
        output.maxfinalscore(i) = output.minfinalscore(-checkedbefore(i));
        output.minstruct(i) = output.maxstruct(-checkedbefore(i));
        output.maxstruct(i) = output.minstruct(-checkedbefore(i));
    end
end

vs = zeros(length(isratio), 2);
for i = 1:length(isratio)
    validindex = output.minfinalscore(i,:) < max_score + feasibilityTolerance;
    if any(validindex)
        vs(i,1) = min(output.minv(i,validindex));
    else
        vs(i,1) = 222;
    end
    validindex = output.maxfinalscore(i,:) < max_score + feasibilityTolerance;
    if any(validindex)
        vs(i,2) = max(output.maxv(i,validindex));
    else
        vs(i,2) = -222;
    end
end

elapsed_time = etime(clock, t_start)
return;

% function [index] = getIndex(rxn, direction, point, numpoints)
% % point goes from 1 .. NUMPOINTS
% % rxn goes from 1 .. NUMRXNS
% % direction is -1 or 1
% % index goes from 1 to NUMPOINTS*NUMRXNS*2
%
% rxn = rxn - 1;
% point = point -1;
% direction = (direction + 1)/2;
%
%
% index = rxn*numpoints*2 + direction*numpoints + point;
% index = index+1;
%
% return;

function [rxn, direction, point] = getValues(index, numpoints)
% point goes from 1 .. NUMPOINTS
% rxn goes from 1 .. NUMRXNS
% direction is -1 or 1
% index goes from 1 to NUMPOINTS*NUMRXNS*2
index = index - 1;
point = mod(index, numpoints);
index = index - point;
index = index/numpoints;
direction = mod(index,2);
index = index - direction;
index = index/2;
rxn = index;

point = point +1; % remap to 1..NUMPOINTS
direction = direction*2-1; % remap to -1,1
rxn = rxn+1; % remap to 1 .. number of rxns;
return;

% function that returns the proper objective coefficient for each reaction
% takes into account the reversibility of reactinos etc.
function [d] =  objective_coefficient(i, model)
d = zeros(length(model.lb),1);
d(i) = 1;
if (model.match(i))
    d(model.match(i)) = -1;
end
d = (d'*model.N)'; % transform to null space;
return

author: opencobra

repo: cobratoolbox


text:
function b = r8gd_vxm ( n, ndiag, offset, a, x )

%*****************************************************************************80
%
%% R8GD_VXM multiplies a vector by a R8GD matrix.
%
%  Discussion:
%
%    The R8GD storage format is suitable for matrices whose only nonzero entries
%    occur along a few diagonals, but for which these diagonals are not all
%    close enough to the main diagonal for band storage to be efficient.
%
%    In that case, we assign the main diagonal the offset value 0.
%    Each successive superdiagonal gets an offset value 1 higher, until
%    the highest superdiagonal (the A(1,N) entry) is assigned the offset N-1.
%    Similarly, the subdiagonals are assigned offsets of -1 through -(N-1).
%
%    Now, assuming that only a few of these diagonals contain nonzeros,
%    then for the I-th diagonal to be saved, we stored its offset in
%    OFFSET(I), and its entries in column I of the matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer NDIAG, the number of diagonals of the matrix
%    that are stored in the array.
%    NDIAG must be at least 1, and no more than 2 * N - 1.
%
%    Input, integer OFFSET(NDIAG), the offsets for the diagonal storage.
%
%    Input, real A(N,NDIAG), the R8GD matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(N), the product X*A.
%
  b(1:n) = 0.0;

  for i = 1 : n
    for diag = 1 : ndiag
      j = i + offset(diag);
      if ( 1 <= j & j <= n )
        b(j) = b(j) + x(i) * a(i,diag);
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function b = c8mat_copy ( m, n, a )

%*****************************************************************************80
%
%% C8MAT_COPY copies a C8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the order of the matrix.
%
%    Input, complex A(M,N), the matrix.
%
%    Output, complex B(M,N), a copy of the matrix.
%
  b(1:m,1:n) = a(1:m,1:n);

  return
end


text:
function [um] = m2um(m)
% Convert length from meters to micrometers or microns.
% Chad A. Greene 2012
um = m*1000000;


text:
function [codebook,scheme] = code_ECOC(m,dist,distfct)
% Generate the codebook for multiclass classification with Error Correcting Output encoding if feasible.
%
% function coding the multiple classes of this classification
% model, using the Error Correcting Output Coding;
%
%   codebook = code_ECOC(m)
%   codebook = code_ECOC(m, dist)
%   codebook = code_ECOC(m, dist, distfct)
%
% a codebook is found such that the minimal distances
% between the 'm' different classes is larger than 'dist' according
% to the distance measure of 'distfct'. The default is 'dist' 2
% for the 'codedist_hamming' distance. Besides the minimal distance
% between class representations, similar binary classifiers are
% also avoided as these do not add reliability in the context of
% deterministic binary classifiers.
%
% A recursive backtracking implementation looks for a
% representation which fullfills the constraint. It can decide
% exhaustively if such a representation is feasable. This can take
% lots of memory and time when 'm' becomes large (>50).
%
%
%  see also:
%    code, code_OneVsOne, code_OneVsAll, code_MOC, codedist_hamming

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

% default
eval('distfct;','distfct=''codedist_hamming'';');
eval('dist;','dist=2'';');
nb = ceil(log2(m*dist));
codebook =[];

candidates = eps.*ones(nb,1);
while isempty(codebook),
  disp(['number of bits ' num2str(nb)]);
  if nb>2^(m-1), error('No such code feasable'); end
  [codebook,sc] = create_code(candidates, m, dist, distfct,[]);
  if isempty(codebook),
    nb=nb+1;
    candidates = eps.*ones(nb,1);
  else
    hd=inf;
    hdM = 0;
    for t1=1:size(codebook,1),    for t2=(t1+1):size(codebook,1),
	hd = min(hd,feval(distfct,codebook(t1,:), codebook(t2,:)));
	hdM = max(hdM,feval(distfct,codebook(t1,:), codebook(t2,:)));
    end; end

    if hd==0|hdM==size(codebook,2),
      candidates = sc;
      codebook=[]; disp('retry');
    end
  end
end

%
% output format, where 'b' stands for binary discriminator
% see also 'code' and 'codelssvm'
scheme = []; for i=1:nb, scheme = [scheme 'b']; end

function [code,shrunkcandidate,rc] = create_code(candidates, m, dist, distfct,foundcand)
%
% recursive called function
%

% base case
if isempty(candidates), code=[]; shrunkcandidate=[]; rc=0; return; end

% pick a candidate
[nb,nc] = size(candidates);
rc=ceil(rand*nc);
acode = candidates(:,rc);

% initate this candidate
% and remove from the candidate list
acode = (acode~=eps).*acode;
aicode = acode +(acode==0).*sign(rand(nb,1)-.5);
if sum(acode==0)==0,
  candidates = candidates(:,[1:(rc-1) (rc+1):nc]);
else
  while(acode==aicode),
    aicode = acode + (acode==0).*sign(rand(nb,1));
  end
end
aicode = aicode+(aicode==0).*eps;
acode = acode+(acode==0).*eps;

candidates = shrink(candidates, aicode, dist, distfct);
shrunkcandidate = shrink(acode, aicode, dist, distfct);

% recursion
if m-1>0,
  shrunkc = candidates;

  fprintf('R;');
  [newcode,shrunkcandidate2,cc] = create_code(candidates,m-1, dist, distfct,[foundcand aicode]);
  fprintf('O;');
  while isempty(newcode),
    if isempty(find(shrunkcandidate2)), code=[]; return; end
    disp('retry with left candidates');
    shrunkc = [shrunkc(:,1:(cc-1)) shrunkcandidate2  shrunkc(:,(cc+1):end)];
    [newcode,shrunkcandidate2,cc] = create_code(shrunkc, m, dist, distfct,foundcand);
  end
  code = [aicode newcode];
 else
  code = aicode;
end

shrunkcandidate = candidates;

function shrunkcandidates = shrinkr(candidates, aicode, dist, distfct)
% refine candidates according to dist
% and shrink list of candidates
%
% recursive algorithm: TAKE CARE many recursions needed

fprintf('r');
% end of recursion
if isempty(candidates),shrunkcandidates=[]; return; end
if size(candidates,2)==1 &sum(candidates==eps)==0,shrunkcandidates=[]; return; end

% recursive step
cand = candidates(:,1);
if feval(distfct, aicode', cand)<dist,
  %zi = find(cand==eps & aicode~=eps);
  zi = find(cand==eps);
  if ~isempty(zi),
    ncandn = [cand(1:(zi-1)); -1; cand((zi+1):end)];
    ncandp = [cand(1:(zi-1)); 1; cand((zi+1):end)];
    candidates = [candidates(:,2:end) ncandp ncandn];
  else
    candidates = candidates(:,2:end);
  end
  shrunkcandidates = shrink(candidates,aicode,dist,distfct);
else
  shrunkcandidates = [cand shrink(candidates(:,2:end),aicode,dist,distfct)];
end
fprintf('o');

function shrunkcandidates = shrink(candidates, aicode, dist, distfct)
% refine candidates according to dist
% and shrink list of candidates
%
% iteration with dynamical list

%aicode
%candidates
i =1;
nb = size(candidates,2);
while i<=nb,
  cand = candidates(:,i);
  if feval(distfct, aicode', cand)<dist,
    zi = find(cand==eps);
    if ~isempty(zi),
      ncandn = [cand(1:(zi-1)); -1; cand((zi+1):end)];
      ncandp = [cand(1:(zi-1)); 1; cand((zi+1):end)];
      [candidates(:,[1:(i-1) (i+1):end]) ncandp ncandn];
      candidates = [candidates(:,[1:(i-1) (i+1):end]) ncandp ncandn];
    else
      candidates(:,[1:(i-1) (i+1):end]);
      candidates = candidates(:,[1:(i-1) (i+1):end]);
    end
  else
    i=i+1;
  end
  nb = size(candidates,2);
end
shrunkcandidates = candidates;

author: peterhcharlton

repo: RRest


text:
function d=meansqtf(b,a)
%AVEPSPEC calculates the mean square transfer function for a filter D=(B,A)
%
% Inputs: B,A         Numerator and denominator filter coefficients.
%
% Output: D           The mean square transfer function of the filter B/A. This equals
%                     the average otuput power when the filter is fed with unit variance
%                     white noise.
%
%                     D may be obtained approximately by:
%                         N=1024; D=sum(filter(B,A,[1 zeros(1,N)]).^2)

% Since the power spectrum is the fourier transform of the autocorrelation, we can calculate
% the average value of pb/pa by taking the 0'th order term of the convolution of the autocorrelation
% functions associated with b and 1/a. Since b is an FIR filter, this convolution is
% a finite sum even though the autocorrelation function of 1/a is infinite in extent.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: meansqtf.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if length(a)==1
   d=(b(:)')*b(:);
else
   m=lpcar2ra(b(:)');
   m(1)=m(1)*0.5;
   d=2*lpcar2rr(a(:)',length(m)-1)*m';
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [stop reason] = stoppingcriterion(problem, x, options, info, last)
% Checks for standard stopping criteria, as a helper to solvers.
%
% function [stop reason] = stoppingcriterion(problem, x, options, info, last)
%
% Executes standard stopping criterion checks, based on what is defined in
% the info(last) stats structure and in the options structure.
%
% The returned number 'stop' is 0 if none of the stopping criteria
% triggered, and a (strictly) positive integer otherwise. The integer
% identifies which criterion triggered:
%  0 : Nothing triggered;
%  1 : Cost tolerance reached;
%  2 : Gradient norm tolerance reached;
%  3 : Max time exceeded;
%  4 : Max iteration count reached;
%  5 : Maximum number of cost evaluations reached;
%  6 : User defined stopfun criterion triggered.
%
% The output 'reason' is a string describing the triggered event.

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    stop = 0;
    reason = '';

    stats = info(last);

    % Target cost attained
    if isfield(stats, 'cost') && isfield(options, 'tolcost') && ...
       stats.cost <= options.tolcost
        reason = 'Cost tolerance reached. See options.tolcost.';
        stop = 1;
        return;
    end

    % Target gradient norm attained
    if isfield(stats, 'gradnorm') && isfield(options, 'tolgradnorm') && ...
       stats.gradnorm < options.tolgradnorm
        reason = 'Gradient norm tolerance reached. See options.tolgradnorm.';
        stop = 2;
        return;
    end

    % Alloted time exceeded
    if isfield(stats, 'time') && isfield(options, 'maxtime') && ...
       stats.time >= options.maxtime
        reason = 'Max time exceeded. See options.maxtime.';
        stop = 3;
        return;
    end

    % Alloted iteration count exceeded
    if isfield(stats, 'iter') && isfield(options, 'maxiter') && ...
       stats.iter >= options.maxiter
        reason = 'Max iteration count reached. See options.maxiter.';
        stop = 4;
        return;
    end

    % Alloted function evaluation count exceeded
    if isfield(stats, 'costevals') && isfield(options, 'maxcostevals') && ...
       stats.costevals >= options.maxcostevals
        reason = 'Maximum number of cost evaluations reached. See options.maxcostevals.';
        stop = 5;
    end

    % Check whether the possibly user defined stopping criterion
    % triggers or not.
    if isfield(options, 'stopfun')
        userstop = options.stopfun(problem, x, info, last);
        if userstop
            reason = 'User defined stopfun criterion triggered. See options.stopfun.';
            stop = 6;
            return;
        end
    end

end

author: andrewssobral

repo: lrslibrary


text:
function obj = updateTernaryPlotPro(obj, ternaryIndex)

    %-AXIS INDEX-%
    axIndex = obj.getAxisIndex(obj.State.Plot(ternaryIndex).AssociatedAxis);

    %-GET DATA STRUCTURES-%
    ternaryData = get(obj.State.Plot(ternaryIndex).Handle);
    axisData = get(obj.State.Plot(ternaryIndex).AssociatedAxis);
    figureData = get(obj.State.Figure.Handle);

    %-CHECK FOR MULTIPLE AXES-%
    [xsource, ysource] = findSourceAxis(obj, axIndex);

    %=========================================================================%
    %
    %-UPDATE TRACE PLOT-%
    %
    %=========================================================================%

    %-get plot data-%
    xData = ternaryData.XData;
    yData = ternaryData.YData;

    %-------------------------------------------------------------------------%

    %-set trace-%
    for t = 1:size(xData,2)

        %-------------------------------------------------------------------------%

        %-get new ternaryIndex-%
        if t > 1
            obj.PlotOptions.nPlots = obj.PlotOptions.nPlots + 1;
            ternaryIndex = obj.PlotOptions.nPlots;
        end

        %-------------------------------------------------------------------------%

        %-trace type-%
        obj.data{ternaryIndex}.type = 'scatterternary';
        obj.data{ternaryIndex}.subplot = sprintf('ternary%d', xsource+1);

        %-------------------------------------------------------------------------%

        %-set mode and properties for trace-%
        if ~strcmpi('none', ternaryData.Marker) && ~strcmpi('none', ternaryData.LineStyle)
            obj.data{ternaryIndex}.mode = 'lines+markers';
            obj.data{ternaryIndex}.marker = extractPatchMarker(ternaryData);

        elseif ~strcmpi('none', ternaryData.Marker)
            obj.data{ternaryIndex}.mode = 'markers';
            obj.data{ternaryIndex}.marker = extractPatchMarker(ternaryData);

        elseif ~strcmpi('none', ternaryData.LineStyle)
            obj.data{ternaryIndex}.mode = 'lines';
            obj.data{ternaryIndex}.line = extractPatchLine(ternaryData);

        else
            obj.data{ternaryIndex}.mode = 'none';
        end

        %-------------------------------------------------------------------------%

        %-convert from cartesian coordinates to trenary points-%
        yTernData = yData(:,t); yTernData(end+1) = yData(1,t);
        xTernData = xData(:,t); xTernData(end+1) = xData(1,t);

        aData = yTernData/sin(deg2rad(60));
        bData = 1 - xTernData - yTernData*cot(deg2rad(60));

        %-------------------------------------------------------------------------%

        %-set plot data-%
        obj.data{ternaryIndex}.a = aData;
        obj.data{ternaryIndex}.b = bData;

        %-------------------------------------------------------------------------%

        %-some trace properties-%
        obj.data{ternaryIndex}.name = ternaryData.DisplayName;
        obj.data{ternaryIndex}.showscale = false;
        obj.data{ternaryIndex}.visible = strcmp(ternaryData.Visible,'on');

        %-trace coloring-%
        faceColor = ternaryData.FaceColor;

        if isnumeric(faceColor)
            fillColor = sprintf('rgb(%f,%f,%f)', 255*faceColor);

        else
            cMap = figureData.Colormap;
            nColors = size(cMap,1);

            switch faceColor

                case 'none'
                    fillColor = 'rgba(0,0,0,0)';

                case {'flat', 'interp'}

                    switch ternaryData.CDataMapping

                        case 'scaled'
                            cMin = axisData.CLim(1);
                            cMax = axisData.CLim(2);

                            if strcmpi(faceColor, 'flat')
                                cData = ternaryData.ZData(1,t);
                            elseif strcmpi(faceColor, 'interp')
                                cData = max(ternaryData.ZData(:,t));
                            end

                            cData = max(min(cData, cMax), cMin);
                            cData = (cData - cMin)/diff(axisData.CLim);
                            cData = 1 + floor( cData*(nColors-1) );

                            fillColor = sprintf('rgb(%f,%f,%f)', 255*cMap(cData,:));

                        case 'direct'
                            fillColor = sprintf('rgb(%f,%f,%f)', 255*cMap(ternary(1,t),:));
                    end
            end
        end

        obj.data{ternaryIndex}.fillcolor = fillColor;
        obj.data{ternaryIndex}.fill = 'toself';

        %---------------------------------------------------------------------%

        %-trace legend-%
        obj.data{ternaryIndex}.showlegend = false;

        %-------------------------------------------------------------------------%
    end

    %=========================================================================%
    %
    %-UPDATE TERNARY AXES-%
    %
    %=========================================================================%

    %-set domain plot-%
    xo = axisData.Position(1);
    yo = axisData.Position(2);
    w = axisData.Position(3);
    h = axisData.Position(4);

    ternary.domain.x = min([xo xo + w],1);
    ternary.domain.y = min([yo yo + h],1);

    %-----------------------------------------------------------------------------%

    %-label settings-%
    l = 1; t = 1;
    labelLetter = {'b', 'a', 'c'};

    for n = 1:length(axisData.Children)
        if strcmpi(axisData.Children(n).Type, 'text')
            stringText = axisData.Children(n).String;

            if any(isletter(stringText))
                labelIndex(l) = n;
                l = l + 1;
            else
                tickIndex(t) = n;
                t = t + 1;
            end
        end
    end

    for l = 1:length(labelIndex)
        n = labelIndex(l);
        patterText = sprintf('ternary.%saxis.title', labelLetter{l});

        labelText = axisData.Children(n).String;
        labelFontColor = sprintf('rgb(%f,%f,%f)', axisData.Children(n).Color);
        labelFontSize = 1.5 * axisData.Children(n).FontSize;
        labelFontFamily = matlab2plotlyfont(axisData.Children(n).FontName);

        eval(sprintf('%s.text = labelText;', patterText));
        eval(sprintf('%s.font.color = labelFontColor;', patterText));
        eval(sprintf('%s.font.size = labelFontColor;', patterText));
        eval(sprintf('%s.font.family = labelFontFamily;', patterText));
    end

    %-----------------------------------------------------------------------------%

    %-tick settings-%
    t0 = tickIndex(1); t1 = tickIndex(2);
    tick0 = str2num(axisData.Children(t0).String);
    tick1 = str2num(axisData.Children(t1).String);
    dtick = tick1 - tick0;

    tickFontColor = sprintf('rgb(%f,%f,%f)', axisData.Children(t0).Color);
    tickFontSize = 1.0 * axisData.Children(t0).FontSize;
    tickFontFamily = matlab2plotlyfont(axisData.Children(t0).FontName);

    for l = 1:3
        patterText = sprintf('ternary.%saxis', labelLetter{l});

        eval(sprintf('%s.tick0 = tick0;', patterText));
        eval(sprintf('%s.dtick = dtick;', patterText));
        eval(sprintf('%s.tickfont.color = tickFontColor;', patterText));
        eval(sprintf('%s.tickfont.size = tickFontSize;', patterText));
        eval(sprintf('%s.tickfont.family = tickFontFamily;', patterText));
    end

    %-----------------------------------------------------------------------------%

    %-set ternary axes to layout-%
    obj.layout = setfield(obj.layout, sprintf('ternary%d', xsource+1), ternary);

    %-----------------------------------------------------------------------------%

    obj.PlotlyDefaults.isTernary = true;

    %-----------------------------------------------------------------------------%
end

function rad = deg2rad(deg)
    rad = deg / 180 * pi;
end

author: plotly


text: %POWER Mapping overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% RES = pointOp(IM, LUT, ORIGIN, INCREMENT, WARNINGS)
%
% Apply a point operation, specified by lookup table LUT, to image IM.
% LUT must be a row or column vector, and is assumed to contain
% (equi-spaced) samples of the function.  ORIGIN specifies the
% abscissa associated with the first sample, and INCREMENT specifies the
% spacing between samples.  Between-sample values are estimated via
% linear interpolation.  If WARNINGS is non-zero, the function prints
% a warning whenever the lookup table is extrapolated.
%
% This function is much faster than MatLab's interp1, and allows
% extrapolation beyond the lookup table domain.  The drawbacks are
% that the lookup table must be equi-spaced, and the interpolation is
% linear.

% Eero Simoncelli, 8/96.

function res = pointOp(im, lut, origin, increment, warnings)

%% NOTE: THIS CODE IS NOT ACTUALLY USED! (MEX FILE IS CALLED INSTEAD)

fprintf(1,'WARNING: You should compile the MEX version of "pointOp.c",\n         found in the MEX subdirectory of matlabPyrTools, and put it in your matlab path.  It is MUCH faster.\n');

X = origin + increment*[0:size(lut(:),1)-1];
Y = lut(:);

res = reshape(interp1(X, Y, im(:), 'linear', 'extrap'),size(im));

repo: SelfExSR


text:
% Fast Inter-Harmonic Reconstruction using the zero-phase version of the
% speech signal, as a pre-process for LPC-based spectral envelope estimation.
%
%
% Description
%  This technique reconstructs inter-harmonics of the voice signal,
% as a pre-process for an efficient spectral envelope estimation in
% high-pitched voices. The technique is fully described in [1].

%
% Inputs
%  wave            : [samples] [Nx1] input signal (speech signal)
%  Fs              : [Hz]      [1x1] sampling frequency
%  f0              : F0 estimate (a 0 value is for an unvoiced frame)
%  order           : Order used for the LP analysis
%
%
% Outputs
%  LP              : LP coefficients
%  Energy          : Energy of the frames
%
%
% Example
%  Please see the HOWTO_envelope.m example file.
%  Please see http://tcts.fpms.ac.be/~drugman/Toolbox/ for more details.
%
% References
%  [1] T.Drugman, Y. Stylianou, "Fast Inter-Harmonic Reconstruction for
%      Spectral Envelope Estimation in High-Pitched Voices", vol. 21, pp.
%      1418-1422, IEEE Signal Processing Letters, 2014.
%      http://tcts.fpms.ac.be/~drugman/files/SPL-FIHR.pdf
%
% Copyright (c) 2014 Toshiba Cambridge Research Laboratory
%
% License
%  This code will be part of the GLOAT toolbox with the following
%  licence:
%  This program is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
% This function will also be part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Thomas Drugman thomas.drugman@umons.ac.be

function [LP, Energy] = env_fihrzp(Seg,Fs,F0_local,order)

F0_target=100;
% F0* in the paper. Here set to 100 Hz.

Energy=sum(Seg.^2);

if F0_local>0

    %Zero-phase version of the input signal
    Nfft=2048;
    ZP=real(ifft(abs(fft(Seg',Nfft))));
    ZP=[ZP(Nfft/2+1:end) ZP(1:Nfft/2)]';
    Seg=ZP(round(length(ZP)/2)-round(length(Seg)/2):round(length(ZP)/2)+round(length(Seg)/2)-1);

    % The rest is similar to the FIHR technique
    TmpSignal=ones(1,length(Seg));
    t_tmp=1:length(TmpSignal);

    I=ceil(F0_local/F0_target);
    for k=1:I-1
        W=(I-k)/I;
        % We use here a linear weighting function. Other functions are
        % possible, as long as they meet the properties mentioned in
        % the paper.
        TmpSignal=TmpSignal+2*W*cos(2*pi*(k/I)*F0_local/Fs*t_tmp);

        % Equation (2) in the paper
    end

    Seg2=Seg.*TmpSignal';
else
    Seg2=Seg;
end

LP=lpc(Seg2,order);

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [distMat_cl,cluster_Groupi,ord,handles] = BF_ClusterDown(distMat,varargin)
% BF_ClusterDown    Reduce a pairwise distance matrix into smaller clusters.
%
% Yields a visualization of a pairwise distance matrix, including a set of
% clusters of objects showing highly correlated behavior.
%
%---INPUTS:
% distMat, a pairwise distance vector or matrix (e.g., generated from pdist)
% [EXTRA OPTIONS]:
% 'clusterThreshold', the threshold in the dendrogram for forming clusters
% 'objectLabels', labels for axes
% 'whatDistance', the type of distance metric used:
%        (i) 'corr' (default): provide 1-R, where R is correlation coefficient
%        (ii) 'abscorr': provide correlation distances, 1-R, but clustering is
%                        done on abscorr distances (sign of R is ignored), and
%                        correlations, R, are plotted
%        (iii) 'abscorr_ii': provide correlation distances, 1-R, but clustering
%                        is done on abscorr distances (sign of R is ignored), and
%                        absolute correlations, |R|, are plotted
%        (iv) 'general': provide any distance matrix, plots the distance matrix
% 'linkageMeth', the linkage method to use (for Matlab's linkage function)
%
%---OUTPUTS:
% distMat_cl, the distance matrix re-ordered by clustering
% cluster_Groupi, the assignment of nodes to clusters (cell of indices)
% ord, the ordering of objects in the dendrogram
% handles, structure containing handles to plotting elements, f (figure),
%               ax1 (dendrogram), ax2 (pairwise correlation matrix)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% Check Inputs
%-------------------------------------------------------------------------------
% Custom plotting options using an input parser
inputP = inputParser;

% Cluster threshold:
default_clusterThreshold = 0.2;
addParameter(inputP,'clusterThreshold',default_clusterThreshold,@isnumeric);

% Object labels for axes
default_objectLabels = {};
addParameter(inputP,'objectLabels',default_objectLabels,@iscell);

% Default input is a distance matrix based on (Spearman) correlation coefficients
default_whatDistance = 'corr';
addParameter(inputP,'whatDistance',default_whatDistance,@ischar);

% errth Set error threshold instead of just top number
% default_errTh = 30;
% addParameter(inputP,'errTh',default_errTh,@isnumeric);

% miOrCorr
% default_miOrCorr = 'corr'; % 'mi','corr'
% addParameter(inputP,'miOrCorr',default_miOrCorr,@ischar);

% clusterMeth, what method to use to do the clustering
% default_clusterMeth = 'linkage'; % 'linkage', 'kmedoids'
% addParameter(inputP,'clusterMeth',default_clusterMeth,@ischar);

% linkageMeth, what method to use to do the linkage clustering
default_linkageMeth = 'average';
addParameter(inputP,'linkageMeth',default_linkageMeth,@ischar);

% plotBig, whether to make a big plot
% default_plotBig = 0;
% check_plotBig = @(x) x==1 || x==0;
% addParameter(inputP,'plotBig',default_plotBig,check_plotBig);

%% Parse inputs:
parse(inputP,varargin{:});

% Make variables from results of input parser:
clusterThreshold = inputP.Results.clusterThreshold;
objectLabels = inputP.Results.objectLabels;
whatDistance = inputP.Results.whatDistance;
% errTh = inputP.Results.errTh;
% miOrCorr = inputP.Results.miOrCorr;
% clusterMeth = inputP.Results.clusterMeth;
linkageMeth = inputP.Results.linkageMeth;
% plotBig = inputP.Results.plotBig;
clear('inputP');

%-------------------------------------------------------------------------------
%% Process input distance data
%-------------------------------------------------------------------------------
[distVec,links,distVec0,numItems] = BF_DoLinkage(distMat,whatDistance,linkageMeth);

% ------------------------------------------------------------------------------
%% Compute dendrogram
% ------------------------------------------------------------------------------
f = figure('color','w');
f.Position(3:4) = [1200,800];

% Get the dendrogram reordering:
ax1 = subplot(1,6,6);
ord = BF_LinkageOrdering(distVec,links);
h_dend = dendrogram(links,0,'Orientation','right','Reorder',ord);
ax1.YDir = 'reverse'; % needs to be reversed to match the reversed y-axis of the imagesc plot

% Save the clustered distance matrix for output
distMat_cl = distMat(ord,ord);

% Make the dendrogram look aight
set(h_dend,'color','k','LineWidth',1)
ax1.YTickLabel = {};
xlabel('Distance');

%-------------------------------------------------------------------------------
%% Cluster into groups:
%-------------------------------------------------------------------------------
% Cluster the dendrogram:
T = cluster(links,'cutoff',clusterThreshold,'criterion','distance');
numClusters = max(T);

fprintf(1,'%s distance-based clustering with %u clusters.\n',whatDistance,numClusters);

% Reorder members of each cluster by distance to other members of the cluster:
cluster_Groupi = cell(numClusters,1);
for j = 1:numClusters
    cluster_Groupi{j} = find(T==j);
    if sum(T==j) > 1
        % Sort by increasing sum of distances to other members of the cluster
        [~,ix] = sort(sum(distMat(cluster_Groupi{j},cluster_Groupi{j})),'ascend');
        cluster_Groupi{j} = cluster_Groupi{j}(ix);
    end
end

% Reorder by decreasing cluster size
[~,ix] = sort(cellfun(@length,cluster_Groupi),'descend');
cluster_Groupi = cluster_Groupi(ix);

% Select the closest to cluster centre in each group
clusterCenters = cellfun(@(x) x(1),cluster_Groupi);

cluster_Groupi_cl = cellfun(@(x) arrayfun(@(y)find(ord==y),x),cluster_Groupi,...
                                                'UniformOutput',false);
clusterCenters_cl = arrayfun(@(y) find(ord==y),clusterCenters);

%-------------------------------------------------------------------------------
%% Plot it
%-------------------------------------------------------------------------------
% Make a square distance matrix, and plot as a pairwise similarity matrix
% (reordered by ord determined above)
ax2 = subplot(1,6,2:5);
switch whatDistance
case {'corr','abscorr_ii','corr_fast'}
    % Input is a correlation distance matrix, plot correlation coefficients:
    distMat0 = squareform(distVec0);
    BF_PlotCorrMat(1 - distMat0(ord,ord),'n1to1');
    rectangleColors = BF_GetColorMap('accent',5,true);
case 'abscorr'
    distMat0 = squareform(distVec0);
    % Input is correlation distance matrix, plot absolute correlation coefficients:
    BF_PlotCorrMat(abs(1 - distMat0(ord,ord)),'0to1');
    % (red colormap)
    rectangleColors = {[0,0,0],ones(1,3)*0.5};
case 'general'
    % Input is a general distance matrix (grayscale colormap):
    BF_PlotCorrMat(distMat_cl);
    rectangleColors = BF_GetColorMap('accent',5,true);
otherwise
    warning('No special plotting options for distance metric: ''%s''',whatDistance);
    BF_PlotCorrMat(distMat_cl);
    rectangleColors = BF_GetColorMap('accent',5,true);
end

%-------------------------------------------------------------------------------
%% Add rectangles to group highly correlated clusters:
%-------------------------------------------------------------------------------
for i = 1:numClusters
    % Label cluster borders:
    rectangle('Position',[min(cluster_Groupi_cl{i})-0.5,min(cluster_Groupi_cl{i})-0.5, ...
                            length(cluster_Groupi_cl{i}),length(cluster_Groupi_cl{i})], ...
                    'EdgeColor',rectangleColors{1},'LineWidth',3);

    % Label cluster centers:
    rectangle('Position',[clusterCenters_cl(i)-0.5,clusterCenters_cl(i)-0.5,1,1], ...
                                'EdgeColor',rectangleColors{2},'LineWidth',3);
end

%-------------------------------------------------------------------------------
%% Tweak plot details:
%-------------------------------------------------------------------------------

% Add object labels on y-axis if provided
if ~isempty(objectLabels)
    ax2.YTick = 1:numItems;
    ax2.YTickLabel = objectLabels(ord);
    ax2.TickLabelInterpreter = 'none';
else
    ax2.YTick = [];
end

% Remove X-ticks
ax2.XTick = [];

% Add a color bar:
cB = colorbar('southoutside');
switch whatDistance
case {'corr','abscorr_ii'}
    cB.Label.String = 'Spearman correlation, \rho';
case 'abscorr'
    cB.Label.String = 'Magnitude of Spearman correlation, |\rho|';
case 'general'
    cB.Label.String = 'distance';
end

% Link axes:
linkaxes([ax1,ax2],'y');

% Fix limits
ax2.YLim = [0.5,numItems+0.5];

% Align with dendrogram
ax1.Position(1) = ax2.Position(1) + ax2.Position(3);
ax1.Position(4) = ax2.Position(4);
ax1.Position(2) = ax2.Position(2);

% Handles
if nargout > 3
    handles = struct();
    handles.f = f;
    handles.ax1 = ax1;
    handles.ax2 = ax2;
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
classdef SettingsMonitoringDocker < AbstractSettings

    properties (Access = protected)
        defaultParamsName = 'paramsMonitoringDocker.json'
    end

    properties (Access = public)
        showOptParams
        refreshInterval
        shallDisplayDesignVar
        shallShowBoundaryConditions

        problemID
        costFuncNames
        costWeights
        constraintFuncs
        dim
        scale

        designVariable
        dualVariable
        optimizerNames
        cost
        constraint

        convergenceVars
        boundaryConditions

        mesh
    end

    methods (Access = public)

        function obj = SettingsMonitoringDocker(varargin)
            if nargin == 1
                obj.loadParams(varargin{1});
            end
            obj.init();
        end

    end

    methods (Access = private)

        function init(obj)
            if ischar(obj)
            end
        end

    end

end

author: SwanLab

repo: Swan