text:
function rez = learnTemplates2(rez, iorder)
% This is the main optimization. Takes the longest time and uses the GPU heavily.

rez.ops.fig = getOr(rez.ops, 'fig', 1); % whether to show plots every N batches

% Turn on sorting of spikes before subtracting and averaging in mpnu8
rez.ops.useStableMode = getOr(rez.ops, 'useStableMode', 1);
useStableMode = rez.ops.useStableMode;

NrankPC = 6; % this one is the rank of the PCs, used to detect spikes with threshold crossings
Nrank = 3; % this one is the rank of the templates

rez.ops.LTseed = getOr(rez.ops, 'LTseed', 1);
rng('default'); rng(rez.ops.LTseed);

ops = rez.ops;

% we need PC waveforms, as well as template waveforms
wTEMP = rez.wTEMP;
wPCA  = rez.wPCA;

% move these to the GPU
wPCA = gpuArray(wPCA(:, 1:Nrank));
wTEMP = gpuArray(wTEMP);
wPCAd = double(wPCA); % convert to double for extra precision
ops.wPCA = gather(wPCA);
ops.wTEMP = gather(wTEMP);
nt0 = ops.nt0;
nt0min  = rez.ops.nt0min;
rez.ops = ops;
nBatches  = rez.temp.Nbatch;
NT  	= ops.NT;
Nfilt 	= ops.Nfilt;
Nchan 	= ops.Nchan;

% two variables for the same thing? number of nearest channels to each primary channel
NchanNear   = min(ops.Nchan, 32);
Nnearest    = min(ops.Nchan, 32);

% decay of gaussian spatial mask centered on a channel
sigmaMask  = ops.sigmaMask;

% spike threshold for finding missed spikes in residuals
ops.spkTh = -6; % why am I overwriting this here?

batchstart = 0:NT:NT*nBatches;

% find the closest NchanNear channels, and the masks for those channels
[iC, mask, C2C] = getClosestChannels(rez, sigmaMask, NchanNear);

niter   = numel(iorder);

% this is the absolute temporal offset in seconds corresponding to the start of the
% spike sorted time segment
t0 = ceil(rez.ops.trange(1) * ops.fs);

nInnerIter  = 60; % this is for SVD for the power iteration

% schedule of learning rates for the model fitting part
% starts small and goes high, it corresponds approximately to the number of spikes
% from the past that were averaged to give rise to the current template
pmi = exp(-1./linspace(ops.momentum(2), ops.momentum(2), niter));

Nsum = min(Nchan,7); % how many channels to extend out the waveform in mexgetspikes
% lots of parameters passed into the CUDA scripts

% initialize the list of channels each template lives on
iList = int32(gpuArray(zeros(Nnearest, Nfilt)));

dWU = gpuArray(rez.dWU);
W = gpuArray(rez.W);
Nfilt = size(W,2); % update the number of filters/templates
% initialize average number of spikes per batch for each template
nsp = gpuArray.zeros(Nfilt,1, 'double');

Params     = double([NT Nfilt ops.Th(1) nInnerIter nt0 Nnearest ...
    Nrank ops.lam pmi(1) Nchan NchanNear ops.nt0min 1 Nsum NrankPC ops.Th(1) useStableMode]);
Params(2) = Nfilt; % update in the CUDA parameters
% this flag starts 0, is set to 1 later
Params(13) = 0;

% kernels for subsample alignment
[Ka, Kb] = getKernels(ops, 10, 1);

p1 = .95; % decay of nsp estimate in each batch

fprintf('Time %3.0fs. Optimizing templates ...\n', toc)

fid = fopen(ops.fproc, 'r');

%%

