text:
function val = viewGetVolume(vw,param,varargin)
% Get data from various view structures
%
% This function is wrapped by viewGet. It should not be called by anything
% else other than viewGet.
%
% This function retrieves information from the view that relates to a
% specific component of the application.
%
% We assume that input comes to us already fixed and does not need to be
% formatted again.

if notDefined('vw'), vw = getCurView; end
if notDefined('param'), error('No parameter defined'); end

mrGlobals;
val = [];

switch param

    case 'nodes'
        % Return the array of nodes. Only gray views have nodes. See help
        % for mrManDist.m for a description of the node structure. In
        % brief, nodes are 8 x nvoxels. The first 3 rows correspond to the
        % voxel location and the next 5 correspond to gray graph-related
        % data.
        %   nodes = viewGet(vw, 'nodes');
        if isfield(vw, 'nodes'), val = vw.nodes;
        else
            val = [];
            warning('vista:viewError', 'Nodes not found.');
        end
    case 'xyznodes'
        % Return the xyz coordinates of the gray voxels as found in nodes
        % array. Assumes a Gray view. See case 'nodes' and help for
        % mrManDist for more information.
        %
        % Must call this sagittal, axial coronal or whatever the mapping is
        % ras, 06/07 -- I believe it's [cor axi sag]. coords is [axi cor sag].
        %
        % xyzNodes = viewGet(vw, 'xyz nodes');
        nodes = viewGet(vw,'nodes');
        val = nodes(1:3,:);
    case 'nodegraylevel'
        % Return the gray level of each voxel as determined by the nodes
        % array. Assumes a Gray view. See case 'nodes' and help for
        % mrManDist for more information.
        %   nodeGrayLevel = viewGet(vw, 'gray level');
        nodes = viewGet(vw,'nodes');
        val = nodes(6,:);
    case 'nnodes'
        % Return the number of nodes. Assumes a Gray view. See case 'nodes'
        % and help for mrManDist for more information.
        %   nNodes = viewGet(vw, 'number of nodes');
        val = size( viewGet(vw, 'nodes'), 2 );
    case 'edges'
        % Return the edge structure of the gray graph. Assumes a Gray view.
        % See help for mrManDist for more information.
        %   edges = viewGet(vw, 'edges');
        if isfield(vw, 'edges'), val = vw.edges;
        else val = []; warning('vista:viewError', 'Edges not found.'); end
    case 'nedges'
        % Return the number of edges in the gray graph. Assumes a Gray
        % view. See case 'edges' and help for mrManDist for more
        % information.
        %   nEdges = viewGet(vw, 'number of edges');
        val = length( viewGet(vw, 'edges') );
    case 'allleftnodes'
        % Return the subset of nodes in the Gray graph that are in the left
        % hemisphere. See mrgGrowGray and mrManDist.
        %   allLeftNodes = viewGet(vw, 'all left nodes');
        val = vw.allLeftNodes;
    case 'allleftedges'
        % Return the subset of edges in the Gray graph that are in the left
        % hemisphere. See mrgGrowGray and mrManDist.
        %   allLeftEdges = viewGet(vw, 'all left edges');
        val = vw.allLeftEdges;
        if isempty(val) && Check4File('Gray/coords')
            % Try laoding from the Gray/coords file
            load('Gray/coords', 'allLeftEdges')
            val = allLeftEdges;
        end
    case 'allrightnodes'
        % Return the subset of nodes in the Gray graph that are in the
        % right hemisphere. See mrgGrowGray and mrManDist.
        %   allRightNodes = viewGet(vw, 'all right nodes');
        val = vw.allRightNodes;
        if isempty(val) && Check4File('Gray/coords')
            % Try laoding from the Gray/coords file
            load('Gray/coords', 'allRightNodes')
            val = allRightNodes;
        end

    case 'allrightedges'
        % Return the subset of edges in the Gray graph that are in the
        % right hemisphere. See mrgGrowGray and mrManDist.
        %   allRightEdges = viewGet(vw, 'all right edges');
        val = vw.allRightEdges;

    case 'allnodes'
        % Return all nodes from Gray graph by taking union of allLeftNodes
        % and allRightNodes.
        %
        % This is NOT necessarily the same as simply returning 'vw.nodes'.
        % When we install a segmentation, we can either keep all the nodes
        % in the gray graph, or only those that fall within the functional
        % field of view (to save space). When we do the latter, the fields
        % vw.coords, vw.nodes, and vw.edges contain only the coords, nodes,
        % and eges within the functional field of view. However the fields
        % vw.allLeftNodes, vw.allLeftEdges, vw.allRightNodes, and
        % vw.allRightEdges contain the edges and nodes for the entire
        % hemisphere
        %
        % Example: nodes = viewGet(vw, 'all nodes');
        val = [vw.allLeftNodes'; vw.allRightNodes']';

    case 'alledges'
        % Return all edges from Gray graph by taking union of allLeftEdges
        % and allRightEdges. See 'allnodes' for explanation.
        %
        % Example: edges = viewGet(vw, 'all edges');
        val = [vw.allLeftEdges vw.allRightEdges];

    case 'coords'
        % Return all the coordinates in the current view. If in Flat view,
        % return the coordinates for a particular slice (slice specified in
        % varargin{1}). If in Inplane view, slice specification is
        % optional. If in Gray or Volume view, slice specification is
        % ignored.
        %   <gray, volume or inplane>
        %       coords = viewGet(vw, 'coords');
        %   <flat or inplane>
        %       slice  = viewGet(vw, 'current slice');
        %       coords = viewGet(vw,'coords', slice);
        try
            switch lower(viewGet(vw, 'viewType'))
                case 'flat'
                    %% separate coords for each flat hemisphere
                    if length(varargin) ~= 1
                        error('You must specify which hemisphere.');
                    end
                    hname = varargin{1};
                    switch hname
                        case 'left'
                            val = vw.coords{1};
                        case 'right'
                            val = vw.coords{2};
                        otherwise
                            error('Bad hemisphere name');
                    end
                case {'gray', 'volume', 'hiddengray'}
                    val = vw.coords;
                case 'inplane'
                    % These coords are for inplane anat. Functional coords
                    % may have different values (if
                    % upSampleFactor(vw,scan) ~= 1)
                    dims = viewGet(vw, 'anatomysize');
                    if length(varargin) >= 1  % then we want coords for just one slice
                        slice = varargin{1};
                        indices = 1+prod([dims(1:2) slice-1]):prod([dims(1:2) slice]);
                        val=indices2Coords(indices,dims);
                    else
                        indices = 1:prod(dims);
                        val=indices2Coords(indices,dims);
                    end
            end
        catch ME
            val=[];
            warning(ME.identifier, ME.message);
            fprintf('[%s]: Coords not found.', mfilename);
        end

    case 'allcoords'
        % Return all coords from Gray graph, including those that are not
        % included in the functional field of view. See 'allnodes' for
        % explanation. If session was initialized with the option
        % 'keepAllNodes' == true, then this call will be identical to
        % viewGet(vw.coords).
        %
        % Optional input: hemi ('right' 'left' 'both') for right, left or both
        %
        % Example: coords = viewGet(vw, 'allcoords', 'right');
        nodes = viewGet(vw, 'nodes');

        if length(varargin) ~= 1
            hname = 'both';
        else
            hname = varargin{1};
        end

        switch hname
            case 'left'
                subnodes  = viewGet(vw, 'all left nodes');
                [~, ia] = intersectCols(nodes(1:3,:), subnodes(1:3,:));
            case 'right'
                subnodes  = viewGet(vw, 'all right nodes');
                [~, ia] = intersectCols(nodes(1:3,:), subnodes(1:3,:));
            otherwise
                ia = 1:size(nodes,2);
        end

        val = nodes([2 1 3], ia);

    case 'coordsfilename'
        % Return the path to the file in which coordinates are stored.
        % Assumes that a gray view has been created (though current view
        % can be any type).
        %   coordsFileName = viewGet(vw, 'coords file name');
        homeDir = viewGet(vw, 'homedir');
        if isempty(homeDir), val = ['Gray' filesep 'coords.mat'];
        else val = [homeDir filesep 'Gray' filesep 'coords.mat'];
        end
    case 'ncoords'
        % Return the number of coordinates in the current view. See case
        % 'coords'.
        %   nCoords = viewGet(vw, 'number of coords');
        val = size( viewGet(vw, 'Coords'), 2 );

    case 'classfilename'
        % Return the path to either the left or the right gray/white
        % classification file.
        %   fname = viewGet(vw, 'class file name', 'left');
        %   fname = viewGet(vw, 'class file name', 'right');
        if (length(varargin) == 1), hemisphere = varargin{1};
        else  hemisphere = 'left';
        end
        switch lower(hemisphere)
            case {'left' 'l'}
                if ~checkfields(vw,'leftClassFile') || isempty(vw.leftClassFile);
                    [~,val] = GetClassFile(vw, 0, 1);
                else
                    val = vw.leftClassFile;
                end
            case {'right' 'r'}
                if ~checkfields(vw,'rightClassFile') || isempty(vw.rightClassFile)
                    [~,val] = GetClassFile(vw, 1, 1);
                else
                    val = vw.rightClassFile;
                end
            otherwise
                error('Unknown hemisphere');
        end

    case {'classdata','class','classification','whitematter'}
        % classFileRight = viewGet(vw,'class data','right');
        if length(varargin) == 1, hemisphere = varargin{1};
        else error('You must specify right/left hemisphere.');
        end
        switch lower(hemisphere)
            case 'left'
                val = GetClassFile(vw, 0);
            case 'right'
                val = GetClassFile(vw, 1);
            otherwise
                error('Unknown hemisphere');
        end

    case {'graymatterfilename','graypath','grayfilename','grayfile'}
        % grayFile = viewGet(vw,'Gray matter filename','right');
        if length(varargin) == 1, hemisphere = varargin{1};
        else error('You must specify right/left hemisphere.');
        end
        switch lower(hemisphere)
            case 'left'
                if checkfields(vw,'leftPath')
                    val = vw.leftPath;
                end
            case 'right'
                if checkfields(vw,'rightPath')
                    val = vw.rightPath;
                end
            otherwise
                error('Unknown hemisphere');
        end

    otherwise
        error('Unknown viewGet parameter');

end

return

author: vistalab

repo: vistasoft


text:
function  [cct,  ta1, ta2] = srcint(la, k2, ct4, pp1, pp2, ctc0, nprimtilt, lmax, K, stfname, logstep)
% srcpint.m integrates the source including polarisation over conformal time using the approach of Seljak :
% integration over the photon past light cone, see M Zaldarriaga et al., ApJ nr 494, 491 (1998)
% the function needs an l-range (la), a k-range (k2) a ctime-range (ct4),
% ctc0 : the present conformal time, a series of splines of the source in conformal time
% for the temperature (pp1) and for the polarisation (pp2),
% the index of the primary power spectrum (nprimtilt), the maximum la of the anistropy spectrum (lmax)
% and the curvature energy content from the Friedman equation (K), the startfile name of the
% ultra-spherical function (stfname), an option for logaritmic variation of the
% k2 vector : default logstep = 0 (constant steps in k2). k2 has 10-log steps, if logstep = 1.
%
% the result (cct) is the square of the temperature spectrum (1) resp the E-polarisation spectrum (2)
% resp the cross-correlation of temperature and polarisation (3).
%
% D Vangheluwe 4 april 2005
% remark 1: we use ultra-spherical bessel values which have two parameters and must be found by integration:
%   for the calculation of the ultra-spherical bessel functions see: cmb/usphint.m
% remark 2: for the ode equation of the u-function, see equation (36) of Zaldarriaga & Seljak.
% remark 3: take attention the value of k2 should lie within the range of kb1 values : there is no check!!
%   see spline of start  values x10 and pbd0.
% remark 4: how do we get the ultra-spherical bessel function values? see development and test routine usphpar.m

lla = size(la, 2);
lk2 = size(k2, 2);
k2min = k2(1);
k2max = k2(lk2);
lct4 = size(ct4, 2);
kk0 = sign(K);

% load the table and find the start values for the integration of the ode for u-functions
% the startfile can be obtained by running usphst.m
%  stv1 = load('usphst.dat');
stv1 = load(stfname);
lmax0 = stv1.ust.la(end);

ek2step = 1;
if logstep == 1, ek2step = log(k2(2)/k2(1))/log(10); end
nrsteps = 10000;   % default 10000
%lmax = 1500;
xmax = k2(lk2) * ctc0;   % =default 3000
xstep = 2*lmax0/nrsteps;
kx = 1e-12 : xstep : xmax;

% allocate the data in order to speed up the routine
kctc0 = zeros(1, lct4);
kctc = zeros(1, lct4);
jl = zeros(lct4, 1);
src = zeros(lct4, 1);
htable = zeros(lct4, 1);
tetatl = zeros(1, lk2);
tetael = zeros(1, lk2);
cmbtable = zeros(1, lk2);

% calculate the curvature length times wavenumber : a parameter for the ultra_spherical bessel function
kb = sqrt(-K) ./k2;

