text:
function xEst=DopplerOnlyInit6D(rDot,lRx,algorithm,opts,AbsTol,scratchFolderPath,execPath)
%%DOPPLERONLYINIT6D Determine the 3D position and velocity of a target
%          using six simultaneous Doppler measurements. The measurements
%          are assumed to be one-way (from a moving emitter). The receivers
%          are stationary. This function makes use of the function
%          solvePolySysWithExtProg to call an external program to solve the
%          simultaneous multivariate polynomials that define the solution.
%
%INPUTS: rDot A 6X1 or 1X6 vector of Doppler measurements of the emitter.
%     lRx A 3X6 set of the six Cartesian locations of the receivers
%         taking the Doppler measurements in rDot.
% algorithm An optional parameter specifying which algorithm to use. The
%         solvers for simultaneous multivariate polynomials are external
%         programs called via solvePolySysWithExtProg. Possible values
%         are:
%         0 (The default if omitted or an empty matrix is passed) Use
%           Bertini.
%         1 Use PHCpack.
%         2 Use the certified homotopy algorithm that is built into
%           Macaulay2 (in NAG4M2). This only uses the normalized solver
%           with the default options.
%    opts An optional input specifying options for the solver in the
%         function solvePolySysWithExtProg. This is described in more
%         detail in solvePolySysWithExtProg. Omitting this parameter or
%         passing am empty matrices uses the default values, except if
%         Bertini is used, then SecurityMaxNorm, EndpointFiniteThreshold
%         and PathTruncationThreshold are set to 1e9 by default.
%  AbsTol An absolute tolerance on the imaginary part of a solution to the
%         multivariate polynomials used here real. The default if this
%         parameter is omitted or an empty matrix is passed is 1e-8.
% scratchFolderPath An optional parameter specifying the folder in which
%         temporary files will be written. This is needed, because all
%         solvers only works through files. If this parameter is omitted or
%         an empty matrix is passed, then a folder named temp in the folder
%         enclosing this function is used. Note that if an error occurs
%         while executing this function, then temporary files might be left
%         in the temp folder after this function ends.
%  execPath The command line command to use to execute the solver. The
%         default if this parameter is omitted or an empty matrix is
%         passed is just the standard name on the command line: bertini,
%         phc and M2 for each of the algorithms. The default assumes that
%         the program is already in the default search path. For example,
%         on *NIX systems, one can usually put the executable in
%         /usr/local/bin for it to be in the default search path.
%         However, for some reason, that directory is often not
%         included in Matlab's search path. Matlab 2016a's help under
%         "Run External Commands, Scripts, and Programs" specifically
%         says how to add that folder to the default search path.
%
%OUTPUTS: xEst A 6XnumEst set of solutions.
%
%The algorithm solves a system of simultaneous multivariate polynomials as
%formulated in [1] (Equation numbers in the code refer to [1]). However,
%there is a typo in [1], which is corrected in the same derivation in [2].
%
%EXAMPLE:
% lRx=[1000,  500, 1100, 2500,     0, 1000;
%      3000, 2500, 2500,    0,     0, 1000;
%         0,    0,    0,  400,  8000, 8000];%(Stationary) sensor locations
% xTrue=[1e3;5e3;4e3;0;300;0];
% rDot=zeros(6,1);
% for curMeas=1:6
%     diff=xTrue(1:3)-lRx(:,curMeas);
%     rDot(curMeas)=(diff'/norm(diff))*xTrue(4:6);
% end
% DopplerOnlyInit6D(rDot,lRx,0)
%
%REFERENCES:
%[1] T.-L. Lee, S.-S. Lin, W.-W. Lin, S.-T. Yau, and J. Zhu, "Polynomial
%    calculations in Doppler tracking," Communications in Information and
%    Systems, vol. 12, no. 2, pp. 157-184, 2012.
%[2] Y.-C. Kuo, W.-W. Lin, and S.-T. Yau, "A novel efficient homotopy
%    continuation method in tracking," Communications in Information and
%    Systems, vol. 14, no. 1, pp. 57-78, 2014.
%
%March 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%We subtract out the locations of the sensors and then put them back in at
%the end. This should help with numerical sensitivity issues when
%localizing things on the ground but using ECEF coordinates. The offset is
%removed in the end.
centerOfRegion=mean(lRx,2);
lRx=bsxfun(@minus,lRx,centerOfRegion);

%Now, we scale the distances so that the farthest one is 1. This should
%help deal with convergence problems when the scale of the probolem
%changes. The scaling will be reversed after the estimation.
scalFactor=sqrt(max(sum(lRx.*lRx,1)));
lRx=lRx/scalFactor;
rDot=rDot/scalFactor;

if(nargin<3||isempty(algorithm))
    algorithm=0;
end

if(nargin<4)
    opts=[];
end

if(nargin<5||isempty(AbsTol))
    AbsTol=1e-8;
end

if(nargin<5)
    scratchFolderPath=[];
end

if(nargin<6)
    execPath=[];
end

if(algorithm==0&&isempty(opts))%If using Bertini, set default tolerances
    %Default options when using Bertini
    opts=struct('SecurityMaxNorm',1e9,'EndpointFiniteThreshold',1e9,'PathTruncationThreshold',1e9);
elseif(algorithm==0)
    if(~isfield(opts,'SecurityMaxNorm'))
       opts.SecurityMaxNorm=1e9;
    end

    if(~isfield(opts,'EndpointFiniteThreshold'))
        opts.EndpointFiniteThreshold=1e9;
    end

    if(~isfield(opts,'PathTruncationThreshold'))
        opts.PathTruncationThreshold=1e9;
    end
end

%The following are from Equation 4.
V0=lRx';
n=sum(V0.*V0,2);

u1=lRx(:,1);
u2=lRx(:,2);
u3=lRx(:,3);
u4=lRx(:,4);
u5=lRx(:,5);
u6=lRx(:,6);

%From Equation 5
RDot=diag(rDot);

%From Equation 6
C=[-1, 1, 0, 0, 0, 0;
   -1, 0, 1, 0, 0, 0;
   -1, 0, 0, 1, 0, 0;
   -1, 0, 0, 0, 1, 0;
   -1, 0, 0, 0, 0, 1];

%From Equation 10
A=(V0'*(C'*C)*V0)\(V0'*(C'*C));

%From Equation 15
CHat=[0, -1, 1, 0, 0, 0;
      0, -1, 0, 1, 0, 0;
      0, -1, 0, 0, 1, 0;
      0, -1, 0, 0, 0, 1];