for ibatch = 1:niter
    k = iorder(ibatch); % k is the index of the batch in absolute terms

    % obtained pm for this batch
    Params(9) = pmi(ibatch);
    pm = pmi(ibatch);

    % loading a single batch (same as everywhere)
    offset = 2 * ops.Nchan*batchstart(k);
    fseek(fid, offset, 'bof');
    dat = fread(fid, [ops.Nchan NT+ops.ntbuff], '*int16');
    dat = dat';
    dataRAW = single(gpuArray(dat))/ ops.scaleproc;
    Params(1) = size(dataRAW,1);

    % resort the order of the templates according to best peak channel
    % this is important in order to have cohesive memory requests from the GPU RAM
    [~, iW] = max(abs(dWU(nt0min, :, :)), [], 2); % max channel (either positive or negative peak)
    iW = int32(squeeze(iW));

    % decompose dWU by svd of time and space (via covariance matrix of 61 by 61 samples)
    % this uses a "warm start" by remembering the W from the previous iteration
    [W, U, mu] = mexSVDsmall2(Params, dWU, W, iC-1, iW-1, Ka, Kb);

    % UtU is the gram matrix of the spatial components of the low-rank SVDs
    % it tells us which pairs of templates are likely to "interfere" with each other
    % such as when we subtract off a template
    [UtU, maskU] = getMeUtU(iW, iC, mask, Nnearest, Nchan); % this needs to change (but I don't know why!)

    % main CUDA function in the whole codebase. does the iterative template matching
    % based on the current templates, gets features for these templates if requested (featW, featPC),
    % gets scores for the template fits to each spike (vexp), outputs the average of
    % waveforms assigned to each cluster (dWU0),
    % and probably a few more things I forget about
    [st0, id0, x0, featW, dWU0, drez, nsp0, featPC, vexp, errmsg] = ...
        mexMPnu8(Params, dataRAW, single(U), single(W), single(mu), iC-1, iW-1, UtU, iList-1, ...
        wPCA);

    % errmsg returns 1 if caller requested "stableMode" but mexMPnu8 was
    % compiled without the sorter enabled (i.e. STABLEMODE_ENABLE = false
    % in mexGPUAll). Send an error message to the console just once if this
    % is the case:
    if (ibatch == 1)
        if( (useStableMode == 1) && (errmsg == 1) )
            fprintf( 'useStableMode selected but STABLEMODE not enabled in compiled mexMPnu8.\n' );
        end
    end

    % Sometimes nsp can get transposed (think this has to do with it being
    % a single element in one iteration, to which elements are added
    % nsp, nsp0, and pm must all be row vectors (Nfilt x 1), so force nsp
    % to be a row vector.
    [nsprow, nspcol] = size(nsp);
    if nsprow<nspcol
        nsp = nsp';
    end

    % updates the templates as a running average weighted by recency
    % since some clusters have different number of spikes, we need to apply the
    % exp(pm) factor several times, and fexp is the resulting update factor
    % for each template
    fexp = exp(double(nsp0).*log(pm));
    fexp = reshape(fexp, 1,1,[]);
    dWU = dWU .* fexp + (1-fexp) .* (dWU0./reshape(max(1, double(nsp0)), 1,1, []));

    % nsp just gets updated according to the fixed factor p1
    nsp = nsp * p1 + (1-p1) * double(nsp0);
    % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
    if (rem(ibatch, 100)==1)
        % this is some of the relevant diagnostic information to be printed during training
        fprintf('%2.2f sec, %d / %d batches, %d units, nspks: %2.4f, mu: %2.4f, nst0: %d \n', ...
            toc, ibatch, niter, Nfilt, sum(nsp), median(mu), numel(st0))

        % these diagnostic figures should be mostly self-explanatory
        if ops.fig
            try
                if ibatch==1
                    figHand = figure;
                else
                    figure(figHand);
                end
                make_fig(W, U, mu, nsp)
            catch ME
               warning('Error making figure was: %s',ME.message);
            end
        end
    end
end
fclose(fid);
toc

%%
% final covariance matrix between all templates
[WtW, iList] = getMeWtW(single(W), single(U), Nnearest);
% iW is the final channel assigned to each template
[~, iW] = max(abs(dWU(nt0min, :, :)), [], 2);
iW = int32(squeeze(iW));
% the similarity score between templates is simply the correlation,
% taken as the max over several consecutive time delays
rez.simScore = gather(max(WtW, [], 3));

rez.iNeighPC    = gather(iC(:, iW));

% the neihboring templates idnices are stored in iNeigh
rez.iNeigh   = gather(iList);

rez = memorizeW(rez, W, dWU, U, mu); % memorize the state of the templates
rez.ops = ops; % update these (only rez comes out of this script)
rez.nsp = nsp;

% save('rez_mid.mat', 'rez');

fprintf('Finished learning templates \n')
%%

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function [blastStructure,diamondReport]=getDiamond(organismID,fastaFile,...
    modelIDs,refFastaFiles,developMode,hideVerbose)
% getDiamond
%   Uses DIAMOND to perform a bidirectional BLAST between the organism
%   of interest and a set of template organisms
%
%   Input:
%   organismID      the id of the organism of interest. This should also
%                   match with the id supplied to getModelFromHomology
%   fastaFile       a FASTA file with the protein sequences for the
%                   organism of interest
%   modelIDs        a cell array of model ids. These must match the
%                   "model.id" fields in the "models" structure if the
%                   output is to be used with getModelFromHomology
%   refFastaFiles   a cell array with the paths to the corresponding FASTA
%                   files
%   developMode     true if blastReport should be generated that is used
%                   in the unit testing function for DIAMOND (opt, default
%                   false)
%   hideVerbose     true if no status messages should be printed (opt,
%                   default false)
%
%   Output:
%   blastStructure  structure containing the bidirectional homology
%                   measurements which are used by getModelFromHomology
%   diamondReport   structure containing MD5 hashes for FASTA database
%                   files and non-parsed BLAST output data. Will be blank
%                   if developMode is false.
%
%   NOTE: This function calls DIAMOND to perform a bidirectional homology
%   search between the organism of interest and a set of other organisms
%   using the '--more-sensitive' setting from DIAMOND. For the most
%   sensitive results, the use of getBlast() is adviced, however,
%   getDiamond() is a fast alternative (>15x faster). The blastStructure
%   generated is in the same format as those obtained from getBlast().
%
%   Usage: [blastStructure,diamondReport]=getDiamond(organismID,fastaFile,...
%    modelIDs,refFastaFiles,developMode,hideVerbose)

if nargin<5
    developMode = false;
end
if nargin<6
    hideVerbose = false;
end

%Everything should be cell arrays
organismID=convertCharArray(organismID);
fastaFile=convertCharArray(fastaFile);
modelIDs=convertCharArray(modelIDs);
refFastaFiles=convertCharArray(refFastaFiles);

%Create blank structures for results
blastStructure=[];
diamondReport.dbHashes={};
diamondReport.diamondTxtOutput={};

%Get the directory for RAVEN Toolbox.
ravenPath=findRAVENroot();

%Generate temporary names for DIAMOND databases and output files
tmpDB=tempname;
outFile=tempname;

%Check for existence of files. If no full path is specified for a file,
%assume that it is in the current folder
if isrow(refFastaFiles)
    files=horzcat(fastaFile,refFastaFiles);
else
    files=vertcat(fastaFile,refFastaFiles);
end

files=checkFileExistence(files,2); %Copy files to temp dir
fastaFile = files(1);
refFastaFiles = files(2:end);

%Identify the operating system
if isunix
    if ismac
        binEnd='.mac';
    else
        binEnd='';
    end
elseif ispc
    binEnd='.exe';
else
    dispEM('Unknown OS, exiting.')
    return
end

%Run DIAMOND multi-threaded to use all logical cores assigned to MATLAB.
cores = evalc('feature(''numcores'')');
cores = strsplit(cores, 'MATLAB was assigned: ');
cores = regexp(cores{2},'^\d*','match');
cores = cores{1};

%Create a database for the new organism and blast each of the refFastaFiles
%against it
[status, message]=system(['"' fullfile(ravenPath,'software','diamond',['diamond' binEnd]) '" makedb --in "' fastaFile{1} '" --db "' fullfile(tmpDB) '"']);
if developMode
    diamondReport.dbHashes{numel(diamondReport.dbHashes)+1} = char(regexp(message,'[a-f0-9]{32}','match'));
end
if status~=0
    EM=['DIAMOND makedb did not run successfully, error: ', num2str(status)];
    dispEM(EM,true);
end

for i=1:numel(refFastaFiles)
    if ~hideVerbose
        fprintf(['Running DIAMOND blastp with "' modelIDs{i} '" against "' organismID{1} '"..\n']);
    end
    [status, ~]=system(['"' fullfile(ravenPath,'software','diamond',['diamond' binEnd]) '" blastp --query "' refFastaFiles{i} '" --out "' outFile '_' num2str(i) '" --db "' fullfile(tmpDB) '" --more-sensitive --outfmt 6 qseqid sseqid evalue pident length bitscore ppos --threads ' cores ]);
    if developMode
        diamondReport.diamondTxtOutput{numel(diamondReport.diamondTxtOutput)+1}=importdata([outFile '_' num2str(i)]);
    end
    if status~=0
        EM=['DIAMOND blastp did not run successfully, error: ', num2str(status)];
        dispEM(EM,true);
    end
end
delete([tmpDB filesep 'tmpDB*']);

%Then create a database for each of the reference organisms and blast the
%new organism against them
for i=1:numel(refFastaFiles)
    if ~hideVerbose
        fprintf(['Running DIAMOND blastp with "' organismID{1} '" against "' modelIDs{i} '"..\n']);
    end
    [status, message]=system(['"' fullfile(ravenPath,'software','diamond',['diamond' binEnd]) '" makedb --in "' refFastaFiles{i} '" --db "' fullfile(tmpDB) '"']);
    if status~=0
        EM=['DIAMOND makedb did not run successfully, error: ', num2str(status)];
        dispEM(EM,true);
    end
    [status, ~]=system(['"' fullfile(ravenPath,'software','diamond',['diamond' binEnd]) '" blastp --query "' fastaFile{1} '" --out "' outFile '_r' num2str(i) '" --db "' fullfile(tmpDB) '" --more-sensitive --outfmt 6 qseqid sseqid evalue pident length bitscore ppos --threads ' cores]);
    if developMode
        diamondReport.dbHashes{numel(diamondReport.dbHashes)+1} = char(regexp(message,'[a-f0-9]{32}','match'));
        diamondReport.diamondTxtOutput{numel(diamondReport.diamondTxtOutput)+1}=importdata([outFile '_r' num2str(i)]);
    end
    if status~=0
        EM=['DIAMOND blastp did not run successfully, error: ', num2str(status)];
        dispEM(EM,true);
    end
    delete([tmpDB filesep 'tmpDB*']);
end

%Done with the DIAMOND blastp, do the parsing of the text files
for i=1:numel(refFastaFiles)*2
    tempStruct=[];
    if i<=numel(refFastaFiles)
        tempStruct.fromId=modelIDs{i};
        tempStruct.toId=organismID{1};
        A=readtable([outFile '_' num2str(i)],'Delimiter','\t','Format','%s%s%f%f%f%f%f');
    else
        tempStruct.fromId=organismID{1};
        tempStruct.toId=modelIDs{i-numel(refFastaFiles)};
        A=readtable([outFile '_r' num2str(i-numel(refFastaFiles))],'Delimiter','\t','Format','%s%s%f%f%f%f%f');
    end
    tempStruct.fromGenes=A{:,1};
    tempStruct.toGenes=A{:,2};
    tempStruct.evalue=table2array(A(:,3));
    tempStruct.identity=table2array(A(:,4));
    tempStruct.aligLen=table2array(A(:,5));
    tempStruct.bitscore=table2array(A(:,6));
    tempStruct.ppos=table2array(A(:,7));
    blastStructure=[blastStructure tempStruct];
end

%Remove the old tempfiles
delete([outFile '*']);
%Remove the temp fasta files
delete(files{:})
end

author: SysBioChalmers

repo: RAVEN


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

m = size(X,1)

for i = 1:m
    distance_array = zeros(1,K);
    for j = 1:K
        distance_array(1,j) = sqrt(sum(power((X(i,:)-centroids(j,:)),2)));
    end
    [d, d_idx] = min(distance_array);
    idx(i,1) = d_idx;
end

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

end

author: Borye


text:
classdef DensityPlotterForPerimeter < handle

    properties (Access = private)
        inputFile
        mesh
        scale
        density
        densityVariable
        plotter
    end

    methods (Access = public)

        function obj = DensityPlotterForPerimeter(cParams)
            obj.init(cParams)
        end

        function plot(obj)
            obj.createDesignVariable();
            obj.createPlotter();
            obj.plotter.refresh();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh       = cParams.mesh;
            obj.inputFile  = cParams.inputFile;
            obj.scale      = cParams.scale;
            obj.density    = cParams.density;
        end

        function createDesignVariable(obj)
            s.mesh      = obj.mesh;
            s.inputFile = obj.inputFile;
            s.scale     = obj.scale;
            s.type      = 'Density';
            d = DesignVariableCreatorSettings(s);
            s = d.create();
            obj.densityVariable       = DesignVariable.create(s);
            obj.densityVariable.update(obj.density);
            obj.densityVariable.rho   = obj.density;
        end

        function createPlotter(obj)
            sD.shallDisplay   = true;
            sD.showBC         = false;
            sD.bc             = false;
            sD.designVariable = obj.densityVariable;
            sD.optimizerNames = '';
            sD.dim            = '2D';
            sD.scale          = obj.scale;
            sD.mesh           = Mesh_Total(obj.mesh);
            f = DesignVarMonitorFactory;
            obj.plotter = f.create(sD);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function SO3F = smooth(SO3F,varargin)
% smooth SO3FunComposition
%
% Input
%  SO3F - @SO3FunComposition
%  psi - @SO3Kernel (smoothing kernel)
%
% Options
%  halfwidth
%
% Output
%  SO3F - smoothed @SO3Fun
%
% See also
% SO3Fun/smooth

% smooth components
for i = 1:length(SO3F.components)
  SO3F.components{i} = SO3F.components{i}.smooth(varargin{:});
end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function out = draw_arrow(startpoint,endpoint,headsize)
%by Ryan Molecke
% accepts two [x y] coords and one double headsize

   v1 = headsize*(startpoint-endpoint)/2.5;

   theta = 22.5*pi/180;
   theta1 = -1*22.5*pi/180;
   rotMatrix = [cos(theta)  -sin(theta) ; sin(theta)  cos(theta)];
   rotMatrix1 = [cos(theta1)  -sin(theta1) ; sin(theta1)  cos(theta1)];

   v2 = v1*rotMatrix;
   v3 = v1*rotMatrix1;
   x1 = endpoint;
   x2 = x1 + v2;
   x3 = x1 + v3;
hold on;
   fill([x1(1) x2(1) x3(1)],[x1(2) x2(2) x3(2)],[0 0 0]);     % this fills the arrowhead (black)
   plot([startpoint(1) endpoint(1)],[startpoint(2) endpoint(2)],'linewidth',2,'color',[0 0 0]);

author: Sable

repo: mcbench-benchmarks


text:
function pattern = elementaryCellularAutomata(rule, n, width, randfrac)
%elementaryCellularAutomata Elementary 1D cellular automaton patterns
%
%   PATTERN = elementaryCellularAutomata(RULE, NITER), where NITER is a
%   scalar, returns an NITER x 2*NITER+1 matrix whose entries are all 0 or
%   1. The I'th row of the matrix contains the state of the elementary 1D
%   cellular automaton at iteration I-1, counting the initial state as
%   iteration 0. The integer RULE specifies the rule to use as set out at
%   http://mathworld.wolfram.com/ElementaryCellularAutomaton.html.
%
%   The initial state is all zeros except for a 1 at element NITER+1 (the
%   central element) of the CA.
%
%   PATTERN = elementaryCellularAutomata(RULE, NITER, WID), where WID is a
%   scalar, returns an NITER x WID matrix. The array of cells is taken to
%   be circular, so that PATTERN(I+1,I) depends on PATTERN{(I,WID), (I,1)
%   and (I,2)}. Similarly, PATTERN(I+1,WID) depends on PATTERN{(I,WID-1),
%   (I,WID) and (I,1)}. This only matters if WID < 2*NITER+1 and RULE is
%   such that the pattern propagates outwards, so reaching the boundaries
%   of the array.
%
%   This wraparound allows long, thin patterns to be generated if an
%   appropriate rule is chosen. All patterns with a fixed width will be
%   periodic, unless some random noise is added.
%
%   The state on the first iteration is all zeros except for a 1 at element
%   floor((WID+1)/2) of the CA.
%
%   PATTERN = elementaryCellularAutomata(RULE, NITER, START) where START is
%   a 1 x WID row vector containing only the values 0 and 1, is as above
%   except that the initial state is given by the entries in START. Thus on
%   exit, PATTERN(1,:) is equal to START.
%
%   PATTERN = elementaryCellularAutomata(RULE, NITER, WIDSTART, FNOISE) is
%   as above except that noise is added to the process. WIDSTART can be
%   either a scalar giving the width or a vector giving the start state; an
%   empty matrix is equivalent to 2*NITER+1. FNOISE is a number from 0 to 1
%   giving the probability that any given cell will be set to the wrong
%   state (the complement of the state given by the rule) on any one
%   iteration.
%
%   Example
%   -------
%       % show 50 rows of each pattern
%       for rule = 0:255
%           pattern = elementaryCellularAutomata(rule, 50);
%           imshow(pattern); pause;
%       end

%   Copyright 2010 David Young

% check arguments and supply defaults
error(nargchk(2, 4, nargin));
validateattributes(rule, {'numeric'}, {'scalar' 'integer' 'nonnegative' '<=' 255}, ...
    'elementaryCellularAutomata', 'RULE');
validateattributes(n, {'numeric'}, {'scalar' 'integer' 'positive'}, ...
    'elementaryCellularAutomata', 'N');
if nargin < 3 || isempty(width)
    width = 2*n-1;
elseif isscalar(width)
    validateattributes(width, {'numeric'}, {'integer' 'positive'}, ...
        'elementaryCellularAutomata', 'WIDTH');
else
    validateattributes(width, {'numeric' 'logical'}, {'binary' 'row'}, ...
        'elementaryCellularAutomata', 'START');
end
if nargin < 4 || isempty(randfrac)
    dorand = false;
else
    validateattributes(randfrac, {'double' 'single'}, {'scalar' 'nonnegative' '<=' 1}, ...
        'elementaryCellularAutomata', 'FNOISE');
    dorand = true;
end

% set up machine
if isscalar(width)
    patt = ones(1, width);
    patt(floor((width+1)/2)) = 2;
else
    patt = width + 1;  % change 0,1 to 1,2 so can use sub2ind
    width = length(patt);
end

% unpack rule
rulearr = (bitget(rule, 1:8) + 1);

% initialise output array
pattern = zeros(n, width);

% iterate to generate rest of pattern
for i = 1:n
    pattern(i, :) = patt;   % record current state in output array

    % core step: apply CA rules to propagate to next 1D pattern
    ind = sub2ind([2 2 2], ...
        [patt(2:end) patt(1)], patt, [patt(end) patt(1:end-1)]);
    patt = rulearr(ind);

    %optional randomisation
    if dorand
        flip = rand(1, width) < randfrac;
        patt(flip) = 3 - patt(flip);
    end
end

% change symbols from 1 and 2 to 0 and 1
pattern = pattern-1;
end


text:
function pde = hyperIntf(am,ap,bm,bp,r,x0,y0,z0)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,...
    'exactu1',@exactu1,'exactu2',@exactu2,'exactu3',@exactu3,...
    'um1',@um1,'um2',@um2,'um3',@um3,'up1',@up1,'up2',@up2,'up3',@up3,...
    'Dxu',@Dxu,'Dxum',@Dxum,'Dxup',@Dxup,'Dyu',@Dyu,...
    'Dyum',@Dyum,'Dyup',@Dyup,'Dzu',@Dzu,'Dzum',@Dzum,'Dzup',@Dzup,...
    'f1',@f1,'f2',@f2,'f3',@f3,...
    'fm1',@fm1,'fm2',@fm2,'fm3',@fm3,...
    'fp1',@fp1,'fp2',@fp2,'fp3',@fp3,...
    'A',@A,'Am',@Am,'Ap',@Ap,'one',@one,...
    'B',@B,'Bm',@Bm,'Bp',@Bp);

pde.am = am;
pde.ap = ap;
pde.bm = bm;
pde.bp = bp;
%% interface function
    function u = intf(x,y,z)
        u = ((x-x0).^2 + (y-y0).^2 - (z-z0).^2)-r;
    end

%% exact solution
    function u = exactu1(x,y,z)
        u = um1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up1(x(id),y(id),z(id));
    end
    function u = exactu2(x,y,z)
        u = um2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up2(x(id),y(id),z(id));
    end
    function u = exactu3(x,y,z)
        u = um3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up3(x(id),y(id),z(id));
    end
    coef1 = 1; coef2 = 0; coef3 = 1;
    function u = um1(x,y,z)
        u = coef1*(z-z0) + coef3*(x-x0)/bm + coef2*intf(x,y,z).*(x-x0)/am;
    end
    function u = um2(x,y,z)
        u =                coef3*(y-y0)/bm + coef2*intf(x,y,z).*(y-y0)/am;
    end
    function u = um3(x,y,z)
        u = coef1*(x-x0) - coef3*(z-z0)/bm + coef2*intf(x,y,z).*(z-z0)/am;
    end
    function u = up1(x,y,z)
        u = coef1*(z-z0) + coef3*(x-x0)/bp + coef2*intf(x,y,z).*(x-x0)/ap;
    end
    function u = up2(x,y,z)
        u =                coef3*(y-y0)/bp + coef2*intf(x,y,z).*(y-y0)/ap;
    end
    function u = up3(x,y,z)
        u = coef1*(x-x0) - coef3*(z-z0)/bp + coef2*intf(x,y,z).*(z-z0)/ap;
    end
%% Boundary Function
    function u = gD1(x,y,z)
        u = exactu1(x,y,z);
    end
    function u = gD2(x,y,z)
        u = exactu2(x,y,z);
    end
    function u = gD3(x,y,z)
        u = exactu3(x,y,z);
    end
%% Derivative of the exact solution
    function u = Dxu(x,y,z)
        u = Dxum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dxup(x(id),y(id),z(id));
    end
    function u = Dyu(x,y,z)
        u = Dyum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dyup(x(id),y(id),z(id));
    end
    function u = Dzu(x,y,z)
        u = Dzum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dzup(x(id),y(id),z(id));
    end
    function u = Dxum(x,y,z)
        u = 4*(y-y0).*(z-z0)/am*coef2;
    end
    function u = Dyum(x,y,z)
        u = -4*(x-x0).*(z-z0)/am*coef2;
    end
    function u = Dzum(x,y,z)
        u = zeros(size(x));
    end
    function u = Dxup(x,y,z)
        u = 4*(y-y0).*(z-z0)/ap*coef2;
    end
    function u = Dyup(x,y,z)
        u = -4*(x-x0).*(z-z0)/ap*coef2;
    end
    function u = Dzup(x,y,z)
        u = zeros(size(x));
    end

%% right hand side function
    function u = f1(x,y,z)
        u = fm1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp1(x(id),y(id),z(id));
    end
    function u = f2(x,y,z)
        u = fm2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp2(x(id),y(id),z(id));
    end
    function u = f3(x,y,z)
        u = fm3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp3(x(id),y(id),z(id));
    end

    function u = fm1(x,y,z)
        u = 4*(x-x0)*coef2 + bm*um1(x,y,z);
    end
    function u = fm2(x,y,z)
        u = 4*(y-y0)*coef2 + bm*um2(x,y,z);
    end
    function u = fm3(x,y,z)
        u = -8*(z-z0)*coef2 + bm*um3(x,y,z);
    end
    function u = fp1(x,y,z)
        u = 4*(x-x0)*coef2 + bp*up1(x,y,z);
    end
    function u = fp2(x,y,z)
        u = 4*(y-y0)*coef2 + bp*up2(x,y,z);
    end
    function u = fp3(x,y,z)
        u = -8*(z-z0)*coef2 + bp*up3(x,y,z);
    end

%% Diffusion coefficient function
    function u = A(x,y,z)
        u = Am(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Ap(x(id),y(id),z(id));
    end
    function u = Am(x,y,z)
        u = am*ones(size(x));
    end
    function u = Ap(x,y,z)
        u = ap*ones(size(x));
    end

%% Mass coefficient function
    function u = B(x,y,z)
        u = Bm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Bp(x(id),y(id),z(id));
    end
    function u = Bm(x,y,z)
        u = bm*ones(size(x));
    end
    function u = Bp(x,y,z)
        u = bp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end

repo: ifem


text:
function outputImage = remove_clusters(this, varargin)
% Removes all voxel clusters within a given range of voxel counts, and/or
% a given range of percent of the area that is (approximately) filled;
% either slice-wise or volume-wise (8 or 26 neighbours) using Matlab's
% bwareaopen
%
%   Y = MrImage()
%   clusterRemovedImage = Y.remove_clusters('nPixelRange', nPixelRange, ...
%           'pAreaRange', pAreaRange, 'applicationDimension', applicationDimension)
%
% This is a method of class MrImage.
%
% IN
%   nPixelsPerClusterRange
%   percentAreaFilledRange
%           percent area of the cluster filled estimated using the
%           Solidity property
%   applicationDimension
%           dimensionality to perform operation
%           '2d' = slicewise application, separate 2d images (cluster
%           neighbours only considered within slice)
%           '3d' = neighbourhood considered as volume
%
% OUT
%   outputImage
%           MrImage where data matrix does not contain removed clusters
%
% EXAMPLE
%   Y = MrImage();
%   % remove all pixel clusters with 15 or less pixels (3D neighbourhood)
%   Y.remove_clusters('nPixelRange', [0 15], 'applicationDimension', '3D');
%   Y.remove_clusters('areaRange', [0 40], 'applicationDimension', '2D');
%
%   See also MrImage imdilate MrImage.imerode perform_unary_operation bwareaopen

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

defaults.nPixelRange = [0 0];
defaults.pAreaRange = [0 0];
defaults.applicationDimension = '3D';

args = tapas_uniqc_propval(varargin, defaults);
tapas_uniqc_strip_fields(args);

areaRange = pAreaRange/100;
is3D = strcmpi(applicationDimension, '3D');

BW = this.copyobj.binarize(0, 'exclude');

if is3D
    conn = 26; % use 26 neighbours in 3D (all faces and corners connected)
else
    conn = 8; % use 8 neighbours (in 2D)
end

%% step 1 - remove cluster based on voxel count

% only crop, if required (i.e., connectivity components greater than 1)
if nPixelRange(1) > 1
    % since we want to reinclude this mask later on (range!), we do not add
    % the +1 usually needed for bwareaopen, i.e., we only remove smaller
    % clusters here
    BW1 = BW.perform_unary_operation(@(x) bwareaopen(x, nPixelRange(1), conn), ...
        applicationDimension);
else
    BW1 = BW.copyobj;
end

if nPixelRange(2) > 0
    % since we indeed use this as a subtractive mask, we have to add +1 to
    % remove the clusters at the upper end of the cluster range as well
    BW2 = BW.perform_unary_operation(@(x) bwareaopen(x, nPixelRange(2) + 1, conn), ...
        applicationDimension);
else
    BW2 = BW.copyobj;
end

% by the following subtraction, excluded voxels in both BW1 and BW2 will be
% reincluded, hence, an exclusive or is reached, only removing clusters in
% the range between BW1 and BW2
BW = abs(BW - BW1 - BW2);

%% step 2 - remove cluster based on percent area filled
% extract region properties
if is3D
    stats = bwconncomp(BW.data, conn);
    props = regionprops3(stats, 'Solidity');
    % threshold solidity
    solidity = [props.Solidity];
    removeClusters = (solidity > areaRange(1) & solidity < areaRange(2));
    removeVoxels = vertcat(stats.PixelIdxList{removeClusters});
    % remove voxel
    BW.data(removeVoxels) = 0;
else
    % regionprops does not automatically take into account the defined
    % connectivity, so we have to do the split by hand
    nSamplesZ = BW.dimInfo.nSamples('z');
    for z = 1:nSamplesZ
        tempData = BW.select('z', z).data;
        stats = bwconncomp(tempData, conn);
        props = regionprops(stats, 'Solidity');
        % threshold solidity
        solidity = [props.Solidity];
        removeClusters = (solidity > areaRange(1) & solidity < areaRange(2));
        removeVoxels = vertcat(stats.PixelIdxList{removeClusters});
        % remove voxel
        tempData(removeVoxels) = 0;
        BW.data(:,:,z) = tempData;
    end
end

outputImage = this.*BW;

end

author: translationalneuromodeling

repo: tapas


text:
function varargout = pdftops(cmd)
%PDFTOPS  Calls a local pdftops executable with the input command
%
% Example:
%   [status result] = pdftops(cmd)
%
% Attempts to locate a pdftops executable, finally asking the user to
% specify the directory pdftops was installed into. The resulting path is
% stored for future reference.
%
% Once found, the executable is called with the input command string.
%
% This function requires that you have pdftops (from the Xpdf package)
% installed on your system. You can download this from:
% http://www.foolabs.com/xpdf
%
% IN:
%   cmd - Command string to be passed into pdftops.
%
% OUT:
%   status - 0 iff command ran without problem.
%   result - Output from pdftops.

% Copyright: Oliver Woodford, 2009-2010

% Thanks to Jonas Dorn for the fix for the title of the uigetdir window on
% Mac OS.
% Thanks to Christoph Hertel for pointing out a bug in check_xpdf_path
% under linux.

% Call pdftops
[varargout{1:nargout}] = system(sprintf('"%s" %s', xpdf_path, cmd));
return

function path = xpdf_path
% Return a valid path
% Start with the currently set path
path = user_string('pdftops');
% Check the path works
if check_xpdf_path(path)
    return
end
% Check whether the binary is on the path
if ispc
    bin = 'pdftops.exe';
else
    bin = 'pdftops';
end
if check_store_xpdf_path(bin)
    path = bin;
    return
end
% Search the obvious places
if ispc
    path = 'C:\Program Files\xpdf\pdftops.exe';
else
    path = '/usr/local/bin/pdftops';
end
if check_store_xpdf_path(path)
    return
end
% Ask the user to enter the path
while 1
    if strncmp(computer,'MAC',3) % Is a Mac
        % Give separate warning as the uigetdir dialogue box doesn't have a
        % title
        uiwait(warndlg('Pdftops not found. Please locate the program, or install xpdf-tools from http://users.phg-online.de/tk/MOSXS/.'))
    end
    base = uigetdir('/', 'Pdftops not found. Please locate the program.');
    if isequal(base, 0)
        % User hit cancel or closed window
        break;
    end
    base = [base filesep];
    bin_dir = {'', ['bin' filesep], ['lib' filesep]};
    for a = 1:numel(bin_dir)
        path = [base bin_dir{a} bin];
        if exist(path, 'file') == 2
            break;
        end
    end
    if check_store_xpdf_path(path)
        return
    end
end
error('pdftops executable not found.');

function good = check_store_xpdf_path(path)
% Check the path is valid
good = check_xpdf_path(path);
if ~good
    return
end
% Update the current default path to the path found
if ~user_string('pdftops', path)
    warning('Path to pdftops executable could not be saved. Enter it manually in pdftops.txt.');
    return
end
return

function good = check_xpdf_path(path)
% Check the path is valid
[good message] = system(sprintf('"%s" -h', path));
% system returns good = 1 even when the command runs
% Look for something distinct in the help text
good = ~isempty(strfind(message, 'PostScript'));
return

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% StackExchange Signal Processing Q52760
% https://dsp.stackexchange.com/questions/52760
% Convolution Strategy / Method for the Fastest 1D Convolution
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     27/04/2020
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

CONVOLUTION_METHOD_DFT          = 2;
CONVOLUTION_METHOD_OVERLAP_SAVE = 3;

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

vConvMethod = [CONVOLUTION_METHOD_DFT, CONVOLUTION_METHOD_OVERLAP_SAVE];
vConvShape  = [CONVOLUTION_SHAPE_FULL, CONVOLUTION_SHAPE_SAME, CONVOLUTION_SHAPE_VALID];

%% Simulation Parameters

vNumSamplesSignal = 2 * (10 .^ [2:0.5:6]);
vNumSamplesKernel = 10 .^ [1:0.5:6];
convShape   = CONVOLUTION_SHAPE_SAME;

%% Generate / Load Data

numSamplesSteps = length(vNumSamples);
numMethods      = length(vConvMethod);

vS = randn(vNumSamples(numSamplesSteps), 1);
vK = randn(vNumSamples(numSamplesSteps), 1);

%% Run Time Analysis

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        convString = 'full';
    case(CONVOLUTION_SHAPE_SAME)
        convString = 'same';
    case(CONVOLUTION_SHAPE_VALID)
        convString = 'valid';
end

mRunTime = zeros(numSamplesSteps, numSamplesSteps, numMethods);

for jj = 1:numSamplesSteps
    % Kernel Length
    numSamplesKernel = vNumSamples(jj);
    vKK = vK(1:numSamplesKernel);
    for ii = jj:numSamplesSteps
        % Signal Length
        numSamplesSignal = vNumSamples(ii);
        vSS = vS(1:numSamplesSignal);
        for kk = 1:numMethods
            switch(vConvMethod(kk))
                case(CONVOLUTION_METHOD_DIRECT)
                    hF = @() conv2(vSS, vKK, convString);
                case(CONVOLUTION_METHOD_OVERLAP_SAVE)
                    hF = @() ConvolutionOverlapSave(vSS, vKK, convShape);
                case(CONVOLUTION_METHOD_DFT)
                    hF = @() ConvolutionDft(vSS, vKK, convShape);
            end

            mRunTime(ii, jj, kk) = timeit(hF);

        end
    end
end

[~, mBestMethod] = min(mRunTime, [], 3);
for jj = 1:numSamplesSteps
    for ii = 1:numSamplesSteps
        if(jj > ii)
            mBestMethod(ii, jj) = 0;
        end
    end
end

%% Display Results

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge); %<! [x, y, width, height]
hAxes       = axes(); %<! [x, y, width, height]
hImageObj   = imagesc(vNumSamples, vNumSamples, mBestMethod);
set(hAxes, 'DataAspectRatio', [1, 1, 1]);
set(get(hAxes, 'Title'), 'String', {['Linear Convolution Run Time Comparison']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Kernel Length']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Signal Length']}, ...
    'FontSize', fontSizeTitle);
