text:
function DeltaR=solidTideShift(rStation,rSun,rMoon,Jul1,Jul2,addPermTide)
%%SOLIDTIDESHIFT Compute the vector offset of a point on the surface of the
%             Earth due to solid Earth tidal effects. These effects can be
%             over 31cm. If a Julian date is given, third-order solid Earth
%             tidal effects will be taken into account. Note that site
%             displacement due to ocean loading, which is not included in
%             the solid Earth tide offset, can also be quite significant,
%             with the IERS Conventions 2010 listing it up to 10cm.
%
%INPUTS: rStation A vector from the geocenter to the (non-offset) location
%                 of a station on the surface of the Earth in WGS-84 Earth-
%                 centered Earth-fixed (ECEF) coordinates in meters. This
%                 value is converted into spherical coordinates by this
%                 function, and only the angles in spherical coordinates,
%                 not the radius, are important for this algorithm.
%                 However, since locations on the ground are generally
%                 given in WGS-84 ellipsoidal coordinates, the ellipsoidal
%                 height of the ground plays a role in the conversion to
%                 spherical coordinates. However, variations in the
%                 ellipsoidal height provided on the order of 2km result in
%                 differences in the output DeltaR on the order of microns,
%                 so in practice, one could just set rStation to the
%                 Cartesian position corresponding to a particular
%                 ellipsoidal latitude and longitude with zero ellipsoidal
%                 height. Ideally, however, rStation is the Cartesian
%                 location of the land in a mean-tide model (if
%                 addPermTide=false), or in a tide-free model if
%                 addPermTide=true.
%            rSun A vector from the geocenter to the sun in ITRS
%                 coordinates in meters. This can be obtained using the
%                 readJPLEphem and GCRS2ITRS functions.
%           rMoon A vector from the geocenter to the moon in ITRS
%                 coordinates in meters. This can be obtained using the
%                 readJPLEphem and GCRS2ITRS functions.
%      Jul1, Jul2 Two parts of a Julian date given in terrestrial time
%                 (TT). The units of the date are days. The full date is
%                 the sum of both terms. The date is broken into two
%                 parts to provide more bits of precision. It does not
%                 matter how the date is split. These parameters must be
%                 given for the third-order tidal components to be taken
%                 into account.
%     addPermTide The deltaR values are for a conventional tide-free model.
%                 If addPermTide is true, then a constant (depending on
%                 latitude) permanent tide offset will be added to put the
%                 point into a mean-tide model. If this parameter is
%                 omitted, the default value is "false".
%
%OUTPUTS: DeltaR The offset of a point on the surface of the Earth due to
%                solid Earth tides. rStation+DeltaR is the location of the
%                station (point on the ground) at the terrestial time given
%                by Jul1 and Jul2 taking into account solid Earth tides.
%                Note that the addPermTide term must be consistent with the
%                coordinate system of rStation. If rStation already
%                includes the permanent tides (is in a mean-tide model),
%                then addPermTide should be false. Otherwise, addPermTide
%                should be true.
%
%The Sun and Moon cause the crust of the Earth to warp and points on the
%Earth in WGS-84 coordinates to move over time. The formulae for computing
%tidal shits of the crust are given in Section 7.1 of [1].
%
%Section 7.1.1 discusses solid Earth tides. While the second-order
%corrections are well documneted, the third order corrections can not be
%implemented directly from the standard. However, the standard specifies
%FORTRAN routines, whose algorithms do not resemble the work given in the
%standard. These Fortran routines, have been converted to Matlab and are
%called for the third-order corrections if Julian dates are given. The
%converted routines are separate functions at the bottom of this file.
%
%Note that third-order model has an implied set of ephemerides built into
%it that might not be perfectly consistent with whatever model was used to
%obtain rSun and rMoon.
%
%One time component is supposed to be in Julian centuries since J2000.0. As
%J2000.0 is defined in TDB and the time is given in TT, the TDB2TT function
%is used. The extra parameters needed for the TDB2TT function are not
%requested as inputs to this function, since it is assumed that the
%difference between the tides is too small to matter for the model.
%
%To test this function, one can use the same values that are given in the
%IERS's implementation, DEHANTTIDEINEL.F Those are
% rStation = [4075578.385;931852.890;4801570.154];
% rSun=[137859926952.015;54228127881.4350;23509422341.6960];
% rMoon=[-179996231.920342;-312468450.131567;-169288918.592160];
% addPermTide=false;
% [Jul1,Jul2]=Cal2TT(2009,4,13,0,0,0);
% DeltaR=solidTideShift(rStation,rSun,rMoon,addPermTide,Jul1,Jul2);
%
%The value of DeltaR truncated to 9 places should be about
%DeltaR = 0.0770042035
%         0.0630405632
%         0.0551656815
%
%REFERENCES:
%[1] G. Petit and B. Luzum, IERS Conventions (2010), International Earth
%    Rotation and Reference Systems Service Std. 36, 2010.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The convention for indexing arrays is that 1 refers to a value with
%respect to lunar parameters and 2 refers to a value with respect to solar
%parameters.

if(nargin<6||isempty(addPermTide))
    addPermTide=false;
end

%%%%%%General constants that are used%%%%%%
Re=Constants.EarthEqRadius;%meters
GMGMe(1)=Constants.MoonEarthMassRatio;
GMGMe(2)=Constants.SunEarthMassRatio;

%%%%%%Convert the input parameters into different coordinate systems%%%%%%

%Obtain the spherical coordinates of the station.
pointSpherStation=Cart2Sphere(rStation);
lambda=pointSpherStation(2);%Geocentric longitude
phi=pointSpherStation(3);%Geocentric latitude.
r=pointSpherStation(1);
rHat=rStation/r;%A unit vector in the direction of the station.

%Get East-North-Up Axes using a local SPHERICAL Earth model. A flattening
%factor of zero makes the Earth flat.
uENU=getENUAxes([phi;lambda;0],false,r,0);
uE=uENU(:,1);%Spherical East unit vector
uN=uENU(:,2);%Spherical North unit vector.
uU=uENU(:,3);%Spherical radial *up) unit vector

%Extract the x,y, and z coordinates of the moon.
Xj(1)=rMoon(1);
Yj(1)=rMoon(2);
Zj(1)=rMoon(3);
%Obtain the spherical coordinates of the moon.
pointSpherMoon=Cart2Sphere(rMoon);
Lambdaj(1)=pointSpherMoon(2);
Phij(1)=pointSpherMoon(3);
Rj(1)=pointSpherMoon(1);
RjHat(:,1)=rMoon/Rj(1);%A unit vector in the direction of the moon.

%Extract the x,y, and z coordinates of the sun.
Xj(2)=rSun(1);
Yj(2)=rSun(2);
Zj(2)=rSun(3);
%Obtain the spherical coordinates of the sun.
pointSpherSun=Cart2Sphere(rSun);
Lambdaj(2)=pointSpherSun(2);
Phij(2)=pointSpherSun(3);
Rj(2)=pointSpherSun(1);
RjHat(:,2)=rSun/Rj(2);%A unit vector in the direction of the sun.

%%%%%%Determine the solid Earth tidal displacement.%%%%%%
%This is a three step process. The first two steps are described in the
%table on page 103 of the 2010 IERS conventions.

%%%Step 1: time-domain contributions
%%First, the in-phase contribution.

%Determine the second degree nomial Love (h) and Shida (l) numbers using
%the latitude correction in the second paragraph on page 105.
h0=0.6078;
h2=-0.0006;
l0=0.0847;
l2=0.0002;

h2=h0+h2*(3*sin(phi)^2-1)/2;
l2=l0+l2*(3*sin(phi)^2-1)/2;

%Equation 7.5 for the degree 2 tidal offset
DeltaR=0;
for n=1:2
    RjHatDotrHat=dot(RjHat(:,n),rHat);
    DeltaR=DeltaR+GMGMe(n)*(Re^4/Rj(n)^3)*(h2*rHat*(3*RjHatDotrHat^2-1)/2+3*l2*RjHatDotrHat*(RjHat(:,n)-RjHatDotrHat*rHat));
end

%Equation 7.6 for the degree 3 tidal offset
h3=0.292;
l3=0.015;
for n=1:2
    RjHatDotrHat=dot(RjHat(:,n),rHat);
    DeltaR=DeltaR+GMGMe(n)*(Re^5/Rj(n)^4)*(h3*rHat*((5/2)*RjHatDotrHat^3-(3/2)*RjHatDotrHat)+l3*((15/2)*RjHatDotrHat^2-3/2)*(RjHat(:,n)-RjHatDotrHat*rHat));
end

%Equation 7.8 for the diurnal latitude dependence of the Love numbers.
l1Diurn=0.0012;
deltaT=0;
for n=1:2
    P12Cos=3*Xj(n)*Zj(n)/Rj(n)^2;%Equation 7.7b
    P12Sin=3*Yj(n)*Zj(n)/Rj(n)^2;%Equation 7.7b

    %This is a term in the sum of Equation 7.8. However, it does not look
    %like Equation 7.8, since angle difference formulae had to be used to
    %get terms suitable for the values in Equation 7.7b to be used.
    %Specifically,
    %sin(lambda-Lambdaj(n))=sin(lambda)*cos(Lambdaj(n))-cos(lambda)*sin(Lambdaj(n));
    %cos(lambda-Lambdaj(n))=cos(lambda)*cos(Lambdaj(n))+sin(lambda)*sin(Lambdaj(n));
    deltaT=deltaT+GMGMe(n)*(Re^4/Rj(n)^3)*(uN*sin(phi)*(cos(lambda)*P12Cos+sin(lambda)*P12Sin)-uE*cos(2*phi)*(sin(lambda)*P12Cos-cos(lambda)*P12Sin));
end
deltaT=-l1Diurn*sin(phi)*deltaT;

DeltaR=DeltaR+deltaT;

%Equation 7.9 for the semi-diurnal latitude dependence of the Love numbers.
l1SemiDiurn=0.0024;
deltaT=0;
for n=1:2
    P22Cos=(3/Rj(n)^2)*(Xj(n)^2-Yj(n)^2);%Equation 7.7c
    P22Sin=(6/Rj(n)^2)*Xj(n)*Yj(n);%Equation 7.7c

    %This is a term in the sum of Equation 7.9. However, it does not look
    %like Equation 7.9, since angle difference formulae had to be used to
    %get terms suitable for the values in Equation 7.7c to be used.
    %Specifically,
    %sin(2*(lambda-Lambdaj(n)))=sin(2*lambda)*cos(2*Lambdaj(n))-cos(2*lambda)*sin(2*Lambdaj(n));
    %cos(2*(lambda-Lambdaj(n)))=cos(2*lambda)*cos(2*Lambdaj(n))+sin(2*lambda)*sin(2*Lambdaj(n));
    deltaT=deltaT+GMGMe(n)*(Re^4/Rj(n)^3)*(uN*(cos(2*lambda)*P22Cos+sin(2*lambda)*P22Sin)+uE*sin(phi)*(sin(2*lambda)*P22Cos-cos(2*lambda)*P22Sin));
end
deltaT=-(1/2)*l1SemiDiurn*sin(phi)*cos(phi)*deltaT;

DeltaR=DeltaR+deltaT;

%%Next, the out-of-phase contributions for the second degree terms

%First, the diurnal tides in Equations 7.10a and 7.10b
hIDiurn=-0.0025;
lIDirun=-0.0007;

deltaR=0;
deltaT=0;
for n=1:2
    deltaR=deltaR+GMGMe(n)*(Re^4/Rj(n)^3)*sin(2*Phij(n))*sin(2*phi)*sin(lambda-Lambdaj(n));
    deltaT=deltaT+GMGMe(n)*(Re^4/Rj(n)^3)*sin(2*Phij(n))*(cos(2*phi)*sin(lambda-Lambdaj(n))*uN+sin(phi)*cos(lambda-Lambdaj(n))*uE);
end
deltaR=-(3/4)*hIDiurn*deltaR;
deltaT=-(3/2)*lIDirun*deltaT;

DeltaR=DeltaR+deltaR*uU+deltaT;

%Next, the semidiurnal tides in Equations 7.11a and 7.11b.
hISemiDiurn=-0.0022;
lISemiDiurn=-0.0007;

