text:
function [ Ep ] = tapas_rdcm_empty_par(DCM)
% [ Ep ] = tapas_rdcm_empty_par(DCM)
%
% Creates an empty valid parameter structure for a given DCM
%
%   Input:
%   	DCM         - model structure
%
%   Output:
%   	Ep          - parameter structure
%

% ----------------------------------------------------------------------
%
% Authors: Stefan Fraessle (stefanf@biomed.ee.ethz.ch), Ekaterina I. Lomakina
%
% Copyright (C) 2016-2022 Translational Neuromodeling Unit
%                         Institute for Biomedical Engineering
%                         University of Zurich & ETH Zurich
%
% This file is part of the TAPAS rDCM Toolbox, which is released under the
% terms of the GNU General Public License (GPL), version 3.0 or later. You
% can redistribute and/or modify the code under the terms of the GPL. For
% further see COPYING or <http://www.gnu.org/licenses/>.
%
% Please note that this toolbox is in an early stage of development. Changes
% are likely to occur in future releases.
%
% ----------------------------------------------------------------------

% number of regions and inputs
[nr, nu] = size(DCM.c);

% empty parameter structure
Ep.A        = zeros(nr,nr);
Ep.B        = zeros(nr,nr,nu);
Ep.C        = zeros(nr,nu);
Ep.D        = zeros(nr,nr,0);
Ep.transit  = zeros(nr,1);
Ep.decay    = zeros(nr,1);
Ep.epsilon  = 0;

end

author: translationalneuromodeling

repo: tapas


text:
function [Xs,Ys] = ml_get_grid(X,dims,num_pts )
%ML_GET_GRID  Returns Xs and Ys meshgrid with limits obtained from X.
%
%   input -----------------------------------------------------------------
%
%       o X : (N x D), data
%
%       o dims : (1 x 2), dimensions to create grid from
%

D = size(X,2);

if D == 1
    max_d1  = max(X);
    min_d1  = min(X);
    Xs      = linspace(min_d1,max_d1,num_pts);
    Ys      = [];
else
    max_d1  = max(X(:,dims(1)));
    min_d1  = min(X(:,dims(1)));
    max_d2  = max(X(:,dims(2)));
    min_d2  = min(X(:,dims(2)));
    [Xs,Ys] = meshgrid(linspace(min_d1,max_d1,num_pts),linspace(min_d2,max_d2,num_pts));
end

end

author: epfl-lasa


text:
function [x, xdot, x2dot] = modifiednewmarkint(M, C, K, R, x0, xdot0, ...
    t, varargin)
%Newmark's Direct Integration Method
%--------------------------------------------------------------------------
% Code written by : - Siva Srinivas Kolukula                              |
%                     Senior Research Fellow                              |
%                     Structural Mechanics Laboratory                     |
%                     Indira Gandhi Center for Atomic Research            |
%                     INDIA
%                   - Hamed Nokhostin                                     |
%                     B.S. Student                                        |
%                     Mechanical Engineering Faculty                      |
%                     K.N.Toosi Uiversity of Technoloygy                  |
%                     I.R.Iran                                            |
% E-mail : allwayzitzme@gmail.com
%          h_nokhostin@yahoo.com                                          |
%-------------------------------------------------------------------------
% PURPOSE
%        ???
% SYNTAX
%        [x, xdot, x2dot] = newmarkint(M, C, K, R, x0, xdot0, t, varargin)
% INPUT
%        [M] :       System Mass              [n,n]
%        [C] :       System Damping           [n,n]
%        [K] :       System Stiffness         [n,n]
%        [R] :       Externally Applied Load  [n,nt]
%        [x0] :      Initial Position         [n,1]
%        [xdot0] :   Initial Velocity         [n,1]
%        [t] :       Time Vector              [1,nt]
%        [varargin]: Options
%
% OUTPUT
%       [x]:        Displacemente Response   [n,nt]
%       [xdot]:     Velocity                 [n,nt]
%       [x2dot]:    Acceleration             [n,nt]
%
%
%  nt = number of time steps
%  n = number of nodes
% The options include changing the value of the "gamma" and "beta"
% coefficient which appear in the formulation of the method. By default
% these values are set to gamma = 1/2 and beta = 1/4.
%
% EXAMPLE
% To change nemark's coefficients, say to gamma = 1/3 and beta = 1/5,
% the syntax is:
%       [u, udot, u2dot] = newmark_int(t,p,u0,udot0,m,k,xi, 1/3, 1/5)
%
%-------------------------------------------------------------------------
if nargin == 7
    disp('Using default values:');
    disp('    gamma = 1/2');
    disp('    beta  = 1/4');
    gamma = 1 / 2;
    beta = 1 / 4;
elseif nargin == 9
    gamma = varargin{1};
    beta = varargin{2};
    disp('Using user''s values:');
    disp(['    gamma = ', num2str(alpha)]);
    disp(['    beta  = ', num2str(delta)]);
else
    error('Incorrect number of imput arguments');
end

dt = t(2) - t(1);
nt = fix((t(length(t) )- t(1)) / dt);
n = length(M);

% Constants used in Newmark's integration
a1 = gamma / (beta * dt);
a2 = 1 / (beta * dt ^ 2);
a3 = 1 / (beta * dt);
a4 = gamma / beta;
a5 = 1/(2 * beta);
a6 = (gamma / (2 * beta) - 1) * dt;

x = zeros(n,nt);
xdot = zeros(n,nt);
x2dot = zeros(n,nt);

% Initial Conditions
x(:, 1) = x0;
xdot(:, 1) = xdot0;
%R0 = zeros(n,1);
x2dot(:,1) = M \ (R(:, 1) - C * xdot(:, 1) - K * x(:, 1)) ;

Kcap = K + a1 * C + a2 * M;

a = a3 * M + a4 * C;
b = a5 * M + a6 * C;

% Tme step starts
for i = 1 : nt
    delR = R(:, i) + a * xdot(:, i) + b * x2dot(:, i);
    delx = Kcap \ delR ;
    delxdot = a1 * delx - a4 * xdot(:, i) - a6 * x2dot(:, i);
    delx2dot = a2 * delx - a3 * xdot(:, i) - a5 * x2dot(:, i);
    x(:, i + 1) = x(:, i) + delx;
    xdot(:, i + 1) = xdot(:, i) + delxdot;
    x2dot(:, i + 1) = x2dot(:, i) + delx2dot;
end

author: Sable

repo: mcbench-benchmarks


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% TLD is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with TLD.  If not, see <http://www.gnu.org/licenses/>.

path = 'X:\PhD\Publications\2009-06_OLCV2009\Web\MultiFace\';
clf;
files = img_dir(path);
load([path 'gt.mat']);
for i = 1:length(files)
    img =  imread([path files(i).name]);
    for ii = 1:length(gt)
        img = bb_burn(img,gt{ii}.bb(:,1763+i),255);
    end
    imshow(img);
    drawnow;
end

repo: OpenTLD


text:
function f = poissoncontunuos(x,lambda)
% function f = poissoncontunuos(x,lambda)
f = lambda.^x*exp(-lambda)./gamma(x+1);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function quadrule_test40 ( )