set(hAxes, 'LooseInset', [0.05, 0.05, 0.05, 0.05]);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [Price, ProbaITM , CI] = GetOptionPrice(Paths,Exercise,TimeToExpiry,RiskFreeRate,Optiontype)

switch Optiontype
    case 'Asian'

        MeanPrices = mean(Paths,1) ;
        OptionPricesAtMaturity = max(MeanPrices- Exercise,0);
           case 'Vanilla'
        OptionPricesAtMaturity = max(Paths(end,:)- Exercise,0);

end;
[MeanPriceAtMaturity, dummy,CI] = normfit(OptionPricesAtMaturity,0.01);
CI                     = CI* exp(-TimeToExpiry * RiskFreeRate);
Price                  = MeanPriceAtMaturity * exp(-TimeToExpiry * RiskFreeRate);
ProbaITM               = 100 * sum(OptionPricesAtMaturity > 0) ./ length(OptionPricesAtMaturity);


text:
%%*****************************************************************
%% linsysolve: solve linear system to get dy, and direction
%%             corresponding to unrestricted variables.
%%
%% [xx,coeff,L,resnrm] = linsysolve(schur,UU,EE,Bmat,rhs);
%%
%% child functions: mybicgstable.m
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [xx,coeff,L,resnrm] = HSDlinsysolve(par,schur,UU,EE,Bmat,rhs);

    global solve_ok msg
    global nnzmat nnzmatold matfct_options matfct_options_old use_LU

    spdensity  = par.spdensity;
    printlevel = par.printlevel;
    iter       = par.iter;

    m = length(schur);
    if (iter==1); use_LU = 0; matfct_options_old = ''; end
    if isempty(nnzmatold); nnzmatold = 0; end