deltaR=0;
deltaT=0;
for n=1:2
    deltaR=deltaR+GMGMe(n)*(Re^4/Rj(n)^3)*cos(Phij(n))^2*cos(phi)^2*sin(2*(lambda-Lambdaj(n)));
    deltaT=deltaT+GMGMe(n)*(Re^4/Rj(n)^3)*cos(Phij(n))^2*(sin(2*phi)*sin(2*(lambda-Lambdaj(n)))*uN-2*cos(phi)*cos(2*(lambda-Lambdaj(n)))*uE);
end
deltaR=-(3/4)*hISemiDiurn*deltaR;
deltaT=(3/4)*lISemiDiurn*deltaT;

DeltaR=DeltaR+deltaR*uU+deltaT;

%%%Step 2: frequency-domain contributions

%Only apply the corrections if a Julian date has been provided.
if(nargin>4)
    %%First, the contributions for the diurnal band
      %Convert Terrestrial Time to Julian centuries since J2000.0. A Julian
      %century is defined to have 36525 days in it. The offset of 2451545.0
      %days is the Julian date at J2000.0 in TDB. In TT, a difference of a
      %few milliseconds exists.
      [TDB1,TDB2]=TT2TDB(Jul1,Jul2);
      T=((TDB1-2451545.0)+TDB2)/36525;

      %Convert the Julian date to UTC and determine the fractional number
      %of hours passed in the day assuming precisely 24 hours in a say.
      [Jul1,Jul2]=TT2UTC(Jul1,Jul2);
      [~,~,~,dayFrac]=UTC2Cal(Jul1,Jul2,true);
      FHR=24*dayFrac;

    DeltaR=DeltaR+diurBandCorr(rStation,FHR,T);

    %%Next, the contributions from the long-period band
    DeltaR=DeltaR+longBandCorr(rStation,T);
end

if(addPermTide~=false)
    P2=(3*sin(phi)^2-1)/2;
    %Equation 7.14a in the IERS conventions
    deltaR=(-0.1206+0.0001*P2)*P2;
    %Equation 7.14b in the IERS conventions
    deltaT=(-0.0252-0.0001*P2)*sin(2*phi);
    DeltaR=DeltaR+deltaR*uU+deltaT*uN;
end
end

function XCORSTA=diurBandCorr(XSTA,FHR,T)
%%DIURBANDCORR  This subroutine is a Matlab translation of the subroutine
%               STEP2DIU that is provided by the IERS at
%               ftp://tai.bipm.org/iers/convupdt/chapter7/
%               The algorithm is not fully documented by the IERS 2010
%               conventions. For example, most of the constants are not
%               listed anywhere in the documentation.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%
%SUBROUTINE STEP2DIU (XSTA,FHR,T,XCORSTA)
%
%  - - - - - - - - - - -
%   S T E P 2 D I U
%  - - - - - - - - - - -
%
%  This routine is part of the International Earth Rotation and
%  Reference Systems Service (IERS) Conventions software collection.
%
%  This subroutine gives the in-phase and out-of-phase corrections
%  induced by mantle anelasticity in the diurnal band.
%
%  In general, Class 1, 2, and 3 models represent physical effects that
%  act on geodetic parameters while canonical models provide lower-level
%  representations or basic computations that are used by Class 1, 2, or
%  3 models.
%
%  Status: Class 1
%
%     Class 1 models are those recommended to be used a priori in the
%     reduction of raw space geodetic data in order to determine
%     geodetic parameter estimates.
%     Class 2 models are those that eliminate an observational
%     singularity and are purely conventional in nature.
%     Class 3 models are those that are not required as either Class
%     1 or 2.
%     Canonical models are accepted as is and cannot be classified as a
%     Class 1, 2, or 3 model.
%
%  Given:
%     XSTA          d(3)   Geocentric position of the IGS station (Note 1)
%     FHR           d      Fractional hours in the day (Note 2)
%     T             d      Centuries since J2000
%
%  Returned:
%     XCORSTA       d(3)   In phase and out of phase station corrections
%                          for diurnal band (Note 4)
%
%  Notes:
%
%  1) The IGS station is in ITRF co-rotating frame.  All coordinates are
%     expressed in meters.
%
%  2) The fractional hours in the day is computed as the hour + minutes/60.0
%     + sec/3600.0.  The unit is expressed in Universal Time (UT).
%
%  4) All coordinates are expressed in meters.
%
%  Test case:
%     given input: XSTA(1) = 4075578.385D0 meters
%                  XSTA(2) =  931852.890D0 meters
%                  XSTA(3) = 4801570.154D0 meters
%                  FHR     = 0.00D0 hours
%                  T       = 0.1059411362080767D0 Julian centuries
%
%     expected output:  XCORSTA(1) = 0.4193085327321284701D-02 meters
%                       XCORSTA(2) = 0.1456681241014607395D-02 meters
%                       XCORSTA(3) = 0.5123366597450316508D-02 meters
%
%  References:
%
%     Mathews, P. M., Dehant, V., and Gipson, J. M., 1997, ''Tidal station
%     displacements," J. Geophys. Res., 102(B9), pp. 20,469-20,477
%
%     Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
%     IERS Technical Note No. 36, BKG (2010)
%
%  Revisions:
%  1996 March    23 V. Dehant      Original code
%  2009 July     31 B.E. Stetzler  Initial standardization of code
%  2009 August   06 B.E. Stetzler  Provided a test case
%  2009 August   06 B.E. Stetzler  Capitalized all variables for
%                                  Fortran 77 compatibility
%  2010 October  20 B.E. Stetzler  Input T corrected to be number of
%                                  centuries since J2000
%-----------------------------------------------------------------------

      D2PI = 6.283185307179586476925287;

      DATDI=[-3,    0,  2,   0,  0, -0.01,   0,      0,      0;
             -3,    2,  0,   0,  0, -0.01,   0,      0,      0;
             -2,    0,  1,  -1,  0, -0.02,   0,      0,      0;
             -2,    0,  1,   0,  0, -0.08,   0,     -0.01,   0.01;
             -2,    2, -1,   0,  0, -0.02,   0,      0,      0;
             -1,    0,  0,  -1,  0, -0.10,   0,      0,      0;
             -1,    0,  0,   0,  0, -0.51,   0,     -0.02,   0.03;
             -1,    2,  0,   0,  0,  0.01,   0,      0,      0;
              0,   -2,  1,   0,  0,  0.01,   0,      0,      0;
              0,    0, -1,   0,  0,  0.02,   0,      0,      0;
              0,    0,  1,   0,  0,  0.06,   0,      0,      0;
              0,    0,  1,   1,  0,  0.01,   0,      0,      0;
              0,    2, -1,   0,  0,  0.01,   0,      0,      0;
              1,   -3,  0,   0,  1, -0.06,   0,      0,      0;
              1,   -2,  0,  -1,  0,  0.01,   0,      0,      0;
              1,   -2,  0,   0,  0, -1.23,  -0.07,   0.06,   0.01;
              1,   -1,  0,   0, -1,  0.02,   0,      0,      0;
              1,   -1,  0,   0,  1,  0.04,   0,      0,      0;
              1,    0,  0,  -1,  0, -0.22,   0.01,   0.01,   0;
              1,    0,  0,   0,  0, 12.00,  -0.80,  -0.67,  -0.03;
              1,    0,  0,   1,  0,  1.73,  -0.12,  -0.10,   0;
              1,    0,  0,   2,  0, -0.04,   0,      0,      0;
              1,    1,  0,   0, -1, -0.50,  -0.01,   0.03,   0;
              1,    1,  0,   0,  1,  0.01,   0,      0,      0;
              0,    1,  0,   1, -1, -0.01,   0,      0,      0;
              1,    2, -2,   0,  0, -0.01,   0,      0,      0;
              1,    2,  0,   0,  0, -0.11,   0.01,   0.01,   0;
              2,   -2,  1,   0,  0, -0.01,   0,      0,      0;
              2,    0, -1,   0,  0, -0.02,   0,      0,      0;
              3,    0,  0,   0,  0,  0,      0,      0,      0;
              3,    0,  0,   1,  0,  0,      0,      0,      0];

    DEG2RAD = D2PI/360;

%  Compute the phase angles in degrees.
    S = 218.31664563+(481267.88194+(-0.0014663889+(0.00000185139)*T)*T)*T;

    TAU = FHR*15+280.4606184+(36000.7700536+(0.00038793+(-0.0000000258)*T)*T)*T+(-S);

    PR = (1.396971278+(0.000308889+(0.000000021+(0.000000007)*T)*T)*T)*T;

    S = S + PR;

    H = 280.46645+(36000.7697489+(0.00030322222+(0.000000020+(-0.00000000654)*T)*T)*T)*T;

    P = 83.35324312+(4069.01363525+(-0.01032172222+(-0.0000124991+(0.00000005263)*T)*T)*T)*T;

    ZNS = 234.95544499+(1934.13626197+(-0.00207561111+(-0.00000213944+(0.00000001650)*T)*T)*T)*T;

    PS = 282.93734098+(1.71945766667+(0.00045688889+(-0.00000001778+(-0.00000000334)*T)*T)*T)*T;

% Reduce angles to between the range 0 and 360.
    S =  mod(S,360);
    TAU = mod(TAU,360);
    H =  mod(H,360);
    P =  mod(P,360);
    ZNS = mod(ZNS,360);
    PS = mod(PS,360);

    RSTA = norm(XSTA);
    SINPHI = XSTA(3)/RSTA;
    COSPHI = norm(XSTA(1:2))/RSTA;

    COSLA = XSTA(1)/COSPHI/RSTA;
    SINLA = XSTA(2)/COSPHI/RSTA;
    ZLA = atan2(XSTA(2),XSTA(1));

% Initialize.
    XCORSTA=zeros(3,1);

    for J=1:31
        % Convert from degrees to radians.
        THETAF=(TAU+DATDI(J,1)*S+DATDI(J,2)*H+DATDI(J,3)*P+DATDI(J,4)*ZNS+DATDI(J,5)*PS)*DEG2RAD;

        DR=DATDI(J,6)*2*SINPHI*COSPHI*sin(THETAF+ZLA)+DATDI(J,7)*2*SINPHI*COSPHI*cos(THETAF+ZLA);

        DN=DATDI(J,8)*(COSPHI^2-SINPHI^2)*sin(THETAF+ZLA)+DATDI(J,9)*(COSPHI^2-SINPHI^2)*cos(THETAF+ZLA);
        %      DE=DATDI(8,J)*SINPHI*COS(THETAF+ZLA)+
        %     Modified 20 June 2007

        DE=DATDI(J,8)*SINPHI*cos(THETAF+ZLA)-DATDI(J,9)*SINPHI*sin(THETAF+ZLA);

        XCORSTA(1)=XCORSTA(1)+DR*COSLA*COSPHI-DE*SINLA-DN*SINPHI*COSLA;
        XCORSTA(2)=XCORSTA(2)+DR*SINLA*COSPHI+DE*COSLA-DN*SINPHI*SINLA;
        XCORSTA(3)=XCORSTA(3)+DR*SINPHI+DN*COSPHI;
    end

    XCORSTA=XCORSTA/1000;

%  Finished.

%+----------------------------------------------------------------------
%
%  Copyright (C) 2008
%  IERS Conventions Center
%
%  ==================================
%  IERS Conventions Software License
%  ==================================
%
%  NOTICE TO USER:
%
%  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING TERMS AND CONDITIONS
%  WHICH APPLY TO ITS USE.
%
%  1. The Software is provided by the IERS Conventions Center ("the
%     Center").
%
%  2. Permission is granted to anyone to use the Software for any
%     purpose, including commercial applications, free of charge,
%     subject to the conditions and restrictions listed below.
%
%  3. You (the user) may adapt the Software and its algorithms for your
%     own purposes and you may distribute the resulting "derived work"
%     to others, provided that the derived work complies with the
%     following requirements:
%
%     a) Your work shall be clearly identified so that it cannot be
%        mistaken for IERS Conventions software and that it has been
%        neither distributed by nor endorsed by the Center.
%
%     b) Your work (including source code) must contain descriptions of
%        how the derived work is based upon and/or differs from the
%        original Software.
%
%     c) The name(s) of all modified routine(s) that you distribute
%        shall be changed.
%
%     d) The origin of the IERS Conventions components of your derived
%        work must not be misrepresented; you must not claim that you
%        wrote the original Software.
%
%     e) The source code must be included for all routine(s) that you
%        distribute.  This notice must be reproduced intact in any
%        source distribution.
%
%  4. In any published work produced by the user and which includes
%     results achieved by using the Software, you shall acknowledge
%     that the Software was used in obtaining those results.
%
%  5. The Software is provided to the user "as is" and the Center makes
%     no warranty as to its use or performance.   The Center does not
%     and cannot warrant the performance or results which the user may
%     obtain by using the Software.  The Center makes no warranties,
%     express or implied, as to non-infringement of third party rights,
%     merchantability, or fitness for any particular purpose.  In no
%     event will the Center be liable to the user for any consequential,
%     incidental, or special damages, including any lost profits or lost
%     savings, even if a Center representative has been advised of such
%     damages, or for any claim by any third party.
%
%  Correspondence concerning IERS Conventions software should be
%  addressed as follows:
%
%                     Gerard Petit
%     Internet email: gpetit[at]bipm.org
%     Postal address: IERS Conventions Center
%                     Time, frequency and gravimetry section, BIPM
%                     Pavillon de Breteuil
%                     92312 Sevres  FRANCE
%
%     or
%
%                     Brian Luzum
%     Internet email: brian.luzum[at]usno.navy.mil
%     Postal address: IERS Conventions Center
%                     Earth Orientation Department
%                     3450 Massachusetts Ave, NW
%                     Washington, DC 20392
%
%
%-----------------------------------------------------------------------
end

