text:
function result=rsimpls(x,y,varargin)

%RSIMPLS is a 'Robust method for Partial Least Squares Regression based on the
% SIMPLS algorithm'. It can be applied to both low and high-dimensional predictor variables x
% and to one or multiple response variables y. It is resistant to outliers in the data.
% The RSIMPLS algorithm is built on two main stages. First, a matrix of scores is derived
% based on a robust covariance criterion (see robpca.m),
% and secondly a robust regression is performed based on the results from ROBPCA.
%
% The RSIMPLS method is described in:
%    Hubert, M., and Vanden Branden, K. (2003),
%    "Robust Methods for Partial Least Squares Regression",
%    Journal of Chemometrics, 17, 537-549.
%
% To select the number of components in the regression model, a robust RMSECV (root mean squared
% error of cross validation) curve is drawn, based on a fast algorithm for
% cross-validation. This approach is described in:
%
%    Engelen, S., Hubert, M. (2005),
%    "Fast model selection for robust calibration methods",
%    Analytica Chimica Acta, 544, 219-228.
%
% Required input arguments:
%            x : Data matrix of the explanatory variables
%                (n observations in rows, p variables in columns)
%            y : Data matrix of the response variables
%                (n observations in rows, q variables in columns)
%
% Optional input arguments:
%            k : Number of components to be used.
%                (default = min(rank([x,y]),kmax)). If k is not specified,
%                it can be selected using the option 'rmsecv'.
%         kmax : Maximal number of components to be used (default = 9).
%                If k is provided, kmax does not need to be specified, unless k is larger
%                than 9.
%        alpha : (1-alpha) measures the fraction of outliers the algorithm should
%                resist. Any value between 0.5 and 1 may be specified. (default = 0.75)
%            h : (n-h+1) measures the number of observations the algorithm should
%                resist. Any value between n/2 and n may be specified. (default = 0.75*n)
%                Alpha and h may not both be specified.
%       rmsecv : If equal to zero and k is not specified, a robust R-squared curve
%                is plotted and an optimal k value can be chosen. (default)
%                If equal to one and k is not specified, a robust component selection-curve is plotted and
%                an optimal k value can be chosen. This curve computes a combination of
%                the cross-validated root mean squared error and the robust residual
%                sum of squares for k=1 to kmax. Rmsecv and k may not both
%                be specified.
%        rmsep : If equal to one, the robust RMSEP-value (root mean squared error of
%                prediction) for the model with k components (default = 0). This value
%                is automatically given if rmsecv = 1.
%        plots : If equal to one, a menu is shown which allows to draw several plots,
%                such as a robust score outlier map and a regression
%                outlier map are drawn. (default)
%                If the input argument 'classic' is equal to one, the classical
%                diagnostic plots are drawn as well.
%                If 'plots' is equal to zero, all plots are suppressed.
%                See also makeplot.m
%        labsd : The 'labsd' observations with largest score distance are
%                labeled on the outlier map (default = 3)
%        labod : The 'labod' observations with largest orthogonal distance are
%                labeled on the outlier map (default = 3)
%      labresd : The 'labresd' observations with largest residual distance are
%                labeled on the outlier map (default = 3)
%      classic : If equal to one, the classical SIMPLS analysis will be performed as well
%                (see also csimpls.m). (default = 0)
%
% Options for advanced users:
%           kr : Total number of components used by the ROBPCA method.
%                We advise to use kr=k+q. (default)
%        kmaxr : Maximal number of components used by the ROBPCA method.
%                default = min(kmax+q,rank([x,y]))
%  plotsrobpca : If equal to one, a robust score outlier map from ROBPCA is drawn.
%                If the input argument 'classic' is equal to one, the classical
%                outlier map is drawn as well.
%                If 'plotsrobpca' is equal to zero, all plots are suppressed. (default)
%           st : Indicates the current stage of the algorithm for cross-validation (RMSECV)
%                default = 0 -> performs all the stages of the algorithm and computes all the
%                parameters for the full model.
%                st=1, robpca is still performed, but when
%                st=2, the algorithm proceeds based on the previous knownledge of the output from robpca.
%          out : Is an empty structure, but is constructed while performing the cross-validation.
%                (see rrmse.m)
%
% I/O: result=rsimpls(x,y,'k',k,'kmax',10,'alpha',0.75,'h',h,'rmsecv',0,'rmsep',0,...
%       'plots',1,'labsd',3,'labod',3,'labresd',3,'classic',1,'kr',kr,...
%       'kmaxr',kmaxr,'plotsrobpca',0,'st',0,'out',[]);
%  The user should only give the input arguments that have to change their default value.
%  The name of the input arguments needs to be followed by their value.
%  The order of the input arguments is of no importance.
%
% Examples:
%   rsimpls(x,y,'k',5,'plots',1);
%   rsimpls(x,y,'classic',1,'rmsecv',1);
%
% The output of RSIMPLS is a structure containing:
%
%   result.slope     : Robust slope estimate
%   result.int       : Robust intercept estimate
%   result.fitted    : Robust fitted values
%   result.res       : Robust residuals
%   result.cov       : Estimated variance-covariance matrix of the residuals
%   result.T         : Robust scores
%   result.weights.r : Robust simpls weights
%   result.weights.p : Robust simpls weights
%   result.Tcenter   : Robust center of the scores
%   result.Tcov      : Robust covariance matrix of the scores
%   result.rsquared  : Robust R-squared value for the optimal k
%   result.rcs       : Robust Component Selection Criterion:
%                      This is a matrix with kmax columns. The first row contains the approximate
%                      R-squared values (for k=1,...,kmax) and the second row the square root of the weighted
%                      residuals sum of squares. This is equal to the RCS-value with
%                      gamma = 0 (see Engelen and Hubert (2005) for the definition).
%                      If the input argument rmsecv = 1, the third row contains the RCS-value for
%                      gamma = 0.5 and the fourth for gamma = 1. The last one is equal to the robust
%                      cross-validated RMSE values.
%                      Note that all the entries in this matrix depend on the choice of kmax.
%   result.rmsep     : Robust RMSEP value
%   result.k         : Number of components used in the regression
%   result.h         : The quantile h used throughout the algorithm
%   result.sd        : Robust score distances
%   result.od        : Robust orthogonal distances
%   result.resd      : Residual distances (when there are several response variables).
%                      If univariate regression is performed, it contains the standardized residuals.
%   result.cutoff    : Cutoff values for the score (result.cutoff.sd), orthogonal
%                      (result.cutoff.od) and residual distances (result.cutoff.resd).
%                      We use 0.975 quantiles of the chi-squared distribution.
%   result.flag      : The observations whose score distance is larger than
%                      'result.cutoff.sd' receive a flag 'result.flag.sd' equal
%                      to zero (good leverage points). Otherwise 'result.flag.sd'
%                      is equal to one.
%                      The components 'result.flag.od' and 'result.flag.resd' are
%                      defined analogously, and determine the orthogonal outliers,
%                      resp. the bad leverage points/vertical outliers.
%                      The observations with 'result.flag.od' and 'result.flag.resd'
%                      equal to zero, can be considered as calibration outliers and receive
%                      'result.flag.all' equal to zero. The regular observations and the good leverage
%                      points have 'result.flag.all' equal to one.
%   result.class     : 'RSIMPLS'
%   result.classic   : If the inputargument 'classic' is equal to 1, this structure
%                      contains results of the classical SIMPLS analysis. (see also csimpls.m)
%   results.robpca   : The results of robpca on [X,Y].
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Karlien Vanden Branden
% Version date: 07/04/2004
% Last update: 04/08/2006

%
%initialization with defaults
%
if rem(nargin,2)~=0
    error('Number of input arguments must be even!');
end
[n,p1]=size(x);
[n2,q1]=size(y);
z=[x,y];
rx=rank(x);
rz=rank(z);
if n~=n2
    error('The response variables and the predictor variables have a different number of observations.')
end
niter=100;counter=1;mcd=0;alfa=0.75;
kmax=min([9,rx,floor(n/2),p1]);
kmaxr=min([kmax+q1,rz]);
h=floor(2*floor((n+kmaxr+1)/2)-n+2*(n-floor((n+kmaxr+1)/2))*alfa);
labsd=3;labod=3;labresd=3;
plotsrobpca=0;plots=1;k=0;
kr=k+q1;
st=0;rmsecv=0;
out=[];classic=0;rmsep=0;rmsep_value=nan;rmsecv_value = nan;rsquared_value = nan;rss_value = nan;
default=struct('alpha',alfa,'h',h,'labsd',labsd,'labod',labod,'labresd',labresd,'k',k,'kr',kr,...
    'plotsrobpca',plotsrobpca,'plots',plots,'kmax',kmax,'st',st,...
    'out',out,'rmsecv',rmsecv,'classic',classic,'rmsep',rmsep,...
    'kmaxr',kmaxr,'rmsep_value',rmsep_value,'rmsecv_value',rmsecv_value,'rsquared_value',...
    rsquared_value,'rss_value',rss_value);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
%
if nargin>2
    %
    %placing inputfields in array of strings
    %
    for j=1:nargin-3
        if rem(j,2)~=0
            chklist{i}=varargin{j};
            i=i+1;
        end
    end
    dummy=sum(strcmp(chklist,'h')+2*strcmp(chklist,'alpha'));
    switch dummy
        case 0 %Take on default values
            options.alpha=alfa;%0.75;
            options.h=h;
        case 3
            error('Both input arguments alpha and h are provided. Only one is required.')
    end
    %
    %Checking which default parameters have to be changed
    % and keep them in the structure 'options'.
    %
    while counter<=IN
        index=strmatch(list(counter,:),chklist,'exact');%contains the users input one by one
        if ~isempty(index) %in case of similarity
            for j=1:nargin-3 %searching the index of the accompanying field
                if rem(j,2)~=0 %fieldnames are placed on odd index
                    if strcmp(chklist{index},varargin{j})
                        I=j;
                    end
                end
            end
            options=setfield(options,chklist{index},varargin{I+1});
            index=[];
        end
        counter=counter+1;
    end
    options.h=floor(options.h);
    options.kmax=floor(options.kmax);
    options.k=floor(options.k);
    options.kmaxr=floor(options.kmaxr);
    options.kr=floor(options.kr);
    kmax=min([options.kmax,floor(n/2),rz,p1]);
    kmaxr=max([min([options.kmaxr,rz]),kmax+q1]);
    k=min(options.k,kmax);
    while k<0
        k=input(['The number of components can not be negative.\n'...
        'How many principal components would you like to retain?\n']);
    end
    if any(strcmp(chklist,'kr'))
        kr=floor(max([options.kr,k+q1]));
    else
        kr=k+q1;
    end
    if dummy==1 %checking inputvariable h
        if options.h-floor(options.h)~=0
            mess=sprintf('Attention (rsimpls.m): h must be an integer. \n');
            disp(mess)
        end
        if kr==0
            if options.h<floor((n+kmaxr+1)/2)
                options.h=floor((n+kmaxr+1)/2);
                mess=sprintf(['Attention (rsimpls.m): h should be larger than (n+kmaxr+1)/2.\n',...
                        'It is set to its minimum value',num2str(options.h)]);
                disp(mess)
            end
            options.alpha=options.h/n;
        else
            if options.h<floor((n+kr+1)/2)
                options.h=floor((n+kr+1)/2);
                mess=sprintf(['Attention (rsimpls.m): h should be larger than (n+kr+1)/2.\n',...
                        'It is set to its minimum value',num2str(options.h)]);
                disp(mess)
            end
            options.alpha=options.h/n;
        end
        if options.h>n
            options.alpha=0.75;
            if kr==0
                options.h=floor(2*floor((n+kmaxr+1)/2)-n+2*(n-floor((n+kmaxr+1)/2))*options.alpha);
            else
                options.h=floor(2*floor((n+kr+1)/2)-n+2*(n-floor((n+kr+1)/2))*options.alpha);
            end
            mess=sprintf(['Attention (rsimpls.m): h should be smaller than n. \n',...
                    'It is set to its default value ',num2str(options.h)]);
            disp(mess)
        end
    elseif dummy==2
        if options.alpha < 0.5
            options.alpha=0.5;
            mess=sprintf(['Attention (rsimpls.m): Alpha should be larger than 0.5. \n',...
                    'It is set to 0.5.']);
            disp(mess)
        end
        if options.alpha > 1
            options.alpha=0.75;
            mess=sprintf(['Attention (rsimpls.m): Alpha should be smaller than 1.\n',...
                    'It is set to 0.75.']);
            disp(mess)
        end
        if kr==0
            options.h=floor(2*floor((n+kmaxr+1)/2)-n+2*(n-floor((n+kmaxr+1)/2))*options.alpha);
        else
            options.h=floor(2*floor((n+kr+1)/2)-n+2*(n-floor((n+kr+1)/2))*options.alpha);
        end
    end
    h=options.h;alfa=options.alpha;labsd=max(0,min(floor(options.labsd),n));
    dummyh = strcmp(chklist,'h');
    dummykmax = strcmp(chklist,'kmax');
    if all(dummyh == 0) && any(dummykmax)
        h = floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*alfa);
    end
    labod=max(0,min(floor(options.labod),n));labresd=max(0,min(floor(options.labresd),n));
    plotsrobpca=options.plotsrobpca;plots=options.plots;
    st=options.st;
    out=options.out;
    rmsecv=options.rmsecv;
    classic=options.classic;
    rmsep=options.rmsep;
    rmsep_value = options.rmsep_value;
    rmsecv_value = options.rmsecv_value;
    rmsecv = options.rmsecv;
    rsquared_value = options.rsquared_value;
    rss_value = options.rss_value;