%%
%% diagonal perturbation
%% old: pertdiag = 1e-15*max(1,diagschur);
%%
    diagschur = abs(full(diag(schur)));
    const = 1e-2/max(1,norm(par.dy2));
    alpha = max(1e-14,min(1e-10,const*norm(par.rp))/(1+norm(diagschur.*par.dy2)));
    pertdiag = alpha*max(1e-8,diagschur); %% Note: alpha is close to 1e-15.
    mexschurfun(schur,pertdiag);
    %%if (printlevel); fprintf(' %3.1e ',alpha); end
    if (par.depconstr) | (min(diagschur) < min([1e-20*max(diagschur), 1e-4]))
       lambda = 0.1*min(1e-14,const*norm(par.rp)/(1+norm(par.diagAAt.*par.dy2)));
       mexschurfun(schur,lambda*par.diagAAt);
       %%if (printlevel); fprintf('*'); end
    end
    if (max(diagschur)/min(diagschur) > 1e14) & (par.blkdim(2) == 0) ...
       & (iter > 10)
       tol = 1e-6;
       idx = find(diagschur < tol); len = length(idx);
       pertdiagschur = zeros(m,1);
       if (len > 0 & len < 5) & (norm(rhs(idx)) < tol)
          pertdiagschur(idx) = 1*ones(length(idx),1);
          mexschurfun(schur,pertdiagschur);
          if (printlevel); fprintf('#'); end
       end
    end
%%
%%
%%
    UU = [UU, Bmat];
    if ~isempty(EE)
       len = max(max(EE(:,1)),max(EE(:,2)));
    else
       len = 0;
    end
    tmp = [len+1,len+3,-1; len+2,len+4,1; len+3,len+1,1; len+4,len+2,-1;
           len+2,len+2,par.addschur];  %% this is the -inverse
    EE = [EE; tmp];
    ncolU = size(UU,2);
%%
%% assemble coefficient matrix
%%
    if isempty(EE)
       coeff.mat22 = [];
    else
       coeff.mat22 = spconvert(EE);
    end
    coeff.mat12 = UU;
    coeff.mat11 = schur; %% important to use perturbed schur matrix
%%
%% pad rhs with zero vector
%% decide which solution methods to use
%%
    rhs = [rhs; zeros(m+ncolU-length(rhs),1)];
    if (ncolU > 300); use_LU = 1; end
%%
%% Cholesky factorization
%%
    L = []; resnrm = []; xx = inf*ones(m,1);
    if (~use_LU)
       nnzmat = mexnnz(coeff.mat11);
       nnzmatdiff = (nnzmat ~= nnzmatold);
       solve_ok = 1;  solvesys = 1;
       if (nnzmat > spdensity*m^2) | (m < 500)
          matfct_options = 'chol';
       else
          matfct_options = 'spchol';
       end
       if (printlevel > 2); fprintf(' %s',matfct_options); end
       L.matdim = length(schur);
       if strcmp(matfct_options,'chol')
          if issparse(schur); schur = full(schur); end;
          if (iter<=5); %%--- to fix strange anonmaly in Matlab
             mexschurfun(schur,1e-20,2);
          end
          L.matfct_options = 'chol';
          [L.R,indef] = chol(schur);
          L.perm = [1:m];
       elseif strcmp(matfct_options,'spchol')
          if ~issparse(schur); schur = sparse(schur); end;
          L.matfct_options = 'spchol';
          [L.R,indef,L.perm] = chol(schur,'vector');
          L.Rt = L.R';
       end
       if (indef)
          solve_ok = -2; solvesys = 0;
          msg = 'HSDlinsysolve: Schur complement matrix not positive definite';
          if (printlevel); fprintf('\n  %s',msg); end
       end
       if (solvesys)
          if (ncolU)
             tmp = coeff.mat12'*linsysolvefun(L,coeff.mat12)-coeff.mat22;
	     if issparse(tmp); tmp = full(tmp); end
             [L.Ml,L.Mu,L.Mp] = lu(tmp);
             tol = 1e-16;
             condest = max(abs(diag(L.Mu)))/min(abs(diag(L.Mu)));
             idx = find(abs(diag(L.Mu)) < tol);
             if ~isempty(idx) | (condest > 1e50*sqrt(norm(par.diagAAt))); %% old: 1e30
                solvesys = 0; solve_ok = -4;
                use_LU = 1;
                msg = 'SMW too ill-conditioned, switch to LU factor';
                if (printlevel); fprintf('\n  %s, %2.1e.',msg,condest); end
             end
          end
          if (solvesys)
             [xx,resnrm,solve_ok] = HSDbicgstab(coeff,rhs,L,[],[],printlevel);
             if (solve_ok<=0) & (printlevel)
                fprintf('\n  warning: HSDbicgstab fails: %3.1f.',solve_ok);
             end
          end
       end
       if (solve_ok < 0)
          if (m < 6000 & strcmp(matfct_options,'chol')) | ...
             (m < 1e5 & strcmp(matfct_options,'spchol'))
             use_LU = 1;
             if (printlevel); fprintf('\n  switch to LU factor'); end
          end
       end
    end
%%
%% LU factorization
%%
    if (use_LU)
       nnzmat = mexnnz(coeff.mat11)+mexnnz(coeff.mat12);
       nnzmatdiff = (nnzmat ~= nnzmatold);
       solve_ok = 1;
       if ~isempty(coeff.mat22)
          raugmat = [coeff.mat11, coeff.mat12; coeff.mat12', coeff.mat22];
       else
          raugmat = coeff.mat11;
       end
       if (nnzmat > spdensity*m^2) | (m+ncolU < 500)
          matfct_options = 'lu';
       else
          matfct_options = 'splu';
       end
       if (printlevel > 2); fprintf(' %s ',matfct_options); end
       L.matdim = length(raugmat);
       if strcmp(matfct_options,'lu')
          if issparse(raugmat); raugmat = full(raugmat); end
          L.matfct_options = 'lu';
          [L.L,L.U,L.p] = lu(raugmat,'vector');
       elseif strcmp(matfct_options,'splu')
          if ~issparse(raugmat); raugmat = sparse(raugmat); end
          L.matfct_options = 'splu';
          [L.L,L.U,L.p,L.q,L.s] = lu(raugmat,'vector');
          L.s = full(diag(L.s));
       elseif strcmp(matfct_options,'ldl')
          if issparse(raugmat); raugmat = full(raugmat); end
          L.matfct_options = 'ldl';
          [L.L,L.D,L.p] = ldl(raugmat,'vector');
          L.D = sparse(L.D);
       elseif strcmp(matfct_options,'spldl')
          if ~issparse(raugmat); raugmat = sparse(raugmat); end
          L.matfct_options = 'spldl';
          [L.L,L.D,L.p,L.s] = ldl(raugmat,'vector');
          L.s  = full(diag(L.s));
          L.Lt = L.L';
       end
       [xx,resnrm,solve_ok] = HSDbicgstab(coeff,rhs,L,[],[],printlevel);
       if (solve_ok<=0) & (printlevel)
          fprintf('\n  warning: HSDbicgstab fails: %3.1f,',solve_ok);
       end
    end
    if (printlevel>=3); fprintf('%2.0f ',length(resnrm)-1); end
%%
    nnzmatold = nnzmat; matfct_options_old = matfct_options;
%%***************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [clust_num temp auto_sort] = find_temp(tree,clu,par) %,spikes,ipermut)

min_clus = par.min_clus;

c_ov = par.c_ov;
elbow_min = par.elbow_min;
clu = clu(1:end-1,3:end)+1; %first dim temp

y = tree(1:end,5);
maxdiff = max(diff(tree(1:end,6:end)),[],2);
maxdiff(maxdiff<0)=0;
prop = (y(2:end)+maxdiff.*(maxdiff>0))./y(1:end-1);
aux = find(prop<elbow_min,1,'first')+1; %percentaje of the rest

% The next if removes the particular case where just a class is found at the
% lowest temperature and with just a small change the rest appears
% all together at the next temperature
if ~isempty(aux) && par.mintemp==0 && aux==2
    aux = find(prop(2:end)<elbow_min,1,'first')+2; %percentaje of the rest
end

tree = tree(1:end-1,5:end);
clus = zeros(size(tree));
clus(tree(:,:) >= min_clus)=1; %only check the ones that cross the thr

dt = diff(tree);
clus = clus & [ones(size(clus(1,:)));dt(1:end,:)>min_clus];

for ii = 1:size(clus,1)
    detect = find(clus(ii,:),1,'last');
    if ~isempty(detect)
        clus(ii,1:detect)=1;
    end
end

auto_sort.elbow = size(tree,1);
if ~isempty(aux)
    clus(aux:end,1:end)=0;
    auto_sort.elbow = aux;
end
auto_sort.peaks = clus;

for ti = size(clus,1):-1:1
    detect = find(clus(ti,:));
    for ci = 1:length(detect) %the clusters removed aren't detected
        cl = (clu(ti,:) == detect(ci));
        for tj = ti-1:-1:1
            toremove = find(clus(tj,:));
            for j = 1:length(toremove)
                totest = (clu(tj,:) == toremove(j));
                if nnz(cl & totest)/min(nnz(totest),nnz(cl)) >= c_ov
                    clus(tj,toremove(j))=0;
                end
            end
        end
    end
end

[temp clust_num]=find(clus);

author: csn-le


text:
function p = pathdef
%PATHDEF Search path defaults.
%   PATHDEF returns a string that can be used as input to MATLABPATH
%   in order to set the path.

%   Copyright 1984-2007 The MathWorks, Inc.
%   $Revision: 1.4.2.2 $ $Date: 2007/06/07 14:45:14 $

% DO NOT MODIFY THIS FILE.  IT IS AN AUTOGENERATED FILE.
% EDITING MAY CAUSE THE FILE TO BECOME UNREADABLE TO
% THE PATHTOOL AND THE INSTALLER.

