text:
function area=hypersphereArea(r,numDim)
%%HYPERSPHEREAREA Determine the surface area of a hypersphere of a given
%                 radius in a given number of dimensions.
%
%INPUTS:     r The radius of the hypersphere. r>0.
%       numDim The number of dimensions of the hypersphere. numDim>=1. One
%              dimensions is a line, two a circle, three a sphere, etc.
%
%OUTPUTS: area The area of the hypersphere.
%
%The surface area of a hypersphere is given in terms of the gamma function
%in [1].
%
%REFERENCES:
%[1] S. Li, "Concise formulas for the area and volume of a hyperspherical
%    cap," Asian Journal of Mathematics and Statistics, vol. 4, no. 1, pp.
%    66-70, 2011.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

area=2*pi^(numDim/2)/gamma(numDim/2)*r^(numDim-1);

end

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


text:
function [Wmat, destGridAtlas, destGridLoc, destGridOrient] = tess_interp_mixed( ResultsMat, WmatSurf, srcSurfMat, destSurfMat, sMriSrc, sMriDest, isInteractive )
% TESS_INTERP_MIXED: Compute an interpolation matrix between two surfaces for a mixed source model.
%
% WARNING: All the source files must be flattened first (GridAtlas.Grid2Source is an identity matrix).

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2016

% Parse inputs
if (nargin < 7) || isempty(isInteractive)
    isInteractive = 1;
end

% Check the source model
if isfield(ResultsMat, 'nComponents') && (ResultsMat.nComponents ~= 0)
    error('This function is aimed to process only mixed source models.');
end
% Check that the field GridAtlas is available
if ~isfield(ResultsMat, 'GridAtlas') || ~isfield(ResultsMat.GridAtlas, 'Scouts') || isempty(ResultsMat.GridAtlas.Scouts)
    error('The field GridAtlas is missing.');
end
% Get Structures atlas in the destination surface
iStructDest = find(strcmpi({destSurfMat.Atlas.Name}, 'Structures'));
if isempty(iStructDest)
    error('Destination surface must include an atlas "Structures".');
end
% % Get "Source model" atlas in the destination surface
% iModelDest = find(strcmpi({destSurfMat.Atlas.Name}, 'Source model'));

% GridAtlas: Must be updated to fit the destination surface
destGridAtlas = ResultsMat.GridAtlas;
destGridAtlas.Scouts(:) = [];

% Initialize interpolation matrices
Wmat = sparse(1, size(ResultsMat.GridLoc,1));
destGridLoc    = [];
destGridOrient = [];
iVert2Grid = [];
isFirstMniWarning = 1;
hFig1 = [];

% Process each region separately
for iScoutSrc = 1:length(ResultsMat.GridAtlas.Scouts)
    sScoutSrc = ResultsMat.GridAtlas.Scouts(iScoutSrc);
    sScoutSrc.Vertices = sScoutSrc.Vertices(:)';
    % Check if the source model has been flattened
    if (sScoutSrc.Region(3) ~= 'C')
        error('Source maps must be flattened before being projected with this function.');
    end
    % Get region in destination surface
    iScoutDest = find(strcmpi(sScoutSrc.Label, {destSurfMat.Atlas(iStructDest).Scouts.Label}));
    if isempty(iScoutDest)
        disp(['PROJECT> Warning: Structure not found in destination surface: ' sScoutSrc.Label]);
        continue;
    end
%     % If there is already an atlas "Source model", make sure the constraints match
%     if ~isempty(iModelDest)
%         iScoutDestModel = find(strcmpi(sScoutSrc.Label, {destSurfMat.Atlas(iModelDest).Scouts.Label}));
%         if isempty(iModelDest) || ~isequal(destSurfMat.Atlas(iModelDest).Scouts(iScoutDestModel).Region, sScoutSrc.Region)
%             error(['The atlas "Source model" has different constraints in the destination surface for region "' sScoutSrc.Label '".']);
%         end
%     end
    % Update scout structure
    sScoutDest = sScoutSrc;
    sScoutDest.Vertices = destSurfMat.Atlas(iStructDest).Scouts(iScoutDest).Vertices(:)';
    sScoutDest.Seed     = destSurfMat.Atlas(iStructDest).Scouts(iScoutDest).Seed;

    % Process depends on region type
    switch (sScoutSrc.Region(2))
        % Volume: Compute a grid
        case 'V'
            % Source grid: Get from the source file
            srcLoc = ResultsMat.GridLoc(sScoutSrc.GridRows,:);
            % Compute a grid based on the destination region
            destLoc = dba_anatmodel(sScoutDest.Vertices, sScoutDest, destSurfMat, 'vol');
            destOri = 0 * destLoc;
            % Convert to MNI coordinates
            srcLocMni  = cs_convert(sMriSrc,  'scs', 'mni', srcLoc);
            destLocMni = cs_convert(sMriDest, 'scs', 'mni', destLoc);
            if (isempty(srcLocMni) || isempty(destLocMni))
                if isFirstMniWarning
                    strWarning = 'For accurate results, compute the MNI transformation for both subjects before running this interpolation.';
                    if isInteractive
                        java_dialog('warning', strWarning);
                    else
                        disp(['PROJECT> Warning: ' strWarning]);
                        bst_report('Warning', 'process_project_sources', [], strWarning);
                    end
                    isFirstMniWarning = 0;
                end
                srcLocMni  = srcLoc;
                destLocMni = destLoc;
            end
            % Set the indices of the regions in the Atlas
            sScoutDest.GridRows = size(destGridLoc,1) + (1:size(destLoc,1));
            % Compute interpolation
            Wmat(sScoutDest.GridRows, sScoutSrc.GridRows) = bst_shepards(destLocMni, srcLocMni, 8, 0);

            % === DISPLAY ALIGNMENT ===
            if isInteractive && ~isequal(sMriSrc, sMriDest)
                % Close previous figures
                if ~isempty(hFig1) && ishandle(hFig1)
                    close(hFig1);
                end
                % Compute convex envelopes
                facesSrc  = convhull(srcLoc);
                facesDest = convhull(destLoc);
                % Display figure
                [hFig1, iDS, iFig, hPatch] = view_surface_matrix(destLoc, facesDest, .4, [1 0 0]);
                set(hPatch, 'EdgeColor', 'r');
                [hFig1, iDS, iFig, hPatch] = view_surface_matrix(srcLoc, facesSrc, .4, [], hFig1);
                set(hPatch, 'EdgeColor', [.6 .6 .6]);
                line(destLoc(:,1), destLoc(:,2), destLoc(:,3), ...
                     'LineStyle',   'none', ...
                     'Color',       [1 1 0], ...
                     'MarkerSize',  2, ...
                     'Marker',      'o');
                line(srcLoc(:,1), srcLoc(:,2), srcLoc(:,3), ...
                     'LineStyle',   'none', ...
                     'Color',       [0 0 1], ...
                     'MarkerSize',  2, ...
                     'Marker',      'o');
                set(hFig1, 'Name', [sScoutSrc.Label ' (source grid)']);
                drawnow;
            end

        % Surface: Use the surface-surface interpolation computed previously
        case 'S'
            destLoc = destSurfMat.Vertices(sScoutDest.Vertices, :);
            destOri = destSurfMat.VertNormals(sScoutDest.Vertices, :);
            % Set the indices of the regions in the Atlas
            sScoutDest.GridRows = size(destGridLoc,1) + (1:size(destLoc,1));
            % Re-use surface inteprolation
            Wmat(sScoutDest.GridRows, sScoutSrc.GridRows) = WmatSurf(sScoutDest.Vertices, sScoutSrc.Vertices);
            % Add the match of the vertices in the cortex surface and the GridLoc matrix
            iVert2Grid = [iVert2Grid; sScoutDest.Vertices', sScoutDest.GridRows'];

        otherwise
            disp(['PROJECT> Warning: Structure ignored: ' sScoutSrc.Label]);
            continue;
    end
    % Concatenate GridLoc and GridOrient
    destGridLoc    = [destGridLoc;    destLoc];
    destGridOrient = [destGridOrient; destOri];
    % Add scout to destination atlas
    destGridAtlas.Scouts(end+1) = sScoutDest;
end

% Create sparse conversion matrices between indices
% Vert2Grid: Correspondance between the vertices of the cortex surface and the GridLoc field
destGridAtlas.Vert2Grid   = logical(sparse(iVert2Grid(:,2), iVert2Grid(:,1), ones(size(iVert2Grid,1),1)));
% Grid2Source: Correspondance between the GridLoc and ImageGridamp matrices (identity matrix because the source maps must be flatten before projecting)
destGridAtlas.Grid2Source = logical(speye(size(destGridLoc,1)));

author: brainstorm-tools


text:
function IM = updateVMCInfluence(IM, structROIV, sampleRateV, updateM)
%"updateVMCInfluence"
%   Uses the VMC++ MC engine to populate the beamlet fields of an IM structure.
%   StructROIV is the list of structures, sampleRateV is a vector of sample
%   rates, one for each structROIV.
%   The doses are stored in sparse format.
%   The stored index is with respect to the structure mask registered to the
%   uniformized CT scan.
%
%
%      Rewrite June 06, JOD, to be a 'drop-in' alternative to the pencil
%      beam calculations.
%
%      Based on code by P. Lindsay, J. Alaly, and J. Deasy
%
%      Last modified: JJW, 21 June 2006
%                     JJW, 26 June 2006
%                     JC,  12 July 2006
%                           Implement to use 'IM.params.cutoffDistance'.
%                           Remove the assumption of SID = 100 cm.
%           JC, 26 Feb 2007
%               Fix bug, "the dose was shifted by one voxel". to.
%               offset(1:2)=offset(1:2)+yres/2;
%
% 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/>.

%  Flow:
%
%  Loop over beams
%    Loop over beamlets
%      Get MC2CT beamlet matrix.
%      Loop over structures
%            Get infl contribution from that beamlet to that structure
%                at that sample rate.
%            Store
%      end structure loop
%    End beamlet loop
%  End beam loop

%This is a different organizational structure compared to the PB generate influence
%flow.  In particular, we loop over beams first, not structures, because the MC
%calculation is always done to the entire CT matrix, and it is by far the slowest
%computation.

%==========Set up MC params=======================%

global planC
global stateS
indexS=planC{end};

%obtain associated scanNum for structures. It is assumed that all the
%structures are associated to same scan (which is checked in IMRTP.m)
scanNum = getStructureAssociatedScan(structROIV(1));

numBeams = length(IM.beams);

IMRTPGui('statusbar', 'Preparing VMC++ calculation');
percentDone = 0;
IMRTPGui('waitbar', percentDone);

%Create all path parameters.
IMRTPdir = fileparts(which('IMRTP'));
VMCPath         = fullfile(IMRTPdir , 'vmc++', '');
runsPath        = fullfile(VMCPath, 'runs', '')
phantomPath     = fullfile('.', 'phantoms', '');
energyPath      = fullfile('.', 'spectra', '');

%Create phantom default filename.
phantomFilename = fullfile(phantomPath, 'CERR_IMRT.ct');

%Set environment variables for PC.  Unix set at runtime.
if ispc
    dos(['"' fullfile(VMCPath, 'setx" ') 'vmc_home "' VMCPath '" -m']);
    dos(['"' fullfile(VMCPath, 'setx" ') 'vmc_dir "' runsPath '" -m']);
    dos(['"' fullfile(VMCPath, 'setx" ') 'xvmc_dir "' VMCPath '" -m']);

    [jnk, out1] = dos('set vmc_home');
    [jnk, out2] = dos('set vmc_dir');
    [jnk, out3] = dos('set xvmc_dir');

    %(1:end-1) to avoid CR at end of output1-3.
    if ~strcmpi(out1(1:end-1), ['vmc_home=' VMCPath]) | ~strcmpi(out2(1:end-1), ['vmc_dir=' runsPath ]) | ~strcmpi(out3(1:end-1), ['xvmc_dir=' VMCPath])
        error('Environment Variables were not properly set for VMC++.  CERR has set them--please restart Matlab and begin the calculation again.');
    end
end

%Calculation will be based on the resolution of uniformized CT dataset.
currDir = cd;
cd(VMCPath);
fillWater = 0;
[offset, bbox] = generateCT_uniform(getUniformizedCTScan(0,scanNum), phantomFilename, scanNum, fillWater);
cd(currDir);

scandr = bbox(2) - bbox(1) + 1;
scandc = bbox(4) - bbox(3) + 1;
scands = bbox(6) - bbox(5) + 1;

%Get x,y,z resolution.
[xV, yV, zV] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));
xres = abs(xV(2) - xV(1));
yres = abs(yV(2) - yV(1));
zres = abs(zV(2) - zV(1));

%Left hand side versus centre of each voxel
%offset(1:2)=offset(1:2)-yres/2;
% This is the fix, as in <generateCT_uniform.m> %THIS IS THE FIX::: Y=Y(end:-1:1);
offset(1:2)=offset(1:2)+yres/2;
offset(3:4)=offset(3:4)-xres/2;
offset(5:6)=offset(5:6)-zres/2;