end

if q1==1 && k>=(h-2)
    mess=sprintf(['Attention (rsimpls.m): The number of components, k = ',num2str(k),...
            '\n is larger than our recommended maximum value of k = ',num2str(h-2)-1,'.']);
    disp(mess)
elseif q1>1 && k>=((h/q1)-(q1/2)-0.5)
    mess=sprintf(['Attention (rsimpls.m): The number of components, k = ',num2str(k),...
            '\n is larger than our recommended maximum value of k = ',num2str(floor((h/q1)-(q1/2)-1.5)),'.']);
    disp(mess)
end

%
%MAIN PART
%
% selection of number of components
if k == 0
    if rmsecv == 0
        [R2,final]=rsquared(x,y,kmax,'RSIMPLS',options.h);
        rss = final.rss;
        k = final.k;
        result=rsimpls(x,y,'k',k,'kr',k+q1,'h',h,'rmsecv',0,'rmsep',options.rmsep,'plots',0,'classic',classic,'rsquared_value',R2,'rss_value',rss);
    else
        out=rrmse(x,y,h,kmax,'RSIMPLS',1);
        R2 = out.R2;
        rss = out.rss;
        k=out.k;
        rmsecv_value = out.rmsecv;
        pred = rrmse(x,y,h,kmax,'RSIMPLS',0,k,out.weight,out.res);
        result=rsimpls(x,y,'k',k,'kr',k+q1,'h',h,'rmsecv',0,'rmsep',0,'plots',0,'classic',classic,'rmsep_value',pred.rmsep,'rmsecv_value',rmsecv_value,'rsquared_value',R2,'rss_value',rss);
    end
else
    if rmsecv
        error(['Both RMSECV and k were given.', ...
            'Please rerun your analysis with one of these inputs.  (see help file)'])
    end