%#########
kb(1)
kb(end)
if all(kb == 0), kbzero = 1; x10 = 1e-12 * ones(1, lk2);
else  % not all kb are zero

  kbzero = 0;
  ikb = find(kb > 10);
  if ~isempty(ikb),
      kb(ikb) = 10 * ones(1, size(ikb, 2));
      message('boundary kb > 10 reached')
  end

% split the structure stv1 from the startfile : x1: start values where pbi_beta = 1e-6,
% pbd : dphi_beta/dx at the start values, all data exact within 1e-6.
% the la-values resp kb-vector should be the same as in the program usphst.m
  la1 = stv1.ust.la;
  lla1 = size(la1, 2);
  kb1 = stv1.ust.kb;
  lkb1 = size(kb1, 2);
  x1v1 = stv1.ust.x1;
  pbdv1 = stv1.ust.pbd;

% find all start values (hs) for the ode integration step as a function of la and kb (is a surface):
% the long sought magic formula : la1 is the la-vector from usphst.m!!
  hs = sqrt(0.6* x1v1 ./ repmat(((la1+250) .^0.905)', 1, lkb1));
  hs = hs.* repmat((kb1 + 1) ./ (4*kb1 + 1), lla1, 1);

end %all(kb == 0)

%############################# start
for il = 1:lla

    l = la(il)

    if ~kbzero
       il1 = find(la1 == l);
       if isempty(il1)
           error('la value not found, try another value')
           break;
       else
           x10 = spline(kb1, x1v1(il1,:), kb);
           pbd0 = spline(kb1, pbdv1(il1,:), kb);
       end

% calculate the start values for the u-function : u0 = [u, du/dx] with u = r(x) * phi_beta(x) and
%  du/dx = phi_beta(x) * dr(x)/dx + r(x) * dphi_beta(x)/dx :
       r10 = sinh(kb .* x10) ./kb;
       rd10 = cosh(kb .* x10);
       u0 = [r10 *1e-6; r10 .* pbd0 + 1e-6 * rd10];

% define the step and set 'odeint' to a constant number of steps for all cases (kb)
       hstart = spline(kb1, hs(il1,:), kb);
       toi = 1e-4;
       hmax = 0.3;   %we take hmax = 0.3 as the default value
       maxstp = 150;   % default maxstp = 150
       x1e = ones(1, lk2) * xmax;
% solve the u-function instead of the phi_beta function : xs2 does not have a constant step
       [xs2, u2, nok, nbad, nfev] = odeintp(u0, x10, x1e, toi, hstart, 0, hmax, maxstp, @uspheq1, l, kb, kk0);

% prepare a spline of the solution for a constant step in xs (the values are found with ppval1):
       y2der = splinep(xs2, u2);
% calculate the number of constant steps, xstep in the solution xs2: ilast is the last one
       ilast = ceil((xs2(end,:) - x10)/xstep);
       ncsteps = 8;
       ilast = ilast - ncsteps;
% set the number of overlapping steps for the matching to the asymptotic solution (default=40):
       noverflow = 40;
%       noverflow = 80;  %used for smaller xstep
    end  % kbzero

% calculate also the spherical bessel function as we may need it for kb < 1-5
    table_bv = sphbes(l, kx);
% set some constants needed for the polarisation formula
    gl = sqrt((l + 2) * (l + 1) * l * (l - 1));

% interpolate and integrate (sum with the Simpson rule) over conformal time :
    for i = 1:lk2
%    for i = 967:967

       if kb(i) > 1e-5 * (1500/l)

% define a vector of argument values for the ultra-spherical bessel function (x1): reverse ct4
          clear('xs', 'x1', 'x2', 'yspl2','iover','iasym', 'yas', 'table_pb')
          x1 = x10(i) : xstep : xmax;
          lx1 = size(x1, 2);

% in the next 10 lines make a table over x1 (table_pb) of ultra-spherical bessel values:

% make the interpolation of the ode solution for a constant step xstep:
          [xs, yspl2] = ppval1(xs2(:,i), u2(:,i), y2der(:,i), xstep);

% define the overflow index in x1 (iover):
%default 40 steps should at least include one zero and one maximum of u2
          iover = (ilast(i) - noverflow) : ilast(i);
% define the index of x1 where the asymptotic appr is applied (iasym)
          iasym = ilast(i) : lx1;

% calculate the ultra-spherical bessel values
          kctc0 = ctc0 * k2(i) * ones(1, lct4);
          kctc = k2(i) * ct4;
          xctc = kctc0 - kctc;
          ixode = find(xctc > x1(1)  &  xctc < x1(ilast(i)));
          ul = zeros(1, lct4);
          if  ~isempty(ixode)
             r1 = sinh(kb(i) * xs)/kb(i);
% interpolate the table of ultra-spherical bessel function in the range of k values of the ode solution
             ul(ixode) = interpl(xs, yspl2 ./ r1, xctc(ixode));
% calculate the necessry values of the asymptotic expansion of the ultra-spherical bessel function
             ixas = find(xctc >= x1(ilast(i)));
% find the asymptotic values and the phase correction in the overflow region:
             if ~isempty(ixas)
                 dphi = phdif(x1(iover), usphas(l, kb(i), x1(iover), 0, 0, kk0), yspl2(iover));
                 ul(ixas) = usphas1(l, kb(i), xctc(ixas), dphi, 1, kk0);
             end
% integrate the source function over ct4, following formula (40) of Zaldarriaga et al.(1998)
             src = ppval(pp1, k2(i));

figure(1)
ix5= ixode;
plot(ct4(ix5), src(ix5)', ct4(ix5), ul(ix5),'--')
% as we have an integration over tau (not k*tau) and the ultra-spherical bessel function is normalised on k2,
% we have to multiply the integrant with k2 (then it agrees with the flat case, see kb < 1e-5)
             htable = src .* ul';
             tetatl(i) = simpsint(ct4(1), ct4(end), htable);
%plot(ct4, src', ct4, ul,'--', ct4, htable)
%a1= simpsint(ct4(1), ct4(end), src)
%tetatl(i)
             src = ppval(pp2, k2(i));
             htable = src .* ul';
             tetael(i) = simpsint(ct4(1), ct4(end), htable);
          else
             tetatl(i) = 0;
             tetael(i) = 0;
          end

       else  %if kb(i) <= 1e-5 *(1500/l)

  % interpolate the spherical bessel function
          kctc0 = ctc0 * k2(i) * ones(1, lct4);
          kctc = k2(i) * ct4;
          xctc = kctc0 - kctc;
          jl = interpl(kx, table_bv, xctc)';
% integrate the source function following (12) and (13) of Seljak resp (18) of Zaldarriaga and Seljak
% the integration interval is limited to ct4 where the source <> 0
          src = ppval(pp1, k2(i));
          htable = src .* jl;
          tetatl(i) = simpsint(ct4(1), ct4(end), htable);
          src = ppval(pp2, k2(i));
          htable = src .* jl;
          tetael(i) = simpsint(ct4(1), ct4(end), htable);
       end  %if kb(i)

   end  % forloop k2

ta1 = tetatl;
ta2 = tetael;
factor = 1;
f2 = 1;
%if (K ~= 0),  f2 = coth(pi*k2/sqrt(-K)); end
if  logstep == 1, factor = k2 * log(10) * ek2step * (lk2-1)/(k2max - k2min); end
% perform the integration over k2min-k2max, following (9) of Seljak resp (19) of Zaldarriaga and Seljak
%   cmbtable = (tetatl .^2) ./ (k2 .^ (2-nprimtilt));
   cmbtable = factor .* f2 .* (tetatl .^2) .* (k2 ./ (k2 .^2 - K)) .^ (2 - nprimtilt);
   cct.ctt(il) = l*(l + 1) * simpsint(k2min, k2max, cmbtable');
%   cmbtable = (tetael .^2) ./ (k2 .^ (2-nprimtilt));
   cmbtable = factor .* (tetael .^2) .* (k2 ./ (k2 .^2 - K)) .^ (2 - nprimtilt);
   cct.cee(il) = gl^2 * l*(l + 1) * simpsint(k2min, k2max, cmbtable');
%   cmbtable = (tetael .* tetatl) ./ (k2 .^ (2-nprimtilt));
   cmbtable = factor .* (tetael .* tetatl) .* (k2 ./ (k2 .^2 - K)) .^ (2 - nprimtilt);
   cct.cte(il) = l*(l + 1) * gl * simpsint(k2min, k2max, cmbtable');

end  % forloop la

author: Sable

repo: mcbench-benchmarks


text:
function k = whitefixedKernDiagCompute(kern, x)

% WHITEFIXEDKERNDIAGCOMPUTE Compute diagonal of WHITEFIXED kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the fixed parameter white noise kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : whitefixedKernParamInit, kernDiagCompute, kernCreate, whitefixedKernCompute
%
% COPYRIGHT : Nathaniel J. King, 2006

% KERN

k = whiteKernDiagCompute(kern, x);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% teststat - EEGLAB statistical testing function
%
% Statistics are critical for inference testing in Science. It is thus
% primordial to make sure than all the statistics implemented are
% robust and at least bug free. Statistical function using complex
% formulas are inherently prone to bugs. EEGLAB functions are all the
% more prone to bugs given that they only use complex Matlab code to
% avoid loops and speed up computation.
%
% This test function does not garantee that EEGLAB statistical functions
% are bug free. It does assure though that bugs are unlikely and minor
% if they are present.
%
% This function test 3 things.
%
% * First, it checks that for vector inputs the EEGLAB functions return
%   the same output as other reference functions from the Matlab statistical
%   toolbox or from other packages tested against the SPSS software for
%   repeated measure ANOVA (rm_anova2 function).
%
% * Second, it checks that array inputs with different number of dimensions
%   (from 1 to 3) the EEGLAB function return the same output.
%
% * Third, it checks that the permutation and bootstrap methods shuffle
%   the data properly by running multiple tests.

function teststat;

% testing paired t-test
% ---------------------
a = { rand(1,10) rand(1,10)+0.5 };
[t df pvals surog] = statcond(a, 'mode', 'param', 'verbose', 'off', 'paired', 'on');
[h p tmp stats] = ttest(a{1}, a{2});
fprintf('Statistics paired statcond    t-value %2.2f df=%d p=%0.4f\n', t, df, pvals);
fprintf('Statistics paired ttest func. t-value %2.2f df=%d p=%0.4f\n', stats.tstat, stats.df, p);
assertsame([t stats.tstat], [df stats.df], [pvals p]);
disp('--------------------');

% testing unpaired t-test
% -----------------------
[t df pvals surog] = statcond(a, 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[h p tmp stats] = ttest2(a{1}, a{2});
fprintf('Statistics paired statcond     t-value %2.2f df=%d p=%0.4f\n', t, df, pvals);
fprintf('Statistics paired ttest2 func. t-value %2.2f df=%d p=%0.4f\n', stats.tstat, stats.df, p);
assertsame([t stats.tstat], [df stats.df], [pvals p]);
disp('--------------------');

% testing paired 1-way ANOVA
% --------------------------
a = { rand(1,10) rand(1,10) rand(1,10)+0.2; rand(1,10) rand(1,10)+0.2 rand(1,10) };
[F df pvals surog] = statcond(a(1,:), 'mode', 'param', 'verbose', 'off', 'paired', 'on');
z = zeros(10,1); o = ones(10,1); t = ones(10,1)*2;
stats = rm_anova2(  [ a{1,1}';a{1,2}';a{1,3}'], repmat([1:10]', [3 1]), [o;o;o], [z;o;t], {'a','b'});
fprintf('Statistics 1-way paired statcond        F-value %2.2f df1=%d df2=%d p=%0.4f\n', F, df(1), df(2), pvals);
fprintf('Statistics 1-way paired rm_avova2 func. F-value %2.2f df1=%d df2=%d p=%0.4f\n', stats{3,5}, stats{3,3}, stats{6,3}, stats{3,6});
assertsame([F stats{3,5}], [df(1) stats{3,3}], [df(2) stats{6,3}], [pvals stats{3,6}]);
disp('--------------------');

% testing paired 2-way ANOVA
% --------------------------
[F df pvals surog] = statcond(a, 'mode', 'param', 'verbose', 'off', 'paired', 'on');
z = zeros(10,1); o = ones(10,1); t = ones(10,1)*2;
stats = rm_anova2(  [ a{1,1}';a{1,2}';a{1,3}';a{2,1}';a{2,2}';a{2,3}' ], ...
            repmat([1:10]', [6 1]), [o;o;o;z;z;z], [z;o;t;z;o;t], {'a','b'});
fprintf('Statistics 2-way paired statcond        F-value %2.2f df1=%d df2=%d p=%0.4f\n', F{3}, df{3}(1), df{3}(2), pvals{3});
fprintf('Statistics 2-way paired rm_avova2 func. F-value %2.2f df1=%d df2=%d p=%0.4f\n', stats{4,5}, stats{4,3}, stats{7,3}, stats{4,6});
assertsame([F{3} stats{4,5}], [df{3}(1) stats{4,3}], [df{3}(2) stats{7,3}], [pvals{3} stats{4,6}]);
disp('--------------------');

% testing 1-way unpaired ANOVA
% ----------------------------
[F df pvals surog] = statcond(a(1,:), 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[p stats] = anova1( [ a{1,1}' a{1,2}' a{1,3}' ],{}, 'off');
fprintf('Statistics 1-way unpaired statcond     F-value %2.2f df1=%d df2=%d p=%0.4f\n', F, df(1), df(2), pvals);
fprintf('Statistics 1-way unpaired anova1 func. F-value %2.2f df1=%d df2=%d p=%0.4f\n', stats{2,5}, stats{2,3}, stats{3,3}, stats{2,6});
assertsame([F stats{2,5}], [df(1) stats{2,3}], [df(2) stats{3,3}], [pvals stats{2,6}]);
disp('--------------------');

% testing 2-way unpaired ANOVA
% ----------------------------
[F df pvals surog] = statcond(a, 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[p stats] = anova2( [ a{1,1}' a{1,2}' a{1,3}'; a{2,1}' a{2,2}' a{2,3}' ], 10, 'off');
fprintf('Statistics 2-way paired statcond       F-value %2.2f df1=%d df2=%d p=%0.4f\n', F{3}, df{3}(1), df{3}(2), pvals{3});
fprintf('Statistics 1-way unpaired anova2 func. F-value %2.2f df1=%d df2=%d p=%0.4f\n', stats{4,5}, stats{4,3}, stats{5,3}, stats{4,6});
assertsame([F{3} stats{4,5}], [df{3}(1) stats{4,3}], [df{3}(2) stats{5,3}], [pvals{3} stats{4,6}]);
disp('--------------------');

% testing different dimensions in statcond
% ----------------------------------------
a = { rand(1,10)      rand(1,10)+0.5      rand(1,10)};
b = { rand(10,10)     rand(10,10)+0.5     rand(10,10)};     b{1}(4,:)     = a{1}; b{2}(4,:)     = a{2}; b{3}(4,:)     = a{3};
c = { rand(5,10,10)   rand(5,10,10)+0.5   rand(5,10,10)};   c{1}(2,4,:)   = a{1}; c{2}(2,4,:)   = a{2}; c{3}(2,4,:)   = a{3};
d = { rand(2,5,10,10) rand(2,5,10,10)+0.5 rand(2,5,10,10)}; d{1}(1,2,4,:) = a{1}; d{2}(1,2,4,:) = a{2}; d{3}(1,2,4,:) = a{3};
[t1 df1 pvals1] = statcond(a(1:2), 'mode', 'param', 'verbose', 'off', 'paired', 'on');
[t2 df2 pvals2] = statcond(b(1:2), 'mode', 'param', 'verbose', 'off', 'paired', 'on');
[t3 df3 pvals3] = statcond(c(1:2), 'mode', 'param', 'verbose', 'off', 'paired', 'on');
[t4 df4 pvals4] = statcond(d(1:2), 'mode', 'param', 'verbose', 'off', 'paired', 'on');
fprintf('Statistics paired statcond t-test dim1 t-value %2.2f df=%d p=%0.4f\n', t1, df1, pvals1);
fprintf('Statistics paired statcond t-test dim2 t-value %2.2f df=%d p=%0.4f\n', t2(4), df2, pvals2(4));
fprintf('Statistics paired statcond t-test dim3 t-value %2.2f df=%d p=%0.4f\n', t3(2,4), df3, pvals3(2,4));
fprintf('Statistics paired statcond t-test dim4 t-value %2.2f df=%d p=%0.4f\n', t4(1,2,4), df4, pvals4(1,2,4));
assertsame([t1 t2(4) t3(2,4) t4(1,2,4)], [df1 df2 df3 df4], [pvals1 pvals2(4) pvals3(2,4) pvals4(1,2,4)]);
disp('--------------------');
[t1 df1 pvals1] = statcond(a(1:2), 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[t2 df2 pvals2] = statcond(b(1:2), 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[t3 df3 pvals3] = statcond(c(1:2), 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[t4 df4 pvals4] = statcond(d(1:2), 'mode', 'param', 'verbose', 'off', 'paired', 'off');
fprintf('Statistics unpaired statcond t-test dim1 t-value %2.2f df=%d p=%0.4f\n', t1, df1, pvals1);
fprintf('Statistics unpaired statcond t-test dim2 t-value %2.2f df=%d p=%0.4f\n', t2(4), df2, pvals2(4));
fprintf('Statistics unpaired statcond t-test dim3 t-value %2.2f df=%d p=%0.4f\n', t3(2,4), df3, pvals3(2,4));
fprintf('Statistics unpaired statcond t-test dim4 t-value %2.2f df=%d p=%0.4f\n', t4(1,2,4), df4, pvals4(1,2,4));
assertsame([t1 t2(4) t3(2,4) t4(1,2,4)], [df1 df2 df3 df4], [pvals1 pvals2(4) pvals3(2,4) pvals4(1,2,4)]);
disp('--------------------');
[t1 df1 pvals1] = statcond(a, 'mode', 'param', 'verbose', 'off', 'paired', 'on');
[t2 df2 pvals2] = statcond(b, 'mode', 'param', 'verbose', 'off', 'paired', 'on');
[t3 df3 pvals3] = statcond(c, 'mode', 'param', 'verbose', 'off', 'paired', 'on');
[t4 df4 pvals4] = statcond(d, 'mode', 'param', 'verbose', 'off', 'paired', 'on');
fprintf('Statistics paired statcond anova 1-way dim1 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t1, df1(1), df1(2), pvals1);
fprintf('Statistics paired statcond anova 1-way dim2 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t2(4), df2(1), df2(2), pvals2(4));
fprintf('Statistics paired statcond anova 1-way dim3 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t3(2,4), df3(1), df3(2), pvals3(2,4));
fprintf('Statistics paired statcond anova 1-way dim4 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t4(1,2,4), df4(1), df4(2), pvals4(1,2,4));
assertsame([t1 t2(4) t3(2,4) t4(1,2,4)], [df1 df2 df3 df4], [pvals1 pvals2(4) pvals3(2,4) pvals4(1,2,4)]);
disp('--------------------');
[t1 df1 pvals1] = statcond(a, 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[t2 df2 pvals2] = statcond(b, 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[t3 df3 pvals3] = statcond(c, 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[t4 df4 pvals4] = statcond(d, 'mode', 'param', 'verbose', 'off', 'paired', 'off');
fprintf('Statistics unpaired statcond anova 1-way dim1 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t1, df1(1), df1(2), pvals1);
fprintf('Statistics unpaired statcond anova 1-way dim2 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t2(4), df2(1), df2(2), pvals2(4));
fprintf('Statistics unpaired statcond anova 1-way dim3 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t3(2,4), df3(1), df3(2), pvals3(2,4));
fprintf('Statistics unpaired statcond anova 1-way dim4 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t4(1,2,4), df4(1), df4(2), pvals4(1,2,4));
assertsame([t1 t2(4) t3(2,4) t4(1,2,4)], [df1 df2 df3 df4], [pvals1 pvals2(4) pvals3(2,4) pvals4(1,2,4)]);
disp('--------------------');
a(2,:) = a; a{1} = a{1}/2;
b(2,:) = b; b{1} = b{1}/2;
c(2,:) = c; c{1} = c{1}/2;
d(2,:) = d; d{1} = d{1}/2;
[t1 df1 pvals1] = statcond(a, 'mode', 'param', 'verbose', 'off', 'paired', 'on');
[t2 df2 pvals2] = statcond(b, 'mode', 'param', 'verbose', 'off', 'paired', 'on');
[t3 df3 pvals3] = statcond(c, 'mode', 'param', 'verbose', 'off', 'paired', 'on');
[t4 df4 pvals4] = statcond(d, 'mode', 'param', 'verbose', 'off', 'paired', 'on');
fprintf('Statistics paired statcond anova 2-way dim1 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t1{3}, df1{3}(1), df1{3}(2), pvals1{3});
fprintf('Statistics paired statcond anova 2-way dim2 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t2{3}(4), df2{3}(1), df2{3}(2), pvals2{3}(4));
fprintf('Statistics paired statcond anova 2-way dim3 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t3{3}(2,4), df3{3}(1), df3{3}(2), pvals3{3}(2,4));
fprintf('Statistics paired statcond anova 2-way dim4 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t4{3}(1,2,4), df4{3}(1), df4{3}(2), pvals4{3}(1,2,4));
assertsame([t1{3} t2{3}(4) t3{3}(2,4) t4{3}(1,2,4)], [df1{3}(1) df2{3}(1) df3{3}(1) df4{3}(1)], [df1{3}(2) df2{3}(2) df3{3}(2) df4{3}(2)], [pvals1{3} pvals2{3}(4) pvals3{3}(2,4) pvals4{3}(1,2,4)]);
disp('--------------------');
[t1 df1 pvals1] = statcond(a, 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[t2 df2 pvals2] = statcond(b, 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[t3 df3 pvals3] = statcond(c, 'mode', 'param', 'verbose', 'off', 'paired', 'off');
[t4 df4 pvals4] = statcond(d, 'mode', 'param', 'verbose', 'off', 'paired', 'off');
fprintf('Statistics unpaired statcond anova 2-way dim1 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t1{3}, df1{3}(1), df1{3}(2), pvals1{3});
fprintf('Statistics unpaired statcond anova 2-way dim2 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t2{3}(4), df2{3}(1), df2{3}(2), pvals2{3}(4));
fprintf('Statistics unpaired statcond anova 2-way dim3 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t3{3}(2,4), df3{3}(1), df3{3}(2), pvals3{3}(2,4));
fprintf('Statistics unpaired statcond anova 2-way dim4 t-value %2.2f df1=%d df2=%d p=%0.4f\n', t4{3}(1,2,4), df4{3}(1), df4{3}(2), pvals4{3}(1,2,4));
assertsame([t1{3} t2{3}(4) t3{3}(2,4) t4{3}(1,2,4)], [df1{3}(1) df2{3}(1) df3{3}(1) df4{3}(1)], [df1{3}(2) df2{3}(2) df3{3}(2) df4{3}(2)], [pvals1{3} pvals2{3}(4) pvals3{3}(2,4) pvals4{3}(1,2,4)]);
disp('--------------------');

% testing shuffling and permutation for bootstrap
% -----------------------------------------------
clear a;
m1 = [1:10];
m2 = [1:10]+100;
m3 = [1:10]+1000;
a{1} = { m1 m2 };
a{2} = { m1 m2 m3 };
a{3} = { [ zeros(9,10); m1] [ zeros(9,10); m2] };
a{4} = { [ zeros(9,10); m1] [ zeros(9,10); m2] [ zeros(9,10); m3] };
tmpa = zeros(9,8,10); tmpa(end,end,:) = m1;
tmpb = zeros(9,8,10); tmpb(end,end,:) = m2;
tmpc = zeros(9,8,10); tmpc(end,end,:) = m3;
a{5} = { tmpa tmpb };
a{6} = { tmpa tmpb tmpc };

for method = 1:2
    if method == 2, opt = {'arraycomp', 'off'}; else opt = {}; end;
    for dim = 1:length(a)
        [sa1] = statcond(a{dim}, 'mode', 'bootstrap', 'verbose', 'off', 'paired', 'on',  'returnresamplingarray', 'on', opt{:}, 'naccu', 10);
        [sa2] = statcond(a{dim}, 'mode', 'perm'     , 'verbose', 'off', 'paired', 'on',  'returnresamplingarray', 'on', opt{:}, 'naccu', 10);
        [sa3] = statcond(a{dim}, 'mode', 'bootstrap', 'verbose', 'off', 'paired', 'off', 'returnresamplingarray', 'on', opt{:}, 'naccu', 10);
        [sa4] = statcond(a{dim}, 'mode', 'perm'     , 'verbose', 'off', 'paired', 'off', 'returnresamplingarray', 'on', opt{:}, 'naccu', 10);

        % select data
        nd = ndims(sa1{1});
        if nd == 2 && size(sa1{1},2) > 1
            for t=1:length(sa1),
                sa1{t} = sa1{t}(end,:);
                sa2{t} = sa2{t}(end,:);
                sa3{t} = sa3{t}(end,:);
                sa4{t} = sa4{t}(end,:);
            end;
        elseif nd == 3
            for t=1:length(sa1),
                sa1{t} = squeeze(sa1{t}(end,end,:));
                sa2{t} = squeeze(sa2{t}(end,end,:));
                sa3{t} = squeeze(sa3{t}(end,end,:));
                sa4{t} = squeeze(sa4{t}(end,end,:));
            end;
        elseif nd == 4
            for t=1:length(sa1),
                sa1{t} = squeeze(sa1{t}(end,end,end,:));
                sa2{t} = squeeze(sa2{t}(end,end,end,:));
                sa3{t} = squeeze(sa3{t}(end,end,end,:));
                sa4{t} = squeeze(sa4{t}(end,end,end,:));
            end;
        end;

        % for paired bootstrap, we make sure that the resampling has only shuffled between conditions
        % for instance [101 2 1003 104 ...] is an acceptable sequence
        if all(rem(sa1{1}(:)',10) == [1:9 0]) && all(rem(sa1{2}(:)',10) == [1:9 0])
            fprintf('Bootstrap paired dim%d resampling method %d Pass\n', dim, method);
        else  error('Bootstrap paired resampling Error');
        end;
        % for paired permutation, in addition, we make sure that the sum accross condition is constant
        % which is not true for bootstrap
        msa = meansa(sa2); msa = msa(:)-msa(1);
        if all(rem(sa1{1}(:)',10) == [1:9 0]) && all(rem(sa1{2}(:)',10) == [1:9 0]) && ...
            all(round(msa) == [0:9]') && length(unique(sa2{1})) == 10 && length(unique(sa2{2})) == 10
            fprintf('Permutation paired dim%d resampling method %d Pass\n', dim, method);
        else  error('Permutation paired resampling Error');
        end;
        % for unpaired bootstrap, only make sure there are enough unique
        % values
        if length(unique(sa3{1})) > 3 && length(unique(sa3{2})) > 3
             fprintf('Bootstrap unpaired dim%d reampling method %d Pass\n', dim, method);
        else   error('Bootstrap unpaired reampling Error');
        end;
        % for unpaired permutation, the number of unique values must be 10
        % and the sum must be constant (not true for bootstrap)
        if length(unique(sa4{1})) == 10 && length(unique(sa4{2})) == 10 && ( floor(mean(meansa(sa4))) == 55 || floor(mean(meansa(sa4))) == 372 )
             fprintf('Permutation unpaired dim%d reampling method %d Pass\n', dim, method);
        else   error('Permutation unpaired reampling Error');
        end;

        disp('------------------------');
    end;
end;

% function to check
function assertsame(varargin)

for ind = 1:length(varargin)
    if length(varargin{1}) > 2
        for tmpi = 1:length(varargin{1})-1
            assertsame(varargin{1}(tmpi:tmpi+1));
        end;
        return;
    else
        if (varargin{ind}(1)-varargin{ind}(2)) > abs(mean(varargin{ind}))*0.01
            error('Test failed');
        end;
    end;
end;
disp('Test pass');

function [meanmat] = meansa(mat)

meanmat = zeros(size(mat{1}));
for index = 1:length(mat)
    meanmat = meanmat+mat{index}/length(mat);
end;

function stats = rm_anova2(Y,S,F1,F2,FACTNAMES)
%
% function stats = rm_anova2(Y,S,F1,F2,FACTNAMES)
%
% Two-factor, within-subject repeated measures ANOVA.
% For designs with two within-subject factors.
%
% Parameters:
%    Y          dependent variable (numeric) in a column vector
%    S          grouping variable for SUBJECT
%    F1         grouping variable for factor #1
%    F2         grouping variable for factor #2
%    F1name     name (character array) of factor #1
%    F2name     name (character array) of factor #2
%
%    Y should be a 1-d column vector with all of your data (numeric).
%    The grouping variables should also be 1-d numeric, each with same
%    length as Y. Each entry in each of the grouping vectors indicates the
%    level # (or subject #) of the corresponding entry in Y.
%
% Returns:
%    stats is a cell array with the usual ANOVA table:
%      Source / ss / df / ms / F / p
%
% Notes:
%    Program does not do any input validation, so it is up to you to make
%    sure that you have passed in the parameters in the correct form:
%
%       Y, S, F1, and F2 must be numeric vectors all of the same length.
%
%       There must be at least one value in Y for each possible combination
%       of S, F1, and F2 (i.e. there must be at least one measurement per
%       subject per condition).
%
%       If there is more than one measurement per subject X condition, then
%       the program will take the mean of those measurements.
%
% Aaron Schurger (2005.02.04)
%   Derived from Keppel & Wickens (2004) "Design and Analysis" ch. 18
%

%
% Revision history...
%
% 11 December 2009 (Aaron Schurger)
%
% Fixed error under "bracket terms"
% was: expY = sum(Y.^2);
% now: expY = sum(sum(sum(MEANS.^2)));
%

stats = cell(4,5);

F1_lvls = unique(F1);
F2_lvls = unique(F2);
Subjs = unique(S);

a = length(F1_lvls); % # of levels in factor 1
b = length(F2_lvls); % # of levels in factor 2
n = length(Subjs); % # of subjects

INDS = cell(a,b,n); % this will hold arrays of indices
CELLS = cell(a,b,n); % this will hold the data for each subject X condition
MEANS = zeros(a,b,n); % this will hold the means for each subj X condition

% Calculate means for each subject X condition.
% Keep data in CELLS, because in future we may want to allow options for
% how to compute the means (e.g. leaving out outliers > 3stdev, etc...).
for i=1:a % F1
    for j=1:b % F2
        for k=1:n % Subjs
            INDS{i,j,k} = find(F1==F1_lvls(i) & F2==F2_lvls(j) & S==Subjs(k));
            CELLS{i,j,k} = Y(INDS{i,j,k});
            MEANS(i,j,k) = mean(CELLS{i,j,k});
        end
    end
end

% make tables (see table 18.1, p. 402)
AB = reshape(sum(MEANS,3),a,b); % across subjects
AS = reshape(sum(MEANS,2),a,n); % across factor 2
BS = reshape(sum(MEANS,1),b,n); % across factor 1

A = sum(AB,2); % sum across columns, so result is ax1 column vector
B = sum(AB,1); % sum across rows, so result is 1xb row vector
S = sum(AS,1); % sum across columns, so result is 1xs row vector
T = sum(sum(A)); % could sum either A or B or S, choice is arbitrary

% degrees of freedom
dfA = a-1;
dfB = b-1;
dfAB = (a-1)*(b-1);
dfS = n-1;
dfAS = (a-1)*(n-1);
dfBS = (b-1)*(n-1);
dfABS = (a-1)*(b-1)*(n-1);

% bracket terms (expected value)
expA = sum(A.^2)./(b*n);
expB = sum(B.^2)./(a*n);
expAB = sum(sum(AB.^2))./n;
expS = sum(S.^2)./(a*b);
expAS = sum(sum(AS.^2))./b;
expBS = sum(sum(BS.^2))./a;
expY = sum(sum(sum(MEANS.^2))); %sum(Y.^2);
expT = T^2 / (a*b*n);

% sums of squares
ssA = expA - expT;
ssB = expB - expT;
ssAB = expAB - expA - expB + expT;
ssS = expS - expT;
ssAS = expAS - expA - expS + expT;
ssBS = expBS - expB - expS + expT;
ssABS = expY - expAB - expAS - expBS + expA + expB + expS - expT;
ssTot = expY - expT;

% mean squares
msA = ssA / dfA;
msB = ssB / dfB;
msAB = ssAB / dfAB;
msS = ssS / dfS;
msAS = ssAS / dfAS;
msBS = ssBS / dfBS;
msABS = ssABS / dfABS;

% f statistic
fA = msA / msAS;
fB = msB / msBS;
fAB = msAB / msABS;

% p values
pA = 1-fcdf(fA,dfA,dfAS);
pB = 1-fcdf(fB,dfB,dfBS);
pAB = 1-fcdf(fAB,dfAB,dfABS);

% return values
stats = {'Source','SS','df','MS','F','p';...
         FACTNAMES{1}, ssA, dfA, msA, fA, pA;...
         FACTNAMES{2}, ssB, dfB, msB, fB, pB;...
         [FACTNAMES{1} ' x ' FACTNAMES{2}], ssAB, dfAB, msAB, fAB, pAB;...
         [FACTNAMES{1} ' x Subj'], ssAS, dfAS, msAS, [], [];...
         [FACTNAMES{1} ' x Subj'], ssBS, dfBS, msBS, [], [];...
         [FACTNAMES{1} ' x ' FACTNAMES{2} ' x Subj'], ssABS, dfABS, msABS, [], []};

 return

author: PatternRecognition

repo: OpenBMI


text:
function [mg,nu,sig,info] = spm_rice_mixture(h,x,K)
% Fit a mixture of Ricians to a histogram
% FORMAT [mg,nu,sig] = spm_rice_mixture(h,x,K)
% h    - histogram counts
% x    - bin positions (plot(x,h) to see the histogram)
% K    - number of Ricians
% mg   - integral under each Rician
% nu   - "mean" parameter of each Rician
% sig  - "standard deviation" parameter of each Rician
% info - This struct can be used for plotting the fit as:
%            plot(info.x(:),info.p,'--',info.x(:), ...
%                 info.h/sum(info.h)/info.md,'b.', ...
%                 info.x(:),info.sp,'r');
%
% An EM algorithm is used, which involves alternating between computing
% belonging probabilities, and then the parameters of the Ricians.
% The Koay inversion technique is used to compute the Rician parameters
% from the sample means and standard deviations. This is described at
% https://en.wikipedia.org/wiki/Rician_distribution
%__________________________________________________________________________
% Copyright (C) 2012-2019 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_rice_mixture.m 7595 2019-05-23 13:48:53Z mikael $

mg  = ones(K,1)/K;
nu  = (0:(K-1))'*max(x)/(K+1);
sig = ones(K,1)*max(x)/K/10;
lam = (sum(x.*h)/sum(h)/K).^2;

m0 = zeros(K,1);
m1 = zeros(K,1);
m2 = zeros(K,1);
ll = -Inf;
for iter=1:10000
    p  = zeros(numel(x),K);
    for k=1:K
        % Product Rule
        % p(class=k, x | mg, nu, sig) = p(class=k|mg) p(x | nu, sig, class=k)
        p(:,k) = mg(k)*ricepdf(x(:),nu(k),sig(k)^2) + eps;
    end

    % Sum Rule
    % p(x | mg, nu, sig) = \sum_k p(class=k, x | mg, nu, sig)
    sp  = sum(p,2);
    oll = ll;
    ll  = sum(log(sp).*h(:)); % Log-likelihood
    if ll-oll<1e-8*sum(h), break; end

    %fprintf('%g\n',ll);
    %md = mean(diff(x));
    %plot(x(:),p,'--',x(:),h/sum(h)/md,'b.',x(:),sp,'r'); drawnow

    % Bayes Rule
    % p(class=k | x, mg, nu, sig) = p(class=k, x | mg, nu, sig) / p(x | mg, nu, sig)
    p = bsxfun(@rdivide,p,sp);

    % Compute moments from the histograms, weighted by the responsibilities (p).
    for k=1:K
        m0(k) = sum(p(:,k).*h(:));              % Number of voxels in class k
        m1(k) = sum(p(:,k).*h(:).*x(:));        % Sum of the intensities in class k
        m2(k) = sum(p(:,k).*h(:).*x(:).*x(:));  % Sum of squares of intensities in class k
    end

    mg = m0/sum(m0); % Mixing proportions
    for k=1:K
        mu1 = m1(k)./m0(k);                                % Mean
        mu2 = (m2(k)-m1(k)*m1(k)/m0(k)+lam*1e-3)/(m0(k)+1e-3); % Variance

        % Compute nu & sig from mean and variance
        [nu(k),sig(k)] = moments2param(mu1,mu2);
    end
    %disp([nu'; sig'])
end

if nargout >= 4
    % This info can be used for plotting the fit
    info    = struct;
    info.x  = x;
    info.h  = h;
    info.p  = p;
    info.sp = sp;
    info.md = mean(diff(x));
end
%__________________________________________________________________________

%__________________________________________________________________________
function [nu,sig] = moments2param(mu1,mu2)
% Rician parameter estimation (nu & sig) from mean (mu1) and variance
% (mu2) via the Koay inversion technique.
% This follows the scheme at
% https://en.wikipedia.org/wiki/Rice_distribution#Parameter_estimation_.28the_Koay_inversion_technique.29
% This Wikipedia description is based on:
% Koay, C.G. and Basser, P. J., Analytically exact correction scheme
% for signal extraction from noisy magnitude MR signals,
% Journal of Magnetic Resonance, Volume 179, Issue = 2, p. 317–322, (2006)

r     = mu1/sqrt(mu2);
theta = sqrt(pi/(4-pi));
if r>theta
    for i=1:256
        xi    = 2+theta^2-pi/8*exp(-theta^2/2)*((2+theta^2)*besseli(0,theta^2/4)+theta^2*besseli(1,theta^2/4))^2;
        g     = sqrt(xi*(1+r^2)-2);
        if abs(theta-g)<1e-6, break; end
        theta = g;
    end
    if ~isfinite(xi), xi = 1; end
    sig = sqrt(mu2)/sqrt(xi);
    nu  = sqrt(mu1^2+(xi-2)*sig^2);
else
    nu  = 0;
    sig = (2^(1/2)*(mu1^2 + mu2)^(1/2))/2;
end
%__________________________________________________________________________

%__________________________________________________________________________
function p = ricepdf(x,nu,sig2)
% Rician PDF
% p = ricepdf(x,nu,sig2)
% https://en.wikipedia.org/wiki/Rice_distribution#Characterization
p       = zeros(size(x));
tmp     = -(x.^2+nu.^2)./(2*sig2);
msk     = (tmp > -95) & (x*(nu/sig2) < 85) ; % Identify where Rice probability can be computed
p(msk)  = (x(msk)./sig2).*exp(tmp(msk)).*besseli(0,x(msk)*(nu/sig2)); % Use Rician distribution
p(~msk) = (1./sqrt(2*pi*sig2))*exp((-0.5/sig2)*(x(~msk)-nu).^2);      % Use Gaussian distribution

author: spm


text:
function [V, policy, iter, cpu_time] = mdp_policy_iteration(P, R, discount, policy0, max_iter, eval_type)

% max_iter=2000;
% eval_type=1;

% mdp_policy_iteration  Resolution of discounted MDP
%                       with policy iteration algorithm
% Arguments ---------------------------------------------------------------
% Let S = number of states, A = number of actions
%   P(SxSxA) = transition matrix
%              P could be an array with 3 dimensions or
%              a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
%   R(SxSxA) or (SxA) = reward matrix
%              R could be an array with 3 dimensions (SxSxA) or
%              a cell array (1xA), each cell containing a sparse matrix (SxS) or
%              a 2D array(SxA) possibly sparse
%   discount = discount rate, in ]0, 1[
%   policy0(S) = starting policy, optional
%   max_iter = maximum number of iteration to be done, upper than 0,
%              optional (default 1000)
%   eval_type = type of function used to evaluate policy:
%              0 for mdp_eval_policy_matrix, else mdp_eval_policy_iterative
%              optional (default 0)
% Evaluation --------------------------------------------------------------
%   V(S)   = value function
%   policy(S) = optimal policy
%   iter     = number of done iterations
%   cpu_time = used CPU time
%--------------------------------------------------------------------------
% In verbose mode, at each iteration, displays the number
% of differents actions between policy n-1 and n

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

cpu_time = cputime;

global mdp_VERBOSE;

% check of arguments
if iscell(P); S = size(P{1},1); else S = size(P,1); end;
if iscell(P); A = length(P); else A = size(P,3); end

if discount <= 0 || discount >= 1
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: Discount rate must be in ]0; 1[')
    disp('--------------------------------------------------------')
elseif nargin > 3 && (size(policy0,1)~=S || any(mod(policy0,1)) || any(policy0<1)|| any(policy0>A) )
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: policy0 must a (Sx1) vector with integer from 1 to A')
    disp('--------------------------------------------------------')
elseif nargin > 4 && max_iter <= 0
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: The maximum number of iteration must be upper than 0')
    disp('--------------------------------------------------------')
else

    PR = mdp_computePR(P,R);

    % initialization of optional arguments
    if nargin < 6; eval_type = 0; end;
    if nargin < 5; max_iter = 1000; end;
    if nargin < 4;
        % initialization of policy:
        % the one wich maximizes the expected immediate reward
        [nil, policy0] = mdp_bellman_operator(P,PR,discount,zeros(S,1));
    end;

%     if mdp_VERBOSE; disp('  Iteration  Number_of_different_actions'); end;
        disp('  Iteration  Number_of_different_actions');

    iter = 0;
    policy = policy0;
    is_done = false;
    while ~is_done
        iter = iter + 1;
        if  (eval_type==0)
            V = mdp_eval_policy_matrix(P,PR,discount,policy);
        else
            V = mdp_eval_policy_iterative(P,PR,discount,policy);
        end;
        [nil, policy_next] = mdp_bellman_operator(P,PR,discount,V);

	n_different = sum(policy_next ~= policy);
%         if mdp_VERBOSE; disp(['       ' num2str(iter) '                 '  num2str(n_different)]); end;
        disp(['       ' num2str(iter) '                 '  num2str(n_different)]);

        if all(policy_next==policy) || iter == max_iter || (iter > 20 && n_different <=5)
            is_done = true;
        else
            policy = policy_next;
        end;
    end;

end;

cpu_time = cputime - cpu_time;

author: matthieukomorowski


text:
function sub = ind2subv(siz,index)
%IND2SUBV   Subscript vector from linear index.
% IND2SUBV(SIZ,IND) returns a vector of the equivalent subscript values
% corresponding to a single index into an array of size SIZ.
% If IND is a vector, then the result is a matrix, with subscript vectors
% as rows.

% Written by Tom Minka
% Part of Tom Minka's lightspeed package.
% (c) Microsoft Corporation. All rights reserved.

n = length(siz);
cum_size = cumprod(siz(:)');
prev_cum_size = [1 cum_size(1:end-1)];
index = index(:) - 1;
sub = rem(repmat(index,1,n),repmat(cum_size,length(index),1));
sub = floor(sub ./ repmat(prev_cum_size,length(index),1))+1;

% slow way
%for dim = n:-1:1
%  sub(:,dim) = floor(index/cum_size(dim))+1;
%  index = rem(index,cum_size(dim));
%end

author: andrewssobral

repo: lrslibrary


text:
function [LS,LD,I] = isolines(V,F,S,iso)
  % ISOLINES compute a list of isolines for a scalar field S defined on the
  % mesh (V,F)
  %
  % [LS,LD,I] = isolines(V,F,S,iso)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of triangle indices
  %   S  #V list of scalar values defined on V
  %   iso #iso list of iso values
  % Outputs:
  %   LS  #L by dim list of isoline start positions
  %   LD  #L by dim list of isoline destination positions
  %   I  #L list of indices into iso revealing corresponding iso value
  %
  % alternative: tracing isolines should result in smaller plot data (every
  % vertex only appears once
  %
  % Example:
  %   iso = linspace(min(S),max(S),10);
  %   [LS,LD] = isolines(V,F,S,iso);
  %   colormap(jet(numel(iso)-1));
  %   tsurf(F,V,'CData',S,fphong);
  %   hold on;
  %   plot([LS(:,1) LD(:,1)]',[LS(:,2) LD(:,2)]', ...
  %     'Color','k','LineWidth',10);
  %   hold off;
  %

  % make sure iso is a ROW vector
  iso = iso(:)';
  % number of isolines
  niso = numel(iso);

  % number of domain positions
  n = size(V,1);
  % number of dimensions
  dim = size(V,2);

  % number of faces
  m = size(F,1);

  % Rename for convenience
  S1 = S(F(:,1),:);
  S2 = S(F(:,2),:);
  S3 = S(F(:,3),:);

  % t12(i,j) reveals parameter t where isovalue j falls on the line from
  % corner 1 to corner 2 of face i
  t12 = bsxfun(@rdivide,bsxfun(@minus,iso,S1),S2-S1);
  t23 = bsxfun(@rdivide,bsxfun(@minus,iso,S2),S3-S2);
  t31 = bsxfun(@rdivide,bsxfun(@minus,iso,S3),S1-S3);

  % replace values outside [0,1] with NaNs
  t12( (t12<-eps)|(t12>(1+eps)) ) = NaN;
  t23( (t23<-eps)|(t23>(1+eps)) ) = NaN;
  t31( (t31<-eps)|(t31>(1+eps)) ) = NaN;

  % masks for line "parallel" to 12
  l12 = ~isnan(t23) & ~isnan(t31);
  l23 = ~isnan(t31) & ~isnan(t12);
  l31 = ~isnan(t12) & ~isnan(t23);

  % find non-zeros (lines) from t23 to t31
  [F12,I12,~] = find(l12);
  [F23,I23,~] = find(l23);
  [F31,I31,~] = find(l31);
  % indices directly into t23 and t31 corresponding to F12 and I12
  %ti12 = sub2ind(size(l12),F12,I12);
  %ti23 = sub2ind(size(l23),F23,I23);
  %ti31 = sub2ind(size(l31),F31,I31);
  % faster sub2ind
  ti12 = F12+(I12-1)*size(l12,1);
  ti23 = F23+(I23-1)*size(l23,1);
  ti31 = F31+(I31-1)*size(l31,1);

  % compute actual position values
  LS = [ ...
    ... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t23(ti12), V(F(F12,2),:)) + ...
    bsxfun(@times,  t23(ti12), V(F(F12,3),:)) ...
    ; ... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t31(ti23), V(F(F23,3),:)) + ...
    bsxfun(@times,  t31(ti23), V(F(F23,1),:)) ...
    ;... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t12(ti31), V(F(F31,1),:)) + ...
    bsxfun(@times,  t12(ti31), V(F(F31,2),:)) ...
    ];
  LD = [...
    ... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t31(ti12), V(F(F12,3),:)) + ...
    bsxfun(@times,  t31(ti12), V(F(F12,1),:)) ...
    ;... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t12(ti23), V(F(F23,1),:)) + ...
    bsxfun(@times,  t12(ti23), V(F(F23,2),:)) ...
    ;... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t23(ti31), V(F(F31,2),:)) + ...
    bsxfun(@times,  t23(ti31), V(F(F31,3),:)) ...
    ];

  % I is just concatenation of each I12
  I = [I12;I23;I31];

end

author: alecjacobson

repo: gptoolbox


text:
function a = r8row_sort_heap_a ( m, n, a )

%*****************************************************************************80
%
%% R8ROW_SORT_HEAP_A ascending heapsorts an R8ROW.
%
%  Discussion:
%
%    In lexicographic order, the statement "X < Y", applied to two real
%    vectors X and Y of length M, means that there is some index I, with
%    1 <= I <= M, with the property that
%
%      X(J) = Y(J) for J < I,
%    and
%      X(I) < Y(I).
%
%    In other words, the first time they differ, X is smaller.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), the array of M rows of N-vectors.
%
%    Output, real A(M,N), the rows of A have been sorted.
%
  if ( m <= 0 )
    return
  end

  if ( n <= 1 )
    return
  end
%
%  Initialize.
%
  i = 0;
  indx = 0;
  isgn = 0;
  j = 0;
%
%  Call the external heap sorter.
%
  while ( 1 )

    [ indx, i, j ] = sort_heap_external ( m, indx, isgn );
%
%  Interchange the I and J objects.
%
    if ( 0 < indx )

      a = r8row_swap ( m, n, a, i, j );
%
%  Compare the I and J objects.
%
    elseif ( indx < 0 )

      isgn = r8row_compare ( m, n, a, i, j );

    elseif ( indx == 0 )

      break

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [nodes, edges, faces] = boundedVoronoi2d(box, germs)
%BOUNDEDVORONOI2D Return a bounded voronoi diagram as a graph structure
%
%   [NODES, EDGES, FACES] = boundedVoronoi2d(BOX, GERMS)
%   GERMS an array of points with dimension 2
%   NODES, EDGES, FACES: usual graph representation, FACES as cell array
%
%   Example
%   [n, e, f] = boundedVoronoi2d([0 100 0 100], rand(100, 2)*100);
%   drawGraph(n, e);
%
%   See also
%     graphs, boundedCentroidalVoronoi2d, clipGraph, clipGraphPolygon

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2007-01-12
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% uniformize input for box.
box = box';
box = box(:);

% add points far enough
boxSizeX  = box(2) - box(1);
boxSizeY  = box(4) - box(3);
farPoints = [...
    box(2)+2*boxSizeX  box(4)+3*boxSizeY;...
    box(1)-3*boxSizeX  box(4)+2*boxSizeY;...
    box(1)-2*boxSizeX  box(3)-3*boxSizeY;...
    box(2)+3*boxSizeX  box(3)-2*boxSizeY;...
    ];

% extract voronoi vertices and face structure
[V, C] = voronoin([germs ; farPoints]);

% initialize graph structure, without edges and without faces
nodes = V(2:end, :);
edges = zeros(0, 2);
faces = {};

for i = 1:size(germs, 1)
    cell = C{i};
    cell = cell-1;
    edges = [edges; sort([cell' cell([2:end 1])'], 2)]; %#ok<AGROW>
    faces{length(faces)+1} = cell; %#ok<AGROW>
end

edges = unique(edges, 'rows');

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function arg = calc_acosInput(a, b)
% calculate normalized dot product and enforce result to be within limits
% [-1, 1] to avoid acos function throwing an error caused by numerical
% inaccuracy

    expr = dot(a,b)/(norm(a)*norm(b));

    % fix to avoid error due to numerical inaccuracy in combination with
    % acos calculation
    arg = max(min(expr, 1.0), -1.0);

end

author: TUMFTM


text:
function [coordVERTICES,varargout] = READ_stl(stlFILENAME,varargin)
% READ_stlascii  Read mesh data in the form of an <*.stl> file
%==========================================================================
% FILENAME:          READ_stl.m
% AUTHOR:            Adam H. Aitkenhead
% INSTITUTION:       The Christie NHS Foundation Trust
% CONTACT:           adam.aitkenhead@physics.cr.man.ac.uk
% DATE:              29th March 2010
% PURPOSE:           Read mesh data in the form of an <*.stl> file.
%
% USAGE:
%
%     [coordVERTICES,coordNORMALS,stlNAME] = READ_stl(stlFILENAME,stlFORMAT)
%
% INPUT PARAMETERS:
%
%     stlFILENAME   - String - Mandatory - The filename of the STL file.
%
%     stlFORMAT     - String - Optional  - The format of the STL file:
%                                        'ascii' or 'binary'
%
% OUTPUT PARAMETERS:
%
%     coordVERTICES - Nx3x3 array - Mandatory
%                                 - An array defining the vertex positions
%                                   for each of the N facets, with:
%                                   1 row for each facet
%                                   3 cols for the x,y,z coordinates
%                                   3 pages for the three vertices
%
%     coordNORMALS  - Nx3 array   - Optional
%                                 - An array defining the normal vector for
%                                   each of the N facets, with:
%                                   1 row for each facet
%                                   3 cols for the x,y,z components of the vector
%
%     stlNAME       - String      - Optional  - The name of the STL object.
%
%==========================================================================

%==========================================================================
% VERSION:  USER:  CHANGES:
% -------   -----  --------
% 100329    AHA    Original version
% 100513    AHA    Totally reworked the code.  Now use textscan to read the
%                  file all at once, rather than one line at a time with
%                  fgetl.  Major speed improvment.
% 100623    AHA    Combined code which reads ascii STLS and code which
%                  reads binary STLs into a single function.
%==========================================================================

%==========================================================================
% STL ascii file format
%======================
% ASCII STL files have the following structure.  Technically each facet
% could be any 2D shape, but in practice only triangular facets tend to be
% used.  The present code ONLY works for meshes composed of triangular
% facets.
%
% solid object_name
% facet normal x y z
%   outer loop
%     vertex x y z
%     vertex x y z
%     vertex x y z
%   endloop
% endfacet
%
% <Repeat for all facets...>
%
% endsolid object_name
%==========================================================================

%==========================================================================
% STL binary file format
%=======================
% Binary STL files have an 84 byte header followed by 50-byte records, each
% describing a single facet of the mesh.  Technically each facet could be
% any 2D shape, but that would screw up the 50-byte-per-facet structure, so
% in practice only triangular facets are used.  The present code ONLY works
% for meshes composed of triangular facets.
%
% HEADER:
% 80 bytes:  Header text
% 4 bytes:   (int) The number of facets in the STL mesh
%
% DATA:
% 4 bytes:  (float) normal x
% 4 bytes:  (float) normal y
% 4 bytes:  (float) normal z
% 4 bytes:  (float) vertex1 x
% 4 bytes:  (float) vertex1 y
% 4 bytes:  (float) vertex1 z
% 4 bytes:  (float) vertex2 x
% 4 bytes:  (float) vertex2 y
% 4 bytes:  (float) vertex2 z
% 4 bytes:  (float) vertex3 x
% 4 bytes:  (float) vertex3 y
% 4 bytes:  (float) vertex3 z
% 2 bytes:  Padding to make the data for each facet 50-bytes in length
%   ...and repeat for next facet...
%==========================================================================

if nargin==2
  stlFORMAT = lower(varargin{1});
else
  stlFORMAT = 'auto';
end

%If necessary, identify whether the STL is ascii or binary:
if strcmp(stlFORMAT,'ascii')==0 && strcmp(stlFORMAT,'binary')==0
  stlFORMAT = IDENTIFY_stl_format(stlFILENAME);
end

%Load the STL file:
if strcmp(stlFORMAT,'ascii')
  [coordVERTICES,coordNORMALS,stlNAME] = READ_stlascii(stlFILENAME);
elseif strcmp(stlFORMAT,'binary')
  [coordVERTICES,coordNORMALS] = READ_stlbinary(stlFILENAME);
  stlNAME = 'unnamed_object';
end %if

%Prepare the output arguments
if nargout == 2
  varargout(1) = {coordNORMALS};
elseif nargout == 3
  varargout(1) = {coordNORMALS};
  varargout(2) = {stlNAME};
end

end %function

%==========================================================================
function [stlFORMAT] = IDENTIFY_stl_format(stlFILENAME)
% IDENTIFY_stl_format  Test whether an stl file is ascii or binary

% Open the file:
fidIN = fopen(stlFILENAME);

% Check the file size first, since binary files MUST have a size of 84+(50*n)
fseek(fidIN,0,1);         % Go to the end of the file
fidSIZE = ftell(fidIN);   % Check the size of the file

if rem(fidSIZE-84,50) > 0

  stlFORMAT = 'ascii';

else

  % Files with a size of 84+(50*n), might be either ascii or binary...

  % Read first 80 characters of the file.
  % For an ASCII file, the data should begin immediately (give or take a few
  % blank lines or spaces) and the first word must be 'solid'.
  % For a binary file, the first 80 characters contains the header.
  % It is bad practice to begin the header of a binary file with the word
  % 'solid', so it can be used to identify whether the file is ASCII or
  % binary.
  fseek(fidIN,0,-1);        % Go to the start of the file
  firsteighty = char(fread(fidIN,80,'uchar')');

  % Trim leading and trailing spaces:
  firsteighty = strtrim(firsteighty);

  % Take the first five remaining characters, and check if these are 'solid':
  firstfive = firsteighty(1:min(5,length(firsteighty)));

  % Double check by reading the last 80 characters of the file.
  % For an ASCII file, the data should end (give or take a few
  % blank lines or spaces) with 'endsolid <object_name>'.
  % If the last 80 characters contains the word 'endsolid' then this
  % confirms that the file is indeed ASCII.
  if strcmp(firstfive,'solid')

    fseek(fidIN,-80,1);     % Go to the end of the file minus 80 characters
    lasteighty = char(fread(fidIN,80,'uchar')');

    if findstr(lasteighty,'endsolid')
      stlFORMAT = 'ascii';
    else
      stlFORMAT = 'binary';
    end

  else
    stlFORMAT = 'binary';
  end

end

% Close the file
fclose(fidIN);

end %function
%==========================================================================

%==========================================================================
function [coordVERTICES,coordNORMALS,stlNAME] = READ_stlascii(stlFILENAME)
% READ_stlascii  Read mesh data in the form of an ascii <*.stl> file

% Read the ascii STL file
fidIN = fopen(stlFILENAME);
fidCONTENTcell = textscan(fidIN,'%s','delimiter','\n');                  %Read all the file
fidCONTENT = fidCONTENTcell{:}(logical(~strcmp(fidCONTENTcell{:},'')));  %Remove all blank lines
fclose(fidIN);

% Read the STL name
if nargout == 3
  line1 = char(fidCONTENT(1));
  if (size(line1,2) >= 7)
    stlNAME = line1(7:end);
  else
    stlNAME = 'unnamed_object';
  end
end

% Read the vector normals
if nargout >= 2
  stringNORMALS = char(fidCONTENT(logical(strncmp(fidCONTENT,'facet normal',12))));
  coordNORMALS  = str2num(stringNORMALS(:,13:end));
end

% Read the vertex coordinates
facetTOTAL       = sum(strcmp(fidCONTENT,'endfacet'));
stringVERTICES   = char(fidCONTENT(logical(strncmp(fidCONTENT,'vertex',6))));
coordVERTICESall = str2num(stringVERTICES(:,7:end));
cotemp           = zeros(3,facetTOTAL,3);
cotemp(:)        = coordVERTICESall;
coordVERTICES    = shiftdim(cotemp,1);

end %function
%==========================================================================

%==========================================================================
function [coordVERTICES,coordNORMALS] = READ_stlbinary(stlFILENAME)
% READ_stlbinary  Read mesh data in the form of an binary <*.stl> file

% Open the binary STL file
fidIN = fopen(stlFILENAME);

% Read the header
fseek(fidIN,80,-1);                   % Move to the last 4 bytes of the header
facetcount = fread(fidIN,1,'int32');  % Read the number of facets in the stl file

% Initialise arrays into which the STL data will be loaded:
coordNORMALS  = zeros(facetcount,3);
coordVERTICES = zeros(facetcount,3,3);

% Read the data for each facet:
for loopF = 1:facetcount

  tempIN = fread(fidIN,3*4,'float');

  coordNORMALS(loopF,1:3)    = tempIN(1:3);    % x,y,z components of the facet's normal vector
  coordVERTICES(loopF,1:3,1) = tempIN(4:6);    % x,y,z coordinates of vertex 1
  coordVERTICES(loopF,1:3,2) = tempIN(7:9);    % x,y,z coordinates of vertex 2
  coordVERTICES(loopF,1:3,3) = tempIN(10:12);  % x,y,z coordinates of vertex 3

  fread(fidIN,1,'int16');   % Move to the start of the next facet.  Using fread is much quicker than using fseek(fidIN,2,0);

end %for

% Close the binary STL file
fclose(fidIN);

end %function
%==========================================================================

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function data = outlier(N, r, dist, outliers, noise)

    if nargin < 1
        N = 600;
    end
    if nargin < 2
        r = 20;
    end
    if nargin < 3
        dist = 30;
    end
    if nargin < 4
        outliers = 0.04;
    end
    if nargin < 5
        noise = 5;
    end

    N1 = round(N * (.5-outliers));
    N2 = N1;
    N3 = round(N * outliers);
    N4 = N-N1-N2-N3;

    phi1 = rand(N1,1) * pi;
    r1 = sqrt(rand(N1,1))*r;
    P1 = [-dist + r1.*sin(phi1) r1.*cos(phi1) zeros(N1,1)];

    phi2 = rand(N2,1) * pi;
    r2 = sqrt(rand(N2,1))*r;
    P2 = [dist - r2.*sin(phi2) r2.*cos(phi2) 3*ones(N2,1)];

    P3 = [rand(N3,1)*noise dist+rand(N3,1)*noise 2*ones(N3,1)];

    P4 = [rand(N4,1)*noise -dist+rand(N4,1)*noise ones(N4,1)];

    data = [P1; P2; P3; P4];

end

repo: gspbox


text:
% Filename: givepd.m
% Author: Shashank G. Sawant
% email: sgsawant@gmail.com
% Description: Given 2 sinusoidal signals of the
% same frequency, the function gives the "phase difference" between the
% 2 given signals
% The phase difference is in RADIANS!!!!!!
% The output is limited to pd={-pi,pi}radians
% Time Stamp: 2010 October 19 2043hrs

% pd - output phase difference (in radians)
% v - first sinusoidal signal
% i - second sinusoidal signal
% Note: v and i should have the same frequency
function pd=givepd(v,i)
L=length(v);
if(L~=length(i))
    error('The length of the 2 sinusoidal input vectors is not same!!!');
end

% The following block calculates the FFT
NFFT = 2^nextpow2(L);
V = fft(v,NFFT)/L;                  %Fourier Transform of signal v

% The following block calculates the phase of the most significant
% frequency component
[value,index]=max(2*abs(V(1:NFFT/2+1)));
pV = angle(V(index));

% The following block calculates the FFT
NFFT = 2^nextpow2(L);
I = fft(i,NFFT)/L;                  %Fourier Transform of signal i

% The following block calculates the phase of the most significant
% frequency component
[value,index]=max(2*abs(I(1:NFFT/2+1)));
pI = angle(I(index));

% The following is the phase difference between the 2 signals
pd=pV-pI;

% The code below limits the output to pd={-pi,pi}radians
while 1
    if pd>pi
        pd=pd-2*pi;
    elseif pd<-pi
        pd=pd+2*pi;
    else
        break;
    end
end

return;
end


text:
function letter=readLetter(snap)
%READLETTER reads the character fromthe character's binary image.
%   LETTER=READLETTER(SNAP) outputs the character in class 'char' from the
%   input binary image SNAP.

load NewTemplates % Loads the templates of characters in the memory.
snap=imresize(snap,[42 24]); % Resize the input image so it can be compared with the template's images.
comp=[ ];
for n=1:length(NewTemplates)
    sem=corr2(NewTemplates{1,n},snap); % Correlation the input image with every image in the template for best matching.
    comp=[comp sem]; % Record the value of correlation for each template's character.
end
vd=find(comp==max(comp)); % Find the index which correspond to the highest matched character.
%*-*-*-*-*-*-*-*-*-*-*-*-*-
% Accodrding to the index assign to 'letter'.
% Alphabets listings.
if vd==1 || vd==2
    letter='A';
elseif vd==3 || vd==4
    letter='B';
elseif vd==5
    letter='C';
elseif vd==6 || vd==7
    letter='D';
elseif vd==8
    letter='E';
elseif vd==9
    letter='F';
elseif vd==10
    letter='G';
elseif vd==11
    letter='H';
elseif vd==12
    letter='I';
elseif vd==13
    letter='J';
elseif vd==14
    letter='K';
elseif vd==15
    letter='L';
elseif vd==16
    letter='M';
elseif vd==17
    letter='N';
elseif vd==18 || vd==19
    letter='O';
elseif vd==20 || vd==21
    letter='P';
elseif vd==22 || vd==23
    letter='Q';
elseif vd==24 || vd==25
    letter='R';
elseif vd==26
    letter='S';
elseif vd==27
    letter='T';
elseif vd==28
    letter='U';
elseif vd==29
    letter='V';
elseif vd==30
    letter='W';
elseif vd==31
    letter='X';
elseif vd==32
    letter='Y';
elseif vd==33
    letter='Z';
    %*-*-*-*-*
% Numerals listings.
elseif vd==34
    letter='1';
elseif vd==35
    letter='2';
elseif vd==36
    letter='3';
elseif vd==37 || vd==38
    letter='4';
elseif vd==39
    letter='5';
elseif vd==40 || vd==41 || vd==42
    letter='6';
elseif vd==43
    letter='7';
elseif vd==44 || vd==45
    letter='8';
elseif vd==46 || vd==47 || vd==48
    letter='9';
else
    letter='0';
end
end


text:
function value = SPEED_OF_LIGHT( )
%SPEED_OF_LIGHT Constant for speed of light in meters per second
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

value = 299792458;

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [handle]  = ml_plot_rvm_boundary( data, labels, model, varargin)
% ML_PLOT_RVM_BOUNDARY plots the training data
%   and decision boundary, given a model produced by LIBSVM
%   input ----------------------------------------------------------------
%
%       o data      : (1 x 1), number of data points to generate.
%
%       o labels    : (1 x 1), dimension of the data. [2 or 3]
%
%       o model     : (1 x 1), number of classes.
%
%       o options   : struct
%
%       o varargin  : string, if 'draw' draw contours otherwise, don't
%
%%  Plot function
handle = figure('Color', [1 1 1]);
labels(find(labels==-1)) = 0;

% Parse model params
weights = model.weights;
kernel_ = model.kernel_;
width   = model.width;
bias    = model.bias;
RVs     = model.RVs;
acc     = model.predict_accuracy;
N_RVs   = model.N_RVs;

% Make classification predictions over a grid of values
xplot = linspace(min(data(:,1)), max(data(:,1)), 100)';
yplot = linspace(min(data(:,2)), max(data(:,2)), 100)';
[grid1, grid2] = meshgrid(xplot, yplot);
Xgrid = [grid1(:),grid2(:)];

% Evaluate RVM
PHI		= SB1_KernelFunction(Xgrid,RVs,kernel_,width);
y_grid  = PHI*weights + bias;

% apply sigmoid for probabilities (option)
p_grid	= 1./(1+exp(-y_grid));

% Plot the RVM Contours
colormap hot
if (size(varargin, 2) == 1) && strcmp(varargin{1}, 'draw')
    contourf(grid1,grid2, reshape(p_grid,size(grid1)), 50, 'LineStyle', 'none');
    legend_names = {'Decision values (p(x))', 'Class 1','Class 0', 'p=0.5', 'p=0.75', 'p=0.25'};
    colorbar
elseif (size(varargin, 2) == 1) && strcmp(varargin{1}, 'surf')
    surf(grid1,grid2, reshape(p_grid,size(grid1))); shading interp;
    legend_names = {'Decision Values f(x)','f(x)=0', 'f(x)= +1','f(x)=-1','Class 1','Class -1'};
    colorbar
else
    legend_names = {'Class 1','Class 0', 'p=0.5', 'p=0.75', 'p=0.25'};
end
hold on;

% Plot the training data on top of the boundary
pos = find(labels == 1);
neg = find(labels == 0);
plot(data(pos,1), data(pos,2), 'ko', 'MarkerFaceColor', 'r'); hold on;
plot(data(neg,1), data(neg,2), 'ko', 'MarkerFaceColor', 'g'); hold on;

% Show decision boundaries p=0.5 (y=0) and illustrate p=0.75(y=1) and p=0.25(y=-1)
contour(grid1,grid2,reshape(y_grid,size(grid1)),[0.5 0.5],'-', 'LineWidth', 2,  'Color', 'k');
contour(grid1,grid2,reshape(y_grid,size(grid1)),[0.75 0.75],'-', 'LineWidth',0.5, 'Color', 'r');
contour(grid1,grid2,reshape(y_grid,size(grid1)),[0.25 0.25],'-', 'LineWidth',0.5, 'Color', 'g');

% Show relevance vectors
scatter(RVs(:,1),RVs(:,2),70,'o','MarkerEdgeColor', [1 1 1], 'LineWidth', 2);
legend(legend_names,'Location','NorthWest')
hold off
title(sprintf('RVM \\sigma = %g, RV = %d, Acc = %g', width, N_RVs, acc), 'Interpreter','tex','FontSize', 20, 'FontWeight', 'Normal');

axis tight
grid on
box on

end

author: epfl-lasa


text:
function ar=lpcrand(p,n,bw)
% generate n random stable polynomials of order p with a minimum pole
% bandwidth of bw*fs where fs is the sampling fequency.
% To limit the pole radius to r set bw=-log(r)/pi
% bw may be a vector specifying a different max bandwidth for each row

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

if nargin < 3
   bw=0;
   if nargin < 2
      n=1;
   end
end
if p
   if ~bw
      ar=lpcrf2ar(2*rand(n,p+1)-1);
   else
      k=exp(-pi*bw(:)*(0:p));
      if size(k,1)==1
         ar=lpcrf2ar(2*rand(n,p+1)-1).*k(ones(n,1),:);
      else
         ar=lpcrf2ar(2*rand(n,p+1)-1).*k;
      end
   end
else
   ar=ones(n,1);
end

author: decouples


text:
function data = readdata(filename, m)
%-------------------------------------------------------------------------------
% data = readdata(filename, m)
% read binary data with m double
% by Dr. Yudan Yi
%-------------------------------------------------------------------------------
data = [];
if (nargin<2) return; end;
if (isempty(m)|m<=0) return; end;
fdat = fopen(filename, 'rb');
if (fdat==-1) return; end;
data = fread(fdat, 'double');
n = size(data,1);
nn = floor(n/m);
data = reshape(data(1:nn*m), m, nn)';
data = sortrows(data, 1);
loc = find(diff(data(:,1))==0);
data(loc, :) = [];
fclose(fdat);
%-------------------------------------------------------------------------------

repo: ekfmonoslam


text:
% eeg_options() - eeglab option script
%
% Note: DO NOT EDIT, instead use pop_editoptions() or the menu
%       /File/Maximize memory in EEGLAB gui

% STUDY and file options (set these checkboxes if you intend to work with studies)
option_storedisk = 0 ; % If set, keep at most one dataset in memory. This allows processing hundreds of datasets within studies.
option_savetwofiles = 1 ; % If set, save not one but two files for each dataset (header and data). This allows faster data loading in studies.
option_saveversion6 = 1 ; % If set, write Matlab files in Matlab v6.5 (max compatibility). If not, write files in Matlab v7.3 (larger than 2Gb).
% Memory options
option_single = 0 ; % If set, use single precision under Matlab 7.x. This saves RAM but can lead to rare numerical imprecisions.
option_memmapdata = 0 ; % If set, use memory mapped array under Matlab 7.x. This may slow down some computation (beta).
option_eegobject = 0 ; % If set, use the EEGLAB EEG object instead of the standard EEG structure (beta).
% ICA options
option_computeica = 0 ; % If set, precompute ICA activations. This requires more RAM but allows faster plotting of component activations.
option_scaleicarms = 1 ; % If set, scale ICA component activities to RMS (Root Mean Square) in microvolt (recommended).
% Folder options
option_rememberfolder = 1 ; % If set, when browsing to open a new dataset assume the folder/directory of previous dataset.

% Toolbox options
option_donotusetoolboxes = 0 ; % If set, do not use Matlab additional toolboxes functions even if they are present (need to restart EEGLAB).
% EEGLAB connectivity and support
option_checkversion = 1 ; % If set, check for new version of EEGLAB at startup.
option_chat = 0 ; % If set, enable EEGLAB chat - currently UCSD only - restart EEGLAB after changing that option.

repo: REST


text:
function grad = B_tmaxpool(prev_layer, curr_layer, future_layers)
input = prev_layer{1}.a;
[D,T,N] = size(input);
max_idx = curr_layer.idx;

if isfield(curr_layer, 'context')
    context = curr_layer.context;
else
    context = 0;
end
if isfield(curr_layer, 'stride')
    stride = curr_layer.stride;
else
    stride = 0;
end

future_grad = GetFutureGrad(future_layers, curr_layer);

if strcmpi(class(future_grad), 'gpuArray')
    grad = gpuArray.zeros(D,T,N);
else
    grad = zeros(D,T,N);
end

if context==0 || stride==0  % global pooling
    if 0
        for n=1:N
            for d = 1:D
                grad(d,max_idx(d,1,n),n) = future_grad(d,1,n);
            end
        end
    elseif 0
        for n = 1:N
            % for new Matlab version, we can use sub2idx.
            % idx2 = sub2idx(size(grad(:,:,1)), 1:D, max_idx(:,1,n));
            offset = (max_idx(:,1,n)-1)*D + D*T*(n-1);
            idx = offset+ [1:D]';
            grad(idx) = future_grad(:,1,n);
        end
    else
        offset = bsxfun(@plus, (squeeze(max_idx)-1)*D, D*T*(0:(N-1)) );
        offset = reshape(offset, numel(future_grad),1);
        idx = offset+ repmat((1:D)', N,1);
        grad(idx) = future_grad;
    end
else
    if 0
        for i=1:size(max_idx,2)
            offset = (i-1)*stride;
            for n=1:N
                for d = 1:D
                    idx = max_idx(d,i,n) + offset;
                    grad(d,idx,n) = future_grad(d,i,n);
                end
            end
        end
    else
        T2 = size(max_idx,2);
        max_idx = reshape(permute(max_idx, [1 3 2]), D*N,T2);
        future_grad = reshape(permute(future_grad, [1 3 2]), D*N,T2);
        grad = reshape(permute(grad, [1 3 2]), D*N,T);

        if 0
            for i=1:T2
                offset = (i-1)*stride;
                offset2 = (offset+max_idx(:,i)-1)*D*N;
                idx = offset2 + [1:(D*N)]';
                grad(idx) = future_grad(:,i);
            end
        else
            offset = (0:(T2-1))*stride;
            offset2 = bsxfun(@plus, (max_idx-1)*D*N, offset*D*N);
            offset2 = reshape(offset2, numel(offset2),1);
            idx = double(offset2) + repmat([1:(D*N)]', T2,1);   % have to use double precision as sometimes the index can be very big

            if stride<context
                % some idx may be redundant due to the fact that the same
                % element may go to multiple pooled elements.
                idx = reshape(idx,D*N,T2);
                for t=1:T2
                    grad(idx(:,t)) = grad(idx(:,t)) + future_grad(:,t);
                end
            else
                grad(idx) = future_grad;
            end
        end
        grad = permute(reshape(grad, D,N,T), [1 3 2]);
    end
end

end

author: singaxiong

repo: SignalGraph


text:
function endianFlag = LittleEndianCheck(mrSESSION);
%
% endianFlag = LittleEndianCheck(mrSESSION);
%
% Code for Initializing data collected at Stanford Lucas Center:
% Checks whether functional data is likely to use big- or
% little-endian format. The issue is this: around March 1,
% 2005, the 3T scanner was updated, and Gary Clover's spiral
% recon code updated, to use Little-Endian format. Prior to
% this, it was big-endian. This code simply checks the
% date stored for the first functional in mrSESSION, compares
% it against March 1, and sets the endian flag appropriately
% for downstream reconning.
%
% However, since there's a chance older data may be reconned
% using the newer little-endian scripts, if the date is
% before 03/01/05, the code prompts the user first.
%
% If endianFlag is 1, downstream code will treat functionals
% as little-endian; otherwise, will treat it as big-endian.
%
% ras, 05/05
if ~checkfields(mrSESSION,'functionals','date')
    endianFlag = 1; % little endian
    return
end

date = mrSESSION.functionals(1).date;
mo = str2num(date(1:2));
da = str2num(date(4:5));
yr = str2num(date(7:8));

if (yr==5 & mo<3) | (yr<5)
    % older data: prompt
    q = ['This data is before March 2005, so is likely ' ...
         'in Big-Endian format. However, it''s possible ' ...
         'that it''s in Little-Endian format, e.g. if the P*.7 ' ...
         'files were reconned later. How should the data ' ...
         'be read?   (type "help LittleEndianCheck" for more info)'];
    resp = questdlg(q,'mrInitRet','LittleEndian','BigEndian',...
                      'Quit','BigEndian');
    switch lower(resp)
        case 'littleendian', endianFlag = 1;
        case 'bigendian', endianFlag = 0;
        otherwise, error('User aborted');
    end
else
    % newer data: little endian
    endianFlag = 1;
end

return

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function ftps = kmph2ftps(kmph)
%KMPH2FTPS Convert speed from kilometers per hour to feet per second
%
%  ftps = KMPH2FTPS(kmph) converts speeds from kilometers per hour to feet
%   per second.
%
%  See also KMPH2KTS, KMPH2MPH, KMPH2MPS, FTPS2KMPH.

% Jonathan Sullivan
% Original: May 2011
% jonathan.sullivan@ll.mit.edu

ftps = kmph/1.09728;


function  prices = SABR_EurBarAmer_func(call, M, T, S0, Kvec, r, CTMCParams, ModParams, contract_type, L)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Pricing Function for European, American, and Barrier Options using
% double Layer CTMC approximation for SABR
%
% Models Supported: SABR
% Returns: price of contract (for vector of strikes)
% Author: Justin Lars Kirkby
%
% References:  (1) General Valuation Framework for SABR and Stochastic Local Volatility
%                   Models. SIAM J. Financial Mathematics, 2018. (w/ Z. Cui
%                   and D. Nguyen)
%
% ----------------------
% Contract/Model Params
% ----------------------
% call = 1 for call option, else Put
% contract_type = type of contact: % 1 = European, 2 = American, 3 = Down and Out Barrier
% Kvec  = strike vector
% S0 = initinal underlying value
% r   = interest rate (e.g. 0.05)
% T   = time remaining until maturity (in years, e.g. T=1)
% ModParams = model parameters: .v0, .alpha, .beta, .rho
% L =  For barrier contract, this is the barrier
%
% ----------------------
% Numerical (CTMC) Params
% ----------------------
% CTMCParams: .m_0 = grid/state size for variance process
%             .N = grid/state stize for underlying
%             .gridMult_v = grid non-uniformity multiplier (for variance)
%             .gridMult_s = grid non-uniformity multiplier (for underlying)
%             .gamma = Grid width param for variance grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m_0        = CTMCParams.m_0;
N          = CTMCParams.N;
gridMult_v = CTMCParams.gridMult_v;
gridMult_s = CTMCParams.gridMult_s;  %Grid mult param for S
gamma      = CTMCParams.gamma;          %Grid width param for variance grid

gridMethod_v = 5;   %%% ALWAYS use 5 for this one (puts v0 on grid)
gridMethod_s = 4;   %%% 5 puts S_0 on grid, but 4 seems better (requires interpolation)

v0    = ModParams.v0;
alpha = ModParams.alpha;
beta  = ModParams.beta;
rho   = ModParams.rho;

%%%%%%%%%%%%%%%%%%%%%%
%%%   Set Asset Grid bounds
%%%%%%%%%%%%%%%%%%%%%%
if S0 < 0.5
    ls = .01*S0;   %lower bound in asset Grid (S_t) space
else
    ls = 0.001*S0;
end

us = max(4.5*S0, S0 + 10*(v0*(S0)^beta)*sqrt(T));  %upper bound in asset Grid (S_t) space

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   Step 1: Variance Grid / Generators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt = T/M;
t = sqrt(T)/2;    %NOTE: this is different than we used to use...

mu_func = @(u) 0*u;
sig_func = @(u) alpha*u;
mu_H = v0;
sig2_H = v0^2*(exp(alpha^2*t) - 1);

lx = max(0.0001,mu_H - gamma*sqrt(sig2_H));
ux = mu_H + gamma*sqrt(sig2_H);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   Step 1: Variance Grid / Generators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
center_v = v0;
[Q,v] = General_Q_Matrix_Newest(m_0,mu_func,sig_func,lx,ux,gridMethod_v,center_v, gridMult_v);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   Step 2: Asset Grid / Grid For Xtilde
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g = @(s) (s ).^(1-beta)/(1-beta);
invOneBet = 1/(1-beta);

center_s = S0;   %center of asset grid (e.g. center points around the strike)
manualPoint_s = center_s;  %manually places S0 on grid
Xgrid = g(getNonUniformGrid(N, ls, us, gridMethod_s, center_s, manualPoint_s, gridMult_s)) - rho/alpha*v0;  %Grid for Xtilde

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   Step 3: Generators (for Xtilde process)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nm = N*m_0;
G = zeros(Nm, Nm);  %Big Generator matrix
I = eye(N,N);

Payoff = zeros(Nm, 1);  %Terminal payoff
sqrtRho = sqrt(1-rho^2);

for j = 1:m_0 %loop through rows of big G matrix
    %%%%%%%%%%
    % Step(1): Find G_j (generator with v(j) fixed)
    %%%%%%%%%%
    nu_j = v(j);
    muX_func_nu  = @(x) -.5*beta*(nu_j)^2.*((1-beta)*(x + rho*nu_j/alpha)).^(-1) ;  %Drift function of Xtilde with v(j) fixed
    sigX_func_nu = @(x) sqrtRho*nu_j*[x>-100];  %Constant function for each fixed variance state

    Gnu = getGenerator_Q_MatrixOnly(Xgrid, muX_func_nu, sigX_func_nu, gridMethod_v);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% FORCE absorbing vs reflecting
    Gnu(1,1) = 0; Gnu(1,2) =0;
    %Gnu(N,N) = 0; Gnu(N,N-1) = 0;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%
    % Step(2): Populate the Generator matrix (recall it is block tridiagonal,
    %%%%%%%%%%
    for k = max(1,j-1):min(m_0,j+1)  %%% NOTE: we skip the ones that are known to be zeros (matrix is block tridiagonal)
        lamjk = Q(j,k);
        if j==k  %diagonal block element of G matrix
            G((j-1)*N + 1:j*N, (k-1)*N + 1:k*N ) = Gnu + lamjk*I;
        else
            G((j-1)*N + 1:j*N, (k-1)*N + 1:k*N ) = lamjk*I;
        end
    end
end

%%%% Find bracketing variance gridpoint
k_0 = 2;
while v0 >= v(k_0) && k_0 < m_0
    k_0 = k_0+1;
end
k_0 = k_0 - 1;  %left bracketing point:  v(k_0) <= v0 < v(k_0 +1)

%%%% Find bracketing Xtilde gridpoint
x0 = g(S0) - rho*v0/alpha;  %initial value on Xtilde grid for (S_0, v0)
j_0 = 2;
while x0 >= Xgrid(j_0) && j_0 < N
    j_0 = j_0+1;
end
j_0 = j_0 - 1;  %left bracketing point:   Xtilde(j_0) <= x0 < Xtilde(j_0 +1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% VALUE : using recursive method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

P = expm(G*dt);   %transition matrix of one dimensional CTMC
initialIndex = (k_0-1)*N + j_0;  %index corresponding to initial conditions

prices = zeros(length(Kvec),1);

for k = 1:length(Kvec)
    K = Kvec(k);
    %%% Calculate Payoff
    if call == 1
        for j = 1:m_0
            Payoff((j-1)*N +1:j*N) = max(0, ((1-beta)*(max(0,Xgrid + rho*v(j)/alpha) )).^invOneBet - K) ;  %Portion of payoff corresponding to v(j)
        end
    else
       for j = 1:m_0
            Payoff((j-1)*N +1:j*N) = max(0, K - ((1-beta)*(max(0,Xgrid + rho*v(j)/alpha) )).^invOneBet );  %Portion of payoff corresponding to v(j)
       end
    end
    %%% Now Price
    if contract_type == 3  % Down and out for now
        %determine which states remain alive
        alive = zeros(Nm,1);
        for j = 1:m_0
            cons = g(L) - v(j)*rho/alpha;
            alive((j-1)*N +1:j*N) = (Xgrid > cons);
        end
        pVec = alive.*Payoff;
        for m = M-1:-1:0
            pVec = exp(-r*dt)*alive.*(P*pVec);
        end
    else %% either American or European
        pVec = exp(-r*dt)*(P*Payoff);  %initialize the value (at last period)
        if contract_type == 2
            for m = M-2:-1:0
                pVec = exp(-r*dt)*(P*pVec);  %continuation value
                pVec = max(pVec, Payoff);   %max of continuation and intrinsic value
            end
        elseif contract_type == 1
            for m = M-2:-1:0
                pVec = exp(-r*dt)*(P*pVec);  %continuation value
            end
        end
    end

    if gridMethod_s == 5  %Both v0 and x0 are on grid  (we assume v0 is a member of vol grid, ie as long as gridMethod_v = 5)
        prices(k) = pVec(initialIndex);
    elseif gridMethod_s == 4 %x0 is not on grid (though we assume v0 is, ie as long as gridMethod_v = 5)
        price1 = pVec(initialIndex);    %corresponds to j_0
        price2 = pVec(initialIndex+1);  %corresponds to j_0 + 1
        prices(k) = price1 + (price2 - price1)*(x0 - Xgrid(j_0))/(Xgrid(j_0+1) - Xgrid(j_0));  %LINEAR INTERPOLATION
    end
end

end


text:
function cnoise_test ( )

%*****************************************************************************80
%
%% CNOISE_TEST tests the CNOISE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CNOISE_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the CNOISE library.\n' );

  m = 10000;
  n = 10000;
  var = 1.0;
  alpha = 1.0;
  cnoise_test01 ( m, n, var, alpha );

  m = 10000;
  n = 10000;
  var = 1.0;
  r = 1.0;
  alpha = 1.0;
  cnoise_test02 ( m, n, var, r, alpha );

  m = 10000;
  n = 10000;
  r = 1.0;
  alpha = 1.0;
  cnoise_test03 ( m, n, r, alpha );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CNOISE_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function el = export_RT_DVH_module_field(args)
%"export_RT_DVH_module_field"
%   Given a planC DVH cell and a dose index, return a properly populated
%   RT_DVH_module field generated for that dose for use with RT Dose IODs.
%
%   For speed, tag must be a decimal representation of the 8 digit
%   hexidecimal DICOM tag desired, ie instead of '00100010', pass
%   hex2dec('00100010');
%
%   Arguments are passed in a structure, arg:
%       arg.tag         = decimal tag of field to fill
%       arg.data        = CERR structure(s) to fill from
%       arg.template    = an empty template of the module created by the
%                         function build_module_template.m
%
%   This function requires arg.data = {doseNum, DVHs};
%
%JRA 06/19/06
%
%Usage:
%   dcmobj = export_RT_DVH_module_field(args)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Init output element to empty.
el = [];

%Unpack input data.
tag         = args.tag;
doseNum     = args.data{1};
DVHs        = args.data{2};
template    = args.template;

switch tag
    case 806092896  %300C0060   Referenced Structure Set Sequence
        templateEl = template.getValue(tag);
        fHandle = @export_referenced_structure_set_sequence;

        %Create new null Sequence
        tmp = org.dcm4che3.data.Attributes;
        el = tmp.newSequence(tag, 0);

        for i=1:1 %Only a single element is allowed in this sequence.
            dcmobj = export_sequence(fHandle, templateEl, {DVHs(1)});
            el.add(i-1, dcmobj);
        end
        %return an attribute
        el = el.getParent();

    case 805568576  %30040040   DVH Normalization Point
        %Currently not implemented.

    case 805568578  %30040042   DVH Normalization Dose Value
        %Currently not implemented.

    case 805568592  %30040050   DVH Sequence
        templateEl = template.getValue(tag);
        fHandle = @export_DVH_sequence;

        %Create new null Sequence
        tmp = org.dcm4che3.data.Attributes;
        el = tmp.newSequence(tag, 0);

        dInd = find([DVHs.doseIndex] == doseNum);

        %Iterate over DVHs matching the passed dose index.
        for i=1:length(dInd)
            %Get the DVH we are exporting this cycle.
            myDVH = DVHs(dInd(i));

            %Build the sequence.
            dcmobj = export_sequence(fHandle, templateEl, {myDVH});

            %Add to sequence element.
            el.add(i-1, dcmobj);
        end
        %return an attribute
        el = el.getParent();

    otherwise
        warning(['No methods exist to populate DICOM RT_DVH module field ' dec2hex(tag,8) '.']);
end


text:
%TR2Q	Convert homogeneous transform to a unit-quaternion
%
%	Q = tr2q(T)
%
%	Return a unit quaternion corresponding to the rotational part of the
%	homogeneous transform T.
%
%	See also Q2TR

%	Copyright (C) 1993 Peter Corke
function q = tr2q(t)
	q = zeros(1,4);
	q(1) = sqrt(trace(t))/2;
	kx = t(3,2) - t(2,3);	% Oz - Ay
	ky = t(1,3) - t(3,1);	% Ax - Nz
	kz = t(2,1) - t(1,2);	% Ny - Ox

	if (t(1,1) >= t(2,2)) & (t(1,1) >= t(3,3))
		kx1 = t(1,1) - t(2,2) - t(3,3) + 1;	% Nx - Oy - Az + 1
		ky1 = t(2,1) + t(1,2);			% Ny + Ox
		kz1 = t(3,1) + t(1,3);			% Nz + Ax
		add = (kx >= 0);
	elseif (t(2,2) >= t(3,3))
		kx1 = t(2,1) + t(1,2);			% Ny + Ox
		ky1 = t(2,2) - t(1,1) - t(3,3) + 1;	% Oy - Nx - Az + 1
		kz1 = t(3,2) + t(2,3);			% Oz + Ay
		add = (ky >= 0);
	else
		kx1 = t(3,1) + t(1,3);			% Nz + Ax
		ky1 = t(3,2) + t(2,3);			% Oz + Ay
		kz1 = t(3,3) - t(1,1) - t(2,2) + 1;	% Az - Nx - Oy + 1
		add = (kz >= 0);
	end

	if add
		kx = kx + kx1;
		ky = ky + ky1;
		kz = kz + kz1;
	else
		kx = kx - kx1;
		ky = ky - ky1;
		kz = kz - kz1;
	end
	nm = norm([kx ky kz]);
	if nm == 0,
		q = [1 0 0 0];
	else
		s = sqrt(1 - q(1)^2) / nm;
		q(2:4) = s*[kx ky kz];
	end