text:
function g = lmcKernGradient(kern, X, X2, covGrad)

% LMCKERNGRADIENT Gradient of LMC kernel's parameters.
% FORMAT
% DESC computes the gradient of parameters associated to the LMC kernel.
% As well as the kernel structure and the input positions, the user
% provides a matrix COVGRAD which gives the partial derivatives of the
% function with respect to the relevant elements of the kernel matrix.
% RETURN g:  gradients of the function of interest with respect to the
% kernel parameters. The ordering of the vector should match that
% provided by the function kernExtractParam.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG X : the input locations for which the gradients are being computed.
% ARG covGrad : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes the
% form of a square matrix of dimension  numData, where numData is the
% number of rows in X.
%
% FORMAT
% DESC  computes the derivatives as above, but input locations are now
% provided in two matrices associated with rows and columns of the kernel
% matrix.
% RETURN g : gradients of the function of interest with respect to the
% kernel parameters.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG X : the input locations associated with the rows of the kernel matrix.
% ARG X2 : the input locations associated with the columns of the kernel
% matrix.
% ARG covGrad : matrix of partial derivatives of the function of interest
% with respect to the kernel matrix. The matrix should have the same number
% of rows as X1 and the same number of columns as X2 has rows.
%
% COPYRIGHT : Mauricio A. Alvarez, 2010.

% KERN

fhandle = str2func([kern.basicKernelType 'KernCompute']);
fhandleGrad = str2func([kern.basicKernelType 'KernGradient']);
gBK = zeros(1, kern.nParamsBK);
gPartialB = zeros(kern.nout);

if iscell(X)
    if nargin>3 && ~iscell(X2)
        error('Time course information is not matched in Cell format!');
    end
    % Collate arguments.
    dim1 = zeros(1, kern.nout);
    dim2 = zeros(1, kern.nout);
    for i=1:kern.nout
        dim1(i) = size(X{i}, 1);
        if nargin > 3
            dim2(i) = size(X2{i}, 1);
        else
            dim2(i) = dim1(i);
            covGrad = X2;
        end
    end
    for i = 1:kern.nout
        startOne = sum(dim1(1:(i-1)))+1;
        endOne = sum(dim1(1:i));
        startThree = sum(dim2(1:(i-1))) + 1;
        endThree = sum(dim2(1:i));
        if nargin > 3
            gBK = gBK + kern.B(i,i)*fhandleGrad(kern, X{i}, X2{i}, covGrad(startOne:endOne, ...
                startThree:endThree));
            basicK = fhandle(kern, X{i}, X2{i});
        else
            gBK = gBK + kern.B(i,i)*fhandleGrad(kern, X{i}, covGrad(startOne:endOne,...
                startThree:endThree));
            basicK = fhandle(kern, X{i});
        end
        gPartialB(i,i) = sum(sum(covGrad(startOne:endOne, startThree:endThree).*basicK));
        for j = 1:i-1
            startTwo = sum(dim2(1:(j-1))) + 1;
            endTwo =  sum(dim2(1:j));
            if nargin > 3
                g2  = kern.B(i,j)*fhandleGrad(kern, X{i}, X2{j}, covGrad(startOne:endOne, ...
                    startTwo:endTwo));
                basicK = fhandle(kern, X{i}, X2{j});
            else
                g2  = kern.B(i,j)*fhandleGrad(kern, X{i}, X{j}, covGrad(startOne:endOne, ...
                    startTwo:endTwo));
                basicK = fhandle(kern, X{i}, X{j});
            end
            gBK = gBK + 2*g2;
            gPartialB(i,j) = sum(sum(covGrad(startOne:endOne, startTwo:endTwo).*basicK));
            gPartialB(j,i) = gPartialB(i,j);
        end
    end
else
    if nargin < 4
        covGrad = X2;
        X2 = X;
    end
    basicK = fhandle(kern, X, X2);
    startOne = 1;
    endOne = 0;
    startThree = 1;
    endThree = 0;
    for i=1:kern.nout
        endOne = endOne + size(X,1);
        endThree = endThree + size(X2,1);
        gBK = gBK + kern.B(i,i)*fhandleGrad(kern, X, X2, ...
            covGrad(startOne:endOne, startThree:endThree));
        gPartialB(i,i) = sum(sum(covGrad(startOne:endOne, startThree:endThree).*basicK));
        startTwo = 1;
        endTwo = 0;
        startFour = 1;
        endFour = 0;
        for j=1:i-1
            endTwo = endTwo + size(X2,1);
            g2  = kern.B(i,j)*fhandleGrad(kern, X, X2, covGrad(startOne:endOne, ...
                startTwo:endTwo));
            gBK = gBK + g2;
            gPartialB(i,j) = sum(sum(covGrad(startOne:endOne, startTwo:endTwo).*basicK));
            if nargin < 3
                gBK = gBK + g2;
                gPartialB(j,i) = gPartialB(i,j);
            else
                endFour = endFour + size(X,1);
                g3 = kern.B(j,i)*fhandleGrad(kern, X, X2, covGrad(startFour:endFour, ...
                    startThree:endThree));
                gBK = gBK + g3;
                gPartialB(j,i) = sum(sum(covGrad(startFour:endFour, startThree:endThree).*basicK));
                startFour = endFour + 1;
            end
            startTwo = endTwo + 1;
        end
        startOne = endOne + 1;
        startThree = endThree + 1;
    end
end