%*****************************************************************************80
%
%% QUADRULE_TEST40 tests NCO_COMPUTE and SUM_SUB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2005
%
%  Author:
%
%    John Burkardt
%
  order_max = 9;

  nfunc = func_set ( 'COUNT', 'DUMMY' );

  a = 0.0;
  b = 1.0;

  nsub = 1;
  xlo = -1.0;
  xhi = +1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADRULE_TEST40\n' );
  fprintf ( 1, '  NCO_COMPUTE sets up an open Newton-Cotes rule;\n' );
  fprintf ( 1, '  SUM_SUB carries it out.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Integration interval is [%f, %f]\n', a, b );
  fprintf ( 1, '  Number of subintervals is %d\n', nsub );
  fprintf ( 1, '  Quadrature order will vary.\n' );
  fprintf ( 1, '  Integrand will vary.\n' );
  fprintf ( 1, '\n' );

  for ilo = 1 : 5 : nfunc

    ihi = min ( ilo + 4, nfunc );

    fprintf ( 1, '\n' );
    fprintf ( 1, '    ' );
    for i = ilo : ihi
      fprintf ( '%14s', fname(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );

    for n = 1 : order_max

      if ( n == 8 )
        continue
      end

      fprintf ( 1, '  %2d', n );

      for i = ilo : ihi

        func_set ( 'SET', i );

        [ x, w ] = nco_compute ( n );

        result(i) = sum_sub ( @func, a, b, nsub, n, xlo, xhi, x, w );

        fprintf ( 1, '  %12f', result(i) );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [any_scanevent_repaired, ons, dur, index] ...
    = tapas_physio_repair_scan_events_PHILIPS(ons, sqpar, verbose)
% repairs the temporal structure of the sequence of scan events
% (time-stamps when slice was acquired) by filling up holes of missing
% slice acquisition or volume acquisition starts
%
% USAGE
%   [any_scanevent_repaired, ons, dur, index] ...
%       = tapas_physio_repair_scan_events_PHILIPS(ons, sqpar, verbose)
%
% -------------------------------------------------------------------------
% INPUT:
%   ons         - ons_samples, as given by tapas_physio_read_physlog_PHILIPS_guess_scanevents.m
%   sqpar       - sequence timing parameters
%           .Nslices        - number of slices per volume in fMRI scan
%           .NslicesPerBeat - usually equals Nslices, unless you trigger
%                             with the heart beat
%           .TR             - repetition time in seconds
%           .Ndummies       - number of dummy volumes
%           .Nscans         - number of full volumes saved (volumes in nifti file,
%                             usually rows in your design matrix)
%
%   verbose   - create informative plots (1= yes, 0 = no)
%
% -------------------------------------------------------------------------
% OUTPUT:
%   any_scanevent_repaired - flag, if function changed ons (1=yes, 0=no)
%   ons         - ons_samples, corrected for missing slice and volume
%                 scan events
%
% Note: The temporal unit of the output ons & dur corresponds to the input,
% i.e. if acq_sliceevents is given in samples, the output unit will be 2 ms
% usually, NOT 1 second!
%
% -------------------------------------------------------------------------
% Lars Kasper, August 2011
% Copyright (C) 2013, Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

thrmin = 1.5; % gap between slices, if spacing >thrmin*min(slicegap)
thrmax = 0.5; % gap is large, if gap > thrmax*max(slicegap)

%% repair scan events: some scan event triggers were missed so we have to
% fill up the holes
ons.acq_slice = ons.acq_slice_all(:);
dur.acq_slice = diff(ons.acq_slice);

%% =======================================================================
%% Find and remove gaps of one missing slice trigger within one volume scan of Nslices
index.hugegaps = find(dur.acq_slice>thrmax*max(dur.acq_slice));

%double gaps, think about it later...
index.doublegaps = setdiff(find(dur.acq_slice > (thrmin+1)*min(dur.acq_slice)),index.hugegaps);
if ~isempty(index.doublegaps)
    disp('There are double gaps or worse in scan events. I will not correct that!');
end

%one slice trigger missing in volume of NslicesPerBeat slices
index.gaps                  = find(dur.acq_slice>thrmin*min(dur.acq_slice) ...
                                & dur.acq_slice<(thrmin+1)*min(dur.acq_slice));
ons.acq_slice_gaps          = floor((ons.acq_slice(index.gaps)+ons.acq_slice(index.gaps+1))/2);
ons.acq_slice_filled_gaps   = sort([ons.acq_slice; ons.acq_slice_gaps]);

%% =======================================================================
%% gaps at beginning or end of Nslices-volume block still possible
dur.acq_most_slices = diff(ons.acq_slice_filled_gaps);

%indices of ons.acq_slice which belong to start of new volume
%ons.acq_vol=[1; find(dur.acq_most_slices>thrmin*min(dur.acq_most_slices))+1];

ind_vol = 1:sqpar.Nslices:length(ons.acq_slice_filled_gaps);%(length(ons.acq_slice_filled_gaps)-sqpar.Nslices*(Ndummies+Nscans)+1):sqpar.Nslices:length(ons.acq_slice_filled_gaps);
ons.acq_vol = ind_vol;%ons.acq_slice_filled_gaps(ind_vol);

%here we have scan volume blocks with e.g. N-1 slices, scan event trigger
%missed, e.g. for 3 slices: |||   ||  |||
index.gaps_startend_vol = find(diff(ons.acq_vol)<sqpar.NslicesPerBeat);

%OK till here
%is it |||   ||  ||| or |||  ||   ||| (larger gap before or after block
%which is too short?)
gapbefore   = ons.acq_slice_filled_gaps(ons.acq_vol(index.gaps_startend_vol))-ons.acq_slice_filled_gaps(ons.acq_vol(index.gaps_startend_vol)-1);
gapafter    = ons.acq_slice_filled_gaps(ons.acq_vol(index.gaps_startend_vol+1))-ons.acq_slice_filled_gaps(ons.acq_vol(index.gaps_startend_vol+1)-1);

ons.acq_slice_all = ons.acq_slice_filled_gaps;
ons.acq_slice_larger_gaps = [];
for i = 1:length(gapbefore)
    if (gapbefore(i) > gapafter(i))
        targetons=(ons.acq_vol(index.gaps_startend_vol(i)));
        ons.acq_slice_larger_gaps(i)=2*ons.acq_slice_filled_gaps(targetons)-ons.acq_slice_filled_gaps(targetons+1);
    else
        targetons=(ons.acq_vol(index.gaps_startend_vol(i)+1)-1);
        ons.acq_slice_larger_gaps(i)=2*ons.acq_slice_filled_gaps(targetons)-ons.acq_slice_filled_gaps(targetons-1);
    end
    ons.acq_slice_all=[ons.acq_slice_all; ons.acq_slice_larger_gaps(i)];
end

ons.acq_slice_all = sort(ons.acq_slice_all);
% END of filling up start/end slice triggers in volume
% update ons.acq_vol
dur.acq_slice_all = diff(ons.acq_slice_all);

%ons.acq_vol_all=[1; find(dur.acq_slice_all>thrmin*min(dur.acq_slice_all))+1];
ons.acq_vol_all = ons.acq_vol;

%1st volume added manually
ons.acq_vol_all = [ons.acq_slice_all(ons.acq_vol_all)];

any_scanevent_repaired = ~(length(ons.acq_slice)==length(ons.acq_slice_all) ...
    && length(ons.acq_vol)==length(ons.acq_vol_all));

if verbose && any_scanevent_repaired
    figh{1} = plot_fix_missing_scan_events(ons, ...
        dur, index);
end

end %repair_scan_events

%%========================================================================
%% auxiliary plot functions follow;

%%=========================================================================
%% plot fixing procedure in scantrigger-time-difference view
function fh = plot_fix_missing_scan_events(ons, dur, index)

fh = tapas_physio_get_default_fig_params();
set(fh,'Name','Read-In: Time-difference view between events: Fix missing scan events');
ax(1) = subplot(3,1,1);
hold off;
plot(1:length(dur.acq_slice),dur.acq_slice);
hold on;
plot(index.hugegaps,dur.acq_slice(index.hugegaps),'bo')
plot(index.doublegaps,dur.acq_slice(index.doublegaps),'go')
plot(index.gaps,dur.acq_slice(index.gaps),'kx')

title('As extracted from the logfile');
xlabel('scan event + 1');
ylabel('\Delta t (ms) from last event')
warning('off', 'MATLAB:legend:IgnoringExtraEntries');
legend(...
    'dur.acq_slice (time diff between scan events', ...
    'index.hugegaps (dur.acq\_slice>thrmax*max(dur.acq\_slice))', ...
    'index.doublegaps (dur.acq\_slice > (thrmin+1)*min(dur.acq\_slice), but no hugegap)', ...
    'index.gaps (roughly size of 2*sliceTR)'...
    );

ax(2) = subplot(3,1,2);
hold off;
plot(1:length(dur.acq_most_slices),dur.acq_most_slices);
hold on;
plot(ons.acq_vol,dur.acq_most_slices(ons.acq_vol),'bx')

title('After fixing the slice-trigger gaps within a volume');
xlabel('scan event + 1');
ylabel('\Delta t (ms) from last event')
legend(...
    'dur.acq\_most\_slices', ...
    'ons.acq\_vol');

ax(3) = subplot(3,1,3);
hold off;
plot(1:length(dur.acq_slice_all),dur.acq_slice_all);
hold on;
plot(ons.acq_vol_all,dur.acq_slice_all(ons.acq_vol),'bx')

legend(...
    'dur.acq\_slice\_all', ...
    'ons.acq\_vol\_all');

title('After inserting the slice triggers at the beginning or end of a volume');
xlabel('scan event + 1');
ylabel('\Delta t (ms) from last event')

linkaxes(ax,'x');
end


text:
classdef MimClearDataset < MimGuiPlugin
    % MimClearDataset. Gui Plugin for displaying an "about box" dialog
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimClearDataset is a Gui Plugin for the MIM Toolkit. The gui will
    %     create a button to run this plugin. Running this plugin will cause a
    %     splash screen dialog to be displayed.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'No Dataset'
        SelectedText = 'No Dataset'
        ToolTip = ''
        Category = 'Developer tools'
        Visibility = 'Developer'
        Mode = 'Segment'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Icon = 'clear_overlay.png'
        Location = 101
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ClearDataset;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.DeveloperMode && gui_app.IsDatasetLoaded;
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% close all, clc
%
% USAGE
%  cc
%
% INPUTS
%
% OUTPUTS
%
% EXAMPLE
%
% See also C, CCC
%
% Piotr's Image&Video Toolbox      Version 1.5
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

close all; clc

author: zouchuhang

repo: LayoutNet


text:
function status = test_binaural_synthesis(modus)
%TEST_BINAURAL_SYNTHESIS tests the binaural synthesis functions
%
%   Usage: status = test_binaural_synthesis(modus)
%
%   Input parameters:
%       modus   - 0: numerical
%                 1: visual (not available)
%
%   Output parameters:
%       status  - true or false

%   TEST_BINAURAL_SYNTHESIS(modus) tests the ir_wfs function for different
%   loudspeaker arrays and source models.

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

status = false;

%% ===== Checking of input parameters ====================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);
if modus
    warning('%s: visual modus not available.',upper(mfilename));
end

%% ===== Settings ========================================================
conf = SFS_config;
conf.c = 343;
conf.fs = 44100;
conf.secondary_sources.x0 = [];
conf.secondary_sources.center = [0 0 0];
conf.secondary_sources.size = 3;
conf.ir.usehcomp = false;
conf.ir.useinterpolation = true;
conf.N = 4096;
conf.dimension = '2.5D';
conf.driving_functions = 'default';
conf.usetapwin = true;
conf.tapwinlen = 0.3;
conf.wfs.usehpre = true;
conf.wfs.hpretype = 'FIR';
conf.wfs.hpreflow = 50;
conf.wfs.hprefhigh = 1200;
conf.delayline.resampling = 'none';
conf.delayline.filter = 'integer';
conf.debug = 0;
conf.showprogress = false;

%% ===== Main ============================================================
% Check if HRTF data set is available, download otherwise
basepath = get_sfs_path();
hrtf_file = [basepath '/data/HRTFs/QU_KEMAR_anechoic_3m.sofa'];
if ~exist(hrtf_file,'file')
    disp('Download');
    url = ['https://github.com/sfstoolbox/data/blob/master/', ...
           'HRTFs/QU_KEMAR_anechoic_3m.sofa?raw=true'];
    download_file(url,hrtf_file);
end
% Load HRTF data set
hrtf = SOFAload(hrtf_file);

%% ===== WFS 2.5D ========================================================
% === Linear array ===
conf.secondary_sources.geometry = 'linear';
conf.secondary_sources.number = 20;
X = [0 -2 0];
phi = pi/2;
conf.xref = X;
% Plane wave
src = 'pw';
xs = [0.5 -1 0];
ir_wfs(X,phi,xs,src,hrtf,conf);
% Point source
src = 'ps';
xs = [0 1 0];
ir_wfs(X,phi,xs,src,hrtf,conf);
% Focused source
src = 'fs';
xs = [0 -1 0 0 -1 0];
ir_wfs(X,phi,xs,src,hrtf,conf);

% === Circular array ===
conf.secondary_sources.geometry = 'circle';
conf.secondary_sources.number = 64;
X = [0 0 0];
phi = pi/2;
conf.xref = X;
% Plane wave
src = 'pw';
xs = [0.5 -1 0];
ir_wfs(X,phi,xs,src,hrtf,conf);
% Point source
src = 'ps';
xs = [0.5 2 0];
ir_wfs(X,phi,xs,src,hrtf,conf);
% Focused source
src = 'fs';
xs = [0.5 0.5 0 -1 -1 0];
ir_wfs(X,phi,xs,src,hrtf,conf);

% === Box shaped array ===
conf.secondary_sources.geometry = 'box';
conf.secondary_sources.number = 80;
X = [0 0 0];
phi = pi/2;
conf.xref = X;
% Plane wave
src = 'pw';
xs = [0.5 -1 0];
ir_wfs(X,phi,xs,src,hrtf,conf);
% Point source
src = 'ps';
xs = [0.5 2 0];
ir_wfs(X,phi,xs,src,hrtf,conf);
% Focused source
src = 'fs';
xs = [0.5 0.5 0 -1 -1 0];
ir_wfs(X,phi,xs,src,hrtf,conf);

status = true;

author: sfstoolbox

repo: sfs-matlab


