text:
function zmap(varargin)
    %% run ZMAP Version 7.X
    %
    % This files start up ZMAP.
    %
    % Options:
    %   -debug    : enables debugging functionality
    %   -restart  : clear all windows and variables, then restart zmap
    %   -initonly : set up zmap paths and prefs, but don't open a window
    %   -RT       : implement real-time mode
    %   -catalog alt_cat: use alt_cat as the prime catalog instead. It may be a variable name, a
    %               ZmapCatalog, or a value that canbe converted by ZmapCatalog
    %
    %
    % Options to be probably implemented
    %   -grid :
    %   -selector :
    %   -shape :
    %   -state :
    %   -bordersonly :
    %
    %
    %
    % MAJOR CHANGES to ZMAP from v 6.0:
    %  - event catalog is now a ZmapCatalog instead of a variable-sized array
    %  - dates/times are handled by MATLAB's built-in datetime and duration data types
    %  - duplicated code has been removed or consolidated
    %  - scripts have been turned into functions and classes (allowing zmap functions to be run programmatically)
    %  - GUI elements have been consolidated, updated
    %  - unreachable code (& scripts not appearing in any menus) have been removed
    %  - added ability to import catalogs from FDSN web services.
    %  - global variables have been minimized. Data sharing occurs through one class.
    %  - only result variables should appear in the base workspace.
    %  - requires minimum of matlab v 2018a,
    %
    %
    %  see README.md and the help for individual parts of ZMAP for other change details
    %

    % The matlab searchpaths are updated, existing windows closed.
    %
    %  Originally created by: Stefan Wiemer  12/94
    %  Modified by: Celso Reyes Spring 2017 - Winter 2018
    %
    % see also: ZmapCatalog, ZmapMainWindow

    global ZG
    % ZG (ZmapGlobal) provides access to all ZMAP's global variables
    % When variables are accessed directly via ZmapGlobal.Data.variablename, they
    % should not modify the original.
    % However, assigning ZmapGlobal.Data to a variable first, provides direct [read/write] access to
    % the variables.
    % This allows assignment.
    %
    %       ZmapGlobal.Data.ra = 100; % value is effectively ignored!
    %       ZG=ZmapGlobal.Data; %provide read/write access to global data.
    %       ZG.ra = 23  % changes ra globally

    % advise matlab where it can find everything zmap
    if ~exist('set_zmap_paths','file')
        p = mfilename('fullpath');
        p(end-length(mfilename) : end) = [];
        addpath(fullfile(p,'src'));
    end
    set_zmap_paths;

    startWindow      = true;
    debugMode        = false;

    possible_options = varargin(cellfun(@ischarlike,varargin));
    options          = lower(possible_options(startsWith(possible_options,'-')));

   setappdata(groot,'ZmapCoordinateSystem', CoordinateSystems.geodetic);

   % set application data

    for option = string(options)
        switch option
            case "-restart"
                msg.infodisp('Restarting ZMAP','initialization')
                restartZmap('restart');
                return

            case "-quit"
                msg.infodisp('Quitting ZMAP','initialization')
                restartZmap('quit');
                return

            case "-debug"
                msg.infodisp('debug mode enabled','initialization')
                debugMode = true;

            case "-initonly"
                msg.infodisp('Initializing Zmap without starting the default main window','initialization')
                startWindow = false;

            case "-rt"
                msg.infodisp('RealTimeMode enabled','initialization');
                setappdata(groot,'ZmapRealTimeMode', true);

            case "-cartesian"
                % must be set prior to acessing ZmapGlobal for the first time
                setappdata(groot,'ZmapCoordinateSystem', CoordinateSystems.cartesian);

                %{

            case {"-catalog"}
                vaa = varargin;
                notstrings =cellfun(@(x)~ischarlike(x),varargin);
                vaa(notstrings)={''};
                idx = find(vaa == "-catalog",1,'last'); %requres all members of vaa be strings
                c = varargin{idx+1};
                if ischarlike(c) && isvarname(c)
                    v=evalin('base',"whos('"+ c +"')");
                    if isempty(v)
                        msg.errordisp("cannot find the variable " + c, 'catalog specified');
                        return
                    end
                    alt_cat = evalin('base',"ZmapCatalog('Name', " + v.name + ")");
                    if isa(alt_cat, 'ZmapCatalog')
                        constructor_options.primeCatalog = alt_cat;
                    end
                    msg.infodisp("set primary catalog to the value from : " + v.name, 'catalog specified');
                else
                    constructor_options.primeCatalog = ZmapCatalog.from(c);
                    msg.infodisp("set primary catalog from a : " + class(c), 'catalog specified');
                end
                    %}

            otherwise
                msg.errordisp('Unknown ZMAP option: ' + option, 'Unknown ZMAP option, stopping');
                return
        end

    end

    switch getappdata(groot,'ZmapCoordinateSystem')
        case CoordinateSystems.geodetic
            setappdata(groot,'ZmapDefaultReferenceEllipsoid',referenceEllipsoid('earth','kilometer'));
        case CoordinateSystems.cartesian
            setappdata(groot,'ZmapDefaultReferenceEllipsoid',nonEllipsoid);
    end

    ZG = ZmapGlobal.Data;
    ZG.debug = debugMode;

    disp(['This is zmap.m - version ', ZmapGlobal.Data.zmap_version])
    % Set up the different computer systems
    sys = computer;

    if verLessThan('matlab',ZG.min_matlab_version)
        baseMsg = 'You are running a version of MATLAB older than %s. ZMAP %s requires MATLAB %s or newer';
        messtext = sprintf(baseMsg, ZG.min_matlab_release, ZG.zmap_version, ZG.min_matlab_version);
        msg.errordisp(messtext,  'Incompatible MATLAB Version');
        errordlg(messtext,'Warning!')
        pause(5)
        return
    end

    tested_systems = {'MAC','PCW'};
    if ~ismember( sys(1:3), tested_systems)
        msg.warndisp('ZMAP has not been tested on this computer type', 'Untested System')
        warndlg('Warning: ZMAP has not been tested on this computer type.','Untested System')
        pause(5)
    end

    assignin('base','ZG',ZmapGlobal.Data);

    % set local preference variables
    ini_zmap

    % set system dependent initial variables
    ini_zmap_sys

    % start the main zmap program
    if startWindow

        %get rid of message box that would exist if zmap was already opened without a catalog
        delete(findall(groot,'Tag','Msgbox_No Active Catalogs'));

        zw = findall(allchild(groot),'Tag','Zmap Main Window');
        s=sprintf('%d ZMAP windows exist\n', numel(zw));
        if ~isempty(zw)
            emptyzw = arrayfun(@(x)~isstruct(x.UserData) || isempty(x.UserData.catalog), zw);
            delete(zw(emptyzw));
            s = s + "... of which " + sum(emptyzw) + " were empty";
        end
        msg.dbdisp(s)
        cw = get(groot,'CurrentFigure');
        if isempty(cw)
            msg.dbdisp('No Figure currently exists');
        else
            switch cw.Tag
                case 'Zmap Main Window'
                    msg.dbdisp('ZMAP Window Exists, and is active')
                otherwise
                    % do nothing
            end
        end
        cw = figure;
        if ~isempty(ZG.primeCatalog) && questdlg('Open previous catalog?','ZMAP','Yes')=="Yes"
            ZmapMainWindow(cw, ZG.primeCatalog);
        else
            ZmapMainWindow(cw, ZmapCatalog);
        end
        if ~isappdata(groot,'ZmapShowTips') || getappdata(groot,'ZmapShowTips')
            show_a_tip();
            setappdata(groot,'ZmapShowTips',false);
        end
    end
end

author: CelsoReyes


text:
function test_two ( )

%*****************************************************************************80
%
%% TEST_TWO performs quadrature with SPINTERP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_TWO:\n' );
  fprintf ( 1, '  An example of the use of SPINTERP\n' );
  fprintf ( 1, '  for quadrature in multiple dimensions.\n' );

  addpath ( '../spinterp' );

  m = 2;
  fprintf ( 1, '  Using spatial dimension M = %d\n', m );
  level_max = 8;
  fprintf ( 1, '  Using maximum level LEVEL_MAX = %d\n', level_max );
%
%  Turn off warnings from SPINTERP that a given integral estimate is poor.
%
  warning ( 'off', 'MATLAB:spinterp:insufficientDepth' )
%
%  Compute the integral estimate for an increasing sequence of levels.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, ' Level    Points        Estimate      Error\n' );
  fprintf ( 1, '\n' );

  for level = 0 : level_max
%
%  Indicate that we wish to use exactly the sparse grid of level LEVEL.
%
    options = spset ( ...
      'MinDepth', level, ...
      'MaxDepth', level, ...
      'GridType', 'Chebyshev', ...
      'FunctionArgType', 'vector' );
%
%  Compute the sparse grid structure.
%
    z = spvals ( @f_two, m, [], options );
%
%  Estimate the integral.
%
    n = z.nPoints;

    q = spquad ( z );

    exact = 1.0;

    e = abs ( q - exact );

    fprintf ( 1, '  %4d  %8d  %14.6g  %14.6g\n', level, n, q, e );

  end

  rmpath ( '../spinterp' )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_TWO:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );
  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [y, P] = mlpClassPred(model, X)