%First stage: Obtain the scores T by first performing ROBPCA on z:
    if st<=1
        out.robpca=robpca(z,'k',kr,'h',h,'plots',plotsrobpca,'kmax',kmaxr,'classic',classic,'mcd',mcd);
        out.h=h;
        out.centerz=out.robpca.M;
        out.sigmaxy=out.robpca.P(1:p1,:)*diag(out.robpca.L)*out.robpca.P(p1+1:p1+q1,:)';
        out.sigmax=out.robpca.P(1:p1,:)*diag(out.robpca.L)*out.robpca.P(1:p1,:)';
        out.xcentr=x-repmat(out.centerz(1:p1),n,1);
        out.ycentr=y-repmat(out.centerz(p1+1:p1+q1),n,1);
        out.weights2=out.robpca.flag.all;
    end
    if st
        i=k;
    else
        i=1;
    end
    while i<=k
        out.sigmayx=out.sigmaxy';
        if q1>p1
            [RR,LL]=eig(out.sigmaxy*out.sigmayx);
            [LL,I]=greatsort(diag(LL));
            rr=RR(:,I(1));
            qq=out.sigmayx*rr;
            qq=qq/norm(qq);
        else
            [QQ,LL]=eig(out.sigmayx*out.sigmaxy);
            [LL,I]=greatsort(diag(LL));
            qq=QQ(:,I(1));
            rr=out.sigmaxy*qq;
            rr=rr/norm(rr);
        end
        tt=out.xcentr*rr;
        uu=out.ycentr*qq;
        pp=out.sigmax*rr/(rr'*out.sigmax*rr);
        vv=pp;
        if i>1
            vv=vv-out.v*(out.v'*pp);
        end
        if vv'*vv==0
            error('The number of components is too large')
        end
        vv=vv./norm(vv);
        out.sigmaxy=out.sigmaxy-vv*(vv'*out.sigmaxy);
        out.v(:,i)=vv;
        out.q(:,i)=qq;
        out.t(:,i)=tt;
        out.u(:,i)=uu;
        out.p(:,i)=pp;
        out.r(:,i)=rr;
        i=i+1;
    end

    %Second Stage : Robust ROBPCA-regression
    robpcareg=robpcaregres(out.t,y,out.weights2);
    breg=robpcareg.coeffs(1:k,:);
    Yhat=out.t*breg+repmat(robpcareg.coeffs(k+1,:),n,1);
    b=out.r*breg;
    int=robpcareg.coeffs(k+1,:)-out.centerz(1:p1)*out.r*breg;

    if rmsep==1
        rmse=rrmse(x,y,h,kmax,'RSIMPLS',0,k);
        out.rmsep=rmse.rmsep;
    end

    % testing several output parameters
    if ~isnan(rmsep_value)
        out.rmsep = rmsep_value;
        options.rmsep = 1;
    end

    if ~isnan(rsquared_value)
        out.rcs = rsquared_value;
    end
    if ~isnan(rss_value)
        out.rcs = [out.rcs;sqrt(rss_value)];
    end
    if ~isnan(rmsecv_value)
        gammahalf = 0.5*sqrt(rss_value) + 0.5*rmsecv_value;
        out.rcs = [out.rcs;gammahalf;rmsecv_value];
        options.rmsecv = 1;
    end
    if any(isnan(rsquared_value)) && any(isnan(rss_value)) && any(isnan(rmsecv_value))
        out.rcs = 0;
    end

   %The output:
   out.T=out.t;
   out.weights.p=out.p;
   out.weights.r=out.r;
   out.kr=kr;
   out.h=h;
   out.alpha=alfa;
   out.slope=b;
   out.int=int;
   out.yhat=x*b+repmat(int,n,1);
   out.x=x;
   out.y=y;
   out.res=y-out.yhat;
   out.class='RSIMPLS';
   out.k=k;
   out.cov=robpcareg.cov;

   if ~st
       %calculation of robust distances
       %Score distance
       out.Tcov=robpcareg.sigma(1:k,1:k);
       out.Tcenter=robpcareg.center(1:k);
       out.sd=sqrt(mahalanobis(out.t,out.Tcenter,'cov',out.Tcov))';
       out.cutoff.sd=sqrt(chi2inv(0.975,k));
       %robust residual distance
       if q1==1
           out.resd=out.res/sqrt(out.cov);
       else
           out.resd=sqrt(mahalanobis(out.res,zeros(1,q1),'cov',out.cov))';
       end
       %robust orthogonal distances
       xtilde=out.t*out.p';
       Cdiff=out.xcentr-xtilde;
       for i=1:n
           out.od(i,1)=norm(Cdiff(i,:));
       end
       r=rank(x);
       if k~=r
           [m,s]=unimcd(out.od.^(2/3),out.h);
           out.cutoff.od = sqrt(norminv(0.975,m,s)^3);
       else
           out.cutoff.od=0;
       end
       out.cutoff.resd=sqrt(chi2inv(0.975,q1));

       %Computing flags
       out.flag.od=out.od<=out.cutoff.od;
       out.flag.resd=abs(out.resd)<=out.cutoff.resd;
       out.flag.all=(out.flag.od & out.flag.resd);

       %Multivariate Rsquared
       Yw=y(out.flag.all==1,:);
       cYw=mcenter(Yw);
       res=out.res(out.flag.all==1,:);
       out.rsquared=1-(det(res'*res)/det(cYw'*cYw));

       %Assigning output
       if options.rmsep~=1 && options.rmsecv~=1
           out.rmsep=0;
       end

       if classic
           resultclassic=csimpls(x,y,'k',k,'plots',0);
       else
           resultclassic=0;
       end
       result=struct('slope',{out.slope}, 'int',{out.int},'fitted',{out.yhat},'res',{out.res}, 'cov',{out.cov},...
           'T',{out.T}, 'weights', {out.weights},'Tcenter',{out.Tcenter},'Tcov',{out.Tcov},'rsquared',{out.rsquared},'rcs',{out.rcs},'rmsep',{out.rmsep},...
           'k',{out.k},'alpha',{out.alpha},'h',{out.h},'sd', {out.sd},'od',{out.od},...
           'resd',{out.resd},'cutoff',{out.cutoff},'flag',{out.flag},'class',{out.class},'classic',{resultclassic},'robpca',{out.robpca});
       if result.rcs==0
           result=rmfield(result,'rcs');
       end
       if result.rmsep==0
           result=rmfield(result,'rmsep');
       end
   else
       result=out;
   end
end

% Plots
try
    if plots && options.classic
        makeplot(result,'classic',1,'labsd',labsd,'labod',labod,'labresd',labresd)
    elseif plots
        makeplot(result,'labsd',labsd,'labod',labod,'labresd',labresd)
    end
catch %output must be given even if plots are interrupted
end

author: marianux

repo: ecg-kit

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


text:
function outW = extract(w, method, startV, endV)
   %EXTRACT creates a waveform with a subset of another's data.
   %   waveform = extract(waveform, 'TIME', startTime, endTime)
   %       returns a waveform with the subset of data from startTime to
   %       endTime.  Both times are matlab formatted (string or datenum)
   %
   %   waveform = extract(waveform, 'INDEX', startIndex, endIndex)
   %       returns a waveform with the subset of data from StartIndex to
   %       EndIndex.  this is roughly equivelent to grabbing the waveform's
   %       data into an array, as in D = get(W,'data'), then returning a
   %       waveform with the subset of data,
   %       ie. waveform = set(waveform,'data', D(startIndex:endIndex));
   %
   %   waveform = extract(waveform, 'INDEX&DURATION', startIndex, duration)
   %       a hybrid method that starts from data index startIndex, and then
   %       returns a specified length of data as indicated by duration.
   %       Duration is a matlab formatted time (string or datenum).
   %
   %   waveform = extract(waveform, 'TIME&SAMPLES', startTime, Samples)
   %       a hybrid method that starts from data index startIndex, and then
   %       returns a specified length of data as indicated by duration.
   %       Duration is a matlab formatted time (string or datenum).
   %
   %   Input Arguments:
   %       WAVEFORM: waveform object        N-DIMENSIONAL
   %       METHOD: 'TIME', or 'INDEX', or 'INDEX&DURATION'
   %           TIME: starttime and endtime are absolute times
   %                   (include the date)
   %           INDEX: startt and endt are the offset (index) within the data
   %           INDEX&DURATION: first value is an offset (index), the next says
   %                           how much data to retrieve...
   %           TIME&SAMPLES: grab first value at time startTime, and grab
   %                         Samplength data points
   %       STARTTIME:  Start time (matlab or text format)
   %       ENDTIME:    End time (matlab or text format)
   %       STARTINDEX: position within data array to begin extraction
   %       ENDINDEX:   final grabbed position within data array
   %       DURATION:   matlab format time indicating duration of data to grab
   %       SAMPLES: the number of data points to grab.
   %
   %   the output waveform will have the new, appropriate start time.
   %   if the times are outside the range of the waveform object, then the
   %   output waveform will contain only the portion of the data that is
   %   appropriate.
   %
   %   *MULTIPLE EXTRACTIONS* can be received if the time values are vectors.
   %   Both starttime/startindex and endtime/endindex/endduration/samples must
   %   have the same number of elements.  In this case the output waveforms
   %   will be reshaped with each waveform represented by row, and each
   %   extracted time represented by column.  that is...
   %
   %  The output of this function, for multiple waveforms and times will be:
   %         t1   t2  t3 ... tn
   %    -----------------------
   %    w1 |
   %    w2 |
   %    w3 |
   %     . |
   %     . |
   %    wn |
   %
   %
   %%   examples:
   %       % say that Win is a waveform that starts 1/5/2007 04:00, and
   %       % contains 1 hour of data at 100 Hz (360000 samples)
   %
   %       % grab samples between 4:15 and 4:20
   %       Wout = extract(Win, 'TIME', '1/5/2007 4:15:00','1/5/2007 4:20:00');
   %
   %       % grab 3 minutes, starting at the 10000th sample
   %       Wout = extract(Win, 'INDEX&DURATION', 10000 , '0/0/0 00:03:00');
   %
   %
   %%     example of multiple extract:
   %  % declare the times we're interested in
   %         firstsnippet = datenum('6/20/2003 00:00:00');
   %         lastsnippet = datenum('6/20/2003 24:00:00');
   %
   %         % divide the day into 1-hour segments.
   %         % note, 25 peices. equivelent to 0:1:24, including both midnights
   %         alltimes = linspace(firstsnippet, lastsnippet, 25);
   %         starttimes = alltimes(1:end-1);
   %         endtimes = alltimes(2:end);
   %
   %         % grab each hour of time, and shove it into wHours
   %         wHours = extract(wDay, 'time',starttimes, endtimes);
   %
   %         scaleFactor = 4 * std(double(wDay));
   %         wHours = wHours ./ scaleFactor;
   % %
   %          for n = 1:length(wHours)
   %            wHours(n) = -wHours(n) + n; %add offset for plotting
   %          end
   %          plot(wHours,'xunit','m','b'); %plot it in blue with at nm scaling
   %          axis ([0 60 1 25])
   %          set(gca,'ytick',[0:2:24],'xgrid', 'on','ydir','reverse');
   %          ylabel('Hour');
   %
   %   See also WAVEFORM/SET -- Sample_Length

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   %% Set up condition variables, and ensure validity of input
   MULTIPLE_WAVES = ~isscalar(w);

   %if either of our times are strings, it's 'cause they're actually dates
   if ischar(startV)
      startV = datenum(startV);
   end
   if ischar(endV)
      endV = datenum(endV);
   end

   if numel(startV) ~= numel(endV)
      error('Waveform:extract:indexMismatch',...
         'Number of start times (or indexes) must equal number of end times')
   end

   % are we getting a series of extractions from each waveform?
   MULTIPLE_EXTRACTION = numel(endV) > 1;

   if MULTIPLE_WAVES && MULTIPLE_EXTRACTION
      w = w(:);
   end

   %%
   if numel(w)==0 || numel(startV) ==0
      warning('Waveform:extract:emptyWaveform','no waveforms to extract');
      return
   end

   % changed by glenn thompson 20161006 as it was losing calibration
   % information
   %outW(numel(w),numel(startV)) = waveform; % removed
   outW = w; % added

   for m = 1: numel(startV) %loop through the number of extractions
      for n=1:numel(w); %loop through the waveforms
         inW = w(n);
         myData = inW.data;

         switch lower(method)
            case 'time'

               % startV and endV are both matlab formated dates
               %sampleTimes = get(inW,'timevector');

               %   ensure the format of our times
               if startV(m) > endV(m)
                  warning('Waveform:extract:reversedValues',...
                     'Start time prior to end time.  Flipping.');
                  [startV(m), endV(m)] = swap(startV(m), endV(m));
               end

               %if requested data is outside the existing waveform, change the
               %start time, and clear out the data.
               startsAfterWave = startV(m) > get(inW,'end') ;
               endsBeforeWave = endV(m) < get(inW,'start');
               if startsAfterWave || endsBeforeWave
                  myStart = startV(m);
                  myData = [];
               else
                  %some aspect of this data must be represented by the waveform
                  [myStartI myStartTime] = time2offset(inW,startV(m));
                  [myEndI] = time2offset(inW,endV(m)) - 1;
                  if isempty(myStartTime)
                     %waveform starts sometime after requested start
                     myStartTime = get(inW,'start');
                     myStartI = 1;
                  end

                  if myEndI > numel(myData)
                     myEndI = numel(myData);
                  end
                  myData = myData(myStartI:myEndI);
                  myStart = myStartTime;
               end
            case 'index'
               %startV and endV are both indexes into the data

               if startV(m) > numel(myData)
                  warning('Waveform:extract:noDataFound',...
                     'no data after start index');
                  return
               end;
               if endV(m) > numel(myData)
                  endV(m) = length(myData);
                  warning('Waveform:extract:truncatingData',...
                     'end index too long, truncating to match data');
               end

               if startV(m) > endV(m)
                  warning('Waveform:extract:reversedValues',...
                     'Start time prior to end time.  Flipping.');
                  [startV(m), endV(m)] = swap(startV(m), endV(m));
               end

               myData = myData(startV(m):endV(m));
               sampTimes = get(inW,'timevector'); % grab individual sample times
               myStart = sampTimes(startV(m));

            case 'index&duration'
               % startV is an index into the data, endV is a matlab date
               myData = myData(startV(m):end); %grab the data starting at our index

               sampTimes = get(inW,'timevector'); % grab individual sample times
               sampTimes = sampTimes(startV(m):end); % truncate to match data

               myStart = sampTimes(1); %grab our starting date before hacking it

               sampTimes = sampTimes - sampTimes(1); %set first time to zero
               count = sum(sampTimes <= endV(m)) -1;
               myData = myData(1:count);

            case 'time&samples'
               % startV is a matlab date, while endV is an index into the data
               sampTimes = get(inW,'timevector'); % grab individual sample times

               index_to_times = sampTimes >= startV(m); %mask of valid times
               goodTimes = find(index_to_times,endV(m));%first howevermany of these good times

               myData = myData(goodTimes); % keep matching samples

               try
                  myStart = sampTimes(goodTimes(1)); %first sample time is new waveform start
               catch
                  warning('Waveform:extract:NoDataFound',...
                     'no data');
                  myStart = startV(1);
               end

            otherwise
               error('Waveform:extract:unknownMethod','unknown method: %s', method);
         end

         if MULTIPLE_EXTRACTION
            outW(n,m) = set(inW,'start',myStart, 'data', myData);
         else
            outW(n) = set(inW,'start',myStart, 'data', myData);
         end
      end % n-loop (looping through waveforms)
   end % m-loop (looping through extractions)
end

function  [B,A] = swap(A,B)
   %do nothing, just flip inputs & outputs
end

author: geoscience-community-codes

repo: GISMO


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MRiLab auto generated file: DO NOT EDIT!     %
% Generated by MRiLab "DoWriteXML2m" Generator %
% MRiLab Version 1.3                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rfAmp,rfPhase,rfFreq,rfCoil,GzAmp,GyAmp,GxAmp,ADC,Ext,uts,ts,flags]=PSD_FIESTA3D
global VCtl
global VVar
CV1=3e-3;
CV10=0;
CV11=0;
CV12=0;
CV13=0;
CV14=0;
CV2=4e-3;
CV3=2*pi/2;
CV4=2;
CV5=0;
CV6=0;
CV7=0;
CV8=0;
CV9=0;
rfAmpAll=[];
rfPhaseAll=[];
rfFreqAll=[];
rfCoilAll=[];
GzAmpAll=[];
GyAmpAll=[];
GxAmpAll=[];
ADCAll=[];
ExtAll=[];
rfTimeAll=[];
GzTimeAll=[];
GyTimeAll=[];
GxTimeAll=[];
ADCTimeAll=[];
ExtTimeAll=[];
SEtAll=[];
uts=[];
ts=[];
flags=[];
if VCtl.PlotSeq == 1
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
ts = [ts tS tE];
end
ts = [0 max(ts)-min(ts)];
return;
end
%==============Pulses 1==============
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
rfTime=[];
GzTime=[];
GyTime=[];
GxTime=[];
ADCTime=[];
ExtTime=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if isempty(tS) | isempty(tE) | (tS>=tE)
error('SE setting is incorrect for Pulses 1!');
end
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{1};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{3};
p.CoilID=1;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='alpha pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=rem(VVar.TRCount-1,CV4)*CV3;
p.tEnd=1.5e-3;
p.tStart=1.2e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp1,rfPhase1,rfFreq1,rfCoil1,rfTime1]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil1(1)
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
else
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{3};
p.CoilID=2;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='alpha pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=rem(VVar.TRCount-1,CV4)*CV3;
p.tEnd=1.5e-3;
p.tStart=1.2e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp2,rfPhase2,rfFreq2,rfCoil2,rfTime2]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil2(1)
rfAmp=[rfAmp rfAmp2];
rfPhase=[rfPhase rfPhase2];
rfFreq=[rfFreq rfFreq2];
rfCoil=[rfCoil rfCoil2];
rfTime=[rfTime rfTime2];
end
else
rfAmp=[rfAmp rfAmp2];
rfPhase=[rfPhase rfPhase2];
rfFreq=[rfFreq rfFreq2];
rfCoil=[rfCoil rfCoil2];
rfTime=[rfTime rfTime2];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{3};
p.CoilID=3;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='alpha pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=rem(VVar.TRCount-1,CV4)*CV3;
p.tEnd=1.5e-3;
p.tStart=1.2e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp3,rfPhase3,rfFreq3,rfCoil3,rfTime3]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil3(1)
rfAmp=[rfAmp rfAmp3];
rfPhase=[rfPhase rfPhase3];
rfFreq=[rfFreq rfFreq3];
rfCoil=[rfCoil rfCoil3];
rfTime=[rfTime rfTime3];
end
else
rfAmp=[rfAmp rfAmp3];
rfPhase=[rfPhase rfPhase3];
rfFreq=[rfFreq rfFreq3];
rfCoil=[rfCoil rfCoil3];
rfTime=[rfTime rfTime3];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{3};
p.CoilID=4;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='alpha pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=rem(VVar.TRCount-1,CV4)*CV3;
p.tEnd=1.5e-3;
p.tStart=1.2e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp4,rfPhase4,rfFreq4,rfCoil4,rfTime4]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil4(1)
rfAmp=[rfAmp rfAmp4];
rfPhase=[rfPhase rfPhase4];
rfFreq=[rfFreq rfFreq4];
rfCoil=[rfCoil rfCoil4];
rfTime=[rfTime rfTime4];
end
else
rfAmp=[rfAmp rfAmp4];
rfPhase=[rfPhase rfPhase4];
rfFreq=[rfFreq rfFreq4];
rfCoil=[rfCoil rfCoil4];
rfTime=[rfTime rfTime4];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{3};
p.CoilID=5;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='alpha pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=rem(VVar.TRCount-1,CV4)*CV3;
p.tEnd=1.5e-3;
p.tStart=1.2e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp5,rfPhase5,rfFreq5,rfCoil5,rfTime5]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil5(1)
rfAmp=[rfAmp rfAmp5];
rfPhase=[rfPhase rfPhase5];
rfFreq=[rfFreq rfFreq5];
rfCoil=[rfCoil rfCoil5];
rfTime=[rfTime rfTime5];
end
else
rfAmp=[rfAmp rfAmp5];
rfPhase=[rfPhase rfPhase5];
rfFreq=[rfFreq rfFreq5];
rfCoil=[rfCoil rfCoil5];
rfTime=[rfTime rfTime5];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{3};
p.CoilID=6;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='alpha pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=rem(VVar.TRCount-1,CV4)*CV3;
p.tEnd=1.5e-3;
p.tStart=1.2e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp6,rfPhase6,rfFreq6,rfCoil6,rfTime6]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil6(1)
rfAmp=[rfAmp rfAmp6];
rfPhase=[rfPhase rfPhase6];
rfFreq=[rfFreq rfFreq6];
rfCoil=[rfCoil rfCoil6];
rfTime=[rfTime rfTime6];
end
else
rfAmp=[rfAmp rfAmp6];
rfPhase=[rfPhase rfPhase6];
rfFreq=[rfFreq rfFreq6];
rfCoil=[rfCoil rfCoil6];
rfTime=[rfTime rfTime6];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{3};
p.CoilID=7;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='alpha pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=rem(VVar.TRCount-1,CV4)*CV3;
p.tEnd=1.5e-3;
p.tStart=1.2e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp7,rfPhase7,rfFreq7,rfCoil7,rfTime7]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil7(1)
rfAmp=[rfAmp rfAmp7];
rfPhase=[rfPhase rfPhase7];
rfFreq=[rfFreq rfFreq7];
rfCoil=[rfCoil rfCoil7];
rfTime=[rfTime rfTime7];
end
else
rfAmp=[rfAmp rfAmp7];
rfPhase=[rfPhase rfPhase7];
rfFreq=[rfFreq rfFreq7];
rfCoil=[rfCoil rfCoil7];
rfTime=[rfTime rfTime7];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{3};
p.CoilID=8;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='alpha pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=rem(VVar.TRCount-1,CV4)*CV3;
p.tEnd=1.5e-3;
p.tStart=1.2e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp8,rfPhase8,rfFreq8,rfCoil8,rfTime8]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil8(1)
rfAmp=[rfAmp rfAmp8];
rfPhase=[rfPhase rfPhase8];
rfFreq=[rfFreq rfFreq8];
rfCoil=[rfCoil rfCoil8];
rfTime=[rfTime rfTime8];
end
else
rfAmp=[rfAmp rfAmp8];
rfPhase=[rfPhase rfPhase8];
rfFreq=[rfFreq rfFreq8];
rfCoil=[rfCoil rfCoil8];
rfTime=[rfTime rfTime8];
end
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gz1Sign=1;
p.Gz2Sign=-1;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=VCtl.TE-CV1;
p.t1Start=VCtl.TE-CV2;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GzAmp1,GzTime1]=GzCartesian(p);
GzAmp=[GzAmp GzAmp1];
GzTime=[GzTime GzTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gy1Sign=1;
p.Gy2Sign=-1;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=VCtl.TE-CV1;
p.t1Start=VCtl.TE-CV2;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GyAmp1,GyTime1]=GyCartesian(p);
GyAmp=[GyAmp GyAmp1];
GyTime=[GyTime GyTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gx1Sign=-1;
p.Gx2Sign=1;
p.Gx3Sign=-1;
p.Notes='cartesian frequency';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1Start=VCtl.TE-CV2;
p.t2Middle=VCtl.TE;
p.t3Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GxAmp1,GxTime1]=GxCartesian(p);
GxAmp=[GxAmp GxAmp1];
GxTime=[GxTime GxTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Notes='cartesian readout';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tMiddle=VCtl.TE;
if strcmp(p.Switch,'on')
[ADC1,ADCTime1]=ADCCartesian(p);
ADC=[ADC ADC1];
ADCTime=[ADCTime ADCTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=5;
p.Notes='calculate remaining scan time';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=0;
if strcmp(p.Switch,'on')
[Ext1,ExtTime1]=ExtBit(p);
Ext=[Ext Ext1];
ExtTime=[ExtTime ExtTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=7;
p.Notes='demodulate signal phase referring to rf phase at ADC';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=1.35e-3;
if strcmp(p.Switch,'on')
[Ext2,ExtTime2]=ExtBit(p);
Ext=[Ext Ext2];
ExtTime=[ExtTime ExtTime2];
end
p=[];
%--------------------
SEt=[tS tE];
rfAmp(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfPhase(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfFreq(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfCoil(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzAmp(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyAmp(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxAmp(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADC(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
Ext(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfTime(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzTime(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyTime(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxTime(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADCTime(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
ExtTime(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfAmp(abs(rfAmp)<eps) = 0;
rfTime = rfTime + SEt(1);
GzTime = GzTime + SEt(1);
GyTime = GyTime + SEt(1);
GxTime = GxTime + SEt(1);
ADCTime = ADCTime + SEt(1);
ExtTime = ExtTime + SEt(1);
rfAmpAll=[rfAmpAll rfAmp];
rfPhaseAll=[rfPhaseAll rfPhase];
rfFreqAll=[rfFreqAll rfFreq];
rfCoilAll=[rfCoilAll rfCoil];
GzAmpAll=[GzAmpAll GzAmp];
GyAmpAll=[GyAmpAll GyAmp];
GxAmpAll=[GxAmpAll GxAmp];
ADCAll=[ADCAll ADC];
ExtAll=[ExtAll Ext];
rfTimeAll=[rfTimeAll rfTime];
GzTimeAll=[GzTimeAll GzTime];
GyTimeAll=[GyTimeAll GyTime];
GxTimeAll=[GxTimeAll GxTime];
ADCTimeAll=[ADCTimeAll ADCTime];
ExtTimeAll=[ExtTimeAll ExtTime];
SEtAll=[SEtAll SEt];
end
%====================================
if isempty(rfTimeAll)
error('rf sequence line can not be empty! Master Tx coil element must be used.');
end
if isempty(GzTimeAll)
error('GzSS sequence line can not be empty!');
end
if isempty(GyTimeAll)
error('GyPE sequence line can not be empty!');
end
if isempty(GxTimeAll)
error('GxR sequence line can not be empty!');
end
if isempty(ADCTimeAll)
error('ADC sequence line can not be empty!');
end
if isempty(ExtTimeAll)
error('Ext sequence line can not be empty!');
end
SEflag=repmat([0 0 0 0 0 0]',[1 2]);
rfflag=repmat([1 0 0 0 0 0]',[1 max(size(rfTimeAll))]);
Gzflag=repmat([0 1 0 0 0 0]',[1 max(size(GzTimeAll))]);
Gyflag=repmat([0 0 1 0 0 0]',[1 max(size(GyTimeAll))]);
Gxflag=repmat([0 0 0 1 0 0]',[1 max(size(GxTimeAll))]);
ADCflag=repmat([0 0 0 0 1 0]',[1 max(size(ADCTimeAll))]);
Extflag=repmat([0 0 0 0 0 1]',[1 max(size(ExtTimeAll))]);
ts=[[min(SEtAll) max(SEtAll)] rfTimeAll GzTimeAll GyTimeAll GxTimeAll ADCTimeAll ExtTimeAll]-min(SEtAll);
flags=[SEflag rfflag Gzflag Gyflag Gxflag ADCflag Extflag];
[ts,ind]=sort(ts);
uts=unique(ts);
flags=flags(:,ind);
[rfTime,ind]=sort(rfTimeAll-min(SEtAll));
rfAmp=rfAmpAll(:,ind);
rfPhase=rfPhaseAll(:,ind);
rfFreq=rfFreqAll(:,ind);
rfCoil=rfCoilAll(:,ind);
[GzTime,ind]=sort(GzTimeAll-min(SEtAll));
GzAmp=GzAmpAll(:,ind);
[GyTime,ind]=sort(GyTimeAll-min(SEtAll));
GyAmp=GyAmpAll(:,ind);
[GxTime,ind]=sort(GxTimeAll-min(SEtAll));
GxAmp=GxAmpAll(:,ind);
[ADCTime,ind]=sort(ADCTimeAll-min(SEtAll));
ADC=ADCAll(:,ind);
[ExtTime,ind]=sort(ExtTimeAll-min(SEtAll));
Ext=ExtAll(:,ind);
rfAmp(1) = 0;
rfPhase(1) = 0;
rfFreq(1) = 0;
GzAmp(1) = 0;
GyAmp(1) = 0;
GxAmp(1) = 0;
ADC(1) = 0;
Ext(1) = 0;
rfAmp(end) = 0;
rfPhase(end) = 0;
rfFreq(end) = 0;
GzAmp(end) = 0;
GyAmp(end) = 0;
GxAmp(end) = 0;
ADC(end) = 0;
Ext(end) = 0;
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [xx, yy] = secdraw(start_angle, sector_angle, radius)

%  usage: [xx, yy]  = secdraw(start_angle, sector_angle, radius, clk_wise)
%
%  Input:
%            start_angle  - starting angle in degrees, [0 is default]
%                           NOTE: this can be positive or negative
%            sector_angle - central angle of sector in degrees
%                           NOTE: if negative, plotting will be clockwise
%                                 direction.
%            radius       - radius of sector (default 1)
%
% Output:
%            returns the corner points (xx, yy) of a patch for a sector,
%            which can plotted using the patch(...) function
%            NOTE:
%              1. if less than 2 or NO output variables are given,
%                 secdraw(...) draws the sector on the current (or a new)
%                 figure
%              2. if one argument is given,
%                 the argument is taken as the central angle
%
% Example 1 > [xx, yy]  = secdraw(45, 90, 10);
%             patch( xx, yy, 'b', 'EdgeColor', 'b', 'EraseMode', 'none') ;
%             set( gca, 'DataAspectRatio', [1 1 1] );
%             - this will draw a sector from 45 to 90 degrees with radius=1
%
% Example 2 > secdraw( 60 )
%             - this will draw a sector from 0 to 60 degrees counter
%               clockwise
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% author: Laine Berhane Kahsay
%
% date  : 19-Feb-2004
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 1 | isempty( nargin )
    %disp '';
    disp '-- error in usage --------------------------------';
    disp '      Not enough input arguments.';
    help secdraw;
    return;
elseif nargin == 1
    sector_angle   = start_angle;
    start_angle = 0;
    radius      = 1;
else nargin == 2
    radius  = 1;
end

theta0 = pi*0/180;              % offset angle in radians
theta1 = pi*start_angle/180;    % starting angle in radians
theta  = pi*sector_angle/180;   % central angle in radians

angle_ratio = theta/(2*pi);     % ratio of the sector to the complete circle
rho = abs(radius);              % abs(...) can be removed
MaxPts = 100;                   % maximum number of points in a whole circle.
                                % if set to small # (e.g 10) the resolution
                                % of the curve will be poor.
                                % generally values greater than 50
                                % will give (very) good resolution

n = abs( ceil( MaxPts*angle_ratio ) );
r = [ 0; rho*ones(n+1,1); 0 ];
theta = theta0 + theta1 + [ 0; angle_ratio*(0:n)'/n; 0 ]*2*pi;

% output
[xx,yy] = pol2cart(theta,r);

% plot if not enough output variable are given
if nargout < 2,
    hh=patch(xx, yy, 'b','EdgeColor','b','EraseMode','none');
end

author: Sable

repo: mcbench-benchmarks


text:
function a = lehmer ( m, n )

%*****************************************************************************80
%
%% LEHMER returns the LEHMER matrix.
%
%  Discussion:
%
%    This matrix is also known as the "Westlake" matrix.
%
%  Formula:
%
%    A(I,J) = min ( I, J ) / max ( I, J )
%
%  Example:
%
%    N = 5
%
%    1/1  1/2  1/3  1/4  1/5
%    1/2  2/2  2/3  2/4  2/5
%    1/3  2/3  3/3  3/4  3/5
%    1/4  2/4  3/4  4/4  4/5
%    1/5  2/5  3/5  4/5  5/5
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is positive definite.
%
%    A is totally nonnegative.
%
%    The inverse of A is tridiagonal.
%
%    The condition number of A lies between N and 4*N*N.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Morris Newman, John Todd,
%    The evaluation of matrix inversion programs,
%    Journal of the Society for Industrial and Applied Mathematics,
%    Volume 6, Number 4, 1958, pages 466-476.
%
%    Solutions to problem E710, proposed by DH Lehmer: The inverse of
%    a matrix.
%    American Mathematical Monthly,
%    Volume 53, Number 9, November 1946, pages 534-535.
%
%    John Todd,
%    Basic Numerical Mathematics, Volume 2: Numerical Algebra,
%    Academic Press, 1977, page 154.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  for i = 1 : m
    for j = 1 : n
      a(i,j) = ( min ( i, j ) ) / ( max ( i, j ) );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef DC2_DTLZ1 < PROBLEM
% <multi/many> <real> <large/none> <constrained> <expensive/none>
% DTLZ1 with constrains in decision space

%------------------------------- Reference --------------------------------
% K. Li, R. Chen, G. Fu, and X. Yao, Two-archive evolutionary algorithm for
% constrained multi-objective optimization, IEEE Transactions on
% Evolutionary Computation, 2018, 23(2): 303-315.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = obj.M + 4; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            g      = 100*(obj.D-obj.M+1+sum((PopDec(:,obj.M:end)-0.5).^2-cos(20.*pi.*(PopDec(:,obj.M:end)-0.5)),2));
            PopObj = 0.5*repmat(1+g,1,obj.M).*fliplr(cumprod([ones(size(PopDec,1),1),PopDec(:,1:obj.M-1)],2)).*[ones(size(PopDec,1),1),1-PopDec(:,obj.M-1:-1:1)];
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            g = 100*(obj.D-obj.M+1+sum((PopDec(:,obj.M:end)-0.5).^2-cos(20.*pi.*(PopDec(:,obj.M:end)-0.5)),2));
            PopCon(:,1) = 0.5 - cos(3*pi*g);
            PopCon(:,2) = 0.5 - exp(-g);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,obj.M)/2;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a = linspace(0,1,10)';
                R = {a*a'/2,a*(1-a')/2,(1-a)*ones(size(a'))/2};
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function env = spx_get_env();
    global spx_env;
    if isempty(spx_env)
        filepath = which(mfilename);
        root = fileparts(filepath);
        spx_settings_cache = fullfile(root, 'spx_local_settings.mat');
        data = load(spx_settings_cache);
        spx_env = data.globals;
    end
    env = spx_env;
end

author: indigits

repo: sparse-plex


text:
function [hCylinder,hPlate1,hPlate2] = Cylinder( orgin,r,h,dir,n,closed )
% This Function plots a cylinder at specified orgin, with specified radius, height, direction(axis), no of points along the circumference
%
%   Typical Call Cylinder( orgin,r,h,dir,n ), Note: there is a matlab function "cylinder"
%
%   orgin: vector of order 3 x 1 specifies orgin
%   h    : Height of the Cylinder
%   dir  : String to specify axis of extrution 'x' or 'y' or 'z' (Only along these axes...!)
%   n    : no of points along the circumference
%

t=linspace(0,2*pi,n)';

x1=r*cos(t);
x2=r*sin(t);
h1=orgin(3);
h2=h1+h;

if dir=='y'
    xx1=[[x1;x1(1)] [x1;x1(1)]]+orgin(1);
    xx2=[repmat(h1,length(x1)+1,1) repmat(h2,length(x1)+1,1)]+orgin(2);
    xx3=[[x2;x2(1)] [x2;x2(1)]]+orgin(3);
elseif dir =='x'
    xx1=[repmat(h1,length(x1)+1,1) repmat(h2,length(x1)+1,1)]+orgin(1);
    xx2=[[x1;x1(1)] [x1;x1(1)]]+orgin(2);
    xx3=[[x2;x2(1)] [x2;x2(1)]]+orgin(3);
else
    xx1=[[x1;x1(1)] [x1;x1(1)]]+orgin(1);
    xx2=[[x2;x2(1)] [x2;x2(1)]]+orgin(2);
    xx3=[repmat(h1,length(x1)+1,1) repmat(h2,length(x1)+1,1)]+orgin(3);
end
hCylinder=surf(xx1,xx2,xx3,repmat(3,size(xx1)));

if strcmp(closed,'closed')==1
    hold on
    hPlate1=fill3(xx1(:,1),xx2(:,1),xx3(:,1),'g');
    hPlate2=fill3(xx1(:,2),xx2(:,2),xx3(:,2),'g');
end


text:
function ned = ecef2ned(ecef, llh_org)
% ecef2ned: converts from ECEF coordinates to NED coordinates.
%
% INPUTS
%   ecef: Nx3 ECEF coordinates [X Y Z] (m, m, m).
%   llh_org: 1x3 system origin [lat, lon, h] (rad, rad, m).
%
% OUTPUTS
%   ned: Nx3 NED coordinates [X Y Z] (m, m, m).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Reference:
%	R. Gonzalez, J. Giribet, and H. Patiño. NaveGo: a
% simulation framework for low-cost integrated navigation systems,
% Journal of Control Engineering and Applied Informatics}, vol. 17,
% issue 2, pp. 110-120, 2015. Inverse process of Eq. 15.
%
% Version: 003
% Date:    2019/01/16
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

lat = llh_org(1);
lon = llh_org(2);

ecef_org = llh2ecef(llh_org)';

slat = sin(lat);
clat = cos(lat);
slon = sin(lon);
clon = cos(lon);

R = [  -slat*clon  -slat*slon   clat; ...
       -slon          clon         0; ...
       -clat*clon  -clat*slon  -slat];

[MAX, N] = size(ecef);
ned = zeros(MAX, N);

for i=1:MAX

    ned_t = R * (ecef(i, :)' - ecef_org) ;
    ned (i, :) = ned_t';
end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


function [ytrend,ycycle]=one_sided_hpfilter_serial(y,lambda,discard)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% I am grateful to Martin Kliem for pointing out this alternative to the
% Kalman variant to me and for sharing his code with me.
%
% Copyright: Alexander Meyer-Gohde
%
% You are free to use/modify/redistribute this program so long as original
% authorship credit is given and you in no way impinge on its free
% distribution
%
%
% This program executes a one-sided HP filter by running the
% standard two-sided HP filter successively with each new observation
%
%  Input: y - a Txn data matrix, where T is the number of observations
%            on n variables (i.e.,data is assumed to be in column format).
%         lambda - a smoothing scalar. If not entered, default is 1600.
%         discard  -  a scalar. The first discard periods will be
%            discarded resulting in output matrices of size
%            (T-discard)xn. Optional: if not entered, default value is 0.
%
% Output: ytrend - a (T-discard)xn matrix of extracted trends for
%                  each of the n variables.
%         ycycle a  (T-discard)xn matrix of deviations from the trends
%                     for  each of the n variables. Optional.
%
%       Usage example:
%       [ytrend]=one_sided_hp_filter_kalman(y)
%                   will yield the Txn matrix of trends using the data in
%                   y with lambda set to 1600
%
%
% This one-side HP filter finds a series {ytrend_t}_{t=1}^T for
% each n by calculating for each t the standard HP filtered trend using
% all data upto that t and setting ytrend_t equal to the resulting trend
% value for period t.
%
% See: Mehra, Y.P. (2004). "The Output Gap, Expected Future Inflation and
%                        Inflation Dynamics: Another Look,"
%           The B.E. Journal of Macroeconomics, Berkeley Electronic Press.
%
% The standard HP filter finds a series {ytrend_t}_{t=1}^T that solves
% the following minimization probelm
%
% min sum_{t=1}^T(y_t-ytrend_t)
% +lambda*sum_{t=2}^{T-1}[(ytrend_{t+1}-ytrend_{t})-(ytrend_{t}-ytrend_{t-1})]
%
% Rearanging the first order conditions yields:  A*ytrend=y   where
% A=[   1+lambda    -2*lambda       lambda        0               0           0 ...
%     [    -2*lambda   1+5*lambda   -4*lambda       lambda       0           0 ...
%     [    lambda       -4*lambda        1+6*lambda   -4*lambda  lambda  0 ...
%     [                                         ...
%     [    0      ...   0   lambda       -4*lambda        1+6*lambda   -4*lambda  lambda 0  ... 0]
%     [                                         ...      ]
%     [  0   ...    lambda       -4*lambda        1+6*lambda  -4*lambda  lambda]
%     [  0   ...           0     lambda  -4*lambda       1+5*lambda      -2*lambda ]
%     [  0   ...           0       0        lambda              -2*lambda    1+lambda]
%
%
%   Hodrick, R.J. and E.C.Prescott (1997), "Postwar U.S. Business Cycles:
%   An Empirical Investigation." Jounal of Money, Credit and Banking.
%   29(1), Feb. pp. 1--16.
%
%  Note that filter will be calculated over and over again with an
%  expanding matrix A. This program uses sparse matrices and exploits the
%  patter in A as t progresses.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 2,  lambda = 1600; end%If the user does not provide lambda, set it to 1600
[T,n] = size (y); % Calculate the number of periods and the number of variables in the series

%%Preliminary calculations
x1=[1+lambda, -2*lambda, lambda]; %The non-zero elements of the first row of A
x2=[ -2*lambda, 1+5*lambda, -4*lambda, lambda];%The non-zero elements of the second row of A
x3=[lambda, -4*lambda, 1+6*lambda, -4*lambda, lambda];%The non-zero elements of thej'th row of A, 2<j<T-2
x2rev=x2(end:-1:1);%The non-zero elements of the second-to-last row of A (just x2 in reverse)
x1rev=x1(end:-1:1); %The non-zero elements of the last row of A (just x1 in reverse)

%The minimization problem for the HP filter is either trivially given by
%the following or is not defined when there are less than three
%observations
for t=1:2
    ytrend(t,:)=y(t,:); %Thus, the trend is set to the observation
end
t=t+1;
%For three observations, the minimization problem yields the following
junk=[x1;-2*lambda, 1+4*lambda, -2*lambda;x1rev]\y(1:t,:);%Build the matrix A and solve the system A*y_standard_hp_trend=y
ytrend(t,:)=junk(t,:);%select the t'th element of the resulting trend
%Starting with four observations, the pattern using the preliminary
%calculations begins
t=t+1;
Ibegin=[1; 1; 1; 2; 2; 2; 2];%make a list (a column vector) containing at position i the row of the i'th non-zero element in the first two rows of A
Jbegin=[1;2;3;1;2;3;4];%make a list (a column vector) containing at position i the column of the i'th non-zero element in the first two rows of A
Xbegin=[x1';x2'];%make a list (a column vector) containing at position i the i'th non-zero element in the first two rows of A
Iend=[3; 3; 3; 3; 4; 4; 4];%make a list (a column vector) containing at position i the row of the i'th non-zero element in the last two rows of A
Jend=[1;2;3;4;2;3;4];%make a list (a column vector) containing at position i the column of the i'th non-zero element in the last two rows of A
Xend=[x2rev';x1rev'];%make a list (a column vector) containing at position i the i'th non-zero element in the last two rows of A
junk=sparse([Ibegin; Iend],[Jbegin; Jend],[Xbegin; Xend])\y(1:t,:);
ytrend(t,:)=junk(t,:);
%For the fifth observation on up, the simple pattern can be exploited by
%adding the new row and column indicies of the new non-zero elements to the Ibegin and Jbegin as well as
%the values of the elements themselves. Note that this requires the
%row and column indicies of the non-zeor elements in the last two rows to
%be advanced
while t<T
    t=t+1;
    Ibegin=[Ibegin;t-2;t-2;t-2;t-2;t-2];%Add to the list the rows of the non-zero elements of the t-2'th row
    Jbegin=[Jbegin;t-4;t-3;t-2;t-1;t];%Add to the list the col.umns of the non-zero elements of the t-2'th row
    Xbegin=[Xbegin;x3'];%Add to the list the non-zero elements of the t-2'th row
    Iend=Iend+1;%Advance the rows of the non-zero elements in the last two rows
    Jend=Jend+1;%Advance the columns of the non-zero elements in the last two rows
    junk=sparse([Ibegin; Iend],[Jbegin; Jend],[Xbegin; Xend],t,t,(t-4)*5+14)\y(1:t,:);%Build the matrix A and solve the system A*y_standard_hp_trend=y
    ytrend(t,:)=junk(t,:);%select the t'th element of the resulting trend
end

if nargout==2%Should the user have requested a second output
    ycycle=y-ytrend;%The second output will be the deviations from the HP trend
end

if nargin==3%If the user provided a discard parameter
    ytrend=ytrend(discard+1:end,:);%Remove the first "discard" periods from the trend series
    if nargout==2%Should the user have requested a second output
        ycycle=ycycle(discard+1:end,:);%The second output will be the deviations from the HP trend, likewise with  first "discard" periods removed
    end
end
end


text:
function Ne=vertex_neighbours_double(Fa,Fb,Fc,Vx,Vy,Vz)

F=[Fa Fb Fc];
V=[Vx Vy Vz];

% Neighbourh cell array
Ne=cell(1,size(V,1));

% Loop through all faces
for i=1:length(F)
    % Add the neighbors of each vertice of a face
    % to his neighbors list.
    Ne{F(i,1)}=[Ne{F(i,1)} [F(i,2) F(i,3)]];
    Ne{F(i,2)}=[Ne{F(i,2)} [F(i,3) F(i,1)]];
    Ne{F(i,3)}=[Ne{F(i,3)} [F(i,1) F(i,2)]];
end

% Loop through all neighbor arrays and sort them (Rotation same as faces)
for i=1:size(V,1)

    Pneighf=Ne{i};
    if(isempty(Pneighf))
        Pneig=[];
    else
        start=1;
        for index1=1:2:length(Pneighf)
            found=false;
            for index2=2:2:length(Pneighf),
                if(Pneighf(index1)==Pneighf(index2))
                    found=true; break
                end
            end
            if(~found)
                start=index1; break
            end
        end
        Pneig=[];
        Pneig(1)=Pneighf(start);
        Pneig(2)=Pneighf(start+1);

        % Add the neighbours with respect to original rotation
        for j=2+double(found):(length(Pneighf)/2)
            found = false;
            for index=1:2:length(Pneighf),
                if(Pneighf(index)==Pneig(end))
                    if(sum(Pneig==Pneighf(index+1))==0)
                        found =true;
                        Pneig=[Pneig Pneighf(index+1)];
                    end
                end
            end
            if(~found) % This only happens with weird edge vertices
                for index=1:2:length(Pneighf),
                    if(sum(Pneig==Pneighf(index))==0)
                        Pneig=[Pneig Pneighf(index)];
                        if(sum(Pneig==Pneighf(index+1))==0)
                            Pneig=[Pneig Pneighf(index+1)];
                        end
                    end
                end
            end
        end
        % Add forgotten neigbours
        if(length(Pneig)<length(Pneighf))
            for j=1:length(Pneighf)
                if(sum(Pneig==Pneighf(j))==0)
                    Pneig=[Pneig Pneighf(j)];
                end
            end
        end
    end
    Ne{i}=Pneig;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function save_to_base(overwrite)
%
% save_to_base() copies all variables in the calling function to the base
% workspace.  This makes it possible to examine a function's internal
% variables from the Matlab command prompt after the calling function
% terminates.
%
% If the optional argument "overwrite" is present and has a non-zero value,
% variables in the workspace of the caller will overwrite variables in the
% base workspace having the same name.  Otherwise, preexisting variables in
% the base workspace will not be overwritten.
%
% Phillip M. Feldman

% Version 1.1, 18 Aug. 2009: Added fix to exclude 'ans' variable.
%
% Version 1.0, 2 June 2009: Initial version.
%
% Note: I'd like to acknowledge assistance from Sean Little of the
% MathWorks.

ws_caller= evalin('caller','who()');

if length(nargin) && overwrite
   % Overwrite preexisting variables:
   variables= ws_caller;
else
   % Do not overwrite preexisting variables:
   ws_base= evalin('base','who()');
   variables= setdiff(ws_caller, ws_base);
end

for i= 1 : length(variables)
   if ~strcmp(variables{i},'ans')
      tmp= evalin('caller',variables{i});
      assignin('base',variables{i},tmp);
   end
end


text:
function [trop_delay,var]=trop_mops(time,pos,azel)

k1=77.604;k2=382000;rd=287.054;gm=9.784;g=9.80665;
persistent pos_ zh zw
sinel=sin(azel(2)); h=pos(3);

if isempty(pos_),pos_=zeros(3,1);end
if isempty(zh),zh=0;end
if isempty(zh),zh=0;end

if pos(3)<-100||pos(3)>10000||azel(2)==0
    trop_delay=0;
    var=0;
    return;
end

if zh==0||abs(pos(1)-pos_(1))>1e-7||abs(pos(2)-pos_(2))>1e-7||abs(pos(3)-pos_(3))>1
    met=getmet(pos(1)*180/pi);
    if pos(1)>=0
        tmp=28;
    else
        tmp=211;
    end
    doy=time2doy(time);
    c=cos(2*pi*(doy-tmp)/365.25);
    for i=1:5
        met(i)=met(i)-met(i+5)*c;
    end
    zh=1E-6*k1*rd*met(1)/gm;
    zw=1E-6*k2*rd/(gm*(met(5)+1.0)-met(4)*rd)*met(3)/met(2);
    zh=zh*(1.0-met(4)*h/met(2))^(g/(rd*met(4)));
    zw=zw*(1.0-met(4)*h/met(2))^((met(5)+1.0)*g/(rd*met(4))-1.0);
    for i=1:3
        pos_(i)=pos(1);
    end
end

m=1.001/sqrt(0.002001+sinel*sinel); %mapping function of GCAT model
trop_delay=(zh+zw)*m;
var=0.12*0.12*m*m;

return

repo: GINav


text:
function out = cf(func,varargin)
%CF Convenience wrapper for cellfun with non-uniform output.
% Usage:
%   out = cf(func, C); % equivalent to out = cellfun(func, C, 'UniformOutput, false)

out = cellfun(func, varargin{:}, 'UniformOutput', false);

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function display(sp)
% sympoly/display: Display a sympoly object, calls disp

name = inputname(1);
if ~isempty(name)
  display([name,' ='])
else
  display('ans =')
end

% is it a scalar or an array?
s=size(sp);
if any(s>1)
  % an array or vector
  disp(['Sympoly array has size = [',num2str(s),']'])
  disp(' ')

  sj = [0,ones(1,length(s)-1)];
  sj(1) = 0;
  for j=1:numel(sp)
    sj(1) = sj(1)+1;
    while any(sj>s)
      k = find(sj>s,1,'first');
      sj(k) = 1;
      sj(k+1) = sj(k+1)+1;
    end

    disp(['Sympoly array element [',num2str(sj),']'])
    disp(sp(j))
  end

else
  % A scalar

  disp(sp)

end


text:
close all; clear all; clc;

rng default;
rho = 1880/6;
%rho = 320/6;
% Ambient space dimension
M = 9;
% Number of subspaces
K = 5;
% common dimension for each subspace
D = 6;
% dimensions of each subspace
Ds = D * ones(1, K);
bases = spx.data.synthetic.subspaces.random_subspaces(M, K, Ds);
% Number of points on each subspace
Sk = rho * D;
cluster_sizes = Sk * ones(1, K);
% total number of points
S = sum(cluster_sizes);
fprintf('Points per subspace: %d, Total points: %d\n', Sk, S);
% Let us generate uniformly distributed points in each subspace
points_result = spx.data.synthetic.subspaces.uniform_points_on_subspaces(bases, cluster_sizes);
X = points_result.X;
tic;
solver = spx.cluster.ssc.SSC_MC_OMP(X, D, K);
clustering_result = solver.solve();
elapsed_time = toc;
cluster_labels = clustering_result.labels;
true_labels = spx.cluster.labels_from_cluster_sizes(cluster_sizes);
% Time to compare the clustering
comparsion_result = spx.cluster.clustering_error(cluster_labels, true_labels, K);
clustering_error_perc = comparsion_result.error_perc;
clustering_acc_perc = 100 - comparsion_result.error_perc;
% Compute the statistics related to subspace preservation
spr_stats = spx.cluster.subspace.subspace_preservation_stats(clustering_result.Z, cluster_sizes);
spr_error = spr_stats.spr_error;
spr_flag = spr_stats.spr_flag;
spr_perc = spr_stats.spr_perc;
fprintf('\nPoint density: %0.2f: , clustering error: %0.2f %% \n, clustering accuracy: %0.2f %%, mean spr error: %0.2f \npreserving : %0.2f %%, connectivity: %0.2f, \n elapsed time: %0.2f sec', rho, clustering_error_perc, clustering_acc_perc, spr_stats.spr_error, spr_stats.spr_perc, clustering_result.connectivity, elapsed_time);
fprintf('\n\n');


text:
function c = gain_cond_check(model_param, B, psi_1, evmax, c1, c2, kx, kv, kr, kw)

m = model_param.mass;
J = model_param.I;

alpha = sqrt(psi_1*(2-psi_1));
W1 = [ c1*kx/m -c1*kv/(2*m)*(1+alpha);-c1*kv/(2*m)*(1+alpha) kv*(1-alpha)-c1 ];
W12 = [kx*evmax+c1/m*B 0;B 0];
W2 = [c2*kr/max(eig(J)) -c2*kw/(2*min(eig(J)));-c2*kw/(2*min(eig(J))) kw-c2];

c = min(eig(W2))-4*norm(W12)^2/min(eig(W1));

author: yorgoon

repo: minimum-snap-geometric-control


text:
 function printm(varargin)
%function printm(varargin)
%|
%| like printf except that it puts the mfilename in front of it
%| so that you know where the message originated.

[caller, line] = caller_name;
if ~isempty(line) && line ~= 0
	caller = [caller sprintf(' %d', line)];
	caller = strrep(caller, 'LiveEditorEvaluationHelper', 'LiveE*');
end
if length(varargin)
	tmp = varargin{1};
	if streq(tmp, '\n', 2) % trick for newline to precede caller name
		disp(' ') % blank line
		varargin{1} = tmp(3:end); % strip '\n' at start
	end
	disp([caller ': ' sprintf(varargin{:})])
else
	disp([caller ': '])
end
if ir_is_octave
	fflush(stdout);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function varargout = wave_clus_aux2(varargin)
% WAVE_CLUS_AUX1 M-file for wave_clus_aux.fig
%      WAVE_CLUS_AUX, by itself, creates a new WAVE_CLUS_AUX or raises the existing
%      singleton*.
%
%      H = WAVE_CLUS_AUX returns the handle to a new WAVE_CLUS_AUX or the handle to
%      the existing singleton*.
%
%      WAVE_CLUS_AUX('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in WAVE_CLUS_AUX.M with the given input arguments.
%
%      WAVE_CLUS_AUX('Property','Value',...) creates a new WAVE_CLUS_AUX or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before wave_clus_aux_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to wave_clus_aux_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help wave_clus_aux

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @wave_clus_aux2_OpeningFcn, ...
                   'gui_OutputFcn',  @wave_clus_aux2_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before wave_clus_aux is made visible.
function wave_clus_aux2_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to wave_clus_aux (see VARARGIN)

% Choose default command line output for wave_clus_aux
handles.output = hObject;

set(handles.isi14_accept_button,'value',1);
set(handles.isi15_accept_button,'value',1);
set(handles.isi16_accept_button,'value',1);
set(handles.isi17_accept_button,'value',1);
set(handles.isi18_accept_button,'value',1);
set(handles.fix14_button,'value',0);
set(handles.fix15_button,'value',0);
set(handles.fix16_button,'value',0);
set(handles.fix17_button,'value',0);
set(handles.fix18_button,'value',0);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes wave_clus_aux wait for user response (see UIRESUME)
% uiwait(handles.wave_clus_aux);

% --- Outputs from this function are returned to the command line.
function varargout = wave_clus_aux2_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

h_figs=get(0,'children');
h_fig = findobj(h_figs,'tag','wave_clus_figure');
USER_DATA = get(h_fig,'UserData');
par = USER_DATA{1};

for i = 14:18
	si = num2str(i);
	set(eval(['handles.isi' si '_accept_button']),'value',1);
	set(eval(['handles.isi' si '_reject_button']),'value',0);
	set(eval(['handles.fix' si '_button']),'value',0);

	eval(['set(handles.isi' si '_nbins,''string'',par.nbins' si ');']);
	eval(['set(handles.isi' si '_bin_step,''string'',par.bin_step' si ');']);
end

plot_spikes_aux(handles,USER_DATA,2)

author: csn-le


text:
function [ tDigitImg, vDigitLabel ] = ParseMnistData( imgFile, labelFile, trimImg, scaleImg )
% ----------------------------------------------------------------------------------------------- %
%[tDigitImg, vDigitLabel] = ParseMnistData( imgFile, labelFile )
% Parses MNIST data set into a tensor of images and a vector of labels.
% Input:
%   - imgFile           -   Full Path to Images File.
%                           Full path to images file in MNIST format.
%                           Structure: Vector.
%                           Type: 'Char'.
%                           Range: NA.
%   - labelFile         -   Full Path to Labels File.
%                           Full path to labels file in MNIST format.
%                           Structure: Vector.
%                           Type: 'Char'.
%                           Range: NA.
%   - trimImg           -   Trim Image.
%                           If set to 1 will trim images into 20x20
%                           (Center).
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {0, 1}.
%   - scaleImg          -   Scale Image.
%                           If set to 1 will trim images into the [0, 1]
%                           range.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {0, 1}.
% Output:
%   - tDigitImg         -   Digit Images Tensor.
%                           A tensor where each 3rd dimension slice is a
%                           digit image.
%                           Structure: Tensor (numRowx x numCols x numImages).
%                           Type: 'Single' / 'Double'.
%                           Range: [0, 255].
%   - vDigitLabel       -   Image Labels Vector.
%                           A vector where its 'ii' element is the label of
%                           the 'ii' image in the images tensor.
%                           Structure: Vector (numImages X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: [0, 9].
% References
%   1.  A
% Remarks:
%   1.  It works on 2 files: Images, Labels.
% Known Issues:
%   1.  C
% TODO:
%   1.  D
% Release Notes:
%   -   1.0.000     30/09/2021  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

arguments
    imgFile char {mustBeFile, mustBeVector}
    labelFile char {mustBeFile, mustBeVector}
    trimImg (1, 1) {mustBeNumeric, mustBeReal, mustBeInteger, mustBeMember(trimImg, [0, 1])} = 0
    scaleImg (1, 1) {mustBeNumeric, mustBeReal, mustBeInteger, mustBeMember(scaleImg, [0, 1])} = 1
end

IMG_FILE_HEADER     = 2051;
LABEL_FILE_HEADER   = 2049;

% Parse the digit images
fileID = fopen(imgFile, 'r', 'b');
fileHeader = fread(fileID, 1, 'int32');
if(fileHeader ~= IMG_FILE_HEADER)
    error('Invalid digits images file header');
end

numImages = fread(fileID, 1, 'int32');

numRows = fread(fileID, 1, 'int32');
numCols = fread(fileID, 1, 'int32');

tDigitImg = zeros(numRows, numCols, numImages);

for kk = 1:numImages
    for ii = 1:numRows
        tDigitImg(ii, :, kk) = fread(fileID, numCols, 'uint8');
    end
end

% Close file, Finished parsing images
fclose(fileID);

% Read digit labels
fileID = fopen(labelFile, 'r', 'b');
fileHeader = fread(fileID, 1, 'int32');
if(fileHeader ~= LABEL_FILE_HEADER)
    error('Invalid label file header');
end
numLabels = fread(fileID, 1, 'int32');

if(numLabels ~= numImages)
    error('The number of images doesn''t match the number of lables')
end

vDigitLabel = fread(fileID, numLabels, 'uint8');
% Close file, Finished parsing labels
fclose(fileID);

if(trimImg)
    % Returns the center 20x20 image
    tDigitImg = tDigitImg(5:24, 5:24, :);
end

if(scaleImg)
    tDigitImg = tDigitImg ./ 255;
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [ x, seed ] = r8_normal_ab ( a, b, seed )

%*****************************************************************************80
%
%% R8_NORMAL_AB returns a scaled pseudonormal R8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the mean of the normal PDF.
%
%    Input, real B, the standard deviation of the normal PDF.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the standard normal PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ x, seed ] = r8_normal_01 ( seed );
  x = a + b * x;

  return
end


text:
function [cm] = km2cm(km)
% Convert length from kilometers to centimeters.
% Chad A. Greene 2012
cm = km*100000;


text:
function [vector] = createUnitVector(vector)
% Normalize a [1xn] or [nx1] vector
%Description: This function will simplify the process of creating a unit
%vector in one direction by dividing a vector by its length.
%Use: [vectorOut] = fcn_createUnitVector(vectorIn[1xn])
%  or [vectorOut] = fcn_createUnitVector(vectorIn[mx1])

%Author: Jim West
%Date: 7/12/2010
%Updated 8/31/2010
%Made a heavy revision to this script based on feedback from John D'Errico
[r c] = size(vector);

%Convert the vector to columns
if (r > c && c == 1);
    vector = vector';
elseif (r < c && r==1);
    %do Nothing
else
    error('Input vector must be [mx1] or [1xn]');
end

%Process the unit transform
vector = vector./norm(vector);

%Check if the result is illogical, or has an error
if not(sqrt(sum(vector.^2)) < 1.00001) && not(sqrt(sum(vector.^2)) > 0.99999)
    error(['There has been a calculation error, because the square root of the sum of the squares does not equal 1']);
end

%Check if the result is NaN
if sum(isnan(vector)) >= 1
    error(['The vector is ill-conditioned or dividing by zero with a result of ', num2str(vector)]);
end

%If originally in rows, convert back.
if r > c; vector = vector'; end


text:
%  Copyright (c) 2014, Karen Simonyan
%  All rights reserved.
%  This code is made available under the terms of the BSD license (see COPYING file).

function [res, extra] = eval_best(config, scores, gt)

    % finds an optimal threshold - the threshold which maximises the accuracy

    % threshold scores and get the sign

    res = -Inf;
    extra.bestThresh = [];

    % thresh loop
    for i=1:numel(scores)

        curThresh = scores(i);
        class = 2 * (scores >= curThresh) - 1;

        % class-n accuracy
        acc = mean(class == gt);

        if acc > res

            res = acc;
            extra.bestThresh = curThresh;
        end
    end

    res = res * 100;

end

author: AlfredXiangWu


text:
function [sx, dsdx, dsdp] = VBA_sparsifyPrior (x, logExponent, smoothness)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [sx, dsdx, dsdp] = VBA_sparsifyPrior (x, varargin)
% parameter transformation that emulates Laplace priors (L1-norm)
%
% IN:
%   - x: input value to be transformed
%   - logExponent: optional sparsity parameter. If P = 0, then the mapping
%     is regular and no sparsity is emulated. If P=log(2), then the mapping
%     emulates L1-norm like priors.
%   - smoothness: optional smoothness parameter of the sign approximation
%     (default = 1)
%
% In the general case, you won't have to change any of the optional
% parameters. Those options allows to explore the properties of the
% transformation and should not be changed.
%
% OUT:
%   - sx: transformed value
%   - dsdx: gradient of sx wrt x
%   - dsdP: gradient of sx wrt logExponent parameter
%
% /////////////////////////////////////////////////////////////////////////

% check parameters
% =========================================================================
if nargin < 2
    logExponent = log(2);
end

if nargin < 3
    smoothness = []; % use VBA_sign default
end

% shortcuts
% =========================================================================
[signX, d_signX] = VBA_sign(x, smoothness);
[absX, d_absX] = VBA_abs (x);
exponent = exp(logExponent);

% compute transformation
% =========================================================================

sx = signX .* (absX .^ exponent);

% derivatives
% =========================================================================

% wrt to x
% -------------------------------------------------------------------------
if nargout < 2
    return;
end

dsdx = d_signX .* (absX .^ exponent) ...
     + exponent .* absX .^ (exponent - 1) .* signX .* d_absX;

% wrt to logExponent
% -------------------------------------------------------------------------
if nargout < 3
    return;
end

dsdp = (absX .^ exponent) .* signX .* log(absX) .* exponent;

author: MBB-team

repo: VBA-toolbox


text:
function [f] = ma_maxTotalPropMassObjFunc(stateLog, eventID, celBodyData, maData)
%ma_maxTotalPropMass Summary of this function goes here
%   Detailed explanation goes here

    if(ischar(eventID) && strcmpi(eventID,'final'))
        eventNum = max(stateLog(:,13));
    else
        [~, eventNum] = getEventByID(eventID, maData.script);
    end

    eventLog = stateLog(stateLog(:,13)==eventNum,:);
    finalEntry = eventLog(end,:);
    f = -sum(finalEntry(10:12));
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function transpose(f) %#ok<*INUSD>
%TRANSPOSE   CHEBTECH objects are not transposable, so this method will throw an
%   error.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

error('CHEBFUN:CHEBTECH:transpose:notPossible', ...
    'CHEBTECH objects are not transposable.')

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function plot_tf(T,M)

% plot_tf - plot a tensorial field.
%
%   plot_tf(T,M);
%
%   Display min/max eigenvalue and eigensystem.
%   Should be used for debug purpose, not fine-tuned display.
%
%   See also: plot_tensor_field.
%
%   Copyright (c) 2004 Gabriel Peyre

if nargin<2
    M = [];
end

[e1,e2,l1,l2] = perform_tensor_decomp(T);

clf;

subplot(2,2,1);
imagesc(l1);
title('Min eigenvalue.');
axis square;
axis off;

subplot(2,2,2);
imagesc(l2);
title('Max eigenvalue.');
axis square;
axis off;

subplot(2,2,3);
plot_vf(e1, M, 0);
title('Max eigenvector.');

subplot(2,2,4);
plot_vf(e2, M, 0);
title('Max eigenvector.');

author: gpeyre

repo: matlab-toolboxes


text:
function [uu,Up,Uc] = depixellise(pix,cal)

% DEPIXELLISE  Pixellic to metric conversion
%   DEPIXELLISE(PIX,CAL) returns the metric coordinates of an image pixel
%   PIX given camera intrinsic parameters CAL = [u0 v0 au av]'
%
%   [u,Up,Uc] = DEPIXELLISE(...) gives the jacobians wrt PIX and CAL.
%
%   See also INVPINHOLE.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

u0 = cal(1);
v0 = cal(2);
au = cal(3);
av = cal(4);

u = pix(1,:);
v = pix(2,:);

uu = [(u-u0)/au
    (v-v0)/av];

if nargout > 1
    Up = [...
        [ 1/au,    0]
        [    0, 1/av]];
    Uc = [...
        [        -1/au,            0, (-u+u0)/au^2,            0]
        [            0,        -1/av,            0, (-v+v0)/av^2]];
end

return

%% jacobians test

syms u0 v0 au av u v real
pix = [u;v];
cal = [u0;v0;au;av];

[u,Up,Uc] = depixellise(pix,cal);

simplify(Up - jacobian(u,pix))
simplify(Uc - jacobian(u,cal))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
% RDIVIDE_LINTEST   Check that correct linearity information is returned for
% various operations with ./

function pass = test_lintest_rdivide()
%% Initialisation
x = chebfun(@(x) x);
u = adchebfun(x) + 2;
v = adchebfun(x) + 3;

%% One variable involved
pass = [];
w = u/2;
pass(length(pass) + 1) = w.linearity == 1;

w = u./(x+2);
pass(length(pass) + 1) = w.linearity == 1;

w = sin(u)./(x+2);
pass(length(pass) + 1) = w.linearity == 0;

w = 2./u;
pass(length(pass) + 1)  = w.linearity == 0;

w = x./u;
pass(length(pass) + 1)  = w.linearity == 0;

w = u./(u+2);
pass(length(pass) + 1)  = w.linearity == 0;

%% Two variables, but still only one in compuations
u = seed(u, 1, 2);
v = seed(v, 2, 2);
w = u/2;
pass(length(pass) + 1)  = all( w.linearity == [1 1]);

w = u./(x+2);
pass(length(pass) + 1)  = all( w.linearity == [1 1]);

w = sin(u)/2;
pass(length(pass) + 1)  = all( w.linearity == [0 1]);

w = 2./u;
pass(length(pass) + 1)  = all( w.linearity == [0 1]);

w = x./u;
pass(length(pass) + 1)  = all( w.linearity == [0 1]);

w = u./(u+2);
pass(length(pass) + 1)  = all( w.linearity == [0 1]);

%% Combination of variables
w = u./v;
pass(length(pass) + 1)  = all( w.linearity == [0 0]);

w = sin(u)./v;
pass(length(pass) + 1)  = all( w.linearity == [0 0]);

w = u./(u.*v);
pass(length(pass) + 1)  = all( w.linearity == [0 0]);

end


text:
% This file is part of the following project:
% Oliver Parson, Siddhartha Ghosh, Mark Weal, Alex Rogers.
% Non-intrusive Load Monitoring using Prior Models of General Appliance Types.
% In: 26th AAAI Conference on Artificial Intelligence. Toronto, Canada. 2012.
% Code available for download: https://sites.google.com/site/oliparson/phd-work/research-files/aaai-2012-code.zip?attredirects=0
% Copyright: Oliver Parson et al., University of Southhampton, 2012.

% Modified by Romano Cicchetti, ETH Zurich, in the context of the NILM-Eval project

function [ bnet ] = make_hmm( init, emit_mean, emit_cov, trans )

%generic hmm
intra = [0 1; 0 0];
inter = [1 0; 0 0];
hstates = size(init,2); % discrete states
ostates = 1; % continuous
node_sizes = [hstates ostates];
observed_nodes = [2]; % per slice. all others are assumed hidden
discrete_nodes = [1]; % per slice. all others are assumed continuous
eclass1 = [1 2];
eclass2 = [3 2];
eclass = [eclass1 eclass2];
bnet = mk_dbn(intra, inter, node_sizes, ...
'discrete', discrete_nodes, 'observed', observed_nodes, ...
'eclass1', eclass1, 'eclass2', eclass2);

%model params
bnet.CPD{1} = tabular_CPD(bnet, 1, 'CPT', init);
bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', emit_mean, 'cov', emit_cov);
bnet.CPD{3} = tabular_CPD(bnet, 3, 'CPT', trans);

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function brand()
    % brand draws x random samples of size N from the current dataset and computes the b-value
    % sw, last modifies 9/2001
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    ar2 = [];
    arm2 = [];
    br2 = [];
    brm2 = [];

    sdlg.prompt='Minimum number of events per sample ?'; sdlg.value=50;
    sdlg(2).prompt='Step width in events ? '; sdlg(2).value=10;
    sdlg(3).prompt='Maximum number of events per sample?'; sdlg(3).value=200;
    sdlg(4).prompt='Number of samples drawn ?'; sdlg(4).value=100;

    [~,~,n1,ns,n2,nr]=smart_inputdlg('Random b-value calculation', sdlg);

    %n1 = str2double(prmptdlg('Minimum number of events per sample','50'));
    %ns = str2double(prmptdlg('Step width in events','10'));
    %n2 = str2double(prmptdlg('Maximum number of events per sample','200'));
    %nr = str2double(prmptdlg('Numer of samples drawn ','100'));
    tic
    niv = n1:ns:n2;

    globalcatalog=ZG.primeCatalog;

    for ni = n1:ns:n2
        ni
        ar = [];
        arm = [];
        br = [];
        brm = [];
        for i = 1:nr
            l = ceil(rand([ni 1])*globalcatalog.Count);
            %[bv magco stan,  av] =  bvalca3(newa(l,:), McAutoEstimate.manual);
            %br = [br bv];
            %ar = [ar av];
            [bv2 stan av2 ] = calc_bmemag(globalcatalog.Magnitude(l),0.1);
            brm = [brm bv2];
            arm = [arm av2];
        end
        %br2 = [br2 ; br];
        brm2 = [brm2 ; brm];
        %ar2 = [ar2 ; ar];
        arm2 = [arm2 ; arm];
    end

    figure
    pl1 =plot(niv,prctile2(brm2',50),'k')
    set(pl1,'LineWidth',2.0)
    set(gca,'NextPlot','add')
    pl2=plot(niv,prctile2(brm2',95),'r--');
    set(pl2,'LineWidth',1.0,'color',[0.3 0.3 0.3])
    pl3=plot(niv,prctile2(brm2',5),'r-.');
    set(pl3,'LineWidth',1.0,'color',[0.3 0.3 0.3])

    legend([pl1 pl2 pl3],'mean','95%','5%');

    set(gca,'box','on',...
        'SortMethod','childorder','TickDir','out','FontWeight',...
        'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
    xlabel('Number of eqs')
    ylabel('Range of b-value')

    figure

    pl1=plot(niv,prctile2(arm2',50),'k');
    set(pl1,'LineWidth',2.0)
    set(gca,'NextPlot','add')
    pl2=plot(niv,prctile2(arm2',95),'r--');
    set(pl2,'LineWidth',1.0,'color',[0.3 0.3 0.3])
    pl3=plot(niv,prctile2(arm2',5),'r-.');
    set(pl3,'LineWidth',1.0,'color',[0.3 0.3 0.3])
    legend([pl1 pl2 pl3],'mean','95%','5%');

    set(gca,'box','on',...
        'SortMethod','childorder','TickDir','out','FontWeight',...
        'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
    xlabel('Number of eqs')
    ylabel('Range of a-value')
    grid

    toc
    %
    return

    % experimental code ...

    A = [];
    for i = 1:1:99
        i
        A = [A ; niv' prctile2(brm2',i)' niv'*0+i];
    end
    % l = A(:,3)>50; A(l,3) = 100 - A(l,3);
    [ X, Y ] = meshgrid(n1:ns:n2,0.5:0.01:1.5);

    Z = griddata(A(:,1),A(:,2),A(:,3),X,Y);

    figure
    contourf(X,Y,Z,[1 5 10 50 90 95 99]);

    g = gray(6);
    g = g(11:-1:2,:);
    colormap(g);
end

author: CelsoReyes


text:
classdef denoising_mppca < AbstractModel
% denoising_mppca :  4d image denoising and noise map estimation
%                      by exploiting data redundancy in the PCA domain using universal
%                      properties of the eigenspectrum or random covariance matrices,
%                      i.e. Marchenko Pastur distribution
%
% Assumptions:
%   Noise follows a rician distribution
%   image bounderies are not processed
%
% Inputs:
%   Data4D              4D data (any modality)
%  (Mask)               Binary mask with region-of-interest. (OPTIONAL)
%
% Outputs:
%   Data4D_denoised     denoised 4D data
%   sigma_g               standard deviation of the rician noise
%
% Options:
%	sampling
%   	'full'          sliding window
%       'fast'          block processing (warning: undersampled noise map will be returned)
%   kernel              window size, typically in order of [5 x 5 x 5]
%
% Example of command line usage:
%   Model = denoising_mppca;  % Create class from model
%   data.Data4D = load_nii_data('Data4D.nii.gz');  % Load data
%   FitResults = FitData(data,Model,1);  % Fit each voxel within mask
%   FitResultsSave_nii(FitResults,'Data4D.nii.gz');  % Save in local folder: FitResults/
%
% Author: Tanguy Duval, 2016
%
% References:
%   Please cite the following if you use this module:
%     Veraart, J.; Fieremans, E. & Novikov, D.S. Diffusion MRI noise mapping using random matrix theory Magn. Res. Med., 2016, early view, doi:10.1002/mrm.26059
%   In addition to citing the package:
%     Karakuzu A., Boudreau M., Duval T.,Boshkovski T., Leppert I.R., Cabana J.F.,
%     Gagnon I., Beliveau P., Pike G.B., Cohen-Adad J., Stikov N. (2020), qMRLab:
%     Quantitative MRI analysis, under one umbrella doi: 10.21105/joss.02343

    properties (Hidden=true)
        onlineData_url = 'https://osf.io/j5239/download?version=2';
    end

    properties
        MRIinputs = {'Data4D','Mask'};
        xnames = {};
        voxelwise = 0;

        % Protocol
        Prot  = struct(); % You can define a default protocol here.

        % Model options
        buttons = {'sampling',{'fast','full'},'kernel',[5 5 5]};
        options= struct(); % structure filled by the buttons. Leave empty in the code

    end

    methods

        function obj = denoising_mppca
            obj.options = button2opts(obj.buttons);
            obj = UpdateFields(obj);
        end

        function obj = UpdateFields(obj)
        end

        function FitResults = fit(obj,data)
            dims=size(data.Data4D);
            kernel = min(dims(1:3),obj.options.kernel);
            [V,Ind] = min(dims(1:3));
            if ~isfield(data,'Mask'), data.Mask = []; end
            if V<7 && kernel(Ind)>1, helpdlg(['your dataset has very few slices. To avoid loosing too many slices, kernel is set to 1 in the dimension #' num2str(Ind) '.']); kernel(Ind)=1; end
            [FitResults.Data4D_denoised, FitResults.sigma_g] = MPdenoising(data.Data4D,data.Mask,kernel, obj.options.sampling);
        end

    end

    methods(Access = protected)
        function obj = qMRpatch(obj,loadedStruct, version)
            obj = qMRpatch@AbstractModel(obj,loadedStruct, version);
            if checkanteriorver(version,[2 0 11])
                obj.buttons{2} = obj.buttons{2}([2 1]); % old: '|G| (T/m)', new Gnorm (T/m)
            end
        end
    end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function toms443_test01 ( )

%*****************************************************************************80
%
%% TOMS443_TEST01 tests WEW_A
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 June 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS443_TEST01\n' );
  fprintf ( 1, '  Test WEW_A to evaluate\n' );
  fprintf ( 1, '  Lambert''s W function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X      Exact    Computed      Error\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, w1 ] = lambert_w_values ( n_data );

    if ( n_data <= 0 )
      break
    end

    if ( x == 0.0 )
      w2 = 0.0;
    else
      [ w2, en ] = wew_a ( x );
    end

    fprintf ( 1, '  %12.4f  %16.8g  %16.8g  %10.2e\n', ...
      x, w1, w2, abs ( w1 - w2 ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m