p = [...
%%% BEGIN ENTRIES %%%
     '/afs/inf.ed.ac.uk/user/s08/s0880377/matlab/statscopy:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/gitformatlab:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/junk64/dip/common/mlv7_4/diplib:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/junk64/dip/common/mlv7_4/dipimage_mex:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/junk64/dip/common/dipimage:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/residualanalysis:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/statfun:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/conjgradfunctions:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/simulationdatatool:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/fit/fitfunctions:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/fit:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/aux:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/rory:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/FastICA_25:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/image_proc:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/ploting:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/qdots:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/PatternAnalysis:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/PatternAnalysis/netlab:', ...
     '/afs/inf.ed.ac.uk/user/s08/s0880377/matlab:', ...
     matlabroot,'/toolbox/matlab/general:', ...
     matlabroot,'/toolbox/matlab/ops:', ...
     matlabroot,'/toolbox/matlab/lang:', ...
     matlabroot,'/toolbox/matlab/elmat:', ...
     matlabroot,'/toolbox/matlab/randfun:', ...
     matlabroot,'/toolbox/matlab/elfun:', ...
     matlabroot,'/toolbox/matlab/specfun:', ...
     matlabroot,'/toolbox/matlab/matfun:', ...
     matlabroot,'/toolbox/matlab/datafun:', ...
     matlabroot,'/toolbox/matlab/polyfun:', ...
     matlabroot,'/toolbox/matlab/funfun:', ...
     matlabroot,'/toolbox/matlab/sparfun:', ...
     matlabroot,'/toolbox/matlab/scribe:', ...
     matlabroot,'/toolbox/matlab/graph2d:', ...
     matlabroot,'/toolbox/matlab/graph3d:', ...
     matlabroot,'/toolbox/matlab/specgraph:', ...
     matlabroot,'/toolbox/matlab/graphics:', ...
     matlabroot,'/toolbox/matlab/uitools:', ...
     matlabroot,'/toolbox/matlab/strfun:', ...
     matlabroot,'/toolbox/matlab/imagesci:', ...
     matlabroot,'/toolbox/matlab/iofun:', ...
     matlabroot,'/toolbox/matlab/audiovideo:', ...
     matlabroot,'/toolbox/matlab/timefun:', ...
     matlabroot,'/toolbox/matlab/datatypes:', ...
     matlabroot,'/toolbox/matlab/verctrl:', ...
     matlabroot,'/toolbox/matlab/codetools:', ...
     matlabroot,'/toolbox/matlab/helptools:', ...
     matlabroot,'/toolbox/matlab/demos:', ...
     matlabroot,'/toolbox/matlab/timeseries:', ...
     matlabroot,'/toolbox/matlab/hds:', ...
     matlabroot,'/toolbox/matlab/guide:', ...
     matlabroot,'/toolbox/matlab/plottools:', ...
     matlabroot,'/toolbox/local:', ...
     matlabroot,'/toolbox/shared/controllib:', ...
     matlabroot,'/toolbox/shared/dastudio:', ...
     matlabroot,'/toolbox/matlab/datamanager:', ...
     matlabroot,'/toolbox/compiler:', ...
     matlabroot,'/toolbox/distcomp:', ...
     matlabroot,'/toolbox/distcomp/user:', ...
     matlabroot,'/toolbox/distcomp/mpi:', ...
     matlabroot,'/toolbox/distcomp/pctdemos:', ...
     matlabroot,'/toolbox/distcomp/parallel:', ...
     matlabroot,'/toolbox/distcomp/parallel/datafun:', ...
     matlabroot,'/toolbox/distcomp/parallel/datatypes:', ...
     matlabroot,'/toolbox/distcomp/parallel/elfun:', ...
     matlabroot,'/toolbox/distcomp/parallel/elmat:', ...
     matlabroot,'/toolbox/distcomp/parallel/lapack:', ...
     matlabroot,'/toolbox/distcomp/parallel/matfun:', ...
     matlabroot,'/toolbox/distcomp/parallel/ops:', ...
     matlabroot,'/toolbox/distcomp/parallel/sparfun:', ...
     matlabroot,'/toolbox/distcomp/parallel/specfun:', ...
     matlabroot,'/toolbox/distcomp/parallel/util:', ...
     matlabroot,'/toolbox/distcomp/lang:', ...
     matlabroot,'/toolbox/gads:', ...
     matlabroot,'/toolbox/gads/gads:', ...
     matlabroot,'/toolbox/gads/gadsdemos:', ...
     matlabroot,'/toolbox/images/colorspaces:', ...
     matlabroot,'/toolbox/images/images:', ...
     matlabroot,'/toolbox/images/imdemos:', ...
     matlabroot,'/toolbox/images/imuitools:', ...
     matlabroot,'/toolbox/images/iptformats:', ...
     matlabroot,'/toolbox/images/iptutils:', ...
     matlabroot,'/toolbox/shared/imageslib:', ...
     matlabroot,'/toolbox/shared/spcuilib:', ...
     matlabroot,'/toolbox/nnet:', ...
     matlabroot,'/toolbox/nnet/nncontrol:', ...
     matlabroot,'/toolbox/nnet/nndemos:', ...
     matlabroot,'/toolbox/nnet/nnet:', ...
     matlabroot,'/toolbox/nnet/nnet/nnanalyze:', ...
     matlabroot,'/toolbox/nnet/nnet/nncustom:', ...
     matlabroot,'/toolbox/nnet/nnet/nndistance:', ...
     matlabroot,'/toolbox/nnet/nnet/nnformat:', ...
     matlabroot,'/toolbox/nnet/nnet/nninit:', ...
     matlabroot,'/toolbox/nnet/nnet/nnlearn:', ...
     matlabroot,'/toolbox/nnet/nnet/nnnetinput:', ...
     matlabroot,'/toolbox/nnet/nnet/nnnetwork:', ...
     matlabroot,'/toolbox/nnet/nnet/nnperformance:', ...
     matlabroot,'/toolbox/nnet/nnet/nnplot:', ...
     matlabroot,'/toolbox/nnet/nnet/nnprocess:', ...
     matlabroot,'/toolbox/nnet/nnet/nnsearch:', ...
     matlabroot,'/toolbox/nnet/nnet/nntopology:', ...
     matlabroot,'/toolbox/nnet/nnet/nntrain:', ...
     matlabroot,'/toolbox/nnet/nnet/nntransfer:', ...
     matlabroot,'/toolbox/nnet/nnet/nnweight:', ...
     matlabroot,'/toolbox/nnet/nnguis:', ...
     matlabroot,'/toolbox/nnet/nnguis/nftool:', ...
     matlabroot,'/toolbox/nnet/nnguis/nntool:', ...
     matlabroot,'/toolbox/nnet/nnobsolete:', ...
     matlabroot,'/toolbox/nnet/nnresource:', ...
     matlabroot,'/toolbox/nnet/nnutils:', ...
     matlabroot,'/toolbox/nnet/nnguis/nnguiutils:', ...
     matlabroot,'/toolbox/nnet/nndemos/nndatasets:', ...
     matlabroot,'/toolbox/nnet/nnguis/nntraintool:', ...
     matlabroot,'/toolbox/optim/optim:', ...
     matlabroot,'/toolbox/optim/optimdemos:', ...
     matlabroot,'/toolbox/shared/optimlib:', ...
     matlabroot,'/toolbox/signal/signal:', ...
     matlabroot,'/toolbox/signal/sigtools:', ...
     matlabroot,'/toolbox/signal/sptoolgui:', ...
     matlabroot,'/toolbox/signal/sigdemos:', ...
     matlabroot,'/toolbox/shared/siglib:', ...
     matlabroot,'/toolbox/stats:', ...
     matlabroot,'/toolbox/shared/statslib:', ...
     matlabroot,'/toolbox/wavelet/wavelet:', ...
     matlabroot,'/toolbox/wavelet/wmultisig1d:', ...
     matlabroot,'/toolbox/wavelet/wavedemo:', ...
     matlabroot,'/toolbox/wavelet/compression:', ...
     matlabroot,'/toolbox/bioinfo/bioinfo:', ...
     matlabroot,'/toolbox/bioinfo/biolearning:', ...
     matlabroot,'/toolbox/bioinfo/microarray:', ...
     matlabroot,'/toolbox/bioinfo/mass_spec:', ...
     matlabroot,'/toolbox/bioinfo/proteins:', ...
     matlabroot,'/toolbox/bioinfo/biomatrices:', ...
     matlabroot,'/toolbox/bioinfo/biodemos:', ...
     matlabroot,'/toolbox/bioinfo/graphtheory:', ...
     matlabroot,'/toolbox/curvefit/curvefit:', ...
     matlabroot,'/toolbox/curvefit/cftoolgui:', ...
     matlabroot,'/toolbox/simbio/simbio:', ...
     matlabroot,'/toolbox/simbio/simbiodemos:', ...
%%% END ENTRIES %%%
     ...
];

p = [userpath,p];


text:
%--- help for abstvar/variance_decomposition ---
%
%  Compute the variance decompsition of the VAR
%
%  ::
%
%     [vd,retcode] = variance_decomposition(self,params,Rfunc,shock_names,varargin);
%
%  Args:
%     self (var object):
%     params :(optional) parameter values
%     Rfunc (function handle): (optional) transform parameters into dynamics
%     shock_names : names of shocks
%     varargin
%
%  Returns:
%     :
%
%     - **vd** (struct): a struct containing the variance decomposition:
%
%        - **infinity** (struct):
%        - **conditional** (struct):
%
%     - **retcode** (integer): This passes through the retcode given by Rfunc. Currently (2018/07/19), this output always returns 0.
%
%
%    Other functions named variance_decomposition
%
%       generic/variance_decomposition
%

author: jmaih


text:
function [pred_boxes, scores] = rfcn_im_detect(conf, caffe_net, im, boxes, max_rois_num_in_gpu)
% [pred_boxes, scores] = rfcn_im_detect(conf, caffe_net, im, boxes, max_rois_num_in_gpu)
% --------------------------------------------------------
% R-FCN implementation
% Modified from MATLAB Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2016, Jifeng Dai
% Licensed under The MIT License [see LICENSE for details]

    [im_blob, rois_blob, ~] = get_blobs(conf, im, boxes);

    % When mapping from image ROIs to feature map ROIs, there's some aliasing
    % (some distinct image ROIs get mapped to the same feature ROI).
    % Here, we identify duplicate feature ROIs, so we only compute features
    % on the unique subset.
    [~, index, inv_index] = unique(rois_blob, 'rows');
    rois_blob = rois_blob(index, :);
    boxes = boxes(index, :);

    % permute data into caffe c++ memory, thus [num, channels, height, width]
    im_blob = im_blob(:, :, [3, 2, 1], :); % from rgb to brg
    im_blob = permute(im_blob, [2, 1, 3, 4]);
    im_blob = single(im_blob);
    rois_blob = rois_blob - 1; % to c's index (start from 0)
    rois_blob = permute(rois_blob, [3, 4, 2, 1]);
    rois_blob = single(rois_blob);

    total_rois = size(rois_blob, 4);
    total_scores = cell(ceil(total_rois / max_rois_num_in_gpu), 1);
    total_box_deltas = cell(ceil(total_rois / max_rois_num_in_gpu), 1);
    for i = 1:ceil(total_rois / max_rois_num_in_gpu)

        sub_ind_start = 1 + (i-1) * max_rois_num_in_gpu;
        sub_ind_end = min(total_rois, i * max_rois_num_in_gpu);
        sub_rois_blob = rois_blob(:, :, :, sub_ind_start:sub_ind_end);

        net_inputs = {im_blob, sub_rois_blob};

        % Reshape net's input blobs
        caffe_net.reshape_as_input(net_inputs);
        caffe_net.forward(net_inputs);

        if conf.test_binary
            % simulate binary logistic regression
            scores = caffe_net.blobs('cls_score').get_data();
            scores = squeeze(scores)';
            % Return scores as fg - bg
            scores = bsxfun(@minus, scores, scores(:, 1));
        else
            % use softmax estimated probabilities
            scores = caffe_net.blobs('cls_prob').get_data();
            scores = squeeze(scores)';
        end

        % Apply bounding-box regression deltas
        box_deltas = caffe_net.blobs('bbox_pred').get_data();
        box_deltas = squeeze(box_deltas)';

        total_scores{i} = scores;
        total_box_deltas{i} = box_deltas;
    end

    scores = cell2mat(total_scores);
    box_deltas = cell2mat(total_box_deltas);

    pred_boxes = rfcn_bbox_transform_inv(boxes, box_deltas);
    pred_boxes = clip_boxes(pred_boxes, size(im, 2), size(im, 1));

    % Map scores and predictions back to the original set of boxes
    scores = scores(inv_index, :);
    pred_boxes = pred_boxes(inv_index, :);

    % remove scores and boxes for back-ground
    pred_boxes = pred_boxes(:, 5:end);
    scores = scores(:, 2:end);
    if conf.bbox_class_agnostic
        pred_boxes = repmat(pred_boxes, [1, size(scores,2)]);
    end
end

function [data_blob, rois_blob, im_scale_factors] = get_blobs(conf, im, rois)
    [data_blob, im_scale_factors] = get_image_blob(conf, im);
    rois_blob = get_rois_blob(conf, rois, im_scale_factors);
end

function [blob, im_scales] = get_image_blob(conf, im)
    [ims, im_scales] = arrayfun(@(x) prep_im_for_blob(im, conf.image_means, x, conf.test_max_size), conf.test_scales, 'UniformOutput', false);
    im_scales = cell2mat(im_scales);
    blob = im_list_to_blob(ims);
end

function [rois_blob] = get_rois_blob(conf, im_rois, im_scale_factors)
    [feat_rois, levels] = map_im_rois_to_feat_rois(conf, im_rois, im_scale_factors);
    rois_blob = single([levels, feat_rois]);
end