Jij = zeros(kern.nout, kern.rankCorregMatrix);
Jji = zeros(kern.rankCorregMatrix, kern.nout);
gB = zeros(kern.nout, kern.rankCorregMatrix);
for i=1:kern.nout
    for j=1:kern.rankCorregMatrix
        Jij(i,j) = 1; Jji(j,i) = 1;
        gB(i,j) = sum(sum((kern.A*Jij' + Jji'*kern.A').*gPartialB));
        Jij(i,j) = 0; Jji(j,i) = 0;
    end
end

g = [gBK gB(:)'];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%MEDIANBLUR  Blurs an image using the median filter
%
%     dst = cv.medianBlur(src)
%     dst = cv.medianBlur(src, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ Input 1-, 3-, or 4-channel image. When `KSize` is 3 or 5, the
%   image type should be `uint8`, `uint16`, or `single`. For larger aperture
%   sizes, it can only be `uint8`.
%
% ## Output
% * __dst__ Destination array of the same size and type as `src`.
%
% ## Options
% * __KSize__ Aperture linear size. It must be odd and greater than 1, for
%   example 3, 5, 7 ... default 5
%
% The function smooths an image using the median filter with the
% `KSize x KSize` aperture. Each channel of a multi-channel image is
% processed independently.
%
% Note: The median filter uses `BorderType=Replicate` internally to cope with
% border pixels. See cv.copyMakeBorder.
%
% See also: cv.bilateralFilter, cv.blur, cv.boxFilter, cv.GaussianBlur,
%  medfilt2
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% this function detects potential malware in the current folder and subfolders
%
% Author: A. Delorme, Cotober 2013

function detectmalware(currentFolder);

if nargin < 1
    currentFolder = pwd;
end;

folderContent = dir(currentFolder);
folderContent = { folderContent.name };

malwareStrings = { 'eval(' 'evalin(' 'evalc(' 'delete(' 'movefile(' 'rmdir' 'mkdir' 'copyfile' 'system(' '!' };

for iFile = 1:length(folderContent)

    currentFile = folderContent{iFile};
    if length(currentFile) > 2 && strcmpi(currentFile(end-1:end), '.m')

        fid = fopen(fullfile(currentFolder, currentFile), 'r');
        countLine = 0;
        prevstr = '';
        while ~feof(fid)
            str = fgetl(fid);
            countLine = countLine+1;

            if length(str) > 1 && str(1) ~= '%'
                res = cellfun(@(x)~isempty(findstr(x, str)), malwareStrings(1:end-1));
                if str(1) == '!', res(end+1) = 1; end;
                if any(res)
                    pos = find(res); pos = pos(1);
                    disp('************************************')
                    fprintf('Potential malware command detected containing "%s" in\n %s line %d\n', malwareStrings{pos}, fullfile(currentFolder, currentFile), countLine);
                    fprintf('%d: %s\n%d: %s\n%d: %s\n', countLine-1, prevstr, countLine, str, countLine+1, fgetl(fid));
                    countLine = countLine+1;
                end;
            end;

            prevstr = str;
        end;
        fclose(fid);

    elseif exist(currentFile) == 7 && ~strcmpi(currentFile, '..') && ~strcmpi(currentFile, '.')
        detectmalware(fullfile(currentFolder, currentFile));
    end;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ alpha1, alpha2, beta, beta1, beta2, czero, difeta, etazero, fzero, ...
  gamma1, gamma2, kappa, lambda1, lambda2, m, nu1, nu2, vzero, xl, xr ] = ...
  bio_constants ( )

%*****************************************************************************80
%
%% BIO_CONSTANTS sets biological constants.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2005
%
%  Parameters:
%
%    Output, real XL, XR, the left and right endpoints of the region.
%

%
%  The extent of the domain.
%
  xl = 0.0;
  xr = 1.0;
%
%  The initial conditions.
%
  vzero = 15.0;
  kappa = 0.00554245;
  m = 10;
  czero = 0.0;
  fzero = 1.0;
  etazero = 1.0;
%
%  Parameters occurring in equations.
%
  lambda1 = 73.0;
  lambda2 = 19.0;
  nu1 = 0.007;
  nu2 = 1.28;
  difeta = 0.0144;
  beta = 0.222;
%
%  Parameters for calculation of TAU.
%
  alpha1 = 0.001;
  alpha2 = 1.0;
  beta1 = 1.0;
  beta2 = 0.001;
  gamma1 = 1.2;
  gamma2 = 1.2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear; close all; clc;

%%

filePath=mfilename('fullpath');
toolboxPath=fileparts(fileparts(fileparts(filePath)));

libFolder=fullfile(toolboxPath,'lib');
helpFolder=fullfile(toolboxPath,'docs');
outputFolder=fullfile(toolboxPath,'docs');

%%

%Get list of M-files in library
fileListLib = dir(fullfile(libFolder,'*.m'));
fileListLib={fileListLib(1:end).name};
fileListLib=sort(fileListLib(:));
NumberOfLibFiles=numel(fileListLib);

%Get list of M-files in help folder
fileListHelp = dir(fullfile(helpFolder,'*.m'));
fileListHelp={fileListHelp(1:end).name};
fileListHelp=sort(fileListHelp(:));
NumberOfHelpFiles=numel(fileListHelp);

%%

licenseBoilerPlate=fullfile(toolboxPath,'licenseBoilerPlate.txt');
[T_license]=txtfile2cell(licenseBoilerPlate);

footerTargetText='% _*GIBBON footer text*_ ';
licenseLink='https://github.com/gibbonCode/GIBBON/blob/master/LICENSE';

%Add comment symbols
for q=1:1:numel(T_license)
    T_license{q}=['% ',T_license{q}];
end

%Add target header
T_license=[{'%% '};{footerTargetText};{'% '};{['% License: <',licenseLink,'>']};T_license(1:end)];

%%

if exist(outputFolder,'dir')==0
   mkdir(outputFolder);
end
addpath(outputFolder);

hw = waitbar(0,'Creating documentation templates...');
for q=1:1:NumberOfLibFiles
   libFileName=fileListLib{q};
   libFileNameFull=fullfile(libFolder,libFileName);
   helpFileName=fullfile(helpFolder,['HELP_',libFileName]);
   helpFileNameNew=fullfile(outputFolder,['HELP_',libFileName]);

   if exist(helpFileName,'file')==0
       funcName=libFileName(1:end-2);
       [T_func]=txtfile2cell(libFileNameFull); %A cell containing the M-file text
       funcLine=T_func{1}; %First line containing function statement
       funcLine = strtrim(funcLine); %First line with leading/trailing spaces removed
       funcLine = regexprep(funcLine,' +',' '); %Remove possible extra spaces
       codeSyntaxSimple=funcLine(10:end);

       %Build M-file text
       T={};
       T{1,1}=['%% ',funcName];
       T{end+1,1}=['% Below is a demonstration of the features of the |',funcName,'| function'];
       T{end+1,1}='';
       T{end+1,1}='%%';
       T{end+1,1}='clear; close all; clc;';
       T{end+1,1}='';
       T{end+1,1}='%% Syntax';
       T{end+1,1}=['% |',codeSyntaxSimple,';|'];
       T{end+1,1}='';
       T{end+1,1}='%% Description ';
       T{end+1,1}='% UNDOCUMENTED ';
       T{end+1,1}='%% Examples ';
       T{end+1,1}='% ';
       T{end+1,1}='%%';
       T{end+1,1}='% ';
       T{end+1,1}='% <<gibbVerySmall.gif>>';
       T{end+1,1}='% ';
       T{end+1,1}='% _*GIBBON*_ ';
       T{end+1,1}='% <www.gibboncode.org>';
       T{end+1,1}='% ';
       T{end+1,1}='% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>';
       T(end+1:end+numel(T_license),1)=T_license;

       cell2txtfile(helpFileNameNew,T,0);

   end
   waitbar(q/NumberOfLibFiles,hw,sprintf('%6.2f',q/NumberOfLibFiles))
end
delete(hw);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%EoM_Single_Pendulum.m

% This script derives the equations of motion for a single pendulum using
% the lagrange equations and the matlab symbolic toolbox.

% Model: Single Pendulum, with a point mass bob. Theta is measured from the
% -j axis. Gravity points in the -j direction. The pendulum is a fixed
% length, and there is no friction.

% The lagrangian (L) is defined as:
%
%       L = T - U
%
% where
%       T = system's kinetic energy
%       U = system's potential energy

% How we go about expressing the equations of motion:
%
%       dL      d  / dL \         * Note that some of those 'd' should be
%       ---  =  -- | -- |           curvy 'd' to represent partial
%       dth     dt \ dw /           derivatives

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                        set up variables                                 %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

syms th w a 'real';
% th == angle from -j axis
% w  == angular rate w = d(th)/dt
% a  == angular acceleration a = d(w)/dt

% Note that in the standard notation we would have:
%    th == "q"  --  generalized coordinate
%    w  == "q dot"
%    a  == "q double dot"

syms m g l 'real';
% m == (point) mass of the pendulum bob
% g == acceleration due to gravity
% l == length of the pendulum

T = (1/2)*m*(l^2)*(w^2);  % kinetic energy of the pendulum
U = m*g*l*(1-cos(th));  % potential energy of the pendulum
L = T - U;  % lagrangian

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                  evaluate partial derivatives                           %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

%              dL
%  dL_dth  ==  ---      Note that 'd' is partial derivative here
%              dth
%
dL_dth = jacobian(L,th);

%              dL
%  dL_dw   ==  ---      Note that 'd' is partial derivative here
%              dw
%
dL_dw = jacobian(L,w);

%               d  / dL \         * Note that some of those 'd' should be
% ddL_dtdw  ==  -- | -- |           curvy 'd' to represent partial
%               dt \ dw /           derivatives
%
% Note the application of the chain rule:
%
ddL_dtdw = jacobian(dL_dw,[th, w]) * [w, a]';

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                  solve equations and write files                        %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Equations of motion:
eqn = ddL_dtdw - dL_dth;

% Solve for acceleration:
soln = solve(eqn,a);

% Write dynamics function:
filename = 'singlePendulumDynamics.m';
matlabFunction(soln,'file',filename,'vars',{th,g,l},'output',{'ddth'});

% Write energy function:
filename = 'singlePendulumEnergy.m';
totalEnergy = T+U;
energy = [totalEnergy;T;U];
matlabFunction(totalEnergy,T,U,'file',filename,'vars',{th,w,m,g,l},'output',{'energy','kinetic','potential'});

author: MatthewPeterKelly

repo: dscTutorials


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% shanChen.m: Multi-component fluid, using a LB method,
%   based on the Shan-Chen model
% [X.Shan and H.Chen, http://dx.doi.org/10.1103/PhysRevE.47.1815].
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample, written in Matlab
% Copyright (C) 2008 Orestis Malaspinas, Andrea Parmigiani, Jonas Latt
% Address: EPFL-STI-LIN Station 9
% E-mail: orestis.malaspinas@epfl.ch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% You should have received a copy of the GNU General Public
% License along with this program; if not, write to the Free
% Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
% Boston, MA  02110-1301, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clf

% GENERAL FLOW CONSTANTS

ly           = 201;
lx           = 201;

G = -1.2;  % Amplitude of the molecular interaction force

omega1 = 1.;  % Relaxation parameter for fluid 1
omega2 = 1.;  % Relaxation parameter for fluid 2

maxT   = 80000;    % total number of iterations
tPlot  = 40;       % iterations between successive graphical outputs

% D2Q9 LATTICE CONSTANTS
tNS   = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cxNS  = [  0,   1,  0, -1,  0,    1,  -1,  -1,   1];
cyNS  = [  0,   0,  1,  0, -1,    1,   1,  -1,  -1];
oppNS = [  1,   4,  5,  2,  3,    8,   9,   6,   7];

[y,x] = meshgrid(1:ly,1:lx);

drho = 0.001;
delta_rho = -drho*(1-2.0*rand(lx));

% INITIAL CONDITION FOR BOTH DISTRIBUTION FUNCTIONS: (T=0) ==> TIn(i) = t(i)
for i=1:9
    fIn(i,1:lx,1:ly) = tNS(i).*(1.0 + delta_rho);
    gIn(i,1:lx,1:ly) = tNS(i).*(1.0 - delta_rho);
end

rho1  = reshape(sum(fIn),lx,ly);
imagesc(rho1');
colorbar
title('Fluid 1 density');
axis equal off; drawnow

% MAIN LOOP (TIME CYCLES)
Gomega1 = G/omega1;
Gomega2 = G/omega2;
for cycle = 1:maxT
    % MACROSCOPIC VARIABLES
    rho1 = sum(fIn);
    rho2 = sum(gIn);
    jx1  = reshape ( (cxNS * reshape(fIn,9,lx*ly)), 1,lx,ly);
    jy1  = reshape ( (cyNS * reshape(fIn,9,lx*ly)), 1,lx,ly);
    jx2  = reshape ( (cxNS * reshape(gIn,9,lx*ly)), 1,lx,ly);
    jy2  = reshape ( (cyNS * reshape(gIn,9,lx*ly)), 1,lx,ly);

    rhoTot_OMEGA = rho1*omega1 + rho2*omega2;
    uTotX = (jx1*omega1+jx2*omega2) ./ rhoTot_OMEGA;
    uTotY = (jy1*omega1+jy2*omega2) ./ rhoTot_OMEGA;

    rhoContrib1x = 0.0;
    rhoContrib2x = 0.0;

    rhoContrib1y = 0.0;
    rhoContrib2y = 0.0;
    for i=2:9
        rhoContrib1x = rhoContrib1x + circshift(rho1*tNS(i), [0,cxNS(i),cyNS(i)])*cxNS(i);
        rhoContrib1y = rhoContrib1y + circshift(rho1*tNS(i), [0,cxNS(i),cyNS(i)])*cyNS(i);

        rhoContrib2x = rhoContrib2x + circshift(rho2*tNS(i), [0,cxNS(i),cyNS(i)])*cxNS(i);
        rhoContrib2y = rhoContrib2y + circshift(rho2*tNS(i), [0,cxNS(i),cyNS(i)])*cyNS(i);
    end

    uTotX1 = uTotX - Gomega1.*rhoContrib2x; %POTENTIAL CONTRIBUTION OF FLUID 2 ON 1
    uTotY1 = uTotY - Gomega1.*rhoContrib2y;

    uTotX2 = uTotX - Gomega2.*rhoContrib1x; %POTENTIAL CONTRIBUTION OF FLUID 2 ON 1
    uTotY2 = uTotY - Gomega2.*rhoContrib1y;

   % COLLISION STEP FLUID 1 AND 2
   for i=1:9
      cuNS1        = 3*(cxNS(i)*uTotX1+cyNS(i)*uTotY1);
      cuNS2        = 3*(cxNS(i)*uTotX2+cyNS(i)*uTotY2);

      fEq(i,:,:)   = rho1 .* tNS(i) .* ...
                       ( 1 + cuNS1 + 0.5*(cuNS1.*cuNS1) - 1.5*(uTotX1.^2+uTotY1.^2) );

      gEq(i,:,:)   = rho2 .* tNS(i) .* ...
                       ( 1 + cuNS2 + 0.5*(cuNS2.*cuNS2) - 1.5*(uTotX2.^2+uTotY2.^2) );

      fOut(i,:,:)  = fIn(i,:,:) - omega1 .* (fIn(i,:,:)-fEq(i,:,:));
      gOut(i,:,:)  = gIn(i,:,:) - omega2 .* (gIn(i,:,:)-gEq(i,:,:));
   end

   % STREAMING STEP FLUID 1 AND 2
   for i=1:9
      fIn(i,:,:) = circshift(fOut(i,:,:), [0,cxNS(i),cyNS(i)]);
      gIn(i,:,:) = circshift(gOut(i,:,:), [0,cxNS(i),cyNS(i)]);
   end

   % VISUALIZATION
   if(mod(cycle,tPlot)==0)
       rho1     = reshape(rho1,lx,ly);
       imagesc(rho1'); colorbar
       title('Fluid 1 density');
       axis equal off; drawnow
   end
end

repo: Aero-matlab


text:
function S = sum(disc)
%SUM    Integration functional for TRIGCOLLOC.

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

[ignored, S] = functionPoints(disc);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% test_all_emission.m

list = {
	'eml_curvature test'
	'eml_osem_example'
	'eml_em_test'
	'eml_osem_test'
	'eml_sps_os_test'
	'eml_psca_test'
	'eql_os_emdp_test'
	'eql_sps_os_test'
	'hole_example1'
	'psf_mismatch_example1'
	'psf_mismatch_example2'
};

run_mfile_local(list)
%run_mfile_local(list, 'pause', 1)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function planC = findAndSetMinCTSpacing(planC, minSpacing, maxSpacing, alternateSpacing, scanNumV)
%"findAndSetMinCTSpacing"
%   Finds the minimum spacing in a CT scan and records the superior and inferior slice numbers.
%   minSpacing should be the same as optS.smallestUniformCTSliceSpacing.
%
%VHC 23 Aug 02
%
%Latest modifications:
% 25 Sept 02, JOD.
% 19 Feb 03, JOD, Change the specification of smallest CT such that a range is given:
%                 optS.lowerLimitUniformCTSliceSpacing to optS.upperLimitUniformCTSliceSpacing.
%                 Given a set of CT values, we choose the largest block which CT slice spacing which
%                 falls within these limits.  Otherwise we create a set with a spacing equal to
%                 optS.alternateUniformCTSliceSpacing.  Rewritten using blockS structure.
% 29 Apr 03, JOD, If slice thicknesses are not available, they are assigned by using the zValues.
%  8 May 03, JOD, only deduce slice thicknesses if thicknesses are not already present.
%  9 May O3, JOD, corrected faulty indexing in determining slice block to keep.
% 20 May 03, ES,  changed < to <= in block recognition loop.
% 20 May 03, JOD, corrected bug in assigning slice thicknesses to blocks of slices.
%
%Usage:
%   function planC = findAndSetMinCTSpacing(planC, minSpacing, maxSpacing, alternateSpacing)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

accuracy = 0.001; %this value seems to work best for CT scans.  This is how much different the slice spacing
                     %can be and still be considered as approximately the same slice spacing.

%get scan info
indexS = planC{end};

% Get scan indices
if ~exist('scanNumV','var')
    scanNumV = 1:length(planC{indexS.scan});
end

for scanNum = scanNumV

	thicknessV = deduceSliceWidths(planC,scanNum);
	CERRStatusString('Using zValues to compute voxel thicknesses.')

    for i = 1 : length(thicknessV)  %put back into planC
        planC{indexS.scan}(scanNum).scanInfo(i).sliceThickness = thicknessV(i);
    end

	%Get all blocks of the CT scan which have constant slice spacing
	%and store their characteristics.
	slice = 1;
	blockNum = 1;
    blockS = struct('');
	while slice <= length(thicknessV) - 1

        start = slice;
        finish  = slice;

        same = 1;

      go = 0;
      while same && (slice <= length(thicknessV) - 1)      %changed < to <=, ES.
        go = 1;
        if (abs(thicknessV(slice + 1) - thicknessV(slice)) < accuracy)
          slice = slice + 1;
          finish = slice;
        else
          slice = slice + 1;
          same = 0;
        end
      end

      %store block information
      blockS(blockNum).start = start;
      blockS(blockNum).finish = finish;
      blockS(blockNum).width = finish - start + 1;
      blockS(blockNum).spacing = thicknessV(start);
      blockNum = blockNum + 1;
      if go == 0
        slice = slice + 1;
      end

	end

	numSlices = length(thicknessV);

	%Choose which blocks have thicknesses fulfilling the criteria:

	okV = [];
	for i = 1 : length(blockS)
      if blockS(i).spacing <= maxSpacing &&  blockS(i).spacing >= minSpacing
        okV = [okV, i];
      end
	end

	if ~isempty(okV)
      %Of those, which has the largest number of consecutive slices?
      widthV = [blockS(okV).width];

      indV = find(max(widthV) == widthV);
      ind = indV(1);

      ind2 = okV(ind);

      sliceThickness = blockS(ind2).spacing;
      sliceNumSup = blockS(ind2).start;
      sliceNumInf = blockS(ind2).finish;

	else

      sliceThickness = alternateSpacing;
      sliceNumMiddle = floor(numSlices/2);
      sliceNumSup = sliceNumMiddle;
      sliceNumInf = sliceNumMiddle;

	end

	uniformScanInfo = planC{indexS.scan}(scanNum).scanInfo(1);
	uniformScanInfo.sliceNumSup = sliceNumSup;
	uniformScanInfo.sliceNumInf = sliceNumInf;
	uniformScanInfo.sliceThickness = sliceThickness;
	uniformScanInfo.supInfScansCreated = 0; %scans not created yet, but the field is initialized
    scanArray = getScanArray(planC{indexS.scan}(scanNum));
	uniformScanInfo.minCTValue = double(min(scanArray(:)));
	uniformScanInfo.maxCTValue = double(max(scanArray(:)));

	optS = planC{indexS.CERROptions};

    if ~isfield(optS, 'uniformizedDataType')
        try
            if planC{indexS.scan}(scanNum).uniformScanInfo.bytesPerPixel == 1
                optS.uniformizedDataType = 'uint8';
            else
                optS.uniformizedDataType = 'uint16';
            end
        catch
            optFromFile = opts4Exe([getCERRPath,'CERROptions.json']);
            optS.uniformizedDataType = optFromFile.uniformizedDataType;
        end
        planC{indexS.CERROptions} = optS;
    end

	switch optS.uniformizedDataType

      case 'uint8'

        uniformScanInfo.bytesPerPixel = 1;

      case 'uint16'

        uniformScanInfo.bytesPerPixel = 2;

      otherwise

        error('Error determining CT uniformized data type.')

	end

	%Remove unneeded holdover fields copied from scanInfo(1).
	if isfield(uniformScanInfo, 'imageNumber')
        uniformScanInfo = rmfield(uniformScanInfo, 'imageNumber');
	end
	if isfield(uniformScanInfo, 'zValue')
        uniformScanInfo = rmfield(uniformScanInfo, 'zValue');
	end
	if isfield(uniformScanInfo, 'voxelThickness')
        uniformScanInfo = rmfield(uniformScanInfo, 'voxelThickness');
	end

	planC{indexS.scan}(scanNum).uniformScanInfo = uniformScanInfo;

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [model]  = setRxnSubSystems(model,reactions,subSystems)
% Sets the subSystems of the specified reactions to the specified
% subSystems
%
% USAGE:
%
%    [model]  = addSubSystemToReaction(model,reaction,subSystem)
%
% INPUT:
%    model:                 A COBRA model struct with at least rxns and
%                           subSystems fields
%    reactions:             Either a string identifying a reaction, or a
%                           cell array of strings identifying multiple
%                           reactions, or a double vector or boolean vector
%                           identifying positions.
%    subSystems:            A String identifying a subsystem, or a cell
%                           array of strings to assign multiple subSystems.
%                           If multiple strings are present, an empty
%                           string will be removed
%
% OUTPUT:
%    model:                 The model with the subSystems set to the
%                           respective reactions (old values will be removed).
%
% EXAMPLE:
%    % Set TCA as subSystem for a set of reactions
%    [model]  = setRxnSubSystems(model, {'ACONTa';'ACONTb';'AKGDH';'CS';'FUM';'ICDHyr';'MDH';'SUCOAS'}, 'TCA')
%    % Set Glycolysis/Gluconeogenesis to reaction 18 of the model
%    [model]  = setRxnSubSystems(model, 18, {'Glycolysis','Gluconeogenesis'})
%
% .. Author: - Thomas Pfau Nov 2017

if ischar(reactions) || iscell(reactions)
    [tempreactions] = ismember(model.rxns,reactions);
    if (iscell(reactions) && (sum(tempreactions) ~= length(reactions))) || (ischar(reactions) && sum(tempreactions) ~= 1)
        %If there are an unequal amount of reactions
        warning('There were reaction IDs not present in the model');
    end
    reactions = tempreactions;
end

if ~isfield(model,'subSystems')
    model.subSystems = cell(numel(model.rxns),1);
    model.subSystems(:) = {{''}};
end

if ischar(subSystems)
    subSystems = {subSystems};
end

if numel(subSystems) > 1
    subSystems = setdiff(subSystems,{''});
end

model.subSystems(reactions) = {subSystems};

author: opencobra

repo: cobratoolbox


text:
function varargout=paramsetch(varargin)
%
% For calling details please see v_paramsetch.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_paramsetch.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: paramsetch.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_paramsetch(varargin{:});
else
    v_paramsetch(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function sep = graph_separated(G, X, Y, S)

G2 = G;
G2(S,:) = 0;
G2(:,S) = 0;
conn = reachability_graph(G2);
conn2 = conn(X,Y);
sep = all(conn2(:)==0);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function value = gamma_inc ( p, x )

%*****************************************************************************80
%
%% GAMMA_INC computes the incomplete Gamma function.
%
%  Discussion:
%
%    GAMMA_INC(P,X) = Integral ( 0 <= T <= X ) T**(P-1) EXP(-T) DT / GAMMA(P).
%
%    GAMMA_INC(P,       0) = 0,
%    GAMMA_INC(P,Infinity) = 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 September 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    B L Shea,
%    Chi-squared and Incomplete Gamma Integral,
%    Algorithm AS239,
%    Applied Statistics,
%    Volume 37, Number 3, 1988, pages 466-473.
%
%  Parameters:
%
%    Input, real P, the exponent parameter.
%    0.0 < P.
%
%    Input, real X, the integral limit parameter.
%    If X is less than or equal to 0, GAMMA_INC is returned as 0.
%
%    Output, real VALUE, the value of the function.
%
  exp_arg_min = -88.0;
  overflow = 1.0E+37;
  plimit = 1000.0;
  tol = 1.0E-07;
  xbig = 1.0E+08;

  value = 0.0;

  if ( p <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GAMMA_INC - Fatal error!\n' );
    fprintf ( 1, '  Parameter P <= 0.\n' );
    error ( 'GAMMA_INC - Fatal error!' );
  end

  if ( x <= 0.0 )
    value = 0.0;
    return
  end
%
%  Use a normal approximation if PLIMIT < P.
%
  if ( plimit < p )
    pn1 = 3.0 * sqrt ( p ) * ( ( x / p )^( 1.0 / 3.0 ) + 1.0 / ( 9.0 * p ) - 1.0 );
    cdf = normal_01_cdf ( pn1 );
    value = cdf;
    return
  end
%
%  Is X extremely large compared to P?
%
  if ( xbig < x )
    value = 1.0;
    return
  end
%
%  Use Pearson's series expansion.
%  (P is not large enough to force overflow in the log of Gamma.
%
  if ( x <= 1.0 | x < p )

    arg = p * log ( x ) - x - gammaln ( p + 1.0 );
    c = 1.0;
    value = 1.0;
    a = p;

    while ( 1 )

      a = a + 1.0;
      c = c * x / a;
      value = value + c;

      if ( c <= tol )
        break
      end

    end

    arg = arg + log ( value );

    if ( exp_arg_min <= arg )
      value = exp ( arg );
    else
      value = 0.0;
    end

  else
%
%  Use a continued fraction expansion.
%
    arg = p * log ( x ) - x - gammaln ( p );
    a = 1.0 - p;
    b = a + x + 1.0;
    c = 0.0;
    pn1 = 1.0;
    pn2 = x;
    pn3 = x + 1.0;
    pn4 = x * b;
    value = pn3 / pn4;

    while ( 1 )

      a = a + 1.0;
      b = b + 2.0;
      c = c + 1.0;
      pn5 = b * pn3 - a * c * pn1;
      pn6 = b * pn4 - a * c * pn2;

      if ( 0.0 < abs ( pn6 ) )

        rn = pn5 / pn6;

        if ( abs ( value - rn ) <= min ( tol, tol * rn ) )

          arg = arg + log ( value );

          if ( exp_arg_min <= arg )
            value = 1.0 - exp ( arg );
          else
            value = 1.0;
          end

          return

        end

        value = rn;

      end

      pn1 = pn3;
      pn2 = pn4;
      pn3 = pn5;
      pn4 = pn6;
%
%  Rescale terms in continued fraction if terms are large.
%
      if ( overflow <= abs ( pn5 ) )
        pn1 = pn1 / overflow;
        pn2 = pn2 / overflow;
        pn3 = pn3 / overflow;
        pn4 = pn4 / overflow;
      end

    end

  end

  return
end


text:
% NOTE: DO NOT REMOVE THIS LINE XMAKEFILE_TOOL_CHAIN_CONFIGURATION
function toolChainConfiguration = arduinouno()
%ARDUINOUNO Defines a tool chain configuration

% Copyright 2012 The MathWorks, Inc.

try
hSetup = realtime.setup.PackageInfo;
    spPkg = hSetup.getSpPkgInfo('Arduino Uno');
    arduinoRootDir = hSetup.getTpPkgRootDir('ARDUINO', spPkg);
catch me
    toolChainConfiguration = {};
    return;
end

% General
toolChainConfiguration.Configuration = 'Arduinouno';
toolChainConfiguration.Version = '2.0';
toolChainConfiguration.Description = 'MS Visual Studio';
toolChainConfiguration.Operational = true;
toolChainConfiguration.Decorator = 'linkfoundation.xmakefile.decorator.eclipseDecorator';
% Make
toolChainConfiguration.MakePath = @(src) fullfile(arduinoRootDir, 'hardware', 'tools', 'avr', 'utils', 'bin', 'make');
toolChainConfiguration.MakeFlags = '-f "[|||MW_XMK_GENERATED_FILE_NAME[R]|||]" [|||MW_XMK_ACTIVE_BUILD_ACTION_REF|||]';
toolChainConfiguration.MakeInclude = '';
% Compiler
toolChainConfiguration.CompilerPath = @(src) fullfile(arduinoRootDir, 'hardware', 'tools', 'avr', 'bin', 'avr-gcc');
toolChainConfiguration.CompilerFlags = '-c -x none';
toolChainConfiguration.SourceExtensions = '.c,.cpp';
toolChainConfiguration.HeaderExtensions = '.h';
toolChainConfiguration.ObjectExtension = '.o';
% Linker
toolChainConfiguration.LinkerPath = @(src) fullfile(arduinoRootDir, 'hardware', 'tools', 'avr', 'bin', 'avr-gcc');
toolChainConfiguration.LinkerFlags = '-o [|||MW_XMK_GENERATED_TARGET_REF|||]';
toolChainConfiguration.LibraryExtensions = '.lib,.a';
toolChainConfiguration.TargetExtension = '.elf';
toolChainConfiguration.TargetNamePrefix = '';
toolChainConfiguration.TargetNamePostfix = '';
% Archiver
toolChainConfiguration.ArchiverPath = @(src) fullfile(arduinoRootDir, 'hardware', 'tools', 'avr', 'bin', 'avr-ar');
toolChainConfiguration.ArchiverFlags = '-crs $(TARGET_FILE)';
toolChainConfiguration.ArchiveExtension = '.a';
toolChainConfiguration.ArchiveNamePrefix = 'lib';
toolChainConfiguration.ArchiveNamePostfix = '';
% Pre-build
toolChainConfiguration.PrebuildEnable = false;
toolChainConfiguration.PrebuildToolPath = '';
toolChainConfiguration.PrebuildFlags = '';
% Post-build
toolChainConfiguration.PostbuildEnable = true;
toolChainConfiguration.PostbuildToolPath = @(src) fullfile(arduinoRootDir, 'hardware', 'tools', 'avr', 'bin', 'avr-objcopy');
toolChainConfiguration.PostbuildFlags = '-O ihex -R .eeprom $(TARGET)  [|||MW_XMK_MODEL_NAME|||].hex';
% Execute
toolChainConfiguration.ExecuteDefault = false;
toolChainConfiguration.ExecuteToolPath = '';
toolChainConfiguration.ExecuteFlags = '';
% Directories
toolChainConfiguration.DerivedPath = '[|||MW_XMK_SOURCE_PATH_REF|||]';
toolChainConfiguration.OutputPath = '';
% Custom
toolChainConfiguration.Custom1 = '';
toolChainConfiguration.Custom2 = '';
toolChainConfiguration.Custom3 = '';
toolChainConfiguration.Custom4 = '';
toolChainConfiguration.Custom5 = '';
end

% LocalWords:  XMAKEFILE Uno linkfoundation xmakefile avr utils XMK crs ihex
% LocalWords:  eeprom

author: Sable

repo: mcbench-benchmarks


text:
function [xh,pf, P, L] = RBPF(x_true, y, LRS, u, a, pf, Nsamples, resampling_strategy)

    % x_true: true position of robot
    % y: True laser scan
    % LRS: laser scan inforamtion
    % u: Movement control
    % a: Odometry parameters found via Motion_Model_Velocity_Test
    % pf: pf class with all particle filter information
    % Nsamples: number of samples for the Particle Filter
    % resampling_startegy: DO NOT CHANGE

    %% Variable Initialization

    k = pf.k;
    NPC = pf.NPC;
    q_past = pf.q(:,k-1);
    x_past= pf.particles(:,:,k-1);

    % Initialise  new variables
    x_initial_estimate=zeros(3, NPC);           % Initial estimate from previous Xt-1 and Ut
    yh=zeros(2, 361, NPC);                      % Laser scan estimate from x_initial_estimate
    x_new_sm=zeros(3, NPC);                     % New position after scan matching (ICP)
    x_sampled=zeros(3,Nsamples, NPC);           % Smaple position around every particle position
    p_xt=zeros(Nsamples, NPC);                  % Density function for odometry (Motion_Model_velocity)
    p_zt=zeros(Nsamples, NPC);                  % Density function for scans (Scan Matching)
    tau_xt=zeros(Nsamples, NPC);                % Optimal Proposal Distribution
    n=zeros(1,NPC);                             % Sample normalizer
    mu=zeros(3, NPC);                           % Mean for Gaussian distribution from all particle samples
    cov=zeros(1, NPC);                          % Covariance for Gaussian distribution from all particle samples
    q=zeros(size(q_past));                      % Particle weights
    x_final=zeros(3, NPC);                      % final estimate particle position after RBPF

    %% RBPF for every particle
    for i = 1:NPC

        %Estimate particle position from x_(t-1) after Ut
        x_initial_estimate(:,i) = Kinematics(x_past(:,i),u+sqrt(pf.R1)*randn(length(u),1));

        if  x_initial_estimate(3,i)>2*pi
            x_initial_estimate(3,i)=x_initial_estimate(3,i)-2*pi;
        elseif x_initial_estimate(3,i)<-2*pi
            x_initial_estimate(3,i)=x_initial_estimate(3,i)+2*pi;
        end

        %Scan Matcher returns estimate of laser given an initial particle position estimate
        yh(:,:,i) = Measurement_Estimate_From_Grid(x_initial_estimate(:,i),pf.P(:,:, i),LRS.Resolution,LRS.FoV, LRS.MaxDistance, pf.UsableArea);

        % Scan pre-processing prior to ICP algortithm
            y_help=y;
            yh_help=yh(:,:,i);

            % Use data inside Usable Area only
            y_help(y_help>pf.UsableArea)=nan;
            yh_help(yh_help>pf.UsableArea)=nan;

            % Chose data that whose angle is usable in both scans
            find_y_help=isnan(y_help(2,:));
            find_yh_help=isnan(yh_help(2,:));
            del_cols=or(find_y_help, find_yh_help);
            y_help(:,(del_cols(1,:)==1))=[];
            yh_help(:,(del_cols(1,:)==1))=[];

            % Transform to Polar coordinates and World frame
            y_help=Polar2Cart(y_help);
            yh_help=Polar2Cart(yh_help);
            y_help=Rotate_Data(y_help, x_true);
            yh_help=Robot2World([x_initial_estimate(1,i); x_initial_estimate(2,i); x_initial_estimate(3,i)], yh_help);

            % Add the particles position so that its also transformed
            y_help=horzcat(y_help, x_true(1:2));
            yh_help=horzcat(yh_help, x_initial_estimate(1:2,i));

        failure=false;                               % Initialse ICP failure to FALSE

        try
        [TR,TT, data] = ICP(y_help, yh_help);        % ICP algortithm returns the translational and rotational matrices and the resulting data
        catch
            failure=true;
        end
        % Check if ICP failed
        if abs(TT(1))>50  || abs(TT(2))>50
            failure=true;
        end

        if failure
            disp('ICP failed')
            x_final(:,i)=x_initial_estimate(:,i);    % If ICP failed used last known position estimate
            q(i)=q_past(i);                          % do not change particle weights
        else

            theta_new_sm=x_initial_estimate(3,i)+asin(TR(2,1));     % Modify theta angle according to the rotation given by ICP
            x_new_sm(:,i)=[data(1:2,end); theta_new_sm];            % Use rotated+translated initila_estimate as new (and better) positon

            %Plot grid-map, all positions and all scans (Black: True postion, Blue: Position and scan initial estimates, red:
            %Converged position and laser data
                figure(42)
                set(gcf,'units','normalized','outerposition',[1 0 1 1]);
                clf;
                caxis([0, 1]);
                colormap gray;
                colormap(flipud(colormap));
                hold on;
                imagesc(flipud(pf.P(:,:, i)));

                % Plot scans
                plot(y_help(1,:), y_help(2,:), '.g')
                plot(yh_help(1,:), yh_help(2,:), '.b')
                plot(data(1,:), data(2,:), '.r')

                % Plot positions
                plot(x_true(1), x_true(2), 'dk')
                plot(x_initial_estimate(1,i), x_initial_estimate(2, i), 'db')
                plot(x_new_sm(1,i), x_new_sm(2,i), 'dr')
                hold off;

                % Save each iteration map into as png
                % str=sprintf('path_images/fig%d.png', pf.k);
                % saveas(gcf, str)

            % Randomize samples around good estimate (x_new_sm) of particle position
            theta_samples=x_new_sm(3,i);
            if theta_samples > 2*pi
                theta_samples=theta_samples-2*pi;
            elseif theta_samples < -2*pi
                theta_samples=theta_samples-2*pi;
            end

            x_sampled(:,:, i) = [ normrnd([x_new_sm(1,i).*ones(1, Nsamples); x_new_sm(2,i).*ones(1, Nsamples)],0.5, [2 Nsamples]); theta_samples.*ones(1, Nsamples)];

            for h=1:Nsamples

                % Calculate odometry probability for that sample
                p_xt(h,i)= Motion_Model_Velocity(x_sampled(:,h,i), u, x_past(:,i), a);

                % Calculate scan matching probability for that sample
                ytru=y;
                yest = Measurement_Estimate_From_Grid(x_sampled(:,h,i),pf.P(:,:, i),LRS.Resolution,LRS.FoV, LRS.MaxDistance, pf.UsableArea);
                yest(yest>pf.UsableArea)=nan;
                find_yest=isnan(yest(2,:));
                yest(:,(find_yest(1,:)==1))=[];
                ytru(:,(find_yest(1,:)==1))=[];
                p_zt(h,i) = Measurement_Model(ytru,yest,[1 1]);

                % Calculate optimal proposal distribution
                tau_xt(h,i)=abs(p_xt(h,i)*p_zt(h,i));

            end

            n(i)=sum(tau_xt(:,i));                                  % Calculate sample normalizer
            mu(1, i)=(1/n(i))*(x_sampled(1,:,i)*tau_xt(:,i));       % Calculate mu for gaussian distribution from samples
            mu(2, i)=(1/n(i))*(x_sampled(2,:,i)*tau_xt(:,i));
            mu(3, i)=theta_samples;                                 % Orientation of samples is constant (the one from ICP)

            % Calculate covariance for gaussian distribution from samples
            cov_sum=0;
            for h=1:Nsamples
                cov_sum=cov_sum+(x_sampled(1:2,h,i)-mu(1:2,i))'*(x_sampled(1:2,h,i)-mu(1:2,i))*tau_xt(h,i);
            end
            cov(i)=(1/n(i))*cov_sum;

            %Calculate new weights
            q(i)=q_past(i)*n(i);

            %Sample gaussian for final position of particle
            x_final(:,i)=[normrnd([mu(1,i); mu(2,i)],cov(i), [2 1]); mu(3, i)];
        end

        %Update GridMap for taht particle
        [pf.P(:,:,i), pf.L(:,:,i)] = Occupancy_Grid_Mapping(pf.L(:,:,i), x_final(:,i), y, pf.UsableArea, pf.L0, pf.gridsize, LRS);

    end

    q=q/norm(q,1);                      % Normalize weights
    Neff = 1/sum(q.^2);                 % Calcuklate Neff for possible resampling

    %% Resampling step

    resample_percentage = 0.3;          % 0.5 is the best resample percentage proposed by Cyrill Stachniss
    Nt = resample_percentage*NPC;
    if Neff < Nt
        disp('Resampling');
        [x_final(:,:), q, idx] = Resample(x_final(:,:), q, resampling_strategy);
        P_save_resamp=pf.P(:,:,:);
        L_save_resamp=pf.L(:,:,:);
        for n=1:NPC
           pf.P(:,:,n)=P_save_resamp(:,:,idx(n));
           pf.L(:,:,n)=L_save_resamp(:,:,idx(n));
        end
    end

    pf.q(:,k) = q;                      % Save particle weights (might have changed during resampling)

    %% Calculate estimate of true position (xh) given all weighted particles
    xh = zeros(3,1);
    for i=1:NPC
            xh = xh +(q(i).*x_final(:,i));
    end
    pf.xh(:, k)= xh;

    %% Store final position of all particles
    pf.particles(:,:,k) = x_final(:,:);

    % Make sure no errors have occured. If so probably becasue of p_xt or p_zt being equal to 0
    if isnan(x_final(1,i)) || isnan(x_final(2,i))
                1; %debugger
                x_final=x_initial_estimate;
    end

    %% Visualize results: Plot Particles
%         figure(89)
%         set(gcf,'units','normalized','outerposition',[-1 0 1 1]);
%         clf;
%         hold on;
%         for g = 1:size(map,2)
%             plot([map(1,g) map(3,g)],[map(2,g) map(4,g)],'k')
%         end
%         plot(squeeze(pf.particles(1,:,pf.k-1))',squeeze(pf.particles(2,:,pf.k-1))','bo');
%         plot(squeeze(pf.particles(1,:,pf.k))',squeeze(pf.particles(2,:,pf.k))','k+');
%         plot(squeeze(pf.xh(1,pf.k))',squeeze(pf.xh(2,pf.k))','ro');
%         plot(x_true(1), x_true(2), 'go')
%         hold off;

end

author: Allopart

repo: rbpf-gmapping


text:
function r = subsasgn(r,s,b)
%SUBSASGN     Implements subscripted assignments for hessians
%
%  example  r(2,:) = b
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 05/09/07     S.M. Rump  assignment r(:)=...
% modified 08/26/12     S.M. Rump  global variables removed
% modified 09/26/12     S.M. Rump  index handling (thanks to Matthew Weinstein)
%

  N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');

  if length(s)>1
    error('multiple indexing for hessian assignment not allowed')
  end

  if strcmp(s.type,'()')     % assignment r(i) = b

    rEmpty = isempty(r);     % assignment r(i) = b for empty r
    if isempty(b)
      % does not work in Matlab 5.3 for sparse r
      r.x(s.subs{:}) = [];
      value = zeros(size(r.x));
      value(s.subs{:}) = 1;
      index = find(value);
      r.dx( : , index ) = [];
      r.hx( : , index ) = [];
      return
    end

    if ~isa(b,'hessian')
      b = hessian(b);
    end

    if ~rEmpty
      resultIsintval = isa(r.x,'intval');
      if ~resultIsintval & isa(b.x,'intval')
        r = intval(r);
        resultIsintval = 1;
      end
      sizeincreased = 0;
      if length(s.subs)==1              % single index
        if ~isequal(s.subs{1},':')      % not r(:)=...
          if size(r.x,1)==1             % row vector
            sizeincreased = ( s.subs{1}>size(r.x,2) );
          else
            sizeincreased = ( s.subs{1} > prod(size(r.x)) );
          end
          if sizeincreased
            srx = size(r.x);
            if length(srx)==2
              if all(prod(srx)~=srx)
                error('matrix cannot be resized by assignment a(I) = b')
              end
            else
              error('attempt to grow size of array along ambiguous dimension')
            end
          end
        end
      else                            % multiple index
        for i=1:length(s.subs)
          if ~isequal(s.subs{i},':')
            sizeincreased = sizeincreased | ( s.subs{i} > size(r.x,i) );
          end
        end
      end
      if sizeincreased                % size increased, adapt .x and .dx and .hx
        rx = r.x;
        rdx = r.dx;
        rhx = r.hx;
        value = ones(size(r.x));
        value( s.subs{:} ) = 0;
        r.x = zeros(size(value));
        INTLAB_HESSIAN_SPARSE = getappdata(0,'INTLAB_HESSIAN_SPARSE');
        if N<INTLAB_HESSIAN_SPARSE
          r.dx = zeros(N,prod(size(value)));
          r.hx = zeros(N*N,prod(size(value)));
        else
          r.dx = sparse([],[],[],N,prod(size(value)));
          r.hx = sparse([],[],[],N*N,prod(size(value)));
      end
        if resultIsintval
          r.x = intval(r.x);
          r.dx = intval(r.dx);
          r.hx = intval(r.hx);
        end
        index = find(value);
        r.x( index ) = rx;
        r.dx( : , index ) = rdx;
        r.hx( : , index ) = rhx;
      end
    else                     % assignment r(i) = b for empty r
      resultIsintval = isa(b.x,'intval');
      r.x(s.subs{:}) = 0;
      INTLAB_HESSIAN_SPARSE = getappdata(0,'INTLAB_HESSIAN_SPARSE');
      if N<INTLAB_HESSIAN_SPARSE
        r.dx = zeros(N,prod(size(r.x)));
        r.hx = zeros(N*N,prod(size(r.x)));
      else
        r.dx = sparse([],[],[],N,prod(size(r.x)));
        r.hx = sparse([],[],[],N*N,prod(size(r.x)));
      end
      if resultIsintval
        r.x = intval(r.x);
        r.dx = intval(r.dx);
        r.hx = intval(r.hx);
      end
    end
    r.x( s.subs{:} ) = b.x;
    value = reshape(1:prod(size(r.x)),size(r.x));
    index = value( s.subs{:} );
    if ( prod(size(b.x))==1 ) & ( length(index)~=1 )
      r.dx(:,index) = b.dx(:,ones(1,length(index)));
      r.hx(:,index) = b.hx(:,ones(1,length(index)));
    else
      r.dx(:,index) = b.dx;
      r.hx(:,index) = b.hx;
    end
    if rEmpty
      r = class(r,'hessian');
    end
  else
    error('invalid index reference for hessian')
  end

  % avoid Matlab 6.5f bug:
  % a = sparse([],[],[],1,1); a = reshape(a,1,1); abs(a)
  % produces  9.6721e-317  or similar number in underflow range
  if prod(size(r.x))==1
    r.x = full(r.x);
  end
  if prod(size(r.dx))==1
    r.dx = full(r.dx);
  end
  if prod(size(r.hx))==1
    r.hx = full(r.hx);
  end

repo: OpenMAS


text:
function [coherogram,phase,t,f] = bz_MTCoherogram(lfp1,lfp2,varargin)

%MTCoherogram - Compute LFP coherogram by multi-taper estimation.
%
%  USAGE
%
%    [coherogram,phase,t,f] = MTCoherogram(lfp1,lfp2,<options>)
%
%    lfp1,lfp2      wide-band LFPs (one channel each).
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'frequency'   sampling rate (in Hz) (default = from timestamps if
%                   available, otherwise 1250Hz)
%     'range'       frequency range (in Hz) (default = all)
%     'window'      duration (in s) of the time window (default = 5)
%     'overlap'     overlap between successive windows (default = window/2)
%     'step'        step between successive windows (default = window/2)
%     'tapers'      relative resolution and order of the tapers [NW K]
%                   (default = [3 5])
%     'pad'         FFT padding (see help for cohgramc) (default = 0)
%     'show'        plot results (default = 'off')
%     'cutoffs'     cutoff values for color plot (default = [0 1])
%    =========================================================================
%
%  NOTES
%
%    The LFP can be provided either as a time stamped matrix (list of time-voltage
%    pairs), or as a voltage vector - in which case the frequency must be specified.
%
%    The time displacement between successive short time coherences can be supplied
%    either as a 'step' (explicit time difference) or as an 'overlap' (between
%    successive time windows).
%
%  OUTPUT
%
%    coherogram     coherogram magnitude
%    phase          coherogram phase
%    t              time bins
%    f              frequency bins
%
%  DEPENDENCIES
%
%    This function requires the chronux toolbox.
%
%  SEE
%
%    See also MTCoherence, MTSpectrum, MTSpectrogram, PlotColorMap.

% Copyright (C) 2010-2014 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Make sure chronux is installed and functional
CheckChronux('cohgramc');

% Defaults
f = 1250;
frequency = [];
window = 5;
range = [];
overlap = [];
step = [];
show = 'off';
tapers = [3 5];
pad = 0;
cutoffs = [0 1];

% Check number of parameters
if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help MTCoherogram'' for details).');
end

% Check parameter sizes
if size(lfp1,2) ~= 1 && size(lfp1,2) ~= 2,
	error('Parameter ''lfp1'' is not a vector or a Nx2 matrix (type ''help MTCoherogram'' for details).');
end
if size(lfp2,2) ~= 1 && size(lfp2,2) ~= 2,
	error('Parameter ''lfp2'' is not a vector or a Nx2 matrix (type ''help MTCoherogram'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help MTCoherogram'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'frequency',
			frequency = varargin{i+1};
			if ~isdscalar(frequency,'>0'),
				error('Incorrect value for property ''frequency'' (type ''help MTCoherogram'' for details).');
			end
		case 'range',
			range = varargin{i+1};
			if ~isdvector(range,'#2','<','>=0'),
				error('Incorrect value for property ''range'' (type ''help MTCoherogram'' for details).');
			end
		case 'window',
			window = varargin{i+1};
			if ~isdscalar(window,'>0'),
				error('Incorrect value for property ''window'' (type ''help MTCoherogram'' for details).');
			end
		case 'overlap',
			overlap = varargin{i+1};
			if ~isdscalar(overlap,'>0'),
				error('Incorrect value for property ''overlap'' (type ''help MTCoherogram'' for details).');
			end
		case 'step',
			step = varargin{i+1};
			if ~isdscalar(step,'>0'),
				error('Incorrect value for property ''step'' (type ''help MTCoherogram'' for details).');
			end
		case 'tapers',
			tapers = varargin{i+1};
			if ~isivector(tapers,'#2','>0'),
				error('Incorrect value for property ''tapers'' (type ''help MTCoherogram'' for details).');
			end
		case 'pad',
			pad = varargin{i+1};
			if ~isiscalar(pad,'>-1'),
				error('Incorrect value for property ''pad'' (type ''help MTCoherogram'' for details).');
			end
		case 'show',
			show = varargin{i+1};
			if ~isstring_FMAT(show,'on','off'),
				error('Incorrect value for property ''show'' (type ''help MTCoherogram'' for details).');
			end
		case 'cutoffs',
			cutoffs = varargin{i+1};
			if ~isdvector(cutoffs,'#2','>=0','<'),
				error('Incorrect value for property ''cutoffs'' (type ''help MTCoherogram'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help MTCoherogram'' for details).']);
	end
end

% Determine LFP frequency
if isempty(frequency),
	if size(lfp1,2) == 2,
		frequency = 1/median(diff(lfp1(:,1)));
	else
		frequency = f;
	end
end

% Determine step/overlap
if isempty(step),
	if isempty(overlap),
		overlap = window/2;
	end
else
	if isempty(overlap),
		overlap = window-step;
	elseif overlap ~= window-step,
		error('Incompatible ''step'' and ''overlap'' parameters (type ''help MTCoherogram'' for details).');
	end
end

% Compute and plot coherogram
parameters.Fs = frequency;
if ~isempty(range), parameters.fpass = range; end
parameters.tapers = tapers;
parameters.pad = pad;
[coherogram,phase,~,~,~,t,f] = cohgramc(lfp1,lfp2,[window window-overlap],parameters);
% t = t'+lfp1(1,1);
f = f';
coherogram = coherogram';
% coherogram = permute(coherogram,[2 1 3]);  % Previous code by Gabrielle Girardeau, keep it around just in case
phase = phase';
if strcmp(lower(show),'on'),
  figure;hold on;
  subplot(2,1,1);
  PlotColorMap(coherogram,'x',t,'y',f,'cutoffs',cutoffs,'newfig','off');
  xlabel('Time (s)');
  ylabel('Frequency (Hz)');
  title('Coherogram Amplitude');
  subplot(2,1,2);
  PlotColorMap(phase,'x',t,'y',f,'cutoffs',[-pi pi],'newfig','off');
  xlabel('Time (s)');
  ylabel('Frequency (Hz)');
  title('Coherogram Phase');
end


text:
function [params, names] = rbfard2KernExtractParam(kern)

% RBFARD2KERNEXTRACTPARAM Extract parameters from the RBFARD2 kernel structure.
% FORMAT
% DESC Extract parameters from the automatic relevance determination
% radial basis function kernel structure into a vector of parameters for
% optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% DESC Extract parameters and parameter names from the automatic
% relevance determination radial basis function kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings containg the parameter names.
%
% SEEALSO rbfard2KernParamInit, rbfard2KernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% COPYRIGHT : Michalis K. Titsias, 2009
%
% KERN

params = [kern.variance kern.inputScales];
if nargout > 1
  names = {'variance'};
  for i = 1:length(kern.inputScales)
    names{1+i} = ['input scale ' num2str(i)];
  end
end


text:
function obj=geom_jitter(obj,varargin)
% geom_jitter Display data as jittered points
%
% Example syntax (default arguments): gramm_object.geom_jitter('width',0.2,'height',0.2)
% In case datapoints are grouped together and are hard to see,
% it's possible to randomly jitter them in an area of width
% 'width' and height 'height' using this function.

p=inputParser;
my_addParameter(p,'width',0.2);
my_addParameter(p,'height',0);
my_addParameter(p,'dodge',0);
my_addParameter(p,'alpha',1);
parse(p,varargin{:});

obj.geom=vertcat(obj.geom,{@(dobj,dd)my_jitter(dobj,dd,p.Results)});
obj.results.geom_jitter_handle={};
end

function  hndl=my_jitter(obj,draw_data,params)

draw_data.x=comb(draw_data.x);
draw_data.y=comb(draw_data.y);

draw_data.x=dodger(draw_data.x,draw_data,params.dodge);

draw_data.x=draw_data.x+rand(size(draw_data.x))*params.width-params.width/2;
draw_data.y=draw_data.y+rand(size(draw_data.y))*params.height-params.height/2;

%We adjust axes limits to accomodate for the jittering
if max(draw_data.x)>obj.plot_lim.maxx(obj.current_row,obj.current_column);
    obj.plot_lim.maxx(obj.current_row,obj.current_column)=max(draw_data.x);
end
if min(draw_data.x)<obj.plot_lim.minx(obj.current_row,obj.current_column);
    obj.plot_lim.minx(obj.current_row,obj.current_column)=min(draw_data.x);
end
if max(draw_data.y)>obj.plot_lim.maxy(obj.current_row,obj.current_column);
    obj.plot_lim.maxy(obj.current_row,obj.current_column)=max(draw_data.y);
end
if min(draw_data.y)<obj.plot_lim.miny(obj.current_row,obj.current_column);
    obj.plot_lim.miny(obj.current_row,obj.current_column)=min(draw_data.y);
end

%hndl=my_point(obj,draw_data);
hndl=line(draw_data.x,draw_data.y,'LineStyle','none','Marker',draw_data.marker,'MarkerEdgeColor','none','markerSize',draw_data.point_size,'MarkerFaceColor',draw_data.color);

set_alpha(hndl,1,params.alpha);

obj.results.geom_jitter_handle{obj.result_ind,1}=hndl;
end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
function Mij = get_block(M, i, j, row_range, col_range)
	%% ================== File info ==========================
	% Author		: Tiep Vu (http://www.personal.psu.edu/thv102/)
	% Time created	: 1/27/2016 2:46:04 AM
	% Last modified	: 1/27/2016 2:46:09 AM
	% Description	: Get a block matrix of a big matrix X.
	% 	INPUT
	%		M        : the big matrix.
	%			M = [	M11, M12, ..., M1m]
	% 					M21, M22, ..., M2m]
	%					..................
	% 					Mn1, Mn2, ..., Mnm]
	%		i        : row block index
	%		j        : column block index
	%		row_range: a vector storing the last index of each block. row_range(1) = 0.
	%					i-th block is indexed by row_range(i)+1: row_range(i+1).
	%		col_range: a vector storing the last index of each block. row_range(1) = 0.
	%					i-th block is indexed by col_range(i)+1: col_range(i+1).
	% 	OUTPUT
	%		Mi: output block matrix
	%
	%% ================== end File info ==========================
	range1 = row_range(i) + 1: row_range(i + 1);
	range2 = col_range(j) + 1: col_range(j + 1);
	Mij = M(range1, range2);
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
%
%   FILE NAME:    addManuallyCuratedMetAssoc.m
%
%   PURPOSE: Incorporate manual curation results of met association
%            with Recon/BiGG/MetaNetX ids to ihumanMets2MNX_v2.mat
%

% Load met association to MNX/BiGG
load('ihumanMets2MNX_v2.mat');

% Get the list of empty elements of Recon3D/BiGG
emptyList=find(cellfun(@isempty, ihuman.metRecon3DID));
MNXID=reformatElements(ihuman.metMNXID,'cell2str');
for i=1:length(emptyList)
		m=emptyList(i);
		fprintf('%s\t%s\t%s\n',ihuman.mets{m},ihuman.metNames{m},MNXID{m});
end

% Locate the ones with multiple Recon3D associations
HMRmetsNoComp=regexprep(ihuman.mets,'\w$','');
multi_index=find(cellfun(@numel, ihuman.metRecon3DID)>1);
% Obtain the non-unique assoc
metString=reformatElements(ihuman.metRecon3DID,'cell2str');
[multiMetAssoc.Recon3DID, I, ~]=unique(metString(multi_index));
multiMetAssoc.HMRID=HMRmetsNoComp(multi_index(I));
multiMetAssoc.MNXID=reformatElements(ihuman.metMNXID(multi_index(I)),'cell2str');
for i=1:length(multiMetAssoc.HMRID)
		fprintf('%s\t%s\t%s\n',multiMetAssoc.HMRID{i},multiMetAssoc.MNXID{i},multiMetAssoc.Recon3DID{i});
end

% These ones with missing and multiple Recon3D/BiGG ids were organized into
% an excel file metaboliteCuration_20180605.xlsx for manual cuartion

% The following updates are based on the manual curation results
% 1. For the elements without Recon3D/BiGG assoc
ihuman.metRecon3DID{find(strcmp('m00077p',ihuman.mets))}={'CE2416'};
ihuman.metRecon3DID{find(strcmp('m00678m',ihuman.mets))}={'dec24dicoa'};
ihuman.metRecon3DID{find(strcmp('m00678p',ihuman.mets))}={'dec24dicoa'};
ihuman.metRecon3DID{find(strcmp('m03035m',ihuman.mets))}={'dece3coa'};
ihuman.metRecon3DID{find(strcmp('m03035p',ihuman.mets))}={'dece3coa'};
ihuman.metRecon3DID{find(strcmp('m00980m',ihuman.mets))}={'dece4coa'};
ihuman.metRecon3DID{find(strcmp('m00980p',ihuman.mets))}={'dece4coa'};
ihuman.metRecon3DID{find(strcmp('m01422c',ihuman.mets))}={'cbtnCCP'};
ihuman.metRecon3DID{find(strcmp('m01942g',ihuman.mets))}={'gd1a_hs'};

ihuman.metMNXID{find(strcmp('m01942g',ihuman.mets))}={'MNXM11644','MNXM8637'};
ihuman.metMNXID{find(strcmp('m02149c',ihuman.mets))}={'MNXM56888'};
ihuman.metMNXID{find(strcmp('m02147m',ihuman.mets))}={'MNXM527231'};
ihuman.metMNXID{find(strcmp('m02147x',ihuman.mets))}={'MNXM527231'};
ihuman.metMNXID{find(strcmp('m02147s',ihuman.mets))}={'MNXM527231'};
ihuman.metMNXID{find(strcmp('m02147c',ihuman.mets))}={'MNXM527231'};

% 2. For the elements with multiple Recon3D/BiGG association
ihuman.metRecon3DID{find(strcmp('m00490c',ihuman.mets))}={'mag_hs'};
ihuman.metRecon3DID{find(strcmp('m02733c',ihuman.mets))}={'pa_hs'};
ihuman.metRecon3DID{find(strcmp('m02733r',ihuman.mets))}={'pa_hs'};
ihuman.metRecon3DID{find(strcmp('m00240c',ihuman.mets))}={'dag_hs'};
ihuman.metRecon3DID{find(strcmp('m00240g',ihuman.mets))}={'dag_hs'};
ihuman.metRecon3DID{find(strcmp('m00240n',ihuman.mets))}={'dag_hs'};
ihuman.metRecon3DID{find(strcmp('m01818c',ihuman.mets))}={'11_cis_retfa'};
ihuman.metRecon3DID{find(strcmp('m01818s',ihuman.mets))}={'11_cis_retfa'};
ihuman.metRecon3DID{find(strcmp('m01818x',ihuman.mets))}={'11_cis_retfa'};
ihuman.metRecon3DID{find(strcmp('m02959s',ihuman.mets))}={'tag_hs'};

save('ihumanMets2MNX_v2.mat','ihuman');  % 2018-06-12
% Some MNX IDs also need refinement (to be continued)

% 3. Generate the array structure (metAssocHMR2Recon3D) of metabolite
% association for further curation and model integration, by trimming
% off duplicate mets in multiple compartments
HMRmets=regexprep(ihuman.mets,'\w$','');   % HMR met ids without comp id
[metAssocHMR2Recon3.metHMRID, I, ~]=unique(HMRmets,'stable');
metAssocHMR2Recon3.metBiGGID=ihuman.metBiGGID(I);
metAssocHMR2Recon3.metMNXID=ihuman.metMNXID(I);
metAssocHMR2Recon3.metRecon3DID=ihuman.metRecon3DID(I);
metAssocHMR2Recon3.metNames=ihuman.metNames(I);
save('metAssocHMR2Recon3.mat','metAssocHMR2Recon3');  % 2018-06-17

% 4. Detect mets associated from one Recon3D id to multiple HMR ids
% get the array of non-empty Recon3D met assocations
Recon3DID=reformatElements(metAssocHMR2Recon3.metRecon3DID,'cell2str');
Recon3DID_nonEmpty=Recon3DID(getNonEmptyList(Recon3DID));

% get the array of unique Recon3D met id and the occurrences
check=countFrequency(Recon3DID_nonEmpty);
list=find([check.frequency{:}] > 1);   % Recon3D mets with multiple occurrences
for i=1:numel(list)
		m=list(i);
		ind=find(strcmp(check.uniqueList{m},Recon3DID));
		fprintf('%s\t%s\n',check.uniqueList{m},strjoin(metAssocHMR2Recon3.metHMRID(ind),';'));
end
% These Recon3D mets associated to multiple HMR ids were subjected to
% manual cuartion 2018-06-18

% 5. The following updates of HMR-Recon3D assoc are based on the curation results
% in excel file metaboliteCuration_20180618_HW.xlsx
metAssocHMR2Recon3.metRecon3DID{find(strcmp('m00591',metAssocHMR2Recon3.metHMRID))}='';
metAssocHMR2Recon3.metRecon3DID{find(strcmp('m00352',metAssocHMR2Recon3.metHMRID))}='';
metAssocHMR2Recon3.metRecon3DID{find(strcmp('m00379',metAssocHMR2Recon3.metHMRID))}='';
metAssocHMR2Recon3.metRecon3DID{find(strcmp('m00555',metAssocHMR2Recon3.metHMRID))}='';
metAssocHMR2Recon3.metRecon3DID{find(strcmp('m01123',metAssocHMR2Recon3.metHMRID))}={'CE7097'};
metAssocHMR2Recon3.metRecon3DID{find(strcmp('m01911',metAssocHMR2Recon3.metHMRID))}='';
metAssocHMR2Recon3.metRecon3DID{find(strcmp('m01942',metAssocHMR2Recon3.metHMRID))}='';
metAssocHMR2Recon3.metRecon3DID{find(strcmp('m02410',metAssocHMR2Recon3.metHMRID))}='';
metAssocHMR2Recon3.metRecon3DID{find(strcmp('m02487',metAssocHMR2Recon3.metHMRID))}='';
metAssocHMR2Recon3.metRecon3DID{find(strcmp('m02578',metAssocHMR2Recon3.metHMRID))}='';

% 6. Associate Recon3D mets 'M00196' and 'protein' to HMR2 met 'm00196'
ind=find(strcmp(metAssocHMR2Recon3.metHMRID,'m00196'));
metAssocHMR2Recon3.metRecon3DID{ind}{2}='protein';
save('metAssocHMR2Recon3.mat','metAssocHMR2Recon3');  % 2018-08-03

% 7. Synchronize curtated metabolite association previously in % 5 and % 6
% from 'metAssocHMR2Recon3.mat' to 'ihumanMets2MNX_v2.mat'
load('ihumanMets2MNX_v2.mat');
ihuman.metRecon3DID{find(strcmp('m00591c',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m00591s',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m00352r',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m00379c',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m00555c',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m00555g',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m00555r',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m01123c',ihuman.mets))}{1}='CE7097';
ihuman.metRecon3DID{find(strcmp('m01911c',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m01942g',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02410c',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02410m',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02410r',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02487c',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02487m',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02578c',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02578m',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02578n',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02578p',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02578r',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02578s',ihuman.mets))}='';
ihuman.metRecon3DID{find(strcmp('m02578x',ihuman.mets))}='';
% Associate Recon3D mets 'M00196' and 'protein' to HMR2 met 'm00196'
ihuman.metRecon3DID{find(strcmp('m00196c',ihuman.mets))}{2}='protein';
save('ihumanMets2MNX_v2.mat','ihuman');  % 2018-09-03

author: SysBioChalmers

repo: Human-GEM


text:
function [Nimages, NimagesFullyAnnotated, NumberOfAnnotatedObjects, NumberOfAnnotatedObjectsWithParts, NumberOfParts] = SUNstats(D)
%
%
% Number of images
% Number of fully annotated images
% Number of annotated objects
% Number of annotated root objects with parts
% Number of parts

% Number of images
Nimages = length(D);

% Number of fully annotated images
disp('Counting fully labeled images')
relativearea = LMlabeledarea(D);
NimagesFullyAnnotated = sum(relativearea>.9);

% Number of annotated objects

% Number of annotated root objects with parts
[counts, countrootwithparts, countparts] = LMcountparts(D);

NumberOfAnnotatedObjects = sum(counts);
NumberOfAnnotatedObjectsWithParts = sum(countrootwithparts);
NumberOfParts = sum(countparts);

fprintf('Number of images= %d \n', Nimages)
fprintf('Number of images (90%%) labeled = %d \n', NimagesFullyAnnotated)
fprintf('Number of annotated objects = %d \n', NumberOfAnnotatedObjects)
fprintf('Number of annotated objects with parts = %d \n', NumberOfAnnotatedObjectsWithParts)
fprintf('Number of annotated parts = %d \n', NumberOfParts)

author: CSAILVision

repo: LabelMeToolbox


text:
% Exact solution of Riemann problem; plot adapted to Osher scheme

% Theory in Section 10.2 of:

% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% This program is called by Euler schemes

% Functions called: f

global  PRL  CRL MACHLEFT  gamma  pleft  pright  rholeft  rhoright  uleft...
	uright  tend  lambda		% lambda = dt/dx

gammab = 1/(gamma - 1); gam1 = gamma-1;

%  Assumed structure of exact solution
%
%    \         /      |con |       |s|
%     \   f   /       |tact|       |h|
% left \  a  /  state |disc| state |o| right
% state \ n /    2    |cont|   3   |c| state
%   1    \ /          |tinu|       |k|   4
%         |           |ity |       | |

PRL = pright/pleft;
cright = sqrt(gamma*pright/rhoright); cleft = sqrt(gamma*pleft/rholeft);
CRL = cright/cleft;
MACHLEFT = (uleft - uright)/cleft;

p34 = fzero('f',3);		% p34 = p3/p4

tend = n*dt			% Exact time

p3 = p34*pright; alpha = (gamma+1)/(gamma-1);
rho3 = rhoright*(1+alpha*p34)/(alpha+p34);
rho2 = rholeft*(p34*pright/pleft)^(1/gamma);
u2 = uleft-uright+(2/(gamma-1))*cleft*...
  (1-(p34*pright/pleft)^((gamma-1)/(2*gamma)));
c2 = sqrt(gamma*p3/rho2);

spos = 0.5 + ...		% Shock position
 tend*cright*sqrt((gamma-1)/(2*gamma) + (gamma+1)/(2*gamma)*p34) + tend*uright;
conpos = 0.5 + u2*tend + tend*uright;	% Position of contact discontinuity
pos1 = 0.5 + (uleft - cleft)*tend;	% Start of expansion fan
pos2 = 0.5 + (u2+uright-c2)*tend;	% End of expansion fan

xx = 0:0.002:1;			% Preallocation
pexact = zeros(size(xx)); uexact= zeros(size(xx)); rhoexact = zeros(size(xx));
machexact = zeros(size(xx));  cexact = zeros(size(xx));
for i = 1:length(xx)
  if xx(i) <= pos1
    pexact(i) = pleft;    rhoexact(i) = rholeft;
    uexact(i) = uleft;    cexact(i) = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  elseif xx(i) <= pos2
    pexact(i) = pleft*(1+(pos1-xx(i))/(cleft*alpha*tend))^(2*gamma/(gamma-1));
    rhoexact(i) = rholeft*(1+(pos1-xx(i))/(cleft*alpha*tend))^(2/(gamma-1));
    uexact(i) = uleft + (2/(gamma+1))*(xx(i)-pos1)/tend;
    cexact(i) = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  elseif xx(i) <= conpos
    pexact(i) = p3;    		rhoexact(i) = rho2;
    uexact(i) = u2+uright;    	cexact(i) = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  elseif xx(i) <= spos
    pexact(i) = p3;     	rhoexact(i) = rho3;
    uexact(i) = u2+uright;      cexact(i) = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  else
    pexact(i) = pright;    	rhoexact(i) = rhoright;
    uexact(i) = uright;    	cexact(i) = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  end
end
entroexact = log(pexact./rhoexact.^gamma);

subplot(2,3,1), 	plot(xx,rhoexact,'-');
subplot(2,3,2), 	plot(xx,uexact,'-');
subplot(2,3,3),		plot(xx,pexact,'-');
subplot(2,3,4),		plot(xx,machexact,'-');
subplot(2,3,5),		enthexact = (gamma/(gamma-1))*pexact./rhoexact;
			plot(xx,enthexact,'-');
subplot(2,3,6),		plot(xx,entroexact,'-');


text:
function [] = v2structDemo1()
%  V2STRUCTDEMO1 is a demo function showing some examples for using the function V2STRUCT
%
%% Last update:
%    30.06.2011, Adi N.

%%
fprintf(['\nRunning demo function with examples for using the function vs2truct:\n',...
        'Press any key between breaks.\n']);
pause()
fprintf(['\nInitialize demo variables:\n',...
         'x = zeros(3); x2 = ones(3); y = ''Testing123''; z = cell(2,3);\n',...
         'fieldNames1 = {''fieldNames'',''x'',''y'',''z''};\n',...
         'fieldNames2 = {''fieldNames'',''a'',''b'',''c''};\n',...
         'fieldNames3 = {''fieldNames'',''x''};\n',...
         'nameOfStruct2Update = ''S'';\n']);
pause()
x = zeros(3); x2 = ones(3); y = 'Testing123'; z = cell(2,3);
fieldNames1 = {'fieldNames','x','y','z'};
fieldNames2 = {'fieldNames','a','b','c'};
fieldNames3 = {'fieldNames','x'};
nameOfStruct2Update = 'S';

fprintf(['\nExamples for packing and unpacking:\n',...
         '\nPacking:\n']);
pause()
fprintf('\nExample 1: Variables x,y and z are packed to structure S\n');
fprintf('S = v2struct(x,y,z) ');
pause()
S = v2struct(x,y,z)

pause()
fprintf(['Example 2: Field names defined according to the cell array fieldNames.\n',...
         'variables with the names in fieldNames must exist.\n']);
clear S
fprintf('S = v2struct(fieldNames1)');
pause()
S = v2struct(fieldNames1)

pause()
fprintf('Example 3: Same as the first but arguments are passed explicitly\n');
clear S
fprintf('S = v2struct(zeros(3),''Testing123'',cell(2,3),fieldNames1)');
pause()
S = v2struct(zeros(3),'Testing123',cell(2,3),fieldNames1)

pause()
fprintf(['Example 4: Field names defined by content of fieldNames2 while the values are\n',...
         'set according to the passed arguments. In this case the structure S returned would\n'...
         'be: S.a=x, S.b=y, S.c=z\n']);
clear S
fprintf('S = v2struct(x,y,z,fieldNames2)');
pause()
S = v2struct(x,y,z,fieldNames2)

pause()
fprintf(['Example 5: Update structure S. The fields that would be updated are according to content\n',...
         'of fieldNames3. Note that you must pass a variable with the name\n',...
         '''nameOfStruct2Update'' placed before ''fieldNames3''. This variable should\n',...
         'contain the name of the structure you want to update as a string. Also note\n',...
         'that if you set an output structure name which is different than the one stated\n'...
         'in nameOfStruct2Update a new structure would be created and the structure that\n'...
         'was meant to be updated would not get updated.']);
fprintf(['S.oldField = ''field to be saved for future use''\n',...
        'S = v2struct(x2,nameOfStruct2Update,fieldNames3)\n']);
pause()
clear S
S.oldField = 'field to be saved for future use'
S = v2struct(x2,nameOfStruct2Update,fieldNames3)

pause()
fprintf('Example 6: Pack all variables in caller workspace.\n');
fprintf('S = v2struct');
pause()
S = v2struct

%%
pause()
fprintf(['\nExamples for Unpacking:\n']);
pause()

fprintf(['\nFisrt clear variables from ''packing'' and initialize variables for\n',...
         '''unpacking'':\n',...
         'clear all\n',...
         'S.x = zeros(3); S.y = ''Testing123''; S.z = cell(2,3);\n',...
         'fieldNames3 = {''fieldNames'',''x'',''z''};\n'])

clear all
S.x = zeros(3); S.y = 'Testing123'; S.z = cell(2,3); fieldNames3 = {'fieldNames','x','z'};

pause()
fprintf(['\nExample 1: Create or overwrite variables x, y, z in the caller environment\n',...
         'with the contents of the corresponding named fields.\n']);
pause()
fprintf(['v2struct(S)\n']);
pause()

v2struct(S)

fprintf(['\nExample 2: Assign the contents of the fields of the scalar structure S to\n',...
         'the variables a,b,c rather than overwriting variables in the caller. If there\n',...
         'are fewer output variables than there are fields in S, the remaining fields\n',...
         'are not extracted.\n']);
pause()
fprintf(['[a,b,c] = v2struct(S)\n']);

pause()

[a,b,c] = v2struct(S)

fprintf(['\nExample 3: Create or overwrite variables x and z in the caller with the\n',...
         'contents of the corresponding named fields.\n']);
pause()
fprintf(['v2struct(S,fieldNames3)\n']);
pause()
v2struct(S,fieldNames3)

fprintf(['\nExample 4: Assign the contents of the fields ''x'' and ''z'' defined by\n',...
         'fieldNames3 of the scalar structure S to the variables a and b rather than\n',...
         'overwriting variables in the caller. If there are fewer output variables than\n',...
         'there are fields in S, the remaining fields are not extracted.\n']);
pause()
fprintf(['[a,b] = v2struct(S,fieldNames3)\n']);
pause()
[a,b] = v2struct(S,fieldNames3)

fprintf(['\nExample 5: Unpack variable ''z'' only without overwriting variable ''x''. Note\n',...
         'the addition of the field named ''avoidOverWrite'' to the structure to be\n',...
         'unpacked. This is mandatory in order to make this functionality work. The\n',...
         'contents of this field can be anything, it does not matter.\n']);

pause()
fprintf(['S.avoidOverWrite = ''foo(contents does not matter)'';\n',...
         'x = ''do not overwrite me'';\n',...
         'clear y\n']);
pause()
S.avoidOverWrite = 'foo(contents does not matter)';
x = 'do not overwrite me';
v2struct(S)
S

fprintf('\nend of v2structDemo1\n');

repo: buzcode


text:
% FORMAT:
%
% sampval = time2sample(timeop, time, fs, rounding, offset)
%
% INPUTS:
%
% timeop      - second or millisecond; 0 or 1
% timeval      - time value
% fs          - sample rate
%
% Optional inputs:
%
% rounding    - 1 means round the result {default}, 0 means do not round
% offset      - offset time in samples (you can use time2sample recursively here). Default 0
%
% OUTPUT:
%
% sampval     - time in samples
%
%
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2013

function sampval = time2sample(timeop, timeval, fs, rounding, offset)
if nargin<1
        help time2sample
        return
end
if nargin<5
        offset = 0;
end
if nargin<4
        rounding = 1;
end
if nargin<3
        fs = 1;
end
if nargin<2
        error('Two inputs are requiered at least.')
end
if timeop==0 % sec
    kktime=1;
elseif timeop==1 % msec
    kktime=1000;
else
    error('unknow timeop input.')
end
sampval = offset + timeval*fs/kktime;
if rounding
        sampval = round(sampval);
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function X = caffe_init(X, solver_file)

if ~exist(solver_file,'file') || isempty(solver_file)
    error('You need a solver prototxt definition.')
end
X.Solver_ = caffe.Solver(solver_file); % to cpp
if isfield(X, 'state_file') && ~isempty(X.state_file) && exist(X.state_file, 'file')
    X.Solver_.restore(X.state_file);
    fprintf('resume from snapshot file %s...\n', X.state_file);
elseif isfield(X, 'model_file') && ~isempty(X.model_file)
    X.Solver_.net.copy_from(X.model_file);
    fprintf('resume from caffemodel file %s...\n', X.model_file);
end

% set caffe mode
if strcmp(X.solver_mode,'GPU')
    caffe.set_mode_gpu();
    caffe.set_device(X.device_id);
else
    caffe.set_mode_cpu();
end
fprintf('Done with init %s\n', solver_file);
end

author: rwenqi

repo: GFN-dehazing

save_path: github-repos/MATLAB/rwenqi-GFN-dehazing


text:
function rs=sonomeasure(segArray,rpos,polygon,ns)
rx=rpos(1);
ry=rpos(2);
n=size(segArray,1);
if nargin<4
    ns=24;
    if nargin<3
        polygon=false;
    end
end
if polygon
    segArray=[segArray(1:end-1,:),segArray(2:end)];
end
a=2*pi/ns;%angle between two adjacent sensors
% Simulate the ultrasonic data
% Measure: Point Discretization %to be improved
dl=.1;
segadd=[];
for k=1:n
    lx=segArray(k,1)-segArray(k,3);
    ly=segArray(k,2)-segArray(k,4);
    l=sqrt(lx^2+ly^2);
    if lx
        xadd=segArray(k,1):-(dl*lx/l):segArray(k,3);
        yadd=(xadd-segArray(k,1))*ly/lx+segArray(k,2);
    else
        yadd=segArray(k,2):-(dl*ly/l):segArray(k,4);
        xadd=segArray(k,1)*ones(1,length(yadd));
    end
    segadd=[segadd;[xadd',yadd']];
end
segadd=[segadd;segArray(:,1:2);segArray(:,3:4)];

%determine sections
nadd=size(segadd,1);
nsect=zeros(1,nadd);
ssect=zeros(ns,nadd+1);
for k=1:nadd
    ak=vectorangle([1,0],segadd(k,:)-[rx,ry],1)/a;
    nsect(k)=ceil(ak)+(ak==0);
    sk=nsect(k);
    ssect(sk,1)=ssect(sk,1)+1;
    ssect(sk,ssect(sk,1)+1)=k;
end
%choose min
rs=zeros(1,ns);%distance radius of probing
for k=1:ns
    if ssect(k,1)
        pts=segadd(ssect(k,2:1+ssect(k,1)),:);
        rs(k)=min((pts(:,1)-rx).^2+(pts(:,2)-ry).^2);rs(k)=sqrt(rs(k));
    else
        rs(k)=0;%inf;
    end
end
end


text:
%% t_meshFromClass
%
%   Convert itkGray Class file to mesh using Matlab tools.
%
% See also: t_meshShow
%
% BW (c) Stanford VISTA Team

%% You might first check whether you can visualize a VISTASOFT mesh
% Either run: t_meshShow
% Or, just run this code
%  load(fullfile(mrvDataRootPath,'anatomy','anatomyNIFTI','leftMesh.mat'));
%  meshVisualize(msh);

%% Build a matlab mesh from the itkGray class file

niCFile = fullfile(mrvDataRootPath,'anatomy','anatomyV','t1_class.nii.gz');
niClass = niftiRead(niCFile);
Ds = uint8(niClass.data);

% The ITKGRAY class labels are
%    0: unlabeled
%    1: CSF
%    2: Subcortical
%    3: left white matter
%    5: left gray matter
%    4: right white matter
%    6: right gray matter

% Set the labels for everything that is not left gray or white to
% unlabeled.
Ds(Ds == 1) = 0;
Ds(Ds == 2) = 0;
Ds(Ds == 4) = 0;
Ds(Ds == 5) = 0;
Ds(Ds == 6) = 0;
% showMontage(double(Ds))

%% Matlab calculations reducing the patches and computing normals

% Make a picture of the cortical mesh.  We have all the information.
fv = isosurface(Ds,1);

% Reduce to about 30K faces.  That's enough
nFaces = size(fv.faces,1);
maxFaces = 4*10^4;    % 40,000 faces?
if nFaces > maxFaces
    p = maxFaces/nFaces;
    fv = reducepatch(fv.faces,fv.vertices, p);
end

% Smooth
smoothMode = 1; nIter = 3;
fv = smoothpatch(fv,smoothMode,nIter);

% mrvNewGraphWin;
% patch(fv, 'FaceColor','red','EdgeColor','none');
% view(3); daspect([1,1,1]); axis tight
% camlight; camlight(-80,-10); lighting phong;
% title('Iso surface via isonormals')

%% Now, take the Matlab information and put it into VISTASOFT mesh
mmPerVox = [1 1 1];
windowID = 1000;
actor = 33;
msh = meshFV2msh(fv,mmPerVox,windowID, actor);
meshVisualize(msh);

%% End

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function M = DiagM(n)
% function M = DiagM(n)
% M = TraceM(n, 'lo')
% Return "Diagonal-matrix", when applied on vectorized matrix (of order n)
% It gives the diagonal:
%
% Ouput are sparse
%
% DiagM(size(A))*A(:) == diag(A); % <- all true
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% Date: 21/March/2009

if isscalar(n)
    n = [n n];
end
I = idiag(n,0);

% Result
M = sparse((1:length(I)).', I, 1, length(I), prod(n));

end


text:
function Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

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

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%

for i = 1:size(X, 1)
    x = X(i, :)';
    Z(i, :) = x' * U(:, 1:K);
end

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

end

author: AvaisP


text:
Z = W*X;cz = [Z(:,ind)];

if exist('ZV2')
    [ZV,ZV2{2}] = ZV_gen2(Z,cz,ind,ZV2{1});
else
    [ZV,ZV2] = ZV_gen(Z,cz,Y,ind);
end;

CalcG;

RegMat = GG(ind,:) + Ib*RegEpsilon;

alpha = (Hhat + lambda*RegMat)\hhat;
Is_tmp = IsCalc(Hhat,hhat,alpha);

author: BehnoodRasti


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

% problem 3 -  Polar form and real and imaginary parts of the DFT of cos(n*pi/3)

n=0:8;
x=cos((1/3)*pi*n);
Xk=fft(x);
Xk.'

A=abs(Xk).*exp(j*angle(Xk));
A.'

B=real(Xk)+j*imag(Xk);
B.'


text:
function [newProj, rectGrid] = gridding(cg, proj )

indexT = 1:cg.nt;
indexS = 1:cg.ns;

ws = (cg.ns+1)/2;
wt = (cg.nt+1)/2;

coordT = (indexT - wt) * cg.dt;
coordS = (indexS - ws) * cg.ds;
d = cg.dsd;
r = cg.dso;

matrix = zeros(cg.nt, cg.ns);
%create matrix
for i = 1:1:cg.ns
	matrix(:,i) = coordT .* sqrt(r^2-coordS(i)^2) / d;
end

%find the maximum value
upperLimit = max(max(matrix));
lowerLimit = min(min(matrix)); % should just be -max

%create new grid in the vertical direction
new_dt = (upperLimit-lowerLimit)/cg.nt; % original
%new_dt = cg.dt;
rectGrid = lowerLimit+(upperLimit-lowerLimit)/cg.nt:new_dt:upperLimit;
rectGrid = rectGrid'; % 1d set of "t" samples we want

%must also loop over the projection angles
newProj = zeros(cg.ns, numel(rectGrid), cg.na);
for j = 1:1:cg.na
	for i = 1:1:cg.ns
		newProj(i,:,j) = ...
		interp1(matrix(:,i), proj(i,:,j), rectGrid, 'linear', 'extrap');
	end
end

end % gridding()


text:
function vesselness_wrapper = PTKComputeVesselnessFromHessianeigenvalues(hessian_eigs_wrapper, voxel_size)
    % PTKComputeVesselnessFromHessianeigenvalues. Vesselness filter for detecting blood vessels
    %
    %     PTKComputeVesselnessFromHessianeigenvalues computes a mutiscale
    %     vesselness filter based on Frangi et al., 1998. "Multiscale Vessel
    %     Enhancement Filtering". The filter returns a value at each point which
    %     in some sense representes the probability of that point belonging to a
    %     blood vessel.
    %
    %     This function takes in a PTKWraper object which can either contain a nx6
    %     matrix containing the 3 Hessian matrix eigenvalues for each of n
    %     points, or it can be an ixjxkx3 matrix representing the 3 Hessian
    %     matrix eigenvalues for an image of dimension ixjxk.
    %
    %     The output is a single vesselness value for each input point.
    %
    %     See the PTKVesselness plugin for example usage.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    % lam1 = smallest eigenvalue, lam3 = largest eigenvalue

    vesselness_wrapper = CoreWrapper;

    % The input matrix could be a linear set of points, or an image matrix
    if ndims(hessian_eigs_wrapper.RawImage) == 2
        lam1 = hessian_eigs_wrapper.RawImage(:, 1); % smallest
        lam2 = hessian_eigs_wrapper.RawImage(:, 2);
        lam3 = hessian_eigs_wrapper.RawImage(:, 3); % biggest
    else
        lam1 = hessian_eigs_wrapper.RawImage(:,:,:,1); % smallest
        lam2 = hessian_eigs_wrapper.RawImage(:,:,:,2);
        lam3 = hessian_eigs_wrapper.RawImage(:,:,:,3); % biggest
    end

    term_1 = abs(lam2./lam3); % Ra
    alpha = 0.5;
    term_1 = 1 - exp((-term_1.^2)./(2*alpha^2));

    term_2 = abs(lam1)./sqrt(abs(lam2.*lam3)); % Rb
    beta = 0.5;
    term_2 = exp((-term_2.^2)./(2*beta^2));

    term_3 = sqrt(lam1.^2 + lam2.^2 + lam3.^2); % S

    % Frangi et al. choose a noise threshold based on the Hessian mean.
    % However, we find a fixed experimentally-chosen threshold works better, but
    % this must be scaled by the voxel size
    multiple_matrix = voxel_size'*voxel_size;
    c_scaling = 1/mean(multiple_matrix(:));
    c = 400*c_scaling;

    term_3 = 1 - exp((-term_3.^2)./(2*c^2));

    check_signs = (lam2 <= 0) & (lam3 <= 0);

    vesselness_wrapper.RawImage = term_1.*term_2.*term_3.*check_signs;
end

author: tomdoel

repo: pulmonarytoolkit


text:
function y = gaussian_response(rect_size, sigma)
%GAUSSIAN_RESPONSE create the (fixed) target response of the correlation filter response
    [i1, i2] = circ_grid(rect_size(1), rect_size(2));
    y = exp(-(i1.^2 + i2.^2) / (2 * sigma^2));
end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function pth = mrtInstallSampleData(sourceFolder, projectName, ...
    dFolder, forceOverwrite, varargin)
%MRTINSTALLSAMPLEDATA Install sample data set on local path.
%
% Syntax:
%   pth = MRTINSTALLSAMPLEDATA(sourceFolder, projectName, ...
%      [dFolder], [forceOverwrite], varargin)
%
% Description:
%   A data set is stored on the OSF website (https://osf.io/t8pb2/). If
%   found locally in fullfile(vistaRootPath, 'local'), the local version is
%   used. If not, it is downloaded from the OSF page.
%
%
% Inputs
%     sourceFolder:     Name of folder where project is stored
%                         Examples: 'functional' | 'anatomy' | 'diffusion'
%     projectName:      Name of project
%                         Examples: 'mrBOLD_01' | 'prfInplane'
%     dFolder:          Destination folder
%                         [default = fullfile(vistaRootPath, 'local')]
%     forceOverwrite:   If true, unzip the project even if project folder
%                           is found. If false, do not unzip if project
%                           folder is found.
%                          [default = true]
%     varargin:         Pairs of parameters, values
%                           'filetype', {'zip' 'dat' 'mat' etc}
% Outputs
%      pth:             Full path to installed project folder or file
%
% Example:
%    pth = mrtInstallSampleData('functional', 'mrBOLD_01');
%
% See also:
%   MRVTEST

% Examples:
%{
  srcFolder = 'anatomy';
  project = 'anatomyNIFTI';
  dFolder = fullfile(vistaRootPath,'local');
  mrtInstallSampleData('anatomy','anatomyNIFTI',dFolder);

%}

% Check inputs
if notDefined('forceOverwrite'), forceOverwrite = true; end
if notDefined('dFolder'), dFolder = fullfile(vistaRootPath, 'local'); end
if exist('varargin', 'var')
    for ii = 1:2:length(varargin)
        switch varargin{ii}
            case 'filetype', filetype = varargin{ii+1};
            otherwise, error('%s parameter unrecognized.', varargin{ii});
        end
    end
end

% Check whether test data is installed locally
pth = fullfile(vistaRootPath, 'local', 'testData');
if isfolder(pth) % go on, local copy found
else
    fprintf('Test data not found. Downloading from OSF. This may take several minutes.');
    testdata = websave(fullfile(vistaRootPath, 'local', 'testData.zip'), ...
        'https://osf.io/t8pb2/download');
    unzip(testdata, fullfile(vistaRootPath, 'local'));
end

% By default, assume we are loading a zip file
if notDefined('filetype'), filetype = 'zip'; end

% Find the zip file
pth = fullfile(vistaRootPath, 'local', 'testData', sourceFolder);
d = dir (fullfile(pth, sprintf('%s*.%s', projectName, filetype)));

% Return the directory containing the unzipped data
pth = fullfile(dFolder, projectName);

% If the filetype was not a zip file, we are done. Otherwise unzip.
if ~strcmpi(filetype, 'zip')
    pth = fullfile(d.folder, d.name);
    return;
end

% Unzip
if exist(pth, 'dir') && ~forceOverwrite
    % If the project directory already exists and user has requested NOT to
    % force overwrite, then we do not unzip
    fprintf('Data directory %s already exists. Skipping unzip.\n', pth);
else
    zipfile = fullfile(d.folder, d.name);
    unzip(zipfile, dFolder);
end

if ~exist(pth, 'dir')
    warning('Cannot locate unzipped data directory %s', pth)
end

return


text:
function [fhat,xhat,fcount,retcode] = sims_csminit(fcn,x0,f0,g0,badg,H0,varargin)
% [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,...
%                                       P1,P2,P3,P4,P5,P6,P7,P8)
% retcodes: 0, normal step.  5, largest step still improves too fast.
% 4,2 back and forth adjustment of stepsize didn't finish.  3, smallest
% stepsize still improves too slow.  6, no improvement found.  1, zero
% gradient.
%---------------------
% Modified 7/22/96 to omit variable-length P list, for efficiency and compilation.
% Places where the number of P's need to be altered or the code could be returned to
% its old form are marked with ARGLIST comments.
%
% Fixed 7/17/93 to use inverse-hessian instead of hessian itself in bfgs
% update.
%
% Fixed 7/19/93 to flip eigenvalues of H to get better performance when
% it's not psd.
%
%tailstr = ')';
%for i=nargin-6:-1:1
%   tailstr=[ ',P' num2str(i)  tailstr];
%end
%ANGLE = .03;
ANGLE = .005;
%THETA = .03;
THETA = .3; %(0<THETA<.5) THETA near .5 makes long line searches, possibly fewer iterations.
FCHANGE = 1000;
MINLAMB = 1e-9;
% fixed 7/15/94
% MINDX = .0001;
% MINDX = 1e-6;
MINDFAC = .01;
fcount=0;
lambda=1;
xhat=x0;
f=f0;
fhat=f0;
g = g0;
gnorm = norm(g);
%
if (gnorm < 1.e-12) & ~badg % put ~badg 8/4/94
   retcode =1;
   dxnorm=0;
   % gradient convergence
else
   % with badg true, we don't try to match rate of improvement to directional
   % derivative.  We're satisfied just to get some improvement in f.
   %
   %if(badg)
   %   dx = -g*FCHANGE/(gnorm*gnorm);
   %  dxnorm = norm(dx);
   %  if dxnorm > 1e12
   %     disp('Bad, small gradient problem.')
   %     dx = dx*FCHANGE/dxnorm;
   %   end
   %else
   % Gauss-Newton step;
   %---------- Start of 7/19/93 mod ---------------
   %[v d] = eig(H0);
   %toc
   %d=max(1e-10,abs(diag(d)));
   %d=abs(diag(d));
   %dx = -(v.*(ones(size(v,1),1)*d'))*(v'*g);
%      toc
   dx = -H0*g;
%      toc
   dxnorm = norm(dx);
   if dxnorm > 1e12
      dx = dx*FCHANGE/dxnorm;
   end
   dfhat = dx'*g0;
   %end
   %
   %
   if ~badg
      % test for alignment of dx with gradient and fix if necessary
      a = -dfhat/(gnorm*dxnorm);
      if a<ANGLE
         dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g;
         % suggested alternate code:  ---------------------
         dx = dx*dxnorm/norm(dx)    % This keeps scale invariant to the angle correction
         % ------------------------------------------------
         dfhat = dx'*g;
         % dxnorm = norm(dx);  % this line unnecessary with modification that keeps scale invariant
      end
   end
   %
   % Have OK dx, now adjust length of step (lambda) until min and
   % max improvement rate criteria are met.
   done=0;
   factor=3;
   shrink=1;
   lambdaMin=0;
   lambdaMax=inf;
   lambdaPeak=0;
   fPeak=f0;
   lambdahat=0;
   while ~done
      if size(x0,2)>1
         dxtest=x0+dx'*lambda;
      else
         dxtest=x0+dx*lambda;
      end
      % home
      f = feval(fcn,dxtest,varargin{:});
      if f<fhat
         fhat=f;
         xhat=dxtest;
         lambdahat = lambda;
      end
      fcount=fcount+1;
      shrinkSignal = (~badg & (f0-f < max([-THETA*dfhat*lambda 0]))) | (badg & (f0-f) < 0) ;
      growSignal = ~badg & ( (lambda > 0)  &  (f0-f > -(1-THETA)*dfhat*lambda) );
      if  shrinkSignal  &   ( (lambda>lambdaPeak) | (lambda<0) )
         if (lambda>0) & ((~shrink) | (lambda/factor <= lambdaPeak))
            shrink=1;
            factor=factor^.6;
            while lambda/factor <= lambdaPeak
               factor=factor^.6;
            end
            %if (abs(lambda)*(factor-1)*dxnorm < MINDX) | (abs(lambda)*(factor-1) < MINLAMB)
            if abs(factor-1)<MINDFAC
               if abs(lambda)<4
                  retcode=2;
               else
                  retcode=7;
               end
               done=1;
            end
         end
         if (lambda<lambdaMax) & (lambda>lambdaPeak)
            lambdaMax=lambda;
         end
         lambda=lambda/factor;
         if abs(lambda) < MINLAMB
            if (lambda > 0) & (f0 <= fhat)
               % try going against gradient, which may be inaccurate
               lambda = -lambda*factor^6;
            else
               if lambda < 0
                  retcode = 6;
               else
                  retcode = 3;
               end
               done = 1;
            end
         end
      elseif  (growSignal & lambda>0) |  (shrinkSignal & ((lambda <= lambdaPeak) & (lambda>0)))
         if shrink
            shrink=0;
            factor = factor^.6;
            %if ( abs(lambda)*(factor-1)*dxnorm< MINDX ) | ( abs(lambda)*(factor-1)< MINLAMB)
            if abs(factor-1)<MINDFAC
               if abs(lambda)<4
                  retcode=4;
               else
                  retcode=7;
               end
               done=1;
            end
         end
         if ( f<fPeak ) & (lambda>0)
            fPeak=f;
            lambdaPeak=lambda;
            if lambdaMax<=lambdaPeak
               lambdaMax=lambdaPeak*factor*factor;
            end
         end
         lambda=lambda*factor;
         if abs(lambda) > 1e20;
            retcode = 5;
            done =1;
         end
      else
         done=1;
         if factor < 1.2
            retcode=7;
         else
            retcode=0;
         end
      end
   end
end

author: european-central-bank

repo: BEAR-toolbox