function XCORSTA=longBandCorr(XSTA,T)
%%DIURBANDCORR  This subroutine is a Matlab translation of the subroutine
%               STEP2LON that is provided by the IERS at
%               ftp://tai.bipm.org/iers/convupdt/chapter7/
%               The algorithm is not fully documented by the IERS 2010
%               conventions. For example, most of the constants are not
%               listed anywhere in the documentation.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%
%SUBROUTINE STEP2LON (XSTA,T,XCORSTA)
%
%  - - - - - - - - - - -
%   S T E P 2 L O N
%  - - - - - - - - - - -
%
%  This routine is part of the International Earth Rotation and
%  Reference Systems Service (IERS) Conventions software collection.
%
%  This subroutine gives the in-phase and out-of-phase corrections
%  induced by mantle anelasticity in the long period band.
%
%  In general, Class 1, 2, and 3 models represent physical effects that
%  act on geodetic parameters while canonical models provide lower-level
%  representations or basic computations that are used by Class 1, 2, or
%  3 models.
%
%  Status: Class 1
%
%     Class 1 models are those recommended to be used a priori in the
%     reduction of raw space geodetic data in order to determine
%     geodetic parameter estimates.
%     Class 2 models are those that eliminate an observational
%     singularity and are purely conventional in nature.
%     Class 3 models are those that are not required as either Class
%     1 or 2.
%     Canonical models are accepted as is and cannot be classified as a
%     Class 1, 2, or 3 model.
%
%  Given:
%     XSTA          d(3)   Geocentric position of the IGS station (Note 1)
%     T             d      Centuries since J2000
%
%  Returned:
%     XCORSTA       d(3)   In phase and out of phase station corrections
%                          for diurnal band (Note 2)
%
%  Notes:
%
%  1) The IGS station is in ITRF co-rotating frame.  All coordinates are
%     expressed in meters.
%
%  2) All coordinates are expressed in meters.
%
%  Test case:
%     given input: XSTA(1) = 4075578.385D0 meters
%                  XSTA(2) =  931852.890D0 meters
%                  XSTA(3) = 4801570.154D0 meters
%                  T       = 0.1059411362080767D0 Julian centuries
%
%     expected output:  XCORSTA(1) = -0.9780962849562107762D-04 meters
%                       XCORSTA(2) = -0.2236349699932734273D-04 meters
%                       XCORSTA(3) =  0.3561945821351565926D-03 meters
%
%  References:
%
%     Mathews, P. M., Dehant, V., and Gipson, J. M., 1997, ''Tidal station
%     displacements," J. Geophys. Res., 102(B9), pp. 20,469-20,477
%
%     Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
%     IERS Technical Note No. 36, BKG (2010)
%
%  Revisions:
%  1996 March    23 V. Dehant      Original code
%  2009 August   07 B.E. Stetzler  Initial standardization of code
%                                  and found unnecessary variables tau
%                                  and fhr
%  2009 August   07 B.E. Stetzler  Provided a test case
%  2009 August   07 B.E. Stetzler  Capitalized all variables for
%                                  Fortran 77 compatibility
%  2010 October  20 B.E. Stetzler  Input T corrected to be number of
%                                  centuries since J2000
%-----------------------------------------------------------------------

    D2PI = 6.283185307179586476925287;

    DATDI=[0, 0, 0, 1, 0,   0.47, 0.23, 0.16, 0.07;
         0, 2, 0, 0, 0,  -0.20,-0.12,-0.11,-0.05;
         1, 0,-1, 0, 0,  -0.11,-0.08,-0.09,-0.04;
         2, 0, 0, 0, 0,  -0.13,-0.11,-0.15,-0.07;
         2, 0, 0, 1, 0,  -0.05,-0.05,-0.06,-0.03];

    DEG2RAD = D2PI/360;

%  Compute the phase angles in degrees.
    S = 218.31664563+(481267.88194+(-0.0014663889+(0.00000185139)*T)*T)*T;

    PR = (1.396971278+(0.000308889+(0.000000021+(0.000000007)*T)*T)*T)*T;

    S = S + PR;

    H = 280.46645+(36000.7697489+(0.00030322222+(0.000000020+(-0.00000000654)*T)*T)*T)*T;

    P = 83.35324312+(4069.01363525+(-0.01032172222+(-0.0000124991+(0.00000005263)*T)*T)*T)*T;

    ZNS = 234.95544499+(1934.13626197+(-0.00207561111+(-0.00000213944+(0.00000001650)*T)*T)*T)*T;

    PS = 282.93734098+(1.71945766667+(0.00045688889+(-0.00000001778+(-0.00000000334)*T)*T)*T)*T;

    RSTA=norm(XSTA);
    SINPHI=XSTA(3)/RSTA;
    COSPHI=norm(XSTA(1:2))/RSTA;

    COSLA=XSTA(1)/COSPHI/RSTA;
    SINLA=XSTA(2)/COSPHI/RSTA;

% Reduce angles to between the range 0 and 360.
    S =  mod(S,360);
    %      TAU = DMOD(TAU,360D0)
    H =  mod(H,360);
    P =  mod(P,360);
    ZNS = mod(ZNS,360);
    PS = mod(PS,360);

    DR_TOT = 0;
    DN_TOT = 0;

    XCORSTA=zeros(3,1);

    for J=1:5
        THETAF=(DATDI(J,1)*S+DATDI(J,2)*H+DATDI(J,3)*P+DATDI(J,4)*ZNS+DATDI(J,5)*PS)*DEG2RAD;

        DR=DATDI(J,6)*(3D0*SINPHI^2-1)/2*cos(THETAF)+DATDI(J,8)*(3D0*SINPHI^2-1)/2*sin(THETAF);

        DN=DATDI(J,7)*(COSPHI*SINPHI*2)*cos(THETAF)+DATDI(J,9)*(COSPHI*SINPHI*2)*sin(THETAF);

        DE = 0;
        DR_TOT = DR_TOT+DR;
        DN_TOT = DN_TOT+DN;

        XCORSTA(1)=XCORSTA(1)+DR*COSLA*COSPHI-DE*SINLA-DN*SINPHI*COSLA;
        XCORSTA(2)=XCORSTA(2)+DR*SINLA*COSPHI+DE*COSLA-DN*SINPHI*SINLA;
        XCORSTA(3)=XCORSTA(3)+DR*SINPHI+DN*COSPHI;
    end

    XCORSTA=XCORSTA/1000;

%  Finished.

%+----------------------------------------------------------------------
%
%  Copyright (C) 2008
%  IERS Conventions Center
%
%  ==================================
%  IERS Conventions Software License
%  ==================================
%
%  NOTICE TO USER:
%
%  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING TERMS AND CONDITIONS
%  WHICH APPLY TO ITS USE.
%
%  1. The Software is provided by the IERS Conventions Center ("the
%     Center").
%
%  2. Permission is granted to anyone to use the Software for any
%     purpose, including commercial applications, free of charge,
%     subject to the conditions and restrictions listed below.
%
%  3. You (the user) may adapt the Software and its algorithms for your
%     own purposes and you may distribute the resulting "derived work"
%     to others, provided that the derived work complies with the
%     following requirements:
%
%     a) Your work shall be clearly identified so that it cannot be
%        mistaken for IERS Conventions software and that it has been
%        neither distributed by nor endorsed by the Center.
%
%     b) Your work (including source code) must contain descriptions of
%        how the derived work is based upon and/or differs from the
%        original Software.
%
%     c) The name(s) of all modified routine(s) that you distribute
%        shall be changed.
%
%     d) The origin of the IERS Conventions components of your derived
%        work must not be misrepresented; you must not claim that you
%        wrote the original Software.
%
%     e) The source code must be included for all routine(s) that you
%        distribute.  This notice must be reproduced intact in any
%        source distribution.
%
%  4. In any published work produced by the user and which includes
%     results achieved by using the Software, you shall acknowledge
%     that the Software was used in obtaining those results.
%
%  5. The Software is provided to the user "as is" and the Center makes
%     no warranty as to its use or performance.   The Center does not
%     and cannot warrant the performance or results which the user may
%     obtain by using the Software.  The Center makes no warranties,
%     express or implied, as to non-infringement of third party rights,
%     merchantability, or fitness for any particular purpose.  In no
%     event will the Center be liable to the user for any consequential,
%     incidental, or special damages, including any lost profits or lost
%     savings, even if a Center representative has been advised of such
%     damages, or for any claim by any third party.
%
%  Correspondence concerning IERS Conventions software should be
%  addressed as follows:
%
%                     Gerard Petit
%     Internet email: gpetit[at]bipm.org
%     Postal address: IERS Conventions Center
%                     Time, frequency and gravimetry section, BIPM
%                     Pavillon de Breteuil
%                     92312 Sevres  FRANCE
%
%     or
%
%                     Brian Luzum
%     Internet email: brian.luzum[at]usno.navy.mil
%     Postal address: IERS Conventions Center
%                     Earth Orientation Department
%                     3450 Massachusetts Ave, NW
%                     Washington, DC 20392
%
%
%-----------------------------------------------------------------------
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 D = cluster_discrim(cl,beh,varargin)
% D = cluster_discrim(cl,beh,[covs of no interest])
%
%x = cat(2,pain_matrix_early(:).timeseries);
%x = cat(2,pain_matrix_peak(:).timeseries);

diary Cluster_Discrim_Output.txt

disp('cluster_discrim.m: Predictions of behavior with contrast scores');
fprintf(1,'\n-----------------------------------------------------------------\n');

disp(['Input clusters: ' inputname(1)])
disp('Loading data from cl.timeseries');

% Get values from clusters
% ---------------------------------------------
x = cat(2,cl(:).timeseries);
if any(isnan(sum(x))),
    whnan = find(isnan(sum(x)));
    cl(whnan) = []; x = cat(2,cl(:).timeseries);
    disp(['Found NaNs! Removing clusters ' num2str(whnan)]);
end

% Remove covariates of no interest
% ---------------------------------------------
covs = [];
if length(varargin) > 0
    covs = varargin{1};
    D.covs = covs;
    disp('Found covariate(s) of no interest; Removing them from behavior and brain data.')
        % partialcorr gets adjusted x, y, and correls; here just use it to
        % get adjusted x
        beh  = partialcor([beh covs],ones(size(beh,1),1),1);

    for i = 1:size(x,2)
        x(:,i) = partialcor([x(:,i) covs],ones(size(x,1),1),1);
    end
end

D.data = x;

figure;imagesc(x); colorbar; title('Data for discriminant analysis'); drawnow

D.xyzmm = cat(1,cl(:).mm_center);

y = beh;
D.beh = beh;

y2 = mediansplit(y);
D.beh_medsplit = y2;

if isfield(cl,'imnames'),D.imnames = cl(1).imnames;,end