%Setup the options for the VMC++ input.
VMCOpt = VMCOptInit;
phantomFilename(strfind(phantomFilename, '\')) = '/'; %VMC requires '/'.
VMCOpt.startGeometry.startXYZGeometry.phantomFile=phantomFilename;

%Set options from IM.params.VMC
VMCOpt.startMCControl.NCase                     = IM.params.VMC.NumParticles;
VMCOpt.startMCControl.NBatch                    = IM.params.VMC.NumBatches;
VMCOpt.startScoring.startDoseOptions.scoreDoseToWater = IM.params.VMC.scoreDoseToWater;
VMCOpt.startBeamletSource.monoEnergy            = IM.params.VMC.monoEnergy;
VMCOpt.startVarianceReduction.repeatHistory     = IM.params.VMC.repeatHistory;
VMCOpt.startVarianceReduction.splitPhotons      = IM.params.VMC.splitPhotons;
VMCOpt.startVarianceReduction.photonSplitFactor = IM.params.VMC.photonSplitFactor;
VMCOpt.startQuasi.base                          = IM.params.VMC.base;
VMCOpt.startQuasi.dimension                     = IM.params.VMC.dimension;
VMCOpt.startQuasi.skip                          = IM.params.VMC.skip;
switch lower(IM.params.VMC.includeError)
    case 'yes'
        VMCOpt.startScoring.startOutputOptions.dumpDose = 1;
    case 'no'
        VMCOpt.startScoring.startOutputOptions.dumpDose = 2;
    otherwise
        error('Invalid value for includeError property.');
end

% Precalc voxel indices for all structures
scanIndV=cell(max(structROIV),1);
for i = 1 : length(structROIV)
    if sampleRateV(i) ~= 1
        %Ensures that interpolative downsampling won't miss edge points:
        maskSingle3D = getSurfaceExpand(structROIV(i),0.5,1);
        if rem(log2(sampleRateV(i)),1) ~= 0
            error('Sample factor must (currently) be a power of 2.')
        end
        maskSample3D = getDown3Mask(maskSingle3D, sampleRateV(i), 1);
        %Get a mask of where to sample points
        tmpM = logical(maskSample3D) & maskSingle3D;
        clear maskSample3D maskSingle3D;
    else
        tmpM = getUniformStr(structROIV(i));
    end
    scanIndV{i} = find(tmpM);  %Indices with respect to the uniformized scan.
    clear tmpM
end

%===========Loop over beams=============%

% To generate independent datasets, need to set
% VMCOpt.startQuasi.skip to skip the numbers of histories
% used in previous calculations

beamletCounter = 1;  %Keep up with how many beamlets have been computed.
errorBeamlets = [];

indV = 1:(scandr*scandc*scands);
[r,c,s] = ind2sub([scandr,scandc,scands], indV);
r = r + bbox(1) - 1;
c = c + bbox(3) - 1;
s = s + bbox(5) - 1;
%indV = sub2ind(getUniformizedSize(planC), r, c, s);
indV = sub2ind(getUniformScanSize(planC{indexS.scan}(scanNum)), r, c, s);
clear r c s

% Added by JC July 11, 2006
% Used later in chop doseV/doseM off by 'IM.params.cutoffDistance'
[xM yM zM] = meshgrid(xV, yV, zV);
pM = [xM(indV)', yM(indV)', zM(indV)'];
clear xM yM zM

for beamIndex=1:numBeams

    if ~any(updateM(:,beamIndex)) % no goal changed for this beam
        continue
    elseif ~all(updateM(:,beamIndex))
        indstructV = find(updateM(:,beamIndex));
        if ~ismember('skin',lower({IM.goals(indstructV).structName})) && ...
                ismember('skin',lower({IM.goals.structName}))
            skinInd = find(strcmpi('skin',{IM.goals.structName}));
            recomputeDose = 0;
        else
            recomputeDose = 1;
        end
    else
        recomputeDose = 1;
    end

    if recomputeDose

        pb(:,1)=IM.beams(beamIndex).RTOGPBVectorsM_MC(:,1)+IM.beams(beamIndex).x;
        pb(:,2)=IM.beams(beamIndex).RTOGPBVectorsM_MC(:,2)+IM.beams(beamIndex).y;
        pb(:,3)=IM.beams(beamIndex).RTOGPBVectorsM_MC(:,3)+IM.beams(beamIndex).z;

        dx=IM.beams(beamIndex).beamletDelta_x;
        dy=IM.beams(beamIndex).beamletDelta_y;

        IC=[IM.beams(beamIndex).isocenter.x IM.beams(beamIndex).isocenter.y  IM.beams(beamIndex).isocenter.z];
        th=IM.beams(beamIndex).gantryAngle*pi/180;

        M=[cos(th) sin(th) 0; -sin(th) cos(th) 0; 0 0 1];
        Minv=[cos(th) -sin(th) 0; sin(th) cos(th) 0; 0 0 1];

        % translate to IC frame of reference
        pb=pb-repmat(IC, [length(pb), 1]);

        % anti rotate (i.e., to gantry angle of 0)
        pb=(Minv*pb')';

        % assuming SID=100
        % back project beams up by -50 cm
        %     s1=0.5;
        %     s2=50;
        %     s3=s1*0.5;
        % take SID from the IM struct. back project beams up by -50 cm
        s2 = 50;
        s1 = s2/IM.beams(beamIndex).isodistance;
        s3 = s1 * s1;

        pb(:,[1 3])=pb(:, [1 3])*s1;
        pb(:,2)=pb(:,2)+s2;

        for i=1:length(pb),
            pb_new{i}=[pb(i, 1)+dx*s3, pb(i, 2),  pb(i, 3)+dy*s3
                pb(i, 1)-dx*s3,  pb(i, 2),  pb(i, 3)+dy*s3
                pb(i, 1)-dx*s3,  pb(i, 2),  pb(i, 3)-dy*s3];
        end

        % rotate by the gantry angle
        for i=1:length(pb),
            pb_rot{i}=(M*pb_new{i}')' + repmat(IC, [3 1]);
        end

        % include the offset from the CT image
        virtualSource=[IM.beams(beamIndex).x-offset(3), offset(2)-IM.beams(beamIndex).y, IM.beams(beamIndex).z-offset(5)];

        IC=[IC(1)-offset(3) offset(2)-IC(2) IC(3)-offset(5)];

        for i=1:length(pb),
            pb_rot{i}=pb_rot{i}-repmat([offset(3) 0 offset(5)], [3 1]);
            pb_rot{i}(:,2)=offset(2)-pb_rot{i}(:,2);
        end

        % rounding because VMC++ requires the edges of the beamlet form
        % a square, by testing the dot product      %PUT error function here??
        for i=1:length(pb),
            temp=pb_rot{i};
            pb_rot{i}(:,1)=temp(:,2);
            pb_rot{i}(:,2)=temp(:,1);
            pb_rot{i}=round(pb_rot{i}*100000)/100000;
        end

        VMCOpt.startBeamletSource.virtualPointSourcePosition = [virtualSource(2) virtualSource(1) virtualSource(3)];

        %Write .vmc files for each pencil beam.
        for i=1:length(pb_rot),

            %Set beam spectrum file.
            switch IM.beams(beamIndex).beamEnergy
                case 6
                    energyFile = fullfile(energyPath, 'var_6MV.spectrum');
                case 15
                    energyFile = fullfile(energyPath, 'mohan_15MV.spectrum');
                case 18
                    energyFile = fullfile(energyPath, 'var_18MV.spectrum');
                otherwise
                    error('Invalid beam energy.'); %Eventually add custom beamEnergy files here.
            end

            %if a spectrum was entered in the GUI, use this one instead
            if ~strcmp(IM.params.VMC.spectrum,'') && ~isempty(IM.params.VMC.spectrum)
                if exist(fullfile(VMCPath,'spectra',IM.params.VMC.spectrum), 'file') == 2
                    energyFile = fullfile(energyPath,IM.params.VMC.spectrum);
                else
                    fprintf('\nERROR: Spectrum %s does not exist!\n\n',fullfile(VMCPath,'spectra',IM.params.VMC.spectrum));
                    return
                end
            end

            energyFile(strfind(energyFile, '\')) = '/'; %Must use '/' instead of '\' for VMC.

            VMCOpt.startBeamletSource.spectrum = energyFile;

            switch lower(IM.beams(beamIndex).beamModality)
                case 'photons'
                    VMCOpt.startBeamletSource.charge = 0;
                case 'electrons'
                    VMCOpt.startBeamletSource.charge = -1;
                otherwise
                    error('Invalid beamModality.');
            end

            %Create different rands 1:30000 for each PB.
            VMCOpt.startMCControl.RNGSeeds = [round(rand*30000)+1 round(rand*30000)+1];

            %Set pencil beam edges.
            VMCOpt.startBeamletSource.beamletEdges= [pb_rot{i}(1, :), pb_rot{i}(2,:),pb_rot{i}(3, :)];

            %Write .vmc file.
            outfile=['MCpencilbeam_', int2str(beamIndex), '_', int2str(i), '.vmc'];
            VMCInputC = makeVMCInput(VMCOpt, fullfile(runsPath, outfile));
        end

        npb = length(pb_rot);

        clear pb pb_new pb_rot

        switch VMCOpt.startScoring.startOutputOptions.dumpDose
            case 1
                precision = 4;
            case 2
                precision = 2;
            otherwise
                error('Incorrect value for dumpDose in VMCOpt. Try 1 or 2.')
        end

        sourceM = repmat([IM.beams(beamIndex).x, IM.beams(beamIndex).y, IM.beams(beamIndex).z], length(pM), 1);
        pRelM = pM - sourceM;

        %=========Loop over beamlets=============%
        for pbNum = 1:npb

            if pbNum > 1
                disp(['Computed ' int2str(pbNum-1) ' out of ' num2str(npb)]); pause(0.003);
                try
                    IMRTPGui('status', beamIndex, numBeams, ' CT scan ', pbNum - 1, npb);
                end
            end

            outfile=['MCpencilbeam_', int2str(beamIndex), '_', int2str(pbNum)];
            if ispc
                currPath = cd;
                cd(VMCPath);
                dos(['start /low /B /wait ' fullfile('.', 'bin', 'vmc_Windows.exe') ' ' outfile '']);

                %if exist('VMCOpt.denoising')
                %  if strcmpi(VMCOpt.denoising,'yes')
                %    niter=2;
                %    cerror=0.05;
                %    str=fullfile(VMCPath, 'runs', [outfile '_' VMCOpt.startScoring.startDoseOptions.scoreInGeometries '.dos'])
                %    dose = readDoseFile(fullfile(VMCPath, 'runs', [outfile '_' VMCOpt.startScoring.startDoseOptions.scoreInGeometries '.dos']), precision, scandr, scandc, scands);
                %    dose = reshape(dose, [scandr, scandc, scands]);
                %    maxValue = max(dose(:));
                %    mask = dose < (IM.params.Scatter.Threshold * maxValue);
                %%    lowDoses = zeros(size(dose));
                %    lowDoses(mask) = dose(mask);
                %    dose(mask) = 0;
                %    smoothLowDoses=anisodiff3d_miao(lowDoses, cerror, niter);  %denoising via anisotropic diffusion.
                %      clear mask
                %      clear lowDoses
                %    dose  = dose + smoothLowDoses;
                %    doseV = dose(:);
                %  end
                %end
                doseV = readDoseFile(fullfile('.', 'runs', [outfile '_' VMCOpt.startScoring.startDoseOptions.scoreInGeometries '.dos']), precision, scandr, scandc, scands);

            elseif isunix
                currPath = cd;
                cd(VMCPath);
                dos(['chmod u+x ' fullfile(VMCPath, 'vmc_wrapper')]);
                dos(['chmod u+x ' fullfile(VMCPath, 'bin', 'vmc_Linux.exe')]);
                dos(['dos2unix ' fullfile(VMCPath, 'vmc_wrapper')]); % to make sure endOfLine characters are correct
                dos([fullfile(VMCPath,'vmc_wrapper') ' ' VMCPath ' ' runsPath ' ' VMCPath ' ' fullfile(VMCPath, 'bin', 'vmc_Linux.exe') ' ' outfile]);
                doseV = readDoseFile(fullfile(VMCPath, 'runs', [outfile '_' VMCOpt.startScoring.startDoseOptions.scoreInGeometries '.dos']), precision, scandr, scandc, scands);

            end
            cd(currPath);

            % Added by JC July 11, 2006
            % Used later in chop doseV/doseM off by 'IM.params.cutoffDistance'
            RTOGPBVectorsM = IM.beams(beamIndex).RTOGPBVectorsM;
            PBV = RTOGPBVectorsM(pbNum,:);
            distV = pRelM * PBV';    %Each row is the dot product of (the un-normalized) source-to-calc point
            %direction and the PB unit vector.
            %Hence, this is the distance along the PB ray line for each dose calc point.

            qM = sourceM + distV * PBV;  %The last term gives the vector i, j, and k components of
            clear distV
            %distance to the depth of closest approach.
            %Each row of sourceM is the i,j,k position of the source.
            %So vector sum locates the positions of closest approach
            %to the dose calc points.

            rM  = pM - qM;    %Each row of pM is the x, y, z location of a dose calc point.
            clear qM
            %rM is then the vector pointing from the dose calc point to the point
            %at which the PB vector makes closest approach.

            %Trap points which are too far away to need dose calcs
            %Uses sepsq if possible, much faster.
            try
                rDistVSquared = sepsq([0 0 0]', rM');  %mex file speedup
                goV = rDistVSquared < (IM.params.cutoffDistance^2);
                %clear rDistVSquared
            catch
                rDistV = (rM(:,1).^2 + rM(:,2).^2 + rM(:,3).^2);    %if mex not there
                goV = rDistV < IM.params.cutoffDistance^2;         %Mod by JOD to speedup, Dec 05.
                %clear rDistV
            end

            % Cut off the dose by IM.params.cutoffDistance.
            doseV(~goV) = 0;
            clear goV rM

            % apply compression
            doseV = applyIMRTCompression(IM.params, doseV);

            % build 3D dose matrix
            doseM = zeros(getUniformScanSize(planC{indexS.scan}(scanNum)));
            if length(indV)==length(doseV)
                doseM(indV) = doseV;
            else
                errorBeamlets = [errorBeamlets beamletCounter];
            end

            %-----------Insert dose data into correct beamlet indices---------------------%
            for i = 1 : length(structROIV)
                beamlet = createIMBeamlet(doseM(scanIndV{i}), scanIndV{i}, beamIndex, 0);
                if isempty(IM.beams(beamIndex).beamlets)
                    IM.beams(beamIndex).beamlets = beamlet;
                else
                    IM.beams(beamIndex).beamlets(i,pbNum) = beamlet;
                end
                IM.beams(beamIndex).beamlets(i,pbNum).structureName = planC{planC{end}.structures}(structROIV(i)).structureName;
                IM.beams(beamIndex).beamlets(i,pbNum).sampleRate    = sampleRateV(i);
                IM.beams(beamIndex).beamlets(i,pbNum).strUID = planC{planC{end}.structures}(structROIV(i)).strUID;
            end
            beamletCounter = beamletCounter + 1;

        end

    else % obtain doseV from skin

        beamletsSkin = IM.beams(beamIndex).beamlets;
        npb = size(beamletsSkin,2);
        for PBNum = 1:npb
            dosebmtV  = double(beamletsSkin(skinInd,PBNum).influence);
            indV      = beamletsSkin(skinInd,PBNum).indexV;
            maxVal    = beamletsSkin(skinInd,PBNum).maxInfluenceVal;
            sizeParam = beamletsSkin(skinInd,PBNum).fullLength;
            if isfield(beamletsSkin, 'lowDosePoints')
                lowDosePoints = unpackLogicals(beamletsSkin(skinInd,PBNum).lowDosePoints, size(indV));
                doseScaledV(~lowDosePoints) = dosebmtV(~lowDosePoints) * (maxVal / (2^8 -1));
                doseScaledV(lowDosePoints) = dosebmtV(lowDosePoints) * (maxVal / (2^8 -1) / (2^8 -1));
            else
                doseScaledV = dosebmtV * (maxVal / (2^8 -1));
            end
            % build 3D dose matrix
            doseM = zeros(getUniformScanSize(planC{indexS.scan}(scanNum)));
            doseM(indV) = doseScaledV;
            doseScaledV = [];

            %-----------Insert dose data into correct beamlet indices---------------------%
            for i = 1 : length(structROIV)
                beamlet = createIMBeamlet(doseM(scanIndV{i}), scanIndV{i}, beamIndex, 0);
                if isempty(IM.beams(beamIndex).beamlets)
                    IM.beams(beamIndex).beamlets = beamlet;
                else
                    IM.beams(beamIndex).beamlets(i,PBNum) = beamlet;
                end
                IM.beams(beamIndex).beamlets(i,PBNum).structureName = planC{planC{end}.structures}(structROIV(i)).structureName;
                IM.beams(beamIndex).beamlets(i,PBNum).sampleRate    = sampleRateV(i);
                IM.beams(beamIndex).beamlets(i,PBNum).strUID = planC{planC{end}.structures}(structROIV(i)).strUID;
            end

        end

    end

    IMRTPGui('status', beamIndex, numBeams, ' CT scan ', npb, npb);

end

if length(errorBeamlets)>0
    disp(' ');
    disp('Error in generateVMCInfluence.m!');
    disp('The following beamlets were set to 0:');
    disp(errorBeamlets);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function out = sampleInSeconds(stimList,ISI,varargin)
% out = sampleInSeconds(stimList,ISI,varargin)
% input: stimList, output: stimlist sampled in .1 seconds
%                          OR sampled at your specified frequency

scale = ceil(ISI*10);

if nargin > 2
    scale = round(ISI/varargin{1});
end

numstim = size(stimList,1);

out = zeros(numstim*scale,1);

for i = 0:numstim-1

    out(i*scale+1,1) = stimList(i+1,1);

end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function visualize_approximate_dynamic_images(images)
% VISUALIZE_DYNAMIC_IMAGES

di = compute_approximate_dynamic_images(images) ;

di = di - min(di(:)) ;
di = 255 * di ./ max(di(:)) ;
image(uint8(di)) ;

author: hbilen

repo: dynamic-image-nets


text:
function [sigma q_a out_counts score lv_score word_score] = tamEStepNoMex(x,beta,alpha,a_log_prior,sigma)
%function [sigma q_a out_counts score lv_score word_score] = tamEStepNoMex(x,beta,alpha,a_log_prior,sigma)
%
%this contains a lot of optimizations that make it kind of complicated.
%maybe i should release a simplified version?
[W K A] = size(beta);
if nargin < 5
    if isscalar(alpha), sigma = ones(1,K) * (numel(x)  / K + alpha);
    else sigma = ones(1,K) * numel(x) / K + alpha; end
end
mydudes = x>0;
q_a = normalize_vec(ones(1,A));
myphi = 1/K * ones(sum(mydudes),K);
myx = x(mydudes);
mybeta = beta(mydudes,:,:);

%new_counts = spalloc(K,W,sum(mydudes)*K);
xrep = full(repmat(myx',1,K));
iter = newIterator(5,'thresh',1e-4,'debug',false);
%iter_delta = newDeltaIterator(1000,'thresh',1e-3);
iter_ctr = 0;
while ~iter.done
    sigma = (myx * myphi + alpha);
    dig_sig = digamma(sigma) - digamma(sum(sigma));
    warning off;
    myphi = logToSimplex2(tprod(mybeta,[1 2 -1],q_a,[3 -1],'n') + repmat(dig_sig,size(myx,2),1));
    [q_a log_q_a] = logToSimplex2(squeeze(tprod(mybeta,[-1 -2 1],myphi,[-1 -2],'n'))' + a_log_prior);
    warning on;
    %score it. this is slow, so we do it less often.
    if rem(iter_ctr,5) == 0
        word_score = 0;
        %new_counts(:,mydudes) = transpose(myphi.*xrep);
        my_new_counts = transpose(myphi.*xrep);
        %for a = 1:A
        %new_counts = new_counts * q_a(a);
        %    word_score = word_score + q_a(a) * scoreWords(my_new_counts',mybeta(:,:,a));
        %end
        %same as above, but faster
        word_score = q_a * tprod(my_new_counts,[-1 -2],mybeta,[-2 -1 1],'n');
        %using mybeta(:,:,1) looks weird, but actually this part is being ignored, we're only using the latent variable score. so this is legit.
        [ig ig2 lv_score] = scoreDoc(my_new_counts',mybeta(:,:,1),myphi,myx,sigma,alpha,dig_sig);
        lv_score = lv_score + q_a * (a_log_prior - log_q_a)';
        score = word_score + lv_score;

        iter = updateIterator(iter,score);
    end
    iter_ctr = iter_ctr + 1;
end
%compute output counts for real
if nargout >= 3
    out_counts = zeros(W,K,A);
    for k=1:K, for a=1:A, out_counts(mydudes,k,a) = myphi(:,k) .* myx' .* q_a(a); end; end
end
end

author: jacobeisenstein

repo: SAGE


text:
function [test_failed_s, test_failed_b, test_failed_l]=test_waveletfilters
%TEST_WAVELETFILTERS  Test the erbfilters filter generator

[f, fs] = gspi;
f=f(1:1000);
Ls = length(f);
scales = linspace(10,0.1,100);
scales = flip(scales);
fmin = 250;
fmax = 20000;
bins = 8;
M=100;
alpha = 1-2/(1+sqrt(5)); % 1-1/(goldenratio) delay sequence
delays = @(n,a) a*(mod(n*alpha+.5,1)-.5);

wavelettypes = {{'cauchy', 300},{'fbsp', 4, 3}, {'morse'} };

for ii = 1:numel(wavelettypes)
    %first call includes the option to set a starting frequency for the wavelet
    %frequency range (can be applied to others too)
    [g_scales,a_scales,~,L_scales, info_scales]=waveletfilters(Ls,scales, 'delay', delays, wavelettypes{ii}, 'single', 'uniform', 'redtar', 4, 'complex');
    [g_bins,a_bins,~,L_bins,info_bins] = waveletfilters(Ls,'bins', fs,fmin, fmax, bins,'delay', delays, wavelettypes{ii}, 'repeat', 'uniform', 'startfreq', 800);
    [g_linear,a_linear,~,L_linear,info_linear] = waveletfilters(Ls,'linear', fs,fmin, fmax, M,'delay', delays, wavelettypes{ii}, 'uniform');

    Lscales = filterbanklength(L_scales, a_scales);
    Lbins = filterbanklength(L_bins, a_bins);
    Llinear = filterbanklength(L_linear, a_linear);

    gd_scales=filterbankdual(g_scales,a_scales,Lscales, 'asfreqfilter');
    gd_bins_asf=filterbankrealdual(g_bins,a_bins,Lbins, 'asfreqfilter');
    %gd_bins_e=filterbankrealdual(g_bins,a_bins,Lbins, 'econ');
    gd_linear=filterbankrealdual(g_linear,a_linear,Llinear, 'asfreqfilter');

    if 0
        % Inspect it: Dual windows, frame bounds and the response
        disp('Frame bounds scales:')
        [A,B]=filterbankbounds(g_scales,a_scales,Lscales);
        A
        B
        B/A
        filterbankresponse(g_scales,a_scales,Lscales,'real','plot');
        disp('Frame bounds bins:')
        [A,B]=filterbankrealbounds(g_bins,a_bins,Lbins);
        A
        B
        B/A
        filterbankresponse(g_bins,a_bins,Lbins,'real','plot');
        disp('Frame bounds linear:')
        [A,B]=filterbankrealbounds(g_linear,a_linear,Llinear);
        A
        B
        B/A
        filterbankresponse(g_scales,a_scales,Lscales,'real','plot');
        hold on
        filterbankresponse(g_bins,a_bins,Lbins,'real','plot');
        filterbankresponse(g_linear,a_linear,Llinear,'real','plot');
        figure; filterbankfreqz(g_scales,a_scales,Ls,fs,'plot','linabs','posfreq');
        figure; filterbankfreqz(g_bins,a_bins,Ls,fs,'plot','linabs','posfreq');
        figure; filterbankfreqz(g_linear,a_linear,Ls,fs,'plot','linabs','posfreq');
    end

    %reconstruct them all and calculate the norm to input
    c_scales = filterbank(f, g_scales, a_scales);
    fhat_scales = 2*real(ifilterbank(c_scales, gd_scales, a_scales));

    c_bins = filterbank(f, g_bins, a_bins);
    fhat_bins = 2*real(ifilterbank(c_bins, gd_bins_asf, a_bins));

    c_linear = filterbank(f, g_linear, a_linear);
    fhat_linear = 2*real(ifilterbank(c_linear, gd_linear, a_linear));

    if length(fhat_scales) > length(f)
        res_scales=norm(fhat_scales(1:length(f)) - f);
    else
        res_scales=norm(f(1:length(fhat_scales)) - fhat_scales);
    end
    if length(fhat_bins) > length(f)
        res_bins=norm(fhat_bins(1:length(f)) - f);
    else
        res_bins=norm(f(1:length(fhat_bins)) - fhat_bins);
    end
    if length(fhat_linear) > length(f)
        res_linear=norm(fhat_linear(1:length(f)) - f);
    else
        res_linear=norm(f(1:length(fhat_linear)) - fhat_linear);
    end

    %[test_failed_s,fail_s]=ltfatdiditfail(res_scales,0, 400);
    test_failed_s = 0;
    [test_failed_b,fail_b]=ltfatdiditfail(res_bins,0, 0.00001);
    [test_failed_l,fail_l]=ltfatdiditfail(res_linear,0, 0.00001);
    %s_scales=sprintf(['WAVELETFILTER DUAL SCALES:%3i %0.5g %s'],res_scales,fail_s);
    s_bins=sprintf(['WAVELETFILTER DUAL BINS:%3i %0.5g %s'],res_bins,fail_b);
    s_linear=sprintf(['WAVELETFILTER DUAL LINEAR:%3i %0.5g %s'],res_linear,fail_l);
    %disp(s_scales);
    disp(s_bins);
    disp(s_linear);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ U, g ] = cgrid_uv2rho(ncRef, uname, vname, hname, aname,...
    itime, klev, jj, ii)
% CGRID_UV2RHO Get 2D slice of velocity from a C-grid model (e.g. ROMS)
% Usage:
%    [U,g] = cgrid_uv2rho(ncRef, uname, vname, hvar, avar, itime, klev, [jj], [ii]);
% Inputs:
%   ncRef = OpenDAP Data URL, Local file, or existing 'ncgeodataset' object
%   uname - u variable to slice
%   vname - v variable to slice
%   hname - rho-centered variable to slice (e.g. 'h')
%   aname - rho-centered angle variable (e.g. 'angle')
  %   itime = time step to extract data. Use -1 for last time step. (ignored for static variables)
  %   klev  = vertical level  Use -1 for last level. (ignored for 2d variables)
%   [jj] =  jindex range (subset) [optional, default is entire dimension]
%   [ii] =  iindex range (subset) [optional, default is entire dimension]
%
% Outputs:
%   U = 2D array of complex velocity
%   g  = grid structure containing lon,lat,z,time (Matlab datenum)
%
% Example:
% url='http://geoport.whoi.edu/thredds/dodsC/examples/bora_feb.nc';
% hname='h'; uname='u'; vname='v'; aname='angle';
% nc=ncgeodataset(url);
% itime=3; % 3rd time step
% klev=-1; % last (top) layer
%  % Whole domain:
% [ U,g ] = cgrid_uv2rho(nc,uname,vname,hname,aname,itime,klev);
% % Subset to jj,ii range:
% [ U,g ] = cgrid_uv2rho(nc,uname,vname,hname,aname,itime,klev,1:58,1:70);
% % 3D variable:
% hname='h'; uname='ubar'; vname='vbar'; aname='angle';
% [ U,g ] = cgrid_uv2rho(nc,uname,vname,hname,aname,itime,klev,1:58,1:70);
% pcolorjw(g.lon,g.lat,abs(U));colorbar;dasp(44);
% arrows(g.lon(1:2:end,1:2:end),g.lat(1:2:end,1:2:end),...
%     U(1:2:end,1:2:end),0.08,'black');
%  title(datestr(g.time));dasp(44);

% Note: this routine assumes that the C grid has variables arranged:
%
%       ---------------------------------
%       rho | u | rho | u | rho | u | rho
%       ---------------------------------
%        v  |   | {v} |   | {v} |   |  v
%       ---------------------------------
%       rho |{u}|(rho)|{u}|(rho)|{u}| rho
%       ---------------------------------
%        v  |   | {v} |   | {v} |   |  v
%       ---------------------------------
%       rho | u | rho | u | rho | u | rho
%       ---------------------------------
%
% .. so that size(rho)=ny,nx;  size(u)=ny,nx-1, size(v)=ny-1,nx)
% Also assumes coordinate dimensions are arranged (t, z, y, x)

if (isa(ncRef, 'ncgeodataset')) %check for ncgeodataset object
    nc = ncRef;
else
    % open CF-compliant dataset as ncgeodataset object
    nc = ncgeodataset(ncRef);
end

uvar = nc.geovariable(uname);
vvar = nc.geovariable(vname);
hvar = nc.geovariable(hname);
avar = nc.geovariable(aname);

% depth from uvar
if nargin < 7
    usize = size(uvar);
  if ( length(usize) == 4)
    klev = usize(2);
  else
    klev = 1;
  end
end

% get lon,lat size from hvar
if nargin < 8
    hsize = hvar.size;
    horder = hvar.getaxesorder;
    if horder.lat == horder.lon
        jj = 1:hsize(horder.lat);
        ii = 1:hsize(horder.lat  + 1);
    else
        %         ii = hsize(order.lat);
        %         jj = hsize(order.lon);
        error('order.lat != order.lon');
    end
else
    %pass
end
ujj = (jj(1) + 1):(jj(end) - 1);
uii = ii(1):ii(end) - 1;
vjj = jj(1):jj(end) - 1;
vii = (ii(1) + 1):(ii(end) - 1);

if (length(size(uvar))==4)
   u = double(squeeze(uvar.data(itime, klev, ujj, uii))); % get u
   v = double(squeeze(vvar.data(itime, klev, vjj, vii))); % get v
elseif (length(size(uvar))==3)
   u = double(squeeze(uvar.data(itime, ujj, uii))); % get u
   v = double(squeeze(vvar.data(itime, vjj, vii))); % get v
elseif (length(size(uvar))==2)
   u = double(squeeze(uvar.data(ujj, uii))); % get u
   v = double(squeeze(vvar.data(vjj, vii))); % get v
end
g = hvar.grid_interop(jj, ii); % get lon,lat from rho-point variable (like 'h' in ROMS)
ang = avar.data(jj, ii); % get angle
U = ones(size(g.lon)) * nan; % template for U at rho points
U(2:end-1, 2:end-1) = complex(av2(u.').', av2(v)); %average u,v to rho
U = U .* exp(sqrt(-1) * ang); % rotate

if nargout == 2
    tim = uvar.timewindowij( double(itime) );
    g.time = tim.time;
    g.klevel = klev;
    g.itime = itime;
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
%--- help for generic/pull_objective ---
%
%  Pulls the objective function to optimize
%
%  ::
%
%    [ff,lb,ub]=pull_objective(obj)
%    [ff,lb,ub]=pull_objective(obj,varargin)
%    [ff,lb,ub,x0]=pull_objective(obj,varargin)
%    [ff,lb,ub,x0,vcov]=pull_objective(obj,varargin)
%    [ff,lb,ub,x0,vcov,obj]=pull_objective(obj,varargin)
%
%  Args:
%     obj (rise | dsge | svar | rfvar): initial model object
%     varargin (pairwise addional inputs): usual RISE arguments
%
%  Returns:
%     :
%
%     - **ff** [function handle]: for computing "minus log posterior kernel"
%     - **lb** [d x 1 vector]: lower bound of the parameters to optimize
%     - **ub** [d x 1 vector]: upper bound of the parameters to optimize
%     - **x0** [d x 1 vector]: posterior mode if available
%     - **vcov** [d x d matrix]: covariance matrix at the posterior mode if
%       available.
%     - **obj** [rise|dsge|svar|rfvar]: updated model object
%
%  Note:
%
%     - The function can be used for :
%
%        - optimization,
%        - gradient computation,
%        - hessian computation,
%        - posterior simulation
%
%     - The updated object should be used for doing various exercises (irfs,
%       simulations, etc.) if the posterior mode is not computed.
%
%     - Using this function is potentially costly, one could alternatively
%       simply use log_posterior_kernel. However, if there are restrictions, they
%       will not be enforced. Nevertheless it is an interesting proposition that
%       should be investigated further.
%
%
%    Other functions named pull_objective
%
%       abstvar/pull_objective    dsge/pull_objective
%

author: jmaih


text:
function [T, M, S, DISTR, df] = dtiTTestStat(g1, g2, Y, mask)

% Computes voxel-wise T-test statistics for two groups from a data array.
%
%   [T, M, S, DISTR, df] = dtiTTestStat(g1, g2, DT_ARRAY, [MASK])
%
% Input:
%   g1, g2      List of indices that correspond to each group out of 1:N
%                   E.g: g1 = 1:7, g2 = 8:14, N = 14
%   DT_ARRAY    Data array of size XxYxZxN (or nxN), where X, Y, Z are the volume
%                   dimensions and N is the number of subjects.
%                   (n is the number of voxels).
%   MASK        Optional XxYxZ binary array. Values of M and S are computed
%                   where mask = 1; in other voxels, M and S are set to 0.
%                   Default is entire volume.
%
% Output:
%   T           XxYxZx1 array of test statistics (0 where mask = 0)
%   M           XxYxZx2 array of means for both groups (0 where mask = 0)
%   S           XxYxZx1 array of pooled standard deviations (0 where mask = 0)
%   DISTR       The string 't'
%   df          The number of degrees of freedom = N-2
%
% Utilities:    ndfun.m, dtiSplitTensor.m, dti33to6.m
%
% WARNING: If using Pentium 4, eliminate NaN's from array before running
% (processor bug).
%
% Copyright by Armin Schwartzman, 2005

% HISTORY:
%   2004.06.23 ASH (armins@stanford.edu) wrote it.
%

% Check inputs
if (ndims(Y)==2 | ndims(Y)==3),
    Ind = 1;    % Data in indexed nxN format
    Y = shiftdim(Y, -2);
else
    Ind = 0;    % Data in XxYxZxN format
end
if (ndims(Y)<4 | ndims(Y)>5),
    error('Wrong input format');
end
if (~exist('mask')),
    mask = ones([size(Y,1) size(Y,2) size(Y,3)]);
end

% Computations
N1 = length(g1);
N2 = length(g2);
N  = N1 + N2;

M = cat(4, mean(Y(:,:,:,g1), 4), mean(Y(:,:,:,g2), 4));
Ystd1 = std(Y(:,:,:,g1), 1, 4); Ystd1(~mask) = 1;
Ystd2 = std(Y(:,:,:,g2), 1, 4); Ystd2(~mask) = 1;
S = sqrt((N1*Ystd1.^2 + N2*Ystd2.^2)/(N-2) .* (1/N1 + 1/N2));
T = (M(:,:,:,1) - M(:,:,:,2)) ./ S;

% Adjust output
if Ind,
    T = shiftdim(T, 2);
    M = shiftdim(M, 2);
    S = shiftdim(S, 2);
end

DISTR = 't';
df = [N-2 N-2];     % Second entry is dummy

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ text ] = tapas_BLT_initVisuals(params)
% tapas_BLT_initVisuals initializes Textures for all task relevant images
%
% Inputs:
%   params: data storage file
%
% Outputs:
%	text: image textures
%

% Set prediction cue images according to experiment mode and cue type
switch params.expMode
    case {'train_4', 'train_6'}
        if params.cuetype ==  1
            params.stimuli.cue1         = fullfile(params.path.imgfolder,'cue3_predYL_img.png');
            params.stimuli.cue1Scale    = fullfile(params.path.imgfolder,'cue3_predScale.png');
            params.stimuli.cue1SL       = fullfile(params.path.imgfolder,'cue3_predYL_selL_img.png');
            params.stimuli.cue1SR       = fullfile(params.path.imgfolder,'cue3_predYL_selR_img.png');
            params.stimuli.cue2         = fullfile(params.path.imgfolder,'cue4_predYL_img.png');
            params.stimuli.cue2Scale    = fullfile(params.path.imgfolder,'cue4_predScale.png');
            params.stimuli.cue2SL       = fullfile(params.path.imgfolder,'cue4_predYL_selL_img.png');
            params.stimuli.cue2SR       = fullfile(params.path.imgfolder,'cue4_predYL_selR_img.png');
        elseif params.cuetype ==  2
            params.stimuli.cue1         = fullfile(params.path.imgfolder,'cue3_predYR_img.png');
            params.stimuli.cue1Scale    = fullfile(params.path.imgfolder,'cue3_predScale.png');
            params.stimuli.cue1SL       = fullfile(params.path.imgfolder,'cue3_predYR_selL_img.png');
            params.stimuli.cue1SR       = fullfile(params.path.imgfolder,'cue3_predYR_selR_img.png');
            params.stimuli.cue2         = fullfile(params.path.imgfolder,'cue4_predYR_img.png');
            params.stimuli.cue2Scale    = fullfile(params.path.imgfolder,'cue4_predScale.png');
            params.stimuli.cue2SL       = fullfile(params.path.imgfolder,'cue4_predYR_selL_img.png');
            params.stimuli.cue2SR       = fullfile(params.path.imgfolder,'cue4_predYR_selR_img.png');
        elseif params.cuetype ==  3
            params.stimuli.cue1         = fullfile(params.path.imgfolder,'cue4_predYL_img.png');
            params.stimuli.cue1Scale    = fullfile(params.path.imgfolder,'cue4_predScale.png');
            params.stimuli.cue1SL       = fullfile(params.path.imgfolder,'cue4_predYL_selL_img.png');
            params.stimuli.cue1SR       = fullfile(params.path.imgfolder,'cue4_predYL_selR_img.png');
            params.stimuli.cue2         = fullfile(params.path.imgfolder,'cue3_predYL_img.png');
            params.stimuli.cue2Scale    = fullfile(params.path.imgfolder,'cue3_predScale.png');
            params.stimuli.cue2SL       = fullfile(params.path.imgfolder,'cue3_predYL_selL_img.png');
            params.stimuli.cue2SR       = fullfile(params.path.imgfolder,'cue3_predYL_selR_img.png');
        elseif params.cuetype ==  4
            params.stimuli.cue1         = fullfile(params.path.imgfolder,'cue4_predYR_img.png');
            params.stimuli.cue1Scale    = fullfile(params.path.imgfolder,'cue4_predScale.png');
            params.stimuli.cue1SL       = fullfile(params.path.imgfolder,'cue4_predYR_selL_img.png');
            params.stimuli.cue1SR       = fullfile(params.path.imgfolder,'cue4_predYR_selR_img.png');
            params.stimuli.cue2         = fullfile(params.path.imgfolder,'cue3_predYR_img.png');
            params.stimuli.cue2Scale    = fullfile(params.path.imgfolder,'cue3_predScale.png');
            params.stimuli.cue2SL       = fullfile(params.path.imgfolder,'cue3_predYR_selL_img.png');
            params.stimuli.cue2SR       = fullfile(params.path.imgfolder,'cue3_predYR_selR_img.png');
        end
    case {'task', 'debug'}
        if params.cuetype ==  1
            params.stimuli.cue1         = fullfile(params.path.imgfolder,'cue1_predYL_img.png');
            params.stimuli.cue1Scale    = fullfile(params.path.imgfolder,'cue1_predScale.png');
            params.stimuli.cue1SL       = fullfile(params.path.imgfolder,'cue1_predYL_selL_img.png');
            params.stimuli.cue1SR       = fullfile(params.path.imgfolder,'cue1_predYL_selR_img.png');
            params.stimuli.cue2         = fullfile(params.path.imgfolder,'cue2_predYL_img.png');
            params.stimuli.cue2Scale    = fullfile(params.path.imgfolder,'cue2_predScale.png');
            params.stimuli.cue2SL       = fullfile(params.path.imgfolder,'cue2_predYL_selL_img.png');
            params.stimuli.cue2SR       = fullfile(params.path.imgfolder,'cue2_predYL_selR_img.png');
        elseif params.cuetype ==  2
            params.stimuli.cue1         = fullfile(params.path.imgfolder,'cue1_predYR_img.png');
            params.stimuli.cue1Scale    = fullfile(params.path.imgfolder,'cue1_predScale.png');
            params.stimuli.cue1SL       = fullfile(params.path.imgfolder,'cue1_predYR_selL_img.png');
            params.stimuli.cue1SR       = fullfile(params.path.imgfolder,'cue1_predYR_selR_img.png');
            params.stimuli.cue2         = fullfile(params.path.imgfolder,'cue2_predYR_img.png');
            params.stimuli.cue2Scale    = fullfile(params.path.imgfolder,'cue2_predScale.png');
            params.stimuli.cue2SL       = fullfile(params.path.imgfolder,'cue2_predYR_selL_img.png');
            params.stimuli.cue2SR       = fullfile(params.path.imgfolder,'cue2_predYR_selR_img.png');
        elseif params.cuetype ==  3
            params.stimuli.cue1         = fullfile(params.path.imgfolder,'cue2_predYL_img.png');
            params.stimuli.cue1Scale    = fullfile(params.path.imgfolder,'cue2_predScale.png');
            params.stimuli.cue1SL       = fullfile(params.path.imgfolder,'cue2_predYL_selL_img.png');
            params.stimuli.cue1SR       = fullfile(params.path.imgfolder,'cue2_predYL_selR_img.png');
            params.stimuli.cue2         = fullfile(params.path.imgfolder,'cue1_predYL_img.png');
            params.stimuli.cue2Scale    = fullfile(params.path.imgfolder,'cue1_predScale.png');
            params.stimuli.cue2SL       = fullfile(params.path.imgfolder,'cue1_predYL_selL_img.png');
            params.stimuli.cue2SR       = fullfile(params.path.imgfolder,'cue1_predYL_selR_img.png');
        elseif params.cuetype ==  4
            params.stimuli.cue1         = fullfile(params.path.imgfolder,'cue2_predYR_img.png');
            params.stimuli.cue1Scale    = fullfile(params.path.imgfolder,'cue2_predScale.png');
            params.stimuli.cue1SL       = fullfile(params.path.imgfolder,'cue2_predYR_selL_img.png');
            params.stimuli.cue1SR       = fullfile(params.path.imgfolder,'cue2_predYR_selR_img.png');
            params.stimuli.cue2         = fullfile(params.path.imgfolder,'cue1_predYR_img.png');
            params.stimuli.cue2Scale    = fullfile(params.path.imgfolder,'cue1_predScale.png');
            params.stimuli.cue2SL       = fullfile(params.path.imgfolder,'cue1_predYR_selL_img.png');
            params.stimuli.cue2SR       = fullfile(params.path.imgfolder,'cue1_predYR_selR_img.png');
        end
end

% Set stimulus answer images according to cue type
if params.cuetype ==  1 || params.cuetype ==  3
    params.stimuli.stimAnswer           = fullfile(params.path.imgfolder,'stimAnsYL_img.png');
    params.stimuli.stimAnswerSL         = fullfile(params.path.imgfolder,'stimAnsYL_selL_img.png');
    params.stimuli.stimAnswerSR         = fullfile(params.path.imgfolder,'stimAnsYL_selR_img.png');
elseif params.cuetype == 2 || params.cuetype ==  4
    params.stimuli.stimAnswer           = fullfile(params.path.imgfolder,'stimAnsYR_img.png');
    params.stimuli.stimAnswerSL         = fullfile(params.path.imgfolder,'stimAnsYR_selL_img.png');
    params.stimuli.stimAnswerSR         = fullfile(params.path.imgfolder,'stimAnsYR_selR_img.png');
end

% Add extra images
params.stimuli.predtimeout      = fullfile(params.path.imgfolder,'pred_timeout_img.png');
params.stimuli.prednoresponse   = fullfile(params.path.imgfolder,'pred_noResponse.png');
params.stimuli.stimulation      = fullfile(params.path.imgfolder,'stimulation.png');
params.stimuli.iti              = fullfile(params.path.imgfolder,'iti.png');
params.stimuli.break            = fullfile(params.path.imgfolder,'break.png');
params.stimuli.ratetimeout      = fullfile(params.path.imgfolder,'rate_timeout.png');
params.stimuli.endfieldmap      = fullfile(params.path.imgfolder,'endfieldmap.png');
params.stimuli.endfinal         = fullfile(params.path.imgfolder,'endfinal.png');

for i = 1:params.nIntroTrain
    params.stimuli.introTrain{i}  = fullfile(params.path.imgfolder,['Slide',num2str(i),'.png']);
    imgIntroTrain = imread(params.stimuli.introTrain{i});
    text.IntroTrain{i} = Screen('MakeTexture', params.screen.window, imgIntroTrain);
end

for i2 = 1:params.nIntroMain
    params.stimuli.introMain{i2}  = fullfile(params.path.imgfolder,['Slide',num2str(i+i2),'.png']);
    imgIntroMain = imread(params.stimuli.introMain{i2});
    text.IntroMain{i2} = Screen('MakeTexture', params.screen.window, imgIntroMain);
end

% read images
imgcue1                         = imread(params.stimuli.cue1);
imgcue2                         = imread(params.stimuli.cue2);
imgcue1Scale                    = imread(params.stimuli.cue1Scale);
imgcue2Scale                    = imread(params.stimuli.cue2Scale);
imgcue1SL                       = imread(params.stimuli.cue1SL);
imgcue1SR                       = imread(params.stimuli.cue1SR);
imgcue2SL                       = imread(params.stimuli.cue2SL);
imgcue2SR                       = imread(params.stimuli.cue2SR);
imgpredtimeout                  = imread(params.stimuli.predtimeout);
imgprednoresponse               = imread(params.stimuli.prednoresponse);
imgstimulation                  = imread(params.stimuli.stimulation);
imgstimanswer                   = imread(params.stimuli.stimAnswer);
imgstimanswerSL                 = imread(params.stimuli.stimAnswerSL);
imgstimanswerSR                 = imread(params.stimuli.stimAnswerSR);
imgiti                          = imread(params.stimuli.iti);
imgbreak                        = imread(params.stimuli.break);
imgratetimeout                  = imread(params.stimuli.ratetimeout);
imgendfieldmap                  = imread(params.stimuli.endfieldmap);
imgendfinal                     = imread(params.stimuli.endfinal);

% Make the images into a textures that can be drawn to the screen
text.imgcue1              = Screen('MakeTexture', params.screen.window, imgcue1);
text.imgcue2              = Screen('MakeTexture', params.screen.window, imgcue2);
text.imgcue1Scale         = Screen('MakeTexture', params.screen.window, imgcue1Scale);
text.imgcue2Scale         = Screen('MakeTexture', params.screen.window, imgcue2Scale);
text.imgcue1SL            = Screen('MakeTexture', params.screen.window, imgcue1SL);
text.imgcue1SR            = Screen('MakeTexture', params.screen.window, imgcue1SR);
text.imgcue2SL            = Screen('MakeTexture', params.screen.window, imgcue2SL);
text.imgcue2SR            = Screen('MakeTexture', params.screen.window, imgcue2SR);
text.imgpredtimeout       = Screen('MakeTexture', params.screen.window, imgpredtimeout);
text.imgprednoresponse    = Screen('MakeTexture', params.screen.window, imgprednoresponse);
text.imgstimulation       = Screen('MakeTexture', params.screen.window, imgstimulation);
text.imgstimanswer        = Screen('MakeTexture', params.screen.window, imgstimanswer);
text.imgstimanswerSL      = Screen('MakeTexture', params.screen.window, imgstimanswerSL);
text.imgstimanswerSR      = Screen('MakeTexture', params.screen.window, imgstimanswerSR);
text.imgiti               = Screen('MakeTexture', params.screen.window, imgiti);
text.imgbreak             = Screen('MakeTexture', params.screen.window, imgbreak);
text.imgratetimeout       = Screen('MakeTexture', params.screen.window, imgratetimeout);
text.endfieldmap          = Screen('MakeTexture', params.screen.window, imgendfieldmap);
text.endfinal             = Screen('MakeTexture', params.screen.window, imgendfinal);
end

author: translationalneuromodeling

repo: tapas


text:
% Q = vbfa(D, Y, W_module, X_module, noise_module, ...)
%
% Variational Bayesian (VB) factor analysis (FA) learning algorithm with
% changeable modules for the latent variables.
%
% Optional parameters:
% 'maxiter'             100
% 'update_x'            1
% 'update_w'            1
% 'update_noise'        1
% 'rotate'              1
% 'rotation_checkgrad'  false
% 'rotation_show'       false
% 'rotation_maxiter'    10
% 'debug'               false
% 'autosave'            false
% 'autosavefile'        'vbfa_autosave'

function Q = vbfa(D, Y, W_module, X_module, noise_module, varargin)

[M,N] = size(Y);

options = struct('maxiter', 100, ...
                 'update_x', 1, ...
                 'update_w', 1, ...
                 'update_noise', 1, ...
                 'initialize_x', true, ...
                 'initialize_w', true, ...
                 'initialize_noise', true, ...
                 'rotate', 1, ...
                 'rotation_checkgrad', false, ...
                 'rotation_show', false, ...
                 'rotation_maxiter', 10, ...
                 'loglikelihood', true, ...
                 'debug', false, ...
                 'autosave', false,...
                 'autosavefile', 'vbfa_autosave');

[options, errmsg] = argparse(options, varargin{:});
error(errmsg);

% Initialize
disp('Initializing variables..')
if options.initialize_w
  [Q.W,Q.CovW,Q.rhow] = W_module.initialize(D,M);
else
  QW = W_module.get_struct();
  Q.W = QW.X;
  Q.CovW = QW.CovX;
  Q.rhow = QW.rho;
end
if options.initialize_x
  [Q.X,Q.CovX,Q.rhox] = X_module.initialize(D,N);
else
  QX = W_module.get_struct();
  Q.X = QX.X;
  Q.CovX = QX.CovX;
  Q.rhox = QX.rho;
end
if options.initialize_noise
  [Q.Tau] = noise_module.initialize(M,N);
else
  QTau = noise_module.get_struct();
  Q.Tau = QTau.Tau;
end
%Q.W_module = W_module;
%Q.X_module = X_module;
%Q.noise_module = noise_module;

% Observed values
Obs = ~isnan(Y);
MN_obs = sum(Obs(:));

%%
%% VB learning

loglikelihood_old = -inf;
Q.loglikelihood = nan(options.maxiter,1);

Q.Y = Y;

% Replace missing values with zeros for computational reasons
Y(~Obs) = 0;

KL_W = -inf;
KL_X = -inf;
KL_Tau = -inf;

% Debugging stuff:
klw = nan(options.maxiter,1);
klx = nan(options.maxiter,1);
kltau = nan(options.maxiter,1);
loglike = nan(options.maxiter,1);

%f = zeros(4,options.maxiter);

disp('Starting the VB iteration..')
for ind=1:options.maxiter
  Q.ind = ind;

  startitercpu = cputime;

  %
  % Update variables
  %

  Tau = Q.Tau;
  Tau(~Obs) = 0;

  if index_selected(Q.ind, options.update_x)
    if index_selected(Q.ind, options.debug)
      disp('Update X')
    end
    %[U,Z] = projection(Y,Obs,W,CovW,Tau);
    [Q.X,Q.CovX,Q.rhox,logrhox,KL_X] = X_module.update(Q.ind, Y, Obs, Q.W, ...
                                                      Q.CovW, Q.rhow, Tau);
  end

  if index_selected(Q.ind, options.update_w)
    if index_selected(Q.ind, options.debug)
      disp('Update W')
    end
    %[U,Z] = projection(Y',Obs',X,CovX,Tau');
    [Q.W,Q.CovW,Q.rhow,logrhow,KL_W] = W_module.update(Q.ind, Y', Obs', Q.X, ...
                                                      Q.CovX, Q.rhox, Tau');
  end

  % Compute squared errors <(y_mn-w_m*x_n)^2>
  % (used by noise update and loglikelihood)
  E2 = zeros(size(Y));
  Yh = Q.W'*Q.X;
  E2 = Y.*Y;
  E2 = E2 - 2*Y.*Yh;
  E2 = spdiag(Q.rhow)*E2*spdiag(Q.rhox);
  if ndims(Q.CovW) == 2 && ndims(Q.CovX) == 2
    E2 = E2 ...
         + spdiag(Q.rhow)*(Yh.^2)*spdiag(Q.rhox) ...
         + spdiag(Q.rhow)*(Q.W'.^2)*Q.CovX ...
         + Q.CovW'*(Q.X.^2)*spdiag(Q.rhox) ...
         + Q.CovW'*Q.CovX;
  elseif ndims(Q.CovW) == 3 && ndims(Q.CovX) == 3
    xx = bsxfun(@times, reshape(Q.X*spdiag(Q.rhox),[D,1,N]), reshape(Q.X,[1,D,N])) ...
         + Q.CovX;
    xx = reshape(xx, [D*D,N]);
    ww = bsxfun(@times, reshape(Q.W*spdiag(Q.rhow),[D,1,M]), reshape(Q.W,[1,D,M])) ...
         + Q.CovW;
    ww = reshape(ww, [D*D,M]);
    E2 = E2 + ww'*xx;
  else
    % TODO: Optimize this..
    warning('Optimize this.. and is rho properly here??');
    for m=1:M
      for n=1:N
        if Obs(m,n)
          if ndims(Q.CovW) == 2
            ww = Q.W(:,m)*Q.W(:,m)' + diag(Q.CovW(:,m));
          else
            ww = Q.W(:,m)*Q.W(:,m)' + Q.CovW(:,:,m);
          end
          if ndims(Q.CovX) == 2
            xx = Q.X(:,n)*Q.X(:,n)' + diag(Q.CovX(:,n));
          else
            xx = Q.X(:,n)*Q.X(:,n)' + Q.CovX(:,:,n);
          end
          %WX_WX = WX_WX + traceprod(ww, xx);
          E2(m,n) = E2(m,n) + traceprod(ww, xx);
        end
      end
    end
  end
  E2(~Obs) = 0;

  if index_selected(Q.ind, options.update_noise)
    if index_selected(Q.ind, options.debug)
      disp('Update Tau')
    end

    [Q.Tau,LogTau,KL_Tau] = noise_module.update(Q.ind, E2, Obs);
%    [Tau,lowerbound] = noise_module.update(Y, Obs, v_W, W, CovW, v_X, X, CovX);
  end

  %
  % Rotate
  %

  if index_selected(Q.ind, options.rotate)

    % TODO: You could optimize the hyperparameters at the same time?
    disp('Rotating..')
    A = eye(D);

    if index_selected(Q.ind, options.rotation_checkgrad)
      mycheckgrad(@rotation_cost, A(:) + 0.5*randn(D^2,1), 1e-3, W_module, ...
                  X_module);
    end
    A = minimize(A(:), @rotation_cost, options.rotation_maxiter, W_module, ...
                 X_module);
    A = reshape(A, [D D]);
    if index_selected(Q.ind, options.rotation_show)
      A
    end
    [Q.W, Q.CovW] = W_module.rotate(A);
    [Q.X, Q.CovX] = X_module.rotate(inv(A)');

  end

  %
  % Evaluate VB lower bound
  %

  if index_selected(Q.ind, options.loglikelihood)
  % Likelihood part: <log p(Y|...)>
  logpdf_y = gaussian_logpdf(Q.Tau(Obs)'*E2(Obs), ...
                             0, ...
                             0, ...
                             -sum(LogTau(Obs))-logrhow'*sum(Obs,2)-sum(Obs,1)*logrhox, ...
                             MN_obs);

  % Debugging stuff
  klw(Q.ind) = -KL_W;
  klx(Q.ind) = -KL_X;
  kltau(Q.ind) = -KL_Tau;
  loglike(Q.ind) = logpdf_y;

  % Lower bound
  Q.loglikelihood(Q.ind) = logpdf_y - KL_W - KL_X - KL_Tau;

  if Q.loglikelihood(Q.ind) < loglikelihood_old
    warning(sprintf('Loglikelihood lower bound decreased relatively %e!', ...
                    (loglikelihood_old - Q.loglikelihood(Q.ind)) / ...
                    loglikelihood_old));
    % plot([Q.loglikelihood, klw, klx, kltau, loglike]);
  end

  fprintf('Iteration step %d: loglikelihood=%e (%.2f seconds)\n', Q.ind, ...
          Q.loglikelihood(Q.ind), cputime()-startitercpu);

  loglikelihood_old = Q.loglikelihood(Q.ind);
  end

  if index_selected(Q.ind, options.autosave)
    Q.W_struct = W_module.get_struct();
    Q.X_struct = X_module.get_struct();
    Q.Tau_struct = noise_module.get_struct();
    fprintf('Saving results to %s...', options.autosavefile);
    save(options.autosavefile, '-struct', 'Q');
    fprintf(' done.\n');
  end

end

function [c, dc] = rotation_cost(a, W_module, X_module)
N = sqrt(length(a));
A = reshape(a, N, N);
[U,S,V] = svd(A);
invS = diag(1./diag(S));
invAt = U*invS*V';
[c_w, dc_w] = W_module.rotation_cost(A, U, S, V);
[c_x, dc_x] = X_module.rotation_cost(invAt, U, invS, V);
dc_x = -invAt*dc_x'*invAt;
c = c_w + c_x;
dc = dc_w(:) + dc_x(:);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function result = simplex_nd ( func, n, v )

%*****************************************************************************80
%
%% SIMPLEX_ND approximates an integral inside a simplex in ND.
%
%  Integration region:
%
%    The simplex bounded by the origin and a convex combination of N points.
%
%  Discussion:
%
%    An N+1 point second degree formula is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function which evaluates F(X) at the N-dimensional point
%    X, of the form
%      function value = func ( n, x )
%
%    Input, integer N, the dimension of the space.
%
%    Input, real V(N+1,N), each of the
%    N+1 rows of V contains the N coordinates of one of the
%    vertices of the simplex.
%
%    Output, real RESULT, the approximate integral of the function.
%
  c = 1.0E+00 / sqrt ( n + 2 );
  w = 1.0E+00 / ( n + 1 );

  for i = 1 : n
    x(i) = w * ( 1.0E+00 - c ) * sum ( v(1:n+1,i) );
  end

  quad = 0.0E+00;

  for i = 1 : n+1

    x(1:n) = x(1:n) + c * v(i,1:n);

    quad = quad + w * feval ( func, n, x );

    x(1:n) = x(1:n) - c * v(i,1:n);

  end

  volume = simplex_volume_nd ( n, v );
  result = quad * volume;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
basefolder=fileparts( mfilename( 'fullpath' ) );
elapsedtime=zeros(1,5);

%Training setup
setup.scaling=3; % MAGI-ADAPT
setup.nclusters=2048; %default configuration nclusters=2048
setup.overlap=setup.scaling;%min=1 (sliding window), max=psize (no overlap)
setup.psize=3*setup.scaling;
setup.model='ibp';%can be replaced by 'bicubic', etc.
setup.ntrees=16; %for training
setup.nscales=16; %for training
setup.force_trees=true; %force training - trees
setup.force_regs=true; %force training - regressors
setup.nsamples=1500; %per leaf, to compute regressors
setup.trainingfolder=[ basefolder '/data/training' ];
setup.trainingext='*.bmp';
time=tic;

for i = 1:5
setup=nbsrf_training( setup );
elapsedtime(i,1)=toc(time);

end

average=sum(elapsedtime)/5;

author: thomas-koehler

repo: SupER


text:
classdef PTKSegmentLungs < MimGuiPlugin
    % PTKSegmentLungs. Gui Plugin for activating lung segmentation
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the Pulmonary Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Lungs'
        SelectedText = 'Lungs'
        ToolTip = 'Segment the lungs'
        Category = 'Segment region'
        Visibility = 'Dataset'
        Mode = 'Segment'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Icon = 'icon_lungs.png'
        Location = 2
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.RunPluginCallback('PTKLeftAndRightLungs');
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && (gui_app.IsCT || gui_app.IsMR);
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.IsDatasetLoaded && strcmp(gui_app.GetCurrentPluginName, 'PTKLeftAndRightLungs');
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function value = r4_besk1 ( x )

%*****************************************************************************80
%
%% R4_BESK1 evaluates the Bessel function K of order 1 of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the Bessel function K of order 1 of X.
%
  persistent bk1cs
  persistent ntk1
  persistent xmax
  persistent xmin
  persistent xsml

  if ( isempty ( ntk1 ) )

    bk1cs = [ ...
       0.0253002273389477705, ...
      -0.353155960776544876, ...
      -0.122611180822657148, ...
      -0.0069757238596398643, ...
      -0.0001730288957513052, ...
      -0.0000024334061415659, ...
      -0.0000000221338763073, ...
      -0.0000000001411488392, ...
      -0.0000000000006666901, ...
      -0.0000000000000024274, ...
      -0.0000000000000000070 ]';

    ntk1 = r4_inits ( bk1cs, 11, 0.1 * r4_mach ( 3 ) );
    xmin = exp ( max ( log ( r4_mach ( 1 ) ), ...
      - log ( r4_mach ( 2 ) ) ) + 0.01 );
    xsml = sqrt ( 4.0 * r4_mach ( 3 ) );
    xmax = - log ( r4_mach ( 1 ) );
    xmax = xmax - 0.5 * xmax * log ( xmax ) ...
      / ( xmax + 0.5 ) - 0.01;

  end

  if ( x <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_BESK1 = Fatal error!\n' );
    fprintf ( 1, '  X <= 0.\n' );
    error ( 'R4_BESK1 = Fatal error!' )
  elseif ( x <= xsml )
    y = 0.0;
    value = log ( 0.5 * x ) * r4_besi1 ( x ) + ( 0.75 ...
      + r4_csevl ( 0.5 * y - 1.0, bk1cs, ntk1 ) ) / x;
  elseif ( x <= 2.0 )
    y = x * x;
    value = log ( 0.5 * x ) * r4_besi1 ( x ) + ( 0.75 ...
      + r4_csevl ( 0.5 * y - 1.0, bk1cs, ntk1 ) ) / x;
  elseif ( x <= xmax )
    value = exp ( - x ) * r4_besk1e ( x );
  else
    value = 0.0;
  end

  return
end


text:
% Zero-crossing Frequency Estimation
%
% Estimates the instantaneous frequency of the input signal using the
% zero-crossing instantaneous frequency estimation algorithm.
%
%
% Usage:
%
%     ife = zce( signal, window_length )
%
% Parameters:
%
%
%    signal
%
%	  Input one dimensional signal to be analysed.
%
%    window_length
%
%	  Moving window to count the average number of zero crossings each
%	  time it moves.
%
%
% TFSAP 7.0
% Copyright Prof. B. Boashash
% Qatar University, Doha
% email: tfsap.research@gmail.com

author: HeLiangHIT


text:
% This function calculates the fractional derivative of order d for the
% given function r(t). It is assumed that the vector r contains the
% samples of the continuous signal r(t) which we are going to calculate its
% fractional derivative. h is a constant and represents the sampling
% period of r(t) (the time period between two samples). h must be small
% enough in the sense of Nyquist sampling theorem.
% y is the result achieved by applying the fractional differentiation
% operator on the input r. This contains the samples of the real output
% y(t) with the same sampling period used for r.
% It makes use of the Grünwald-Letnikov definition. The first element of
% the vector "r", i.e. r(1), is always zero.
%
% d        :        the order of fractional differentiation
% r        :        samples of the signal to be differentiated
% h        :        sampling poriod

function [y] = fderiv(d,r,h)

temp = 0;
for i=1:length(r)
    for j=0:i-1
        temp = temp+(-1)^j*(gamma(d+1)/(gamma(j+1)*gamma(d-j+1)))*r(i-j);
    end
    y(i) = temp;
    temp = 0;
end
y = y/(h^d);

author: Sable

repo: mcbench-benchmarks


text:
function [x,y,X,Y] = simplex(func, x, kernel, opts, varargin)
%
% SIMPLEX - multidimensional unconstrained non-linear optimsiation
%
%    X = SIMPLEX(FUNC,X) finds a local minumum of a function, via a function
%    handle FUNC, starting from an initial point X.  The local minimum is
%    located via the Nelder-Mead simplex algorithm [1], which does not require
%    any gradient information.
%
%    [X,Y] = SIMPLEX(FUNC,X) also returns the value of the function, Y, at
%    the local minimum, X.
%
%    X = SIMPLEX(FUNC,X,OPTS) allows the optimisation parameters to be
%    specified via a structure, OPTS, with members
%
%       opts.Chi         - Parameter governing expansion steps
%       opts.Delta       - Parameter governing size of initial simplex.
%       opts.Gamma       - Parameter governing contraction steps.
%       opts.Rho         - Parameter governing reflection steps.
%       opts.Sigma       - Parameter governing shrinkage steps.
%       opts.MaxIter     - Maximum number of optimisation steps.
%       opts.MaxFunEvals - Maximum number of function evaluations.
%       opts.TolFun      - Stopping criterion based on the relative change in
%                          value of the function in each step.
%       opts.TolX        - Stopping criterion based on the change in the
%                          minimiser in each step.
%
%    OPTS = SIMPLEX() returns a structure containing the default optimisation
%    parameters, with the following values:
%
%       opts.Chi         = 2
%       opts.Delta       = 0.01
%       opts.Gamma       = 0.5
%       opts.Rho         = 1
%       opts.Sigma       = 0.5
%       opts.MaxIter     = 200
%       opts.MaxFunEvals = 1000
%       opts.TolFun      = 1e-6
%       opts.TolX        = 1e-6
%
%    X = SIMPLEX(FUNC,X,OPTS, P1, P2, ...) allows addinal parameters to be
%    passed to the function to be minimised.
%
%    [X,Y,XX,YY] = SIMPLEX(FUNC, X) also returns in XX all of the values of
%    X evaluated during the optimisation process and in YY the corresponding
%    values of the function.
%
%    References:
%
%       [1] J. A. Nelder and R. Mead, "A simplex method for function
%           minimization", Computer Journal, 7:308-313, 1965.
%
% Description : Simple implementation of the Nelder-Mead simplex optimisation
%               algorithm [1].  Similar to the fminsearch routine from the
%               MATLAB optimisation toolbox.
%
% References  : [1] J. A. Nelder and R. Mead, "A simplex method for function
%                   minimization", Computer Journal, 7:308-313, 1965.
%
% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

% use default optimisation parameters if none given

if nargin < 4

    opts.Chi         = 2;
    opts.Delta       = 1.2;%0.01;
    opts.Gamma       = 0.5;
    opts.Rho         = 1;
    opts.Sigma       = 0.5;
    opts.MaxIter     = 15;%200;
    opts.MaxFunEvals = 25;
    opts.TolFun      = 1e-6;
    opts.TolX        = 1e-6;

end

% return structure containing default optimisation parameters

if nargin == 0

    x = opts;

    return

end

% get initial parameters

x = x(:);
n = length(x);
x = repmat(x', n+1, 1);
y = zeros(n+1, 1);

% form initial simplex

for i=1:n

    x(i,i) = x(i,i) + opts.Delta;
    y(i)   = func(x(i,:), varargin{:});

end

y(n+1) = func(x(n+1,:), varargin{:});
X      = x;
Y      = y;
count  = n+1;

[~, idx] = min(y);
xx        = x(idx,:);

if strcmp(kernel,'RBF_kernel') || strcmp(kernel,'wav_kernel')|| strcmp(kernel,'RBF4_kernel')
    k = 1;
    format = '  % 4d        % 4d     %e     %.4f        %.4f      %s \n';
    fprintf(1, '\n Iteration   Func-count    min f(x)    log(gamma)    log(sig2)    Procedure\n\n');
    fprintf(1, format, 1, count, min(y),xx(1),xx(2),'initial');
elseif strcmp(kernel,'lin_kernel')
    k = 2;
    format = '  % 4d        % 4d     %e     %.4f           %s \n';
    fprintf(1, '\n Iteration   Func-count    min f(x)    log(gamma)      Procedure\n\n');
    fprintf(1, format, 1, count, min(y),xx(1),'initial');
elseif strcmp(kernel,'poly_kernel')
    k = 3;
    format = '  % 4d        % 4d     %e     %.4f         %.4f         %s   \n';
    fprintf(1, '\n Iteration   Func-count    min f(x)    log(gamma)      log(t)        Procedure\n\n');
    fprintf(1, format, 1, count, min(y),xx(1),xx(2),'initial');
end

% iterative improvement

for i=2:opts.MaxIter

    % order

    [y,idx] = sort(y);
    x       = x(idx,:);

    % reflect

    centroid = mean(x(1:end-1,:));
    x_r      = centroid + opts.Rho*(centroid - x(end,:));
    y_r      = func(x_r, varargin{:});
    count    = count + 1;
    X        = [X ; x_r];
    Y        = [Y ; y_r];

    if y_r >= y(1) && y_r < y(end-1)

        % accept reflection point

        x(end,:) = x_r;
        y(end)   = y_r;
        [~, idx] = min(y);
        xx = x(idx,:);
        if k==1 || k==3
            fprintf(1, format, i, count, min(y),xx(1),xx(2),'reflect');
        elseif k==2
            fprintf(1, format, i, count, min(y),xx(1),'reflect');
        end

    else

        if y_r < y(1)

            % expand

            x_e   = centroid + opts.Chi*(x_r - centroid);
            y_e   = func(x_e, varargin{:});
            count = count + 1;
            X     = [X ; x_e];
            Y     = [Y ; y_e];

            if y_e < y_r

                % accept expansion point

                x(end,:) = x_e;
                y(end)   = y_e;
                [~, idx] = min(y);
                xx = x(idx,:);
                if k==1 || k==3
                   fprintf(1, format, i, count, min(y),xx(1),xx(2), 'expand');
                elseif k==2
                   fprintf(1, format, i, count, min(y),xx(1), 'expand');
                end

            else

                % accept reflection point

                x(end,:) = x_r;
                y(end)   = y_r;
                [~, idx] = min(y);
                xx = x(idx,:);
                if k==1 || k==3
                   fprintf(1, format, i, count, min(y),xx(1),xx(2), 'reflect');
                elseif k==2
                   fprintf(1, format, i, count, min(y),xx(1), 'reflect');
                end

            end

        else

            % contract

            shrink = 0;

            if y(end-1) <= y_r && y_r < y(end)

                % contract outside

                x_c   = centroid + opts.Gamma*(x_r - centroid);
                y_c   = func(x_c, varargin{:});
                count = count + 1;
                X     = [X ; x_c];
                Y     = [Y ; y_c];

                if y_c <= y_r

                    % accept contraction point

                    x(end,:) = x_c;
                    y(end)   = y_c;
                    [~, idx] = min(y);
                    xx = x(idx,:);
                    if k==1 || k==3
                        fprintf(1, format, i, count, min(y),xx(1),xx(2),'contract outside');
                    elseif k==2
                        fprintf(1, format, i, count, min(y),xx(1),'contract outside');
                    end

                else

                    shrink = 1;

                end

            else

                % contract inside

                x_c   = centroid + opts.Gamma*(centroid - x(end,:));
                y_c   = func(x_c, varargin{:});
                count = count + 1;
                X     = [X ; x_c];
                Y     = [Y ; y_c];

                if y_c <= y(end)

                    % accept contraction point

                    x(end,:) = x_c;
                    y(end)   = y_c;
                    [~, idx] = min(y);
                    xx = x(idx,:);
                    if k==1 || k==3
                        fprintf(1, format, i, count, min(y),xx(1),xx(2),'contract inside');
                    elseif k==2
                        fprintf(1, format, i, count, min(y),xx(1),'contract inside');
                    end

                else

                    shrink = 1;

                end

            end

            if shrink

                % shrink

                for j=2:n+1

                    x(j,:) = x(1,:) + opts.Sigma*(x(j,:) - x(1,:));
                    y(j)   = func(x(j,:), varargin{:});
                    count  = count + 1;
                    X      = [X ; x(j,:)];
                    Y      = [Y ; y(j)];

                end
                [~, idx] = min(y);
                xx = x(idx,:);
                if k==1 || k==3
                    fprintf(1, format, i, count, min(y),xx(1),xx(2),'shrink');
                elseif k==2
                    fprintf(1, format, i, count, min(y),xx(1),'shrink');
                end

            end

        end

    end

    % evaluate stopping criterion

    if max(abs(min(x) - max(x))) < opts.TolX

        fprintf(1, 'optimisation terminated sucessfully (TolX criterion)\n');

        break;

    end

    if abs(max(y) - min(y))/max(abs(y))  < opts.TolFun

        fprintf(1, 'optimisation terminated sucessfully (TolFun criterion)\n\n');

        break;

    end
    if count  > opts.MaxFunEvals

        fprintf(1, 'optimisation terminated sucessfully (MaxFunEvals criterion)\n\n');

        break;

    end

end

if i == opts.MaxIter

    fprintf(1, 'Warning : maximim number of iterations exceeded\n\n');

end

% update model structure

[y, idx] = min(y);
x        = x(idx,:);

% bye bye...

author: peterhcharlton

repo: RRest


text:
function plot_square_hhmm(ev)
% Plot the square shape implicit in the evidence.
% ev{i,t} is the value of node i in slice t.
% The observed node contains a velocity (delta increment), which is converted
% into a position.
% The Q2 node specifies which model is used, and hence which color
% to use: 1=red, 2=green, 3=blue, 4=black.

Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;

delta = cell2num(ev(Onode,:)); % delta(:,t)
Q2label = cell2num(ev(Q2,:));

T = size(delta, 2);
pos = zeros(2,T+1);
hold on
cols = {'r', 'g', 'b', 'k'};
for t=2:T+1
  pos(:,t) = pos(:,t-1) + delta(:,t-1);
  plot(pos(1,t), pos(2,t), sprintf('%c.', cols{Q2label(t-1)}));
  if (t==2)
    text(pos(1,t-1),pos(2,t-1),sprintf('%d',t))
  elseif (mod(t,20)==0)
    text(pos(1,t),pos(2,t),sprintf('%d',t))
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [interp_value] = grid_bilin_interp(X_approx, Y_approx, grid, ncols, nrows, cellsize, Xll, Yll, nodata)

% SYNTAX:
%   [interp_value] = grid_bilin_interp(X_approx, Y_approx, grid, ncols, nrows, cellsize, Xll, Yll, nodata);
%
% INPUT:
%   X_approx = X coordinate of the interpolation point
%   Y_approx = Y coordinate of the interpolation point
%   grid = matrix containing the grid
%   ncols = number of columns of the grid
%   nrows = number of rows of the grid
%   cellsize = ground size of a cell
%   Xll = X coordinate of the center of the lower left cell
%   Yll = Y coordinate of the center of the lower left cell
%   nodata = value used for cells not containing data
%
% OUTPUT:
%   interp_value = interpolated value
%
% DESCRIPTION:
%   Function that applies a bilinear interpolation of the four nearest nodes
%   of a georeferenced grid in correspondence of a point of given coordinates.

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Mirko Reguzzoni
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

%preparation of the grid axes
X = (Xll : cellsize : Xll + (ncols - 1) * cellsize)';
Y = (Yll : cellsize : Yll + (nrows - 1) * cellsize)';

if (X_approx <= X(1) | X_approx >= X(end) | Y_approx <= Y(1) | Y_approx >= Y(end))
    interp_value = nodata;
    return
end

%detection of the grid node nearest to the interpolation point
[mX, posX] = min(abs(X - X_approx));
[mY, posY] = min(abs(Y - Y_approx));

%definition of the four grid nodes that sorround the interpolation point
% (i,j) image coordinates (upper-left origin)
% (X,Y) ground coordinates (bottom-left origin)
if (X(posX) > X_approx) | (mX ==0)
    j_left = posX - 1;
    j_right = posX;
    X_left = X(posX - 1);
else
    j_left = posX;
    j_right = posX + 1;
    X_left = X(posX);
end

if (Y(posY) > Y_approx) | (mY ==0)
    i_up = nrows + 1 - posY;
    i_down = i_up + 1;
    Y_down = Y(posY - 1);
else
    i_down = nrows + 1 - posY;
    i_up = i_down - 1;
    Y_down = Y(posY);
end

%if one of the interp_value values of the four sorrounding points is a nodata value, do not interpolate and return nodata value
if (grid(i_up,j_left) == nodata | grid(i_up,j_right) == nodata | grid(i_down,j_left) == nodata | grid(i_down,j_right) == nodata)
    interp_value = nodata;
    return
end

%computation of the parameters of the bilinear function
%f(X, Y) = a*X*Y + b*X + c*Y + d

A = [0 0 cellsize 1; ...
     cellsize^2 cellsize cellsize 1; ...
     0 0 0 1; ...
     0 cellsize 0 1];

B = [grid(i_up,j_left);...
     grid(i_up,j_right);...
     grid(i_down,j_left);...
     grid(i_down,j_right)];

bilin_param = A\B;

i_approx = Y_approx - Y_down;
j_approx = X_approx - X_left;

%computation of the interpolated value
interp_value = bilin_param(1) * j_approx * i_approx + bilin_param(2) * j_approx + bilin_param(3) * i_approx + bilin_param(4);
interp_value = double(interp_value);

author: goGPS-Project


text:
function J = computeCost(X, y, theta)
%COMPUTECOST Compute cost for linear regression
%   J = COMPUTECOST(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.
J = (1/(2*m))*sum(((X*theta) - y).^2)

% =========================================================================

end

author: anirudhjayaraman

repo: Machine-Learning


text:
[K u] = sphere_wp;
h=trisurf(K,u(:,1),u(:,2),u(:,3));
set(h,'FaceColor',[0.7 0.7 0.7],'EdgeColor','none','FaceLighting','gouraud');
axis equal;
light('Position',[1 -1 0],'Style','infinite');
% light('Position',[10 10 10],'Style','infinite');
xlim([-2,10]);
ylim([-2,2]);
zlim([-2,2]);

% animation loop:
u0=u;
for t=0:0.02:5
    R=1+0.2*sin(10*t); % change radius
    x=0+4+3*sin(8.1234*t); % change x position
    u=R*u0;
    u(:,1)=u(:,1)+x;
    set(h,'Vertices',u);
    drawnow;
    pause(1/50);
end


text:
function status = getKaleido()
status=0;
kDir = fullfile(fileparts(mfilename('fullpath')),'..','kaleido');
arch = computer('arch');
plotlyJS = 'https://cdn.plot.ly/plotly-latest.min.js';

if isunix()
    if ismac()
        if strcmp(arch,'maci64')
            url = 'https://github.com/plotly/Kaleido/releases/download/v0.2.1/kaleido_mac_x64.zip';
        else
            url = 'https://github.com/plotly/Kaleido/releases/download/v0.2.1/kaleido_mac_arm64.zip';
        end
    else
        if strcmp(arch,'glnxa64')
            url = 'https://github.com/plotly/Kaleido/releases/download/v0.2.1/kaleido_linux_x64.zip';
        else
            url = 'https://github.com/plotly/Kaleido/releases/download/v0.2.1/kaleido_linux_arm64.zip';
        end
    end
elseif ispc()
    if strcmp(arch,'win64')
        url = 'https://github.com/plotly/Kaleido/releases/download/v0.2.1/kaleido_win_x64.zip';
    else
        url = 'https://github.com/plotly/Kaleido/releases/download/v0.2.1/kaleido_win_x86.zip';
    end
else
    fprintf('\nCouldn''t find a suitable version of "Kaleido" required for exporting plots, for your system\n');
    fprintf('https://github.com/plotly/Kaleido/tags\n');
    fprintf('Please download a suitable Kaleido version for your system, and unzip the files in plotly/Kaleido folder.\n\n');
    return;
end

fprintf('\nTrying to download Kaleido executable and plotly javascript (one-time download)... please wait...');
try
    file = websave(fullfile(kDir,'kaleido.zip'),url);
    websave(fullfile(kDir,'plotly-latest.min.js'),plotlyJS);
catch
    fprintf('\nFailed to download Kaleido.\n');
    fprintf('https://github.com/plotly/Kaleido/tags\n');
    fprintf('Please download a suitable Kaleido version for your system, and unzip the files in plotly/Kaleido folder.\n\n');
    return;
end
unzip(file,kDir);
status=1;
fprintf(' Done\n\n');

author: plotly


text:
function rMTAsaveInExcel(xls_filename, TSscore, deletedGenes, alpha_values, varargin)
% Save in Excel all the results calculated using rMTA.
%
% USAGE:
%
%    rMTAsaveInExcel(output_file, TSscore, deletedGenes, alpha_values, varargin)
%
% INPUTS:
%    xls_filename:       name of the resulting excel file
%    TSscore:            Transformation score by each transformation
%    deletedGenes:       The list of genes/reactions removed in each knock-out
%    alpha_values:       Numeric value or array. Parameter of the quadratic
%                        problem (default = alpha_1, alpha_2,...)
%
% OPTIONAL INPUTS:
%    gene_info:          Table with the gene names and additional
%                        information. It needs a column called 'gene',
%                        which contains the same information as the
%                        deletedGenes variable.
%    differ_genes:       Table with the gene differential expression
%                        information. It needs a column called 'gene',
%                        which contains the same information as the
%                        deletedGenes variable.
%    RankingGeneID:      Caracter array with the ID of the gene which score
%                        is desired to know.
%
% .. Authors:
%       - Luis V. Valcarcel, 08/03/2021, University of Navarra, CIMA & TECNUN School of Engineering.

p = inputParser;
% check required arguments
addRequired(p, 'xls_filename', @(x)ischar(x));
addRequired(p, 'TSscore');
addRequired(p, 'deletedGenes');
% Check optional arguments
addOptional(p, 'alpha_values', -1, @isnumeric);
% Add optional name-value pair argument
addParameter(p, 'gene_info', []);
addParameter(p, 'differ_genes', []);
addParameter(p, 'RankingGeneID', [], @(x)ischar(x));
% extract variables from parser
parse(p, xls_filename, TSscore, deletedGenes, alpha_values, varargin{:});
gene_info = p.Results.gene_info;
differ_genes = p.Results.differ_genes;
RankingGeneID = p.Results.RankingGeneID;
num_alphas = numel(alpha_values);
if alpha_values == -1
    alpha_values = 1:size(TSscore.bTS,2);
end

% In order to store the results of all the cases, a cell array is going to
% store a table with the results for each alpha
results = cell(size(alpha_values));

% threshold for viable results
th = -1e20; % not (-INF) == infeasible

for i = 1:num_alphas
    fprintf('\tProcess results for alpha = %1.2f \n',alpha_values(i));
    T = table(deletedGenes, TSscore.bTS(:,i), TSscore.mTS, TSscore.wTS(:,i),TSscore.rTS(:,i));
    T.Properties.VariableNames = {'gene','bTS','mTS','wTS','rTS'};

    if ~isempty(gene_info)
        % join to extra information of genes
        T = join(gene_info,T,'Keys','gene');
    end

    %sort by MOMA results
    T = sortrows(T,'rTS','descend');

    % store only viable genes
    T = T(T.bTS>th,:);
    results{i} = T;
end

fprintf('\tSelected the "%u" best solutions\n',sum(T.mTS>(th)));

%export to excel
LETTERS = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ';
% Export results for the different alphas
for i = 1:num_alphas
    T = results{i};
    fprintf('\tWrite xls for alpha = %1.2f \n',alpha_values(i));
    sheet_name = ['alpha = ',num2str(alpha_values(i))];

    % write results
    writetable(T,xls_filename,'FileType','spreadsheet','Sheet',sheet_name,'WriteVariableNames',true);

    LETTER = LETTERS(size(T,2)+3);

    if exist('RankingGeneID','var') && ~isempty(RankingGeneID)
        % ranking of gene
        xlswrite(xls_filename,{['position ' RankingGeneID]},sheet_name,[LETTER '2']);
        xlswrite(xls_filename,{['=MATCH(' RankingGeneID ',A:A,0)-1']},sheet_name,[LETTER '3']);
        xlswrite(xls_filename,{'ranking genes with KO'},sheet_name,[LETTER '4']);
        xlswrite(xls_filename,{['=' LETTER '3/' num2str(size(T,1))]},sheet_name,[LETTER '5']);
        xlswrite(xls_filename,{'ranking genes model'},sheet_name, [LETTER '6']);
        xlswrite(xls_filename,{['=' LETTER '3/' num2str(numel(deletedGenes))]}, sheet_name,[LETTER '7']);
    end
end

fprintf('\tWrite selection of best genes\n');
num = numel(results{i}.gene);
TopGenes_ID = cell(num,num_alphas);
titles = cell(1,num_alphas);
for i = 1:num_alphas
    num = sum(results{i}.rTS > 0);
    TopGenes_ID(1:num,i) = results{i}.gene(1:num);
    titles{i} = ['alpha = ',num2str(alpha_values(i))];
end
xlswrite(xls_filename,titles,'TOP_GENES_ID','A1');
xlswrite(xls_filename,TopGenes_ID,'TOP_GENES_ID','A2');

TopGenes_join = TopGenes_ID(~cellfun(@isempty,TopGenes_ID(:,1)),1);
for i = 1:num_alphas
    TopGenes_join = intersect(TopGenes_join,TopGenes_ID(~cellfun(@isempty,TopGenes_ID(:,i)),i));
end

if ~isempty(differ_genes)
    differ_genes = differ_genes(~strcmp(differ_genes.gene,''),:);
    [~, idx, ~] = unique(differ_genes.gene);
    differ_genes = differ_genes(idx,:);
    differ_genes.Properties.RowNames = differ_genes.gene;

    idx = intersect(deletedGenes,differ_genes.gene);

    xlswrite(xls_filename,differ_genes.Properties.VariableNames,'TOP_GENES_diff_expr','A1');
    xlswrite(xls_filename,table2cell(differ_genes(idx,:)),'TOP_GENES_diff_expr','A2');
end

if ~isempty(differ_genes) &&  ~isempty(RankingGeneID)

    fprintf('\tSUMMARY\n');

    differ_genes_aux = differ_genes(differ_genes.logFC>0,:);
    genes_FC_pos = differ_genes_aux.Properties.RowNames;

    SUMMARY = cell(7,2+5*num_alphas);
    table_aux_ranking = array2table(zeros(4,num_alphas));
    table_aux_ranking.Properties.RowNames = {'bTS','mTS','wTS','rTS'};
    table_aux_ranking.Properties.VariableNames = regexprep(strcat('alpha_',cellfun(@num2str,num2cell(alpha_values'),'UniformOutput',0)),'0.','0_');
    table_aux_ranking_FC_pos = table_aux_ranking;
    table_aux_sign = array2table(char(zeros(4,num_alphas)));
    table_aux_sign.Properties.RowNames = {'bTS','mTS','wTS','rTS'};
    table_aux_sign.Properties.VariableNames = regexprep(strcat('alpha_',cellfun(@num2str,num2cell(alpha_values'),'UniformOutput',0)),'0.','0_');
    num_KO = zeros(1,num_alphas);
    num_KO_FC_pos = zeros(1,num_alphas);
    for i = 1:num_alphas
        T = results{i};
        T.Properties.RowNames=T.gene;
        T2 = T(intersect(T.Properties.RowNames, genes_FC_pos),:);
        for j=1:size(table_aux_ranking,1)
            k = table_aux_ranking.Properties.RowNames{j};
            T = sortrows(T,k,'descend');
            T2 = sortrows(T2,k,'descend');
            try RankingGeneID = num2str(RankingGeneID); end
            table_aux_ranking{k,i} = find(strcmp(T.gene,RankingGeneID));
            try table_aux_ranking_FC_pos{k,i} = find(strcmp(T2.gene,RankingGeneID)); end
        end
        T = sortrows(T,'wTS','ascend');
        T2 = sortrows(T2,'wTS','ascend');
        table_aux_ranking{'wTS',i} = find(strcmp(T.gene,RankingGeneID));
        try table_aux_ranking_FC_pos{'wTS',i} = find(strcmp(T2.gene,RankingGeneID)); end
        aux = sign(T{RankingGeneID,table_aux_sign.Properties.RowNames})';
        aux2 = char(aux);
        aux2(:) = '+';
        aux2(aux<0) = '-';
        table_aux_sign{:,i} = aux2;
        num_KO(i) = size(T,1);
        num_KO_FC_pos(i) = size(T2,1);
    end

    table_aux_ranking_numKO = table2array(table_aux_ranking) ./ num_KO;
    table_aux_ranking_numDE = table2array(table_aux_ranking_FC_pos) ./ num_KO_FC_pos;
    table_aux_ranking_numgenes = table2array(table_aux_ranking) ./ length(deletedGenes);

    SUMMARY{1,2+0*num_alphas+1} = 'ranking';
    SUMMARY{1,2+1*num_alphas+1} = 'sign';
    SUMMARY{1,2+2*num_alphas+1} = 'ranking genes with KO';
    SUMMARY{1,2+3*num_alphas+1} = 'ranking genes with DE FC>0';
    SUMMARY{1,2+4*num_alphas+1} = 'ranking genes Recon1';

    SUMMARY(2,:) = [{'Case'} {'Score'} repmat(table_aux_ranking.Properties.VariableNames,1,5)];
    SUMMARY(3:6,2) = table_aux_ranking.Properties.RowNames;
    SUMMARY(3:6,3:end) = [table2cell(table_aux_ranking), table2cell(table_aux_sign), num2cell(table_aux_ranking_numKO), num2cell(table_aux_ranking_numDE), num2cell(table_aux_ranking_numgenes)];
    SUMMARY(end,2+2*num_alphas+(1:num_alphas)) = num2cell(num_KO);
    SUMMARY(end,2+3*num_alphas+(1:num_alphas)) = num2cell(num_KO_FC_pos);
    SUMMARY(end,2+4*num_alphas+(1:num_alphas)) = num2cell(length(deletedGenes));
    SUMMARY{1,1} = date();
    SUMMARY{end,1} = xls_filename;

    xlswrite(xls_filename,SUMMARY,'SUMMARY');
end

fprintf('DONE\n');

end

author: opencobra

repo: cobratoolbox