%From Equation 16 --Note that the RDot matrix has been removed. This is
%because including RDot makes the solution inconsistent with Equation 22.
%However, the RDot term needs to be included in the QR decomposition.
AHat=([(u2-u3)';
      (u2-u4)';
      (u2-u5)';
      (u2-u6)']*A+CHat);

%Equation 18
u1Hat=(1/2)*A*n-u1;

%Get the T values used in Equation 17 (needs the RDot term with AHat.
[~,R]=qr(AHat*RDot);
T1=R(1:2,1:2);
T2=R(1:2,3:end);

%This is from Equation 17
r12Transform=-T1\T2;
%Extract the elements of the transformation
at11=r12Transform(1,1);
at12=r12Transform(1,2);
at13=r12Transform(1,3);
at14=r12Transform(1,4);
at21=r12Transform(2,1);
at22=r12Transform(2,2);
at23=r12Transform(2,3);
at24=r12Transform(2,4);

%A term in the first formula in Equation 23.
M=(1/2)*(A'*A)*RDot;
%Extract the elements
m1=M(1);
m2=M(2);
m3=M(3);
m4=M(4);
m5=M(5);
m6=M(6);
m7=M(7);
m8=M(8);
m9=M(9);
m10=M(10);
m11=M(11);
m12=M(12);
m13=M(13);
m14=M(14);
m15=M(15);
m16=M(16);
m17=M(17);
m18=M(18);
m19=M(19);
m20=M(20);
m21=M(21);
m22=M(22);
m23=M(23);
m24=M(24);
m25=M(25);
m26=M(26);
m27=M(27);
m28=M(28);
m29=M(29);
m30=M(30);
m31=M(31);
m32=M(32);
m33=M(33);
m34=M(34);
m35=M(35);
m36=M(36);

%A vector term in the first formula in Equation 23.
mVec=u1Hat'*A*RDot;
%Extract the elements
mv1=mVec(1);
mv2=mVec(2);
mv3=mVec(3);
mv4=mVec(4);
mv5=mVec(5);
mv6=mVec(6);

%The first polynomial in Equation 23. The components are ordered
%[r3;r4;r5;r6].
terms1=zeros(5,24);%Allocate space for the terms.

terms1(1,1)=(-at11*(rDot(1)+mv1)-at21*mv2-mv3);
terms1(2:end,1)=[1;0;0;0];%r3

terms1(1,2)=(at11^3*m1+m15+at11*(at21^2*m2+m3)+at11^2*(m13+at21*m7)+at21*(at21*(m14+at21*m8)+m9));
terms1(2:end,2)=[3;0;0;0];%r3^3

terms1(1,3)=(-at12*(rDot(1)+mv1)-at22*mv2-mv4);
terms1(2:end,3)=[0;1;0;0];%r4

terms1(1,4)=(2*at21*at22*m14+m21+at12*m3+at11^2*(3*at12*m1+m19+at22*m7)+2*at11*(at21*at22*m2+at12*(m13+at21*m7))+at21^2*(at12*m2+m20+3*at22*m8)+at22*m9);
terms1(2:end,4)=[2;1;0;0];%r3^2*r4

terms1(1,5)=(at12^2*m13+at22^2*m14+m16+at11*(3*at12^2*m1+at22^2*m2+m4+2*at12*(m19+at22*m7))+at21*(m10+2*at12*at22*m2+2*at22*m20+at12^2*m7+3*at22^2*m8));
terms1(2:end,5)=[1;2;0;0];%r3*r4^2

terms1(1,6)=(at12^3*m1+m22+at12*(at22^2*m2+m4)+at12^2*(m19+at22*m7)+at22*(m10+at22*(m20+at22*m8)));
terms1(2:end,6)=[0;3;0;0];%r4^3

terms1(1,7)=(-at13*(rDot(1)+mv1)-at23*mv2-mv5);
terms1(2:end,7)=[0;0;1;0];%r5

terms1(1,8)=(2*at21*at23*m14+m27+at13*m3+at11^2*(3*at13*m1+m25+at23*m7)+2*at11*(at21*at23*m2+at13*(m13+at21*m7))+at21^2*(at13*m2+m26+3*at23*m8)+at23*m9);
terms1(2:end,8)=[2;0;1;0];%r3^2*r5

terms1(1,9)=2*(at22*at23*m14+at12*(at21*at23*m2+at13*(m13+at21*m7))+at11*(at13*m19+at22*at23*m2+at13*at22*m7+at12*(3*at13*m1+m25+at23*m7))+at21*(at23*m20+at22*(at13*m2+m26+3*at23*m8)));
terms1(2:end,9)=[1;1;1;0];%r3*r4*r5

terms1(1,10)=(at22^2*(at13*m2+m26)+m28+at13*m4+at12^2*(3*at13*m1+m25+at23*m7)+2*at12*(at22*at23*m2+at13*(m19+at22*m7))+at23*(m10+2*at22*m20+3*at22^2*m8));
terms1(2:end,10)=[0;2;1;0];%r4^2*r5

terms1(1,11)=(at13^2*m13+at23^2*m14+m17+at11*(3*at13^2*m1+at23^2*m2+m5+2*at13*(m25+at23*m7))+at21*(m11+2*at13*at23*m2+2*at23*m26+at13^2*m7+3*at23^2*m8));
terms1(2:end,11)=[1;0;2;0];%r3*r5^2

terms1(1,12)=(at22*m11+2*at13*at22*at23*m2+m23+at13^2*(m19+at22*m7)+at12*(3*at13^2*m1+at23^2*m2+m5+2*at13*(m25+at23*m7))+at23*(2*at22*m26+at23*(m20+3*at22*m8)));
terms1(2:end,12)=[0;1;2;0];%r4*r5^2;

terms1(1,13)=(at13^3*m1+m29+at13*(at23^2*m2+m5)+at13^2*(m25+at23*m7)+at23*(m11+at23*(m26+at23*m8)));
terms1(2:end,13)=[0;0;3;0];%r5^3

terms1(1,14)=(-at14*(rDot(1)+mv1)-at24*mv2-mv6);
terms1(2:end,14)=[0;0;0;1];%r6

terms1(1,15)=(2*at21*at24*m14+at14*m3+m33+at11^2*(3*at14*m1+m31+at24*m7)+2*at11*(at21*at24*m2+at14*(m13+at21*m7))+at21^2*(at14*m2+m32+3*at24*m8)+at24*m9);
terms1(2:end,15)=[2;0;0;1];%r3^2*r6

terms1(1,16)=2*(at22*at24*m14+at12*(at21*at24*m2+at14*(m13+at21*m7))+at11*(at14*m19+at22*at24*m2+at14*at22*m7+at12*(3*at14*m1+m31+at24*m7))+at21*(at24*m20+at22*(at14*m2+m32+3*at24*m8)));
terms1(2:end,16)=[1;1;0;1];%r3*r4*r6

terms1(1,17)=(at22^2*(at14*m2+m32)+m34+at14*m4+at12^2*(3*at14*m1+m31+at24*m7)+2*at12*(at22*at24*m2+at14*(m19+at22*m7))+at24*(m10+2*at22*m20+3*at22^2*m8));
terms1(2:end,17)=[0;2;0;1];%r4^2*r6

terms1(1,18)=2*(at23*at24*m14+at13*(at21*at24*m2+at14*(m13+at21*m7))+at11*(at23*at24*m2+at14*m25+at14*at23*m7+at13*(3*at14*m1+m31+at24*m7))+at21*(at24*m26+at23*(at14*m2+m32+3*at24*m8)));
terms1(2:end,18)=[1;0;1;1];%r3*r5*r6

terms1(1,19)=2*(at14*at22*at23*m2+at23*at24*m20+at22*at24*m26+at22*at23*m32+at13*(at14*m19+at22*at24*m2+at14*at22*m7)+at12*(at23*at24*m2+at14*m25+at14*at23*m7+at13*(3*at14*m1+m31+at24*m7))+3*at22*at23*at24*m8);
terms1(2:end,19)=[0;1;1;1];%r4*r5*r6

terms1(1,20)=(at23^2*(at14*m2+m32)+m35+at14*m5+2*at13*(at23*at24*m2+at14*m25+at14*at23*m7)+at13^2*(3*at14*m1+m31+at24*m7)+at24*(m11+2*at23*m26+3*at23^2*m8));
terms1(2:end,20)=[0;0;2;1];%r5^2*r6

terms1(1,21)=(at14^2*m13+at24^2*m14+m18+at11*(3*at14^2*m1+at24^2*m2+m6+2*at14*(m31+at24*m7))+at21*(m12+2*at14*at24*m2+2*at24*m32+at14^2*m7+3*at24^2*m8));
terms1(2:end,21)=[1;0;0;2];%r3*r6^2

terms1(1,22)=(at22*m12+2*at14*at22*at24*m2+m24+at14^2*(m19+at22*m7)+at12*(3*at14^2*m1+at24^2*m2+m6+2*at14*(m31+at24*m7))+at24*(2*at22*m32+at24*(m20+3*at22*m8))) ;
terms1(2:end,22)=[0;1;0;2];%r4*r6^2

terms1(1,23)=(at23*m12+2*at14*at23*at24*m2+m30+at14^2*(m25+at23*m7)+at13*(3*at14^2*m1+at24^2*m2+m6+2*at14*(m31+at24*m7))+at24*(2*at23*m32+at24*(m26+3*at23*m8)));
terms1(2:end,23)=[0;0;1;2];%r5*r6^2

terms1(1,24)=(at14^3*m1+m36+at14*(at24^2*m2+m6)+at14^2*(m31+at24*m7)+at24*(m12+at24*(m32+at24*m8)));
terms1(2:end,24)=[0;0;0;3];%r6^3

%Scale the coefficients
terms1(1,:)=terms1(1,:)/max(abs(terms1(1,:)));

%cHat is defined after Equation 20.
cHat=(1/4)*n'*(A'*A)*n-u1'*A*n+u1'*u1;

%A term in the second formula in Equation 23.
Mt=(1/4)*(A'*A);
%Extract the elements
mt1=Mt(1);
mt2=Mt(2);
mt3=Mt(3);
mt4=Mt(4);
mt5=Mt(5);
mt6=Mt(6);
mt7=Mt(7);
mt8=Mt(8);
mt9=Mt(9);
mt10=Mt(10);
mt11=Mt(11);
mt12=Mt(12);
mt13=Mt(13);
mt14=Mt(14);
mt15=Mt(15);
mt16=Mt(16);
mt17=Mt(17);
mt18=Mt(18);
mt19=Mt(19);
mt20=Mt(20);
mt21=Mt(21);
mt22=Mt(22);
mt23=Mt(23);
mt24=Mt(24);
mt25=Mt(25);
mt26=Mt(26);
mt27=Mt(27);
mt28=Mt(28);
mt29=Mt(29);
mt30=Mt(30);
mt31=Mt(31);
mt32=Mt(32);
mt33=Mt(33);
mt34=Mt(34);
mt35=Mt(35);
mt36=Mt(36);

%A vector term in the second formula in Equation 23.
mVec2=u1Hat'*A;
%Extract the elements
mvt1=mVec2(1);
mvt2=mVec2(2);
mvt3=mVec2(3);
mvt4=mVec2(4);
mvt5=mVec2(5);
mvt6=mVec2(6);

%The second polynomial in Equation 23. The components are ordered
%[r3;r4;r5;r6].
terms2=zeros(5,46);%Allocate space for the terms.

terms2(1,1)=cHat;
terms2(2:end,1)=[0;0;0;0];%Constant term

terms2(1,2)=(at11^4*mt1+mt15+at11^2*(mt13+mt3+at21^2*(mt2+mt7))+at21^4*mt8+at21^2*(mt14+mt9));
terms2(2:end,2)=[4;0;0;0];%r3^4

terms2(1,3)=2*(2*at11^3*at12*mt1+at11^2*at21*at22*(mt2+mt7)+at11*at12*(mt13+mt3+at21^2*(mt2+mt7))+at21*at22*(mt14+2*at21^2*mt8+mt9));
terms2(2:end,3)=[3;1;0;0];%r3^3*r4

terms2(1,4)=(mt16+mt21+at12^2*(mt13+mt3)+4*at11*at12*at21*at22*(mt2+mt7)+at11^2*(6*at12^2*mt1+mt19+mt4+at22^2*(mt2+mt7))+at21^2*(mt10+mt20+at12^2*(mt2+mt7)+6*at22^2*mt8)+at22^2*(mt14+mt9));
terms2(2:end,4)=[2;2;0;0];%r3^2*r4^2

terms2(1,5)=2*(at11*at12*(2*at12^2*mt1+mt19+mt4+at22^2*(mt2+mt7))+at21*at22*(mt10+mt20+at12^2*(mt2+mt7)+2*at22^2*mt8));
terms2(2:end,5)=[1;3;0;0];%r3*r4^3

terms2(1,6)=(at12^4*mt1+at22^2*(mt10+mt20)+mt22+at12^2*(mt19+mt4+at22^2*(mt2+mt7))+at22^4*mt8);
terms2(2:end,6)=[0;4;0;0];%r4^4

terms2(1,7)=2*(2*at11^3*at13*mt1+at11^2*at21*at23*(mt2+mt7)+at11*at13*(mt13+mt3+at21^2*(mt2+mt7))+at21*at23*(mt14+2*at21^2*mt8+mt9));
terms2(2:end,7)=[3;0;1;0];%r3^3*r5

terms2(1,8)=2*(2*at12^3*at13*mt1+at12^2*at22*at23*(mt2+mt7)+at12*at13*(mt19+mt4+at22^2*(mt2+mt7))+at22*at23*(mt10+mt20+2*at22^2*mt8));
terms2(2:end,8)=[0;3;1;0];%r4^3*r5

terms2(1,9)=(mt17+mt27+at13^2*(mt13+mt3)+4*at11*at13*at21*at23*(mt2+mt7)+at11^2*(6*at13^2*mt1+mt25+mt5+at23^2*(mt2+mt7))+at21^2*(mt11+mt26+at13^2*(mt2+mt7)+6*at23^2*mt8)+at23^2*(mt14+mt9));
terms2(2:end,9)=[2;0;2;0];%r3^2*r5^2

terms2(1,10)=(mt23+mt28+at13^2*(mt19+mt4)+4*at12*at13*at22*at23*(mt2+mt7)+at22^2*(mt11+mt26+at13^2*(mt2+mt7))+at12^2*(6*at13^2*mt1+mt25+mt5+at23^2*(mt2+mt7))+at23^2*(mt10+mt20+6*at22^2*mt8));
terms2(2:end,10)=[0;2;2;0];%r4^2*r5^2

terms2(1,11)=2*(at11*at13*(2*at13^2*mt1+mt25+mt5+at23^2*(mt2+mt7))+at21*at23*(mt11+mt26+at13^2*(mt2+mt7)+2*at23^2*mt8));
terms2(2:end,11)=[1;0;3;0];%r3*r5^3

terms2(1,12)=2*(at12*at13*(2*at13^2*mt1+mt25+mt5+at23^2*(mt2+mt7))+at22*at23*(mt11+mt26+at13^2*(mt2+mt7)+2*at23^2*mt8));
terms2(2:end,12)=[0;1;3;0];%r4*r5^3

terms2(1,13)=(at13^4*mt1+at23^2*(mt11+mt26)+mt29+at13^2*(mt25+mt5+at23^2*(mt2+mt7))+at23^4*mt8);
terms2(2:end,13)=[0;0;4;0];%r5^4

terms2(1,14)=-2*(at11*at14*(1+mvt1)+at21*at24*mvt2);
terms2(2:end,14)=[1;0;0;1];%r3*r6

terms2(1,15)=2*(2*at11^3*at14*mt1+at11^2*at21*at24*(mt2+mt7)+at11*at14*(mt13+mt3+at21^2*(mt2+mt7))+at21*at24*(mt14+2*at21^2*mt8+mt9));
terms2(2:end,15)=[3;0;0;1];%r3^3*r6

terms2(1,16)=-2*(at12*at14*(1+mvt1)+at22*at24*mvt2);
terms2(2:end,16)=[0;1;0;1];%r4*r6

terms2(1,17)=2*(2*at12^3*at14*mt1+at12^2*at22*at24*(mt2+mt7)+at12*at14*(mt19+mt4+at22^2*(mt2+mt7))+at22*at24*(mt10+mt20+2*at22^2*mt8));
terms2(2:end,17)=[0;3;0;1];%r4^3*r6

terms2(1,18)=-2*(at13*at14*(1+mvt1)+at23*at24*mvt2);
terms2(2:end,18)=[0;0;1;1];%r5*r6

terms2(1,19)=2*(2*at13^3*at14*mt1+at13^2*at23*at24*(mt2+mt7)+at13*at14*(mt25+mt5+at23^2*(mt2+mt7))+at23*at24*(mt11+mt26+2*at23^2*mt8));
terms2(2:end,19)=[0;0;3;1];%r5^3*r6

terms2(1,20)=(-at14^2*(1+mvt1)-at24^2*mvt2-mvt6);
terms2(2:end,20)=[0;0;0;2];%r6^2

terms2(1,21)=(mt18+at14^2*(mt13+mt3)+mt33+4*at11*at14*at21*at24*(mt2+mt7)+at11^2*(6*at14^2*mt1+mt31+mt6+at24^2*(mt2+mt7))+at21^2*(mt12+mt32+at14^2*(mt2+mt7)+6*at24^2*mt8)+at24^2*(mt14+mt9));
terms2(2:end,21)=[2;0;0;2];%r3^2*r6^2

terms2(1,22)=(mt24+mt34+at14^2*(mt19+mt4)+4*at12*at14*at22*at24*(mt2+mt7)+at22^2*(mt12+mt32+at14^2*(mt2+mt7))+at12^2*(6*at14^2*mt1+mt31+mt6+at24^2*(mt2+mt7))+at24^2*(mt10+mt20+6*at22^2*mt8));
terms2(2:end,22)=[0;2;0;2];%r4^2*r6^2

terms2(1,23)=(mt30+mt35+at14^2*(mt25+mt5)+4*at13*at14*at23*at24*(mt2+mt7)+at23^2*(mt12+mt32+at14^2*(mt2+mt7))+at13^2*(6*at14^2*mt1+mt31+mt6+at24^2*(mt2+mt7))+at24^2*(mt11+mt26+6*at23^2*mt8));
terms2(2:end,23)=[0;0;2;2];%r5^2*r6^2

terms2(1,24)=2*(at11*at14*(2*at14^2*mt1+mt31+mt6+at24^2*(mt2+mt7))+at21*at24*(mt12+mt32+at14^2*(mt2+mt7)+2*at24^2*mt8));
terms2(2:end,24)=[1;0;0;3];%r3*r6^3

terms2(1,25)=2*(at12*at14*(2*at14^2*mt1+mt31+mt6+at24^2*(mt2+mt7))+at22*at24*(mt12+mt32+at14^2*(mt2+mt7)+2*at24^2*mt8));
terms2(2:end,25)=[0;1;0;3];%r4*r6^3

terms2(1,26)=2*(at13*at14*(2*at14^2*mt1+mt31+mt6+at24^2*(mt2+mt7))+at23*at24*(mt12+mt32+at14^2*(mt2+mt7)+2*at24^2*mt8));
terms2(2:end,26)=[0;0;1;3];%r5*r6^3

terms2(1,27)=(at14^4*mt1+at24^2*(mt12+mt32)+mt36+at14^2*(mt31+mt6+at24^2*(mt2+mt7))+at24^4*mt8);
terms2(2:end,27)=[0;0;0;4];%r6^4

terms2(1,28)=(-at13^2*(1+mvt1)-at23^2*mvt2-mvt5);
terms2(2:end,28)=[0;0;2;0];%r5^2

terms2(1,29)=2*(at11*(6*at13^2*at14*mt1+2*at13*at23*at24*(mt2+mt7)+at14*(mt25+mt5+at23^2*(mt2+mt7)))+at21*(2*at13*at14*at23*(mt2+mt7)+at24*(mt11+mt26+at13^2*(mt2+mt7)+6*at23^2*mt8)));
terms2(2:end,29)=[1;0;2;1];%r3*r5^2*r6

terms2(1,30)=2*(at12*(6*at13^2*at14*mt1+2*at13*at23*at24*(mt2+mt7)+at14*(mt25+mt5+at23^2*(mt2+mt7)))+at22*(2*at13*at14*at23*(mt2+mt7)+at24*(mt11+mt26+at13^2*(mt2+mt7)+6*at23^2*mt8)));
terms2(2:end,30)=[0;1;2;1];%r4*r5^2*r6

terms2(1,31)=(-at12^2*(1+mvt1)-at22^2*mvt2-mvt4);
terms2(2:end,31)=[0;2;0;0];%r4^2

terms2(1,32)=2*(at11*(6*at12^2*at13*mt1+2*at12*at22*at23*(mt2+mt7)+at13*(mt19+mt4+at22^2*(mt2+mt7)))+at21*(2*at12*at13*at22*(mt2+mt7)+at23*(mt10+mt20+at12^2*(mt2+mt7)+6*at22^2*mt8)));
terms2(2:end,32)=[1;2;1;0];%r3*r4^2*r5

terms2(1,33)=2*(at11*(6*at12^2*at14*mt1+2*at12*at22*at24*(mt2+mt7)+at14*(mt19+mt4+at22^2*(mt2+mt7)))+at21*(2*at12*at14*at22*(mt2+mt7)+at24*(mt10+mt20+at12^2*(mt2+mt7)+6*at22^2*mt8)));
terms2(2:end,33)=[1;2;0;1];%r3*r4^2*r6

terms2(1,34)=2*(2*at12*at22*(at14*at23+at13*at24)*(mt2+mt7)+at13*at14*(mt19+mt4+at22^2*(mt2+mt7))+at12^2*(6*at13*at14*mt1+at23*at24*(mt2+mt7))+at23*at24*(mt10+mt20+6*at22^2*mt8));
terms2(2:end,34)=[0;2;1;1];%r4^2*r5*r6

terms2(1,35)=(-at11^2*(1+mvt1)-at21^2*mvt2-mvt3);
terms2(2:end,35)=[2;0;0;0];%r3^2

terms2(1,36)=2*(2*at11*at21*(at13*at22+at12*at23)*(mt2+mt7)+at12*at13*(mt13+mt3+at21^2*(mt2+mt7))+at11^2*(6*at12*at13*mt1+at22*at23*(mt2+mt7))+at22*at23*(mt14+6*at21^2*mt8+mt9));
terms2(2:end,36)=[2;1;1;0];%r3^2*r4*r5

terms2(1,37)=2*(2*at11*at21*(at14*at22+at12*at24)*(mt2+mt7)+at12*at14*(mt13+mt3+at21^2*(mt2+mt7))+at11^2*(6*at12*at14*mt1+at22*at24*(mt2+mt7))+at22*at24*(mt14+6*at21^2*mt8+mt9));
terms2(2:end,37)=[2;1;0;1];%r3^2*r4*r6

terms2(1,38)=2*(2*at11*at21*(at14*at23+at13*at24)*(mt2+mt7)+at13*at14*(mt13+mt3+at21^2*(mt2+mt7))+at11^2*(6*at13*at14*mt1+at23*at24*(mt2+mt7))+at23*at24*(mt14+6*at21^2*mt8+mt9));
terms2(2:end,38)=[2;0;1;1];%r3^2*r5*r6

terms2(1,39)=-2*(at11*at12*(1+mvt1)+at21*at22*mvt2);
terms2(2:end,39)=[1;1;0;0];%r3*r4

terms2(1,40)=2*(at11*(2*at13*at22*at23*(mt2+mt7)+at12*(6*at13^2*mt1+mt25+mt5+at23^2*(mt2+mt7)))+at21*(2*at12*at13*at23*(mt2+mt7)+at22*(mt11+mt26+at13^2*(mt2+mt7)+6*at23^2*mt8)));
terms2(2:end,40)=[1;1;2;0];%r3*r4*r5^2

terms2(1,41)=4*(at21*(at13*at14*at22+at12*at14*at23+at12*at13*at24)*(mt2+mt7)+at11*(at22*(at14*at23+at13*at24)*(mt2+mt7)+at12*(6*at13*at14*mt1+at23*at24*(mt2+mt7)))+6*at21*at22*at23*at24*mt8);
terms2(2:end,41)=[1;1;1;1];%r3*r4*r5*r6

terms2(1,42)=2*(at11*(2*at14*at22*at24*(mt2+mt7)+at12*(6*at14^2*mt1+mt31+mt6+at24^2*(mt2+mt7)))+at21*(2*at12*at14*at24*(mt2+mt7)+at22*(mt12+mt32+at14^2*(mt2+mt7)+6*at24^2*mt8)));
terms2(2:end,42)=[1;1;0;2];%r3*r4*r6^2

terms2(1,43)=-2*(at11*at13*(1+mvt1)+at21*at23*mvt2);
terms2(2:end,43)=[1;0;1;0];%r3*r5

terms2(1,44)=2*(at11*(2*at14*at23*at24*(mt2+mt7)+at13*(6*at14^2*mt1+mt31+mt6+at24^2*(mt2+mt7)))+at21*(2*at13*at14*at24*(mt2+mt7)+at23*(mt12+mt32+at14^2*(mt2+mt7)+6*at24^2*mt8)));
terms2(2:end,44)=[1;0;1;2];%r3*r5*r6^2

terms2(1,45)=-2*(at12*at13*(1+mvt1)+at22*at23*mvt2);
terms2(2:end,45)=[0;1;1;0];%r4*r5

terms2(1,46)=2*(at12*(2*at14*at23*at24*(mt2+mt7)+at13*(6*at14^2*mt1+mt31+mt6+at24^2*(mt2+mt7)))+at22*(2*at13*at14*at24*(mt2+mt7)+at23*(mt12+mt32+at14^2*(mt2+mt7)+6*at24^2*mt8)));
terms2(2:end,46)=[0;1;1;2];%r4*r5*r6^2

%Scale the coefficients
terms2(1,:)=terms2(1,:)/max(abs(terms2(1,:)));

%The third polynomial in Equation 23:
%A vector term in the third formula in Equation 23.
mVec3=(u2-u3)'*A;

%Extract the elements
mvtt1=mVec3(1);
mvtt2=mVec3(2);
mvtt3=mVec3(3);
mvtt4=mVec3(4);
mvtt5=mVec3(5);
mvtt6=mVec3(6);

terms3=zeros(5,11);

terms3(1,1)=(-u3'*u3+u2'*u2-(u2-u3)'*A*n);
terms3(2:end,1)=[0;0;0;0];%Constant term

terms3(1,2)=1+at11^2*mvtt1+at21^2*(-1+mvtt2)+mvtt3;
terms3(2:end,2)=[2;0;0;0];%r3^2

terms3(1,3)=2*(at11*at12*mvtt1+at21*at22*(-1+mvtt2));
terms3(2:end,3)=[1;1;0;0];%r3*r4

terms3(1,4)=(at12^2*mvtt1+at22^2*(-1+mvtt2)+mvtt4);
terms3(2:end,4)=[0;2;0;0];%r4^2

terms3(1,5)=2*(at11*at13*mvtt1+at21*at23*(-1+mvtt2));
terms3(2:end,5)=[1;0;1;0];%r3*r5

terms3(1,6)=2*(at12*at13*mvtt1+at22*at23*(-1+mvtt2));
terms3(2:end,6)=[0;1;1;0];%r4*r5

terms3(1,7)=(at13^2*mvtt1+at23^2*(-1+mvtt2)+mvtt5);
terms3(2:end,7)=[0;0;2;0];%r5^2

terms3(1,8)=2*(at11*at14*mvtt1+at21*at24*(-1+mvtt2));
terms3(2:end,8)=[1;0;0;1];%r3*r6

terms3(1,9)=2*(at12*at14*mvtt1+at22*at24*(-1+mvtt2));
terms3(2:end,9)=[0;1;0;1];%r4*r6

terms3(1,10)=2*(at13*at14*mvtt1+at23*at24*(-1+mvtt2));
terms3(2:end,10)=[0;0;1;1];%r5*r6

terms3(1,11)=(at14^2*mvtt1+at24^2*(-1+mvtt2)+mvtt6);
terms3(2:end,11)=[0;0;0;2];%r6^2

%Scale the coefficients
terms3(1,:)=terms3(1,:)/max(abs(terms3(1,:)));

%The fourth polynomial in Equation 23:
mVec4=(u2-u4)'*A;
%Extract the elements
mvtt1=mVec4(1);
mvtt2=mVec4(2);
mvtt3=mVec4(3);
mvtt4=mVec4(4);
mvtt5=mVec4(5);
mvtt6=mVec4(6);

terms4=zeros(5,11);

terms4(1,1)=(-u4'*u4+u2'*u2-(u2-u4)'*A*n);
terms4(2:end,1)=[0;0;0;0];%Constant term

terms4(1,2)=(at11^2*mvtt1 + at21^2*(-1 + mvtt2) + mvtt3) ;
terms4(2:end,2)=[2;0;0;0];%r3^2

terms4(1,3)=2*(at11*at12*mvtt1 + at21*at22*(-1 + mvtt2)) ;
terms4(2:end,3)=[1;1;0;0];%r3*r4

terms4(1,4)=(1 + at12^2*mvtt1 + at22^2*(-1 + mvtt2) + mvtt4);
terms4(2:end,4)=[0;2;0;0];%r4^2

terms4(1,5)=2*(at11*at13*mvtt1 + at21*at23*(-1 + mvtt2));
terms4(2:end,5)=[1;0;1;0];%r3*r5

terms4(1,6)=2*(at12*at13*mvtt1 + at22*at23*(-1 + mvtt2));
terms4(2:end,6)=[0;1;1;0];%r4*r5

terms4(1,7)=(at13^2*mvtt1 + at23^2*(-1 + mvtt2) + mvtt5);
terms4(2:end,7)=[0;0;2;0];%r5^2

terms4(1,8)=2*(at11*at14*mvtt1 + at21*at24*(-1 + mvtt2));
terms4(2:end,8)=[1;0;0;1];%r3*r6

terms4(1,9)=2*(at12*at14*mvtt1 + at22*at24*(-1 + mvtt2));
terms4(2:end,9)=[0;1;0;1];%r4*r6

terms4(1,10)=2*(at13*at14*mvtt1 + at23*at24*(-1 + mvtt2));
terms4(2:end,10)=[0;0;1;1];%r5*r6

terms4(1,11)=(at14^2*mvtt1 + at24^2*(-1 + mvtt2) + mvtt6);
terms4(2:end,11)=[0;0;0;2];%r6^2

%Scale the coefficients
terms4(1,:)=terms4(1,:)/max(abs(terms4(1,:)));

%Now, put the term matrices into formats that can be used by the
%multivariate polynomial solvers.
xPolys=cell(4,1);
xPolys{1}=terms1;
xPolys{2}=terms2;
xPolys{3}=terms3;
xPolys{4}=terms4;

varNames={'rc','rd','re','rf'};

for curPoly=1:4
    termMat=xPolys{curPoly};

    numTerms=size(termMat,2);

    thePoly=[];

    for curTerm=1:numTerms
        curCoeff=num2str(termMat(1,curTerm),16);

        curMonomial=[];

        for curVar=1:4
            if(termMat(curVar+1,curTerm)~=0)
                if(~isempty(curMonomial))
                   curMonomial=[curMonomial,'*'];
                end

                if(termMat(curVar+1,curTerm)==1)
                    curMonomial=[curMonomial,varNames{curVar}];
                else
                    curMonomial=[curMonomial,varNames{curVar},'^',num2str(termMat(curVar+1,curTerm))];
                end
            end
        end

        if(isempty(curMonomial))%If it is a constant term.
            %The constant term (if there is one) should be the first
            %term, so we can just set thePoly to it.
            thePoly=curCoeff;
        elseif(~isempty(thePoly))
            thePoly=[thePoly,'+(',curCoeff,')*',curMonomial];
        else%If there is no constant term.
            thePoly=['(',curCoeff,')*',curMonomial];
        end
    end

    xPolys{curPoly}=thePoly;
end

rEst=solvePolySysWithExtProg(xPolys,varNames,algorithm,opts,scratchFolderPath,execPath);

if(isempty(rEst))
    %The solver failed.
    xEst=[];
    return;
end
%Throw out complex solutions. Solutions are deemed complex if the
%imaginary part exceeds AbsTol in magnitude.
sel=all(abs(imag(rEst))<AbsTol,1);
rEst=real(rEst(:,sel));

%Only positive range values are valid.
sel=all(rEst>0,1);
rEst=rEst(:,sel);

numSol=size(rEst,2);
xEst=zeros(6,numSol);

%Given solutions to the r variables, we need to extract the state.
numAdded=0;
for curSol=1:numSol
    %From Equation 17
    r1r2=r12Transform*rEst(:,curSol);

    %Only positive ranges are valid.
    if(all(r1r2)>0)
        r=[r1r2;rEst(:,curSol)];

        %Equation 8
        uDot=A*(-RDot*r);

        %Equation 9
        u=(1/2)*A*(n-r.*r);

        numAdded=numAdded+1;
        xEst(:,numAdded)=[u;uDot];
    end
end

%Shrink to fit the actual number of solutions added.
xEst=xEst(:,1:numAdded);

%Undo the effects of scaling:
xEst=scalFactor*xEst;

%Add back in the offset that was present due to centering the coordinate
%system around the sensors.
xEst(1:3,:)=bsxfun(@plus,xEst(1:3,:),centerOfRegion);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function test_example_glm_nirs

% MEM 4gb
% WALLTIME 00:10:00

%% Using GLM to analyze NIRS timeseries data
%
% This is an example MATLAB script that demonstrates how to compute a simple GLM on the fingertapping NIRS data that is also used in the tutorial on [preprocessing and averaging of single-channel NIRS data](/tutorial/nirs_singlechannel). The data is available from our [FTP server](ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/tutorial/nirs_singlechannel/).
%
% We start with reading the NIRS data from disk.
%
cfg = [];
cfg.dataset = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/nirs_singlechannel/motor_cortex.oxy3');
data_nirs = ft_preprocessing(cfg);

%% # Construct a number of additional channels with the task/stimulus details
%
% Then we read the events, which are represented as markers or triggers, indicating the samples in the data in which the fingertapping started or stopped. We use these to make some additional continuously represented channels that represent the onset, offset, and the motion. Furthermore, we can add two channels for a constant offset, and for a slope. These can be used to remove the baseline and a constant drift in the signal over time.
%
event = ft_read_event(cfg.dataset);

data_stim = [];
data_stim.time = data_nirs.time;
data_stim.label = {
  'onset'
  'offset'
  'motion'
  'constant'
  'slope'
  };
data_stim.chantype = repmat({'stimulus'}, size(data_stim.label));  % Homer stires the experimental design as stimulus (s)
data_stim.chanunit = repmat({'unknown'}, size(data_stim.label));

move_onset  = [event(strcmp({event.value}, 'A')).sample]; % this indicates the beginning of the movement
move_offset = [event(strcmp({event.value}, 'B')).sample]; % this indicates the end of the movement

nchans   = length(data_stim.label);
nsamples = length(data_stim.time{1}); % it is a continous representation, hence one trial/segment

data_stim.trial{1} = zeros(nchans, nsamples);
data_stim.trial{1}(1,move_onset)  = 1;
data_stim.trial{1}(2,move_offset) = 1;
for i=1:numel(move_onset)
  data_stim.trial{1}(3,move_onset(i):move_offset(i)) = 1;
end
data_stim.trial{1}(4,:) = ones(1,nsamples);
data_stim.trial{1}(5,:) = linspace(0,1,nsamples);

% show the experimental design as a matrix
figure
imagesc(data_stim.trial{1})

%
% You actually have to zoom in a lot to see all details, since there are more samples than horizontal pixels on your screen
%
%% # Combine the NIRS data and the description of the task/stimulus details
%
cfg = [];
data_combined = ft_appenddata(cfg, data_nirs, data_stim);

% also keep the optode structure, we need it further down for fieldtrip2homer and plotting
data_combined.opto = data_nirs.opto;

%% # Perform a GLM analysis
%
% This is explained on http://mri-q.com/general-linear-model.html with an excelent introduction, and
% on https://www.brainvoyager.com/bv/doc/UsersGuide/StatisticalAnalysis/TheGeneralLinearModel.html.
%
close all

% take the data from the FieldTrip data structures
y = cat(2, data_nirs.trial{:})';
x = cat(2, data_stim.trial{:})';

% After further processing and artifact removal you could also take them from the `data_combined` structure, which would ensure that they keep nicely aligned - also when you discard artifact segments. But here we simply include all samples.
%
% use the default canonical haemodynamic response function (HRF) from SPM
ft_hastoolbox('spm12', 1);
h = spm_hrf(1/data_nirs.fsample);

% plot the original regressor
figure
plot(x(:,3), 'b')
hold on
title('model')

% convolve the design with the HRF
xc = convn(x', h')';
xc = xc(1:nsamples,:); % remove the trailing part
% replace the first three columns of the design with the HRF convolved version
x(:,1:3) = xc(:,1:3);

% plot the HRF convolved regressor
plot(x(:,3), 'r-')

%
% Remove the baseline, see also **[ft_preproc_baselinecorrect](https://github.com/fieldtrip/fieldtrip/blob/release/preproc/ft_preproc_baselinecorrect.m)** and **[ft_preproc_detrend](https://github.com/fieldtrip/fieldtrip/blob/release/preproc/ft_preproc_detrend.m)**. This is beneficial here, since the experimental regressors are not orthogonal to the confound regressors.
%
y = ft_preproc_polyremoval(y', 2)';

% fit the general linear model (GLM)
% y = x * b + e
b     = x \ y;
model = x * b;
e     = y - model;

chan = find(strcmp(data_nirs.label, 'Rx4b-Tx5 [860nm]'));

figure
hold on
plot(y(:,chan), 'b')
plot(e(:,chan), 'g')
plot(model(:,chan), 'r', 'linewidth', 2)
%set(h, 'LineWidth', 1)
legend({'data', 'noise', 'model'});
title('fit to data')

%
%% # Convert the fitted model into statistical parameters
%
[n,p] = size(x);
r2 = var(model)./var(y);
f = r2*(n-p) ./ ((1-r2)*(p-1));
c = [0 0 1 0 0]';
t = (c' * b) ./ sqrt(var(e) * (c' * (x'*x) * c));

% It would be possible to make topographic maps of these t-values.
%
%% # Doing the same using the Homer data representation
%
% We can also do this analysis using the Homer data representation. For that we can convert the data from FieldTrip to Homer
%
nirs = fieldtrip2homer(data_combined);

% and then take the data matrix |y| and the experimental design matrix |x|
%
y = nirs.d;
x = nirs.s;

% Subsequently we could continue as above ...

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [gt,relres,iter]=gabfirtight(Lsupport,g,a,M,varargin)
%GABFIRTIGHT  Compute FIR tight window
%   Usage: gt=gabfirtight(Lsupport,g,a,M);
%          gt=gabfirtight(Lsupport,g,a,M, varagin);
%
%   Input parameters:
%     Lsupport  : Length of the tight window
%     g      : Initial window function
%     a      : Time shift
%     M      : Number of Channels
%
%   Output parameters:
%     gt     : Tight window
%
%   `gabfirtight(Lsupport,g,a,M)` computes an FIR window *gd* which is tight.
%   The FIR dual window will be supported on *Lsupport* samples.
%
%   This function solve a convex optimization problem that can be written
%   as:
%
%   .. gd  = argmin_x    || alpha x||_1 +  || beta Fx||_1
%
%   ..                 + || omega (x -g_l) ||_2^2 + delta || x ||_S0
%
%   ..                 + gamma || nabla F x ||_2^2 + mu || nabla x ||_2^2
%
%   ..     such that  x is a tight FIR windows
%
%   .. math:: \begin{split}  \text{gd}  = & \text{arg} \min_x   \| \alpha x \|_1 +  \| \beta \mathcal{F}x\|_1  \\ & + \| \omega (x - g_l) \|_2^2  \\ & \delta \| x \|_{S0}+ \mu \| \nabla x \|_2^2 +\gamma \| \nabla \mathcal{F} x \|_2^2 \\ &    \text{such that } x \text{ is a tight FIR window} \end{split}
%
%   **Note**: This function require the unlocbox. You can download it at
%   `<http://unlocbox.sourceforge.net>`_
%
%   The function uses an iterative algorithm to compute the approximate
%   FIR tight windows. Warning The algorithm solve a non convex problem and
%   might be stack in bad local minima. The algorithm can be controlled by
%   the following flags:
%
%     'alpha',alpha  Weight in time. If it is a scalar, it represent the
%                  weights of the entire L1 function in time. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L1 norm (length: Ldual).
%                  Default value is $\alpha=0$.
%                  **Warning**: this value should not be too big in order to
%                  avoid the the L1 norm proximal operator kill the signal.
%                  No L1-time constraint: $\alpha=0$
%
%     'beta',beta  Weight in frequency. If it is a scalar, it represent the
%                  weights of the entire L1 function in frequency. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L1 norm in frequency. (length: Ldual).
%                  Default value is $\beta=0$.
%                  **Warning**: this value should not be too big in order to
%                  avoid the the L1 norm proximal operator kill the signal.
%                  No L1-frequency constraint: $\beta=0$
%
%     'omega',omega
%                  Weight in time of the L2-norm. If it is a scalar,
%                  it represent the
%                  weights of the entire L2 function in time. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L2 norm (length: Ldual).
%                  Default value is $\omega=0$.
%                  No L2-time constraint: $\omega=0$
%
%     'glike',g_l  $g_l$ is a windows in time. The algorithm try to shape
%                  the dual window like $g_l$. Normalization of $g_l$ is done
%                  automatically. To use option omega should be different
%                  from 0. By default $g_d=0$.
%
%     'mu',mu      Weight of the smooth constraint Default value is 1.
%                  No smooth constraint: $\mu=0$
%
%     'gamma',gamma  Weight of the smooth constraint in frequency.
%                    Default value is 1. No smooth constraint: $\gamma=0$
%
%     'delta',delta  Weight of the S0-norm. Default value is 0.
%                    No S0-norm: $\delta=0$
%
%     'dual'       Look for a dual windows (default)
%
%     'painless'   Construct a starting guess using a painless-case
%                  approximation. This is the default
%
%     'zero'       Choose a starting guess of zero.
%
%     'rand'       Choose a random starting phase.
%
%     'tol',t      Stop if relative residual error is less than the
%                  specified tolerance.
%
%     'maxit',n    Do at most n iterations. default 200
%
%     'print'      Display the progress.
%
%     'debug'      Display all the progresses.
%
%     'quiet'      Don't print anything, this is the default.
%
%     'fast'       Fast algorithm, this is the default.
%
%     'slow'       Safer algorithm, you can try this if the fast algorithm
%                  is not working. Before using this, try to iterate more.
%
%     'printstep',p  If 'print' is specified, then print every p'th
%                    iteration. Default value is p=10;
%
%     'hardconstraint' Force the projection at the end (default)
%
%     'softconstaint' Do not force the projection at the end
%
%   See also: gaboptdual, gabdual, gabtight, gabfirdual, gabconvexopt

% Author: Nathanael Perraudin
% Date  : 18 Feb 2014

if nargin<4
  error('%s: Too few input parameters.',upper(mfilename));
end;

[gt,relres,iter]=gabconvexopt(g,a,M,varargin{:}, 'support',Lsupport,'tight');

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function pwl_interp_2d_scattered_test03 ( prob )

%*****************************************************************************80
%
%% PWL_INTERP_2D_SCATTERED_TEST03 tests PWL_INTERP_2D_SCATTERED_VALUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2012
%
%  Author:
%
%    John Burkardt
%
  ni = 25;
  g = 2;
  nd = g00_size ( g );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PWL_INTERP_2D_SCATTERED_TEST03\n' );
  fprintf ( 1, '  PWL_INTERP_2D_SCATTERED_VALUE evaluates a\n' );
  fprintf ( 1, '  piecewise linear interpolant to scattered data.\n' );
  fprintf ( 1, '  Here, we use grid number %d\n', g );
  fprintf ( 1, '  with %d scattered points in the unit square\n', nd );
  fprintf ( 1, '  on problem %d\n', prob );
%
%  Get the data points and evaluate the function there.
%
  [ xd, yd ] = g00_xy ( g, nd );

  zd = f00_f0 ( prob, nd, xd, yd );

  xyd(1,1:nd) = xd(1:nd);
  xyd(2,1:nd) = yd(1:nd);
%
%  Set up the Delaunay triangulation.
%
  [ element_num, triangle, element_neighbor ] = r8tris2 ( nd, xyd );
%
%  Define the interpolation points.
%
  k = 0;
  for i = 1 : 5
    for j = 1 : 5
      k = k + 1;
      xyi(1,k) = ( 2 * i - 1 ) / 10.0;
      xyi(2,k) = ( 2 * j - 1 ) / 10.0;
    end
  end

  xi(1:ni,1) = xyi(1,1:ni);
  yi(1:ni,1) = xyi(2,1:ni);

  ze = f00_f0 ( prob, ni, xi, yi );
%
%  Evaluate the interpolant.
%
  zi = pwl_interp_2d_scattered_value ( nd, xyd, zd, element_num, ...
    triangle, element_neighbor, ni, xyi );

  rms = 0.0;
  for k = 1 : ni
    rms = rms + ( zi(k) - ze(k) )^2;
  end
  rms = sqrt ( rms / ni );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  RMS error is %g\n', rms );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     K      Xi(K)       Yi(K)       Zi(K)       Z(X,Y)\n' );
  fprintf ( 1, '\n' );

  for k = 1 : ni
    fprintf ( 1, '  %4d  %10f  %10f  %10f  %10f\n', ...
      k, xyi(1,k), xyi(2,k), zi(k), ze(k) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%DEMO_PGAUSS  How to use PGAUSS
%
%   This script illustrates various properties of the Gaussian function.
%
%   .. figure::
%
%      Window+Dual+Tight
%
%      This figure shows an optimally centered Gaussian for a
%      given Gabor system, its canonical dual and tight windows
%      and the DFTs of these windows.
%
%   See also: pgauss

disp('Type "help demo_pgauss" to see a description of how this demo works.');

% A quick test: If the second input parameter to
% pgauss is not specified, the output will be
% invariant under an unitary DFT. Matlabs FFT is does not
% preserve the norm, so it must be scaled a bit.

L=128;
g=pgauss(L);

disp('');
disp('Test of DFT invariance: Should be close to zero.');
norm(g-dft(g))

% Setup parameters and length of signal.
% Note that it must hold that L=M*b=N*a for some integers
% b and N, and that a <= M
L=72;  % Length of signal.
a=6;   % Time shift.
M=9;   % Number of modulations.

% Calculate the frequency shift.
b=L/M;

% For this Gabor system, the optimally concentrated Gaussian
% is given by
g=pgauss(L,a/b);

% This is not invarient with respect to a DFT, but it is still
% real and whole point even
disp('');
disp('The function is WP even. The following should be 1.');
isevenfunction(g)

disp('Therefore, its DFT is real.');
disp('The norm of the imaginary part should be close to zero.');
norm(imag(dft(g)))

% Calculate the canonical dual.
gdual=gabdual(g,a,M);

% Calculate the canonical tight window.
gtight=gabtight(g,a,M);

% Plot them:

% Standard note on plotting:
%
% - All windows have real DFTs, but Matlab does not
%   always recoqnize this, so we have to filter away
%   the small imaginary part by calling REAL(...)
%
% - The windows are all centered around zero, but this
%   is not visually pleasing, so the window must be
%   shifted to the middle by an FFTSHIFT
%

gf_plot      = fftshift(real(dft(g)));
gdual_plot   = fftshift(gdual);
gdualf_plot  = fftshift(real(dft(gdual)));
gtight_plot  = fftshift(gtight);
gtightf_plot = fftshift(real(dft(gtight)));
figure(1);

subplot(3,2,1);
x=(1:L).';
plot(x,fftshift(g),'-',...
     x,circshift(fftshift(g),a),'-',...
     x,circshift(fftshift(g),-a),'-');
title('g=pgauss(72,6/8)');

subplot(3,2,2);
plot(gf_plot);
title('g, frequency domain');

subplot(3,2,3);
plot(gdual_plot);
title('Dual window of g');

subplot(3,2,4);
plot(gdualf_plot);
title('dual window, frequency domain');

subplot(3,2,5);
plot(gtight_plot);
title('Tight window generated from g');

subplot(3,2,6);
plot(gtightf_plot);
title('tight window, frequency domain');


text:
%
% Smooth point-set registration method using neighboring constraints
% -------------------------------------------------------------------
%
% Authors: Gerard Sanromà, René Alquézar and Francesc Serratosa
%
% Contact: gsanorma@gmail.com
% Date: 15/02/2012
%
% Performs cleaning of a matrix with an extra row and column of zeros
% corresponding to the outlier assignments
%
% Input:
%   S: matrix of benefit coefficients (with extra row and column)
%
% Output:
%   Sd: discrete {0,1}-match matrix
%

function Sd = clean_sinkhorn(S)

S = S - min(S(:)); % translate so that min is zero
S(end,end) = 0;
fils = size(S,1)-1;
cols = size(S,2)-1;
Sd = zeros(fils,cols);

while sum(S(:)) > 0
    [v i] = max(S(:));
    [f,c] = ind2sub(size(S),i);
    if f <= fils, S(f,:) = 0; end
    if c <= cols, S(:,c) = 0; end
    if f <= fils && c<= cols, Sd(f,c) = 1; end
end

author: Sable

repo: mcbench-benchmarks


text:
%Coverage
%  Generates XML coverage reports of executed code in 1 format:
%  - Cobertura
%
%  Example usage:
%    report = edu.stanford.covert.test.Coverage(<source dir>, <profile data (optional)>);
%    report.exportCoberturaXML(<fileName>);
%
%  This class was inspired by the mlcovr package by David Hart at Sandia
%  National Labs
%  (http://www.mathworks.com/matlabcentral/fileexchange/24985-mlcovr-package).
%
%  Author: Jonathan Karr
%  Affiliation: Covert Lab, Bioengineering, Stanford University
%  Last Updated: 12/16/2010
classdef Coverage
    properties
        sources
        basePath
        timeStamp
        stats
        packages
    end

    methods
        %- sources is a string indicating a source file directory, or a cell
        %  array containing multiple directories.
        %- basepath is string indicaticating path relative to which fileNames
        %  should be reported
        %- Optionall profData is a struct containing profile data generated by
        %  the MATLAB profiler and constructed by calling
        %    profData = profile('info')
        %  If profData is not provided the current profile data is used.
        function this = Coverage(sources, basePath, profData)
            if ~exist('sources', 'var')
                sources = {'.' 'r'};
            end
            if ~exist('basePath', 'var')
                basePath = '.';
            end
            basePath = absolutepath(basePath);

            if ~iscell(sources)
                sources = {sources 'r'};
            end
            for i = 1:size(sources, 1)
                sources{i, 1} = absolutepath(sources{i,1});

                if ~isdir(sources{i,1})
                    throw(MException('Coverage:invalidSource', 'Sources must be directories'));
                end

                if numel(basePath) > numel(sources{i,1}) || ~strcmp(sources{i,1}(1:numel(basePath)), basePath)
                    throw(MException('Coverage:invalidPath', 'source path must be contained within base path'));
                end
            end

            if ~exist('profData', 'var')
                profData = profile('info');
            end

            this.sources = sources;
            this.basePath = basePath;
            this.timeStamp = datestr(now, 31);
            this = this.buildReportFromDir(profData);
        end
    end

    methods (Access = protected)
        function this = buildReportFromDir(this, profData)
            %% get package names
            packageNames = {};
            for i = 1:size(this.sources, 1)
                packageNames = [packageNames; this.getPackageNamesFromDir(this.sources{i,1})];
            end
            packageNames = unique(packageNames);

            %% build report
            %packages
            packages = repmat(struct('name', [], ...
                'linesCovered', [], 'linesValid', [], 'lineRate', [], ...
                'branchesCovered', [], 'branchesValid',[], 'branchRate', [], ...
                'complexity', [], ...
                'classes', []), [0 1]); %#ok<*PROP>

            this.stats = struct('lineRate', [], 'branchRate', []);
            this.stats.linesCovered = 0;
            this.stats.linesValid = 0;
            this.stats.branchesCovered = 0;
            this.stats.branchesValid = 0;
            this.stats.complexity = 0;

            for i = 1:numel(packageNames)
                package = this.buildPackageReport(meta.package.fromName(packageNames{i}), profData);

                if isempty(package)
                    continue;
                end

                packages = [packages;
                    package];

                this.stats.linesCovered = this.stats.linesCovered + package.linesCovered;
                this.stats.linesValid = this.stats.linesValid + package.linesValid;
                this.stats.branchesCovered = this.stats.branchesCovered + package.branchesCovered;
                this.stats.branchesValid = this.stats.branchesValid + package.branchesValid;
                this.stats.complexity = this.stats.complexity + package.complexity;
            end

            %overall stats
            this.packages = packages;
            this.stats.lineRate = this.stats.linesCovered / this.stats.linesValid;
            this.stats.branchRate = this.stats.branchesCovered / this.stats.branchesValid;
        end

        function package = buildPackageReport(this, metaPackage, profData)
            package = [];

            if isempty(metaPackage.Classes) && isempty(metaPackage.Functions)
                return;
            end

            package.name = metaPackage.Name;
            package.classes = [];
            package.linesCovered = 0;
            package.branchesCovered = 0;
            package.linesValid = 0;
            package.branchesValid = 0;
            package.complexity = 0;
            package.classes = [];

            %classes
            for i = 1:numel(metaPackage.Classes)
                class = this.buildClassReport(metaPackage.Classes{i}, profData);

                if isempty(class); continue; end

                package.classes = [
                    package.classes;
                    class];

                package.linesCovered = package.linesCovered + class.linesCovered;
                package.branchesCovered = package.branchesCovered + class.branchesCovered;
                package.linesValid = package.linesValid + class.linesValid;
                package.branchesValid = package.branchesValid + class.branchesValid;
                package.complexity = package.complexity + class.complexity;
            end

            %package functions
            for i = 1:numel(metaPackage.Functions)
                class = this.buildFunctionReport(metaPackage.Name, metaPackage.Functions{i}, profData);

                if isempty(class); continue; end

                package.classes = [
                    package.classes;
                    class];

                package.linesCovered = package.linesCovered + class.linesCovered;
                package.branchesCovered = package.branchesCovered + class.branchesCovered;
                package.linesValid = package.linesValid + class.linesValid;
                package.branchesValid = package.branchesValid + class.branchesValid;
                package.complexity = package.complexity + class.complexity;
            end

            if isempty(package.classes)
                package = [];
                return;
            end

            %package stats
            package.lineRate = package.linesCovered / package.linesValid;
            package.branchRate = package.branchesCovered / package.branchesValid;
        end

        function class = buildClassReport(this, metaClass, profData)
            class = [];

            %check that metadata is valid
            try
                meta.class.fromName(metaClass.Name);
            catch %#ok<CTCH>
                warning('Coverage:warning', 'Invalid meta data for class %s', metaClass.Name);
            end

            fileName = which(metaClass.Name);
            inSources = false;
            for i = 1:size(this.sources, 1)
                source = this.sources{i, 1};
                if ...
                        (isequal(this.sources{i, 2}, 'r') && numel(fileName) >= numel(source) && strcmp(fileName(1:numel(source)+1), [source filesep])) || ...
                        (~isequal(this.sources{i, 2}, 'r') && numel(fileName) >= numel(source) && strcmp(fileName(1:numel(source)+1), [source filesep]) && ~any(strcmp(fileName(numel(source)+2:end), filesep)))
                    inSources = true;
                    break;
                end
            end

            if ~inSources
                return;
            end

            class.name = metaClass.Name;
            class.fileName = fileName;
            class.linesCovered = 0;
            class.branchesCovered = 0;
            class.linesValid = 0;
            class.branchesValid = 0;
            class.complexity = 0;
            class.methods = [];
            class.lines = [];

            shortClassName = metaClass.Name(find(metaClass.Name == '.', 1, 'last')+1:end);
            methodsCallInfo = getcallinfo(fileName);
            runnableLineIndex = unique(callstats('file_lines', fileName));
            code = {};
            fid = fopen(fileName, 'r');
            while ~feof(fid)
                code{end+1} = fgetl(fid);
            end
            fclose(fid);

            %methods
            metaClass; %#ok<VUNUS> %This line needed, otherwise numel(metaClass.Methods) incorrectly returns 0
            for i = 1:numel(metaClass.Methods)
                metaMethod = metaClass.Methods{i};

                if ~isequal(metaMethod.DefiningClass, metaClass)
                    continue;
                end

                signature = sprintf('[%s] = %s(%s)', ...
                    strjoin(', ', metaMethod.OutputNames{:}),...
                    metaMethod.Name, ...
                    strjoin(', ', metaMethod.InputNames{:}));

                method = this.buildMethodReport(...
                    metaMethod.Name, signature, ...
                    shortClassName, fileName, methodsCallInfo, runnableLineIndex, code, profData);

                if isempty(method); continue; end;

                class.methods = [
                    class.methods;
                    method];

                class.linesCovered = class.linesCovered + sum([method.linesCovered]);
                class.branchesCovered = class.branchesCovered + sum([method.branchesCovered]);
                class.linesValid = class.linesValid + sum([method.linesValid]);
                class.branchesValid = class.branchesValid + sum([method.branchesValid]);
                class.complexity = class.complexity + sum([method.complexity]);
            end

            %getters
            for i = 1:numel(metaClass.Properties)
                if ~isequal(metaClass.Properties{i}.DefiningClass, metaClass)
                    continue;
                end

                method = this.buildMethodReport(...
                    ['get.' metaClass.Properties{i}.Name], ['value = get.' metaClass.Properties{i}.Name '(this)'], ...
                    shortClassName, fileName, methodsCallInfo, runnableLineIndex, code, profData);

                if isempty(method); continue; end;

                class.methods = [
                    class.methods;
                    method];

                class.linesCovered = class.linesCovered + method.linesCovered;
                class.branchesCovered = class.branchesCovered + method.branchesCovered;
                class.linesValid = class.linesValid + method.linesValid;
                class.branchesValid = class.branchesValid + method.branchesValid;
                class.complexity = class.complexity + method.complexity;
            end

            %setters
            for i = 1:numel(metaClass.Properties)
                if ~isequal(metaClass.Properties{i}.DefiningClass, metaClass)
                    continue;
                end

                method = this.buildMethodReport(...
                    ['set.' metaClass.Properties{i}.Name], ['set.' metaClass.Properties{i}.Name '(this, value)'], ...
                    shortClassName, fileName, methodsCallInfo, runnableLineIndex, code, profData);

                if isempty(method); continue; end;

                class.methods = [
                    class.methods;
                    method];

                class.linesCovered = class.linesCovered + method.linesCovered;
                class.branchesCovered = class.branchesCovered + method.branchesCovered;
                class.linesValid = class.linesValid + method.linesValid;
                class.branchesValid = class.branchesValid + method.branchesValid;
                class.complexity = class.complexity + method.complexity;
            end

            %nested functions and subfunctions
            for i = 1:numel(methodsCallInfo)
                methodCallInfo = methodsCallInfo(i);
                if ~ismember(methodCallInfo.type, {'nested-function', 'subfunction'})
                    continue;
                end

                method = this.buildMethodReport(...
                    methodCallInfo.fullname(find(methodCallInfo.fullname == '.',1,'first')+1:end), [], ...
                    shortClassName, fileName, methodsCallInfo, runnableLineIndex, code, profData);

                if isempty(method); continue; end;

                class.methods = [
                    class.methods;
                    method];

                class.linesCovered = class.linesCovered + method.linesCovered;
                class.branchesCovered = class.branchesCovered + method.branchesCovered;
                class.linesValid = class.linesValid + method.linesValid;
                class.branchesValid = class.branchesValid + method.branchesValid;
                class.complexity = class.complexity + method.complexity;
            end

            %lines
            for i = 1:numel(class.methods)
                if ~strcmp(class.methods(i).fileName, fileName); continue; end;

                class.lines = [
                    class.lines;
                    class.methods(i).lines];
            end

            %class stats
            class.lineRate = class.linesCovered / class.linesValid;
            class.branchRate = class.branchesCovered / class.branchesValid;
        end

        function method = buildMethodReport(this, name, signature, shortClassName, fileName, methodsCallInfo, runnableLineIndex, code, profData)
            %initialize output
            method = [];

            methodCallInfo = [
                methodsCallInfo(strcmp({methodsCallInfo.type}, 'class method') & strcmp({methodsCallInfo.name}, name));
                methodsCallInfo(strcmp({methodsCallInfo.type}, 'nested-function') & strcmp({methodsCallInfo.fullname}, [shortClassName '>' shortClassName '.' name]));
                methodsCallInfo(strcmp({methodsCallInfo.type}, 'subfunction') & strcmp({methodsCallInfo.fullname}, [shortClassName '>' name]))];

            %method not defined in main class file
            if isempty(methodCallInfo)
                if ~any(fileName == '@'); return; end;
                methodFileName = [fileName(1:find(fileName == filesep, 1, 'last')) name '.m'];
                if ~exist(methodFileName, 'file'); return; end;

                try
                    methodCallInfo = getcallinfo(methodFileName);
                    runnableLineIndex = unique(callstats('file_lines', methodFileName));
                catch %#ok<CTCH>
                    return;
                end
                code = {};
                fid = fopen(methodFileName, 'r');
                while ~feof(fid)
                    code{end+1} = fgetl(fid);
                end
                fclose(fid);

                method = [];
                for i = 1:numel(methodCallInfo)
                    if i == 1
                        methodProfData = profData.FunctionTable(strcmp({profData.FunctionTable.CompleteName}, [methodFileName '>' name]));
                    else
                        methodProfData = profData.FunctionTable(strcmp({profData.FunctionTable.CompleteName}, [methodFileName '>' name '.' methodCallInfo(i).name]));
                        signature = [];
                    end
                    method = [method;
                        this.buildSubFunctionReport(methodFileName, code, runnableLineIndex, methodCallInfo(i), signature, methodProfData)];
                end
                return;
            end

            %method defined in main class file
            method.name = name;
            method.fileName = fileName;
            method.signature = signature;

            %lines
            methodProfData = profData.FunctionTable(strcmp({profData.FunctionTable.CompleteName}, [fileName '>' shortClassName '.' name]));
            [method.lines, method.linesCovered, method.linesValid, method.lineRate, ...
                method.branchesCovered, method.branchesValid, method.branchRate, method.complexity] = ...
                this.buildLinesReport(code, runnableLineIndex, methodCallInfo, methodProfData);

            %method stats
            if ~isempty(methodProfData)
                method.calls = methodProfData.NumCalls;
                method.time = methodProfData.TotalTime;
            else
                method.calls = 0;
                method.time = 0;
            end
        end

        function class = buildFunctionReport(this, packageName, metaFunction, profData)
            class = [];

            fileName = which([packageName '.' metaFunction.Name]);
            inSources = false;
            for i = 1:size(this.sources, 1)
                source = this.sources{i, 1};
                if ...
                        (isequal(this.sources{i, 2}, 'r') && numel(fileName) >= numel(source) && strcmp(fileName(1:numel(source)+1), [source filesep])) || ...
                        (~isequal(this.sources{i, 2}, 'r') && numel(fileName) >= numel(source) && strcmp(fileName(1:numel(source)+1), [source filesep]) && ~any(strcmp(fileName(numel(source)+2:end), filesep)))
                    inSources = true;
                    break;
                end
            end

            if ~inSources
                return;
            end

            class.name = ['_' metaFunction.Name];
            class.fileName = fileName;
            class.methods = [];
            class.lines = [];
            class.linesCovered = 0;
            class.branchesCovered = 0;
            class.linesValid = 0;
            class.branchesValid = 0;
            class.complexity = 0;

            methodsCallInfo = getcallinfo(fileName);
            runnableLineIndex = unique(callstats('file_lines', fileName));
            code = {};
            fid = fopen(fileName, 'r');
            while ~feof(fid)
                code{end+1} = fgetl(fid);
            end
            fclose(fid);

            for i = 1:numel(methodsCallInfo)
                methodCallInfo = methodsCallInfo(i);
                methodProfData = profData.FunctionTable(strcmp({profData.FunctionTable.CompleteName}, [fileName '>' methodCallInfo.name]));
                method = this.buildSubFunctionReport(fileName, code, runnableLineIndex, methodCallInfo, [], methodProfData);

                class.methods = [
                    class.methods;
                    method];

                class.linesCovered = class.linesCovered + method.linesCovered;
                class.branchesCovered = class.branchesCovered + method.branchesCovered;
                class.linesValid = class.linesValid + method.linesValid;
                class.branchesValid = class.branchesValid +method.branchesValid;
                class.complexity = class.complexity + method.complexity;
            end

            %lines
            class.lines = [];
            for i = 1:numel(class.methods)
                class.lines = [
                    class.lines;
                    class.methods(i).lines];
            end

            %stats
            class.lineRate = class.linesCovered / class.linesValid;
            class.branchRate = class.branchesCovered / class.branchesValid;
        end

        function method = buildSubFunctionReport(this, fileName, code, runnableLineIndex, callInfo, signature, profData)
            method.name = callInfo.name;
            method.fileName = fileName;
            method.signature = signature;

            %lines
            [method.lines, method.linesCovered, method.linesValid, method.lineRate, ...
                method.branchesCovered, method.branchesValid, method.branchRate, method.complexity] = ...
                this.buildLinesReport(code, runnableLineIndex, callInfo, profData);

            %method stats
            if ~isempty(profData)
                method.calls = profData.NumCalls;
                method.time = profData.TotalTime;
            else
                method.calls = 0;
                method.time = 0;
            end
        end

        function [lines, linesCovered, linesValid, lineRate, branchesCovered, branchesValid, branchRate, complexity] = ...
                buildLinesReport(~, code, runLineIndex, callInfo, profData)
            lines = repmat(struct('number', [], 'hits', [], 'branch', [], 'time', [], 'complex', []), ...
                sum(runLineIndex >= callInfo.firstline & runLineIndex <= callInfo.lastline), 1);

            linesCovered = 0;
            branchesCovered = 0;
            linesValid = 0;
            branchesValid = 0;
            complexity = 0;

            l2 = 0;
            for l1 = callInfo.firstline:callInfo.lastline
                if ~any(runLineIndex == l1); continue; end;
                l2 = l2+1;

                hits = 0;
                time = 0;
                if ~isempty(profData)
                    idx = find(profData.ExecutedLines(:, 1) == l1);
                    if ~isempty(idx)
                        hits = profData.ExecutedLines(idx, 2);
                        time = profData.ExecutedLines(idx, 3);
                    end
                end

                branch = ~isempty(regexp(code{l1}, '\<(if|switch|try)\>', 'once', 'tokens'));
                complex = ~isempty(regexp(code{l1},'(\<otherwise\>)|(\<catch\>)|(\<elseif\>)|(\<if\>)|(\<while\>)|(\<for\>)|(\<case\>)|(\<continue\>)|(&&)|(||)','once','tokens'));
                lines(l2, 1).number = l1;
                lines(l2, 1).hits = hits;
                lines(l2, 1).branch = branch;
                lines(l2, 1).time = time;
                lines(l2, 1).complex = complex;

                linesValid = linesValid + 1;
                if branch; branchesValid = branchesValid + 1; end;
                if hits > 0; linesCovered = linesCovered + 1; end;
                if hits > 0 && branch; branchesCovered = branchesCovered + 1; end;
                if complex; complexity = complexity + 1; end;
            end

            lineRate = linesCovered / linesValid;
            branchRate = branchesCovered / branchesCovered;
        end

        function packageNames = getPackageNamesFromDir(this, source)
            if source(end) ~= filesep
                source = [source filesep];
            end
            packageNames = {};

            files = dir(source);
            for i = 1:numel(files)
                if files(i).isdir && files(i).name(1) ~= '.'
                    if files(i).name(1) == '+'
                        metaPackage = meta.package.fromName(files(i).name(2:end));
                        if isempty(metaPackage); continue; end;
                        packageNames = [packageNames;
                            this.getSubPackageNamesFromPackage(metaPackage)];
                    else
                        packageNames = [packageNames;
                            this.getPackageNamesFromDir([source files(i).name])];
                    end
                end
            end
        end

        function packageNames = getSubPackageNamesFromPackage(this, metaPackage)
            packageNames = {metaPackage.Name};

            for i = 1:numel(metaPackage.Packages)
                packageNames = [packageNames; this.getSubPackageNamesFromPackage(metaPackage.Packages{i})];
            end
        end

        function relativePath = relativeToBasePath(this, absolutePath)
            if numel(this.basePath) > numel(absolutePath) || ~strcmp(absolutePath(1:numel(this.basePath)), this.basePath)
                throw(MException('Coverage:invalidPath', 'path must be contained within base path'));
            end
            relativePath = absolutePath(numel(this.basePath)+1:end);
        end
    end

    methods
        %See also http://cobertura.sourceforge.net/xml/coverage-04.dtd
        function xml = exportXML(this, fileName)
            %initialize
            xml = sprintf('<?xml version="1.0" encoding="UTF-8"?>\n');
            xml = [xml sprintf('<!DOCTYPE coverage SYSTEM "http://cobertura.sourceforge.net/xml/coverage-04.dtd">\n')];

            %top-level
            xml = [xml sprintf('<coverage\n')];
            xml = [xml sprintf('\tline-rate="%f"\n', this.stats.lineRate)];
            xml = [xml sprintf('\tbranch-rate="%f"\n', this.stats.branchRate)];
            xml = [xml sprintf('\tlines-covered="%d"\n', this.stats.linesCovered)];
            xml = [xml sprintf('\tlines-valid="%d"\n', this.stats.linesValid)];
            xml = [xml sprintf('\tbranches-covered="%d"\n', this.stats.branchesCovered)];
            xml = [xml sprintf('\tbranches-valid="%d"\n', this.stats.branchesValid)];
            xml = [xml sprintf('\tcomplexity="%f"\n', this.stats.complexity)];
            xml = [xml sprintf('\tversion="%s"\n', '0.1')];
            xml = [xml sprintf('\ttimestamp="%s">\n', this.timeStamp)];
            xml = [xml sprintf('\t<sources>\n')];
            for i = 1:size(this.sources, 1)
                xml = [xml sprintf('\t\t<source><![CDATA[%s]]></source>\n', this.sources{i, 1})];
            end
            xml = [xml sprintf('\t</sources>\n')];
            xml = [xml sprintf('\t<packages>\n')];

            %packages
            for i = 1:numel(this.packages)
                package = this.packages(i);

                xml = [xml sprintf('\t\t<package name="%s" lines-covered="%d" lines-valid="%d" line-rate="%f" branches-covered="%d" branches-valid="%d" branch-rate="%f" complexity="%f">\n', ...
                    package.name, ...
                    package.linesCovered, package.linesValid, package.lineRate, ...
                    package.branchesCovered, package.branchesValid, package.branchRate, ...
                    package.complexity)]; %#ok<*AGROW>
                xml = [xml sprintf('\t\t\t<classes>\n')];
                %classes in package
                for j = 1:numel(package.classes)
                    class = package.classes(j);

                    xml = [xml sprintf('\t\t\t\t<class name="%s" filename="%s" lines-covered="%d" lines-valid="%d" line-rate="%f" branches-covered="%d" branches-valid="%d" branch-rate="%f" complexity="%f">\n', ...
                        class.name, this.relativeToBasePath(class.fileName), ...
                        class.linesCovered, class.linesValid, class.lineRate, ...
                        class.branchesCovered, class.branchesValid, class.branchRate, ...
                        class.complexity)];
                    xml = [xml sprintf('\t\t\t\t\t<methods>\n')];
                    %methods in class
                    for k = 1:numel(class.methods)
                        method = class.methods(k);

                        xml = [xml sprintf('\t\t\t\t\t\t<method name="%s" filename="%s" signature="%s" lines-covered="%d" lines-valid="%d" line-rate="%f" branches-covered="%d" branches-valid="%d" branch-rate="%f" complexity="%f" calls="%d" time="%f">\n', ...
                            method.name, this.relativeToBasePath(method.fileName), method.signature, ...
                            method.linesCovered, method.linesValid, method.lineRate, ...
                            method.branchesCovered, method.branchesValid, method.branchRate, ...
                            method.complexity, ...
                            method.calls, method.time)];
                        xml = [xml sprintf('\t\t\t\t\t\t\t<lines>\n')];
                        for l = 1:numel(method.lines)
                            line = method.lines(l);
                            xml = [xml sprintf('\t\t\t\t\t\t\t\t<line number="%d" hits="%d" branch="%d" time="%f" condition-coverage="0%% (0/0)"/>\n', ...
                                line.number, line.hits, line.branch, line.time)];
                        end
                        xml = [xml sprintf('\t\t\t\t\t\t\t</lines>\n')];
                        xml = [xml sprintf('\t\t\t\t\t\t</method>\n')];
                    end
                    xml = [xml sprintf('\t\t\t\t\t</methods>\n')];

                    %lines in class
                    xml = [xml sprintf('\t\t\t\t\t<lines>\n')];
                    for k = 1:numel(class.lines)
                        line = class.lines(k);
                        xml = [xml sprintf('\t\t\t\t\t\t<line number="%d" hits="%d" branch="%d" time="%f" condition-coverage="0%% (0/0)"/>\n',...
                            line.number, line.hits, line.branch, line.time)];
                    end
                    xml = [xml sprintf('\t\t\t\t\t</lines>\n')];
                    xml = [xml sprintf('\t\t\t\t</class>\n')];
                end
                xml = [xml sprintf('\t\t\t</classes>\n')];
                xml = [xml sprintf('\t\t</package>\n')];
            end
            xml = [xml sprintf('\t</packages>\n')];
            xml = [xml sprintf('</coverage>\n')];

            if exist('fileName', 'var')
                fid = fopen(fileName, 'w');
                fwrite(fid, xml);
                fclose(fid);
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function matlab_figures_to_word_document
%This program converts a Matlab Figure File with a title to Word Document
%through ActiveX Control. The program is designed for a figure with a
%single graph with a title.  It converts the graph title to an array,
%deletes the title, copies the figure to the clipboard, pastes the
%clipboard contents to the Word Document, and writes the title as text
%beneath the figure.  The program can operate in two modes:
%	1.	Single file - Use Lines 12 through 15
%	2.	All files in a directory - Use Lines 16 through 25
%
%The program is currently designed for two lines of text, but more can be
%added by copying and pasting the four lines of code that write text and
%adjusting the subscript numbers.  The program could also be adjusted to
%handle a single-line title by deleting the code that writes the second
%line of the title.  The title in the Word Document will also contain
%today's date, but that can be changed by deleting the appropriate code.
%
%At present, the program only creates new files.

%Operate on just one file
	%[filename pathname] = uigetfile('*.fig', 'Select file to open');
		%full_name = [pathname filename];
	%create_word_figure(full_name);
%Get entire directory
	pathname = uigetdir;
	file_directory = dir([pathname '\*.fig']);
		num_files = length(file_directory);
	%Get individual files
	for k = 3:num_files
		filename = file_directory(k,1).name;
		full_name = [pathname '\' filename];
		create_word_figure(full_name);
	end
%End of main program

function create_word_figure(full_name);
	%Open New Word Document
		Doc = actxserver('Word.Application');
		set(Doc,'Visible',1);
		MS = invoke(Doc.Documents,'Add');
		new = 1;
	%Open file and obtain information to be written to word file
		hndl = open(full_name);								%Open figure and get its handle
		ax_title_hndl = get(gca, 'Title');					%Get handle of title
		Ax_Title_Array = get(ax_title_hndl, 'String');		%Get string of handle
		delete(ax_title_hndl);								%Delete title from figure
		print -dmeta										%Copy figure to clipboard
		close(hndl);										%Close the figure
	%Write to document
		%Put figure in document
			invoke(Doc.Selection, 'Paste');
		%Put first line of title in document
			invoke(Doc.Selection, 'TypeParagraph');
			set(Doc.Selection, 'Text',char(Ax_Title_Array(1)));
			set(Doc.Selection.ParagraphFormat,'Alignment',1);
		%Put second line of title in document
			invoke(Doc.Selection,'MoveDown');
			invoke(Doc.Selection, 'TypeParagraph');
			set(Doc.Selection, 'Text',char(Ax_Title_Array(2)));
			set(Doc.Selection.ParagraphFormat,'Alignment',1);
		%Put today's date in document
			invoke(Doc.Selection,'MoveDown');
			invoke(Doc.Selection, 'TypeParagraph');
			set(Doc.Selection, 'Text',date);
			set(Doc.Selection.ParagraphFormat,'Alignment',1);
	%Save word file
		[fpath,fname,fext] = fileparts(full_name);
		word_name = [fname '.doc'];
		full_name = [fpath '\' word_name];
		invoke(MS,'SaveAs',full_name);
		invoke(Doc,'Quit');
		delete(Doc);
%End of function
%End of program


text:
function area = triangle_area_heron ( s )

%*****************************************************************************80
%
%% TRIANGLE_AREA_HERON computes the area of a triangle using Heron's formula.
%
%  Discussion:
%
%    The formula is valid for any spatial dimension, depending only
%    on the lengths of the sides, and not the coordinates of the vertices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real S(3,1), the lengths of the three sides.
%
%    Output, real AREA, the area of the triangle, or -1.0 if the
%    sides cannot constitute a triangle.
%
  area = (   s(1,1) + s(2,1) + s(3,1) ) ...
       * ( - s(1,1) + s(2,1) + s(3,1) ) ...
       * (   s(1,1) - s(2,1) + s(3,1) ) ...
       * (   s(1,1) + s(2,1) - s(3,1) );

  if ( area < 0.0 )
    area = -1.0;
    return
  end

  area = 0.25 * sqrt ( area );

  return
end


text:
% qpwls_qn_test.m
% Compare Quasi-Newton (QN) and SPS algorithms for QPWLS problem.
% Also shows the benefit of initializing QN with good diagonal Hessian approx.
% Copyright July 2000, Jeff Fessler, The University of Michigan

%
% generate data
%
if ~isvar('yi'), printm 'setup qpwls_qn_test'
	em_wls_test_setup
	W = diag_sp(wi(:));
prompt
end

if ~isvar('R'), printm 'make R'
	f.l2b = -2;
	R = Robject(ig.mask, 'type_denom', 'matlab', 'beta', 2^f.l2b);
prompt
end

if ~isvar('xinit')
	f.niter = 12;
	xinit = xfbp;	% warning: no nonnegativity - QN unconstrained!
end

if ~isvar('P0')
	P0 = G' * (W * sum(G')') + R.denom(R, 0);
	P0 = diag_sp(1 ./ P0(:));
end

%
% PCG
%
if ~isvar('xpcg'), printm 'do pcg'
	xpcg = qpwls_pcg(xinit(ig.mask), G, W, yi, 0, R.C, P0, 1+f.niter);
	xpcg = ig.embed(xpcg);
	im clf, im(xpcg, 'QPWLS-PCG iterations')
prompt
end

%
% test new PCG
%
if ~isvar('xnew')
	if ~isvar('xnew'), printm 'do pcg new'
		xnew = pl_pcg_qs_ls(xinit(ig.mask), G, {yi(:), wi(:)}, ...
			@wls_dercurv, R, 'precon', P0, 'niter', f.niter, ...
			'isave', 'all');
		xnew = ig.embed(xnew);
	end
	max_percent_diff(xpcg, xnew)
prompt
end

%
% QN
%
if ~isvar('xqnp'), printm 'do qnp'
	xqnp = qpwls_qn(xinit(ig.mask), G, W, yi, R.C, P0, 1+f.niter);
	xqnp = ig.embed(xqnp);
	im clf, im(xqnp, 'QPWLS-QNP iterations')
prompt
end
if ~isvar('xqnu'), printm 'do qnu'
	xqnu = qpwls_qn(xinit(ig.mask), G, W, yi, R.C, 1, 1+f.niter);
	xqnu = ig.embed(xqnu);
	im clf, im(xqnu, 'QPWLS-QNU iterations')
prompt
end

%
% SPS
%
if ~isvar('xsps'), printm 'do sps'
	xsps = pwls_sps_os(xinit(ig.mask), yi, wi, G, R, ...
		1+f.niter, [-inf inf], [], [], 1, 0); % disable nonnegativity
	xsps = ig.embed(xsps);
	im clf, im(xsps, 'QPWLS-SPS iterations')
prompt
end

if 1
	im pl 2 2, clim = [0 8];
	im(1, xsps, 'SPS matlab', clim), cbar
	im(2, xqnp, 'QNP matlab', clim), cbar
	im(3, xpcg, 'PCG matlab', clim), cbar
	im(4, xqnu, 'QNU matlab', clim), cbar
%	im(223, (xart2-xmat)/max(xmat(:)), '(aspire-matlab)/max(matlab)')

	cost = @(x) pwls_cost(x, G, W, yi(:), R, ig.mask);
	t1 = cost(xsps);
	t2 = cost(xpcg);
	t3 = cost(xqnp);
	t4 = cost(xqnu);
	if im
		subplot(224)
		plot(	0:f.niter, t1-t1(1), '-o', ...
			0:f.niter, t2-t1(1), '-+', ...
			0:f.niter, t3-t1(1), '-x', ...
			0:f.niter, t4-t1(1), '-v')
		xlabel iteration, ylabel '\Phi change',
		ir_legend({'SPS', 'PCG', 'QNP', 'QNU'})
		title('QPWLS: QN vs SPS vs PCG')
%		axis([0 min(f.niter, 20) -18000 -15000])
		axisy(-5000, 500)
	end

%	ir_savefig 'fig_qpwls_qn_vs_sps'
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function result = part_sf_majorize ( n, nparta, a, npartb, b )

%*****************************************************************************80
%
%% PART_SF_MAJORIZE determines if partition A majorizes partition B.
%
%  Discussion:
%
%    The partitions must be in standard form.
%
%    If A, with NPARTA parts, and B, with NPARTB parts, are both partitions
%    of the same positive integer N, then we say that A majorizes B if,
%    for every index K from 1 to N, it is true that
%
%      sum ( 1 <= I <= K ) B(I) <= sum ( 1 <= I <= K ) A(I)
%
%    where entries of A beyond index NPARTA, and of B beyond BPARTB
%    are assumed to be 0.  We say that A strictly majorizes B if
%    A majorizes B, and for at least one index K the inequality is strict.
%
%    For any two partitions of N, it is possible that A majorizes B,
%    B majorizes A, both partitions majorize each other (in which case
%    they are equal), or that neither majorizes the other.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jack vanLint, Richard Wilson,
%    A Course in Combinatorics,
%    Cambridge, 1992,
%    ISBN: 0-521-42260-4,
%    LC: QA164.L56.
%
%  Parameters:
%
%    Input, integer N, the integer to be partitioned.
%    N must be positive.
%
%    Input, integer NPARTA, the number of parts in partition A.
%    1 <= NPARTA <= N.
%
%    Input, integer A(NPARTA), contains partition A in standard
%    form.  A(1) through A(NPARTA) contain nonzero integers which sum to N.
%
%    Input, integer NPARTB, the number of parts in partition B.
%    1 <= NPARTB <= N.
%
%    Input, integer B(NPARTB), contains partition B in standard
%    form.  B(1) through B(NPARTB) contain nonzero integers which sum to N.
%
%    Output, integer RESULT, the result of the comparison.
%    -2, A and B are incomparable, would have been -1.
%    -1, A < B, (A is strictly majorized by B),
%     0, A = B, (A and B are identical),
%    +1, A > B, (A strictly majorizes B),
%    +2, A and B are incomparable, would have been +1.
%

%
%  Check.
%
  ierror = part_sf_check ( n, nparta, a );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PART_SF_MAJORIZE - Fatal error!\n' );
    fprintf ( 1, '  The input array A is illegal.\n' );
    error ( 'PART_SF_MAJORIZE - Fatal error!' );
  end

  ierror = part_sf_check ( n, npartb, b );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PART_SF_MAJORIZE - Fatal error!\n' );
    fprintf ( 1, '  The input array B is illegal.\n' );
    error ( 'PART_SF_MAJORIZE - Fatal error!' );
  end

  result = 0;
  suma = 0;
  sumb = 0;

  for i = 1 : min ( nparta, npartb )

    if ( i <= nparta )
      suma = suma + a(i);
    end

    if ( i <= npartb )
      sumb = sumb + b(i);
    end

    if ( result == -1 )

      if ( sumb < suma )
        result = -2;
        return
      end

    elseif ( result == 0 )

      if ( suma < sumb )
        result = -1;
      elseif ( sumb < suma )
        result = +1;
      end

    elseif ( result == + 1 )

      if ( suma < sumb )
        result = +2;
        return
      end

    end

  end

  return
end


text:
function display(this)
%DISPLAY   Displays parameters of ALT_BA optical flow algorithm
%
%   This is a member function of the class 'alt_ba_optical_flow'.
%
% Authors: Deqing Sun, Department of Computer Science, Brown University
% Contact: dqsun@cs.brown.edu
% $Date: 2009 $
% $Revision: $
%
% Copyright 2009-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  disp(sprintf('Image sequence (size):          %s', ...
               mat2str(size(this.images))));
  disp(sprintf('Spatial term weighting:         %s', ...
               mat2str(this.lambda)));
  disp(sprintf('Linear solver:                  %s', ...
               this.solver));
  disp(sprintf('# of pyramid levels:            %d', ...
               this.pyramid_levels));
  disp(sprintf('Spacing of pyramid levels:      %d', ...
               this.pyramid_spacing));
  disp(sprintf('Max. # of iterations:           %d', ...
               this.max_iters));
  disp(sprintf('Limit flow updates              %d', ...
               this.limit_update));

author: kristinbranson

repo: JAABA


text:
% EEG_REREJECT - reject the same data portion as the ones defined
%                  in a reference dataset. This function is useful when
%                  have filtered a dataset and want to reject again the
%                  same portions of data.
% Usage:
%   >> OUTEEG = eeg_rereject(EEG, cleanEEG)
%
% Inputs:
%   EEG       - input dataset
%   cleanEEG  - EEGLAB dataset with portion of data rejected (information
%               stored in boundary event types).
%
% Optional inputs:
%   'blank'    - ['on'|'off'] replace by 0 instead of removing the data
%                portions. Default is 'off'. Note that this does not
%                add any boundary events.
%   'rmdiscont' - ['on'|'off'] remove discontinuities (by shifting the
%                data in the different discontinuous portions). Note that
%                this discards the boundary events.
%
% Outputs:
%   OUTEEG     - output dataset
%
% Author: Arnaud Delorme, 2020
%
% Note: uses the RESAMPLE function from the signal processing toolbox
%       if present. Otherwise use griddata interpolation method (it should be
%       reprogrammed using spline interpolation for speed up).
%
% See also: RESAMPLE, EEGLAB

% Copyright (C) 2020 Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [EEG, boundloc]= eeg_rereject(EEG, CLEANEDEEG, varargin)

boundLat = [];
if nargin < 2
    help eeg_rereject;
    return;
end

% check inputs
% ------------
g = finputcheck(varargin, { 'blank'      'string'  {'on','off' }             'off';
                            'rmdiscont'  'string'  {'on','off' }             'off' });

if isstr(g), error(g); end

boundEvents = eeg_findboundaries(CLEANEDEEG);
boundloc = [ CLEANEDEEG.event(boundEvents).latency ];
dur      = [ CLEANEDEEG.event(boundEvents).duration ];
cumdur   = cumsum(dur);
boundloc = boundloc + [0 cumdur(1:end-1) ];
boundloc = [ boundloc; boundloc+dur-1]';
if strcmpi(g.blank, 'off')
    EEG = eeg_eegrej(EEG, ceil(boundloc));
    boundEvents = eeg_findboundaries(EEG.event);
    boundLat    = [ EEG.event(boundEvents).latency ];
else
    boundloc = round(boundloc);
    for iBound = size(boundloc,1):-1:1
        EEG.data(:,boundloc(iBound,1):boundloc(iBound,2)) = 0;
    end
    boundEvents = [];
    boundlocTmp = boundloc;
    boundlocTmp(:,1) = boundlocTmp(:,1)-1; % to compensate round(lat+0.5) below
    boundLat = sort(boundlocTmp(:))';
end

% remake data continuous
if strcmpi(g.rmdiscont, 'on')
    for ind = length(boundLat):-1:1
        if boundLat(ind) > 1
            lat = min(round(boundLat(ind)+0.5), size(EEG.data,2));
            EEG.data(:,lat:end) = EEG.data(:,lat:end) - EEG.data(:,lat) + EEG.data(:,lat-1);
            EEG.data(:,lat) = []; % remove 1 sample
            if ~isempty(boundEvents)
                for ind2 = boundEvents(ind):length(EEG.event)
                    EEG.event(ind2).latency = EEG.event(ind2).latency-1; % shift events by 1 sample
                end
            end
        end
    end
    EEG.event(boundEvents) = [];
    EEG.pnts = size(EEG.data,2);
end
EEG = eeg_checkset(EEG);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% JN Kather, NCT Heidelberg / RWTH Aachen, 2017-2018
% see separate LICENSE
%
% This MATLAB script is associated with the following project
% "Deep learning can predict microsatellite instability directly
% from histology in gastrointestinal cancer". Please consider citing this
% publication if you re-use the code
%
% Step 05:
% this is an important step: we will split all image tiles into train and
% test set. This is done on a patient level, that means each patient will
% be assigned to the train or test group and then all tiles of this patient
% will be moved to the train or test folder

clear variables, close all, clc
addpath(genpath('./subroutines/'));
sq = @(varargin) varargin';
doCopy = true; % if false, the images will not be copied (simulation run)

% define the path to the normalized tiles
cnst.tileInputPath = 'D:\Colorectal_MSI_Project\DX_TILES_NORM\'; % source folder
cnst.tileOutputPath = 'E:\DX_TILES_NORM_SPLIT\'; % where to save the train and test set
mkdir(cnst.tileOutputPath);

% define the path to the clinical data table and to the hypermutation data table
% this requires Microsoft Excel to be installed on the machine
cnst.cliniDataTable = './cliniData/merged_TCGA_TUM_clini_table_v1.xlsx';

% read all tiles, might take a while
tic,
allTiles = dir([cnst.tileInputPath,'*.png']);
toc
tic, allTileNames = sq(allTiles.name); toc
disp('successfully read all tile names');

% read clinical data (merged data)
cliniDat = readtable(cnst.cliniDataTable);
cliniPatient = cliniDat.submitter_id;
cliniTumor = cliniDat.project_id;
cliniMSI = cliniDat.MSIStatus;
cliniHiMut = cliniDat.Tota_Mutation_cnt; % extract hypermutated samples.
% hypermutated tumors will have a number here, all other will have NaN

%% MATCH IMAGE PATCHES WITH PATIENT LEVEL DATA FOR TABLE 1
matches = zeros(size(cliniPatient));
tic
for i = 1:numel(cliniPatient)
    currPat = cliniPatient{i}; % current patient name
    TF = contains(allTileNames,currPat); % match all tiles
    matches(i) = sum(TF); % how many tiles were matched to this patient
    if isempty(cliniMSI{i}) % if no MSI info is present, look for mutation data
        % the reason for this is that we have few MSI patients and we do
        % not want to miss any because the MSI field is empty. This will
        % affect only very few patients (8 STAD and 6 COAD patients)
        if isnan(cliniHiMut(i)) % set status not hypermutated
            msistat{i} = 'NA'; % no MSI data, not hypermutated
        else % set status hypermutated
            msistat{i} = 'MUT'; % no MSI data, but hypermutated (highly probable MSI)
        end
    else % set status MSS or MSI
        msistat{i} = cliniMSI{i};
    end
    if mod(i,25)==1
    disp(['progress: ',num2str(round(i/numel(cliniPatient)*100,1))]);
    end
end
toc

% how many MSI and MSS patients could be matched
imgmatch = matches>0;
msimatch = strcmp(msistat,'MSI-H')'&imgmatch; % these are known MSI-H patients
mutmatch = strcmp(msistat,'MUT')'&imgmatch;   % these patients have no MSI status but are hypermutated
mssmatch = strcmp(msistat,'MSS')'&imgmatch;   % these patients are MSS

disp(['I matched ',num2str(sum(msimatch)),' MSI-H patients']);
disp(['I matched ',num2str(sum(mutmatch)),' non-MSI-H MUT patients']);
disp(['I matched ',num2str(sum(mssmatch)),' MSS patients']);

%% RANDOMIZE PATIENTS
rng(1); % reset the random number generator for reproducibility
figure(1),clf, hold on
traintest = 0.70; % split the data set in a 70/30 ratio (approximate)

msitrain = msimatch & rand(size(msimatch))<traintest;
muttrain = mutmatch & rand(size(mutmatch))<traintest;
msstrain = mssmatch & rand(size(mssmatch))<traintest;

plot(1:numel(imgmatch),msitrain,'r.'),plot(1:numel(imgmatch),muttrain,'b.'),plot(1:numel(imgmatch),msstrain,'g.')

disp(['Training cohort contains ',num2str(sum(msitrain)),' MSI-H patients']);
disp(['Training cohort contains ',num2str(sum(muttrain)),' non-MSI-H MUT patients']);
disp(['Training cohort contains ',num2str(sum(msstrain)),' MSS patients']);

% everybody not in training set goes to test set
msitest = msimatch & ~msitrain;
muttest = mutmatch & ~muttrain;
msstest = mssmatch & ~msstrain;

plot(1:numel(imgmatch),msitest,'rx'),plot(1:numel(imgmatch),muttest,'bx'),plot(1:numel(imgmatch),msstest,'gx')

disp(['Test cohort contains ',num2str(sum(msitest)),' MSI-H patients']);
disp(['Test cohort contains ',num2str(sum(muttest)),' non-MSI-H MUT patients']);
disp(['Test cohort contains ',num2str(sum(msstest)),' MSS patients']);

%% COPY IMAGE TILES TO TRAIN AND TEST FOLDERS

% now copy the files to the actual folder. This might take a while.
if doCopy
moveImagesToFolder(cliniPatient,msitrain,cnst.tileInputPath,cnst.tileOutputPath,'msitrain');
moveImagesToFolder(cliniPatient,muttrain,cnst.tileInputPath,cnst.tileOutputPath,'muttrain');
moveImagesToFolder(cliniPatient,msstrain,cnst.tileInputPath,cnst.tileOutputPath,'msstrain');

moveImagesToFolder(cliniPatient,msitest,cnst.tileInputPath,cnst.tileOutputPath,'msitest');
moveImagesToFolder(cliniPatient,muttest,cnst.tileInputPath,cnst.tileOutputPath,'muttest');
moveImagesToFolder(cliniPatient,msstest,cnst.tileInputPath,cnst.tileOutputPath,'msstest');
end

author: jnkather

repo: MSIfromHE

save_path: github-repos/MATLAB/jnkather-MSIfromHE