% Print univariate correlations with regions and location (t-test)
% ---------------------------------------------
disp('Correlations with behavior, after removing covariates of no interest.')
fprintf(1,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t\n','Cluster','Name','x','y','z','V','r', 'p', 'IRLSr', 'IRLSp','Zavg','Inc/Dec');
warning off
fprintf(1,'\n');
for i = 1:size(x,2)
    [dummy,dummy,r,p,rrob,prob]  = partialcor([beh covs],x(:,i),1);

     % print output
     fprintf(1,'Cluster %3.0f\t',i);
     if isfield(cl,'shorttitle'), fprintf(1,'%s\t',cl(i).shorttitle),end
     if isfield(cl,'mm_center'),fprintf(1,'%3.0f\t%3.0f\t%3.0f\t',cl(i).mm_center);,end
     if isfield(cl,'numVox'),fprintf(1,'%3.0f\t',cl(i).numVox);,end
     fprintf(1,'%3.2f\t%3.2f\t%3.2f\t%3.2f\t',[r p rrob prob])

     % t-test
     [hh,p,ci,stats] = ttest(cl(i).timeseries);
     strs = {'Decrease' '--' 'Increase'};
     if ~hh, strs = strs{2};,elseif stats.tstat >0,strs=strs{3};,else,strs = strs{1};end
     fprintf(1,'%3.2f\t%s\t',stats.tstat,strs)

     fprintf(1,'\n');
end
warning on

% --------------------------------------
% get eigenvalues, choose # dims
% --------------------------------------
disp('Computing significant principal components')
[v,score,e] = princomp(x);
[pc,D.pca_npm] = pca_npm(x,1000);
whsig = find(cumsum(D.pca_npm.sig) >= 1:length(D.pca_npm.sig));
ndim = length(whsig) + 1;

try, saveas(gcf,'Cluster_eigenplot','tif');,end

% --------------------------------------
% predict behavior with component scores
% --------------------------------------

fprintf(1,'\n-----------------------------------------------------------------\n');
fprintf(1,'Stepwise regression: Predictions of behavior with component scores');
fprintf(1,'\n-----------------------------------------------------------------\n');

D.PCR = stepwise_tor(score(:,1:ndim),y);

sig = find(D.PCR.inmodel);
D.PCR.sscore = score(:,sig);    % scores for significant components

% plot all significant (now done in cluster_discrim_montage)
%figure('Color','w');
%for i = 1:size(D.PCR.sscore,2)
%    subplot(1,size(D.PCR.sscore,2),i);
%    discrim_plot(D.PCR.sscore(:,i),y,0);
%end

D.PCR.seigv = v(:,sig);         % eigenvectors for sig components
D.PCR.seig = e(sig);            % significant eigenvalues

tmpc = corrcoef([D.PCR.sscore x]);
tmpc(1:size(D.PCR.sscore,2),:) = []; tmpc = tmpc(:,1:size(D.PCR.sscore,2));
D.PCR.compcorr = tmpc;
D.PCR.compcorr_descrip = 'Regions x sig. comps, correlation with component values';
D.PCR.whichcomps = sign(D.PCR.compcorr) .* real(abs(D.PCR.compcorr) > .4);

tmpc(abs(tmpc) < .4) = 0;
D.PCR.threshcomps = tmpc;

% --------------------------------------
% manova to get categorical canonical discrim functions
% this gives us overall p-value and discrim functions
% --------------------------------------
[manovasig,manovap,D.manova] = manova1(score(:,1:ndim),y2);
D.manova.p = manovap;
fprintf(1,'\n-----------------------------------------------------------------\n');
fprintf(1,'Manova to discriminate high from low behavior');
fprintf(1,'\n-----------------------------------------------------------------\n');
fprintf(1,'Wilks lambda (%3.0f,%3.0f) = %3.2f\tChi-sq (%3.0f) = %3.2f\t, p = %3.4f\t\n',D.manova.dfB,D.manova.dfW,D.manova.lambda,D.manova.chisqdf,D.manova.chisq,D.manova.p);
%discrim_plot(D.manova.canon(:,1),y);
fprintf(1,'Eigenvalues\t')
fprintf(1,'%3.2f\t',D.manova.eigenval')
fprintf(1,'\n')
fprintf(1,'Eigenvector 1\t')
fprintf(1,'%3.2f\t',D.manova.eigenvec(:,1)')
fprintf(1,'\n')

% --------------------------------------
% visualize clusters on brain
% --------------------------------------

cluster_discrim_montage(cl,D.PCR.threshcomps,D.PCR.sscore,D.beh);

% --------------------------------------
% manova on individual regions??
% --------------------------------------

diary off

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [H, inls] = at_prosac_p3p(x1, x2, udists, TN, TLmax, tol, doLO, conf)

num_points = 3;
if nargin < 7
  conf = .95;     % in case, conf is not specified.
end;
len = size(x1,2);    % the maxinum number of points, i.e. "N" in PROSAC paper.
max_i = 3; % the number of inliers.
k_star = TN;  % the maximum number for terminating iteration. "n*" in PROSAC paper.

t = 0;
n_sub = num_points; % the number of subset points corresponds to "n" in PROSAC paper.
Tn = prosac_avsamplen(num_points,len,TN); % compute the average sampling number Tn'.
n_star = len; % the number of sampling=iteration corresponds to "n*" in PROSAC paper.

if ~exist('prosac_ImnP3P.mat','file') %|| len > 20000
  Imn = zeros(1,len);
  %%%%adhoc = round(len*(1-conf));
  adhoc = 0;
  for ii = 1:num_points+adhoc
    Imn(ii) = len;
  end
  for ii = num_points+1+adhoc:len
    Imn(ii) = prosac_nonrandom(num_points,ii,conf,0.1,len);
    fprintf('%d\n',ii)
  end
else
  r = load('prosac_ImnP3P.mat');
  Imn = r.Imn; clear r;
end;

H = [];
inls = [];

[~, sindx] = sort(udists);
tmax = TLmax;

x1h = x1;
x2h = x2;

x1h(end+1,:) = 1 ;
x2h(end+1,:) = 1 ;

while (t < k_star && t < tmax)
  t = t+1; %msg(1,'%d\n',t);
  if (t == Tn(n_sub) && n_sub < n_star)
    n_sub = n_sub + 1;
  end

  p_sample = sindx(prosac_rsample(num_points,Tn,t,n_sub));

  x1in = x1h(:,p_sample);
  x2in = x2h(:,p_sample);

  y = [x1in; x2in];

  pf = P3PSolver(y);

  e = PerspRepErr(P,X,K)

%   y = xf(:,i); % get a sample
%   pf = fitF(y); % fit a model
  r  = resF(pf,xr,cnstr); % eval residuals & constraints
  il = abs(r)<=maxres; % inliers
  if size(il,1)>1 % there are more alternative models returned by fitF
    in = sum(il,2); % inlier #
    [in,ix] = max(in);
    il = il(ix,:); % select the best inliers
    r = r(ix,:); % select the residuals
    pf = pf{ix}; % select the best model
  else
    in = sum(il); % inlier #
    if ~isempty(pf)
      if iscell(pf)
        pf = pf{1};
      end
    end
  end

  S1 = centering(x1in) ;
  S2 = centering(x2in) ;
  x1c = S1 * x1in ;
  x2c = S2 * x2in ;
  M = [x1c, zeros(size(x1c)) ;
    zeros(size(x1c)), x1c ;
    bsxfun(@times, x1c,  -x2c(1,:)), bsxfun(@times, x1c,  -x2c(2,:))] ;
  if any(isnan(M(:))), continue; end;
  if any(isinf(M(:))), continue; end;
  [H21,D] = svd(M,'econ') ;
  H21 = reshape(H21(:,end),3,3)' ;
  H21 = S2 \ H21 * S1 ;
  H21 = H21 ./ H21(end) ;

%   H21 = u2H([x2h(:,p_sample); x1h(:,p_sample)]);
%   if isempty(H21), continue; end;
%   H21 = H21 ./ H21(end);

  x1p = H21 * x1h;
  x1p = [x1p(1,:) ./ x1p(3,:) ; x1p(2,:) ./ x1p(3,:)] ;

%   tol = th * sqrt(det(H21(1:2,1:2)));
  dist2 = sum((x2 - x1p).^2,1) ;
  v = dist2 < tol^2;
  no_i = sum(v);

  if( max_i < no_i)
    if doLO && no_i > 4
      bH = u2H([x2h(:,v); x1h(:,v)]);
      bH = bH ./ bH(end);

      x1p = bH * x1h;
      x1p = [x1p(1,:) ./ x1p(3,:) ; x1p(2,:) ./ x1p(3,:)] ;
%       tol = th * sqrt(det(bH(1:2,1:2)));
      bdist2 = sum((x2 - x1p).^2,1) ;
      bv = bdist2 < tol^2;
      no_ib = sum(bv);
      if no_ib > no_i
        no_i = no_ib;
        H21 = bH;
        v = bv;
%         fprintf(1,'LO succeeded\n');
      end;
    end;
    max_i = no_i;
    H = H21;
    inls = find(v);

    %% STEP 6 maximality of iteration
    for cn = num_points:len
      Icn = sum(v(1:cn));
      kcn = nsamples(Icn, cn, num_points, conf);%msg(1,'kcn=%f\n',kcn);
      %                inlfrc = [inlfrc, Icn/cn];
      if cn < n_sub
        kcn = kcn - (t - Tn(cn));
      end;
      if Icn >= Imn(cn) && kcn < k_star
        n_star = cn;
        k_star = kcn;
      end;
    end;

%     k_star = nsamples(max_i, len, num_points, conf);
%     fprintf(1,'PROSAC(H) iter=%d, #inl=%d, k_star=%2.2f\n',t,max_i,k_star);
  end;
end;
% fprintf(1,'PROSAC(H) iter=%d, #inl=%d, k_star=%2.2f\n',t,max_i,k_star);

function urs_indx = prosac_rsample(m,Tn,t,n)
% prosac sampling
if Tn(n) > t
  tmp = randperm(n-1);
  urs_indx = [tmp(1:m-1),n];
else
  tmp = randperm(n-1);
  urs_indx = tmp(1:m);
  %msg(1,'why?\n');
end

function [SampleCnt, q] = nsamples(ni, ptNum, pf, conf)
%SampleCnt calculates number of samples needed to be done

q  = prod (((ni-pf+1) : ni) ./ ((ptNum-pf+1) : ptNum));
if q < eps
  SampleCnt = Inf;
else
  %   SampleCnt  = log(1 - conf) / log(1 - q);
  if q > conf
    SampleCnt = 1;
  else
    SampleCnt  = log(1 - conf) / log(1 - q);
  end;
end;

function Imin = prosac_nonrandom(m,n,conf,beta,len)
% compute the Imin whch express the non-randomness.

phi = 1-conf;
ii = (m:n);
Rn(ii) = binopdf(ii-m,n-m,beta); % see Eq. (7) in PROSAC paper.
for j = m:n
  if sum(Rn(j:n)) < phi
    Imin = j; break;   % see Eq. (8) in PROSAC paper.
  else
      Imin = len;
  end
end

author: HajimeTaira


text:
classdef BadFixture < TestCase

    methods
        function self = BadFixture(name)
            self = self@TestCase(name);
        end

        function setUp(self)
            throw(MException('setUpError:BadFixture', ...
                'BadFixture setUp method always throws exception'));
        end

        function testMethod(self)
        end
    end
end

author: ricedsp


text:
function PlotMesh(coordinates,nodes)
%--------------------------------------------------------------------------
% Code written by : Siva Srinivas Kolukula                                |
%                   Senior Research Fellow                                |
%                   Structural Mechanics Laboratory                       |
%                   Indira Gandhi Center for Atomic Research              |
%                   India                                                 |
% E-mail : allwayzitzme@gmail.com                                         |
%          http://sites.google.com/site/kolukulasivasrinivas/             |
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Purpose:
%         To plot the Finite Element Method Mesh
% Synopsis :
%           PlotMesh(coordinates,nodes)
% Variable Description:
%           coordinates - The nodal coordinates of the mesh
%           -----> coordinates = [node X Y]
%           nodes - The nodal connectivity of the elements
%           -----> nodes = [node1 node2......]
%--------------------------------------------------------------------------

nel = length(nodes) ;                  % number of elements
nnode = length(coordinates) ;          % total number of nodes in system
nnel = size(nodes,2);                % number of nodes per element
%
% Initialization of the required matrices
X = zeros(nnel,nel) ;
Y = zeros(nnel,nel) ;

for iel=1:nel
     for i=1:nnel
     nd(i)=nodes(iel,i);         % extract connected node for (iel)-th element
     X(i,iel)=coordinates(nd(i),1);    % extract x value of the node
     Y(i,iel)=coordinates(nd(i),2);    % extract y value of the node
     end
end

% Plotting the FEM mesh, diaplay Node numbers and Element numbers
     f1 = figure ;
     set(f1,'name','Mesh','numbertitle','off','Color','w') ;
     fill(X,Y,'w')

     title('Finite Element Mesh') ;
     axis off ;

% To disply the node numbers
%      k = nodes(:,1:end);
%      nd = k' ;
%     for i = 1:nel
%         text(X(:,i),Y(:,i),int2str(nd(:,i)),'fontsize',8,'color','k');
%         text(sum(X(:,i))/4,sum(Y(:,i))/4,int2str(i),'fontsize',10,'color','r') ;
%     end

author: Sable

repo: mcbench-benchmarks


text:
function [dictionary] = build_dictionary(filelist, feature, c)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

if(~exist('c', 'var'))
  c = conf();
end

p = c.feature_config.(feature);
if(~isfield(p, 'dictionary_file'))
  dictionary = [];
  return;
end
p.dictionary_file = sprintf(p.dictionary_file, c.cache, p.dictionary_size);
found_dictionary = 0;
check_building = 0;

while(found_dictionary == 0)
  if(~exist(p.dictionary_file, 'file'))
    hostname = getComputerName();
		make_dir(p.dictionary_file);
    save(p.dictionary_file, 'hostname');

    %check for multiple datasets
    if(iscell(filelist{1}))
      trainlists = filelist;
      num_datasets = length(trainlists);
      images_per_dataset = ceil(p.num_images/num_datasets);
      filelists = cellfun(@(x) x(randperm(length(x), min(length(x), images_per_dataset))), trainlists, 'UniformOutput', false);
      filelist = {};
      for i=1:length(filelists)
				if(size(filelists{i},1)~=1), filelists{i}=filelists{i}'; end
        filelist = [filelist filelists{i}];
      end
    end

		vprintf(c.verbosity, 0, 'Learning dictionary for feature: %s, size %d\n', feature, p.dictionary_size);
    perm = randperm(length(filelist));
    descriptors = cell(min(length(filelist), p.num_images), 1);
    num_images = min(length(filelist), p.num_images);
    parfor i=1:num_images
      vprintf(c.verbosity, 1, 'Dictionary learning (%s): %d of %d\n', feature, i, num_images);
      img = imgread(filelist{perm(i)}, p);
      feat = extract_feature(feature, img, c);
      r = randperm(size(feat, 1));
      descriptors{i} = feat(r(1:min(length(r), p.descPerImage)), :);
    end
    descriptors = cell2mat(descriptors);
    ndata = size(descriptors, 1);
    if(ndata>p.num_desc)
      idx = randperm(ndata);
      descriptors = descriptors(idx(1:p.num_desc), :);
    end
    vprintf(c.verbosity, 0, 'Running k-means, dictionary size %d...', p.dictionary_size);
    dictionary = kmeansFast(descriptors, p.dictionary_size);
    vprintf(c.verbosity, 0, 'done!\n');
    vprintf(c.verbosity, 0, 'Saving dictionary: %s\n', p.dictionary_file);
    save(p.dictionary_file, 'dictionary');
    found_dictionary = 1;
  else
    load(p.dictionary_file);
    if(~exist('dictionary', 'var'))
      if(check_building == 0)
        vprintf(c.verbosity, 0, 'Dictionary building in progress on %s..', hostname);
        check_building = 1;
      end
      vprintf(c.verbosity, 1, '.');
      pause(5);
    else
      found_dictionary = 1;
      if(check_building == 1)
        vprintf(c.verbosity, 0, '\n');
      end
    end
  end
end

author: adikhosla

repo: feature-extraction


text:
function [success,message] = xlsappend(file,data,sheet)

% XLSAPPEND Stores numeric array or cell array to the end of specified Excel sheet.
%
%   REQUIRES ONLY ONE CALL TO THE EXCEL ACTXSERVER, so the overhead is less
%   than for successive xlsread/xlswrite calls.
%
%   [SUCCESS,MESSAGE]=XLSAPPEND(FILE,ARRAY,SHEET) writes ARRAY to the Excel
%   workbook, FILE, into the area beginning at COLUMN A and FIRST UNUSED
%   ROW, in the worksheet specified in SHEET.
%   FILE and ARRAY must be specified. If either FILE or ARRAY is empty, an
%   error is thrown and XLSAPPEND terminates. The first worksheet of the
%   workbook is the default. If SHEET does not exist, a new sheet is added at
%   the end of the worksheet collection. If SHEET is an index larger than the
%   number of worksheets, new sheets are appended until the number of worksheets
%   in the workbook equals SHEET. The success of the operation is
%   returned in SUCCESS and any accompanying message, in MESSAGE. On error,
%   MESSAGE shall be a struct, containing the error message and message ID.
%   See NOTE 1.
%
%   [SUCCESS,MESSAGE]=XLSAPPEND(FILE,ARRAY) writes ARRAY to the Excel
%   workbook, FILE, in the first worksheet.
%
%   [SUCCESS,MESSAGE]=XLSAPPEND(FILE,ARRAY) writes ARRAY to the Excel
%   workbook, FILE, starting at cell A[firstUnusedRow] of the first
%   worksheet. The return values are as for the above example.
%
%   XLSAPPEND ARRAY FILE, is the command line version of the above example.
%
%   INPUT PARAMETERS:
%       file:   string defining the workbook file to write to.
%               Default directory is pwd; default extension 'xls'.
%       array:  m x n numeric array or cell array.
%       sheet:  string defining worksheet name;
%               double, defining worksheet index.
%
%   RETURN PARAMETERS:
%       SUCCESS: logical scalar.
%       MESSAGE: struct containing message field and message_id field.
%
%   EXAMPLES:
%
%   SUCCESS = XLSAPPEND('c:\matlab\work\myworkbook.xls',A) will write A to
%   the sheet 1 of workbook file, myworkbook.xls. On success, SUCCESS will
%   contain true, while on failure, SUCCESS will contain false.
%
%   NOTE 1: The above functionality depends upon Excel as a COM server. In
%   absence of Excel, ARRAY shall be written as a text file in CSV format. In
%   this mode, the SHEET argument shall be ignored.
%
% See also XLSREAD, XLSWRITE.
%

% Written by Brett Shoelson, PhD, 8/30/2010.
% Copyright 1984-2010 The MathWorks, Inc.
%==============================================================================

% Set default values.
lastwarn('');
Sheet1 = 1;

if nargin < 3
    sheet = Sheet1;
end

if nargout > 0
    success = true;
    message = struct('message',{''},'identifier',{''});
end

% Handle input.
try
    % handle requested Excel workbook filename.
    if ~isempty(file)
        if ~ischar(file)
            error('MATLAB:xlswrite:InputClass','Filename must be a string.');
        end
        % check for wildcards in filename
        if any(findstr('*', file))
            error('MATLAB:xlswrite:FileName', 'Filename must not contain *.');
        end
        [Directory,file,ext]=fileparts(file);

        if isempty(ext) % add default Excel extension;
            ext = '.xls';
        end
        file = abspath(fullfile(Directory,[file ext]));
        [a1 a2] = fileattrib(file);
        if a1 && ~(a2.UserWrite == 1)
            error('MATLAB:xlswrite:FileReadOnly', 'File cannot be read-only.');
        end
    else % get workbook filename.
        error('MATLAB:xlswrite:EmptyFileName','Filename is empty.');
    end

    % Check for empty input data
    if isempty(data)
        error('MATLAB:xlswrite:EmptyInput','Input array is empty.');
    end

    % Check for N-D array input data
    if ndims(data)>2
        error('MATLAB:xlswrite:InputDimension',...
            'Dimension of input array cannot be higher than two.');
    end

    % Check class of input data
    if ~(iscell(data) || isnumeric(data) || ischar(data)) && ~islogical(data)
        error('MATLAB:xlswrite:InputClass',...
            'Input data must be a numeric, cell, or logical array.');
    end

    % convert input to cell array of data.
     if iscell(data)
        A=data;
     else
         A=num2cell(data);
     end

    if nargin > 2
        % Verify class of sheet parameter.
        if ~(ischar(sheet) || (isnumeric(sheet) && sheet > 0))
            error('MATLAB:xlswrite:InputClass',...
                'Sheet argument must be a string or a whole number greater than 0.');
        end
        if isempty(sheet)
            sheet = Sheet1;
        end
        % Parse sheet
        if ischar(sheet) && ~isempty(strfind(sheet,':'))
            sheet = Sheet1;% Use default sheet.
        end
    end

catch exception
    if ~isempty(nargchk(2,4,nargin))
        error('MATLAB:xlswrite:InputArguments',nargchk(2,4,nargin));
    else
        success = false;
        message = exceptionHandler(nargout, exception);
    end
    return;
end
%------------------------------------------------------------------------------
% Attempt to start Excel as ActiveX server.
try
    Excel = actxserver('Excel.Application');
    % open workbook
    Excel.DisplayAlerts = 0;

    %Workaround for G313142.  For certain files, unless a workbook is
    %opened prior to opening the file, various COM calls return an error:
    %0x800a9c64.  The line below works around this flaw.  Since we have
    %seen only one example of such a file, we have decided not to incur the
    %time penalty involved here.
    %     aTemp = Excel.workbooks.Add(); aTemp.Close();

    ExcelWorkbook = Excel.workbooks.Open(file);
    %ExcelWorkbook.ReadOnly
    format = ExcelWorkbook.FileFormat;
    if  strcmpi(format, 'xlCurrentPlatformText') == 1
        error('MATLAB:xlsread:FileFormat', 'File %s not in Microsoft Excel Format.', file);
    end

    %Sheets = Excel.ActiveWorkBook.Sheets;
    activate_sheet(Excel,sheet);
    readinfo = get(Excel.Activesheet,'UsedRange');
    %get(Excel.Activesheet,'Name')
    if numel(readinfo.value)== 1 && isnan(readinfo.value)
        m2 = 0;
    else
        [m2,n2] = size(readinfo.value);
    end

catch exception %#ok<NASGU>
    warning('MATLAB:xlswrite:NoCOMServer',...
        ['Could not start Excel server for export.\n' ...
        'XLSWRITE will attempt to write file in CSV format.']);
    if nargout > 0
        [message.message,message.identifier] = lastwarn;
    end
    % write data as CSV file, that is, comma delimited.
    file = regexprep(file,'(\.xls[^.]*+)$','.csv');
    try
        dlmwrite(file,data,','); % write data.
    catch exception
        exceptionNew = MException('MATLAB:xlswrite:dlmwrite', 'An error occurred on data export in CSV format.');
        exceptionNew = exceptionNew.addCause(exception);
        if nargout == 0
            % Throw error.
            throw(exceptionNew);
        else
            success = false;
            message.message = exceptionNew.getReport;
            message.identifier = exceptionNew.identifier;
        end
    end
    return;
end
%------------------------------------------------------------------------------
try
    % Construct range string
    % Range was partly specified or not at all. Calculate range.
    % Select range of occupied cells in active sheet.
    % Activate indicated worksheet.
    message = activate_sheet(Excel,sheet);
    [m,n] = size(A);
    range = calcrange('',m,n,m2);
catch exception
    success = false;
    message = exceptionHandler(nargout, exception);
    return;
end

%------------------------------------------------------------------------------
try
    bCreated = ~exist(file,'file');
    ExecuteWrite;
catch exception
    if (bCreated && exist(file, 'file') == 2)
        delete(file);
    end
    success = false;
    message = exceptionHandler(nargout, exception);
end
    function ExecuteWrite
            cleanUp = onCleanup(@()cleaner(Excel, file));
            if bCreated
                % Create new workbook.
                %This is in place because in the presence of a Google Desktop
                %Search installation, calling Add, and then SaveAs after adding data,
                %to create a new Excel file, will leave an Excel process hanging.
                %This workaround prevents it from happening, by creating a blank file,
                %and saving it.  It can then be opened with Open.

                %ExcelWorkbook = Excel.workbooks.Add;
                switch ext
                    case '.xls' %xlExcel8 or xlWorkbookNormal
                       xlFormat = -4143;
                    case '.xlsb' %xlExcel12
                       xlFormat = 50;
                    case '.xlsx' %xlOpenXMLWorkbook
                       xlFormat = 51;
                    case '.xlsm' %xlOpenXMLWorkbookMacroEnabled
                       xlFormat = 52;
                    otherwise
                       xlFormat = -4143;
                end
                ExcelWorkbook.SaveAs(file, xlFormat);
                ExcelWorkbook.Close(false);
            end

            %Open file
            %ExcelWorkbook = Excel.workbooks.Open(file);
            if ExcelWorkbook.ReadOnly ~= 0
                %This means the file is probably open in another process.
                error('MATLAB:xlswrite:LockedFile', 'The file %s is not writable.  It may be locked by another process.', file);
            end
            try % select region.
                % Select range in worksheet.
                Select(Range(Excel,sprintf('%s',range)));

            catch exceptionInner % Throw data range error.
                throw(MException('MATLAB:xlswrite:SelectDataRange', sprintf('Excel returned: %s.', exceptionInner.message)));
            end

            % Export data to selected region.
            set(Excel.selection,'Value',A);
            ExcelWorkbook.Save
    end
end
function cleaner(Excel, filePath)
    try
        %Turn off dialog boxes as we close the file and quit Excel.
        Excel.DisplayAlerts = 0;
        %Explicitly close the file just in case.  The Excel API expects
        %just the filename and not the path.  This is safe because Excel
        %also does not allow opening two files with the same name in
        %different folders at the same time.
        [~, n, e] = fileparts(filePath);
        fileName = [n e];
        Excel.Workbooks.Item(fileName).Close(false);
    catch exception %#ok<NASGU>
        %If something fails in closing, there is nothing to do but attempt
        %to quit.
    end
    Excel.Quit;
end
%--------------------------------------------------------------------------
function message = activate_sheet(Excel,Sheet)
% Activate specified worksheet in workbook.

% Initialise worksheet object
WorkSheets = Excel.sheets;
message = struct('message',{''},'identifier',{''});

% Get name of specified worksheet from workbook
try
    TargetSheet = get(WorkSheets,'item',Sheet);
catch exception  %#ok<NASGU>
    % Worksheet does not exist. Add worksheet.
    TargetSheet = addsheet(WorkSheets,Sheet);
    warning('MATLAB:xlswrite:AddSheet','Added specified worksheet.');
    if nargout > 0
        [message.message,message.identifier] = lastwarn;
    end
end

% activate worksheet
Activate(TargetSheet);
end
%------------------------------------------------------------------------------
function newsheet = addsheet(WorkSheets,Sheet)
% Add new worksheet, Sheet into worsheet collection, WorkSheets.

if isnumeric(Sheet)
    % iteratively add worksheet by index until number of sheets == Sheet.
    while WorkSheets.Count < Sheet
        % find last sheet in worksheet collection
        lastsheet = WorkSheets.Item(WorkSheets.Count);
        newsheet = WorkSheets.Add([],lastsheet);
    end
else
    % add worksheet by name.
    % find last sheet in worksheet collection
    lastsheet = WorkSheets.Item(WorkSheets.Count);
    newsheet = WorkSheets.Add([],lastsheet);
end
% If Sheet is a string, rename new sheet to this string.
if ischar(Sheet)
    set(newsheet,'Name',Sheet);
end
end
%------------------------------------------------------------------------------
function [absolutepath]=abspath(partialpath)

% parse partial path into path parts
[pathname filename ext] = fileparts(partialpath);
% no path qualification is present in partial path; assume parent is pwd, except
% when path string starts with '~' or is identical to '~'.
if isempty(pathname) && isempty(strmatch('~',partialpath))
    Directory = pwd;
elseif isempty(regexp(partialpath,'(.:|\\\\)','once')) && ...
        isempty(strmatch('/',partialpath)) && ...
        isempty(strmatch('~',partialpath));
    % path did not start with any of drive name, UNC path or '~'.
    Directory = [pwd,filesep,pathname];
else
    % path content present in partial path; assume relative to current directory,
    % or absolute.
    Directory = pathname;
end

% construct absulute filename
absolutepath = fullfile(Directory,[filename,ext]);
end
%------------------------------------------------------------------------------
function range = calcrange(range,m,n,offset)
% Calculate full target range, in Excel A1 notation, to include array of size
% m x n

range = upper(range);
cols = isletter(range);
rows = ~cols;
% Construct first row.
if ~any(rows)
    firstrow = offset+1; % Default row.
else
    firstrow = str2double(range(rows)); % from range input.
end
% Construct first column.
if ~any(cols)
    firstcol = 'A'; % Default column.
else
    firstcol = range(cols); % from range input.
end
try
    lastrow = num2str(firstrow+m-1);   % Construct last row as a string.
    firstrow = num2str(firstrow);      % Convert first row to string image.
    lastcol = dec2base27(base27dec(firstcol)+n-1); % Construct last column.

    range = [firstcol firstrow ':' lastcol lastrow]; % Final range string.
catch exception
    error('MATLAB:xlswrite:CalculateRange', 'Invalid data range: %s.', range);
end
end
%----------------------------------------------------------------------
function string = index_to_string(index, first_in_range, digits)

letters = 'A':'Z';
working_index = index - first_in_range;
outputs = cell(1,digits);
[outputs{1:digits}] = ind2sub(repmat(26,1,digits), working_index);
string = fliplr(letters([outputs{:}]));
end
%----------------------------------------------------------------------
function [digits first_in_range] = calculate_range(num_to_convert)

digits = 1;
first_in_range = 0;
current_sum = 26;
while num_to_convert > current_sum
    digits = digits + 1;
    first_in_range = current_sum;
    current_sum = first_in_range + 26.^digits;
end
end
%------------------------------------------------------------------------------
function s = dec2base27(d)

%   DEC2BASE27(D) returns the representation of D as a string in base 27,
%   expressed as 'A'..'Z', 'AA','AB'...'AZ', and so on. Note, there is no zero
%   digit, so strictly we have hybrid base26, base27 number system.  D must be a
%   negative integer bigger than 0 and smaller than 2^52.
%
%   Examples
%       dec2base(1) returns 'A'
%       dec2base(26) returns 'Z'
%       dec2base(27) returns 'AA'
%-----------------------------------------------------------------------------

d = d(:);
if d ~= floor(d) || any(d(:) < 0) || any(d(:) > 1/eps)
    error('MATLAB:xlswrite:Dec2BaseInput',...
        'D must be an integer, 0 <= D <= 2^52.');
end
[num_digits begin] = calculate_range(d);
s = index_to_string(d, begin, num_digits);
end
%------------------------------------------------------------------------------
function d = base27dec(s)
%   BASE27DEC(S) returns the decimal of string S which represents a number in
%   base 27, expressed as 'A'..'Z', 'AA','AB'...'AZ', and so on. Note, there is
%   no zero so strictly we have hybrid base26, base27 number system.
%
%   Examples
%       base27dec('A') returns 1
%       base27dec('Z') returns 26
%       base27dec('IV') returns 256
%-----------------------------------------------------------------------------

if length(s) == 1
   d = s(1) -'A' + 1;
else
    cumulative = 0;
    for i = 1:numel(s)-1
        cumulative = cumulative + 26.^i;
    end
    indexes_fliped = 1 + s - 'A';
    indexes = fliplr(indexes_fliped);
    indexes_in_cells = mat2cell(indexes, 1, ones(1,numel(indexes)));
    d = cumulative + sub2ind(repmat(26, 1,numel(s)), indexes_in_cells{:});
end
end
%-------------------------------------------------------------------------------

function messageStruct = exceptionHandler(nArgs, exception)
    if nArgs == 0
        throwAsCaller(exception);
    else
        messageStruct.message = exception.message;
        messageStruct.identifier = exception.identifier;
    end
end


text:
function transpose(f) %#ok<*INUSD>
%TRANSPOSE   CLASSICFUN objects are not transposable.

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

error('CHEBFUN:CLASSICFUN:transpose:notpossible', ...
    'CLASSICFUN objects are not transposable.')

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function y = sinc_fn(x)

i=find(x==0);
x(i)= 1;
y = sin(pi*x)./(pi*x);
y(i) = 1;

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function tname = treepart(item, dflag)

% function tname = treepart(item, dflag)
% tree part to search - for cfg_repeat/cfg_choice this is val for filled
% cfg_items and values for default cfg_items.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: treepart.m 4867 2012-08-30 13:04:51Z volkmar $

rev = '$Rev: 4867 $'; %#ok

if dflag
    tname = 'values';
else
    tname = 'val';
end;

author: spm


text:
function igood = get_good_channels(ops, chanMap)
% of the channels indicated by the user as good (chanMap)
% further subset those that have a mean firing rate above a certain value
% (default is ops.minfr_goodchannels = 0.1Hz)
% needs the same filtering parameters in ops as usual
% also needs to know where to start processing batches (twind)
% and how many channels there are in total (NchanTOT)

Nbatch = ops.Nbatch;
twind = ops.twind;
NchanTOT = ops.NchanTOT;
NT = ops.NT;
Nchan = numel(chanMap);

% load data into patches, filter, compute covariance
if isfield(ops,'fslow')&&ops.fslow<ops.fs/2
    [b1, a1] = butter(3, [ops.fshigh/ops.fs,ops.fslow/ops.fs]*2, 'bandpass');
else
    [b1, a1] = butter(3, ops.fshigh/ops.fs*2, 'high');
end

fid = fopen(ops.fbinary, 'r');
% irange = [NT/8:(NT-NT/8)];

ibatch = 1;
ich = gpuArray.zeros(5e4,1, 'int16');
k = 0;
ttime = 0;

% from a subset of batches, count threshold crossings
while ibatch<=Nbatch
    offset = twind + 2*NchanTOT*NT* (ibatch-1);
    fseek(fid, offset, 'bof');
    buff = fread(fid, [NchanTOT NT], '*int16');

    if isempty(buff)
        break;
    end

    datr    = gpufilter(buff, ops, chanMap); % apply filters and median subtraction

    % very basic threshold crossings calculation
    datr = datr./std(datr,1,1); % standardize each channel ( but don't whiten)

    mdat = my_min(datr, 30, 1); % get local minima as min value in +/- 30-sample range
    ind = find(datr<mdat+1e-3 & datr<ops.spkTh); % take local minima that cross the negative threshold
    [xi, xj] = ind2sub(size(datr), ind); % back to two-dimensional indexing
    xj(xi<ops.nt0 | xi>NT-ops.nt0) = []; % filtering may create transients at beginning or end. Remove those.
    if k+numel(xj)>numel(ich)
        ich(2*numel(ich)) = 0; % if necessary, extend the variable which holds the spikes
    end
    ich(k + [1:numel(xj)]) = xj; % collect the channel identities for the detected spikes

    k = k + numel(xj);

    ibatch = ibatch + ceil(Nbatch/100); % skip every 100 batches
    ttime = ttime + size(datr,1)/ops.fs; % keep track of total time where we took spikes from
end
fclose(fid);

ich = ich(1:k);

nc = histcounts(ich, .5 + [0:Nchan]); % count how many spikes each channel got
nc = nc/ttime; % divide by total time to get firing rate

% igood = nc>.1;
igood = nc>=getOr(ops, 'minfr_goodchannels', .1); % keep only those channels above the preset mean firing rate

fprintf('found %d threshold crossings in %2.2f seconds of data \n', k, ttime)
fprintf('found %d bad channels \n', sum(~igood))

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function X=get_X_from_xP_lin(xn,Pcomp)
%get_X_from_xP_lin  Estimation of 3D point from image matches and camera matrices, linear.
%  return a 3D point X (column 4-vector, homogeneous, last element 1)
%   from its projections in K images xn (2-by-K matrix) and camera matrices
% P (K-cell of 3-by-4 matrices) by minimizing algebraic distance.
%   xn 3 x K matrix, homogeneous normalized image coordinates,
% last rows of 1, the other values usu. falls between -0.5 to 0.5
K=size(xn,2);
A = zeros(2*K, 4);
for k = 1:K
  A(2*k+(-1:0),:)=[xn(1,k)*Pcomp(3,:,k)-Pcomp(1,:,k); xn(2,k)*Pcomp(3,:,k)-Pcomp(2,:,k)];
end
% A = normx(A')';
[dummy,dummy,X] = svd(A,0);
X = X(:,end);
X=X/X(4);

repo: ekfmonoslam


text:
function this = init_processing_step(this, module, varargin)
% initializes next processing step by creating folders for version tracking,
% shuffling data, and updating processing parameters
%
%   MrSeries = init_processing_step(MrSeries, module)
%
% This is a method of class MrSeries.
%
% IN
%   module      'realign', 'smooth', ...
%
% OUT
%
%   side effects:
%   new folder (with current data):
%       dirObject/<nProcessingSteps+1>_moduleName
%   parameters.processingLog
%   nProcessingSteps
%
% EXAMPLE
%   init_processing_step
%
%   See also MrSeries

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-01
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% NOTE: for each new processing step added here, it has to be decided which
% (input, raw, unprocessed) files are saved additionally

itemsSave = this.parameters.save.items;
doSave = ~strcmpi(itemsSave, 'none');
doSaveRaw = ismember(itemsSave, {'all'});
doSaveNifti = ismember(itemsSave, {'nii', 'all', 'processed'});
doSaveObject = ismember(itemsSave, {'object', 'all', 'processed'});

% set file-saving behavior of MrImage to keep disk files
this.data.parameters.save.keepCreatedFiles = ...
    1 ; % keeps files here, cleanup will happen in finish_processing_step

pathSaveRoot = this.parameters.save.path;

% save initial, unprocessed data
isFirstProcessingStep = ~this.nProcessingSteps;
if isFirstProcessingStep && doSave
    dirProcessing = sprintf('%03d_%s', this.nProcessingSteps, 'unprocessed');
    pathProcessing = fullfile(pathSaveRoot, dirProcessing);
    mkdir(pathProcessing);
    this.data.parameters.save.path = pathProcessing;
    this.data.parameters.save.fileName = 'data.nii';

    % save data (MrImage file)
    if doSaveNifti
        this.data.save();
    end

    % strip and save object as well
    if doSaveObject
        fileObject = fullfile(pathProcessing, 'MrObject.mat');
        MrObject = this.copyobj('exclude', 'data'); % copies object without data
        save(fileObject, 'MrObject');
    end
end

% specify new directory to save data here
this.nProcessingSteps   = this.nProcessingSteps + 1;
dirProcessing           = sprintf('%03d_%s', this.nProcessingSteps, module);
pathProcessing          = fullfile(pathSaveRoot, dirProcessing);

this.processingLog{end+1,1} = dirProcessing;

% module-specific adaptations, e.g. data copying

% for all matlabbatches where additional spm output files are saved
hasMatlabbatch = ismember(module, this.get_all_matlabbatch_methods());

% for all matlabbatches, where data is needed as raw.nii before job start
doesNeedDataNifti = ismember(module, {'realign', 'smooth', ...
    'specify_and_estimate_1st_level'});

if doSave || hasMatlabbatch
    mkdir(pathProcessing);
end

if doesNeedDataNifti % data has to be written to disk before running spm_jobman, prepare file-names!
    this.data.parameters.save.path = pathProcessing;
end

switch module
    case 'specify_and_estimate_1st_level'
        % save raw data
        this.data.save;

    case 'analyze_rois'
        % dummy image for path transfer
        inputDummyImage = varargin{1};
        inputDummyImage.parameters.save.path = pathProcessing;

    case 'compute_masks'
        inputImages = varargin{1};
        nImages = numel(inputImages);

        % set paths, save raw input files and prefix file names with "mask"... for all input files
        for iImage = 1:nImages
            inputImages{iImage}.parameters.save.path = pathProcessing;
            inputImages{iImage}.parameters.save.fileName = ...
                tapas_uniqc_prefix_files(inputImages{iImage}.parameters.save.fileName, ...
                'mask', 0, 1);
            inputImages{iImage}.save;
        end

    case 'compute_stat_images'

        % determine output file names, e.g. mean.nii and path of processing
        % directory
        [handleImageArray, nameImageArray] = this.get_all_image_objects('stats');
        for iImage = 1:numel(handleImageArray)
            handleImageArray{iImage}.parameters.save.path = pathProcessing;
            handleImageArray{iImage}.parameters.save.fileName = ...
                [nameImageArray{iImage} '.nii'];
        end

    case 'compute_tissue_probability_maps'

        % adjust path of input image to make it save-able
        inputImage = varargin{1};
        inputImage.parameters.save.path = pathProcessing;
        inputImage.parameters.save.keepCreatedFiles = 'all';

    case 'coregister'
        transformedImage = varargin{1};
        equallyTransformedImages = varargin{2};
        inputImages = [{transformedImage};...
            equallyTransformedImages];

        % set save path
        nImages = numel(inputImages);
        for iImage = 1:nImages
            inputImages{iImage}.parameters.save.path = ...
                pathProcessing;
            inputImages{iImage}.parameters.save.keepCreatedFiles = 'all';
        end

    case 'realign'

    case 'smooth'

        % set file names and save path for statistical images
    case 't_filter'
        this.data.parameters.save.path = pathProcessing;
        this.data.parameters.save.keepCreatedFiles = 'all';
end

end

author: translationalneuromodeling

repo: tapas


text:
function val = ft_getopt(opt, key, default, emptymeaningful)

% FT_GETOPT gets the value of a specified option from a configuration structure
% or from a cell-array with key-value pairs.
%
% Use as
%   val = ft_getopt(s, key, default, emptymeaningful)
% where the input values are
%   s               = structure or cell-array
%   key             = string
%   default         = any valid MATLAB data type (optional, default = [])
%   emptymeaningful = boolean value (optional, default = false)
%
% If the key is present as field in the structure, or as key-value pair in the
% cell-array, the corresponding value will be returned.
%
% If the key is not present, ft_getopt will return the default, or an empty array
% when no default was specified.
%
% If the key is present but has an empty value, then the emptymeaningful flag
% specifies whether the empty value or the default value should be returned.
% If emptymeaningful==true, then the empty array will be returned.
% If emptymeaningful==false, then the specified default will be returned.
%
% See also FT_SETOPT, FT_CHECKOPT

% Copyright (C) 2011-2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargin<3
  default = [];
end

if nargin < 4
  emptymeaningful = false;
end

if isa(opt, 'struct') || isa(opt, 'config')
  % get the key-value from the structure
  fn = fieldnames(opt);
  if ~any(strcmp(key, fn))
    val = default;
  else
    val = opt.(key);
  end

elseif isa(opt, 'cell')
  % get the key-value from the cell-array
  if mod(length(opt),2)
    error('optional input arguments should come in key-value pairs, i.e. there should be an even number');
  end

  % the 1st, 3rd, etc. contain the keys, the 2nd, 4th, etc. contain the values
  keys = opt(1:2:end);
  vals = opt(2:2:end);

  % the following may be faster than cellfun(@ischar, keys)
  valid = false(size(keys));
  for i=1:numel(keys)
    valid(i) = ischar(keys{i});
  end

  if ~all(valid)
    error('optional input arguments should come in key-value pairs, the optional input argument %d is invalid (should be a string)', i);
  end

  hit = find(strcmpi(key, keys));
  if isempty(hit)
    % the requested key was not found
    val = default;
  elseif length(hit)==1
    % the requested key was found
    val = vals{hit};
  else
    error('multiple input arguments with the same name');
  end

elseif isempty(opt)
  % no options are specified, return default
  val = default;
end % isstruct or iscell or isempty

if isempty(val) && ~isempty(default) && ~emptymeaningful
  % use the default value instead of the empty input that was specified:
  % this applies for example if you do functionname('key', []), where
  % the empty is meant to indicate that the user does not know or care
  % what the value is
  val = default;
end


text:
function [xMin,fMin,exitCode]=conjugateGradMethod(f,x0,epsilon,deltaTestDist,delta,lineSearchParams,resetEvery,gammaVal,maxIter)
%%CONJUGATEGRADMETHOD Perform unconstrained nonlinear optimization using
%              the conjugate gradient algorithm. The algorithm performs
%              unconstrained minimization of a nonlinear function without
%              having to provide a Hessian matrix. On large problems, this
%              algorithm can be faster than Newton's method, which has a
%              large matrix inversion, and steepest ascent, which usually
%              converges more slowly.
%
%INPUTS: f A handle to the function (and its gradient) over which the
%          minimization is to be performed. The function [fVal,gVal]=f(x)
%          takes the NX1 x vector and returns the real scalar function
%          value fVal and NX1 gradient gVal at the point x.
%       x0 The NX1-dimensional point from which the minimization starts.
%  epsilon The parameter determining the accuracy of the desired solution
%          in terms of the gradient. The function terminates when
%          norm(g) < epsilon*max([1, norm(x)])
%          where g is the gradient. The default if omitted or an empty
%          matrix is passed is 1e-6.
% deltaTestDist The number of iterations back to use to compute the
%          decrease of the objective function if a delta-based convergence
%          test is performed. If zero, then no delta-based convergence
%          testing is done. The default if omitted or an empty matrix is
%          passed is zero.
%    delta The delta for the delta convergence test. This determines the
%          minimum rate of decrease of the objective function. Convergence
%          is determined if (f'-f)<=delta*f, where f' is the value of the
%          objective function f deltaTestDist iterations ago,and f is the
%          current objective function value. The default if this parameter
%          is omitted or an empty matrix is passed is 0.
% lineSearchParams An optional structure whose members specify tolerances
%          for the line search. The parameters are described as in the
%          lineSearch function, except if -1 is passed instead of a
%          parameter structure, then no line search is performed.
% resetEvery In Chapter 2.1 of [1], when handling non-quadratic function,
%          it is suggested that the subgradients be reset every N
%          iterations. By reset, this means that the past direction is
%          discarded and the next step is a steepest descent step. This is
%          the number of iterations after which a reset will occur.
%          resetEvery>=1. The default if omitted or an empty matrix is
%          passed is N.
% gammaVal In Chapter 2.1 of [1], it is suggested that the descent
%          direction algorithm be reset if
%          abs(g'*gPrev)>gammaVal*norm(gPrev)^2
%          where 0<gammaVal<1. The notion is that the previous and next
%          gradient directions should be nearly orthogonal, so if they are
%          not, within some tolerance, then a steepest descent step is
%          taken. The default if this parameter is omitted or an empty
%          matrix is passed is 0.1.
%  maxIter The maximum number of iterations to use for the algorithm. The
%          default if this parameter is omitted or an empty matrix is
%          passed is 4*N.
%
%OUTPUTS: xMin The value of x at the minimum point found. If exitCode is
%              negative, then this might be an empty matrix.
%         fMin The cost function value at the minimum point found. If
%              exitCode is negative, then this might be an empty matrix.
%     exitCode A value indicating the termination condition of the
%              algorithm. Nonnegative values indicate success; negative
%              values indicate some type of failure. Possible values are:
%                  0 The algorithm termiated successfully based on the
%                    gradient criterion.
%                  1 The algorithm terminated successfully based on the
%                    accuracy criterion.
%                 -1 The maximum number of overall iterations was reached.
%                 -2 A non-finite value was encoutered outside the line
%                    search.
%              Other negative values correspond to a failure in lineSearch
%              and correspond to the exitCode returned by the lineSearch
%              function.
%
%The algorithm is implemented based on the description in Chapter 2.1 of
%[1].
%
%EXAMPLE:
%The example is that used in the lineSearch file.
% f=@(x)deal((x(1)+x(2)-3)*(x(1)+x(2))^3*(x(1)+x(2)-6)^4,... %The function
%            [(-6+x(1)+x(2))^3*(x(1)+x(2))^2*(54+8*x(1)^2+x(2)*(-45+8*x(2))+x(1)*(-45+16*x(2)));
%            (-6+x(1)+x(2))^3*(x(1)+x(2))^2*(54+8*x(1)^2+x(2)*(-45+8*x(2))+x(1)*(-45+16*x(2)))]);%And the gradient as the second return.
% %Note that the deal function is used to make an anonymous function have
% %two outputs.
% x0=[0.5;0.25];
% [xMin,fMin,exitCode]=conjugateGradMethod(f,x0)
%The optimum point found is such that sum(xMin) is approximately
%1.73539450 with a minimum function value of approximately -2.1860756.
%
%REFERENCES:
%[1] D. P. Bertsekas, Nonlinear Programming, 3rd ed. Belmont, MA: Athena
%    Science, 2016.
%
%January 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release

    if(nargin<3||isempty(epsilon))
        epsilon=1e-6;
    end

    if(nargin<4||isempty(deltaTestDist))
        deltaTestDist=0;
    end

    if(nargin<5||isempty(delta))
        delta=0;
    end

    %Take any line search parameters given.
    useLineSearch=true;
    if(nargin>=6&&~isempty(lineSearchParams))
        if(isnumeric(lineSearchParams)&&lineSearchParams==-1)
            %If one just wishes to blindly take steps.
            useLineSearch=false;
        elseif(~isstruct(lineSearchParams)||(isnumeric(lineSearchParams)&&lineSearchParams~=-1))
            error('Unknown lineSearchParams value given')
        end
    else
        lineSearchParams=[];
    end

    N=size(x0,1);

    if(nargin<7||isempty(resetEvery))
       resetEvery=N;
    end

    if(nargin<8||isempty(gammaVal))
       gammaVal=0.1;
    end

    if(nargin<9||isempty(maxIter))
        maxIter=4*N;
    end

    x=x0;
    [fVal,gradF]=f(x0);

    %The previous values are not yet set.
    gradFPrev=[];
    dPrev=[];

    if(deltaTestDist>0)
        pastFVals=zeros(deltaTestDist,1);
        pastFVals(1)=fVal;
    end

    stepsSinceReset=0;

    for cutIter=1:maxIter
        if(stepsSinceReset>0)
            if(stepsSinceReset>=resetEvery)
                stepsSinceReset=0;
            else
                gradFPrevMag2=gradFPrev'*gradFPrev;

                if(abs(gradF'*gradFPrev)>gammaVal*gradFPrevMag2)
                    stepsSinceReset=0;
                end
            end
        end

        if(stepsSinceReset>0)
            %If it hasn't just been reset.

            %Equation 2.35
            beta=gradF'*(gradF-gradFPrev)/gradFPrevMag2;

            %Equation 2.34 for the descent direction.
            d=-gradF+beta*dPrev;
        else
            %Steepest descent direction.
            d=-gradF;
        end

        if(useLineSearch)
            %Perform a line search in the given descent direction.
            [xCur,fValCur,gradFCur,~,~,exitCode]=lineSearch(f,x,d,[],[],lineSearchParams);
            if(isempty(xCur))
                xMin=[];
                fMin=[];
                return;
            end
        else%Take a step without a line search.
            xCur=xPrev+d;
        end

        if(any(~isfinite(xCur(:))))
            exitCode=-2;
            xMin=[];
            fMin=[];
            return;
        end

        %Check for convergence based on the gradient.
        if(norm(gradFCur)<epsilon*max([1,norm(xCur)]))
            xMin=xCur;
            fMin=fValCur;
            exitCode=0;
            return;
        end

        %Check for convergence based on the actual function value.
        if(deltaTestDist~=0)
            if(pastFVals(end)-fValCur<=delta*fValCur)
                xMin=xCur;
                fMin=fValCur;
                exitCode=1;
                return;
            end
            pastFVals=circshift(pastFVals,[1,0]);
            pastFVals(1)=fVal;
        end

        x=xCur;
        gradFPrev=gradF;
        gradF=gradFCur;
        dPrev=d;
        stepsSinceReset=stepsSinceReset+1;
    end

    %The maximum number of iterations elapsed without convergence
    xMin=xCur;
    fMin=fValCur;
    exitCode=-1;
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.


text:
% TEST_MAXWELL_CUBE_H_DRCHLT: data function for Dirichlet boundary condition.

function h = test_maxwell_cube_h_drchlt (x, y, z, ind)

  h = zeros ([3, size(x)]);
  switch (ind)
    case 1
      h(2,:,:) = -exp(z) .* cos(x);
      h(3,:,:) = exp(x) .* cos(y);
    case 2
      h(2,:,:) = exp(z) .* cos(x);
      h(3,:,:) = -exp(x) .* cos(y);
    case 3
      h(1,:,:) = exp(z) .* cos(x);
      h(3,:,:) = -sin(y) .* z;
    case 4
      h(1,:,:) = -exp(z) .* cos(x);
      h(3,:,:) = sin(y) .* z;
    case 5
      h(1,:,:) = -exp(x) .* cos(y);
      h(2,:,:) = sin(y) .* z;
    case 6
      h(1,:,:) = exp(x) .* cos(y);
      h(2,:,:) = -sin(y) .* z;
    otherwise
      error ('h_drchlt: unknown reference number')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function fixed_num = p10_fixed_num ( )

%*****************************************************************************80
%
%% P10_FIXED_NUM returns the number of fixed points in problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer FIXED_NUM, the number of fixed points.
%
  fixed_num = 4;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function net = cnn_mnist_init(varargin)
% Initialize a simple TensorNet for MNIST.

net.layers = {} ;
inputModeSize = [4, 8, 8, 4] ;
secondModeSize = [5, 5, 5, 5] ;
ranks = [1, 2, 2, 2, 1] ;
W = tt_rand(secondModeSize.*inputModeSize, length(secondModeSize), ranks, []) ;
W = tt_matrix(W, secondModeSize, inputModeSize) ;
W.core = single(W.core) ;
net.layers{end+1} = struct('type', 'custom', ...
                           'forward', @vl_nntt_forward, ...
                           'backward', @vl_nntt_backward, ...
                           'W', W, ...
                           'weights', {{W.core, zeros(1,1,prod(secondModeSize),'single')}}, ...
                           'learningRate', [1, 2], ...
                           'weightDecay', [1, 0], ...
                           'outHeight', 1, ...
                           'outWidth', 1, ...
                           'outChannels', prod(secondModeSize)) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{0.1*randn(1,1,prod(secondModeSize),10, 'single'), zeros(1, 10, 'single')}}, ...
                           'learningRate', [1, 2], ...
                           'weightDecay', [1, 0], ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'softmaxloss') ;

author: Bihaqo

repo: TensorNet

save_path: github-repos/MATLAB/Bihaqo-TensorNet


text:
function varargout = det(varargin)
%DET (overloaded)
%
% t = DET(X)

switch class(varargin{1})

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        X = varargin{1};
        [n,m] = size(X);
        if n~=m
            error('Matrix must be square.')
        end
        if nargin == 2
            if strcmp(varargin{2},'polynomial')
                varargout{1} = polynomialform(X);
                return
            else
                error('If you use two arguments in @sdpvar/det, the second should be ''polynomial''');
            end
        end
        if n==1
            varargout{1} = X;
            return
        else
            y = yalmip('define','det_internal',reshape(X,[],1));
        end
        varargout{1} = y;

    otherwise
end

function d = polynomialform(X)

n = X.dim(1);
m = X.dim(2);

if n~=m
    error('Matrix must be square.');
else
    switch n
        case 1
            d = X;
        case 2
            % Freakin overloading on multiplication doesn't work. Probalby
            % stupid code...
            Y1.type = '()';
            Y2.type = '()';
            Y3.type = '()';
            Y4.type = '()';
            Y1.subs = {1,1};
            Y2.subs = {2,2};
            Y3.subs = {1,2};
            Y4.subs = {2,1};
            d = subsref(X,Y1)*subsref(X,Y2)-subsref(X,Y3)*subsref(X,Y4);
        otherwise
            d = 0;
            Y.type = '()';
            for i = 1:n
                Y.subs = {i,1};
                xi = subsref(X,Y);
                if ~isequal(xi,0)
                    Y.subs = {[1:1:i-1 i+1:1:n],2:n};
                    subX = subsref(X,Y);
                    if isa(subX,'sdpvar')
                        d = d + (-1)^(i+1)*xi*polynomialform(subX);
                    else
                        d = d + (-1)^(i+1)*xi*det(subX);
                    end
                end
            end
    end
end
% Reset info about conic terms
if isa(d,'sdpvar')
    d.conicinfo = [0 0];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function y=udb10(x);

% y=udb10(x);
%
% Convert from db values to linear power values.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

y=10.^(x/10);

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function [vp,samples,output] = vpsample_vbmc(Ns,Ninit,vp,gp,optimState,options,wide_flag)

if nargin < 7 || isempty(wide_flag); wide_flag = false; end

% Assign default values to OPTIMSTATE
if ~isfield(optimState,'delta'); optimState.delta = 0; end
if ~isfield(optimState,'EntropySwitch'); optimState.EntropySwitch = false; end
if ~isfield(optimState,'Warmup'); optimState.Warmup = ~vp.optimize_weights; end
if ~isfield(optimState,'temperature'); optimState.temperature = 1; end

%% Set up sampling variables and options

% Perform quick sieve to determine good starting point
[vp,~,elcbo_beta,compute_var,NSentK] = ...
    vpsieve_vbmc(Ninit,1,vp,gp,optimState,options);

K = vp.K;
D = vp.D;

% Compute soft bounds for variational parameters optimization
[vp,thetabnd] = vpbounds(vp,gp,options,K);

% Move lower bound on scale - we want *wider* distributions
if wide_flag
    lnscale = bsxfun(@plus,log(vp.sigma(:))',log(vp.lambda(:)));
    if vp.optimize_mu; idx = D*K; else; idx = 0; end
    thetabnd.lb(idx+1:idx+K*D) = lnscale;
end

%% Sample variational posterior starting from current

theta0 = get_vptheta(vp)';
Ntheta = numel(theta0);

% MCMC parameters
Widths = 0.5;
sampleopts.Thin = 1;
sampleopts.Burnin = 0;
sampleopts.Display = 'off';
sampleopts.Diagnostics = false;
LB = -Inf(1,Ntheta);
UB = Inf(1,Ntheta);

idx_fixed = false(size(theta0));
if ~optimState.Warmup && 0
    if vp.optimize_mu; idx_fixed(1:D*K) = true; end
%    idx_fixed = true(size(theta0));
%    idx_fixed(idx+1:idx+K) = false;
end

LB(idx_fixed) = theta0(idx_fixed);
UB(idx_fixed) = theta0(idx_fixed);

% Perform sampling
try
    switch lower(options.VariationalSampler)
        case 'slicesample'
            vpmcmc_fun = @(theta_) -negelcbo_vbmc(theta_,elcbo_beta,vp,gp,NSentK,0,compute_var,0,thetabnd);
            [samples,fvals,exitflag,output] = ...
                slicesample_vbmc(vpmcmc_fun,theta0,Ns,Widths,LB,UB,sampleopts);
        case 'malasample'
            if isfield(optimState,'mcmc_stepsize')
                sampleopts.Stepsize = optimState.mcmc_stepsize;
                output.stepsize = sampleopts.Stepsize;
            end
            vpmcmc_fun = @(theta_) vpmcmcgrad_fun(theta_,elcbo_beta,vp,gp,NSentK,compute_var,thetabnd);
            [samples,fvals,exitflag,output] = ...
                malasample_vbmc(vpmcmc_fun,theta0,Ns,Widths,LB,UB,sampleopts);
            % output.accept_rate
    end
catch
    samples = repmat(theta0,[Ns,1]);
end
vp = rescale_params(vp,samples(end,:));

end

function [logp,dlogp] = vpmcmcgrad_fun(theta,elcbo_beta,vp,gp,NSentK,compute_var,thetabnd)
    [nlogp,ndlogp] = negelcbo_vbmc(theta,elcbo_beta,vp,gp,NSentK,1,compute_var,0,thetabnd);
    logp = -nlogp;
    dlogp = -ndlogp;
end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function [ cnt ] = opt_eegStruct( dat, field )
% opt_eegStruct:
%
% Description:
%     This function converts the eeg data into a struct in OpenBMI format.
%
% Example:
%    marker={'1','right';'2','left';'3','foot'};
%    field={'x','t','fs','y_dec','y_logic','y_class','class', 'chan'};
%    [EEG.data, EEG.marker, EEG.info]=Load_EEG(file,{'device','brainVision';'marker',marker;'fs',fs});
%    CNT=opt_eegStruct({EEG.data, EEG.marker, EEG.info}, field);
%
% Input:
%     dat   - EEG data structure in a cell type
%     field - Selected necessary fields you want to get in a cell
%
% Output:
%     cnt   - Raw EEG data in OpenBMI format with fields of input 'field'
%
% Min-ho Lee
% mhlee@image.korea.ac.kr
%

if isempty(field)
    error('OpenBMI: There''s no field information')
elseif ~prod(ismember(field,{'x','t','fs','y_dec','y_logic','y_class','class', 'chan'}))
    error('OpenBMI: Unacceptable field name')
end

for i=1:length(field)
    [cnt.(field{i})]=[];
end

stc=struct;
if iscell(dat)
    nC=length(dat);
    for i=1:nC
        t_stc=opt_selectField(dat{i},field);
        stc=opt_catStruct(stc,t_stc);
    end
end

for i=1:length(field)
    if ~isempty(stc.(field{i}))
        cnt(:).(field{i})= stc.(field{i});
    else
        warn_str=sprintf('OpenBMI: The field "%s" is empty', field{i});
        disp(warn_str);
    end
end

end

author: PatternRecognition

repo: OpenBMI