function [feat_rois, levels] = map_im_rois_to_feat_rois(conf, im_rois, scales)
    im_rois = single(im_rois);

    if length(scales) > 1
        widths = im_rois(:, 3) - im_rois(:, 1) + 1;
        heights = im_rois(:, 4) - im_rois(:, 2) + 1;

        areas = widths .* heights;
        scaled_areas = bsxfun(@times, areas(:), scales(:)'.^2);
        [~, levels] = min(abs(scaled_areas - 224.^2), [], 2);
    else
        levels = ones(size(im_rois, 1), 1);
    end

    feat_rois = round(bsxfun(@times, im_rois-1, scales(levels))) + 1;
end

function boxes = clip_boxes(boxes, im_width, im_height)
    % x1 >= 1 & <= im_width
    boxes(:, 1:4:end) = max(min(boxes(:, 1:4:end), im_width), 1);
    % y1 >= 1 & <= im_height
    boxes(:, 2:4:end) = max(min(boxes(:, 2:4:end), im_height), 1);
    % x2 >= 1 & <= im_width
    boxes(:, 3:4:end) = max(min(boxes(:, 3:4:end), im_width), 1);
    % y2 >= 1 & <= im_height
    boxes(:, 4:4:end) = max(min(boxes(:, 4:4:end), im_height), 1);
end

repo: R-FCN


text:
function [ n_data, x, fx ] = bessel_y1_values ( n_data )

%*****************************************************************************80
%
%% BESSEL_Y1_VALUES returns some values of the Y1 Bessel function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      BesselY[1,x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 16;

  fx_vec = [ ...
     -0.6458951094702027E+01, ...
     -0.7812128213002887E+00, ...
     -0.1070324315409375E+00, ...
      0.3246744247918000E+00, ...
      0.3979257105571000E+00, ...
      0.1478631433912268E+00, ...
     -0.1750103443003983E+00, ...
     -0.3026672370241849E+00, ...
     -0.1580604617312475E+00, ...
      0.1043145751967159E+00, ...
      0.2490154242069539E+00, ...
      0.1637055374149429E+00, ...
     -0.5709921826089652E-01, ...
     -0.2100814084206935E+00, ...
     -0.1666448418561723E+00, ...
      0.2107362803687351E-01 ];

  x_vec = [ ...
      0.1E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      3.0E+00, ...
      4.0E+00, ...
      5.0E+00, ...
      6.0E+00, ...
      7.0E+00, ...
      8.0E+00, ...
      9.0E+00, ...
     10.0E+00, ...
     11.0E+00, ...
     12.0E+00, ...
     13.0E+00, ...
     14.0E+00, ...
     15.0E+00  ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function smallpot = marginalize_pot(bigpot, keep, maximize, useC)
% MARGINALIZE_POT Marginalize a cpot onto a smaller domain.
% smallpot = marginalize_pot(bigpot, keep, maximize, useC)
%
% The maximize argument is ignored - maxing out a Gaussian is the same as summing it out,
% since the mode and mean are equal.
% The useC argument is ignored.

node_sizes = sparse(1, max(bigpot.domain));
node_sizes(bigpot.domain) = bigpot.sizes;
sum_over = mysetdiff(bigpot.domain, keep);

if sum(node_sizes(sum_over))==0 % isempty(sum_over)
  %smallpot = bigpot;
  smallpot = cpot(keep, node_sizes(keep), bigpot.g, bigpot.h, bigpot.K);
else
  [h1, h2, K11, K12, K21, K22] = partition_matrix_vec(bigpot.h, bigpot.K, sum_over, keep, node_sizes);
  n = length(h1);
  K11inv = inv(K11);
  g = bigpot.g + 0.5*(n*log(2*pi) - log(det(K11)) + h1'*K11inv*h1);
  if length(h2) > 0 % ~isempty(keep) % we are are actually keeping something
    A = K21*K11inv;
    h = h2 - A*h1;
    K = K22 - A*K12;
  else
    h = [];
    K = [];
  end
  smallpot = cpot(keep, node_sizes(keep), g, h, K);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function Z = or(X,Y)
%OR Logical OR (|) for tensors.
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

Z = tenfun(@or,X,Y);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
classdef MexFunctionIDs < double
    %MexFunctionID Enumeration that provides ID number of desired function
    %  Copyright 2017 The MathWorks, Inc.

    enumeration
        OPEN_LIB (0)
        CLOSE_LIB (1)
        GET_JNT_POS (2)
        GET_JNT_VEL (3)
        GET_JNT_TORQUE (4)
        GET_JNT_TEMP (5)
        GET_EE_POSE (6)
        GET_EE_WRENCH (7)
        MOVE_TO_HOME (8)
        SET_POS_CTRL (9)
        SET_DIRECT_TORQUE_CTRL (10)
        RUN_GRAVITY_CALIB (11)
        INIT_FINGERS (12)
        SEND_JNT_POS (13)
        SEND_JNT_FING_POS (14)
        SEND_JNT_VEL (15)
        SEND_JNT_TORQUE (16)
        SEND_FINGER_POS (17)
        SEND_CART_POS (18)
        SEND_CART_VEL (19)
        GET_DOF (20)
        START_FORCE_CONTROL(21)
        STOP_FORCE_CONTROL(22)
        GET_EE_OFFSET(23)
        SET_EE_OFFSET(24)
        GET_PROTECT_ZONE(25)
        ERR_PROTECT_ZONES(26)
        SET_PROTECT_ZONE(27)
        GET_GLOB_TRAJECTORY (28)
    end

    methods
    end

end

author: PatternRecognition

repo: OpenBMI


text:
function [cost, gradf, hessian, hessianinv]=costgradf(px,py,pz,rho,u,param)
% compute the cost function

e=1e-6;

h     = param.h;
dim   = param.dim;
gamma = param.gamma;
A1x   = param.A1x;
A1y   = param.A1y;
A1z   = param.A1z;
A2    = param.A2;
A3    = param.A3;
a     = param.a;
c     = param.c;
F1    = param.F1;
F2    = param.F2;

n  = dim(1);
nx = dim(2);
ny = dim(3);
nz = dim(4);
nt = dim(5);
hx = h(1);
hy = h(2);
hz = h(3);
ht = h(4);

%compute p^2 and u^2

pxsq  = px.^2;
pysq  = py.^2;
pzsq  = pz.^2;
usq   = u.^2;

% inverse of rho
rhoinv  = 1./rho;
rhoinv1 = 1./(F1'*rho(:));
rhoinv2 = 1./(F2'*rho(:));

E=eye(n);
E(n,n)=e;
W1=kron(speye((nx-1)*ny*nz*nt),E);
W2=kron(speye((ny-1)*nx*nz*nt),E);
W3=kron(speye((nz-1)*nx*ny*nt),E);

% cost f
%cost = (A1x*W1*pxsq(:)+A1y*W2*pysq(:)+A1z*W3*pzsq(:))'*(A2*rhoinv(:) + a(:))*hx*hy*hz*ht + usq(:)'*(A3*rhoinv1(:)+A3*rhoinv2(:)+ c(:))*hx*hy*hz*ht*gamma;
cost = (A1x*W1*pxsq(:)+A1y*W2*pysq(:)+A1z*W3*pzsq(:))'*(A2*rhoinv(:) + a(:))*hx*hy*hz*ht + usq(:)'*(A3*rhoinv2(:)+ c(:))*hx*hy*hz*ht*gamma;
% compute gradient of f
if nargout > 1
    % derivative wrt px
    A11xe    = 2*W1*A1x'*(A2*rhoinv(:)+a(:));
    nablapxf = px(:).*A11xe(:);

    % derivative wrt py
    A11ye    = 2*W2*A1y'*(A2*rhoinv(:)+a(:));
    nablapyf = py(:).*A11ye(:);

    % derivative wrt pz
    A11ze    = 2*W3*A1z'*(A2*rhoinv(:)+a(:));
    nablapzf = pz(:).*A11ze(:);

    % derivative wrt rho
    A22e     = A2'*(A1x*W1*pxsq(:) + A1y*W2*pysq(:)+ A1z*W3*pzsq(:));
    nablarhof = -A22e.*rhoinv(:).^2-gamma.*F2*((A3'*usq(:)).*rhoinv2(:).^2);
                %-gamma.*F1*((A3'*usq(:)).*rhoinv1(:).^2);

    % derivative wrt u
    %A33e    = 2*gamma*(A3*(rhoinv2(:)+rhoinv1(:))+c(:));
    A33e    = 2*gamma*(A3*(rhoinv2(:))+c(:));
    nablauf = u(:).*A33e(:);

    % gradient wrt p rho u
    gradf   = [nablapxf; nablapyf; nablapzf; nablarhof; nablauf];

    % compute Hessian
    if nargout > 2
        % Hessian over px
        A11x    = A11xe;

        % Hessian over py
        A11y    = A11ye;

        % Hessian over py
        A11z    = A11ze;

         % hessian over rho
        temp1   = 2.*A22e.*rhoinv(:).^3;
        temp2   = 2.*gamma.*(F2*((A3'*usq(:)).*rhoinv2(:).^3));
        %temp3   = 2.*gamma.*(F1*((A3'*usq(:)).*rhoinv1(:).^3));
        A22     = temp1(:)+temp2(:)...+temp3(:)
        +1e-4*ones(n*nx*ny*nz*(nt-1),1);

        % hessian over u
        A33  = A33e;

        % hessian
        hessian = [A11x(:);A11y(:);A11z(:);A22(:);A33(:)];

        % compute the inverse of Hessian
        if nargout > 3
            % inverse of hessian
            hessianinv = spdiags(1./(hessian),0,numel(hessian),numel(hessian));
        end
    end
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function fe2d_n_fast_test ( )

%*****************************************************************************80
%
%% FE2D_N_FAST_TEST tests the FE2D_N_FAST code.
%
%  Discussion:
%
%    This function sets all parameter values and initial condition information
%    necessary to execute the "fast" version of the fe2d_n algorithm.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    28 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Reference:
%
%    Marcus R Garvie, John Burkardt, Jeff Morgan,
%    Simple Finite Element Methods for Approximating Predator-Prey Dynamics
%    in Two Dimensions using MATLAB,
%    Submitted to Bulletin of Mathematical Biology, 2014.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FE2D_N_FAST_TEST:\n' );
  fprintf ( 1, '  Test the FE2D_N_FAST function\n' );
  fprintf ( 1, '  which applies Neumann boundary conditions as it\n' );
  fprintf ( 1, '  approximates a solution to a predator-prey system.\n' );
%
%  Set the parameters.
%
  alpha = 0.4;
  beta = 2.0;
  gamma = 0.6;
  delta = 1.0;
%
%  Use T=150.0 for normal run.
%  Use T=0.50 for a "quick" run that might take 15 minutes of computing.
%
% T = 150.0;
  T = 0.50;
  delt = 1.0 / 384.0;

  t = tic;
  fe2d_n_fast ( alpha, beta, gamma, delta, T, delt, @u0f, @v0f, @guf, @gvf );
  t = toc ( t );

  fprintf ( 1, '  Execution took %10.2g minutes \n', t / 60.0 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FE2D_N_FAST_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

function value = u0f ( x, y )

%*****************************************************************************80
%
%% U0F evaluates the initial condition for U.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    26 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location in the region.
%
%    Output, real VALUE, the initial condition for U at (X,Y).
%
  value = 6.0 / 35.0 - 2.0E-07 * ( x - 0.1 * y - 225.0 ) * ( x - 0.1 * y - 675.0 );

  return
end

function value = v0f ( x, y )

%*****************************************************************************80
%
%% V0F evaluates the initial condition for V.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    26 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location in the region.
%
%    Output, real VALUE, the initial condition for V at (X,Y).
%
  value = 116.0 / 245.0 - 3.0E-05 * ( x - 450.0 ) - 1.2E-04 * ( y - 150.0 );

  return
end

function value = guf ( x, y, t )

%*****************************************************************************80
%
%% GUF evaluates the Neumann boundary condition for U.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    28 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location on the boundary.
%
%    Input, real T, the time.
%
%    Output, real VALUE, the prescribed value of dU/dn at (X,Y,T).
%
  value = 0.0;

  return
end
function value = gvf ( x, y, t )

%*****************************************************************************80
%
%% GVF evaluates the Neumann boundary condition for V.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    28 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location on the boundary.
%
%    Input, real T, the time.
%
%    Output, real VALUE, the prescribed value of dV/dn at (X,Y,T).
%
  value = 0.0;

  return
end


text:
function [gVarmeans gVarcovs gDynKern] = modelPriorReparamGrads(dynModel, gVarmeansLik, gVarcovsLik,gInd)
% MODELPRIORREPARAMGRADS Wrapper function for the gradients of the various types of the
% variational GPLVM bound when dynamics is used.
% FORMAT
% DESC provides a wrapper function for the variational GP-LVM, when there
% are dynamics. It takes a dynamics model structure and according to its type it
% calls the appropriate function to calculate the gradient for the
% variational bound. The gradients returned are w.r.t mu_bar and lambda, i.e. the parameters used by the optimiser
% (the ones introduced by the reparametrization), not the original ones (mu and S).
% These gradients are obtained in two steps: firstly, another funtion is used to calculate the
% gradients that correspond only to the likelihood part of the bound and
% are w.r.t the original parameters mu and S. These quantities are the arguments
% gVarmenasLik and gVarcovsLik and are calculated using the already
% implemented code for the static vargplvm, which assumes that the original
% parameters are not coupled.
% The current funtions receives these quantities and a) amends
% with partial derivatives because the variational parameters are coupled
% via Kt for the dyn. gplvm b) Computes the whole gradients for both parts
% of the bound, the one corresponding to the likelihood and the one
% corresponding to the prior. This must be done in a single function
% because the final formula contains both parts in a nonlinear form.
% c) Also the kernel hyperparameters for the dynamics kernel are being retured.
%
% See the dyn. vargplvm notes for more details.
%
% ARG dynModel : the dynamics model structure for which the gradients are
% to be computed.
% ARG gVarmeansLik, gVarcovsLik: the gradients for the VAR-GPLVM model computed
% for the original parameters and only for the likelihood term.
% ARG gInd: in case the inducing points are tied to the variational means
% this is the partial derivatives of the likelihood part of the variational
% bound w.r.t the inducing points X_u, otherwise it is just [].
% RETURN gVarmeans, gVarcovs : the gradients for the "reparametrized" means
% and covariances (the ones that are visible to the optimiser, not the
% original ones) for the VAR-GPLVM model.
% RETURN gDynKern: the gradient w.r.t the hyperparameters of the dynamics
% kernel
%
% SEEALSO : vargplvmLogLikeGradients vargpTimeDynamicsPriorReparamGrads
%
% COPYRIGHT : Michalis K. Titsias, 2010-2011
% COPYRIGHT : Neil D. Lawrence, 2010-2011
% COPYRIGHT : Andreas C. Damianou, 2010-2011

% VARGPLVM

fhandle = str2func([dynModel.type 'PriorReparamGrads']);
[gVarmeans gVarcovs gDynKern] = fhandle(dynModel, gVarmeansLik, gVarcovsLik, gInd);

author: SheffieldML

repo: vargplvm


text:
function [km] = in2km(in)
% Convert length from inches to kilometers.
% Chad Greene 2012
km = in*0.0000254;


text:
function [cost,source] = cs_results(px,py,pz,rho,u,param)

e=1e-6;

h     = param.h;
dim   = param.dim;
gamma = param.gamma;
A1x   = param.A1x;
A1y   = param.A1y;
A1z   = param.A1z;
A2    = param.A2;
A3    = param.A3;
a     = param.a;
c     = param.c;
F1    = param.F1;
F2    = param.F2;

n  = dim(1);
nx = dim(2);
ny = dim(3);
nz = dim(4);
nt = dim(5);
hx = h(1);
hy = h(2);
hz = h(3);
ht = h(4);

%compute p^2 and u^2

pxsq  = px.^2;
pysq  = py.^2;
pzsq  = pz.^2;
usq   = u.^2;

% inverse of rho
rhoinv  = 1./rho;
rhoinv2 = 1./(F2'*rho(:));

E=eye(n);
E(n,n)=e;
W1=kron(speye((nx-1)*ny*nz*nt),E);
W2=kron(speye((ny-1)*nx*nz*nt),E);
W3=kron(speye((nz-1)*nx*ny*nt),E);

% cost f
cost = (A1x*W1*pxsq(:)+A1y*W2*pysq(:)+A1z*W3*pzsq(:))'*(A2*rhoinv(:) + a(:))*hx*hy*hz*ht + usq(:)'*(A3*rhoinv2(:)+ c(:))*hx*hy*hz*ht*gamma;
source=usq(:)'*(A3*rhoinv2(:)+ c(:))*hx*hy*hz*ht;
end


text:
function wedge_monte_carlo_test ( )

%*****************************************************************************80
%
%% WEDGE_MONTE_CARLO_TEST tests the WEDGE_MONTE_CARLO library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 August 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WEDGE_MONTE_CARLO_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the WEDGE_MONTE_CARLO library.\n' );

  wedge_monte_carlo_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WEDGE_MONTE_CARLO_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [ResultsMaleFemale] = compareMaleFemale(male,female)
% This function compares basic features of the male and female whole-body
% metabolic models
%
% [ResultsMaleFemale] = compareMaleFemale(male,female)
%
% INPUT
% male                  model structure (male whole-body metabolic model)
% female                model structure (female whole-body metabolic model)
%
% OUTPUT
% ResultsMaleFemale     structure containing the basic differences and
%                       commenalities between male and female model
%
% Ines Thiele 2017

% reactions unique to male
ResultsMaleFemale.MaleOnly = setdiff(male.rxns,female.rxns);
ResultsMaleFemale.FemaleOnly = setdiff(female.rxns,male.rxns);
ResultsMaleFemale.BothGender = intersect(female.rxns,male.rxns);

[maleOrgans]=unique(strtok(male.rxns,'_'));
[femaleOrgans]=unique(strtok(female.rxns,'_'));

for i = 1 : length(maleOrgans)
    ResultsMaleFemale.OrgansNumRxnMale(i,1) = length(strmatch(maleOrgans(i),male.rxns));
    ResultsMaleFemale.OrgansNumRxnMale(i,2) = length(strmatch(maleOrgans(i),ResultsMaleFemale.MaleOnly));
    % fraction
    ResultsMaleFemale.OrgansNumRxnMale(i,3) = ResultsMaleFemale.OrgansNumRxnMale(i,2)/ResultsMaleFemale.OrgansNumRxnMale(i,1);
end

for i = 1 : length(femaleOrgans)
    ResultsMaleFemale.OrgansNumRxnFemale(i,1) = length(strmatch(femaleOrgans(i),female.rxns));
    ResultsMaleFemale.OrgansNumRxnFemale(i,2) = length(strmatch(femaleOrgans(i),ResultsMaleFemale.FemaleOnly));
    % fraction
    ResultsMaleFemale.OrgansNumRxnFemale(i,3) = ResultsMaleFemale.OrgansNumRxnFemale(i,2)/ResultsMaleFemale.OrgansNumRxnFemale(i,1);
end

ResultsMaleFemale.maleOrgans = maleOrgans;
ResultsMaleFemale.femaleOrgans = femaleOrgans;

%get subsystems for gall rxns
FemaleSS = female.subSystems(find(ismember(female.rxns,ResultsMaleFemale.FemaleOnly(strmatch('Gall_',ResultsMaleFemale.FemaleOnly)))));
ResultsMaleFemale.FemaleGallSSEnrich = unique(FemaleSS);
for i = 1 : length(ResultsMaleFemale.FemaleGallSSEnrich)
    ResultsMaleFemale.FemaleGallSSEnrich{i,2} = num2str(length(strmatch(ResultsMaleFemale.FemaleGallSSEnrich{i},FemaleSS,'exact')));
end
MaleSS = male.subSystems(find(ismember(male.rxns,ResultsMaleFemale.MaleOnly(strmatch('Gall_',ResultsMaleFemale.MaleOnly)))));
ResultsMaleFemale.MaleGallSSEnrich = unique(MaleSS);
for i = 1 : length(ResultsMaleFemale.MaleGallSSEnrich)
    ResultsMaleFemale.MaleGallSSEnrich{i,2} = num2str(length(strmatch(ResultsMaleFemale.MaleGallSSEnrich{i},MaleSS,'exact')));
end

% unique biofluid exchange reactions

ResultsMaleFemale.MaleOnlyBiofluid = ResultsMaleFemale.MaleOnly(find(~cellfun(@isempty,strfind(ResultsMaleFemale.MaleOnly,'_EX_'))));
ResultsMaleFemale.FemaleOnlyBiofluid = ResultsMaleFemale.FemaleOnly(find(~cellfun(@isempty,strfind(ResultsMaleFemale.FemaleOnly,'_EX_'))));

author: opencobra

repo: cobratoolbox


text:
function F = lininterp1f(X,Y,XI,Ydefault)

disp('For fast linear interpolation, download lininterp1f from Mathworks')
disp('file exchange (File ID: #8627) or use alternative interpolator')
disp('http://www.mathworks.com/matlabcentral/fileexchange/8627-fast-linear-interpolation')
disp(' ')
disp('Place lininterp1f in CORE\utilities')


text:
function diff_vol = dtiSmoothAnisoPM(vol, num_iter, delta_t, kappa, option, voxel_spacing)
%ANISODIFF2D Conventional anisotropic diffusion
%   DIFF_VOL = ANISODIFF3D(VOL, NUM_ITER, DELTA_T, KAPPA, OPTION, VOXEL_SPACING) perfoms
%   conventional anisotropic diffusion (Perona & Malik) upon a stack of gray scale images.
%   A 3D network structure of 26 neighboring nodes is considered for diffusion conduction.
%
%       ARGUMENT DESCRIPTION:
%               VOL      - gray scale volume data (MxNxP).
%               NUM_ITER - number of iterations.
%               DELTA_T  - integration constant (0 <= delta_t <= 3/44).
%                          Usually, due to numerical stability this
%                          parameter is set to its maximum value.
%               KAPPA    - gradient modulus threshold that controls the conduction.
%               OPTION   - conduction coefficient functions proposed by Perona & Malik:
%                          1 - c(x,y,z,t) = exp(-(nablaI/kappa).^2),
%                              privileges high-contrast edges over low-contrast ones.
%                          2 - c(x,y,z,t) = 1./(1 + (nablaI/kappa).^2),
%                              privileges wide regions over smaller ones.
%          VOXEL_SPACING - 3x1 vector column with the x, y and z dimensions of
%                          the voxel (milimeters). In particular, only cubic and
%                          anisotropic voxels in the z-direction are considered.
%                          When dealing with DICOM images, the voxel spacing
%                          dimensions can be extracted using MATLAB's dicominfo(.).
%
%       OUTPUT DESCRIPTION:
%               DIFF_VOL - (diffused) volume with the largest scale-space parameter.
%
%   Example
%   -------------
%   vol = randn(100,100,100);
%   num_iter = 4;
%   delta_t = 3/44;
%   kappa = 70;
%   option = 2;
%   voxel_spacing = ones(3,1);
%   diff_vol = anisodiff3D(vol, num_iter, delta_t, kappa, option, voxel_spacing);
%   figure, subplot 121, imshow(vol(:,:,50),[]), subplot 122, imshow(diff_vol(:,:,50),[])
%
% See also anisodiff1D, anisodiff2D.

% References:
%   P. Perona and J. Malik.
%   Scale-Space and Edge Detection Using Anisotropic Diffusion.
%   IEEE Transactions on Pattern Analysis and Machine Intelligence,
%   12(7):629-639, July 1990.
%
%   G. Grieg, O. Kubler, R. Kikinis, and F. A. Jolesz.
%   Nonlinear Anisotropic Filtering of MRI Data.
%   IEEE Transactions on Medical Imaging,
%   11(2):221-232, June 1992.
%
%   MATLAB implementation based on Peter Kovesi's anisodiff(.):
%   P. D. Kovesi. MATLAB and Octave Functions for Computer Vision and Image Processing.
%   School of Computer Science & Software Engineering,
%   The University of Western Australia. Available from:
%   <http://www.csse.uwa.edu.au/~pk/research/matlabfns/>.
%
% Credits:
% Daniel Simoes Lopes
% ICIST
% Instituto Superior Tecnico - Universidade Tecnica de Lisboa
% danlopes (at) civil ist utl pt
% http://www.civil.ist.utl.pt/~danlopes
%
% May 2007 original version.

% Convert input volume to double.
vol = double(vol);

% Useful variables.
[rows cols pags] = size(vol);

% PDE (partial differential equation) initial condition.
diff_vol = vol;
clear vol

% Center voxel distances.
x = voxel_spacing(1);
y = voxel_spacing(2);
z = voxel_spacing(3);
dx = 1;
dy = 1;
dz = z/x;
dd = sqrt(dx^2+dy^2);
dh = sqrt(dx^2+dz^2);
dc = sqrt(dd^2+dz^2);

% 3D convolution masks - finite differences.
h1 = zeros(3,3,3); h1(2,2,2) = -1; h1(2,2,1) = 1;
h2 = zeros(3,3,3); h2(2,2,2) = -1; h2(2,2,3) = 1;
h3 = zeros(3,3,3); h3(2,2,2) = -1; h3(2,1,2) = 1;
h4 = zeros(3,3,3); h4(2,2,2) = -1; h4(2,3,2) = 1;
h5 = zeros(3,3,3); h5(2,2,2) = -1; h5(3,2,2) = 1;
h6 = zeros(3,3,3); h6(2,2,2) = -1; h6(1,2,2) = 1;

h7 = zeros(3,3,3); h7(2,2,2) = -1; h7(3,1,1) = 1;
h8 = zeros(3,3,3); h8(2,2,2) = -1; h8(2,1,1) = 1;
h9 = zeros(3,3,3); h9(2,2,2) = -1; h9(1,1,1) = 1;
h10 = zeros(3,3,3); h10(2,2,2) = -1; h10(3,2,1) = 1;
h11 = zeros(3,3,3); h11(2,2,2) = -1; h11(1,2,1) = 1;
h12 = zeros(3,3,3); h12(2,2,2) = -1; h12(3,3,1) = 1;
h13 = zeros(3,3,3); h13(2,2,2) = -1; h13(2,3,1) = 1;
h14 = zeros(3,3,3); h14(2,2,2) = -1; h14(1,3,1) = 1;

h15 = zeros(3,3,3); h15(2,2,2) = -1; h15(3,1,2) = 1;
h16 = zeros(3,3,3); h16(2,2,2) = -1; h16(1,1,2) = 1;
h17 = zeros(3,3,3); h17(2,2,2) = -1; h17(3,3,2) = 1;
h18 = zeros(3,3,3); h18(2,2,2) = -1; h18(1,3,2) = 1;

h19 = zeros(3,3,3); h19(2,2,2) = -1; h19(3,1,3) = 1;
h20 = zeros(3,3,3); h20(2,2,2) = -1; h20(2,1,3) = 1;
h21 = zeros(3,3,3); h21(2,2,2) = -1; h21(1,1,3) = 1;
h22 = zeros(3,3,3); h22(2,2,2) = -1; h22(3,2,3) = 1;
h23 = zeros(3,3,3); h23(2,2,2) = -1; h23(1,2,3) = 1;
h24 = zeros(3,3,3); h24(2,2,2) = -1; h24(3,3,3) = 1;
h25 = zeros(3,3,3); h25(2,2,2) = -1; h25(2,3,3) = 1;
h26 = zeros(3,3,3); h26(2,2,2) = -1; h26(1,3,3) = 1;

% Anisotropic diffusion.
for t = 1:num_iter

    % Finite differences. [imfilter(.,.,'conv') can be replaced by convn(.,.,'same')]
    % Due to possible memory limitations, the diffusion
    % will be calculated at each page/slice of the volume.
    for p = 1:pags-2
        diff3pp = diff_vol(:,:,p:p+2);
        aux = imfilter(diff3pp,h1,'conv'); nabla1 = aux(:,:,2);
        aux = imfilter(diff3pp,h2,'conv'); nabla2 = aux(:,:,2);
        aux = imfilter(diff3pp,h3,'conv'); nabla3 = aux(:,:,2);
        aux = imfilter(diff3pp,h4,'conv'); nabla4 = aux(:,:,2);
        aux = imfilter(diff3pp,h5,'conv'); nabla5 = aux(:,:,2);
        aux = imfilter(diff3pp,h6,'conv'); nabla6 = aux(:,:,2);
        aux = imfilter(diff3pp,h7,'conv'); nabla7 = aux(:,:,2);
        aux = imfilter(diff3pp,h8,'conv'); nabla8 = aux(:,:,2);
        aux = imfilter(diff3pp,h9,'conv'); nabla9 = aux(:,:,2);
        aux = imfilter(diff3pp,h10,'conv'); nabla10 = aux(:,:,2);
        aux = imfilter(diff3pp,h11,'conv'); nabla11 = aux(:,:,2);
        aux = imfilter(diff3pp,h12,'conv'); nabla12 = aux(:,:,2);
        aux = imfilter(diff3pp,h13,'conv'); nabla13 = aux(:,:,2);
        aux = imfilter(diff3pp,h14,'conv'); nabla14 = aux(:,:,2);
        aux = imfilter(diff3pp,h15,'conv'); nabla15 = aux(:,:,2);
        aux = imfilter(diff3pp,h16,'conv'); nabla16 = aux(:,:,2);
        aux = imfilter(diff3pp,h17,'conv'); nabla17 = aux(:,:,2);
        aux = imfilter(diff3pp,h18,'conv'); nabla18 = aux(:,:,2);
        aux = imfilter(diff3pp,h19,'conv'); nabla19 = aux(:,:,2);
        aux = imfilter(diff3pp,h20,'conv'); nabla20 = aux(:,:,2);
        aux = imfilter(diff3pp,h21,'conv'); nabla21 = aux(:,:,2);
        aux = imfilter(diff3pp,h22,'conv'); nabla22 = aux(:,:,2);
        aux = imfilter(diff3pp,h23,'conv'); nabla23 = aux(:,:,2);
        aux = imfilter(diff3pp,h24,'conv'); nabla24 = aux(:,:,2);
        aux = imfilter(diff3pp,h25,'conv'); nabla25 = aux(:,:,2);
        aux = imfilter(diff3pp,h26,'conv'); nabla26 = aux(:,:,2);

        % Diffusion function.
        if option == 1
            c1 = exp(-(nabla1/kappa).^2);
            c2 = exp(-(nabla2/kappa).^2);
            c3 = exp(-(nabla3/kappa).^2);
            c4 = exp(-(nabla4/kappa).^2);
            c5 = exp(-(nabla5/kappa).^2);
            c6 = exp(-(nabla6/kappa).^2);
            c7 = exp(-(nabla7/kappa).^2);
            c8 = exp(-(nabla8/kappa).^2);
            c9 = exp(-(nabla9/kappa).^2);
            c10 = exp(-(nabla10/kappa).^2);
            c11 = exp(-(nabla11/kappa).^2);
            c12 = exp(-(nabla12/kappa).^2);
            c13 = exp(-(nabla13/kappa).^2);
            c14 = exp(-(nabla14/kappa).^2);
            c15 = exp(-(nabla15/kappa).^2);
            c16 = exp(-(nabla16/kappa).^2);
            c17 = exp(-(nabla17/kappa).^2);
            c18 = exp(-(nabla18/kappa).^2);
            c19 = exp(-(nabla19/kappa).^2);
            c20 = exp(-(nabla20/kappa).^2);
            c21 = exp(-(nabla21/kappa).^2);
            c22 = exp(-(nabla22/kappa).^2);
            c23 = exp(-(nabla23/kappa).^2);
            c24 = exp(-(nabla24/kappa).^2);
            c25 = exp(-(nabla25/kappa).^2);
            c26 = exp(-(nabla26/kappa).^2);
        elseif option == 2
            c1 = 1./(1 + (nabla1/kappa).^2);
            c2 = 1./(1 + (nabla2/kappa).^2);
            c3 = 1./(1 + (nabla3/kappa).^2);
            c4 = 1./(1 + (nabla4/kappa).^2);
            c5 = 1./(1 + (nabla5/kappa).^2);
            c6 = 1./(1 + (nabla6/kappa).^2);
            c7 = 1./(1 + (nabla7/kappa).^2);
            c8 = 1./(1 + (nabla8/kappa).^2);
            c9 = 1./(1 + (nabla9/kappa).^2);
            c10 = 1./(1 + (nabla10/kappa).^2);
            c11 = 1./(1 + (nabla11/kappa).^2);
            c12 = 1./(1 + (nabla12/kappa).^2);
            c13 = 1./(1 + (nabla13/kappa).^2);
            c14 = 1./(1 + (nabla14/kappa).^2);
            c15 = 1./(1 + (nabla15/kappa).^2);
            c16 = 1./(1 + (nabla16/kappa).^2);
            c17 = 1./(1 + (nabla17/kappa).^2);
            c18 = 1./(1 + (nabla18/kappa).^2);
            c19 = 1./(1 + (nabla19/kappa).^2);
            c20 = 1./(1 + (nabla20/kappa).^2);
            c21 = 1./(1 + (nabla21/kappa).^2);
            c22 = 1./(1 + (nabla22/kappa).^2);
            c23 = 1./(1 + (nabla23/kappa).^2);
            c24 = 1./(1 + (nabla24/kappa).^2);
            c25 = 1./(1 + (nabla25/kappa).^2);
            c26 = 1./(1 + (nabla26/kappa).^2);
        end

    % Discrete PDE solution.
    diff_vol(:,:,p+1) = diff_vol(:,:,p+1) + ...
                        delta_t*(...
                        (1/(dz^2))*c1.*nabla1 + (1/(dz^2))*c2.*nabla2 + ...
                        (1/(dx^2))*c3.*nabla3 + (1/(dx^2))*c4.*nabla4 + ...
                        (1/(dy^2))*c5.*nabla5 + (1/(dy^2))*c6.*nabla6 + ...
                        ...
                        (1/(dc^2))*c7.*nabla7 + (1/(dh^2))*c8.*nabla8 + ...
                        (1/(dc^2))*c9.*nabla9 + (1/(dh^2))*c10.*nabla10 + ...
                        (1/(dh^2))*c11.*nabla11 + (1/(dc^2))*c12.*nabla12 + ...
                        (1/(dh^2))*c13.*nabla13 + (1/(dc^2))*c14.*nabla14 + ...
                        ...
                        (1/(dd^2))*c15.*nabla15 + (1/(dd^2))*c16.*nabla16 + ...
                        (1/(dd^2))*c17.*nabla17 + (1/(dd^2))*c18.*nabla18 + ...
                        ...
                        (1/(dc^2))*c19.*nabla19 + (1/(dh^2))*c20.*nabla20 + ...
                        (1/(dc^2))*c21.*nabla21 + (1/(dh^2))*c22.*nabla22 + ...
                        (1/(dh^2))*c23.*nabla23 + (1/(dc^2))*c24.*nabla24 + ...
                        (1/(dh^2))*c25.*nabla25 + (1/(dc^2))*c26.*nabla26);
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function data = json_decode(string)
% json_decode Parses a JSON string
%
% This function parses a JSON string and returns a cell array with the
% parsed data. JSON objects are converted to structures and JSON arrays are
% converted to cell arrays.
%
% Input:
% - string (string): A JSON string
%
% Output:
% - data (struct, cell): A decoded Matlab structure
%
% Credits: Based on the code of F. Glineur, 2009 (Available on Matlab File
% Exchange)

pos = 1;
len = length(string);
% String delimiters and escape characters are identified beforehand to improve speed
esc = regexp(string, '["\\]'); index_esc = 1; len_esc = length(esc);

if pos <= len
    switch(next_char)
        case '{'
            data = parse_object;
        case '['
            data = parse_array;
        otherwise
            error_pos('Outer level structure must be an object or an array');
    end
end

    function object = parse_object
        parse_char('{');
        object = [];
        if next_char ~= '}'
            while 1
                str = parse_string;
                if isempty(str)
                    error_pos('Name of value at position %d cannot be empty');
                end
                parse_char(':');
                val = parse_value;
                object.(valid_field(str)) = val;
                if next_char == '}'
                    break;
                end
                parse_char(',');
            end
        end
        parse_char('}');
    end

    function object = parse_array
        parse_char('[');
        object = cell(0, 1);
        if next_char ~= ']'
            while 1
                val = parse_value;
                object{end+1} = val;
                if next_char == ']'
                    break;
                end
                parse_char(',');
            end
        end
        parse_char(']');
    end

    function parse_char(c)
        skip_whitespace;
        if pos > len || string(pos) ~= c
            error_pos(sprintf('Expected %c at position %%d', c));
        else
            pos = pos + 1;
            skip_whitespace;
        end
    end

    function c = next_char
        skip_whitespace;
        if pos > len
            c = [];
        else
            c = string(pos);
        end
    end

    function skip_whitespace
        while pos <= len && isspace(string(pos))
            pos = pos + 1;
        end
    end

     function str = parse_string
        if string(pos) ~= '"'
            error_pos('String starting with " expected at position %d');
        else
            pos = pos + 1;
        end
        str = '';
        while pos <= len
            while index_esc <= len_esc && esc(index_esc) < pos
                index_esc = index_esc + 1;
            end
            if index_esc > len_esc
                str = [str string(pos:end)];
                pos = len + 1;
                break;
            else
                str = [str string(pos:esc(index_esc)-1)];
                pos = esc(index_esc);
            end
            switch string(pos)
                case '"'
                    pos = pos + 1;
                    return;
                case '\'
                    if pos+1 > len
                        error_pos('End of file reached right after escape character');
                    end
                    pos = pos + 1;
                    switch string(pos)
                        case {'"' '\' '/'}
                            str(end+1) = string(pos);
                            pos = pos + 1;
                        case {'b' 'f' 'n' 'r' 't'}
                            str(end+1) = sprintf(['\' string(pos)]);
                            pos = pos + 1;
                        case 'u'
                            if pos+4 > len
                                error_pos('End of file reached in escaped unicode character');
                            end
                            str(end+1:end+6) = string(pos-1:pos+4);
                            pos = pos + 5;
                    end
                otherwise % should never happen
                    str(end+1) = string(pos);
                    pos = pos + 1;
            end
        end
        error_pos('End of file while expecting end of string');
    end

    function num = parse_number
        [num, ~, err, delta] = sscanf(string(pos:min(len,pos+20)), '%f', 1); % TODO : compare with json(pos:end)
        if ~isempty(err)
            error_pos('Error reading number at position %d');
        end
        pos = pos + delta-1;
    end

    function val = parse_value
        switch(string(pos))
            case '"'
                val = parse_string;
                return;
            case '['
                val = parse_array;
                return;
            case '{'
                val = parse_object;
                return;
            case {'-','0','1','2','3','4','5','6','7','8','9'}
                val = parse_number;
                return;
            case 't'
                if pos+3 <= len && strcmpi(string(pos:pos+3), 'true')
                    val = true;
                    pos = pos + 4;
                    return;
                end
            case 'f'
                if pos+4 <= len && strcmpi(string(pos:pos+4), 'false')
                    val = false;
                    pos = pos + 5;
                    return;
                end
            case 'n'
                if pos+3 <= len && strcmpi(string(pos:pos+3), 'null')
                    val = [];
                    pos = pos + 4;
                    return;
                end
        end
        error_pos('Value expected at position %d');
    end

    function error_pos(msg)
        poss = max(min([pos-15 pos-1 pos pos+20],len),1);
        if poss(3) == poss(2)
            poss(3:4) = poss(2)+[0 -1];         % display nothing after
        end
        msg = [sprintf(msg, pos) ' : ... ' string(poss(1):poss(2)) '<error>' string(poss(3):poss(4)) ' ... '];
        ME = MException('JSONparser:invalidFormat', msg);
        throw(ME);
    end

    function str = valid_field(str)
    % From MATLAB doc: field names must begin with a letter, which may be
    % followed by any combination of letters, digits, and underscores.
    % Invalid characters will be converted to underscores, and the prefix
    % "alpha_" will be added if first character is not a letter.
        if ~isletter(str(1))
            str = ['alpha_' str];
        end
        str(~isletter(str) & ~('0' <= str & str <= '9')) = '_';
    end

end

author: votchallenge

repo: toolkit-legacy


text:
function [inv_map,bi_inv_map,logdet,invA] = invchol_or_lu(A)
% INVCHOL_OR_LU
% Does a Cholesky decomposition on A and returns logdet, inverse and
% two function handles that respectively map X to A\X and A\X/A.
%

if nargin==0
    test_this();
    return;
end

if isreal(A)

    R = chol(A); %R'*R = A
    inv_map = @(X) R\(R'\X);

    % inv(A)*X*inv(A)
    bi_inv_map = @(X) (inv_map(X)/R)/(R');

    if nargout>2
        logdet = 2*sum(log(diag(R)));
        if nargout>3
            invA = inv_map(eye(size(A)));
        end
    end
else
    [L,T,p] = lu(A,'vector');
    P = sparse(p,1:length(p),1);
                       % P*A = L*U
                       % L is lower triangular, with unit diagonal and unit determinant
                       % T is upper triangular, det(T) = prod(diag(T), may have negative values on diagonal
                       % P is a permuation matrix: P' = inv(P) and det(P) is +1 or -1
                       %

    inv_map = @(X) T\(L\(P*X));

    % inv(A)*X*inv(A)
    bi_inv_map = @(X) ((inv_map(X)/T)/L)*P;

    if nargout>2
        logdet = sum(log(diag(T)))-log(det(P));
        if nargout>3
            invA = T\(L\P);
        end
    end
end

function test_this()
dim = 3;
A = randn(dim)+sqrt(-1)*randn(dim);
[inv_map,bi_inv_map,logdet,iA] = invchol_or_lu(A);
[logdet,log(det(A))]
X = randn(dim,3);
Y1 = A\X,
Y2 = inv_map(X)

Z1 = (A\X)/A
Z2 = bi_inv_map(X)

iA,
inv(A)

author: nesl