% Multilayer perceptron classification prediction
% logistic activation function is used.
% Input:
%   model: model structure
%   X: d x n data matrix
% Ouput:
%   y: 1 x n label vector
%   P: k x n probability matrix
% Written by Mo Chen (sth4nth@gmail.com).
W = model.W;
b = model.b;
T = length(W);
Z = X;
for t = 1:T-1
    Z = sigmoid(W{t}'*Z+b{t});
end
P = softmax(W{T}'*Z+b{T});
[~,y] = max(P,[],1);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function  g =  lfmGradientSigmaH4AP(gamma1_p, gamma1_m, sigma2, t1, ...
    preFactor, preExp, mode)

% LFMGRADIENTSIGMAH4AP Gradient of the function h_i(z) with respect \sigma.
% FORMAT
% DESC Computes the gradient of the function h_i(z) with respect to the
% length-scale of the input "force", \sigma.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode: indicates in which way the vectors t1 and t2 must be transposed
% RETURN g : Gradient of the function with respect to \sigma.
%
% COPYRIGHT : Mauricio Alvarez, 2010

% KERN

if mode==0
    g =  lfmapGradientSigmaUpsilonVector(gamma1_p,sigma2, t1)*( preExp(:,1)/preFactor(1) - preExp(:,2)/preFactor(2)).' ...
        + lfmapGradientSigmaUpsilonVector(gamma1_m,sigma2, t1)*( preExp(:,2)/preFactor(3) - preExp(:,1)/preFactor(4)).';
else
    g =  lfmGradientSigmaUpsilonVector(gamma1_p,sigma2, t1)*(preExp(:,1)/preFactor(1) - preExp(:,2)/preFactor(2)).' ...
        + lfmGradientSigmaUpsilonVector(gamma1_m,sigma2, t1)*(preExp(:,2)/preFactor(3) - preExp(:,1)/preFactor(4)).';
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_suite = test_mdwt
  disp("mdwt")
  test_mdwt_1D
  test_mdwt_2D
  test_mdwt_compute_L1
  test_mdwt_compute_L2
  test_mdwt_compute_L3

function test_mdwt_1D
  x = makesig('LinChirp', 8);
  h = daubcqf(4, 'min');
  L = 2;  % For 8 values in x we would normally be L=2
  [y, L] = mdwt(x, h, L);
  y_corr = [1.1097 0.8767 0.8204 -0.5201 -0.0339 0.1001 0.2201 -0.1401];
  L_corr = 2;
assertVectorsAlmostEqual(y, y_corr, 'relative', 0.001);
assertEqual(L, L_corr);

function test_mdwt_2D
  x = [1 2 3 4; 5 6 7 8 ; 9 10 11 12; 13 14 15 16];
  h = daubcqf(4);
  y = mdwt(x, h);
  y_corr = [34.0000 -3.4641 0.0000 -2.0000; -13.8564 0.0000 0.0000 -2.0000; -0.0000 0.0000 -0.0000 -0.0000; -8.0000 -8.0000 0.0000 -0.0000];
assertVectorsAlmostEqual(y, y_corr, 'relative', 0.001);

function test_mdwt_compute_L1
  x = [1 2];
  h = daubcqf(4, 'min');
  [y, L] = mdwt(x, h);
assertEqual(L, 1);

function test_mdwt_compute_L2
  x = [1 2 3 4];
  h = daubcqf(4, 'min');
  [y, L] = mdwt(x, h);
assertEqual(L, 2);

function test_mdwt_compute_L3
  x = [1 2 3 4 5 6 7 8];
  h = daubcqf(4, 'min');
  [y, L] = mdwt(x, h);
assertEqual(L, 3);

author: ricedsp


text:
function output=rotateAround(image, pointY, pointX, angle, varargin)
% ROTATEAROUND rotates an image.
%   ROTATED=ROTATEAROUND(IMAGE, POINTY, POINTX, ANGLE) rotates IMAGE around
%   the point [POINTY, POINTX] by ANGLE degrees. To rotate the image
%   clockwise, specify a negative value for ANGLE.
%
%   ROTATED=ROTATEAROUND(IMAGE, POINTY, POINTX, ANGLE, METHOD) rotates the
%   image with specified method:
%       'nearest'       Nearest-neighbor interpolation
%       'bilinear'      Bilinear interpolation
%       'bicubic'       Bicubic interpolation
%    The default is 'nearest'.
%
%   Example
%   -------
%       image=imread('eight.tif');
%       imshow(rotateAround(image, 1, 1, 10));
%
%   See also IMROTATE, PADARRAY.

%   Contributed by Jan Motl (jan@motl.us)
%   $Revision: 1.0 $  $Date: 2013/02/22 16:58:01 $

% Parameter checking.
numvarargs = length(varargin);
if numvarargs > 1
    error('myfuns:somefun2Alt:TooManyInputs', ...
        'requires at most 1 optional input');
end
optargs = {'nearest'};    % Set defaults for optional inputs
optargs(1:numvarargs) = varargin;
[method] = optargs{:};    % Place optional args in memorable variable names

% Initialization.
[imageHeight imageWidth] = size(image);
centerX = floor(imageWidth/2+1);
centerY = floor(imageHeight/2+1);

dy = centerY-pointY;
dx = centerX-pointX;

% How much would the watched point shift if rotated around the img center.
[theta, rho] = cart2pol(-dx,dy);
[newX, newY] = pol2cart(theta+deg2rad(angle), rho);
shiftX = round(pointX-(centerX+newX));
shiftY = round(pointY-(centerY-newY));

% Pad the image to preserve the whole image during the rotation.
padX = abs(shiftX);
padY = abs(shiftY);

padded = padarray(image, [padY padX]);

% Rotate the image around the center.
rot = imrotate(padded, angle, 'crop', method);

% Crop the image.
output = rot(padY+1-shiftY:end-padY-shiftY, padX+1-shiftX:end-padX-shiftX);

author: Sable

repo: mcbench-benchmarks


text:
function imdb = getCifarImdb(opts)
% -------------------------------------------------------------------------
% Preapre the imdb structure, returns image data with mean image subtracted
unpackPath = fullfile(opts.dataDir, 'cifar-10-batches-mat');
files = [arrayfun(@(n) sprintf('data_batch_%d.mat', n), 1:5, 'UniformOutput', false) ...
  {'test_batch.mat'}];
files = cellfun(@(fn) fullfile(unpackPath, fn), files, 'UniformOutput', false);
file_set = uint8([ones(1, 5), 3]);

if any(cellfun(@(fn) ~exist(fn, 'file'), files))
  url = 'http://www.cs.toronto.edu/~kriz/cifar-10-matlab.tar.gz' ;
  fprintf('downloading %s\n', url) ;
  untar(url, opts.dataDir) ;
end

data = cell(1, numel(files));
labels = cell(1, numel(files));
sets = cell(1, numel(files));
for fi = 1:numel(files)
  fd = load(files{fi}) ;
  data{fi} = permute(reshape(fd.data',32,32,3,[]),[2 1 3 4]) ;
  labels{fi} = fd.labels' + 1; % Index from 1
  sets{fi} = repmat(file_set(fi), size(labels{fi}));
end

set = cat(2, sets{:});
data = single(cat(4, data{:}));

% remove mean in any case
dataMean = mean(data(:,:,:,set == 1), 4);
data = bsxfun(@minus, data, dataMean);

% normalize by image mean and std as suggested in `An Analysis of
% Single-Layer Networks in Unsupervised Feature Learning` Adam
% Coates, Honglak Lee, Andrew Y. Ng

if isfield(opts,'contrastNormalization')&&opts.contrastNormalization
  z = reshape(data,[],60000) ;
  z = bsxfun(@minus, z, mean(z,1)) ;
  n = std(z,0,1) ;
  z = bsxfun(@times, z, mean(n) ./ max(n, 40)) ;
  data = reshape(z, 32, 32, 3, []) ;
end

if isfield(opts,'whitenData') &&opts.whitenData
  z = reshape(data,[],60000) ;
  W = z(:,set == 1)*z(:,set == 1)'/60000 ;
  [V,D] = eig(W) ;
  % the scale is selected to approximately preserve the norm of W
  d2 = diag(D) ;
  en = sqrt(mean(d2)) ;
  z = V*diag(en./max(sqrt(d2), 10))*V'*z ;
  data = reshape(z, 32, 32, 3, []) ;
end

clNames = load(fullfile(unpackPath, 'batches.meta.mat'));

imdb.images.data = data ;
imdb.images.labels = single(cat(2, labels{:})) ;
imdb.images.set = set;
imdb.meta.sets = {'train', 'val', 'test'} ;
imdb.meta.classes = clNames.label_names;

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function tests = test_cosamp
  tests = functiontests(localfunctions);
end

function [A, x, b, k] = problem_1()
    m = 100;
    n = 1000;
    k = 4;
    A = spx.dict.simple.gaussian_dict(m, n);
    gen = spx.data.synthetic.SparseSignalGenerator(n, k);
    % create a sparse vector
    x =  gen.biGaussian();
    b = A*x;
end

function [dict, reps, signals, k] = problem_2()
    m = 200;
    n = 1000;
    k = 10;
    s = 500;
    dict = spx.dict.simple.gaussian_dict(m, n);
    gen = spx.data.synthetic.SparseSignalGenerator(n, k, s);
    % create a sparse vector
    reps =  gen.biGaussian();
    signals = dict*reps;
end

function test_cosamp_1(testCase)
    [A, x, b, k] = problem_1();
    solver = spx.pursuit.single.CoSaMP(A, k);
    result = solver.solve(b);
    cmpare = spx.commons.SparseSignalsComparison(x, result.z, k);
    %cmpare.summarize();
    verifyTrue(testCase, cmpare.all_have_matching_supports(1.0));
end

function test_cosamp_2(testCase)
    [dict, reps, signals, k] = problem_2();
    solver = spx.pursuit.single.CoSaMP(dict, k);
    ns = size(signals, 2);
    dd = size(dict, 2);
    recovered = zeros(dd, ns);
    for s=1:ns
        signal = signals(:, s);
        result = solver.solve(signal);
        recovered(:, s) = result.z;
    end
    cmpare = spx.commons.SparseSignalsComparison(reps, recovered, k);
    % cmpare.summarize();
    verifyTrue(testCase, cmpare.all_have_matching_supports(1.0));
end

author: indigits

repo: sparse-plex


text:
function st = spm_cfg_st
% SPM Configuration file for Slice Timing Correction
%__________________________________________________________________________
% Copyright (C) 2005-2016 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_st.m 7591 2019-05-15 10:57:38Z john $

%--------------------------------------------------------------------------
% scans Session
%--------------------------------------------------------------------------
scans         = cfg_files;
scans.tag     = 'scans';
scans.name    = 'Session';
scans.help    = {'Select images to slice-time correct.'};
scans.filter  = 'image';
scans.ufilter = '.*';
scans.num     = [1 Inf];
scans.preview = @(f) spm_check_registration(char(f));

%--------------------------------------------------------------------------
% generic Data
%--------------------------------------------------------------------------
generic         = cfg_repeat;
generic.tag     = 'generic';
generic.name    = 'Data';
generic.help    = {'Subjects or sessions. The same parameters specified below will be applied to all sessions.'};
generic.values  = {scans };
generic.num     = [1 Inf];

%--------------------------------------------------------------------------
% nslices Number of Slices
%--------------------------------------------------------------------------
nslices         = cfg_entry;
nslices.tag     = 'nslices';
nslices.name    = 'Number of Slices';
nslices.help    = {'Enter the number of slices.'};
nslices.strtype = 'n';
nslices.num     = [1 1];

%--------------------------------------------------------------------------
% tr TR
%--------------------------------------------------------------------------
tr         = cfg_entry;
tr.tag     = 'tr';
tr.name    = 'TR';
tr.help    = {'Enter the TR (in seconds).'};
tr.strtype = 'r';
tr.num     = [1 1];

%--------------------------------------------------------------------------
% ta TA
%--------------------------------------------------------------------------
ta         = cfg_entry;
ta.tag     = 'ta';
ta.name    = 'TA';
ta.help    = {
    'Enter the TA (in seconds). It is usually calculated as TR-(TR/nslices).'
    ''
    'You can simply enter this equation with the variables replaced by appropriate numbers.'
    ''
    'If the next two items are entered in milliseconds, this entry will not be used and can be set to 0.'
    }';
ta.strtype = 'r';
ta.num     = [1 1];

%--------------------------------------------------------------------------
% so Slice order
%--------------------------------------------------------------------------
so         = cfg_entry;
so.tag     = 'so';
so.name    = 'Slice order';
so.help    = {
    'Enter the slice order.'
    ''
    'Bottom slice = 1. Sequence types and examples of code to enter are given below:'
    ''
    '  ascending (first slice=bottom): [1:1:nslices]'
    ''
    '  descending (first slice=top): [nslices:-1:1]'
    ''
    '  interleaved (middle-top):'
    '     for k = 1:nslices'
    '         round((nslices-k)/2 + (rem((nslices-k),2) * (nslices - 1)/2)) + 1,'
    '     end'
    ''
    '  interleaved (bottom -> up): [1:2:nslices 2:2:nslices]'
    ''
    '  interleaved (top -> down): [nslices:-2:1, nslices-1:-2:1]'
    ''
    'Alternatively you can enter the slice timing in ms for each slice individually.'
    'If doing so, the next item (Reference Slice) will contain a reference time (in ms) instead of the slice index of the reference slice.'
    'For Siemens scanners, this can be acquired in MATLAB from the dicom header as follows (use any volume after the first one):'
    '   hdr = spm_dicom_headers(''dicom.ima'');'
    '   slice_times = hdr{1}.Private_0019_1029'
    'Note that slice ordering is assumed to be from foot to head. If it is not, enter instead: TR - INTRASCAN_TIME - SLICE_TIMING_VECTOR'
    }';
so.strtype = 'r';
so.num     = [1 Inf];

%--------------------------------------------------------------------------
% refslice Reference Slice
%--------------------------------------------------------------------------
refslice         = cfg_entry;
refslice.tag     = 'refslice';
refslice.name    = 'Reference Slice';
refslice.help    = {'Enter the reference slice.'
                    ''
                    'If slice times are provided instead of slice indices in the previous item, this value should represent a reference time (in ms) instead of the slice index of the reference slice.'};
refslice.strtype = 'r';
refslice.num     = [1 1];

%--------------------------------------------------------------------------
% prefix Filename Prefix
%--------------------------------------------------------------------------
prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the slice-time corrected image file(s).','Default prefix is ''a''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.def     = @(val)spm_get_defaults('slicetiming.prefix', val{:});

%--------------------------------------------------------------------------
% st Slice Timing
%--------------------------------------------------------------------------
st         = cfg_exbranch;
st.tag     = 'st';
st.name    = 'Slice Timing';
st.val     = {generic nslices tr ta so refslice prefix};
st.help    = {
    'Correct differences in image acquisition time between slices.'
    ''
    'Slice-time corrected files are prepended with an ''a''.'
    ''
    'Note: The sliceorder arg that specifies slice acquisition order is a vector of N numbers, where N is the number of slices per volume. Each number refers to the position of a slice within the image file. The order of numbers within the vector is the temporal order in which those slices were acquired. To check the order of slices within an image file, use the SPM Display option and move the cross-hairs to a voxel co-ordinate of z=1.  This corresponds to a point in the first slice of the volume.'
    ''
    'The function corrects differences in slice acquisition times. This routine is intended to correct for the staggered order of slice acquisition that is used during echo-planar scanning. The correction is necessary to make the data on each slice correspond to the same point in time. Without correction, the data on one slice will represent a point in time as far removed as 1/2 the TR from an adjacent slice (in the case of an interleaved sequence).'
    ''
    'This routine "shifts" a signal in time to provide an output vector that represents the same (continuous) signal sampled starting either later or earlier. This is accomplished by a simple shift of the phase of the sines that make up the signal. Recall that a Fourier transform allows for a representation of any signal as the linear combination of sinusoids of different frequencies and phases. Effectively, we will add a constant to the phase of every frequency, shifting the data in time.'
    ''
    'Shifter - This is the filter by which the signal will be convolved to introduce the phase shift. It is constructed explicitly in the Fourier domain. In the time domain, it may be described as an impulse (delta function) that has been shifted in time the amount described by TimeShift. The correction works by lagging (shifting forward) the time-series data on each slice using sinc-interpolation. This results in each time series having the values that would have been obtained had the slice been acquired at the same time as the reference slice. To make this clear, consider a neural event (and ensuing hemodynamic response) that occurs simultaneously on two adjacent slices. Values from slice "A" are acquired starting at time zero, simultaneous to the neural event, while values from slice "B" are acquired one second later. Without correction, the "B" values will describe a hemodynamic response that will appear to have began one second EARLIER on the "B" slice than on slice "A". To correct for this, the "B" values need to be shifted towards the Right, i.e., towards the last value.'
    ''
    'This correction assumes that the data are band-limited (i.e. there is no meaningful information present in the data at a frequency higher than that of the Nyquist). This assumption is support by the study of Josephs et al (1997, Human Brain Mapping)/* \cite{josephs1997event} */ that obtained event-related data at an effective TR of 166 msecs. No physio-logical signal change was present at frequencies higher than our typical Nyquist (0.25 HZ).'
    ''
    'When using the slice timing correction it is very important that you input the correct slice order, and if there is any uncertainty then users are encouraged to work with their physicist to determine the actual slice acquisition order.'
    ''
    'One can also consider augmenting the model by including the temporal derivative in the informed basis set instead of slice timing, which can account for +/- 1 second of changes in timing.'
    ''
    'Written by Darren Gitelman at Northwestern U., 1998.  Based (in large part) on ACQCORRECT.PRO from Geoff Aguirre and Eric Zarahn at U. Penn.'
    }';
st.prog = @spm_run_st;
st.vout = @vout;
st.modality = {'FMRI'};

%==========================================================================
function dep = vout(job)
for k=1:numel(job.scans)
    dep(k)            = cfg_dep;
    dep(k).sname      = sprintf('Slice Timing Corr. Images (Sess %d)', k);
    dep(k).src_output = substruct('()',{k}, '.','files');
    dep(k).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
end

author: spm


text:
function f = fast_specific_rhs(x,y,nel)
%zero_rhs   zero RHS forcing function
%   f = specific_rhs(x,y,nel)
%   input
%          x          x coordinate vector
%          y          y coordinate vector
%          nel        number of elements
%   IFISS function: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
f=zeros(nel,1);
return


text:
function A = spblockdiags(B,d,m,n)
% BLOCKDIAGS : Create sparse block diagonal matrices.
%
% A = spblockdiags(B,d,m,n).
%
%   Blockdiags, which generalizes the function "spdiags",
%   produces a sparse matrix with specified block diagonals.
%
%   A is an m*k-by-n*k matrix, or an m-by-n matrix of k-by-k blocks.
%       The nonzero blocks of A are located on p block diagonals.
%   B is a min(m,n)*k-by-p*k matrix whose k-by-k block columns
%       are the block diagonals of A.
%   (Alternatively, B is k-by-p*k, and then A is block Toeplitz.)
%   d is a vector of p integers in the range -m+1 : n-1,
%       specifying which block diagonals in A are to be nonzero.
%   The values of p and k are determined from the dimensions of B and d.
%
%   For k=1 this is exactly the same as A = spdiags(B,d,m,n); see spdiags
%   for examples of use.  For k>1 this is conceptually the same as spdiags,
%   but k-by-k blocks replace matrix elements everywhere.
%
%   For example, the following code sets A to the n^2-by-n^2 matrix of
%   the Laplacian on an n-by-n square grid; the matrix is block tridiagonal,
%   and the nonzero blocks themselves are tridiagonal or the identity.
%
%        a = spblockdiags ([-1 4 -1], -1:1, n, n);
%        I = speye (n, n);
%        A = spblockdiags ([-I a -I], -1:1, n, n);
%
%  John Gilbert, Xerox PARC, 17 April 1991.
%  Viral Shah, UCSB, 12 April 2006.
% Copyright (c) 1990-1996 by Xerox Corporation.  All rights reserved.
% HELP COPYRIGHT for complete copyright and licensing notice.

if nargin ~= 4
   error ('Usage: A = blockdiags(B,d,m,n)');
end

k = length(d);
[nrB,ncB] = size(B);
p = ncB/k;

% Check for reasonable input.

if any(size(m)>1) | any(size(n)>1)
   error ('blockdiags(B,d,m,n): m or n not scalar');
end
if min(size(d)~=1)
   error ('blockdiags(B,d,m,n): d not a vector');
end
if any(rem(size(B),p))
   error ('blockdiags(B,d,m,n): block size does not divide size of B');
end

A = sparse (m*p, n*p);

for i=1:k
  S = spdiags (ones(m,1), d(i), m, n);
  block = B(:, [(i-1)*p+1:i*p]);
  if isscalar(block)
    A = A + block .* S;
  else
    A = A + kron (S, block);
  end
end

return;

repo: REST


text:
		% Incarcarea continutului fisierului
load patrmag.dat
		% Alocarea variabilei B continutul fisierului
B=patrmag


text:
function ts = spcalib(varargin)
% function t = spcalib()
% function t = spcalib('auto', n) % autocalibration n time, n=3 by default
% function t = spcalib('reset') % Recover the timing set during the session
% function t = spcalib(newT) % Set the new timing
%
% Return the typical running time of each atomic task (mostly sparse
% access) for the specific computer. Make sure no other program do not run
% at the same time
%
% Author Bruno Luong <brunoluong@yahoo.com>
% Last update: 11-Apil-2009
%              13-Nov-2010, move mexing to build_spidxmex

global TREF;

if nargin>=1
    v1 = varargin{1};
    if isstruct(v1)
        TREF = v1;
    elseif ischar(v1) && ~isempty(strmatch(lower(v1),'auto'))

        ntest = 3;
        if nargin>=2
            v2 = varargin{2};
            if isnumeric(v2) && isscalar(v2) && v2>=1
                ntest = floor(v2);
            end
        end

        % Prompt user to save the customize time
        msg='Compile mexfile get/setspvalmex?';
        ButtonName = questdlg(msg, 'SPcalib', 'Yes', 'No', 'No');
        if strcmp(ButtonName,'Yes')
            try
                build_spidxmex();
            catch
                warning('Could not compile mex file');
            end
        end

        msg='Make sure the computer is not loaded by other applications';
        questdlg(msg, 'SPcalib', 'OK', 'OK');
        drawnow;

        fprintf('Auto calibration in progess, please be patient...\n');

        n=1e5;
        nz=1e6;
        i=ceil(n*rand(1,nz));
        j=ceil(n*rand(1,nz));
        iminmax=min(maxlinind([n n]),n*n);
        ilin=ceil(iminmax*rand(1,nz));
        s=rand(size(i));

        % Allocate time arrays
        tfind=zeros(1,ntest);
        tcreate=zeros(1,ntest);
        treadlin=zeros(1,ntest);
        treadfor=zeros(1,ntest);
        twritefor=zeros(1,ntest);
        tadd=zeros(1,ntest);
        tsetval=zeros(1,ntest);
        tgetval=zeros(1,ntest);
        tinsert=zeros(1,ntest);
        tremove=zeros(1,ntest);
        tismember=zeros(1,ntest);

        for l=1:ntest

            % ismember
            tic
            ismember([i j],[i j],'rows');
            tismember(l)=toc/numel(i);
            fprintf('.');

            % create sparse
            tic
            S=sparse(i,j,s,n,n);
            tcreate(l)=toc/nnz(S);
            fprintf('.')

            % find
            tic
            [I J v]=find(S); %#ok
            tfind(l)=toc/nnz(S);
            fprintf('.');

            % read with linear indexes
            tic
            v=S(ilin); %#ok
            treadlin(l)=toc/numel(ilin);
            fprintf('.');

            % read with for loop
            tic
            nloop=min(length(i),Inf);
            for k=1:nloop
                s(k) = S(i(k),j(k));
            end
            treadfor(l)=toc/nloop;
            fprintf('.');

            % write with for loop
            tic
            nloop=min(length(i),10000);
            for k=1:nloop
                S(i(k),j(k)) = s(k);
            end
            twritefor(l)=toc/nloop;
            fprintf('.');

            % add large sparse matrices to a small one
            S=sparse([],[],[],n,n);
            tic
            Sadd = sparse(i,j,s,n,n);
            S = S + Sadd;
            tadd(l)=toc/nnz(S);
            fprintf('.');

            % getspvalmex
            [i j]=find(S);
            tic
            s = getspvalmex(S,i(:),j(:));
            tgetval(l)=toc/numel(i)/max(log2(nnz(S)),7);
            fprintf('.');

            % setspvalmex
            [i j s]=find(S);
            tic
            S = setspvalmex(S,i(:),j(:),s(:));
            tsetval(l)=toc/numel(s)/max(log2(nnz(S)),7);
            fprintf('.');

            % insert new elements
            S=sparse([],[],[],n,n,nz);
            tic
            nloop=min(length(i),10000);
            for k=1:nloop
                S(i(k),j(k)) = s(k);
            end
            tinsert(l)=toc/nloop;
            fprintf('.');

            % remove existing elements
            % Use the same n-loop
            % nloop=min(length(i),10000);
            tic
            for k=1:nloop
                S(i(k),j(k)) = 0;
            end
            tremove(l)=toc/nloop;
            fprintf('.');
        end

        fprintf('\n');

        clear S i j s

        TREF = struct('factory', false, ...
            'tcreate', median(tcreate), ...
            'tfind', median(tfind), ...
            'treadlin',  median(treadlin), ...
            'tismember', median(tismember), ...
            'treadfor',  median(treadfor), ...
            'twritefor',  median(twritefor), ...,
            'tadd',  median(tadd), ...
            'tgetval',  median(tgetval), ...
            'tsetval',  median(tsetval), ...
            'tinsert',  median(tinsert), ...
            'tremove',  median(tremove));

        fprintf('\n');

        disp(TREF);

        % Prompt user to save the customize time
        msg='Save the calibration permanently?';
        ButtonName = questdlg(msg, 'SPcalib', 'Yes', 'No', 'No');

        % Write to a file
        if strcmp(ButtonName,'Yes')
            path = mfilename('fullpath');
            path = fileparts(path);
            fid = fopen([path filesep 'defaultTref.m'], 'w');
            if fid>0
                contains = {
                    'function tref = defaultTref'
                    '%'
                    '% Return the typical running time of each task'
                    '% This must be customized for each specific computer'
                    '% The structure here is the same as output returned by calling'
                    '% SPCALIB(''auto'')'
                    '%'
                    '% Do not change here'
                    '% Reference factory, obtained from 2007 Sony laptop '
                    '%  tref = struct( ...'
                    '%      ''factory'', true, ...'
                    '%      ''tcreate'', 7.8671e-00, ...'
                    '%      ''tfind'', 5.088857e-008, ...'
                    '%      ''treadlin'', 1.1566e-007, ...'
                    '%      ''tismember'', 1.716785e-007, ...'
                    '%      ''treadfor'', 3.4761e-006, ...'
                    '%      ''twritefor'', 7.439304e-006, ...'
                    '%      ''tadd'', 1.0071e-007, ...'
                    '%      ''tgetval'', 1.2649e-008, ...'
                    '%      ''tsetval'', 1.1009e-008, ...'
                    '%      ''tinsert'', 6.796743e-005, ...'
                    '%      ''tremove'', 5.114582e-005);'
                    sprintf('%% Create date: %s', datestr(now))
                    ''
                    ' tref = struct( ...'
                            '     ''factory'', false, ...'
                    sprintf('     ''tcreate'', %e, ...', TREF.tcreate)
                    sprintf('     ''tfind'', %e, ...', TREF.tfind)
                    sprintf('     ''treadlin'', %e, ...', TREF.treadlin)
                    sprintf('     ''tismember'', %e, ...', TREF.tismember)
                    sprintf('     ''treadfor'', %e, ...', TREF.treadfor)
                    sprintf('     ''twritefor'', %e, ...', TREF.twritefor)
                    sprintf('     ''tadd'', %e, ...', TREF.tadd)
                    sprintf('     ''tgetval'', %e, ...', TREF.tgetval)
                    sprintf('     ''tsetval'', %e, ...', TREF.tsetval)
                    sprintf('     ''tinsert'', %e, ...', TREF.tinsert)
                    sprintf('     ''tremove'', %e);', TREF.tremove)
                    ''
                    'end'
                    };
                count = fprintf(fid, '%s\n',contains{:});
                if count<sum(cellfun(@length,contains)+1)
                    error('Could not save defaultTref.m');
                end
            else
                error('Could not save defaultTref.m');
            end
            fclose(fid);
            drawnow;
        end

    else
        TREF = [];
    end

end

% Factory tune
if isempty(TREF)
    TREF = defaultTref();
end

ts = TREF;

end


text:
function plotstations(ax, options)
    %plotstations plots stations on specified axes
    % (Not fully IMPLEMENTED)
    % plotstatiosn(ax) plots the default station list, cropped to the current window
    % to the axis AX.
    % to the specified axes. STATIONS is a structure or table
    % with fields Network, Name, Latitude, Longitude,
    % additionally, it should have Elevation, StartTime, and EndTime.
    % by default, only plots stations active during main catalog dates.
    %
    % if stations is a file, then it should have a variable `'stations' with the
    % same fields.
    % options would be a struct with filtervalues
    %
    % eg.
    %  options=struct('Network','CH'); % show only stations belonging to network CH
    %  options=struct('Network',{'CH','DK'}) % show only stations belonging either to CH or DK
    %  valid options are:
    %  Network, Station, StartBefore, EndAfter, Marker, MarkerSize, MarkerColor
    %
    % FileName:
    % path&name of a station file meeting above criteria. Samples that come with Zmap:
    %  resrc/features/stations_irisall_20170714.mat       %unfiltered IRIS station list
    %  or  resrc/features/stations_sedall_20170714.mat    %unfiltered SED station list

    % TODO Maybe these should be actually go into a get/load stations file.
    % TODO Provide ability to update the stations file?
    % TODO set default network?

    % all swiss stations:
    % http://eida.ethz.ch/fdsnws/station/1/query?format=text&level=station&nodata=404
    unimplemented_error()
    report_this_filefun();

    ZG=ZmapGlobal.Data;
    if ~exist('options','var')
        options=struct();
    end

    if ~isfield(options,'Marker'), options.Marker='s';end
    if ~isfield(options,'MarkerSize'), options.MarkerSize=6;end
    if ~isfield(options,'MarkerEdgeColor'), options.MarkerEdgeColor='gray';end
    if ~isfield(options,'MarkerFaceColor'), options.MarkerFaceColor='gray';end
    if ~isfield(options,'StartBefore')
        %     w   [x   y ]  z   want stations that start some time before y and end some time after x
        options.StartTime=max(ZG.primeCatalog.Date);
    end
    if ~isfield(options,'EndAfter')
        options.EndTime=min(ZG.primeCatalog.Date);
    end
    stations=ZG.features('stations');
    if isempty(stations)
        stations.load();
    end

    % filter the stations down to something managable
    % WARNING: If same station has multiple periods at same place, it gets plotted twice.
    yl = ylim(ax);
    xl = xlim(ax);
    idx = stations.Latitude <= yl(2) & ...
        stations.Latitude >= yl(1) & ...
        stations.Longitude <= xl(2) & ...
        stations.Longitude >= xl(1); %doesn't account for dateline
    if isa(stations,'Table')
        if ismember('StartTime',stations.Properties.VariableNames) && ~isempty(options.StartBefore)
            idx=idx & stations.StartTime <= options.StartBefore;
        end
        if ismember('EndTime',stations.Properties.VariableNames) && ~isempty(options.EndAfter)
            idx=idx & stations.EndTime >= options.EndAfter;
        end
        if isfield(options.Network) && ~isempty(options.Network)
            idx=idx & ismember(stations.Network, options.Network);
        end
        stas = stations(idx,:);
    else % struct
        if isfield(stations,'StartTime') && ~isempty(options.StartBefore)
            idx=idx & stations.StartTime <= options.StartBefore;
        end
        if isfield(stations,'EndTime') && ~isempty(options.EndAfter)
            idx=idx & stations.EndTime >= options.EndAfter;
        end
        if isfield(options.Network) && ~isempty(options.Network)
            idx=idx & ismember(stations.Value.Network, options.Network);
        end
        stas = stations(idx);
    end

    showNames = ~isfield(options,'ShowNames') || options.ShowNames;
    options = rmfield(options, {'StartTime','EndTime','Network', 'ShowNames'});

    set(gca,'NextPlot','add')
    dx = abs(s1_east-s2_west)/130;
    pl = plot(ax, stas.Longitude,stas.Latitude,options.Marker);
    set(pl,'LineWidth',1,options);
    if showNames
        te1 = text(stas.Longitude+dx,lastas.Latitudet,char(stas.Name),'clipping','on');
        set(te1,'FontWeight','bold','Color','k','FontSize',9);
    end
    drawnow
    set(gca,'NextPlot','replace')
end


text:
function [state,Energy,TotalEnergy1MCS,ColorMatrix] = InitializeMatricesQPOTTS(x,MonteCarloSteps,Q);

state=zeros(size(x,1),size(x,1));
Energy=zeros(size(x,1),size(x,1),MonteCarloSteps);
TotalEnergy1MCS=zeros(1,MonteCarloSteps);
ColorMatrix=rand(Q,3);

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34784-monte-carlo-simulation-of-two-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 2D square-lattice - microstructure/InitializeMatricesQPOTTS.m


text:
function out = TreatData(data,treat,vnames,dates,DatesOpt)
% =======================================================================
% Treat data.
% =======================================================================
% [nobs, dates] = CountDate(fo_year,lo_year,frequency,fo_period,lo_period)
% -----------------------------------------------------------------------
% INPUT
%	- data: structure where the the data is stored
%	- treat: vector (of length N) with type of treatment
%	- vnames: vector (of length N) with variable names
% -----------------------------------------------------------------------
% OUTPUT
%	- out: matrix of treated data
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

% Retrieve some parameters
nobs = DatesOpt.nobs;
fo_year = DatesOpt.fo_year;
lo_year = DatesOpt.lo_year;
frequency = DatesOpt.frequency;
fo_period = DatesOpt.fo_period;
lo_period = DatesOpt.lo_period;

% Get first
if strcmp(frequency,'y')
    aux = num2str(fo_year);
else
    aux = [num2str(fo_year) frequency num2str(fo_period)];
end
fo = find(strcmp(aux,dates));

% Get last
if strcmp(frequency,'y')
    aux = num2str(lo_year);
else
    aux = [num2str(lo_year) frequency num2str(lo_period)];
end
lo = find(strcmp(aux,dates));

% Initialize matrix
nvar = length(vnames);
out = nan(nobs,nvar);

for ii=1:nvar
    % No treatment
    if treat(ii)==0
        out(:,ii) = data.(vnames{ii})(fo:lo);
    % Log
    elseif treat(ii)==1
        out(:,ii) = log(data.(vnames{ii})(fo:lo));
    % Log-diff
    elseif treat(ii)==2
        out(:,ii) = XoX(data.(vnames{ii})(fo:lo),1,'logdiff');
    % Diff
    elseif treat(ii)==3
        out(:,ii) = XoX(data.(vnames{ii})(fo:lo),1,'diff');
    end
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
% fig2: anat averages from paired regressors: sensory/motor correlations
% plots a pair of regressors in a graded double colormap, saves as tiff
% stacks, and also saves anat average plot as tiff file.
% (compare to fig2A_correlation_1reg_with_hist for more details)

clear all; close all; clc

outputDir = GetOutputDataDir;

%% Init load
hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

%% set params
thres_reg_const = 0.5;
% M_reg_thres = {0.5,0.5,0.5};

isKeepTopPrct = 0; % init, can override
isMotorseed = 0; % init, can override

for setflag = 9
    % reset
    ResetDisplayParams(hfig);

    switch setflag
        case 1 % Phototaxis
            %             stimrange = 1;
            M_stimrange = GetStimRange('P');
            M_regtypeflag = [1,1,2]; % 0 for VAR, 1 for stim and 2 for motor
            M_reg_name = {'PT','HalfField','PT_SwimLR'};
            M_reg_range = {[3,2],[5,6],[1,3]};
            M_fishrange = {[1:18],[1:7],[1:3,5:18]};

        case 2 % OMR
            %             stimrange = 2;
            M_stimrange = GetStimRange('O');
            M_regtypeflag = [1,1,2]; % 1 for stim and 2 for motor
            M_reg_name = {'OMR_FwBw','OMR_LR','OMR_SwimLR'};
            M_reg_range =  {[7,6],[9,8],[1,3]};
            M_fishrange = {[8:18],[8:18],[8:18]};

        case 3 % Looming
            %             stimrange = 5;
            M_stimrange = GetStimRange('L');
            M_regtypeflag = [1,2]; % 1 for stim and 2 for motor
            M_reg_name = {'Loom_LR','Loom_SwimLR'};
            M_reg_range =  {[11,12],[1,3]};
            M_fishrange = {[9:15,17:18],[9:15,17:18]};

        case 4 % Dark Flash (black-white)
            %             stimrange = 3;
            M_stimrange = GetStimRange('D');
            M_regtypeflag = [1,2]; % 1 for stim and 2 for motor
            M_reg_name = {'DF_BW','DF_SwimLR'};
            M_reg_range =  {[1,4],[1,3]};
            M_fishrange = {[1:5,12:15,17:18],[12:15,17:18]}; % up till 7/10
%             M_fishrange = {[12:15,17:18],[12:15,17:18]}; % up till 7/10

%             % from VAR:

        case 5 % ABN
            M_regtypeflag = [0];
            M_stimrange = GetStimRange();%('2');
%             M_reg_name = {'ABN_top1%_defstimrange'};
            M_reg_name = {'ABN_reg0.5_defstimrange'};
            M_clus_range = {[12,1]};
            range = GetFishRange('e');%[1:8,11,12,14:17]
            M_fishrange = {range};%{[1:12,14:18]};

%             isKeepTopPrct = 1;
%             prct_const = 1;

        case 6 % Fw
            M_regtypeflag = [0];
            M_stimrange = GetStimRange();%('2');
            M_reg_name = {'Fw_seed_reg2%_defS'};
%             M_reg_name = {'Fw_reg0.5_defS'};
            M_clus_range = {[11,4]};
            M_fishrange = {[1:3,5:18]};

            isKeepTopPrct = 1;
            prct_const = 2;

        case 7 % HBO 4 stripes
            M_regtypeflag = [0,0];
            M_stimrange = GetStimRange();%('2');
            M_reg_name = {'HBO-L2_reg1%_tRes_defS','HBO-R2_reg1%_tRes_defS'};
            M_clus_range = {[10,2],[10,3]};
            M_fishrange = {[1:3,5:18],[1:3,5:18]};

            isKeepTopPrct = 1;
            prct_const = 1;
            setappdata(hfig,'isTrialRes',1);

        case 8 % HBO 2 halves
            M_regtypeflag = [0,0];
            M_stimrange = GetStimRange();%('2');
            M_reg_name = {'HBO2_reg1%_tRes_defS'};
            M_clus_range = {[10,5]};
            M_fishrange = {[1:3,5:18]};

            isKeepTopPrct = 1;
            prct_const = 1;
            setappdata(hfig,'isTrialRes',1);

        case 9 % spontaneous
            M_stimrange = GetStimRange('S'); % spontaneous
            M_regtypeflag = [2]; % 1 for stim and 2 for motor
            M_reg_name = {'spt_SwimLR_notmotorseed_top2%'};
            M_reg_range =  {[1,3]}; % {[1,2]} for ismotorseed = 1
            M_fishrange = {[8:15,17:18]};

            isMotorseed = 0; % 1
            isKeepTopPrct = 1;
            prct_const = 2;
    end

    n_reg = length(M_reg_name);
    M_fishrange_im = M_fishrange;

    range = 1:18; % init,M_fishrange below specifies actual range #oldcode
    IM_full = cell(n_reg,18);

    %%

%     [M_fishrange_im,fishrange_load] = CheckIfLoadFish(M_fishrange,M_ClusterIDs,M_stimrange);
% [cIX_load,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,M_ClusterIDs{i_set},M_stimrange{i_set});
    for i_fish = range
        if ismember(i_fish,cell2mat(M_fishrange))
            ClusterIDs = [1,1];%GetClusterIDs('all');
            stimrange = M_stimrange{i_fish};
            % load fish data
            if isempty(stimrange)
                [cIX_load,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs);
            else
                [cIX_load,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs,stimrange);
            end
        end

        %% Load stim/motor
        for i_set = 1:n_reg

            if ~ismember(i_fish,M_fishrange{i_set})
                continue;
            end

            reg_thres = thres_reg_const; %M_reg_thres{i_set};

            %% Get Regressors
            if M_regtypeflag(i_set)==0
                ClusterIDs = M_clus_range{i_set};%[12,1];% GetClusterIDs('all');
                [cIX,gIX] = LoadCluster_Direct(i_fish,ClusterIDs(1),ClusterIDs(2));
                M = UpdateIndices_Manual(hfig,cIX,gIX);
                Reg = FindClustermeans(gIX,M);

            elseif M_regtypeflag(i_set)==1
                fishset = getappdata(hfig,'fishset');
                [~,names,regressors] = GetStimRegressor(stim,fishset,i_fish);

                reg_range = M_reg_range{i_set}; % left/right pair
                Reg = regressors(reg_range,:);

            elseif M_regtypeflag(i_set)==2
%                 isMotorseed = 0;
                setappdata(hfig,'isMotorseed',isMotorseed);
                [~,~,behavior] = UpdateTimeIndex(hfig);
                [~,names,regressors] = GetMotorRegressor(behavior,i_fish);

                reg_range = M_reg_range{i_set}; % left/right pair
                Reg = regressors(reg_range,:);
            end

            %% Regression
            % code adapted from 'best regressor regression' code 'AllRegsRegression'

            if isempty(Reg)
                M_fishrange_im{i_set} = setdiff(M_fishrange_im{i_set} ,i_fish);
                continue;
            end

            Corr = corr(Reg',M_0');

            if isKeepTopPrct
                % keep best regression only
                [Corr_rows,corr_max] = ConvertCorrToBestRegRows(Corr);

                % top x percentile
                nCells_total = size(M_0,1);
%                 prct_const = 1;
                [CIX,RegThres] = ChooseTopPercentOfFish(nCells_total,prct_const,Corr_rows);

                if length(CIX)==1
                    cIX = CIX{1};
                    % get map color
                    reg_thres = 0.25;
                    gIX = MapXto1Dcolormap(corr_max(cIX),[reg_thres,1],64);
                else
                    cIX1 = CIX{1};
                    cIX2 = CIX{2};

                    % get map color
                    clrIX1 = MapXto1Dcolormap(corr_max(cIX1),[reg_thres,1],64);
                    clrIX2 = MapXto1Dcolormap(corr_max(cIX2),[reg_thres,1],64);
                    %     clrIX1 = MapXto1Dcolormap(corr_max(cIX1),[reg_thres1,1],64);
                    %     clrIX2 = MapXto1Dcolormap(corr_max(cIX2),[reg_thres2,1],64);

                    cIX = [cIX1;cIX2];
                    %     if isempty(cIX)
                    %         M_fishrange_im{i_set} = setdiff(M_fishrange_im{i_set} ,i_fish); %#ok<SAGROW>
                    %         continue;
                    %     end

                    clrIX = [clrIX1;clrIX2];
                    gIX_offset = [ones(size(cIX1));2*ones(size(cIX2))];
                    gIX = clrIX+(gIX_offset-1)*64;
                    %     gIX = [clrIX1;64+clrIX2];
                    numK = length(unique(gIX));
                end

            else
                %%
                [corr_max,IX_regtype] = max(Corr,[],1);
                cIX = find(corr_max>reg_thres)';
                gIX_offset = IX_regtype(cIX)';
                clrIX = MapXto1Dcolormap(corr_max(cIX),[reg_thres,1],64);

                if isempty(cIX)
                    M_fishrange_im{i_set} = setdiff(M_fishrange_im{i_set} ,i_fish);
                    continue;
                end

                gIX = clrIX+(gIX_offset-1)*64;
                numK = length(unique(gIX));
            end

            %% make double colormap
            clr1 = [1,0,0];
            clr1_ = [0.7,0.5,0.5];
            clr2 = [0,1,1];
            clr2_ = [0.5,0.7,0.7];
            numC = 64;
            clrmap1 = Make1DColormap([clr1_;clr1],numC);
            clrmap2 = Make1DColormap([clr2_;clr2],numC);
            clrmap = [clrmap1;clrmap2];

            %% make figure
            I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX,clrmap);
            [h,im_full] = DrawCellsOnAnat(I);

            %         % add 2 colorbars
            % %         AddColorbarToAnat(clrmap,cmin,cmax)
            % %         colormap(clrmap1);
            % %         caxis([reg_thres,1])
            % %         colorbar('Location','manual','Position',[0.8,0.7,0.05,0.15],'Units','normalized')
            %         ax = axes('Position',[0.75,0.8,0.05,0.15],'Units','normalized');
            %         DrawCustomColorbar(clrmap1,[reg_thres,1],2,ax);
            %
            %         ax = axes('Position',[0.9,0.8,0.05,0.15],'Units','normalized');
            %         DrawCustomColorbar(clrmap2,[reg_thres,1],2,ax);
            %
            %% save figure
            close(h);
            IM_full{i_set,i_fish} = im_full;

        end
    end

    %% save as tiff stack
    for i_set = 1:n_reg
        range_im = M_fishrange_im{i_set};
        tiffdir = fullfile(outputDir,[M_reg_name{i_set},'_allfish.tiff']);
        IM = IM_full(i_set,range_im);

        SaveImToTiffStack(IM,tiffdir);
    end

    %% [for later] plot from tiff stack
    isPlotfromtiffstack = 0;

    if isPlotfromtiffstack
        IM_full = cell(n_reg,18);
        for i_set = 1:n_reg
            %% get tiff-stack path
            tiffdir = fullfile(outputDir,[M_reg_name{i_set},'_allfish.tiff']);

            %% load
            range_load = 1:17;
            for i_fish = range_load
                im = double(imread(tiffdir,i_fish))./255;
                IM_full{i_set,i_fish} = im;
%                 IM_full{i_set,i_fish} = im(317:1236,1:621,:);
            end
        end
    end

    %% make average plots
    % M_k_scale = {1,1.5,1};
    % M_k_contrast = {1.2,1.5,1.2};

    for i_set = 1:n_reg
        range_im = M_fishrange_im{i_set};%[1:3,5:7];%[1:3,5:18];
        cellarray = IM_full(i_set,range_im);

        % adjust params for visualization
        k_scale = 0.7;%1/1.5;%M_k_scale{i_set};
        k_contrast = 1.1;%M_k_contrast{i_set};

        [h_anat,im_avr] = AverageAnatPlot(cellarray,k_contrast,k_scale);

        tiffdir = fullfile(outputDir,[M_reg_name{i_set},'_avr.tiff']);
        imwrite(im_avr, tiffdir, 'compression','none','writemode','overwrite');
    end

end

author: xiuyechen

repo: FishExplorer


text:
function t = mytrace(A,B)

% This compute the trace inner product of two matrices A, B, of the same
% dimension.

t=sum(sum(A.*B));

author: zzstefan

repo: BrainNetClass


text:
%DEMO_PHASEPLOT  Give demos of nice phaseplots
%
%   This script creates a synthetic signal and then uses |phaseplot| on it,
%   using several of the possible options.
%
%   For real-life signal only small parts should be analyzed. In the chosen
%   demo the fundamental frequency of the speaker can be nicely seen.
%
%   .. figure::
%
%      Synthetic signal
%
%      Compare this to the pictures in reference 2 and 3. In
%      the first two figures a synthetic signal is analyzed. It consists of a
%      sinusoid, a small Delta peak, a periodic triangular function and a
%      Gaussian. In the time-invariant version in the first part the periodicity
%      of the sinusoid can be nicely seen also in the phase coefficients. Also
%      the points of discontinuities can be seen as asymptotic lines approached
%      by parabolic shapes. In the third part both properties, periodicity and
%      discontinuities can be nicely seen. A comparison to the spectogram shows
%      that the rectangular part in the middle of the signal can be seen by the
%      phase plot, but not by the spectogram.
%
%      In the frequency-invariant version, the fundamental frequency of the
%      sinusoid can still be guessed as the position of an horizontal
%      asymptotic line.
%
%   .. figure::
%
%      Synthetic signal, thresholded.
%
%      This figure shows the same as Figure 1, except that values with low
%      magnitude has been removed.
%
%   .. figure::
%
%      Speech signal.
%
%      The figure shows a part of the 'linus' signal. The fundamental
%      frequency of the speaker can be nicely seen.
%
%   References: carmonmultiridge1 Carmona98practical gross1

disp('Type "help demo_phaseplot" to see a description of how this demo works.');

tt=0:98;
f1=sin(2*pi*tt/33); % sinusoid

f2=zeros(1,100);
f2(50)=1; % delta-like

f3=fftshift(firwin('tria',32)).';

f4 = fftshift(pgauss(100)).';
f4 = f4/max(f4);

sig = 0.9*[f1 0 f2 f3 -f3 f3 f4 0 0 0 0];

figure(1);
sgram(sig,'lin','nf');

figure(2);
subplot(3,1,1);
plot(sig);
title('Synthetic signal');

subplot(3,1,2);
phaseplot(sig,'freqinv');
title('Phaseplot of synthetic signal - frequency-invariant phase');

subplot(3,1,3);
phaseplot(sig,'timeinv')
title('Phaseplot of synthetic signal - time-invariant phase');

figure(3);
subplot(3,1,1);
plot(sig);
title('Synthetic signal');

subplot(3,1,2);
phaseplot(sig,'freqinv','thr',0.001)
title('Phaseplot of synthetic signal - thresholded version, freq. inv. phase');

subplot(3,1,3);
phaseplot(sig,'thr',0.001)
title('Phaseplot of synthetic signal - thresholded version, time inv. phase');

figure(4);
f=linus;
f = f(4500:8000);

subplot(3,1,1);
plot(f);
axis tight;
title('Speech signal: linus');

subplot(3,1,2);
phaseplot(f)
title('Phaseplot of linus');

subplot(3,1,3);
phaseplot(f,'thr',.001)
title('Phaseplot of linus - thresholded version');

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%% This is for testing the Plotting functions in the robotics Toolbox
function tests = TransformationsTest
  tests = functiontests(localfunctions);
  clc
end

function teardownOnce(tc)
    close all
end

%    tranimate                  - animate a coordinate frame
function tranimate_test(tc)
    X1 = eye(3,3); X2 = rotx(pi/2);
    tranimate(X1, X2);
    tranimate(X1, X2, 'nsteps', 10);

    clf
    tranimate(X1, X2, 'axis', [-10 10 -20 20 -30 30]);
    ax = gca; v = [ax.XLim ax.YLim ax.ZLim];
    tc.verifyEqual(v, [-10 10 -20 20 -30 30]); % 19b doesnt give

    tranimate(X1, X2, 'noxyz');
    tranimate(X1, X2, 'rgb');
    tranimate(X1, X2, 'retain');
    tranimate(X1, X2, 'fps', 20);

    X1 = eye(4,4); X2 = transl(1,2,3)*trotx(pi/2);
    tranimate(X1, X2);
    tranimate(X1, X2, 'nsteps', 10);

    clf
    tranimate(X1, X2, 'axis', [-10 10 -20 20 -30 30]);
    %v = axis;
    ax = gca; v = [ax.XLim ax.YLim ax.ZLim];
    tc.verifyEqual(v, [-10 10 -20 20 -30 30]);

    tranimate(X1, X2, 'noxyz');
    tranimate(X1, X2, 'rgb');
    tranimate(X1, X2, 'retain');
    tranimate(X1, X2, 'fps', 20);
end

function tranimate2_test(tc)
    X1 = eye(2,2); X2 = rot2(pi/2);
    tranimate2(X1, X2);
    tranimate2(X1, X2, 'nsteps', 10);

    clf
    tranimate2(X1, X2, 'axis', [-10 10 -20 20]);
    v = axis;
    tc.verifyEqual(v, [-10 10 -20 20]);

    tranimate2(X1, X2, 'noxyz');
    tranimate2(X1, X2, 'rgb');
    tranimate2(X1, X2, 'retain');
    tranimate2(X1, X2, 'fps', 20);

    X1 = eye(3,3); X2 = transl2(1,2)*trot2(pi/2);
    tranimate2(X1, X2);
    tranimate2(X1, X2, 'nsteps', 10);

    clf
    tranimate2(X1, X2, 'axis', [-10 10 -20 20]);
    v = axis;
    tc.verifyEqual(v, [-10 10 -20 20]);

    tranimate2(X1, X2, 'noxyz');
    tranimate2(X1, X2, 'rgb');
    tranimate2(X1, X2, 'retain');
    tranimate2(X1, X2, 'fps', 20);
end

%    trplot                     - plot HT as a coordinate frame
function trplot_test(tc)
    %%
    Rt1 = [1.0000         0         0         0
                0    0.8253   -0.5646         0
                0    0.5646    0.8253         0
                0         0         0    1.0000];
    trplot(Rt1);
    clf
    h = trplot(Rt1);
    trplot(Rt1, 'handle', h);
    trplot(Rt1, 'color', 'r');
    trplot(Rt1, 'color', [1 0 1]);
    trplot(Rt1, 'noaxes');
    trplot(Rt1, 'frame', 'bob');
    trplot(Rt1, 'frame', 'A', 'text_opts', {'FontSize', 10, 'FontWeight', 'bold'})
    trplot(Rt1, 'view', [10 20]);
    trplot(Rt1, 'arrow')
    trplot(Rt1, '3d')
    trplot(Rt1, '3d', 'anaglyph', 'mo')
    trplot(Rt1, '3d', 'dispar', 0.3);
end

 function trplot2_test(tc)
     %%
    Rt1 = [
    0.9553   -0.2955    1.0000
    0.2955    0.9553    2.0000
         0         0    1.0000];
    trplot2(Rt1);
    clf
    h = trplot2(Rt1);
    trplot2(Rt1, 'handle', h);
    h = trplot2(Rt1, 'handle', h);
    trplot2(Rt1, 'color', 'r');
    trplot2(Rt1, 'color', [1 0 1]);
    trplot2(Rt1, 'noaxes');
    trplot2(Rt1, 'frame', 'bob');
    trplot2(Rt1, 'frame', 'A', 'text_opts', {'FontSize', 10, 'FontWeight', 'bold'})
    trplot2(Rt1, 'view', [10 20]);
    trplot2(Rt1, 'arrow')
    trplot2(Rt1, 'framelabel', 'A')

    clf
    a = gca;
    trplot2(Rt1, 'axhandle', a);

    h = trplot2(Rt1);
    tc.verifyWarning( @() trplot2(h, Rt1), 'SMTB:trplot2:deprecated')
end

%    plot2                      - plot trajectory
function plot2_test(tc)
    %%
        th = [0:0.1:pi];
        sth = sin(th);
        p = [th' sth'];
        plot2(p);

        p = [th' sth' sth'];
        plot2(p);
end

%    plot_box                   - draw a box
function plot_box_test(tc)
        plot_box(1,1,5,5,'b');
        plot_box(2,2,4,4,'g');
end

%    plot_circle                - draw a circle
function plot_circle_test(tc)
    plot_circle([1 2],2,'g');
    plot_circle([1 2],2,'fillcolor', 'g');
    plot_circle([1 2],2,'fillcolor', 'g', 'alpha', 0.5);
    plot_circle([1 2],2,'edgecolor', 'b');
    plot_circle([1 2],2,'fillcolor', 'g', 'edgecolor', 'b');
end

%    plot_ellipse               - draw an ellipse
function plot_ellipse_test(tc)
    %%
    %2d

    C = diag([1,4]);
    plot_ellipse(C,[2 3],'r');
    plot_ellipse(C,[2 3],'fillcolor', 'g');
    plot_ellipse(C,[2 3],'fillcolor', 'g', 'alpha', 0.5);
    plot_ellipse(C,[2 3],'edgecolor', 'g');
    plot_ellipse(C,[2 3],'fillcolor', 'g', 'edgecolor', 'b');

    % with 3d centre
    plot_ellipse([1 0; 0 4],[2 3 1],'r');

    %3d

    C = diag([1,4,2]);
    plot_ellipse(C,[2 3 1],'r');
    plot_ellipse(C,[2 3 1],'fillcolor', 'g');
    plot_ellipse(C,[2 3 1],'fillcolor', 'g', 'alpha', 0.5);
    plot_ellipse(C,[2 3 1],'edgecolor', 'g');
    plot_ellipse(C,[2 3 1],'fillcolor', 'g', 'edgecolor', 'b');
end

%    plot_homline               - plot homogeneous line
function plot_homline_test(tc)
    plot_homline([1 2 3]','y');
end

%    plot_point                 - plot points
function plot_point_test(tc)
    plot_point([1; 2]);
end

%    plot_poly                  - plot polygon
function plot_poly_test(tc)
    p = [1 2 2 1; 1 1 2 2];
    plot_poly(p,'g');
end

%    plot_sphere                - draw a sphere
function plot_sphere_test(tc)
     plot_sphere([1 2 3],5,'r');
end

function xyzlabel_test(tc)

    plot3(1, 2, 3)
    xyzlabel()
end

function plotvol_test(tc)

    % 2D cases
    plotvol([2 3])
    tc.verifyEqual(axis, [-2 2 -3 3]);
    plotvol([2 3 4 5])
    tc.verifyEqual(axis, [2 3 4 5]);

    % 3D cases
    plotvol(3)
    ax = gca; v = [ax.XLim ax.YLim ax.ZLim];
    tc.verifyEqual(v, 3*[-1 1 -1 1 -1 1]);

    plotvol([1 2 3]);
    ax = gca; v = [ax.XLim ax.YLim ax.ZLim];
    tc.verifyEqual(v, [-1 1 -2 2 -3 3]);

    plotvol([1 2 3 4 5 6]);
    ax = gca; v = [ax.XLim ax.YLim ax.ZLim];
    tc.verifyEqual(v, [1 2 3 4 5 6]);
end

author: petercorke

repo: spatialmath-matlab


text:
function out = spm_shoot_warp(job)
% Register images with template
% format spm_shoot_warp(job)
% Fields of job:
%     job.images{1} first set of images (eg rc1*.nii)
%     job.images{2} second set of images (eg rc2*.nii)
%     etc
%     job.templates template files
% Other settings are defined in spm_shoot_defaults.m
%
% The outputs are flow fields (v*.nii), deformation fields (y*.nii) and
% Jacobian determinants (j*.nii)
%_______________________________________________________________________
% Copyright (C) Wellcome Trust Centre for Neuroimaging (2009)

% John Ashburner
% $Id: spm_shoot_warp.m 7389 2018-08-06 13:35:48Z john $

%_______________________________________________________________________
d       = spm_shoot_defaults;
tname   = d.tname;   % Base file name for templates

cyc_its = d.cyc_its; % No. multigrid cycles and inerations
sched   = d.sched;   % Schedule for coarse to fine
nits    = numel(sched)-1;
rparam  = d.rparam;  % Regularisation parameters for deformation
eul_its = d.eul_its; % Start with fewer steps
scale   = d.scale;   % Fraction of Gauss-Newton update step to use

bs_args = d.bs_args; % B-spline settings for interpolation
%_______________________________________________________________________

spm_diffeo('boundary',0); % Set boundary condition

% Sort out handles to images
n1 = numel(job.images);
n2 = numel(job.images{1});
NF = struct('NI',[],'vn',[1 1]);
NF(n1,n2) = struct('NI',[],'vn',[1 1]);

% Pick out individual volumes within NIfTI files
for i=1:n1
    if numel(job.images{i}) ~= n2
        error('Incompatible number of images');
    end
    for j=1:n2
        [pth,nam,ext,num] = spm_fileparts(job.images{i}{j});
        NF(i,j).NI        = nifti(fullfile(pth,[nam ext]));
        num               = [str2num(num) 1 1];
        NF(i,j).vn        = num(1:2);
    end
end

dm = [size(NF(1,1).NI.dat) 1];
dm = dm(1:3);

% Sort out which template for each iteration
tmpl_no = round(((1:nits)-1)/(nits-1)*(numel(job.templates)-0.51))+1;

ok = true(n2,1);

NU     = nifti;
NU(n2) = nifti;
NY     = nifti;
NY(n2) = nifti;
NJ     = nifti;
NJ(n2) = nifti;

for i=1:n2 % Loop over subjects

    % Generate files for flow fields, deformations and Jacobian determinants.
    [pth,nam,ext]   = fileparts(NF(1,i).NI.dat.fname);
    if ~isempty(tname), nam = [nam '_' tname]; end
    offs  = 352;

    NU(i) = nifti;
    NU(i).dat = file_array(fullfile(pth,['v_' nam '.nii']),[dm 1 3], 'float32-le', offs, 1, 0);
    NU(i).descrip = sprintf('Velocity (%d %.4g %.4g %.4g)', rparam(1), rparam(2), rparam(3), rparam(4));
    NU(i).mat     = NF(1,i).NI.mat;
    NU(i).mat0    = NF(1,i).NI.mat0;
    create(NU(i)); NU(i).dat(:,:,:,:,:) = 0;

    NY(i) = nifti;
    NY(i).dat = file_array(fullfile(pth,['y_' nam '.nii']),[dm 1 3], 'float32-le', offs, 1, 0);
    NY(i).descrip = 'Deformation (templ. to. ind.)';
    NY(i).mat     = NF(1,i).NI.mat;
    create(NY(i)); NY(i).dat(:,:,:,:,:) = reshape(affind(spm_diffeo('Exp',zeros([dm,3],'single'),[0 1]),NU(i).mat0),[dm,1,3]);

    NJ(i) = nifti;
    NJ(i).dat = file_array(fullfile(pth,['j_' nam '.nii']),[dm 1 1], 'float32-le', offs, 1, 0);
    NJ(i).descrip = 'Jacobian det (templ. to. ind.)';
    NJ(i).mat     = NF(1,i).NI.mat;
    create(NJ(i)); NJ(i).dat(:,:,:)     = 1;

    drawnow
end

for i=1:n2 % Loop over subjects. Can replace FOR with PARFOR.

    % Load image data for this subject
    f  = loadimage(NF(:,i));
    u  = squeeze(single(NU(i).dat(:,:,:,:,:)));
    y  = affind(squeeze(single(NY(i).dat(:,:,:,:,:))),inv(NU(i).mat0));
    dt = squeeze(single(NJ(i).dat(:,:,:)));

    % Re-load first template
    [g,vx] = load_template(job.templates{tmpl_no(1)}, n1, bs_args);

    % The actual work
    for it=1:nits

        % Load template appropriate for this iteration
        if (it>1) && (tmpl_no(it)~=tmpl_no(it-1))
            [g,vx] = load_template(job.templates{tmpl_no(it)}, n1, bs_args);
        end

        % More regularisation in the early iterations, as well as a
        % a less accurate approximation in the integration.
        prm      = [vx, rparam*sched(it+1)*prod(vx)];
        int_args = [eul_its(it), cyc_its];
        drawnow

        fprintf(' %-5d %-3d\t| ',i,it);

        % Gauss-Newton iteration to re-estimate deformations for this subject
        u      = spm_shoot_update(g,f,u,y,dt,prm,bs_args,scale); drawnow
        [y,dt] = defdet(u,prm,int_args);

        if any(~isfinite(dt(:)) | dt(:)>100 | dt(:)<1/100)
            ok(i) = false;
            fprintf('Problem with %s (dets: %g .. %g)\n', NU(i).dat.fname, min(dt(:)), max(dt(:)));
           %clear dt
            break
        end

        drawnow
        NU(i).dat(:,:,:,:,:) = reshape(u,[dm 1 3]);
        NY(i).dat(:,:,:,:,:) = reshape(affind(y,NU(i).mat0),[dm 1 3]);
        NJ(i).dat(:,:,:)     = dt;

    end
    fprintf('\n');

    drawnow
end

if any(~ok)
    fprintf('Problems with:\n');
    for i=find(~ok)'
        fprintf('\t%s\n', NU(i).dat.fname);
    end
end

% Finish off
out.vel = cell(n2,1);
out.def = cell(n2,1);
out.jac = cell(n2,1);
for i=1:n2
    out.vel{i} = NU(i).dat.fname;
    out.def{i} = NY(i).dat.fname;
    out.jac{i} = NJ(i).dat.fname;
end
%=======================================================================

%=======================================================================
function y1 = affind(y0,M)
% Affine transform of deformation
y1 = zeros(size(y0),'single');
for d=1:3
    y1(:,:,:,d) = y0(:,:,:,1)*M(d,1) + y0(:,:,:,2)*M(d,2) + y0(:,:,:,3)*M(d,3) + M(d,4);
end
%=======================================================================

%=======================================================================
function f = loadimage(NF)
n1      = size(NF,1);
f       = cell(n1+1,1);
dm      = [NF(1).NI.dat.dim 1 1 1];
dm      = dm(1:3);
f{n1+1} = ones(dm,'single');
for j=1:n1
    vn      = NF(j,1).vn;
    f{j}    = single(NF(j,1).NI.dat(:,:,:,vn(1),vn(2)));
    msk     = ~isfinite(f{j});
    f{j}(msk) = 0;
    f{n1+1} = f{n1+1} - f{j};
    drawnow
end
f{n1+1}(msk) = 0.00001;
%=======================================================================

%=======================================================================
function [g,vx] = load_template(template, n1, bs_args)
g  = cell(n1+1,1);
NG = nifti(template);
if size(NG.dat,4) < n1+1
    error('Not enough tissues in template (%d < %d+1).', size(NG.dat,4),n1);
end

bg = ones([size(NG.dat,1), size(NG.dat,2), size(NG.dat,3)]);
for j=1:n1
    tmp  = NG.dat(:,:,:,j);
    g{j} = spm_bsplinc(log(tmp), bs_args);
    bg   = bg - tmp;
    clear tmp;
end
g{n1+1}  = log(max(bg,eps));
clear bg

vx = sqrt(sum(NG.mat(1:3,1:3).^2));
%=======================================================================

%=======================================================================
function [y,dt] = defdet(u,prm,int_args)
% Generate deformation
[y,J] = spm_shoot3d(u,prm,int_args);
dt    = spm_diffeo('det',J);
%=======================================================================

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


text:
function peakVals = ode_fun ( bVals, kVals )

%*****************************************************************************80
%
%% ODE_FUN finds peak values of an ODE over a mesh of parameter values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2011
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, real bVals(M), a vector of B values.
%
%    Input, real kVals(N), a vector of K values.
%
%    Output, real peakVals(M,N), a vector containing the maximum value of
%    the ODE solution over a fixed time interval, for every combination
%    of B and K values.
%

%
%  Form a grid of all pairs of K and B:
%
  [ kGrid, bGrid ] = meshgrid ( bVals, kVals );
%
%  Define an array to hold the results, and initialize it to NAN.
%
  peakVals = nan ( size ( kGrid ) );
%
%  Solve the ODE for every pair of K and B values.  (M is fixed at 5.)
%
  m = 5.0;

  for ij = 1 : numel(kGrid)
%
%  Solve the ODE over the time interval 0 <= T <= 25, with
%  initial conditions X(0) = 0, X'(0) = 1.
%
    [ T, Y ] = ode45 ( @(t,y) ode_system ( t, y, m, bGrid(ij), kGrid(ij) ), ...
      [0, 25],  [0, 1] );
%
%  Retrieve the maximum value achieved by this solution.
%
    peakVals(ij) = max ( Y(:,1) );

  end

  return
end


text:
% Calculate the number of connected components using the eigenvalues
%                   of the Laplacian - counting the number of zeros
%
% INPUTS: adjacency matrix, nxn
% OUTPUTs: positive integer - number of connected components
%
% Other routines used: graphSpectrum.m
% GB: last updated: September 22, 2012

function nc=numConnComp(adj)

s=graphSpectrum(adj);
nc=numel(find(s<10^(-5)));   % zero eigenvalues are sometimes close to zeros numerically

author: aeolianine

repo: octave-networks-toolbox


text:
function [pzout, intv00f, intv00g, intv00h] =  gskpzget(pzinl,intv0ln);
pzout=[];
pzout(1,1:3) = pzinl(intv0ln,1:3);
pzout(2,1:3) = pzinl(intv0ln,4:6);
pzout(3,1:3) = pzinl(intv0ln,7:9);
intv00f = pzinl(intv0ln,10);
intv00g = pzinl(intv0ln,11);
intv00h = pzinl(intv0ln,12);


text:
function [ vx, vy ] = basisd_serene ( xq, yq, xw, ys, xe, yn, xx, yy )

%*****************************************************************************80
%
%% BASISD_SERENE differentiates the serendipity basis functions.
%
%  Discussion:
%
%    This procedure assumes that a serendipity element has been defined,
%    whose sides are parallel to coordinate axes.
%
%    The local element numbering is
%
%      YN  3--2--1
%       |  |     |
%       |  4     8
%       |  |     |
%      YS  5--6--7
%       |
%       +--XW---XE--
%
%    We note that each basis function can be written as the product of
%    three linear terms, which never result in an x^2y^2 term.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XQ, YQ, the evaluation point.
%
%    Input, real XW, YS, the coordinates of the lower left corner.
%
%    Input, real XE, YN, the coordinates of the upper right corner.
%
%    Input, real XX(8), YY(8), the coordinates of the 8 nodes.
%
%    Output, real VX(8), VY(8), the derivatives of the basis functions
%    at (XQ,YQ) with respect to X and Y.
%
  vx = zeros(8,1);
  vy = zeros(8,1);

  vx(1) = ...
      not1d ( xq, xw, xx(1) ) ...
    * not1 ( yq, ys, yy(1) ) ...
    * not2 ( xq, yq, xx(8), yy(8), xx(2), yy(2), xx(1), yy(1) ) ...
    + not1 ( xq, xw, xx(1) ) ...
    * not1 ( yq, ys, yy(1) ) ...
    * not2dx ( xq, yq, xx(8), yy(8), xx(2), yy(2), xx(1), yy(1) );

  vy(1) = ...
      not1 ( xq, xw, xx(1) ) ...
    * not1d ( yq, ys, yy(1) ) ...
    * not2 ( xq, yq, xx(8), yy(8), xx(2), yy(2), xx(1), yy(1) ) ...
    + not1 ( xq, xw, xx(1) ) ...
    * not1 ( yq, ys, yy(1) ) ...
    * not2dy ( xq, yq, xx(8), yy(8), xx(2), yy(2), xx(1), yy(1) );

  vx(2) = ...
      not1d ( xq, xw, xx(2) ) ...
    * not1 ( xq, xe, xx(2) ) ...
    * not1 ( yq, ys, yy(2) ) ...
    + not1 ( xq, xw, xx(2) ) ...
    * not1d ( xq, xe, xx(2) ) ...
    * not1 ( yq, ys, yy(2) );

  vy(2) = ...
      not1 ( xq, xw, xx(2) ) ...
    * not1 ( xq, xe, xx(2) ) ...
    * not1d ( yq, ys, yy(2) );

  vx(3) = ...
      not1d ( xq, xe, xx(3) ) ...
    * not1 ( yq, ys, yy(3) ) ...
    * not2 ( xq, yq, xx(2), yy(2), xx(4), yy(4), xx(3), yy(3) ) ...
    + not1 ( xq, xe, xx(3) ) ...
    * not1 ( yq, ys, yy(3) ) ...
    * not2dx ( xq, yq, xx(2), yy(2), xx(4), yy(4), xx(3), yy(3) );

  vy(3) = not1 ( xq, xe, xx(3) ) ...
    * not1d ( yq, ys, yy(3) ) ...
    * not2 ( xq, yq, xx(2), yy(2), xx(4), yy(4), xx(3), yy(3) ) ...
    + not1 ( xq, xe, xx(3) ) ...
    * not1 ( yq, ys, yy(3) ) ...
    * not2dy ( xq, yq, xx(2), yy(2), xx(4), yy(4), xx(3), yy(3) );

  vx(4) = ...
      not1d ( xq, xe, xx(4) ) ...
    * not1 ( yq, yn, yy(4) ) ...
    * not1 ( yq, ys, yy(4) );

  vy(4) = ...
      not1 ( xq, xe, xx(4) ) ...
    * not1d ( yq, yn, yy(4) ) ...
    * not1 ( yq, ys, yy(4) ) ...
    + not1 ( xq, xe, xx(4) ) ...
    * not1 ( yq, yn, yy(4) ) ...
    * not1d ( yq, ys, yy(4) );

  vx(5) = ...
      not1d ( xq, xe, xx(5) ) ...
    * not1 ( yq, yn, yy(5) ) ...
    * not2 ( xq, yq, xx(4), yy(4), xx(6), yy(6), xx(5), yy(5) ) ...
    + not1 ( xq, xe, xx(5) ) ...
    * not1 ( yq, yn, yy(5) ) ...
    * not2dx ( xq, yq, xx(4), yy(4), xx(6), yy(6), xx(5), yy(5) );

  vy(5) = ...
      not1 ( xq, xe, xx(5) ) ...
    * not1d ( yq, yn, yy(5) ) ...
    * not2 ( xq, yq, xx(4), yy(4), xx(6), yy(6), xx(5), yy(5) ) ...
    + not1 ( xq, xe, xx(5) ) ...
    * not1 ( yq, yn, yy(5) ) ...
    * not2dy ( xq, yq, xx(4), yy(4), xx(6), yy(6), xx(5), yy(5) );

  vx(6) = ...
      not1d ( xq, xe, xx(6) ) ...
    * not1 ( xq, xw, xx(6) ) ...
    * not1 ( yq, yn, yy(6) ) ...
    + not1 ( xq, xe, xx(6) ) ...
    * not1d ( xq, xw, xx(6) ) ...
    * not1 ( yq, yn, yy(6) );

  vy(6) = ...
      not1 ( xq, xe, xx(6) ) ...
    * not1 ( xq, xw, xx(6) ) ...
    * not1d ( yq, yn, yy(6) );

  vx(7) = ...
      not1d ( xq, xw, xx(7) ) ...
    * not1 ( yq, yn, yy(7) ) ...
    * not2 ( xq, yq, xx(6), yy(6), xx(8), yy(8), xx(7), yy(7) ) ...
    + not1 ( xq, xw, xx(7) ) ...
    * not1 ( yq, yn, yy(7) ) ...
    * not2dx ( xq, yq, xx(6), yy(6), xx(8), yy(8), xx(7), yy(7) );

  vy(7) = ...
      not1 ( xq, xw, xx(7) ) ...
    * not1d ( yq, yn, yy(7) ) ...
    * not2 ( xq, yq, xx(6), yy(6), xx(8), yy(8), xx(7), yy(7) ) ...
    + not1 ( xq, xw, xx(7) ) ...
    * not1 ( yq, yn, yy(7) ) ...
    * not2dy ( xq, yq, xx(6), yy(6), xx(8), yy(8), xx(7), yy(7) );

  vx(8) = ...
      not1 ( yq, ys, yy(8) ) ...
    * not1 ( yq, yn, yy(8) ) ...
    * not1d ( xq, xw, xx(8) );

  vy(8) = ...
      not1d ( yq, ys, yy(8) ) ...
    * not1 ( yq, yn, yy(8) ) ...
    * not1 ( xq, xw, xx(8) ) ...
    + not1 ( yq, ys, yy(8) ) ...
    * not1d ( yq, yn, yy(8) ) ...
    * not1 ( xq, xw, xx(8) );

  return
end


text:
function tests = VehicleTest
  tests = functiontests(localfunctions);
  clc
end

function Bicycle_constructor_test(tc)

    % default constructor
    v = Bicycle();

    % display
    v

    % char
    s = v.char()
    tc.verifyTrue( ischar(s) );

    % all options
    v = Bicycle( ...
        'steermax', 2, ...
        'accelmax', 3, ...
        'covar', [1 2;3 4], ...
        'speedmax', 5, ...
        'L',2.5, ...
        'x0', [1 2 3], ...
        'dt', 0.5, ...
        'rdim', 0.3 ...
        );

    tc.verifyEqual(v.steermax, 2);
    tc.verifyEqual(v.accelmax, 3);
    tc.verifyEqual(v.V, [1 2;3 4]);
    tc.verifyEqual(v.speedmax, 5);
    tc.verifyEqual(v.L, 2.5);
    tc.verifyEqual(v.x0, [1 2 3]');
    tc.verifyEqual(v.dt, 0.5);
    tc.verifyEqual(v.rdim, 0.3);

    v.init()
    tc.verifyEqual(v.x, [1 2 3]');

end

function Bicycle_deriv_test(tc)

    v = Bicycle('steermax', Inf);

    xd =v.deriv([], [0 0 0], [1 0]);
    tc.verifyEqual( xd, [1 0 0]', 'AbsTol', 1e-6);

    xd =v.deriv([], [0 0 pi/2], [1 0]);
    tc.verifyEqual( xd, [0 1 0]', 'AbsTol', 1e-6);

    xd =v.deriv([], [0 0 0], [0 1]);
    tc.verifyEqual( xd, [0 0 0]', 'AbsTol', 1e-6);

    xd =v.deriv([], [0 0 0], [1 pi/4]);
    tc.verifyEqual( xd, [1 0 1]', 'AbsTol', 1e-6);

    v = Bicycle('steermax', pi/4);
    xd =v.deriv([], [0 0 0], [1 100]);
    tc.verifyEqual( xd, [1 0 1]', 'AbsTol', 1e-6);

    v = Bicycle('speedmax', 1, 'steermax', Inf);
    xd =v.deriv([], [0 0 0], [100 pi/4]);
    tc.verifyEqual( xd, [1 0 1]', 'AbsTol', 1e-6);

    v = Bicycle('accelmax', 1);
    xd =v.deriv([], [0 0 0], [100 0]);
    tc.verifyEqual( xd, [v.dt 0 0]', 'AbsTol', 1e-6);
    for i=1:9
            xd =v.deriv([], [0 0 0], [100 0]);
    end
    tc.verifyEqual( xd, [1 0 0]', 'AbsTol', 1e-6);

end

function Bicycle_update_test(tc)

    v = Bicycle('dt', 1, 'speedmax', Inf, 'steermax', Inf);

    tc.verifyEqual( v.update([1 0]), [1 0], 'AbsTol', 1e-6);
    tc.verifyEqual( v.x, [1 0 0]', 'AbsTol', 1e-6);

    tc.verifyEqual( v.update([1 0]), [1 0], 'AbsTol', 1e-6);
    tc.verifyEqual( v.x, [2 0 0]', 'AbsTol', 1e-6);

    tc.verifyEqual( v.update([-2 0]), [2 0], 'AbsTol', 1e-6);
    tc.verifyEqual( v.x, [0 0 0]', 'AbsTol', 1e-6);

    tc.verifyEqual( size(v.x_hist,1), 3);

    v.init();
    tc.verifyEqual( size(v.x_hist,1), 0);

    tc.verifyEqual( v.update([0 1]), [0 0], 'AbsTol', 1e-6);
    tc.verifyEqual( v.update([1 pi/4]), [1 1], 'AbsTol', 1e-6);
end

function Bicycle_f_test(tc)

    v = Bicycle('steermax', Inf);

    tc.verifyEqual( v.f([0 0 0], [1 0], [0 0]), [1 0 0], 'AbsTol', 1e-6);
    tc.verifyEqual( v.f([2 3 0], [1 0], [0 0]), [3 3 0], 'AbsTol', 1e-6);

    tc.verifyEqual( v.f([0 0 0], [0 0], [1 0]), [1 0 0], 'AbsTol', 1e-6);

    tc.verifyEqual( v.f([0 0 0], [0 1], [0 0]), [0 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.f([0 0 0], [sqrt(2) pi/4], [0 0]), [1 1 pi/4], 'AbsTol', 1e-6);
    tc.verifyEqual( v.f([0 0 0], [sqrt(2) 0], [0 pi/4]), [1 1 pi/4], 'AbsTol', 1e-6);
end

function Bicycle_jacobian_test(tc)

    v = Bicycle('steermax', Inf);

    tc.verifyEqual( v.Fx([0 0 0], [0 0]), [1 0 0; 0 1 0; 0 0 1], 'AbsTol', 1e-6);

    tc.verifyEqual( v.Fx([0 0 0], [2 0]), [1 0 0; 0 1 2; 0 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.Fx([2 3 0], [2 0]), [1 0 0; 0 1 2; 0 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.Fx([0 0 0], [3 pi/2]), [1 0 -3; 0 1 0; 0 0 1], 'AbsTol', 1e-6);

    tc.verifyEqual( v.Fv([0 0 0], [0 0]), [1 0; 0 0; 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.Fv([2 3 0], [0 0]), [1 0; 0 0; 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.Fv([0 0 0], [2 3]), [1 0; 0 0; 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.Fv([0 0 pi/2], [0 0]), [0 0; 1 0; 0 1], 'AbsTol', 1e-6);
end

function Unicycle_constructor_test(tc)

    % default constructor
    v = Unicycle();

    % display
    v

    % char
    s = v.char()
    tc.verifyTrue( ischar(s) );

    % all options
    v = Bicycle( ...
        'accelmax', 3, ...
        'covar', [1 2;3 4], ...
        'speedmax', 5, ...
        'L',2.5, ...
        'x0', [1 2 3], ...
        'dt', 0.5, ...
        'rdim', 0.3 ...
        );

    tc.verifyEqual(v.accelmax, 3);
    tc.verifyEqual(v.V, [1 2;3 4]);
    tc.verifyEqual(v.speedmax, 5);
    tc.verifyEqual(v.L, 2.5);
    tc.verifyEqual(v.x0, [1 2 3]');
    tc.verifyEqual(v.dt, 0.5);
    tc.verifyEqual(v.rdim, 0.3);

    v.init()
    tc.verifyEqual(v.x, [1 2 3]');
end

function Unicycle_deriv_test(tc)

    v = Unicycle();

    xd =v.deriv([], [0 0 0], [1 0]);
    tc.verifyEqual( xd, [1 0 0]', 'AbsTol', 1e-6);

    xd =v.deriv([], [0 0 pi/2], [1 0]);
    tc.verifyEqual( xd, [0 1 0]', 'AbsTol', 1e-6);

    xd =v.deriv([], [0 0 0], [0 1]);
    tc.verifyEqual( xd, [0 0 1]', 'AbsTol', 1e-6);

    xd =v.deriv([], [0 0 0], [1 1]);
    tc.verifyEqual( xd, [1 0 1]', 'AbsTol', 1e-6);

    v = Unicycle('speedmax', 1);
    xd =v.deriv([], [0 0 0], [100 0]);
    tc.verifyEqual( xd, [1 0 0]', 'AbsTol', 1e-6);

    v = Unicycle('accelmax', 1);
    xd =v.deriv([], [0 0 0], [100 0]);
    tc.verifyEqual( xd, [v.dt 0 0]', 'AbsTol', 1e-6);
    for i=1:9
            xd =v.deriv([], [0 0 0], [100 0]);
    end
    tc.verifyEqual( xd, [1 0 0]', 'AbsTol', 1e-6);

end

function Unicycle_update_test(tc)

    v = Unicycle('dt', 1, 'speedmax', Inf);

    tc.verifyEqual( v.update([1 0]), [1 0], 'AbsTol', 1e-6);
    tc.verifyEqual( v.x, [1 0 0]', 'AbsTol', 1e-6);

    tc.verifyEqual( v.update([1 0]), [1 0], 'AbsTol', 1e-6);
    tc.verifyEqual( v.x, [2 0 0]', 'AbsTol', 1e-6);

    tc.verifyEqual( v.update([-2 0]), [2 0], 'AbsTol', 1e-6);
    tc.verifyEqual( v.x, [0 0 0]', 'AbsTol', 1e-6);

    tc.verifyEqual( size(v.x_hist,1), 3);

    v.init();
    tc.verifyEqual( size(v.x_hist,1), 0);

    tc.verifyEqual( v.update([0 1]), [0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.update([2 1]), [2 1], 'AbsTol', 1e-6);
end

function Unicycle_f_test(tc)

    v = Unicycle();

    tc.verifyEqual( v.f([0 0 0], [1 0], [0 0]), [1 0 0], 'AbsTol', 1e-6);
    tc.verifyEqual( v.f([2 3 0], [1 0], [0 0]), [3 3 0], 'AbsTol', 1e-6);

    tc.verifyEqual( v.f([0 0 0], [0 0], [1 0]), [1 0 0], 'AbsTol', 1e-6);

    tc.verifyEqual( v.f([0 0 0], [0 1], [0 0]), [0 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.f([0 0 0], [sqrt(2) pi/4], [0 0]), [1 1 pi/4], 'AbsTol', 1e-6);
    tc.verifyEqual( v.f([0 0 0], [sqrt(2) 0], [0 pi/4]), [1 1 pi/4], 'AbsTol', 1e-6);
end

function Unicycle_jacobian_test(tc)

    v = Unicycle();

    tc.verifyEqual( v.Fx([0 0 0], [0 0]), [1 0 0; 0 1 0; 0 0 1], 'AbsTol', 1e-6);

    tc.verifyEqual( v.Fx([0 0 0], [2 0]), [1 0 0; 0 1 2; 0 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.Fx([2 3 0], [2 0]), [1 0 0; 0 1 2; 0 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.Fx([0 0 0], [3 pi/2]), [1 0 -3; 0 1 0; 0 0 1], 'AbsTol', 1e-6);

    tc.verifyEqual( v.Fv([0 0 0], [0 0]), [1 0; 0 0; 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.Fv([2 3 0], [0 0]), [1 0; 0 0; 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.Fv([0 0 0], [2 3]), [1 0; 0 0; 0 1], 'AbsTol', 1e-6);
    tc.verifyEqual( v.Fv([0 0 pi/2], [0 0]), [0 0; 1 0; 0 1], 'AbsTol', 1e-6);
end

repo: robotics-toolbox-matlab


text:
function varargout = max(varargin)
% global, local and pointwise maxima of spherical functions
%
% Syntax
%   [v,pos] = max(sF) % the position where the maximum is atained
%
%   [v,pos] = max(sF,'numLocal',5) % the 5 largest local maxima
%
%   sF = max(sF, c) % maximum of a spherical functions and a constant
%   sF = max(sF1, sF2) % maximum of two spherical functions
%   sF = max(sF1, sF2, 'bandwidth', bw) % specify the new bandwidth
%
%   % compute the maximum of a multivariate function along dim
%   sF = max(sFmulti,[],dim)
%
% Input
%  sF, sF1, sF2 - @S2Fun
%  sFmulti - a multivariate @S2Fun
%  c       - double
%
% Output
%  v - double
%  pos - @vector3d
%
% Options
%  kmax - number of iterations
%  numLocal      - number of peaks to return
%  startingNodes - @vector3d
%  tolerance     - minimum distance between two peaks
%  resolution    - minimum step size
%  maxStepSize   - maximm step size
%

varargin{1} = -varargin{1};
if nargin>1 && ( ~isempty(varargin{2}) && isa(varargin{2},'double') || isa(varargin{2},'S2Fun' ))
  varargin{2} = -varargin{2};
end

[varargout{1:nargout}] = min(varargin{:});

varargout{1} = - varargout{1};

end

author: mtex-toolbox

repo: mtex

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


text:
function Q=cotes(f,a,b,N,nodes)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% cotes.m
%
% 2 to 11 point Summed Newton-Cotes numerical integration formulas.
%
% Inputs:
%
% f - function to be integrated
% a - lower limit of integration
% b - upper limit of integration
% N - Requested number of grid points
% nodes - number of nodes in Newton-Cotes formula
%
% Sample Usage:
%
% >>cotes(@sin,0,pi/2,20,5)
%
% ans =
%         0.999999999501637
%
% Written by: Greg von Winckel
% Contact: gregvw(at)math(dot)unm(dot)edu
% URL: http://math.unm.edu/~gregvw
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = (nodes-1)*ceil(N/(nodes-1));  N1=N+1;
x = linspace(a,b,N1)';
h = x(2)-x(1);    g=f(x);

endpts = g(1)+g(N1);

switch nodes

    case{2} % Trapezoidal Rule
    Q=(h/2)*(endpts+2*sum(g(2:N)));

    case{3} % Simpson's Rule
    Q=(h/3)*(endpts+4*sum(g(2:2:N))+2*sum(g(3:2:N)));

    case{4} % Simpson's 3/8 Rule
    Q=(3*h/8)*(endpts+3*sum(g(2:3:N)+g(3:3:N))+2*sum(g(4:3:N)));

    case{5} % Boole's 4/90 Rule
    Q=(4*h/90)*(7*endpts+32*sum(g(2:4:N))+...
        12*sum(g(3:4:N))+32*sum(g(4:4:N))+14*sum(g(5:4:N)));

    case{6}
    Q=(5*h/288)*(19*endpts+75*sum(g(2:5:N)+g(5:5:N))+...
        50*sum(g(3:5:N)+g(4:5:N))+38*sum(g(6:5:N)));

    case{7}
    Q=(6*h/840)*(41*endpts+216*sum(g(2:6:N)+g(6:6:N))+...
        27*sum(g(3:6:N)+g(5:6:N))+272*sum(g(4:6:N))+...
        +82*sum(g(7:6:N)));

    case{8}
    Q=(7*h/17280)*(751*endpts+3577*sum(g(2:7:N)+g(7:7:N))+...
        1323*sum(g(3:7:N)+g(6:7:N))+2989*sum(g(4:7:N)+g(5:7:N))+...
        +1502*sum(g(8:7:N)));

    case{9}
    Q=(4*h/14175)*(989*endpts+5888*sum(g(2:8:N)+g(8:8:N))-...
       928*sum(g(3:8:N)+g(7:8:N))+10496*sum(g(4:8:N)+g(6:8:N))-...
       4540*sum(g(5:8:N))+1978*sum(g(9:8:N)));

    case{10}
    Q=(9*h/89600)*(2857*endpts+15741*sum(g(2:9:N)+g(9:9:N))+...
       1080*sum(g(3:9:N)+g(8:9:N))+19344*sum(g(4:9:N)+g(7:9:N))+...
       5778*sum(g(5:9:N)+g(6:9:N))+5714*sum(g(10:9:N)));

    case{11}
    Q=(5*h/299376)*(16067*endpts+106300*sum(g(2:10:N)+g(10:10:N))-...
       48525*sum(g(3:10:N)+g(9:10:N))+272400*sum(g(4:10:N)+g(8:10:N))-...
       260550*sum(g(5:10:N)+g(7:10:N))+427368*sum(g(6:10:N))+...
       32134*sum(g(11:10:N)));
    otherwise
        error('Order must be between 2 and 11');
end

repo: Aero-matlab


text:
function [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse,mets)
% contractModel
%   Contracts a model by grouping all identical reactions. Similar to the
%   deleteDuplicates part in simplifyModel but more care is taken here
%   when it comes to gene associations. If the duplicated reactions have
%   '_EXP_*' suffixes (where * is a digit), then the model is assumed to
%   have been passed through expandModel, and these suffixes are removed
%   here.
%
%   model                  a model structure
%   distReverse            distinguish reactions with same metabolites
%                          but different reversibility as different
%                          reactions (opt, default true)
%   mets                   string or cell array of strings with metabolite
%                          identifiers, whose involved reactions should be
%                          checked for duplication (opt, by default all
%                          reactions are considered) (option is used by
%                          replaceMets)
%
%   reducedModel           a model structure without duplicate reactions
%   removedRxns            cell array for the removed duplicate reactions
%   indexedDuplicateRxns   indexed cell array for the removed duplicate
%                          reactions (multiple valuess separated by semicolon)
%
%   NOTE: This code might not work for advanced grRules strings
%         that involve nested expressions of 'and' and 'or'.
%
%   Usage: [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse)

if nargin<2
    distReverse=true;
end

%First sort the model so that reversible reactions are in the same
%direction
modelS=sortModel(model);

%Get a list of duplicate reactions
if distReverse
    x=[modelS.S; model.rev']';
else
    x=modelS.S';
end
if nargin>2
    toCheck = getIndexes(model,mets,'mets');
    [~, toCheck] = find(modelS.S(toCheck,:));
    x = [x, (1:size(x,1))']; % Make all other rxns unique by additional column
    x(toCheck,end) = 0;
end

[~,I,J] = unique(x,'rows','first');

%Initialize cell array of indexedDuplicateRxns
indexedDuplicateRxns=cell(numel(model.rxns),1);
indexedDuplicateRxns(:)={''};

duplicateRxns=setdiff(1:numel(model.rxns),I);
mergeTo=I(J(duplicateRxns));

mergedRxns=unique(mergeTo);

%Now add all the info from this one. Print a warning if they have different
%bounds or objective function coefficients. Uses the widest bounds and
%largest magnitude of objective coefficient
for i=1:numel(mergedRxns)
    duplRxn=transpose([mergedRxns(i),duplicateRxns(mergeTo==mergedRxns(i))]);
    if numel(unique(model.lb(duplRxn)))>1
        EM=['Duplicates of reaction ' model.rxns{mergedRxns(i)} ' have different lower bound. Uses the most negative/smallest lower bound'];
        dispEM(EM,false);
        model.lb(mergedRxns(i))=min(model.lb(duplRxn));
    end
    if numel(unique(model.ub(duplRxn)))>1
        EM=['Duplicates of reaction ' model.rxns{mergedRxns(i)} ' have different upper bound. Uses the most positive/largest upper bound'];
        dispEM(EM,false);
        model.ub(mergedRxns(i))=max(model.ub(duplRxn));
    end
    if numel(unique(model.c(duplRxn)))>1
        EM=['Duplicates of reaction ' model.rxns{mergedRxns(i)} ' has a different objective function coefficient. Uses the largest coefficient'];
        dispEM(EM,false);
        model.c(mergedRxns(i))=max(model.c(duplRxn));
    end
    if isfield(model,'grRules') && any(~isempty(model.grRules(duplRxn)))
        rules=model.grRules(duplRxn);
        allRules={};
        for j=1:numel(rules)
            rules{j}=ignoreORinComplex(rules{j});
            allRules=[allRules regexp(rules{j},' or ','split')];
        end
        allRules=unique(allRules);
        allRules=strrep(allRules,'__OR__',' or ');
        andRules=contains(allRules,' and ');
        allRules(andRules)=strcat('(',allRules(andRules),')');
        if numel(allRules)==1
            model.grRules{mergedRxns(i)}=allRules{1};
        else
            model.grRules{mergedRxns(i)}=strjoin(allRules,' or ');
        end
    end
    if isfield(model,'eccodes') && any(~isempty(model.eccodes(duplRxn)))
        codes=model.eccodes(duplRxn);
        allCodes={};
        for j=1:numel(codes)
            allCodes=[allCodes regexp(codes{j},';','split')];
        end
        allCodes=unique(allCodes);
        if numel(allCodes)==1
            model.eccodes{mergedRxns(i)}=allCodes{1};
        else
            model.eccodes{mergedRxns(i)}=strjoin(allCodes,';');
        end
    end
    %Generate indexedDuplicateRxns cell array
    if numel(duplRxn)==2
        indexedDuplicateRxns{duplRxn(1)}=model.rxns{duplRxn(2)};
    else
        indexedDuplicateRxns{duplRxn(1)}=strjoin(model.rxns(duplRxn(2:end)),';');
    end
    %If all reactions have _EXP_* (where * is digit) as suffix, then the
    %model had passed through expandModel(), and these suffixes are removed.
    if all(~cellfun(@isempty,regexp(model.rxns(duplRxn),'_EXP_\d+$')))
        model.rxns(duplRxn(1))=regexprep(model.rxns(duplRxn(1)),'_EXP_\d+$','');
    end
end

%Delete the duplicate reactions
reducedModel=removeReactions(model,duplicateRxns);
removedRxns=model.rxns(duplicateRxns);
[~, index]=ismember(reducedModel.rxns,model.rxns);
indexedDuplicateRxns=indexedDuplicateRxns(index);

if isfield(reducedModel,'rxnGeneMat')
    %Fix grRules and reconstruct rxnGeneMat
    [grRules,rxnGeneMat] = standardizeGrRules(reducedModel,true);
    reducedModel.grRules = grRules;
    reducedModel.rxnGeneMat = rxnGeneMat;
end
end

function grRule = ignoreORinComplex(grRule)
%In a grRule, if OR relationship is nested in an AND relationship, then
%obfuscate the OR before splitting isoenzymes
grRule=['(' grRule ')'];
brOpen=strfind(grRule,'(');
brClose=strfind(grRule,')');
andPos=strfind(grRule,' and ');
%Find opening bracket closest before AND
stillCapturing = 0;
for i=1:numel(andPos)
    searchPos = andPos(i);
    while stillCapturing == 0
        closestOpen = brOpen(max(find(brOpen<searchPos)));
        inbetweenClose = brClose(brClose<searchPos & brClose>closestOpen);
        if ~isempty(inbetweenClose)
            searchPos=max(inbetweenClose);
        else
            stillCapturing = 1;
            beginPos = closestOpen;
        end
    end
    stillCapturing = 0;
    searchPos = andPos(i);
    while stillCapturing == 0
        closestClose = brClose(min(find(brClose>searchPos)));
        inbetweenOpen = brOpen(brOpen>searchPos & brOpen<closestOpen);
        if ~isempty(inbetweenOpen)
            searchPos=min(closestClose);
        else
            stillCapturing = 1;
            endPos = closestClose;
        end
    end
    replacePart=regexprep(grRule(beginPos:endPos),' or ','__OR__');
    grRule=strrep(grRule,grRule(beginPos:endPos),replacePart);
end
grRule=grRule(2:end-1);
end

author: SysBioChalmers

repo: RAVEN


text:
function he_plot ( a, b, index, filename )

%*****************************************************************************80
%
%% HE_PLOT plots HE(i,x).
%
%  Discussion:
%
%    He(i,x) represents the probabilist's Hermite polynomial.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the plotting range.
%
%    Input, integer INDEX(*), the orders of 1 or more Hermite functions
%    to be plotted together.
%
%    Input, string FILENAME, the name into which the graphics information is
%    to be stored.  Note that the PNG format will be used.
%
  m = 501;
  x = linspace ( a, b, m );
  x = x';
  index_num = length ( index );

  clf
  hold on
  for i = 1 : index_num
    n = index(i);
    y = he_polynomial_value ( m, n, x );
    plot ( x, y(:,n+1), 'LineWidth', 2 );
  end
  grid on
  xlabel ( '<--- X --->', 'Fontsize', 16 )
  ylabel ( '<--- He(i,x) --->', 'Fontsize', 16 )
  title ( 'Hermite polynomials He(n,x)', 'Fontsize', 24 )
  hold off
  print ( '-dpng', filename )

  return
end

repo: jburkardt-m