text:
classdef TestEstimateGlobalMotionRansac
    %TestEstimateGlobalMotionRansac

    methods (Static)
        function test_1
            im1 = cv.imread(fullfile(mexopencv.root(),'test','RubberWhale1.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            im2 = cv.imread(fullfile(mexopencv.root(),'test','RubberWhale2.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            pt1 = cv.goodFeaturesToTrack(im1, 'MaxCorners',200);
            pt2 = cv.goodFeaturesToTrack(im2, 'MaxCorners',200);
            models = {'Translation', 'TranslationAndScale', 'Rotation', ...
                'Rigid', 'Similarity', 'Affine'};
            for i=1:numel(models)
                [M,rmse,ninliers] = cv.estimateGlobalMotionRansac(pt1, pt2, ...
                    'MotionModel',models{i}, 'RansacParams',models{i});
                validateattributes(M, {'numeric'}, {'2d', 'size',[3 3]});
                validateattributes(rmse, {'numeric'}, {'scalar', 'real'});
                validateattributes(ninliers, {'numeric'}, {'scalar', 'integer'});
            end
        end

        function test_error_argnum
            try
                cv.estimateGlobalMotionRansac();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [At,b,c,K,prep,origcoeff] = pretransfo(At,b,c,K,pars)

% [At,b,c,K,prep] = pretransfo(At,b,c,K)
%
% PRETRANSFO  Checks data and then transforms into internal SeDuMi format.
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% Nearly complete rewrite
% Copyright (c) 2013 Michael C. Grant
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

% -------------------------------------------------------------------------
% Make sure that all fields exist in K, and verify that they are valid
% -------------------------------------------------------------------------

if ~isfield(K,'f') || isempty(K.f)
    K.f = 0;
elseif numel(K.f) ~= 1 || K.f ~= floor(K.f) || K.f < 0 || ~isreal(K.f)
    error('K.f should be nonnegative integer')
end
if ~isfield(K,'l') || isempty(K.l)
    K.l = 0;
elseif K.l ~= floor(K.l) || K.l < 0 || ~isreal(K.l)
    error('K.l should be nonnegative integer')
end
if ~isfield(K,'q') || ~nnz(K.q)
    K.q = zeros(1,0);
else
    K.q = K.q(:)';
    if any(K.q ~= floor(K.q)) || any(K.q<2) || ~isreal(K.q)
        error('K.q should contain only integers bigger than 1')
    end
end
if ~isfield(K,'r') || ~nnz(K.r)
    K.r = zeros(1,0);
else
    K.r = K.r(:)';
    if any(K.r ~= floor(K.r)) || any(K.r<3) || ~isreal(K.r)
        error('K.r should contain only integers bigger than 2')
    end
end
if ~isfield(K,'s') || ~nnz(K.s)
    K.s = zeros(1,0);
else
    K.s = K.s(:)';
    if any(K.s ~= floor(K.s)) || any(K.s<1) || ~isreal(K.s)
        error('K.s should contain only positive integers')
    end
end
% As an alternative to the 'scomplex' flag, I've added a 'z' parameter
% containing a list of Hermitian semidefinite cone sizes. For now, this
% is just translated to 'scomplex' for you. In the future, we may use
% this internally *instead* of scomplex or rsdpN.
if ~isfield(K,'z') || ~nnz(K.z)
    K.z = zeros(1,0);
else
    K.z = K.z(:)';
    if any(K.z ~= floor(K.z)) || any(K.z<1) || ~isreal(K.z)
        error('K.z should contain only positive integers')
    end
end

N_f    = K.f;
N_l    = K.l;
N_fl   = N_f + N_l;
L_q    = length(K.q);
N_q    = sum(K.q);
L_r    = length(K.r);
N_r    = sum(K.r);
N_qr   = N_q + N_r;
L_qr   = L_q + L_r;
L_s    = length(K.s);
L_z    = length(K.z);
L_sz   = L_s + L_z;
N_s    = sum((K.s).^2);
N_z    = sum((K.z).^2);
N_sz   = N_s + N_z;
L_qrsz = L_qr + L_sz;
N_flqr = N_fl + N_qr;
N      = N_flqr + N_sz;

if ~isfield(K,'ycomplex') || isempty(K.ycomplex)
    K.ycomplex = zeros(1,0);
else
    K.ycomplex = sort(K.ycomplex(:))';
    K.ycomplex(find(~diff(K.ycomplex))+1) = [];
    if any(K.ycomplex ~= floor(K.ycomplex)) || any(K.ycomplex<1)
        error('K.ycomplex should contain only positive integers');
    elseif any(K.ycomplex > numel(b))
        error('Elements of K.ycomplex are out of range');
    end
end
if ~isfield(K,'xcomplex') || isempty(K.xcomplex)
    K.xcomplex = zeros(1,0);
else
    K.xcomplex = sort(K.xcomplex(:))';
    K.xcomplex(find(~diff(K.xcomplex))+1) = [];
    if any(K.xcomplex ~= floor(K.xcomplex)) || any(K.xcomplex<1)
        error('K.xcomplex should contain only positive integers');
    elseif any(K.xcomplex>N_flqr)
        error('Elements of K.xcomplex are out of range');
    end
end
if ~isfield(K,'scomplex') || isempty(K.scomplex)
    K.scomplex = zeros(1,0);
else
    K.scomplex = sort(K.scomplex(:))';
    K.scomplex(find(~diff(K.scomplex))+1) = [];
    if any(K.scomplex~=floor(K.scomplex)) || any(K.scomplex<1)
        error('K.scomplex should contain only positive integers');
    elseif any(K.scomplex>L_s)
        error('Elements of K.xcomplex are out of range');
    end
end
if L_z
    K.s = [ K.s, K.z ];
    K.scomplex = [ K.scomplex, L_s + 1 : L_sz ];
    K.z = zeros(1,0);
    L_s = L_s + L_z;
    N_s = N_s + N_z;
    L_z = 0; %#ok
    N_z = 0; %#ok
end

% -------------------------------------------------------------------------
% Verify the size and validity of At, b, and c
% N = # variables
% -------------------------------------------------------------------------

% SeDuMi assumes that if At is not consistent with K but At' is, that the
% user supplied the transpose of the coefficient matrix. This introduces a
% rare ambiguity in the case where At happens to be square. Past versions
% of SeDuMi would not have allowed this, instead rejecting the case where
% m >= N; however, with complex variables, the situation is not quite as
% clear, and there may technically be cases where m >= N is acceptable.

if ndims(At) > 2 %#ok
    error('A must be a matrix');
elseif nnz(isnan(At)) || nnz(isinf(At))
    error('A contains NaN or Inf');
elseif size(At,1) == N
    % nothing
elseif size(At,2) == N
    At = At';
else
    error('(At,K) size mismatch');
end
if all(size(b)>1)
    error('Parameter b must be a vector');
elseif any(isnan(b)) || any(isinf(b))
    error('b contains NaN or Inf');
elseif length(b) ~= size(At,2)
    error('(At,b) size mismatch');
else
    b = b(:);
end
if all(size(c)>1)
    error('Parameter c must be a vector');
elseif any(isnan(c)) || any(isinf(c))
    error('c contains NaN or Inf');
elseif length(c) ~= N
    error('(c,K) size mismatch');
else
    c = c(:);
end

% -------------------------------------------------------------------------
% Save the standardized data for further use if needed
% -------------------------------------------------------------------------

if isfield(pars,'errors') && pars.errors==1
    origcoeff.At=At;
    origcoeff.c=c;
    origcoeff.b=b;
    origcoeff.K=K;
else
    origcoeff=[];
end

% -------------------------------------------------------------------------
% Flag diagonal SDP blocks for removal
% -------------------------------------------------------------------------

% There is some serious MATLAB trickery here (if I do say so myself) that
% merits explanation. "spattern" contains the indices of symbolically
% nonzero elements of the dual variable z = c - At * y. This is a single
% vector across all LMI constraints, so "sblk" tells us which indices
% belong to which constraints. The "rem" statement is what determines if
% a particular element is off-diagonal. If there is even one off-diagonal
% element in an SDP, then its corresponding element of "sdiag" is false.
%
% This new version of the analysis replaces the old preprocessSDP(), and
% seems inexpensive enough to apply to all LMIs regardless of size/count.
% The previous version was applied more sparingly.
%
% In theory one could do a more complex analysis of the block structure of
% an LMI, potentially breaking a larger into smaller ones. But this would
% likely be significantly more expensive.

if L_s && ( ~isfield(pars,'sdp') || pars.sdp )
    ssiz        = (K.s).^2;
    strt        = cumsum([1,ssiz(1:end-1)]);
    sblk        = zeros(1,N_s);
    sblk(strt)  = 1;
    sblk        = cumsum(sblk);
    spattern    = find(c(N_flqr+1:N)~=0|any(At(N_flqr+1:N,:),2))';
    sblk        = sblk(spattern);
    sblk        = sblk(rem(spattern-strt(sblk),K.s(sblk)+1)~=0);
    sdiag       = true(1,L_s);
    sdiag(sblk) = false;
else
    % Even if we disable SDP processing, we're going to move 1x1 SDPs into
    % the nonnegative variable block. It just doesn't make sense to deploy
    % all of that SDP machinery for nonnegative variables.
    sdiag = K.s == 1;
end

% -------------------------------------------------------------------------
% Handle K.ycomplex by splitting apart the complex constraints into pairs
% of real constraints
% -------------------------------------------------------------------------

if ~isempty(K.ycomplex)
    b  = [ real(b) ; imag(b(K.ycomplex)) ];
    At = [ At, 1j * At(:,K.ycomplex) ];
else
    b  = real(b);
end

% -------------------------------------------------------------------------
% Find the locations of the the complex data, so we can convert into
% SeDuMi's internal format, which uses only MATLAB's real representation.
% -------------------------------------------------------------------------
% Strictly speaking, nonnegative variables, the first variable in a Lorentz
% cone, and the first two variables in a rotated Lorentz cone are real. But
% But SeDuMi has allowed them all to be specified as complex; the imaginary
% portions are interpreted as free variables. We have kept that behavior.

% This code replaces whichcpx.c in its entirety. It actually fixes a bug in
% rotated Lorentz cone handling that was probably never exercised.
if isempty(K.xcomplex) && isempty(K.scomplex)
    K.fcplx = zeros(1,0);
    K.qcplx = zeros(1,0);
    K.rcplx = zeros(1,0);
    scplx   = false(1,L_s);
    sreal   = ~sdiag;
    K.rsdpN = L_s;
    N_fc = 0;
    K.cdim = 0;
else
    xc = K.xcomplex;
    tt = xc <= N_fl;
    K.fcplx = xc(tt);
    xc = xc(~tt) -  N_fl;
    tt = xc <= N_q;
    K.qcplx = xc(tt);
    xc = xc(~tt) - N_q;
    tt = xc <= N_r;
    K.rcplx = xc(tt);
    if ~isempty(K.qcplx)
        ndxs = cumsum([1,K.q(1:end-1)]);
        t2 = any(bsxfun(@eq,K.qcplx,ndxs'),1);
        K.fcplx = [ K.fcplx, K.qcplx(t2) + N_fl ];
        K.qcplx(t2) = [];
        t2 = sum(bsxfun(@gt,K.qcplx,ndxs'),1);
        K.q = K.q + full(sparse(1,t2,1,1,L_q));
        K.qcplx = K.qcplx + (1:length(K.qcplx));
        N_q = N_q + length(K.qcplx);
    end
    if ~isempty(K.rcplx)
        ndxs = cumsum([1,K.r(1:end-1)]);
        t2 = any(bsxfun(@eq,K.rcplx,[ndxs,ndxs+1]'),1);
        K.fcplx = [ K.fcplx, K.rcplx(t2) + N_fl + N_q ];
        K.rcplx(t2) = [];
        t2 = sum(bsxfun(@gt,K.rcplx,ndxs'),1);
        K.r = K.r + full(sparse(1,t2,1,1,L_r));
        % This 2*t2 offset is required because QR makes two accesses
        % each of the first two elements of a rotated Lorentz cone.
        K.rcplx = K.rcplx + (1:length(K.rcplx)) + 2 * t2;
        N_r = N_r + length(K.rcplx);
    end
    N_fc = length(K.fcplx);
    N_f  = N_f + N_fc;
    scplx = false(1,L_s);
    scplx(K.scomplex&~sdiag) = true;
    sreal = ~scplx & ~sdiag;
    K.rsdpN = nnz(sreal);
    K.cdim = length(K.xcomplex) + sum(K.s(scplx).^2);
end

% -------------------------------------------------------------------------
% We have significantly rewritten this section of the code. This section
% constructs a sparse matrix that represents the following transformations:
%   --- Free variables split into differences of nonnegative variables; OR
%   --- Free variables placed in a Lorentz-cone
%   --- Rotated lorentz cones translated to standard Lorentz cones
%   --- Lorentz cones rearranged to trace block + norm-bound blocks
%   --- Conversion of diagonal SDPs to nonnegative variables
%   --- SDP coefficients moved to the lower triangle for increased sparsity
% This code replaces the rotlorenz, qreshape, and vectril MEX files.
% -------------------------------------------------------------------------

newL = 0;
newQ = zeros(1,0);
ii = {}; jj = {}; vv = {};

% Split free variables into the difference of nonnegatives
if ~isfield( pars, 'free' ) || pars.free == 2 && L_qrsz
    pars.free = 1;
end
if N_f && ~pars.free
    jt = [ 1 : K.f, K.fcplx ; 1 : K.f, K.fcplx ];
    vt = [ ones(1,K.f), -1j*ones(1,N_fc) ; -ones(1,K.f), 1j*ones(1,N_fc) ];
    ii{end+1} = 1 : 2 * N_f;
    jj{end+1} = jt(:)';
    vv{end+1} = vt(:)';
    newL = 2 * N_f;
    prep.freeL = N_f;
end

% Copy nonnegative variables without change
if K.l
    ii{end+1} = newL + 1 : newL + K.l;
    jj{end+1} = K.f + 1 : K.f + K.l;
    vv{end+1} = ones(1,K.l);
    newL = newL + K.l;
end

% Convert diagonal SDPs to nonnegative variables
if any(sdiag)
    dsize = K.s(sdiag);
    sdpL = sum(dsize);
    prep.sdiag = dsize;
    jstrt = cumsum([N_flqr+1,K.s(1:end-1).^2]);
    jstrt = jstrt(sdiag);
    istrt = cumsum([1,dsize(1:end-1)]);
    dsize = dsize + 1;
    dblks = cumsum(full(sparse(1,istrt,1,1,sdpL)));
    ii{end+1} = newL + 1 : newL + sdpL;
    jj{end+1} = jstrt(dblks) + dsize(dblks) .* ( ( 1 : sdpL ) - istrt(dblks) );
    vv{end+1} = ones(1,sdpL);
    newL = newL + sdpL;
end

% Stuff free variables into a Lorentz cone
tr_off = newL;
nb_off = newL + L_qr;
if N_f && pars.free
    tr_off = tr_off + 1;
    nb_off = nb_off + 1;
    ii{end+1} = nb_off + 1 : nb_off + N_f;
    jj{end+1} = [ 1 : K.f, K.fcplx ];
    vv{end+1} = [ ones(1,K.f), -1j*ones(1,N_fc) ];
    nb_off = nb_off + K.f;
    newQ = N_f + 1;
end

% Rearrange Lorentz cones to trace block + norm-bound blocks
if N_q
    ndxs      = cumsum([1,K.q(1:end-1)]);
    it        = zeros(1,N_q);
    it(ndxs)  = tr_off + 1 : tr_off + L_q;
    it(it==0) = nb_off + 1 : nb_off + ( N_q - L_q );
    jt = K.f + K.l + 1 : K.f + K.l + N_q;
    vt = ones(1,N_q);
    if ~isempty(K.qcplx)
        jt = jt - cumsum(full(sparse(1,K.qcplx,1,1,N_q)));
        vt(K.qcplx) = -1j;
    end
    ii{end+1} = it(:)';
    jj{end+1} = jt;
    vv{end+1} = vt;
    tr_off    = tr_off + L_q;
    nb_off    = nb_off + N_q - L_q;
end

% Transform rotated Lorentz cones to standard Lorentz cones, and rearrange
% to trace block + norm-bound blocks.
if N_r
    ndxr       = cumsum([1,K.r(1:end-1)]);
    ndxp       = ndxr + 2*(0:L_r-1);
    it         = zeros(1,N_r+2*L_r);
    it(ndxp)   = tr_off + 1 : tr_off + L_r;
    it(ndxp+1) = -1;
    it(ndxp+2) = it(ndxp);
    it(it==0)  = nb_off + 1 : nb_off + ( N_r - L_r );
    it(ndxp+1) = it(ndxp+3);
    jt = [ K.f + K.l + N_q + 1, ones(1,N_r+2*L_r-1) ];
    jt([ndxp+1,ndxp+3]) = 0;
    vt = ones(1,N_r+2*L_r);
    vt([ndxp,ndxp+1,ndxp+2]) = sqrt(0.5);
    vt(ndxp+3) = -sqrt(0.5);
    if ~isempty(K.rcplx)
        jt(K.rcplx) = 0;
        vt(K.rcplx) = -1j;
    end
    ii{end+1} = it;
    jj{end+1} = cumsum(jt);
    vv{end+1} = vt;
    nb_off = nb_off + N_r - L_r;
end

% Replace non-diagonal real SDP coefficients with tril(X) + tril(X',-1).
% This cuts the number of nonzeros approximately in half.
if K.rsdpN
    dsize = K.s(sreal);
    sdpL  = sum(dsize.^2);
    jstrt = cumsum([N_flqr+1,K.s(1:end-1).^2]);
    jstrt = jstrt(sreal);
    istrt = cumsum([1,dsize(1:end-1).^2]);
    dblks = cumsum(full(sparse(1,istrt,1,1,sdpL)));
    istrt = istrt + nb_off;
    dsize = dsize(dblks);
    istrt = istrt(dblks);
    jndxs = ( nb_off + 1 : nb_off + sdpL ) - istrt;
    cols  = floor(jndxs ./ dsize);
    rows  = jndxs - dsize .* cols;
    ii{end+1} = max(rows,cols) + min(rows,cols) .* dsize + istrt;
    jj{end+1} = jndxs + jstrt(dblks);
    vv{end+1} = ones(1,sdpL);
    nb_off = nb_off + sdpL;
    clear dsize jstrt istrt dblks jndxs rows cols
end

% Replace Hermitian SDP coefficients with tril(X) + tril(X',-1). This one's
% a bit trickier because we have the real and complex values interleaved,
% and the imaginary values along the diagonal are zero.
if K.rsdpN < length(K.s)
    dsize = K.s(scplx);
    jsize = dsize .^ 2;
    sdpL  = 2 * sum(jsize);
    jstrt = cumsum([N_flqr+1,K.s(1:end-1).^2]);
    jstrt = jstrt(scplx);
    bstrt = cumsum([1,2*jsize(1:end-1)]);
    dblks = cumsum(full(sparse(1,bstrt,1,1,sdpL)));
    istrt = bstrt + nb_off;
    dsize = dsize(dblks);
    istrt = istrt(dblks);
    bndxs = ( nb_off + 1 : nb_off + sdpL ) - istrt;
    cols  = floor( bndxs ./ dsize );
    rows  = bndxs - dsize .* cols;
    imgv  = cols >= dsize;
    cols  = cols - imgv .* dsize;
    indxs = max(rows,cols) + min(rows,cols) .* dsize + imgv .* jsize(dblks) + istrt;
    vals  = ( 1 - 2 * ( cols > rows ) ) .* ( 1 - ( 1 + 1j ) .* imgv );
    keep  = ~imgv | ( rows ~= cols );
    jndxs = rows + cols .* dsize + jstrt(dblks);
    ii{end+1} = indxs(keep);
    jj{end+1} = jndxs(keep);
    vv{end+1} = vals(keep);
    clear dsize jsize jstrt bstrt istrt bndxs rows cols vals imgv keep
end

% Update free, nonnegative, and Lorentz variable counts
K.f = 0;
K.l = newL;
K.q = [ newQ, K.q, K.r ];
K.r = zeros(0,1);
K.s = [ K.s(:,~scplx&~sdiag), K.s(scplx&~sdiag) ];
K.rsdpN = nnz(~scplx&~sdiag);
K.N = K.l + sum(K.q) + sum(K.s(1:K.rsdpN).^2)+2*sum(K.s(K.rsdpN+1:end).^2);

% Create the artificial (x0,z0) variable for the self-dual model by
% appending a zero row to At and c. This is accomplished in QR by adding
% 1 to all of the row indices created above.
K.N  = K.N + 1;
K.l  = K.l + 1;
K.m  = length(b);

% Transform At, c
% The transformation matrix QR does not satisfy QR'*QR=I. But it does, in
% fact, serve as a reverse transformation:
%    --- For Lorentz cones, QR applies a permutation; QR' reverses it.
%    --- For rotated Lorentz cones, QR applies a unitary, self-adjoint
%        rotation on the first two variables, so QR' reverses it.
%    --- For split free variables, QR creates the positive and negative
%        parts; QR' combines them back together.
%    --- For free variables placed in a Lorentz cone, QR moves them to the
%        Lorentz cone block, adding an extra epigraph variable. QR' moves
%        them back and drops the extra variable.
%    --- For semidefinte cones, QR adds the strict upper triangle to the
%        strict lower triangle; QR' copies the strict lower triangle to the
%        strict upper triangle, ensuring symmetry.
%    --- QR adds a row for the self-dual variable; QR' removes it.
[dummy,ndxs] = sort(cellfun(@(x) x(1),jj)); %#ok
QR = sparse( horzcat(ii{ndxs})+1, horzcat(jj{ndxs}), horzcat(vv{ndxs}), K.N, length(c) );
At = real( sparse( QR * At ) );
c  = real( sparse( QR * c  ) );
b  = sparse( b );
prep.QR = QR;
clear ii jj vv

% -------------------------------------------------------------------------
% Now K has field K.{l,q,s}
% Generate a more detailed description of cone K:
% Let K.blkstart(i):K.blkstart(i+1)-1 give the index range of block i.
% Compute maxN=max(order), len=sum(order) for LORENTZ, real PSD, herm PSD
% yields: K.{l,q,s,rsdpN,blkstart,rLen,hLen,qMaxn,rMaxn,hMaxn}
% -------------------------------------------------------------------------
Ksr = K.s(1:K.rsdpN);
Ksc = K.s(K.rsdpN+1:end);
K.blkstart = cumsum([K.l+1,length(K.q)+length(K.r),K.q-1,Ksr.^2,2*Ksc.^2]);
K.rLen = sum(Ksr);
K.hLen = sum(Ksc);
K.qMaxn = max([0,K.q]);
K.rMaxn = max([0,Ksr]);
K.hMaxn = max([0,Ksc]);
K.mainblks = K.blkstart(cumsum([1 1 length(K.q)]));
K.qblkstart = K.blkstart(2:2+length(K.q));  % Also include blkend
K.sblkstart = K.blkstart(2+length(K.q):end);
K.lq = K.mainblks(end)-1;

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function yout = oneinsert(removed_voxels, data)
% Insert ones (for p values) back into data series (i.e., obj.p)
% from which they were removed.
%
% :Usage:
% ::
%
%     yout = oneinsert(removed_voxels, data)
%
% removed_voxels is indicator for the locations of removed voxels
%
% This function is a modified version of naninsert.m, but different.
%
% ..
%    Wani Woo, Aug 2012
% ..
%
% :See also:
% naninsert.m, nanremove.m, zeroinsert.m

yout = zeros(size(removed_voxels,1), size(data,2));
yout(removed_voxels,:) = 1;
yout(~removed_voxels,:) = data;

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function prepData = prepINITModel(origRefModel, taskStruct, spontRxnNames, convertGenes, customRxnsToIgnore, extComp)
% prepINITModel
%
% The purpose of this function is to run time-consuming calculation steps that are not
% dependent on the RNA-Seq data.
%
% origRefModel      The model to use. Expected to be something such as Human-GEM,
%                   Mouse-GEM, etc.
% taskStruct        The essential tasks. Can be loaded with for example
%                   taskStruct = parseTaskList('../data/metabolicTasks_Essential.txt');
% spontRxnNames     The spontaneous rxns. (opt, default {})
% convertGenes      If true the genes are converted to gene names (from
%                   ENSEMBL) (opt, default false)
% customRxnsToIgnore These reactions can be ignored in the ignore mask
%                   (specifying b7=1) (opt, default = {})
% extComp           Name of the external compartment, typically 's' or 'e'. This
%                   is used for identifying exch and import rxns (opt, default = 'e')
% prepData          The resulting prepData structure which is used as input to ftINIT
%
% Usage: prepData = prepINITModel(origRefModel, taskStruct, spontRxnNames, convertGenes, customRxnsToIgnore, extComp)

if nargin < 3
    spontRxnNames = {};
end

if nargin < 4
    convertGenes = false;
end

if nargin < 5
    customRxnsToIgnore = {};
end

if nargin < 6
    extComp = 'e';
end

disp('Step 1: Gene rules')
[origRefModel.grRules, origRefModel.rxnGeneMat] = standardizeGrRules(origRefModel, true);

if convertGenes %For mouse we might want to translate in the opposite direction - this has to be done before calling this function in that case.
    [origRefModel.grRules, origRefModel.genes, origRefModel.rxnGeneMat] = translateGrRules(origRefModel.grRules, 'Name');
end

% Get list of all dead-end/inacessible/constrained reactions in the model.
% We don't merge linear dependent reactions here, that will happen later.
% We also don't remove reversibility here, we don't want that in the final
% models produced. Therefore, this is done as a separate step.
% This takes quite some time to run
disp('Step 2: First simplification')
[~,deletedDeadEndRxns] = simplifyModel(origRefModel,true,false,true,true,true);

% get reduced model
cModel = removeReactions(origRefModel,deletedDeadEndRxns,false,true);

%Then, run the tasks to find out which reactions are essential. These rxns
%will be forced to carry flux in the optimization later

% determine reactions essential for tasks - takes ~10 min
% need to add boundary mets first - but only to a temporary model, we don't
% want those for further processing
disp('Step 3: Check tasks (~10 min)')
if ~isempty(taskStruct)
    bModel = closeModel(cModel);
    [taskReport, essentialRxnMat, ~, essentialFluxes] = checkTasks(bModel,[],true,false,true,taskStruct);

    %extract the essential rxns:
    sel = sum(essentialRxnMat,2) > 0;
    selInd = find(sel);
    essentialRxns = bModel.rxns(selInd);%382 rxns for human-GEM

    %Find metabolites present in taskStruct. We want to avoid removing
    %these metabolites from the final model (even though they may no longer
    %participate in any reacitons) so that the final model is still able to
    %complete all of the tasks without any errors.
    taskMets = union(vertcat(taskStruct.inputs),vertcat(taskStruct.outputs));
    taskMets = union(taskMets, parseRxnEqu(vertcat(taskStruct.equations)));
    modelMets = strcat(cModel.metNames,'[',cModel.comps(cModel.metComps),']');
    [inModel,metInd] = ismember(taskMets,modelMets);
    essentialMetsForTasks = cModel.mets(metInd(inModel));%118 mets

    %Remove tasks that cannot be performed
    taskStruct(taskReport.ok == false) = [];

else
    essentialRxns = {};
    essentialMetsForTasks = {};
    taskReport = {};
end

% remove metabolites separately to avoid removing those needed for tasks
unusedMets = cModel.mets(all(cModel.S == 0,2));%1248 in Human-GEM
cModel = removeMets(cModel, setdiff(unusedMets,essentialMetsForTasks));

if ~isempty(taskStruct)

    %Here, we decide on a direction in which the essential reactions should be forced.
    essentialRevDir = false(length(essentialRxns),1);
    pp = zeros(length(essentialRxns),1);
    nn = zeros(length(essentialRxns),1);
    for i = 1:length(selInd)
        pos = sum(essentialFluxes(selInd(i),essentialRxnMat(selInd(i),:)) > 0);
        neg = sum(essentialFluxes(selInd(i),essentialRxnMat(selInd(i),:)) < 0);
        essentialRevDir(i) = pos < neg;
        pp(i) = pos;
        nn(i) = neg;
    end

    %sum(pp>0 & nn>0) %0, so it is not a big problem that different tasks share essential
    %reactions but in different directions, at least no such cases exist now, so we ignore this for now.

    %Now create a minimal model for the MILP - we want to make this as small as we can, but we
    %don't want to use this later when we create the final models - just for the first MILP.
    %We don't want boundary metabolites in here either.

    %We do a trick with the essential rxns here:
    %What we want is to get rid of the reversibility of the reactions, because reversible reactions require
    %an integer each in the MILP, and we can also get the case that it may be cheaper to force flux in the wrong direction, which
    %will probably not give the gap-filling of reactions that we want.
    %So, we simply change the reversible reactions to irreversible, and flip their direction if the direction
    %to force flux in is negative. We don't need to do anything else, runINIT will then handle all essential rxns
    %as irreversible (and will not add any integers etc. for them).

    minModel1 = cModel;
    %for the essential rxns where we should force flux in the positive direction, just make them irreversible
    %for the negative direction, we do the same, but also flip the direction of the reaction

    %first flip the reactions with negative direction
    %constructEquations(minModel1,minModel1.rxns(selInd(essentialRevDir))) %looks good
    minModel1 = reverseRxns(minModel1, minModel1.rxns(selInd(essentialRevDir)));
    %constructEquations(minModel2,minModel2.rxns(selInd(essentialRevDir))) %looks good

    %Then make them irreversible
    minModel1.rev(selInd) = 0;
    minModel1.lb(selInd) = 0;
else
    minModel1 = cModel;
end
%Now, remove reversibility on any reactions that can only be run in one direction - that simplifies the problem
%This takes a long time to run (at least an hour or so).
%This is an important step, we go from 5533 to 3734 reversible rxns, which reduces the size of the MILP,
%since irreversible rxns are easier to handle.
disp('Step 4: Second simplification (~1 hour)')
minModel2 = simplifyModel(minModel1,false,false,false,false,false,false,true);

disp('Step 5: Final work')

%Now, merge all linear dependent rxns - this reduces the size of the model a lot.
%Size goes from 11888 rxns to 7922 rxns for human-GEM, so this is an important step
[minModel3,origRxnIds,groupIds]=mergeLinear(minModel2, {});

%we now need to figure out what happened to the essential rxns in the merge step above.
%Note that it would be ideal to run checkTasks on the result from mergeLinear, but that doesn't work.
%The exchange rxns have in many cases been merged with other reactions, making it difficult to run
%the tasks. So, instead we run checkTasks before, and figure out which merged reactions the
%essential reactions belong to now

%It may be possible to create a minModel2b, run mergeLinear on that (which is quick), and
%then run checkTasks on the minimized model. Would be good to check if this gives the
%same result, I'm not sure. The code would get less complicated.

if ~isempty(taskStruct)

    %find all potential rxns
    rxnIndOrig = ismember(origRxnIds, essentialRxns);
    ids = unique(groupIds(rxnIndOrig));
    ids = ids(ids > 0);%0 means it has not been merged
    rxnCandidates = unique([essentialRxns;origRxnIds(ismember(groupIds, ids))]);%388, essentialRxns is 382
    newEssentialRxns = minModel3.rxns(ismember(minModel3.rxns,rxnCandidates));%196 in Human-GEM
    %setdiff(newEssentialRxns, essentialRxns)%only 2, that is very few
    %setdiff(essentialRxns, newEssentialRxns)%188, that is very many
    %It is probably the case that many essential reactions are linearly dependent on one another, which makes sense.
else
    newEssentialRxns = {};
end
%What is left now is to remove some rxns from the problem, such as exchange rxns and spontaneous rxns.
%We set the rxn scores of these to exactly 0, and make sure no other scores are zero.
%runINIT will then handle this, we do >0 and <0 to get positive and negative scores, the rest
%are not included in the problem.
%But, this has to happen in the scoring - we identify the reactions in the origrxns here, and then fix this
%in the scoring for each sample (the function groupRxnScores handles this and is run for each sample).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%identify the reactions to ignore in the minModel2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Identify reactions that should be ignored by tINIT, i.e. tINIT will not remove these regardless of score

[~,exchRxnInd] = getExchangeRxns(minModel2);
toIgnoreExch = false(numel(minModel2.rxns),1);
toIgnoreExch(exchRxnInd) = true;
toIgnoreImportRxns = false(numel(minModel2.rxns),1);
toIgnoreSimpleTransp = false(numel(minModel2.rxns),1);
toIgnoreAdvTransp = false(numel(minModel2.rxns),1);
toIgnoreS = false(numel(minModel2.rxns),1);

%sum(toIgnore) %1500 for Human-GEM
%find simple transport reactions with no GPRs:
%We really only skip the reactions with two metabolites
%where they are the same but different compartments - only skip the transport into the cell for now
numMets = sum(minModel2.S ~= 0, 1);
scomp = find(strcmp(minModel2.comps,extComp));
for i = 1:length(minModel2.rxns)
    %check that it has no GPR and that the rxn has two metabolites
    if strcmp(minModel2.grRules(i),'') && (numMets(i) == 2)
      metsComps = minModel2.metComps(minModel2.S(:,i) ~= 0);
      metNames = minModel2.metNames(minModel2.S(:,i) ~= 0);
      if metsComps(1) ~= metsComps(2) && (metsComps(1) == scomp || metsComps(2) == scomp) %different compartments, one is the extComp
          if strcmp(metNames{1}, metNames{2}) %the same metabolite, it is a transport reaction
              toIgnoreImportRxns(i) = true;
          end
      elseif metsComps(1) ~= metsComps(2) %different compartments, no check for extComp
          if strcmp(metNames{1}, metNames{2}) %the same metabolite, it is a transport reaction
              toIgnoreSimpleTransp(i) = true;
          end
      end
    else %check for advanced
        metsInd = find(minModel2.S(:,i) ~= 0);
        %check that we have an even number of mets that is larger than 2 and that we don't have a GPR
        if rem(length(metsInd),2) == 0 && length(metsInd) > 2 && isempty(minModel2.grRules{i})
            %Now check that we have pairs of metabolites that are transported
            SVals = full(minModel2.S(metsInd,i));
            metNames = minModel2.metNames(metsInd);
            comps = minModel2.metComps(metsInd);
            success = true;
            while ~isempty(metNames)
                if length(metNames) < 2 %this should not happen
                   error('We should never arrive here!');
                end
                metMatch = find(strcmp(metNames(2:length(metNames)),metNames{1})) + 1;
                if length(metMatch) ~= 1 %we want a match, and only one (if there is some other complex transport rxn, we skip that)
                   success = false;
                   break;
                end
                if (SVals(1) + SVals(metMatch)) ~= 0 %the stoichiometry is not right
                   success = false;
                   break;
                end
                if comps(1) == comps(metMatch) %they must be in different compartments (maybe this need not to be checked)
                   success = false;
                   break;
                end
                %now remove the pair and continue with the next
                SVals([1;metMatch]) = [];
                metNames([1;metMatch]) = [];
                comps([1;metMatch]) = [];
            end
            toIgnoreAdvTransp(i) = success;
        end
    end
end

%Spontaneous reactions:
toIgnoreSpont = ismember(minModel2.rxns, spontRxnNames);%only 10 rxns in human-GEM

%rxns without GPRs in the s compartment (outside the cell)
sCompInd = find(strcmp(minModel2.comps, extComp));
for i = 1:length(minModel2.rxns)
    metsInd = find(minModel2.S(:,i) ~= 0);
    %check that it has more than 1 metabolite (not exch rxn) and no GPR
    if length(metsInd) > 1 && isempty(minModel2.grRules{i})
        if  all(minModel2.metComps(metsInd) == sCompInd)
            toIgnoreS(i) = true;
        end
    end
end

toIgnoreAllWithoutGPRs = cellfun(@isempty,minModel2.grRules);

%sum(toIgnore)%in total 2375 rxns in human-GEM
%sum(toIgnoreAllTransp)%in total 3337 rxns in human-GEM, so 1,000 more

%Now, try to scale the model to become more favorable for the solver.
%In general, we try to make all fluxes as similar as possible.
%      There is room for improvement here, this function is pretty simple.
%      It only rescales reactions, while it would be possible to rescale
%      metabolites as well (ROS => kROS, albumin => millialbumin, etc., but
%      scaled freely with a mathematical method). It could potentially open up
%      for using less strict margins in the MILP, and make it run faster.
scaledMinModel = rescaleModelForINIT(minModel3);
scaledMinModel.ub(scaledMinModel.ub > 0) = 1000;
scaledMinModel.lb(scaledMinModel.lb < 0) = -1000;

% create data structure with pre-processing results
prepData.taskReport = taskReport;
prepData.essentialRxns = newEssentialRxns; %essential rxns in the minModel
prepData.taskStruct = taskStruct;
prepData.refModel = cModel;
prepData.minModel = scaledMinModel;
prepData.refModelWithBM = closeModel(cModel); %do this here so we don't have to do this for each sample
prepData.groupIds = groupIds;
prepData.essentialMetsForTasks = essentialMetsForTasks;

prepData.toIgnoreExch = toIgnoreExch;
prepData.toIgnoreImportRxns = toIgnoreImportRxns;
prepData.toIgnoreSimpleTransp = toIgnoreSimpleTransp;
prepData.toIgnoreAdvTransp = toIgnoreAdvTransp;
prepData.toIgnoreSpont = toIgnoreSpont;
prepData.toIgnoreS = toIgnoreS;
prepData.toIgnoreCustomRxns = ismember(minModel2.rxns, customRxnsToIgnore);
prepData.toIgnoreAllWithoutGPRs = toIgnoreAllWithoutGPRs;

end

author: SysBioChalmers

repo: RAVEN


text:
function output = CheckIn(input,array)
[a_rows,a_cols] = size(array);
output = 0;
for i = 1:a_cols
    if input == array(i)
        output = 1;
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Fig. 8.10   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

numG=1;
denG=[1 0 0];

numD=[1 .2];
denD=[1 2];

num=conv(numG,numD);
den=conv(denG,denD);
poles=roots(den);
zeros=roots(num);

K1=0:.05:1.22;
K2=[1.25 1.28];  % K for break-in and break-away points
K3=1.5:5:100;
K=[K1 K2 K3];
Ko=.81;

r=rlocus(num,den,K);
ro=rlocus(num,den,Ko);

plot(r,'-'),grid
axis('square')
axis([-2.5 .5 -1.5 1.5])
hold on
plot(ro,'k*')
plot(-.2,0,'o')
plot(-2,0,'x')
plot(0,.01,'x')
plot(0,-.01,'x')
title('Fig. 8.10  s-plane locus vs. K')
xlabel('Re(s)')
ylabel('Im(s)')
text(-2.3,-.7,'*  K_c = 0.81')
hold off


text:
function x = callback_inpainting(y,dir,options)

% callback_inpainting - callback for inpainting applications
%
%   x = callback_inpainting(y,dir,options);
%
%   You must provide options.mask.
%   The missing pixels (removed) are mask==1.
%
%   Copyright (c) 2008 Gabriel Peyre

mask = getoptions(options, 'mask', [], 1);

x = y;
x(mask==1) = 0;

author: gpeyre

repo: matlab-toolboxes


text:
function DAT = read_excel(excelfilename, has_header_row)
% Reads each sheet of an Excel input file into a data structure DAT
%
% :Usage:
% ::
%
%     DAT = read_excel(excelfilename, has_header_row [1 or 0])
%
%  - reads multiple sheets
%  - adds fields named with variable names for easy access
%  - uses importdata <matlab internal> to do most of the work
%
% ..
%     Author and copyright information:
%     -------------------------------------------------------------------------
%     Copyright (C) 2015 Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Input:
%
%   **excelfilename:**
%        Char array with full path of input filename
%
% :Outputs:
%
%   **DAT:**
%        Structure with DAT.(sheetname).(varname) data fields
%        DAT.(sheetname).varnames is cell array of var names
%
%   **behdat:**
%        Data structure as read in by importdata
%
% :Examples:
% ::
%
%    excelfilename = fullfile(basedir, 'data', 'SchulzNewcorn_regressor_demograph.xlsx');
%    has_header_row = 1;
%    [DAT, behdat] = read_excel(excelfilename, has_header_row)
%
% :See also:
% importdata, read_database, read_database2, read_edat_output_2008, read_physio_data
%
% ..
%    Programmers' notes:
%    Created by Tor, 7/2015
%    Updated by Tor, Dec 2021 - keep up with changes in importdata behavior
% ..

DAT = struct;

% Note: importdata is not handling blanks in Excel files well. Use xlsread below.
behdat = importdata(excelfilename); % read the original data

has_multiple_sheets = isstruct(behdat.data) || isstruct(behdat.textdata);

if has_multiple_sheets

    % % Sheet names
    % %-------------------
    mysheets = fieldnames(behdat.data);

    for s = 1:length(mysheets)

        [~, ~, raw] = xlsread(excelfilename, s);  % 2021: it doesn't extract data well with blank entries. Use raw and do it ourselves.

        [data_struct, data_table] = process_sheet(raw, has_header_row);

        DAT.(mysheets{s}).raw = raw;
        DAT.(mysheets{s}).data_struct = data_struct;
        DAT.(mysheets{s}).data_table = data_table;

        % DAT.(mysheets{s}).data_table = process_sheet(behdat.data.(mysheets{s}), behdat.textdata.(mysheets{s}), has_header_row);

    end

else
    % One sheet
    [~, ~, raw] = xlsread(excelfilename, 1);

    DAT.raw = raw;
    [DAT.data_struct, DAT.data_table] = process_sheet(raw, has_header_row);

    %     DAT.data_table = process_sheet(behdat.data, behdat.textdata, has_header_row);
end

end % main function

% ------------------------------------------------------------------------
% Subfunctions
% ------------------------------------------------------------------------

function [data_struct, data_table] = process_sheet(raw, has_header_row)
% Create a table object from raw data from xlsread

% % Variable names
if has_header_row

    varnames = raw(1, :);

else
    varnames = cell(1, size(raw, 2));  % textdata only seems to be full-width

    for i = 1:length(varnames)
        varnames{i} = sprintf('Var_%03d', i);
    end

end

varnames = matlab.lang.makeValidName(varnames);

% % Loop through columns
% --------------------

data_struct = struct();

for i = 1:length(varnames)

    fname = varnames{i};

    % figure out if is text and process (replace empty)
    mydat = raw((has_header_row+1):end, i); % skip first row if headers.

    % replace empty, convert to numeric if possible
    mydat = process_data_column(mydat);

    data_struct.(fname) = mydat;

end

data_table = struct2table(data_struct);

end

function mydat = process_data_column(mydat)

% Are some characters text?
is_char = any(cellfun(@ischar, mydat));

if ~is_char
    % already numeric. convert from cell.
    mydat = cat(1, mydat{:});
    return
end

% Below, some entries at least are char arrays. But may be numeric with
% missing values/empty cells

% some empty cells in text are read in as NaN, which can cause errors below
nanvec = false(size(mydat));
for i = 1:length(mydat)
    if isnan(mydat{i}), nanvec(i) = true; else nanvec(i) = false; end
end
mydat(nanvec) = {'NaN'};

% try to replace empty with 'NaN'. We will convert to numeric later.
wh_empty = cellfun(@isempty, mydat);
mydat(wh_empty) = {'NaN'};

% sometimes will read in mixed format: char and numeric.
% if so, we need to convert ALL to char.
% all-numeric is handled above, so here convert if mixed.
is_numeric = cellfun(@isnumeric, mydat);
for i = 1:length(mydat)
    if is_numeric(i)
        mydat{i} = num2str(mydat{i});
    end
end
% mydat{is_numeric} = cellfun(@num2str, mydat(is_numeric));

% try to convert to numeric
%     numdat = cellfun(@str2num, mydat);
numdat = cellfun(@str2num, mydat, 'UniformOutput', false);

% if any empty after replacing blanks, this is text (leave as cell)
istext = any(cellfun(@isempty, numdat));

if istext
    % Leave as cell array, but make NaN empty char
    mydat(wh_empty) = {' '};

else  % numeric
    mydat = cat(1, numdat{:});

end

end % process_data_column

%
%
%
% function data_table = process_sheet(data, textdata, has_header_row)
% % Create a table object from separate numeric data (data) and text data (textdata).
%
% % % Variable names
% if has_header_row
%
%     varnames = textdata(1, :);
%
% else
%     varnames = cell(1, size(textdata, 2));  % textdata only seems to be full-width
%
%     for i = 1:length(varnames)
%         varnames{i} = sprintf('Var_%03d', i);
%     end
%
% end
%
% varnames = matlab.lang.makeValidName(varnames);
%
% % % Get vector of text vs. numeric
% % --------------------
%
%
% % % All data
% % --------------------
%
% data_struct = struct();
%
% numeric_counter = 1;        % Numeric array has reduced number of columns
%
% for i = 1:length(varnames)
%
%     fname = varnames{i};
%
%     % figure out if is text
%     mytext = textdata((has_header_row+1):end, i); % skip first row if headers.
%     istext = ~all(cellfun(@isempty, mytext));
%
%     if istext
%
%         data_struct.(fname) = mytext;
%
%     else  % numeric
%
%         data_struct.(fname) = data(:, numeric_counter);
%         numeric_counter = numeric_counter + 1;
%
%     end
%
% end
%
% end
%


text:
function res = vl_simplenn_globalModel(net, x, dzdy, res, varargin)
% VL_SIMPLENN  Evaluates a simple CNN
%   RES = VL_SIMPLENN(NET, X) evaluates the convnet NET on data X.
%   RES = VL_SIMPLENN(NET, X, DZDY) evaluates the convnent NET and its
%   derivative on data X and output derivative DZDY.
%
%   The network has a simple (linear) topology, i.e. the computational
%   blocks are arranged in a sequence of layers. Please note that
%   there is no need to use this wrapper, which is provided for
%   convenience. Instead, the individual CNN computational blocks can
%   be evaluated directly, making it possible to create significantly
%   more complex topologies, and in general allowing greater
%   flexibility.
%
%   The NET structure contains two fields:
%
%   - net.layers: the CNN layers.
%   - net.normalization: information on how to normalize input data.
%
%   The network expects the data X to be already normalized. This
%   usually involves rescaling the input image(s) and subtracting a
%   mean.
%
%   RES is a structure array with one element per network layer plus
%   one representing the input. So RES(1) refers to the zeroth-layer
%   (input), RES(2) refers to the first layer, etc. Each entry has
%   fields:
%
%   - res(i+1).x: the output of layer i. Hence res(1).x is the network
%     input.
%
%   - res(i+1).aux: auxiliary output data of layer i. For example,
%     dropout uses this field to store the dropout mask.
%
%   - res(i+1).dzdx: the derivative of the network output relative to
%     variable res(i+1).x, i.e. the output of layer i. In particular
%     res(1).dzdx is the derivative of the network output with respect
%     to the network input.
%
%   - res(i+1).dzdw: the derivative of the network output relative to
%     the parameters of layer i. It can be a cell array for multiple
%     parameters.
%
%   net.layers is a cell array of network layers. The following
%   layers, encapsulating corresponding functions in the toolbox, are
%   supported:
%
%   Convolutional layer::
%     The convolutional layer wraps VL_NNCONV(). It has fields:
%
%     - layer.type = 'conv'
%     - layer.weights = {filters, biases}
%     - layer.stride: the sampling stride (usually 1).
%     - layer.padding: the padding (usually 0).
%
%   Max pooling layer::
%     The max pooling layer wraps VL_NNPOOL(). It has fields:
%
%     - layer.type = 'pool'
%     - layer.method: pooling method ('max' or 'avg').
%     - layer.pool: the pooling size.
%     - layer.stride: the sampling stride (usually 1).
%     - layer.padding: the padding (usually 0).
%
%   Normalization layer::
%     The normalization layer wraps VL_NNNORMALIZE(). It has fields
%
%     - layer.type = 'normalize'
%     - layer.param: the normalization parameters.
%
%   Spatial normalization layer:
%     This is similar to the layer above, but wraps VL_NNSPNORM():
%
%     - layer.type = 'spnorm'
%     - layer.param: the normalization parameters.
%
%   Batch normalization layer:
%     This layer wraps VL_NNBNORM(). It has fields:
%
%     - layer.type = 'bnorm'
%     - layer.weights = {multipliers, biases}.
%
%   ReLU and Sigmoid layers::
%     The ReLU layer wraps VL_NNRELU(). It has fields:
%
%     - layer.type = 'relu'
%
%     The sigmoid layer is the same, but for the sigmoid function, with
%     `relu` replaced by `sigmoid`.
%
%   Dropout layer::
%     The dropout layer wraps VL_NNDROPOUT(). It has fields:
%
%     - layer.type = 'dropout'
%     - layer.rate: the dropout rate.
%
%   Softmax layer::
%     The softmax layer wraps VL_NNSOFTMAX(). It has fields
%
%     - layer.type = 'softmax'
%
%   Log-loss layer::
%     The log-loss layer wraps VL_NNLOSS(). It has fields:
%
%     - layer.type = 'loss'
%     - layer.class: the ground-truth class.
%
%   Softmax-log-loss layer::
%     The softmax-log-loss layer wraps VL_NNSOFTMAXLOSS(). It has
%     fields:
%
%     - layer.type = 'softmaxloss'
%     - layer.class: the ground-truth class.
%
%   P-dist layer:
%     The pdist layer wraps VL_NNPDIST(). It has fields:
%
%     - layer.type = 'pdist'
%     - layer.p = P parameter of the P-distance
%     - layer.noRoot = whether to raise the distance to the P-th power
%     - layer.epsilon = regularization parameter for the derivatives
%
%   Custom layer::
%     This can be used to specify custom layers.
%
%     - layer.type = 'custom'
%     - layer.forward: a function handle computing the block.
%     - layer.backward: a function handle computing the block derivative.
%
%     The first function is called as res(i+1) = forward(layer, res(i), res(i+1))
%     where res() is the struct array specified before. The second function is
%     called as res(i) = backward(layer, res(i), res(i+1)). Note that the
%     `layer` structure can contain additional fields if needed.

% Copyright (C) 2014 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.res = [] ;
opts.conserveMemory = false ;
opts.sync = false ;
opts.disableDropout = false ;
opts.freezeDropout = false ;
opts.accumulate = false;
opts.backPropDepth = +inf ;

opts = vl_argparse(opts, varargin);

n = numel(net.layers) ;

if (nargin <= 2) || isempty(dzdy)
    doder = false ;
else
    doder = true ;
end

gpuMode = isa(x, 'gpuArray') ;

if nargin <= 3 || isempty(res)
    res = struct(...
        'x', cell(1,n+1), ...
        'dzdx', cell(1,n+1), ...
        'dzdw', cell(1,n+1), ...
        'aux', cell(1,n+1), ...
        'time', num2cell(zeros(1,n+1)), ...
        'backwardTime', num2cell(zeros(1,n+1))) ;
end
res(1).x = x ;

for i=1:n
    l = net.layers{i} ;
    res(i).time = tic ;
    switch l.type
        case 'conv'
            res(i+1).x = vl_nnconv(res(i).x, l.weights{1}, l.weights{2}, 'pad', l.pad, 'stride', l.stride) ;
        case 'pool'
            res(i+1).x = vl_nnpool(res(i).x, l.pool, 'pad', l.pad, 'stride', l.stride, 'method', l.method) ;
        case 'normalize'
            res(i+1).x = vl_nnnormalize(res(i).x, l.param) ;
        case 'softmax'
            res(i+1).x = vl_nnsoftmax(res(i).x) ;
        case 'loss'
            res(i+1).x = vl_nnloss(res(i).x, l.class) ;
        case 'softmaxloss'
            res(i+1).x = vl_nnsoftmaxloss(res(i).x, l.class) ;
        case 'svmloss_multiclass'
            res(i+1).x = vl_nnsvmloss(res(i).x, l.class) ;
        case 'relu'
            res(i+1).x = vl_nnrelu(res(i).x) ;
        case 'sigmoid'
            res(i+1).x = vl_nnsigmoid(res(i).x) ;
        case 'noffset'
            res(i+1).x = vl_nnnoffset(res(i).x, l.param) ;
        case 'spnorm'
            res(i+1).x = vl_nnspnorm(res(i).x, l.param) ;
        case 'reshape'
            res(i+1).x = vl_nnreshape(res(i).x, [l.numOutput l.numClasses]);
        case 'dropout'
            if opts.disableDropout
                res(i+1).x = res(i).x ;
            elseif opts.freezeDropout
                [res(i+1).x, res(i+1).aux] = vl_nndropout(res(i).x, 'rate', l.rate, 'mask', res(i+1).aux) ;
            else
                [res(i+1).x, res(i+1).aux] = vl_nndropout(res(i).x, 'rate', l.rate) ;
            end
        case 'bnorm'
            if isfield(l, 'weights')
                res(i+1).x = vl_nnbnorm(res(i).x, l.weights{1}, l.weights{2}) ;
            else
                res(i+1).x = vl_nnbnorm(res(i).x, l.filters, l.biases) ;
            end
        case 'pdist'
            res(i+1) = vl_nnpdist(res(i).x, l.p, 'noRoot', l.noRoot, 'epsilon', l.epsilon) ;
        case 'custom'
            res(i+1) = l.forward(l, res(i), res(i+1)) ;
        otherwise
            error('Unknown layer type %s', l.type) ;
    end
    % optionally forget intermediate results
    forget = opts.conserveMemory ;
    forget = forget & (~doder || strcmp(l.type, 'relu')) ;
    forget = forget & ~(strcmp(l.type, 'loss') || strcmp(l.type, 'softmaxloss')) ;
    forget = forget & (~isfield(l, 'rememberOutput') || ~l.rememberOutput) ;
    if forget
        res(i).x = [] ;
    end
    if gpuMode & opts.sync
        % This should make things slower, but on MATLAB 2014a it is necessary
        % for any decent performance.
        wait(gpuDevice) ;
    end
    res(i).time = toc(res(i).time) ;
end

if doder
    res(n+1).dzdx = dzdy ;
    for i=n:-1:max(1, n-opts.backPropDepth+1)
        l = net.layers{i} ;
        res(i).backwardTime = tic ;
        switch l.type
            case 'conv'
                if ~opts.accumulate
                    [res(i).dzdx, res(i).dzdw{1}, res(i).dzdw{2}] = ...
                        vl_nnconv(res(i).x, l.weights{1}, l.weights{2}, ...
                        res(i+1).dzdx, ...
                        'pad', l.pad, 'stride', l.stride) ;
                else
                    dzdw = cell(1,2) ;
                    [res(i).dzdx, dzdw{1}, dzdw{2}] = ...
                        vl_nnconv(res(i).x, l.weights{1}, l.weights{2}, ...
                        res(i+1).dzdx, ...
                        'pad', l.pad, 'stride', l.stride) ;
                    for j=1:2
                        res(i).dzdw{j} = res(i).dzdw{j} + dzdw{j} ;
                    end
                    clear dzdw ;
                end

            case 'pool'
                res(i).dzdx = vl_nnpool(res(i).x, l.pool, res(i+1).dzdx, ...
                    'pad', l.pad, 'stride', l.stride, 'method', l.method) ;
            case 'normalize'
                res(i).dzdx = vl_nnnormalize(res(i).x, l.param, res(i+1).dzdx) ;
            case 'softmax'
                res(i).dzdx = vl_nnsoftmax(res(i).x, res(i+1).dzdx) ;
            case 'loss'
                res(i).dzdx = vl_nnloss(res(i).x, l.class, res(i+1).dzdx) ;
            case 'softmaxloss'
                res(i).dzdx = vl_nnsoftmaxloss(res(i).x, l.class, res(i+1).dzdx) ;
            case 'reshape'
                res(i).dzdx = vl_nnreshape(res(i).x, [l.numOutput l.numClasses], res(i+1).dzdx);
            case 'svmloss_multiclass'
                res(i).dzdx = vl_nnsvmloss(res(i).x, l.class, res(i+1).dzdx) ;
            case 'relu'
                if ~isempty(res(i).x)
                    res(i).dzdx = vl_nnrelu(res(i).x, res(i+1).dzdx) ;
                else
                    % if res(i).x is empty, it has been optimized away, so we use this
                    % hack (which works only for ReLU):
                    res(i).dzdx = vl_nnrelu(res(i+1).x, res(i+1).dzdx) ;
                end
            case 'sigmoid'
                res(i).dzdx = vl_nnsigmoid(res(i).x, res(i+1).dzdx) ;
            case 'noffset'
                res(i).dzdx = vl_nnnoffset(res(i).x, l.param, res(i+1).dzdx) ;
            case 'spnorm'
                res(i).dzdx = vl_nnspnorm(res(i).x, l.param, res(i+1).dzdx) ;
            case 'dropout'
                if opts.disableDropout
                    res(i).dzdx = res(i+1).dzdx ;
                else
                    res(i).dzdx = vl_nndropout(res(i).x, res(i+1).dzdx, 'mask', res(i+1).aux) ;
                end
            case 'bnorm'
                if ~opts.accumulate
                    if isfield(l, 'weights')
                        [res(i).dzdx, res(i).dzdw{1}, res(i).dzdw{2}] = ...
                            vl_nnbnorm(res(i).x, l.weights{1}, l.weights{2}, ...
                            res(i+1).dzdx) ;
                    else
                        [res(i).dzdx, res(i).dzdw{1}, res(i).dzdw{2}] = ...
                            vl_nnbnorm(res(i).x, l.filters, l.biases, ...
                            res(i+1).dzdx) ;
                    end
                else
                    dzdw = cell(1,2) ;
                    if isfield(l, 'weights')
                        [res(i).dzdx, dzdw{1}, dzdw{2}] = ...
                            vl_nnbnorm(res(i).x, l.weights{1}, l.weights{2}, ...
                            res(i+1).dzdx) ;
                    else
                        [res(i).dzdx, dzdw{1}, dzdw{2}] = ...
                            vl_nnbnorm(res(i).x, l.filters, l.biases, ...
                            res(i+1).dzdx) ;
                    end
                    for j=1:2
                        res(i).dzdw{j} = res(i).dzdw{j} + dzdw{j} ;
                    end
                    clear dzdw ;
                end
            case 'pdist'
                res(i).dzdx = vl_nnpdist(res(i).x, l.p, res(i+1).dzdx, ...
                    'noRoot', l.noRoot, 'epsilon', l.epsilon) ;
            case 'custom'
                res(i) = l.backward(l, res(i), res(i+1)) ;
        end
        if opts.conserveMemory
            res(i+1).dzdx = [] ;
        end
        if gpuMode & opts.sync
            wait(gpuDevice) ;
        end
        res(i).backwardTime = toc(res(i).backwardTime) ;
    end
end

author: aosokin


text:
function p = motionProfile(profileType, t, x, ptotal)
% motionProfile  Compute motion profiles for the crane block
%   PROF = motionProfile(PROFILETYPE, T, X, PTOTAL) returns the profile of
%     the crane block.
%     PROFILETYPE : 'acceleration', 'velocity', 'profile'
%     T           : time vector
%     X           : 3-element vector [tp1, tp2, tf]
%     PTOTAL      : total displacement
%
%   The acceleration profile for this problem is a "bang-coast-bang"
%   profile.
%
%  Example:
%    x = [10 5 20];
%    ptotal = 50;
%    t = 0:.01:30;
%
%    acc = motionProfile('accel', t, x, ptotal);
%    vel = motionProfile('vel', t, x, ptotal);
%    pos = motionProfile('pos', t, x, ptotal);

% Copyright 2011 The MathWorks, Inc.

error(nargchk(4, 4, nargin, 'struct'));

profileType = validatestring(profileType, ...
   {'acceleration', 'velocity', 'position'}, mfilename, 'PROFILETYPE');

tp1 = x(1); tp2 = x(2); tf = x(3);

% Pulse heights
p1 = 2*ptotal/(tp1*(2*tf-tp1-tp2));
p2 = tp1/tp2*p1;

switch profileType
   case 'acceleration'
      p = zeros(size(t));

      % Acceleration phase
      p(t>=0 & t<tp1) = p1;

      % Deceleration phase
      p(t>=tf-tp2 & t<tf) = -p2;

   case 'velocity'
      p = zeros(size(t));

      % Acceleration phase
      idx1 = t>=0 & t<tp1;
      p(idx1) = t(idx1)*p1;

      % Coast
      p(t>=tp1 & t<tf-tp2) = tp1*p1;

      % Deceleration phase
      idx2 = t>=tf-tp2 & t<tf;
      p(idx2) = tp1*p1 - (t(idx2)-t(find(idx2, 1, 'first')))*p2;

   case 'position'
      p = ptotal*ones(size(t));

      % Acceleration phase
      idx1 = t>=0 & t<tp1;
      p(idx1) = 0.5*p1*t(idx1).^2;

      % Coast phase
      idx2 = t>=tp1 & t<tf-tp2;
      idx2_f = find(idx2, 1, 'first');
      p(idx2) = 0.5*p1*t(idx2_f).^2 + (t(idx2)-t(idx2_f))*tp1*p1;

      % Deceleration phase
      idx3 = t>=tf-tp2 & t<tf;
      idx3_f = find(idx3, 1, 'first');
      p(idx3) = 0.5*p1*t(idx2_f).^2 + (t(idx3_f)-t(idx2_f))*tp1*p1 +...
         tp1*p1*(t(idx3)-t(idx3_f)) - 0.5*p2*(t(idx3)-t(idx3_f)).^2;

end


text:
function ply_to_tri_surface_test ( )

%*****************************************************************************80
%
%% PLY_TO_TRI_SURFACE_TEST tests the PLY_TO_TRI_SURFACE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PLY_TO_TRI_SURFACE_TEST\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Test the PLY_TO_TRI_SURFACE library.\n' );

  ply_file_name = 'sphere.ply';

  ply_to_tri_surface_test01 ( ply_file_name );

  ply_file_name = 'cube.ply';

  ply_to_tri_surface_test02 ( ply_file_name );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PLY_TO_TRI_SURFACE_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function figure_num = edge_test37 ( figure_num )

%*****************************************************************************80
%
%% EDGE_TEST037 plots a function with a derivative discontinuity.
%
%  Discussion:
%
%    This is example 3.1 in the reference.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Rick Archibald, Anne Gelb, Jungho Yoon,
%    Determining the locations and discontinuities in the derivatives
%    of functions,
%    Applied Numerical Mathematics,
%    Volume 58, 2008, pages 577-592.
%
  if ( nargin < 1 )
    figure_num = 0;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'EDGE_TEST037:\n' );
  fprintf ( 1, '  Plot 2D test function #5, discontinuous medium wave, U(x,t).\n' );

  figure_num = figure_num + 1;
  figure ( figure_num )

  m = 21;
  x = linspace ( -1.0, 0.0, m );
  n = 21;
  y = linspace ( 0.0, 0.1, n );
  [ X, Y ] = meshgrid ( x, y );
  tri = delaunay ( X, Y );
  Z = fxy5 ( m * n, X, Y );
  trimesh ( tri, X, Y, Z, 'FaceColor', 'Interp', 'EdgeColor', 'None' );

  grid on
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  zlabel ( '<--- Z(X,Y) --->' );
  title ( '2D test function #5, discontinuous medium wave, U(x,t).' );
  view ( -30.0, 45.0 );
  colorbar

  filename = 'edge_test037.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created graphics file "%s".\n', filename );

  return
end


text:
%% When is the probability distribution of the data important?
% The null-hypothesis in the significance tests for WT, XWT and WTC is
% normally distributed AR1 noise. The AR1 coefficient and process variance
% is chosen so that it best fits the observed data. It is therefore quite
% important that the data is close to normal and is reasonably well
% modeled by a Gaussian AR1 process. Otherwise we can trivially reject
% the null-hypothesis and the significance level calculated by the program
% is not appropriate. However, the Central Limit Theorem tells us that the
% distribution tends towards normality as we convolute with longer and longer
% wavelets (_in the absence of long-range persistence_). This means that the
% data distribution is only really important on the shortest scales.
% So, if we are primarily looking at longer scales we do not need to
% worry so much about the distribution. However, for the WT and XWT the
% color of the noise is very important and a very non-normal distribution
% will affect the performance of the ar1 estimators (ar1.m & ar1nv.m).
% The WTC is relatively insensitive to the colour of the noise in the
% significance test (see next question).

author: grinsted

repo: wavelet-coherence


text:
%
%       HDR Toolbox demo 3:
%	   1) Load "CS_Warwick.hdr" HDR image
%	   2) Change the mapping from Longitude Latitude to Cube Map
%	   3) Show the new HDR image using gamma encoding
%	   4) Write the cube map on the disk
%	   5) Write on the disk each face of the cube map ready for OpenGL or
%	   Direct3D
%
%       Author: Francesco Banterle
%       Copyright June 2012 (c)
%
%

clear all;

disp('1) Load "CS_Warwick.hdr" HDR image');
img = hdrimread('demos/CS_Warwick.hdr');

disp('2) Change the mapping from Longitude Latitude to Cube Map');
imgOut = ChangeMapping(img, 'LongitudeLatitude', 'CubeMap');

disp('3) Visualization of the new HDR image using gamma encoding');
ldr = GammaTMO(imgOut, 2.2, 0.0, 1);

disp('4) Write the cube map on the disk');
hdrimwrite(imgOut, 'demos/output/CS_Warwick_CUBE_MAP_format.hdr');

disp('5) Write each face of the cube map on the disk');
CrossCutter(imgOut, 'demos/output/CS_Warwick_CM', 'hdr');

author: banterle


text:
%SKEW Create skew-symmetric matrix
%
% S = SKEW(V) is a skew-symmetric matrix formed from V.
%
% If V (1x1) then S =
%
%           | 0  -v |
%           | v   0 |
%
% and if V (1x3) then S =
%
%           |  0  -vz   vy |
%           | vz    0  -vx |
%           |-vy   vx    0 |
%
%
% Notes::
% - This is the inverse of the function VEX().
% - These are the generator matrices for the Lie algebras so(2) and so(3).
%
% References::
% - Robotics, Vision & Control: Second Edition, Chap 2,
%   P. Corke, Springer 2016.
%
% See also SKEWA, VEX.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function S = skew(v)
    if isvec(v,3)
        % SO(3) case
        S = [  0   -v(3)  v(2)
              v(3)  0    -v(1)
             -v(2) v(1)   0];
    elseif isvec(v,1)
        % SO(2) case
        S = [0 -v; v 0];
    else
        error('SMTB:skew:badarg', 'argument must be a 1- or 3-vector');
    end

author: petercorke

repo: spatialmath-matlab


text:
function imgs = read_neurofinder(foldername,output_name)
%% import neurofinder tiff multifile format

format compact
opts_tiff.append = false;
opts_tiff.big = true;
files = subdir(fullfile(foldername,'*.tif*'));
numFiles = length(files)
data_type = class(imread(files(1).name));
imgs = zeros([size(imread(files(1).name)),numFiles],data_type);
for i = 1:length(files)
    fname = files(i).name;
	imgs(:,:,i) = imread(fname);
    if mod(i,100)==0
        disp(i)
    end
end
if ~isempty(output_name)
    disp('saving to tif')
    saveastiff(cast(imgs,data_type),output_name,opts_tiff);
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [scan,frame,mn] = motionCompNearestImage(view, scans, frames, meanImage)

%    [scan,frame,mn] = motionCompNearestImage(view, [scans], [frames], [meanImage])
%
% gb 01/22/05
%
% Finds the image closest to the mean image. Returns the scan and the frame
% of this image and also its Mean Square Error with the mean image
% Default : scans = current scan
%           frames = all frames
%           meanImage = first frame of the first scan

% Initilizes arguments and variables
global dataTYPES

if ieNotDefined('scans')
    scans = viewGet(view,'curScan');
end
if ieNotDefined('frames')
    frames = 1:length(dataTYPES(curType).scanParams(scans(1)).nFrames);
end
if ieNotDefined('meanImage')
    meanImage = motionCompLoadImages(view,scans(1),1);
end

% In order to save memory, it is necessary to search for the image scan after scan.
%       - If the input 'scans' is a single value, it returns the nearest image
%         for this scan.
%       - If the input 'scans' is an array, it calls it recursively scan
%         after scan and finds the nearest image among the different nearest images

if length(scans) == 1

    % If 'scans' is a single value
    % Loads the whole volume
    scan = scans;
    tSeriesAllSlices = motionCompLoadImages(view, scan);

    % Computes the Mean Square Error of the difference between the images
    % and the mean image
    MSE = motionCompMSE(tSeriesAllSlices(frames,:,:), meanImage);

    % Finds the minimum of this mean square error
    [mn,frame] = min(MSE);

    return
else

    % If 'scans' is an array
    % Initializes the variable
    results = zeros(length(scans),3);

    % Calls itself recursively scan after scan
    for scanIndex = 1:length(scans)
        scanNum = scans(scanIndex);

        fprintf('Searching nearest image in scan %d... ',scanNum);
        [results(scanIndex,1),results(scanIndex,2),results(scanIndex,3)] = ...
            motionCompNearestImage(view, scanNum, frames, meanImage);
        fprintf('Done\n');
    end

    % Finds the closest image among the images found in each scan
    [mn, baseScanIndex] = min(results(:,3));
    scan = results(baseScanIndex,1);
    frame = results(baseScanIndex,2);
    return
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% SP_EXTERIOR_DERIVAITVE: computes the exterior derivative as a matrix with size
%  given by the dimension of two consecutive spaces in the De Rham sequence.
%
%   diff_op = sp_exterior_derivative (space1, space2);
%
% INPUT:
%
%   space1:  domain space of the exterior derivative (number of columns)
%   space2:  image space of the exterior derivative (number of rows)
%
% OUTPUT:
%
%   diff_op: sparse matrix representation of the differential operator.
%
% Copyright (C) 2020-2023 Bernard Kapidani, Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function diff_op = sp_exterior_derivative (space1, space2)

  assert (~strcmpi(space1.transform, 'integral-preserving'), ...
    'The first space cannot be the one for integral-preserving splines (or n-forms)')

  ndim = numel (space1.scalar_spaces{1}.knots);
  if (ndim == 1)
    error ('Not implemented. For dimension 1, use the scalar spaces')
  elseif (ndim == 2)
    if (strcmpi (space1.transform, 'curl-preserving'))
      grad_curl = 'grad';
      output_der = 2;
      deg_shift = {[0 1], [1 0]};
      knt_shift = {[0 2], [2 0]};
      degree = space1.scalar_spaces{1}.degree + [1 0];
      knots = {space1.scalar_spaces{2}.knots{1}, space1.scalar_spaces{1}.knots{2}};
    elseif (strcmpi (space1.transform, 'div-preserving'))
      grad_curl = 'curl';
      output_der = 2;
      deg_shift = {[1 0], [0 1]};
      knt_shift = {[2 0], [0 2]};
      degree = space1.scalar_spaces{1}.degree + [0 1];
      knots = {space1.scalar_spaces{1}.knots{1}, space1.scalar_spaces{2}.knots{2}};
    else
      error ('The second space should be either curl-preserving or div-preserving')
    end
    assert (isa(space2, 'sp_scalar'), 'The two spaces are not compatible')
    for idim = 1:ndim
      assert (all(space1.scalar_spaces{idim}.degree == space2.degree+deg_shift{idim}), 'The degrees are not compatible')
      assert (all(cellfun(@numel,space1.scalar_spaces{idim}.knots) == (cellfun(@numel, space2.knots)+knt_shift{idim})), ...
        'The knot vectors are not compatible')
    end
  elseif (ndim == 3)
    grad_curl = 'grad';
    if (strcmpi (space1.transform, 'curl-preserving'))
      output_der = 2;
      deg_shift = {[-1 1 1], [1 -1 1], [1 1 -1]};
      knt_shift = {[-2 2 2], [2 -2 2], [2 2 -2]};
      degree = space1.scalar_spaces{1}.degree + [1 0 0];
      knots = {space1.scalar_spaces{2}.knots{1}, space1.scalar_spaces{1}.knots{2}, space1.scalar_spaces{1}.knots{3}};
      assert (isa(space2, 'sp_vector'), 'The two spaces are not compatible')
      for idim = 1:ndim
        assert (all(space1.scalar_spaces{idim}.degree == space2.scalar_spaces{idim}.degree+deg_shift{idim}), 'The degrees are not compatible')
        assert (all(cellfun(@numel,space1.scalar_spaces{idim}.knots) == (cellfun(@numel, space2.scalar_spaces{idim}.knots)+knt_shift{idim})), ...
          'The knot vectors are not compatible')
      end
    elseif (strcmpi (space1.transform, 'div-preserving'))
      output_der = 3;
      deg_shift = {[1 0 0], [0 1 0], [0 0 1]};
      knt_shift = {[2 0 0], [0 2 0], [0 0 2]};
      degree = space1.scalar_spaces{1}.degree + [0 1 1];
      knots = {space1.scalar_spaces{1}.knots{1}, space1.scalar_spaces{2}.knots{2}, space1.scalar_spaces{3}.knots{3}};
      assert (isa(space2, 'sp_scalar'), 'The two spaces are not compatible')
      for idim = 1:ndim
        assert (all(space1.scalar_spaces{idim}.degree == space2.degree+deg_shift{idim}), 'The degrees are not compatible')
        assert (all(cellfun(@numel,space1.scalar_spaces{idim}.knots) == (cellfun(@numel, space2.knots)+knt_shift{idim})), ...
          'The knot vectors are not compatible')
      end
    else
      error ('Only implemented for curl-preserving and div-preserving transforms')
    end
  end

  diff_ops = op_geom_exterior (knots, degree, grad_curl);
  diff_op = diff_ops{output_der};
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function test_suite = test_dim_rename()
% tests for cosmo_dim_rename
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite

function test_dim_rename_basics
    ds=cosmo_synthetic_dataset();

    % simple rename
    ds2=cosmo_dim_rename(ds,'j','jj');
    assertEqual(ds2.a.fdim.labels{2},'jj');
    assertEqual(ds2.a.fdim.values,ds.a.fdim.values);
    assertEqual(ds2.samples,ds.samples);
    assertEqual(ds2.sa,ds.sa);

    fa=ds.fa;
    fa.jj=fa.j;
    fa=rmfield(fa,'j');
    assertEqual(ds2.fa,fa);

    % if raise=false, input is same as output
    ds3=cosmo_dim_rename(ds,'xxx','jj',false);
    assertEqual(ds,ds3,ds);

    aet=@(varargin)assertExceptionThrown(@()...
                        cosmo_dim_rename(varargin{:}),'');
    aet(ds,'xxx','jj');
    aet(ds,'xxx','jj',true);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function str = gencode_substruct(subs, name)

% GENCODE_SUBSTRUCT  String representation of subscript structure.
% Generate MATLAB code equivalent to subscript structure subs. See help
% on SUBSTRUCT, SUBSASGN and SUBSREF for details how subscript structures
% are used.
%
% str = gencode_substruct(subs, name)
% Input arguments:
%  subs - a subscript structure
%  name - optional: name of variable to be dereferenced
% Output arguments:
%  str  - a one-line cellstr containing a string representation of the
%         subscript structure
% If name is given, it is prepended to the string.
% For '()' and '{}' also pseudo subscripts are allowed: if subs.subs{...}
% is a string, it will be printed literally, even if it is not equal to
% ':'. This way, it is possible create code snippets that contain
% e.g. references to a loop variable by name.
%
% See also GENCODE, GENCODE_RVALUE, GENCODE_SUBSTRUCTCODE.
%
% This code has been developed as part of a batch job configuration
% system for MATLAB. See
%      http://sourceforge.net/projects/matlabbatch
% for details about the original project.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: gencode_substruct.m 3355 2009-09-04 09:37:35Z volkmar $

rev = '$Rev: 3355 $'; %#ok

ind = 1;
if nargin < 2
    name = '';
end

if ~isstruct(subs) || ~all(isfield(subs, {'type','subs'}))
    if any(exist('cfg_message') == 2:6)
        cfg_message('matlabbatch:usage', 'Item is not a substruct.');
    else
        warning('gencode_substruct:usage', 'Item is not a substruct.');
    end
else
    str = {name};
    for k = 1:numel(subs)
        switch subs(k).type
            case '.',
                str{1} = sprintf('%s.%s', str{1}, subs(k).subs);
            case {'()','{}'},
                str{1} = sprintf('%s%s', str{1}, subs(k).type(1));
                for l = 1:numel(subs(k).subs)
                    if ischar(subs(k).subs{l})
                        substr = subs(k).subs{l};
                    else
                        substr = sprintf('%d ', subs(k).subs{l});
                        if numel(subs(k).subs{l}) > 1
                            substr = sprintf('[%s]', substr(1:end-1));
                        else
                            substr = substr(1:end-1);
                        end
                    end
                    str{1} = sprintf('%s%s, ', str{1}, substr);
                end
                str{1} = str{1}(1:end-2);
                str{1} = sprintf('%s%s', str{1}, subs(k).type(2));
        end
    end
end

author: spm