text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   ABB IRB1200.
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   09/01/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function robot = parameters()

robot.name= 'ABB_IRB1200';

robot.DH.theta= '[q(1) q(2)-pi/2 q(3) q(4) q(5) q(6)+pi]';
robot.DH.d='[0.399 0 0 0.351 0 0.082]';
robot.DH.a='[0 0.35 0.042 0 0 0]';
robot.DH.alpha= '[-pi/2 0 -pi/2 pi/2 -pi/2 0]';
robot.J=[];

robot.inversekinematic_fn = 'inversekinematic_irb1200(robot, T)';
robot.directkinematic_fn = 'directkinematic(robot, q)';

%number of degrees of freedom
robot.DOF = 6;

%rotational: 0, translational: 1
robot.kind=['R' 'R' 'R' 'R' 'R' 'R'];

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-170) deg2rad(170) ; %Axis 1, minimum, maximum
                deg2rad(-100) deg2rad(135); %Axis 2, minimum, maximum
                deg2rad(-200) deg2rad(70); %Axis 3
                deg2rad(-270) deg2rad(270); %Axis 4: Unlimited (400� default)
                deg2rad(-130) deg2rad(130); %Axis 5
                deg2rad(-360) deg2rad(360)]; %Axis 6: Really Unlimited to (800� default)

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(288); %Axis 1, rad/s
                deg2rad(240); %Axis 2, rad/s
                deg2rad(300); %Axis 3, rad/s
                deg2rad(400); %Axis 4, rad/s
                deg2rad(405); %Axis 5, rad/s
                deg2rad(600)];%Axis 6, rad/s

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time

% end effectors maximum velocity
robot.linear_velmax = 7.3; %m/s

%base reference system
%robot.T0 = eye(4);
robot.T0 = [1 0 0 0;
             0 1 0 0;
             0 0 1 0;
             0 0 0 1];

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

% GRAPHICS
robot.graphical.has_graphics=1;
robot.graphical.color = [255 102 51]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-0.5 0.75 -0.75 0.75 0 2.25];
%read graphics files
robot = read_graphics(robot);

robot.has_dynamics=0;


text:
function I = getPcontrolCurrentPack(t,t0,tf,x,param,extra)
%	getPcontrolCurrentPack returns the value of the input current as a function of
%	the time, states and parameters
%
%       I = getPcontrolCurrentPack(t,t0,tf,y,param,extra)
%
%       Inputs:
%               - t         : value of the current time step
%               - t0        : initial integration time
%               - tf        : final integration time
%               - x         : contains the array of all the states
%                             (differential and algebraic) at time t
%               - param     : contains the parameters structure
%               - extra     : extra parameters
%       Outputs:
%               - I     : Applied current desnity [A/m^2]

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% This script provides the value of the applied current density as a
% function of the voltage across the battery pack.

% Firstly we extract the variables related to the first and second cell. To
% this aim the field x_indices of the param structure contains the values
% of the absolute positionong of the variables in the overall x array.
% Indeed, when running simulations with multiple cells, the states array
% (x) will contain as many rows - for a given time instant t - as the sum
% of the states of all the cells involved in the simulation. x_indices
% stores, for each cell, what are the indices of the overall x in which the
% variables of a given cell are stored.

% Extract the first cell variables
cell1_variables = x(param{1}.x_index);
% Extract the second cell variables
cell2_variables = x(param{2}.x_index);

% At this stage, since we have put the variables of the 2 cells into
% cell1_variables and cell2_variables, use their relative indices to
% extract the exact values.

% Cell 1 voltage
V_1 = cell1_variables(param{1}.Phis_indices(1))-cell1_variables(param{1}.Phis_indices(end));

% Cell 2 voltage
V_2 = cell2_variables(param{2}.Phis_indices(1))-cell2_variables(param{2}.Phis_indices(end));

% Define the proportional action. Do not put extreme values, the simulator
% could crash
Kp = 100;

% Define the Voltage Setpoint
V_ref = 8.2;
% Define your linear or nonlinear function of t for evaluate the value
% of I.
I = Kp*(V_ref-(V_1+V_2));

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function [estimate] = ft_inverse_music(cov, leadfield, varargin)

% FT_INVERSE_MUSIC
%
% Use as
%  estimate = ft_inverse_music(cov, leadfield, ...)

author: PatternRecognition

repo: OpenBMI


text:
%% triPatchSmoothRayTraced
% Below is a demonstration of the features of the |triPatchSmoothRayTraced| function

%%
clear; close all; clc;

%% Syntax
% |[VV]=triPatchSmoothRayTraced(F,V,optionStruct);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% op_concatAverages.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% out=op_concatAverages(in1,in2);
%
% DESCRIPTION:
% Concatenate two scans along the averages dimension.  Two scans with 50
% averages each will now look like a single scan with 100 averages.  For
% looping, it is possible to enter a blank structure as the first input.
% In this case, the function will simply return the second input.
%
% INPUTS:
% in1    = first input in matlab structure format.
% in2    = second input in matlab structure format.
%
% OUTPUTS:
% out    = Output following concatenation of inputs along the
%           averages dimension.

function out=op_concatAverages(in1,in2);

if isempty(fieldnames(in1)) %Check if the first argument to be a null string
    out=in2; %If so, return in2;
else
    if in1.dims.averages~=0 && in2.dims.averages~=0 && in1.dims.averages ~= in2.dims.averages
        %RULE:  Averages dimension must be the same for both inputs UNLESS one of
        %the inputs has an averages dimension of zero:
        error('averages dimensions must be the same for both inputs');
    end

    if in1.dims.averages==0 && in2.dims.averages==0
        dim=2;
    elseif in1.dims.averages==0
        dim=in2.dims.averages;
    elseif in2.dims.averages==0;
        dim=in1.dims.averages;
    else
        dim=in1.dims.averages;
    end

    fids=cat(dim,in1.fids,in2.fids);
    specs=cat(dim,in1.specs,in2.specs);
    sz=size(fids);

    %FILLING IN DATA STRUCTURE
    out=in1;
    out.fids=fids;
    out.specs=specs;
    out.sz=sz;
    out.averages=in1.averages+in2.averages;
    out.rawAverages=in1.rawAverages+in2.rawAverages;
    if in1.dims.averages==0
        out.dims.averages=dim;
    end

    %FILLING IN THE FLAGS
    out.flags=in1.flags;
    out.flags.writtentostruct=1;
    out.flags.averaged=0;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
close all
clear
clc

%Select the degraded image
[filename, pathname] = uigetfile({'*.*';'*.bmp';'*.jpg'}, 'Select the Fingerprint image file');

%or_im, read the original image to work on
or_im = imread([pathname,filename]);

[fin_im] =  main_function(or_im);

%imshow(fin_im);
imagesc(fin_im);colormap(gray(256));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [PB] = GB2PB(GB)
% Convert computery things from gigabytes to petabytes.
% Chad A. Greene 2012
PB = GB*2^-20 ;

author: Sable

repo: mcbench-benchmarks


text:
function [b, iCv, iC] = gp_looprep(gp, x, y)
%GP_LOOPREP  Calculates help terms needed in Gaussian GP LOO
%
%  Description
%    [b, iCv] = gp_looprep(gp, x, y, varargin) takes GP structure,
%    training input and output and returns following terms
%      b = C\y
%      iCv = diag(inv(C))
%    [b, iCv, iC] = gp_looprep(gp, x, y, varargin) returns additionally
%      iC = inv(C)
%    In case of CS, FIC, PIC and CS+FIC, these terms are computed
%    using appropriate efficient sparse matrix and inverse lemma
%    computations.
%

% Copyright (c) 2012 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  switch gp.type
    case 'FULL'
      % FULL GP (and compact support GP)
      [K, C] = gp_trcov(gp, x);

      if issparse(C)
        if nargout==3
          iC = full(spinv(C));        % sparse inverse
          iCv = full(diag(iC));       % the diagonal of sparse inverse
        else
          iCv = full(diag(spinv(C))); % the diagonal of sparse inverse
        end
        [LD, notpositivedefinite] = ldlchol(C);
        if notpositivedefinite
          b = NaN;
          iCv = NaN;
          iC = NaN;
          return
        end
        b = ldlsolve(LD,y);
      else
        if nargout==3
          iC = inv(C);                % full inverse
          iCv = diag(inv(C));         % the diagonal of full inverse
        else
          iCv = diag(inv(C));         % the diagonal of full inverse
        end
        b = C\y;
      end

    case 'FIC'
      % FIC
      % Use inverse lemma for FIC low rank covariance matrix approximation
      % Code adapated from gp_pred

      u = gp.X_u;
      m = size(u,1);
      % Turn the inducing vector on right direction
      if size(u,2) ~= size(x,2)
        u=u';
      end
      [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
      K_fu = gp_cov(gp, x, u);   % f x u
      K_uu = gp_trcov(gp, u);     % u x u, noiseles covariance K_uu
      [Luu, notpositivedefinite] = chol(K_uu,'lower');
      if notpositivedefinite
        b = NaN;
        iCv = NaN;
        iC = NaN;
        return
      end
      B=Luu\(K_fu');
      Qv_ff=sum(B.^2)';
      Lav = Cv_ff-Qv_ff;   % 1 x f, Vector of diagonal elements

      % iLaKfu = diag(inv(Lav))*K_fu = inv(La)*K_fu
      iLaKfu = zeros(size(K_fu));  % f x u,
      n=size(x,1);
      for i=1:n
        iLaKfu(i,:) = K_fu(i,:)./Lav(i);  % f x u
      end
      A = K_uu+K_fu'*iLaKfu;
      A = (A+A')./2;
      L = iLaKfu/chol(A);

      if nargout==3
        iC = diag(1./Lav) - L*L'; % Cv = inv(C);
        iCv = diag(iC); %           iCv = diag(inv(C));
      else
        iCv = 1./Lav - sum(L.^2,2); % iCv = diag(inv(C));
      end
      b = y./Lav - L*(L'*y);      % b = C\y;

    case {'PIC' 'PIC_BLOCK'}
      % PIC
      % Use inverse lemma for PIC low rank covariance matrix approximation
      % Code adapated from gp_pred

      u = gp.X_u;
      ind = gp.tr_index;
      if size(u,2) ~= size(x,2)
        % Turn the inducing vector on right direction
        u=u';
      end

      % Calculate some help matrices
      [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
      K_fu = gp_cov(gp, x, u);         % f x u
      K_uu = gp_trcov(gp, u);    % u x u, noiseles covariance K_uu
      Luu = chol(K_uu)';

      % Evaluate the Lambda (La) for specific model
      % Q_ff = K_fu*inv(K_uu)*K_fu'
      % Here we need only the diag(Q_ff), which is evaluated below
      B=Luu\K_fu';
      iLaKfu = zeros(size(K_fu));  % f x u
      for i=1:length(ind)
        Qbl_ff = B(:,ind{i})'*B(:,ind{i});
        [Kbl_ff, Cbl_ff] = gp_trcov(gp, x(ind{i},:));
        La{i} = Cbl_ff - Qbl_ff;
        iLaKfu(ind{i},:) = La{i}\K_fu(ind{i},:);
      end
      A = K_uu+K_fu'*iLaKfu;
      A = (A+A')./2;            % Ensure symmetry
      L = iLaKfu/chol(A);

      % From this on evaluate the prediction
      % See Snelson and Ghahramani (2007) for details
      n=size(y,1);
      iCv=zeros(n,1);
      if nargout==3
        iC=zeros(n,n);
      end
      for i=1:length(ind)
        if nargout==3
          iC(ind{i},ind{i}) = inv(La{i});
        else
          iCv(ind{i},:) = diag(inv(La{i}));
        end
        b(ind{i},:) = La{i}\y(ind{i},:);
      end
      if nargout==3
        iC = iC - L*L';          % iC  = inv(C);
        iCv = diag(iC);          % iCv = diag(inv(C));
      else
        iCv = iCv - sum(L.^2,2); % iCv = diag(inv(C));
      end
      b = b - L*(L'*y);        % b = C\y;

    case 'CS+FIC'
      % CS+FIC
      % Use inverse lemma for CS+FIC
      % Code adapated from gp_pred

      u = gp.X_u;
      if size(u,2) ~= size(x,2)
        % Turn the inducing vector on right direction
        u=u';
      end

      n = size(x,1);
      m = size(u,1);
      ncf = length(gp.cf);

      % Indexes to all non-compact support and compact support covariances.
      cf1 = [];
      cf2 = [];

      % Loop through all covariance functions
      for i1 = 1:ncf
        if ~isfield(gp.cf{i1},'cs')
          % Non-CS covariances
          cf1 = [cf1 i1];
        else
          % CS-covariances
          cf2 = [cf2 i1];
        end
      end

      % First evaluate needed covariance matrices
      % v defines that parameter is a vector
      [Kv_ff, Cv_ff] = gp_trvar(gp, x, cf1); % f x 1  vector
      K_fu = gp_cov(gp, x, u, cf1);          % f x u
      K_uu = gp_trcov(gp, u, cf1);    % u x u, noiseles covariance K_uu
      K_uu = (K_uu+K_uu')./2;         % ensure the symmetry of K_uu

      Luu  = chol(K_uu)';
      % Evaluate the Lambda (La)
      % Q_ff = K_fu*inv(K_uu)*K_fu'
      B=Luu\(K_fu');       % u x f
      Qv_ff=sum(B.^2)';
      Lav = Cv_ff-Qv_ff;   % f x 1, Vector of diagonal elements

      K_cs = gp_trcov(gp,x,cf2);
      La = sparse(1:n,1:n,Lav,n,n) + K_cs;

      iLaKfu = La\K_fu;
      A = K_uu+K_fu'*iLaKfu;
      A = (A+A')./2;     % Ensure symmetry
      L = iLaKfu/chol(A);

      if nargout==3
        iC = inv(La) - L*L'; % iCv = diag(inv(C));
        iCv = diag(iC)     ; % iCv = diag(inv(C));
      else
        iCv = diag(inv(La)) - sum(L.^2,2); % iCv = diag(inv(C));
      end
      b = La\y - L*(L'*y);         % b = C\y;

    case 'SSGP' % SSGP
      error('GP_LOOPRED is not yet implemented for SSGP and Gaussian likelihood')

    otherwise
      error('Unknown type of Gaussian process')
  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [cMap]=viridis(varargin)

% function [cMap]=viridis(n)
% ------------------------------------------------------------------------
% Creates the colormap data for the viridis colormap.
%
% This colormap is a MATLAB implementation of the matplotlib colormap by
% Nathaniel J. Smith, Stefan van der Walt, and Eric Firing.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2018/04/10
%------------------------------------------------------------------------

%%
switch nargin
    case 0
        n=[];
    case 1
        n=varargin{1};
end

%%

[cMap]=matplotlibColormap('viridis',n);

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


text:
function fuse = spm_cfg_eeg_fuse
% configuration file for fusing M/EEG files
%_______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_cfg_eeg_fuse.m 5377 2013-04-02 17:07:57Z vladimir $

rev = '$Rev: 5377 $';
D = cfg_files;
D.tag = 'D';
D.name = 'File Names';
D.filter = 'mat';
D.num = [2 Inf];
D.help = {'Select the M/EEG mat files.'};

prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the fused dataset. Default prefix is ''u''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.val     = {'u'};

fuse = cfg_exbranch;
fuse.tag = 'fuse';
fuse.name = 'Fusion';
fuse.val = {D, prefix};
fuse.help = {'Fuse EEG/MEG data.'};
fuse.prog = @eeg_fuse;
fuse.vout = @vout_eeg_fuse;
fuse.modality = {'EEG'};

function out = eeg_fuse(job)
% construct the S struct
S = job;
S.D = char(S.D);

out.D = spm_eeg_fuse(S);
out.Dfname = {fullfile(out.D.path, out.D.fname)};

function dep = vout_eeg_fuse(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'Fused Data';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2) = cfg_dep;
dep(2).sname = 'Fused Datafile';
% reference field "Dfname" from output
dep(2).src_output = substruct('.','Dfname');
% this can be entered into any file selector
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


text:
function a = toeplitz_pds ( m, n, x, y )

%*****************************************************************************80
%
%% TOEPLITZ_PDS returns the TOEPLITZ_PDS matrix.
%
%  Discussion:
%
%    TOEPLITZ is a Toeplitz matrix that is positive definite symmetric.
%
%  Formula:
%
%    A(I,J) = sum ( 1 <= K <= M ) Y(K) * cos ( 2 * PI * X(K) * (I-J) )
%
%  Example:
%
%    N = 5, M = 5,
%    X = ( -0.0625, - 0.03125, 0.0, 0.03125, 0.0625 ),
%    Y = ( 0.2, 0.2, 0.2, 0.2, 0.2)
%
%    1.000000  0.961866  0.852395  0.685661  0.482843
%    0.961866  1.000000  0.961866  0.852395  0.685661
%    0.852395  0.961866  1.000000  0.961866  0.852395
%    0.685661  0.852395  0.961866  1.000000  0.961866
%    0.482843  0.685661  0.852395  0.961866  1.000000
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A is positive definite or positive semi-definite, depending on
%    the values of X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Cybenko, Charles Van Loan,
%    Computing the minimum eigenvalue of a symmetric positive definite
%    Toeplitz matrix,
%    SIAM Journal on Scientific and Statistical Computing,
%    Volume 7, 1986, pages 123-131.
%
%  Parameters:
%
%    Input, integer M, the number of terms of W and X.
%
%    Input, integer N, the order of A.
%
%    Input, real X(M), used to define the matrix.
%
%    Input, real Y(M), a set of positive weights used to define the matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n
      a(i,j) = 0.0;
      for k = 1 : m
        angle = 2.0 * pi * x(k) * ( i - j );
        a(i,j) = a(i,j) + y(k) * cos ( angle );
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function hs = binaryTreeApply( X, tree, maxDepth, minWeight, nThreads )
% Apply learned binary decision tree classifier.
%
% USAGE
%  hs = binaryTreeApply( X, tree, [maxDepth], [minWeight], [nThreads] )
%
% INPUTS
%  X          - [NxF] N length F feature vectors
%  tree       - learned tree classification model
%  maxDepth   - [] maximum depth of tree
%  minWeight  - [] minimum sample weigth to allow split
%  nThreads   - [16] max number of computational threads to use
%
% OUTPUTS
%  hs         - [Nx1] predicted output log ratios
%
% EXAMPLE
%
% See also binaryTreeTrain
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.40
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if(nargin<3 || isempty(maxDepth)), maxDepth=0; end
if(nargin<4 || isempty(minWeight)), minWeight=0; end
if(nargin<5 || isempty(nThreads)), nThreads=16; end
if(maxDepth>0), tree.child(tree.depth>=maxDepth) = 0; end
if(minWeight>0), tree.child(tree.weights<=minWeight) = 0; end
hs = tree.hs(forestInds(X,tree.thrs,tree.fids,tree.child,nThreads));

end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function V = mttkrp(X,U,n)
%MTTKRP Matricized tensor times Khatri-Rao product for ttensor.
%
%   V = MTTKRP(X,U,n) efficiently calculates the matrix product of the
%   n-mode matricization of X with the Khatri-Rao product of all
%   entries in U, a cell array of matrices, except the nth.  How to
%   most efficiently do this computation depends on the type of tensor
%   involved.
%
%   See also TTENSOR, TTENSOR/TTV
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

N = ndims(X);

if (n==1)
    R = size(U{2},2);
else
    R = size(U{1},2);
end

% Compute cell array of weights to multiply into core
W = cell(N,1);
for i = [1:n-1,n+1:N]
  W{i} = (X.u{i}' * U{i});
end
Y = mttkrp(X.core,W,n);

% Find each column of answer by multiplying columns of X.u{n} with weights
V = X.u{n} * Y;


text:
function [sigma,dkm,g,r]=noiseest(dlogCr,dlogr,method)
%Syntax: [sigma,dkm,g,r]=noiseest(dlogCr,dlogr,method)
%_____________________________________________________
%
% Calculates the noise standard deviation from the derivative of the
% Correlation Integral. Requires the auxilary function "dkmminusg".
%
% sigma is the noise standard deviation.
% dkm is the empirical effect of noise on the Correlation Integral.
% g is the theoritical effect of noise on the Correlation Integral.
% r is the range.
% dlogCr is the derivative of the logCr.
% dlogr is the log(range).
% method can take one of the folloing values:
%  'full' for Schreiber's full minimization
%  'mod' for Leontitsis et al. modified minimization
%
%
% References:
%
% Schreiber T (1993): Determination of the noise level of chaotic time
% series. Physical Review E 48(1): R13-R16
%
% Leontitsis A, Pange J., Bountis T. (2003): Large noise level estimation.
% International Journal of Bifurcation and Chaos 13(8): 2309-2313
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% June 15, 2001.

% dlogCr and dlogr must have the same number of rows
if size(dlogCr,1)~=size(dlogr,1)
   error('dlogCr and dlogr must have the same number of rows.');
end

r=10.^dlogr;
options=optimset('Display','off');
for i=2:size(dlogCr,2)
   dkm(:,i-1)=(dlogCr(:,i)-dlogCr(:,1))/(i-1);
   sigma(i-1)=fminbnd(@dkmminusg,r(1),r(end),options,dkm(:,i-1),r,method);
end
z=r./2./sigma(1);
g=2.*z.*exp(-z.^2)./sqrt(pi)./erf(z);

function error=dkmminusg(s,dkm,r,method)
%Syntax: error=dkmminusg(s,dkm,r,method)
%_______________________________________
%
% Auxilary function for "noiseest". Calculates the  error between the dkm
% and function "g" given a noise standard deviation (s) and a vector of
% range values (r).
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% June 15, 2001.

% Calculate the auxilary variable z
z=r./2./s;

% Calpculate the g fumction
g=2*z.*exp(-z.^2)/sqrt(pi)./erf(z);

switch method
    case 'full'
        % The ordinary (low noise) calculation
        error=norm(dkm-g);
    case 'mod'
        % The modified (large noise) estimation
        if any(dkm<g)==1
            i=find(dkm<g);
            error=sum(abs(dkm(i)-g(i)));
        else
            error=norm(dkm-g,-inf);
        end
   otherwise
       error('You should provide another value for method.');
end


text:
function [mask] = polygon_filter(x, y, XI, YI, in_or_out)
    % POLYGON_FILTER used Analytic Geometry to select points inside or outside polygon
    %
    % [mask] = POLYGON_FILTER(x, y, XI, YI, in_or_out)
    %
    % input params:
    %    x, y : polygon vertices, with x(1)==x(end) and y(1)==y(end)
    %           x & y must be vectors of same size, and include
    %           more than 3 points (len(polyX) > 3)
    %
    %    XI, YI : x and y values for points to be tested
    %
    %    in_or_out : either 'inside' or 'outside'
    %               'inside' : returns a mask that is true for points inside polygon
    %               'outside': returns mask that is true for points outside polygon
    %
    %               default is 'inside'
    %
    %    I don't know about edge cases.
    %
    % see also inpolygon, inpoly

    assert(length(XI) == length(YI), 'test-points should have equal numbers of x and y' );
    assert(length(x) == length(y), 'number of polygon x should equal polygon y');
    assert(numel(x) > 3,'not enough points to be a polygon');

    m = length(x)-1;      %  number of coordinates of polygon
    l = zeros(size(XI));
    l2 = l;            %  Algorithm to select points inside a closed
    %  polygon based on Analytic Geometry    R.Z. 4/94
    for i = 1:m

        l= ((y(i)-YI < 0) & (y(i+1)-YI >= 0)) & ...
            (XI-x(i)-(YI-y(i))*(x(i+1)-x(i))/(y(i+1)-y(i)) < 0) | ...
            ((y(i)-YI >= 0) & (y(i+1)-YI < 0)) & ...
            (XI-x(i)-(YI-y(i))*(x(i+1)-x(i))/(y(i+1)-y(i)) < 0);

        if i ~= 1
            l2(l) = 1 - l2(l);
        else
            l2 = l;
        end         % if i

    end

    if ~exist('in_or_out','var')
        in_or_out = 'inside';
    end

    switch(in_or_out)
        case 'inside'
            mask = l2;
        case 'outside'
            mask = ~l2;
        otherwise
            error("unrecognized in_or_out option. either 'inside' or 'outside'");
    end
end

author: CelsoReyes


text:
function [x_star_all, varx_star_all, mini] = ...
    vargplvmPredictLatent(model, Yts, testInd, batchMode, iters, displayTestOpt, mini, initQx, parallel)

% VARGPLVMPREDICTLATENT Given output observation(s) compute the approximate
% posterior
% DESC Given an output observed test point y* (or set of test points Y*)
% compute the approximate posterior p(X* | Y*) for the corresponding latent
% points.
%
% ARG model: The trained vargplvm model
% ARG Yts: The set of outputs observations
% ARG testInd: The indices of Yts to do inference for (default: all)
% ARG batchMode: If set to true, then if N* points are given in Yts, only
% one optimisation will be done in the N*-dimensional space. Otherwise, N*
% 1-dimensional optimisations will be performed. Batch mode is faster but
% possibly less accurate.
% ARG iters: the number of test iterations for the optimisation(s)
% ARG displayTestOpt: If set to true, then the test iterations will be
% displayed
% ARG mini: The indices of the points in the test set for which the test
% latent points are initialised with, ie x*(i,:) will be initialised with
% Xtrain(mini(i), :). If this vector is not given then this vector is
% computed here in the default way.
% ARG initQx: Direct initialisation (instead of mini). Has to contain
% .means and .covars
%
% SEEALSO:
%
% COPYRIGHT : Andreas C. Damianou, 2013
%
% VARGPLVM

if nargin < 2 || isempty(model) || isempty(Yts), error('At least 2 arguments needed!'); end
if nargin < 3 || isempty(testInd), testInd = 1:size(Yts,1); end
if nargin < 4 || isempty(batchMode), batchMode = false; end
if nargin < 5 || isempty(iters), iters = model.globalOpt.reconstrIters; end
if nargin < 6 || isempty(displayTestOpt), displayTestOpt = false; end
if nargin < 7, mini = []; end
if nargin < 8, initQx = []; end
if nargin < 9, parallel = []; end

try
    pool_open = matlabpool('size')>0;
catch e
    pool_open = 0;
end
if pool_open && ~isempty(parallel) && size(Yts,1) > 3 && ~batchMode
    parallel = parallel;
    warning('Parallel predictions!')
else
    parallel = false;
end

miniGiven = ~isempty(mini);
initQxGiven = ~isempty(initQx);

if initQxGiven && miniGiven
    error('Both initial training indices AND initial X was given!')
end

if batchMode
    y_star = Yts(testInd,:);
    if ~initQxGiven
        if ~miniGiven
            mini = NaN(1,length(testInd));
            for i=1:length(testInd)
                dst = dist2(y_star(i,:), model.y);
                [~, mini(i)] = min(dst);
            end
        end
        model.mini = mini;
        vardistx = vardistCreate(model.vardist.means(mini, :), model.q, 'gaussian');
        vardistx.covars = model.vardist.covars(mini, :);
    else
        vardistx = vardistCreate(initQx.means, model.q, 'gaussian');
        vardistx.covars = initQx.covars;
    end
    model.vardistx = vardistx;
    [x_star_all, varx_star_all] = vargplvmOptimisePoint(model, vardistx, y_star, true, iters);
else
    if parallel
        x_star_allCell = cell(length(testInd), 1);
        varx_star_allCell = cell(length(testInd), 1);
        if ~initQxGiven &&  ~miniGiven
            mini = NaN(1,length(testInd));
        end
        for i=1:length(testInd)
             curInd = testInd(i);
             YYts{i} = Yts(curInd,:);
             if initQxGiven
                initQxMeans{i} = initQx.means(i,:);
                initQxCovars{i} = initQx.covars(i,:);
             else
                 initQxMeans{i} = [];
                 initQxCovars{i} = [];
             end
             miniCell{i} = mini(i);
        end
         Ytr = model.y;
        pb = myProgressBar(length(testInd), min(length(testInd),20));
        for i=1:length(testInd)
            pbb{i} = pb;
        end
        parfor i=1:length(testInd)
            pbb{i} = myProgressBar(pbb{i},i);
            y_star =YYts{i};
            if ~initQxGiven
                if ~miniGiven
                    dst = dist2(y_star, Ytr);
                    [~, miniCur] = min(dst);
                    miniCell{i} = miniCur;
                else
                    miniCur = miniCell{i};
                end
                modelCur = model;
                modelCur.mini = miniCur;
                vardistx = vardistCreate(modelCur.vardist.means(miniCur,:), modelCur.q, 'gaussian');
                vardistx.covars = modelCur.vardist.covars(miniCur,:);
            else
                vardistx = vardistCreate(initQxMeans{i}, modelCur.q, 'gaussian');
                vardistx.covars = initQxCovars{i};
            end
            modelCur.vardistx = vardistx;

            % Find p(X_* | Y_*) which is approximated by q(X_*)
            [x_star, varx_star] = vargplvmOptimisePoint(modelCur, vardistx, y_star, displayTestOpt, iters);
            x_star_allCell{i} = x_star;
            varx_star_allCell{i} = varx_star;
        end
        for i=1:length(testInd)
            x_star_all(i,:) = x_star_allCell{i};
            varx_star_all(i,:) = varx_star_allCell{i};
            mini(i) = miniCell{i};
        end
        fprintf(1, '\n');
%         %.. TODO
%         x_star_all = zeros(length(testInd), model.q);
%         varx_star_all = zeros(length(testInd), model.q);
%         if ~initQxGiven &&  ~miniGiven
%             mini = NaN(1,length(testInd));
%         end
%         pb = myProgressBar(length(testInd), min(length(testInd),20));
%         YYts = cell(length(testInd),1);
%         initQxMeans = cell(length(testInd),1);
%         initQxCovars = cell(length(testInd),1);
%         for i=1:length(testInd)
%             curInd = testInd(i);
%             YYts{i} = Yts(curInd,:);
%             initQxMeans{i} = initQx.means(i,:);
%             initQxCovars{i} = initQx.covars(i,:);
%         end
%         parfor i=1:length(testInd)
%             pb = myProgressBar(pb,i);
%             %curInd = testInd(i);
%             y_star = YYts{i};
%             modelTmp = model;
%             if ~initQxGiven
%                 if ~miniGiven
%                     dst = dist2(y_star, model.y);
%                     [~, mini(i)] = min(dst);
%                 end
%                 modelTmp.mini = mini(i);
%                 vardistx = vardistCreate(modelTmp.vardist.means(mini(i),:), modelTmp.q, 'gaussian');
%                 vardistx.covars = modelTmp.vardist.covars(mini(i),:);
%             else
%                 vardistx = vardistCreate(initQxMeans{i}, modelTmp.q, 'gaussian');
%                 vardistx.covars = initQxCovars{i};
%             end
%             modelTmp.vardistx = vardistx;
%
%             % Find p(X_* | Y_*) which is approximated by q(X_*)
%             [x_star, varx_star] = vargplvmOptimisePoint(modelTmp, vardistx, y_star, displayTestOpt, iters);
%             x_star_all(i,:) = x_star;
%             varx_star_all(i,:) = varx_star;
%         end
%         fprintf(1, '\n');
    else
        x_star_all = zeros(length(testInd), model.q);
        varx_star_all = zeros(length(testInd), model.q);
        if ~initQxGiven &&  ~miniGiven
            mini = NaN(1,length(testInd));
        end
        pb = myProgressBar(length(testInd), min(length(testInd),20));
        for i=1:length(testInd)
            pb = myProgressBar(pb,i);
            curInd = testInd(i);
            y_star = Yts(curInd,:);
            if ~initQxGiven
                if ~miniGiven
                    dst = dist2(y_star, model.y);
                    [~, mini(i)] = min(dst);
                end
                model.mini = mini(i);
                vardistx = vardistCreate(model.vardist.means(mini(i),:), model.q, 'gaussian');
                vardistx.covars = model.vardist.covars(mini(i),:);
            else
                vardistx = vardistCreate(initQx.means(i,:), model.q, 'gaussian');
                vardistx.covars = initQx.covars(i,:);
            end
            model.vardistx = vardistx;

            % Find p(X_* | Y_*) which is approximated by q(X_*)
            [x_star, varx_star] = vargplvmOptimisePoint(model, vardistx, y_star, displayTestOpt, iters);
            x_star_all(i,:) = x_star;
            varx_star_all(i,:) = varx_star;
        end
        fprintf(1, '\n');
    end
end

author: SheffieldML

repo: vargplvm


text:
function init(fname, nbytes, opts)
% Initialise binary file on disk
% FORMAT init(fname, nbytes[, opts])
% fname   - filename
% nbytes  - data size {bytes}
% opts    - optional structure with fields:
%   .offset   - file offset {bytes} [default: 0]
%   .wipe     - overwrite exisiting values with 0 [default: false]
%   .truncate - truncate file if larger than requested size [default: true]
%
% This function is normally called by file_array/initialise
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: init.m 7147 2017-08-03 14:07:01Z spm $

%-This is merely the help file for the compiled routine
error('init.c not compiled - see Makefile');


text:
function blurTSeriesPlot(blurFcn)
% Blurs the tseries data in the current graph window
%
%   blurTSeriesPlot(blurFcn)
%
% The blurring function, by default, sums to one to preserve the mean. It
% is [0.2500    0.5000    0.2500].  For more blurring, repeat this call.
%
% The plot is created either from plotting mean tSeries, or multiple
% tSeries.  This routine also runs on Single Cycle Plots.
%
% Re-written so that it might be extended in the future.  We should
% probably add identifiers to the userData structure in the plot windows.
% (BW)
%
% djh, 3/2001
global GRAPHWIN
if (isempty(GRAPHWIN) || GRAPHWIN==0), myErrorDlg('blurPlot: no graph to blur');
else  set(0,'CurrentFigure',GRAPHWIN);
end
if notDefined('blurFcn'), blurFcn = [1 2 1]; blurFcn = blurFcn/sum(blurFcn); end
axisHandles = get(gcf,'Children');
for h = 1:length(axisHandles)
    subplot(axisHandles(h));
    data = get(gca,'userData');

    if isfield(data,'tSeries')
        typeOfPlot = 'tseries';
        data.time = data.frameNumbers;
        data.tSeries = conv2(data.tSeries,blurFcn,'same');
    elseif isfield(data,'x') && isfield(data,'y')
        typeOfPlot = 'singlecycle';
        data.time = data.x;
        data.tSeries = conv2(data.y,blurFcn,'same');
    else
        myErrorDlg('blurPlot: cannot identify plot type.  Edit blurTseriesPlot')
    end

    % get plot properties
    xLim = get(gca,'xLim');    yLim = get(gca,'yLim');
    xTick = get(gca,'xTick');  yTick = get(gca,'yTick');
    xGrid = get(gca,'xGrid');  yGrid = get(gca,'yGrid');
    fontName = get(get(gca,'xlabel'),'fontName');
    fontSize = get(get(gca,'xLabel'),'fontSize');
    xLabel   = get(get(gca,'xLabel'),'String');
    yLabel   = get(get(gca,'yLabel'),'String');
    titleStr = get(get(gca,'title'),'String');

    color = get(get(gca,'Children'),'Color');
    if iscell(color), color = color{1}; end
    lineWidth = get(get(gca,'Children'),'lineWidth');
    if iscell(lineWidth), lineWidth = lineWidth{1}; end

    switch lower(typeOfPlot)
        case 'tseries'
            p = plot(data.time,data.tSeries);
        case 'singlecycle'
            p = errorbar(data.time,data.tSeries,data.e);
        otherwise
    end

    % Replot the data
    set(p,'lineWidth',lineWidth);
    set(p,'Color',color);

    % reset properties
    set(gca,'xLim',xLim);   set(gca,'yLim',yLim);
    set(gca,'xTick',xTick); set(gca,'yTick',yTick);
    set(gca,'xGrid',xGrid); set(gca,'yGrid',yGrid);
    xlabel(xLabel,'fontName',fontName,'fontSize',fontSize);
    ylabel(yLabel,'fontName',fontName,'fontSize',fontSize);
    title(titleStr,'fontName',fontName,'fontSize',fontSize);
    set(gca,'UserData',data);

end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% minimum distance from head of fly to central of any other fly
function [data,units] = compute_dhead2central(trx,n)

larvae = trx.exp2flies{n};
nlarvae = numel(larvae);
data = cell(1,nlarvae);

for i1 = 1:nlarvae,
  larva1 = larvae(i1);
  % access closestfly to ensure that dhead2central is computed
  trx(larva1).closestlarva_head2central;
  data{i1} = trx(larva1).dhead2central;
end
units = parseunits('mm');

author: kristinbranson

repo: JAABA


text:
function x = meicv2(c,N1,N2,ns,nag)
% meicv2 - 2D inverse mirror-extended curvelet transform
% -----------------
% INPUT
% --
% c is a cell array which contains the curvelets coefficients. If
% tp=='ortho', then c{j}{l}(n1,n2) is the coefficient at scale j,
% direction l and spatial index (n1,n2). The directional index l
% iterates through the wedges in the first quadrant. Notice that, for
% the mirror-extended wave atoms, the spatial indices wrap around once.
% --
% N1, N2 are positive integers.
% --
% ns is the number of levels, including the coarsest level. ns =
% ceil(log2(min(N1,N2)) - 3) is commonly used.
% --
% nag is the number of angles used for the second coarsest level.
% nag is required to be a multiple of 4 and nag = 16 is often used.
% -----------------
% OUTPUT
% --
% x is an N1-by-N2 matrix.
% -----------------
% Written by Lexing Ying and Laurent Demanet, 2007

  E1 = ceil(N1/3);  E2 = ceil(N2/3);  %E1 = 0;  E2 = 0;
  A1 = 2*(N1+E1);  A2 = 2*(N2+E2);
  fd = zeros(A1,A2);

  G1 = 4/3*N1;  G2 = 4/3*N2;

  for s=ns:-1:2
    R1 = 2^(s-ns)*G1;
    R2 = 2^(s-ns)*G2;

    idx1 = [ceil(-R1):floor(R1)];
    [wl,wr] = cvwindow((idx1+R1/1)/(R1/2));    tmpa = wl;
    [wl,wr] = cvwindow((idx1-R1/2)/(R1/2));    tmpb = wr;
    coef1 = tmpa.*tmpb;
    idx2 = [ceil(-R2):floor(R2)];
    [wl,wr] = cvwindow((idx2+R2/1)/(R2/2));    tmpa = wl;
    [wl,wr] = cvwindow((idx2-R2/2)/(R2/2));    tmpb = wr;
    coef2 = tmpa.*tmpb;
    lowpass = coef1'*coef2;

    idx1 = [ceil(-R1):floor(R1)];
    [wl,wr] = cvwindow((idx1+R1/2)/(R1/4));    tmpa = wl;
    [wl,wr] = cvwindow((idx1-R1/4)/(R1/4));    tmpb = wr;
    coef1 = tmpa.*tmpb;
    idx2 = [ceil(-R2):floor(R2)];
    [wl,wr] = cvwindow((idx2+R2/2)/(R2/4));    tmpa = wl;
    [wl,wr] = cvwindow((idx2-R2/4)/(R2/4));    tmpb = wr;
    coef2 = tmpa.*tmpb;
    tmppass = coef1'*coef2;
    hghpass = sqrt(1-tmppass.^2);

    pass = lowpass.*hghpass;
    [M1,M2] = size(pass);

    fh = zeros(M1,M2);
    %get fh

    nbangles = nag*2^(ceil((s-2)/2));
    %---------
    [M1,M2] = size(fh);
    nd = nbangles/4;
    cs = c{s};
    W1 = 2*R1/nd;  W2 = 2*R2/nd;

    %take only first quadrant
    cnt = 1;
    for g=nd/2:nd-1
      xs = R1/4-(W1/2)/4;    xe = R1;
      ys = -R2 + (2*g-1)*W2/2;		ye = -R2 + (2*g+3)*W2/2;
      xn = ceil(xe-xs);    yn = ceil(ye-ys);
      if(g==0)
        thts = atan2(-1.0, 1.0-1.0/nd);
        thtm = atan2(-1.0+1.0/nd, 1.0);
        thte = atan2(-1.0+3.0/nd, 1.0);
      elseif(g==nd-1)
        thts = atan2(-1.0+(2.0*g-1.0)/nd, 1.0);
        thtm = atan2(-1.0+(2.0*g+1.0)/nd, 1.0);
        thte = atan2(1.0, 1.0-1.0/nd);
      else
        thts = atan2(-1.0+(2.0*g-1.0)/nd, 1.0);
        thtm = atan2(-1.0+(2.0*g+1.0)/nd, 1.0);
        thte = atan2(-1.0+(2.0*g+3.0)/nd, 1.0);
      end
      %fprintf(1,'%d %d %d\n',thts,thtm,thte);
      R21 = R2/R1;
      wpdata = fft2(cs{cnt}) / sqrt(numel(cs{cnt}));
      cnt = cnt+1;
      for xcur=ceil(xs):xe
        yfm = ceil( max([-R2, R21*xcur*tan(thts)]) );
        yto = floor( min([R2, R21*xcur*tan(thte)]) );
        ycur = yfm:yto;
        thtcur = atan2(ycur/R2,xcur/R1);
        [al,ar] = cvwindow((thtcur-thts)/(thtm-thts));
        [bl,br] = cvwindow((thtcur-thtm)/(thte-thtm));
        pou = al.*br;
        fh(mod(xcur,M1)+1,mod(ycur,M2)+1) = fh(mod(xcur,M1)+1,mod(ycur,M2)+1) + wpdata(mod(xcur,xn)+1,mod(ycur,yn)+1) .* pou;
      end
    end

    for f=nd-1:-1:nd/2
      ys = R2/4-(W2/2)/4;		  ye = R2;
      xs = -R1 + (2*f-1)*W1/2;		  xe = -R1 + (2*f+3)*W1/2;
      xn = ceil(xe-xs);		  yn = ceil(ye-ys);
      if(f==0)
        phis = atan2(-1.0, 1.0-1.0/nd);
        phim = atan2(-1.0+1.0/nd, 1.0);
        phie = atan2(-1.0+3.0/nd, 1.0);
      elseif(f==nd-1)
        phis = atan2(-1.0+(2.0*f-1.0)/nd, 1.0);
        phim = atan2(-1.0+(2.0*f+1.0)/nd, 1.0);
        phie = atan2(1.0, 1.0-1.0/nd);
      else
        phis = atan2(-1.0+(2.0*f-1.0)/nd, 1.0);
        phim = atan2(-1.0+(2.0*f+1.0)/nd, 1.0);
        phie = atan2(-1.0+(2.0*f+3.0)/nd, 1.0);
      end
      %fprintf(1,'%d %d %d\n',phis,phim,phie);
      R12 = R1/R2;
      wpdata = fft2(cs{cnt}) / sqrt(numel(cs{cnt}));
      cnt = cnt+1;
      for ycur=ceil(ys):ye
        xfm = ceil( max([-R1, R12*ycur*tan(phis)]) );
        xto = floor( min([R1, R12*ycur*tan(phie)]) );
        xcur = xfm:xto;
        phicur = atan2(xcur/R1, ycur/R2);
        [al,ar] = cvwindow((phicur-phis)/(phim-phis));
        [bl,br] = cvwindow((phicur-phim)/(phie-phim));
        pou = al.*br;
        fh(mod(xcur,M1)+1,mod(ycur,M2)+1) = fh(mod(xcur,M1)+1,mod(ycur,M2)+1) + wpdata(mod(xcur,xn)+1,mod(ycur,yn)+1) .* pou';
      end
    end

    %put back into fd
    fd(mod(idx1,A1)+1,mod(idx2,A2)+1) = fd(mod(idx1,A1)+1,mod(idx2,A2)+1) + pass .* fh(mod(idx1,M1)+1,mod(idx2,M2)+1);
  end

  if(1)
    s = 1;
    R1 = 2^(s-ns)*G1;
    R2 = 2^(s-ns)*G2;
    idx1 = [ceil(-R1):floor(R1)];
    [wl,wr] = cvwindow((idx1+R1/1)/(R1/2));    tmpa = wl;
    [wl,wr] = cvwindow((idx1-R1/2)/(R1/2));    tmpb = wr;
    coef1 = tmpa.*tmpb;
    idx2 = [ceil(-R2):floor(R2)];
    [wl,wr] = cvwindow((idx2+R2/1)/(R2/2));    tmpa = wl;
    [wl,wr] = cvwindow((idx2-R2/2)/(R2/2));    tmpb = wr;
    coef2 = tmpa.*tmpb;
    pass = coef1'*coef2;
    [M1,M2] = size(pass);

    cs = c{s};
    tmp = fft2(cs{1}) / sqrt(numel(cs{1}));
    tmp = tmp/4;
    tmp = mescatter(tmp,0);
    tmp = mescatter(tmp',0)';
    [K1,K2] = size(tmp);
    fh = zeros(M1,M2);
    fh(mod(idx1,M1)+1,mod(idx2,M2)+1) = tmp(mod(idx1,K1)+1,mod(idx2,K2)+1);

    fd(mod(idx1,A1)+1,mod(idx2,A2)+1) = fd(mod(idx1,A1)+1,mod(idx2,A2)+1) + pass .* fh(mod(idx1,M1)+1,mod(idx2,M2)+1);
  end

  fd = mecombine(fd,E1);
  fd = mecombine(fd',E2)';

  x = idct2(fd);

author: thomaskuestner


text:
function [probas, Isubj,nf]=taldir(dirname, th_pval, DirTable)
%
% Computes the probability of the Talairach transform matrices
%    of all the subjects found in the directory "dirname".
%  Uses the mean vector and covariance matrix obtained with talairachin_table.m from
%     the data set (default data set: /space/neo/2/recon/buckner)
%  Uses th_pval as a threshold for the p-values to detect unlikely transform  matrices
%

%
% talairaching_dir_afd.m
%
% Original Author: Laurence Wastiaux
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.5 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

if (nargin<2 | nargin>3)
    msg=sprintf('USAGE: [probas, Isubj , nf]=taldir(SubjectDir, th_pval, <DirTable>)');
    disp(msg)
end

files=dir(dirname);
nf=files;
probas=[];
Isubj=[];
Pval=[];
count=0;
probas_flag=1;

if(nargin==3)
    fsafdDir=DirTable;
else
    %%% Get the default tables'directory %%%
    if(getenv('FREESURFER_HOME'))
        fsh=getenv('FREESURFER_HOME');
        fsafdDir=strcat(fsh, '/fsafd');
    else
        error(sprintf('Impossible to find FREESURFER_HOME\n'));
    end
end

%%% 1x9 mean vector obtained from the training set %%%
%load('/space/okapi/3/data/laurence/ADF/talairaching/transfo_param_mean.mat');
%file_mu='/space/okapi/3/data/laurence/ADF/talairaching/TalairachingMean.adf'; %Loads mu
file_mu=strcat(fsafdDir, '/TalairachingMean.adf');
fi=fopen(file_mu);
pos=0;
if(fi==-1)
    mess=sprintf('Could not find %s', file_mu);
    error(mess)
else
    while(strfind(fgetl(fi), '#'))  % skip the header
        pos=ftell(fi);
    end
    fseek(fi, pos, 'bof');
    mu=(fscanf(fi, '%g'))';
    fclose(fi);
end

%%% 9x9 covariance matrix obtained from the training set %%%
%load('/space/okapi/3/data/laurence/ADF/talairaching/transfo_param_regularizedCov2.mat'); %loads sigma
%sigma_file='/space/okapi/3/data/laurence/ADF/talairaching/TalairachingCovariance.adf';
sigma_file=strcat(fsafdDir, '/TalairachingCovariance.adf');
fis=fopen(sigma_file);
pos=0;
if(fis==-1)
    mess=sprintf('Could not find %s', sigma_file);
    error(mess)
else
    while(strfind(fgetl(fis), '#'))  % skip the header
        pos=ftell(fis);
    end
    fseek(fis, pos, 'bof');
    sig=fscanf(fis, '%g');
    sigma=reshape(sig, [9,9]);
    fclose(fis);
end

%%% Probabilities of the transform matrices  %%%
%stat_file='/space/okapi/3/data/laurence/ADF/talairaching/TalairachingProbas.adf';
stat_file=strcat(fsafdDir, '/TalairachingProbas.adf');
fid=fopen(stat_file);
pos=0;
if(fid==-1)
    mess=sprintf('Could not find %s... compute probas for the first time', stat_file);
    disp(mess)
    probas_flag=0;
else
    while(strfind(fgetl(fid), '#'))  % skip the header
        pos=ftell(fid);
    end
    fseek(fid, pos, 'bof');
    yy=fscanf(fid, '%g');
    fclose(fid);
end

for i=1:(length(files))
    s=strcat(dirname,'/',files(i).name);
    ttfile=strcat(s,'/mri/transforms/talairach.xfm');
    fid=fopen(ttfile, 'r');
    if ( (fid ~= -1) && ( length(strfind(files(i).name,'0'))>=1 | (length(strfind(files(i).name,'1'))>=1 )))
        while feof(fid) == 0
            linef=fgetl(fid);
            nb=findstr(linef, 'Linear_Transform');
            if nb == 1
                pos=ftell(fid);
                break
            end
        end
        A=(fscanf(fid, '%g',12))';
        A=[A(1) A(2) A(3) A(5) A(6) A(7) A(9) A(10) A(11)]; % if trans=1, the translation parameters are not taken into account
        status=fclose(fid);
        p=mvnpdf(A,mu,sigma);
        probas=[probas p];
        Isubj=[Isubj i];
        nf(length(probas)).name=files(i).name;
        if(probas_flag)
            [pinf]=compute_pval(p, yy);
            Pval=[Pval pinf];
            if (pinf < th_pval)
                %mess=sprintf('Talairach Transform: %s failed (%g)', files(i).name, p);
                mess=sprintf('Talairach Transform: %s failed, (p=%g pval=%g)', files(i).name, p, pinf);
                disp(mess)
                count=count+1;
            else
                mess2=sprintf('Talairach Transform: %s OK, (p=%g pval=%g)', files(i).name, p, pinf);
                disp(mess2)
            end
        end
    else
        i=i+1;
        if (fid ~= -1)
            status=fclose(fid);
        end
    end
end

if (length(probas) == 0)
    messdir=sprintf('No subject found in this directory');
    disp(messdir)
end
messcount=sprintf('Number of unlikely Talairach transforms : %d', count);
disp(messcount)
% if(~probas_flag)
%     outprobas=strcat(DirTable, '/TalairachingProbas.adf');
%     save(outprobas, 'probas', '-ASCII');
% end
%y=probas;
%save('/space/okapi/3/data/laurence/talairaching/transfo_param_probas.mat', 'y');

% subfunction compute_pval() %
function [p_inf]=compute_pval(val,y)
%load('/space/okapi/3/data/laurence/ADF/talairaching/transfo_param_probas.mat'); %loads y
pas=0.05;
x=0:pas:1;
[h] = hist(y,x);
p = h/sum(h);
dinf=find(x<=val);
xinf=x(dinf);
pinf=p(1:length(xinf));
if (val>=0 & length(xinf) >1 )
    p_inf=trapz(xinf,pinf)/pas;
elseif (val>=0 & (length(xinf)<2))
    pas2=pas/10;
    x2=0:pas2:1;
    [h2] = hist(y,x2);
    p2 = h2/sum(h2);
    dinf2=find(x2<=val);
    xinf2=x2(dinf2);
    pinf2=p2(1:length(xinf2));
    if(length(xinf2)>1)
        p_inf=trapz(xinf2,pinf2)/pas2;
    else
        p_inf=0;
    end
else
    p_inf=0;

end


text:
function [R,C] = syncRot(T)

[~,L,Q] = proj_deformable_approx(T');
s = sign(L(find(abs(L)==max(abs(L)),1)));
C = s*L';
R = s*Q';
R(3,:) = cross(R(1,:),R(2,:));

% [R,C] = projectNonrigid(T);
% R(3,:) = cross(R(1,:),R(2,:));

author: geopavlakos


text:
function [status]=seof(HDR)
% SEOF checks for end of signal-file
%    status = seof(HDR)
%
% returns 1 if End-of-EDF-File is reached
% returns 0 otherwise
%
% See also: SOPEN, SREAD, SWRITE, SCLOSE, SSEEK, SREWIND, STELL, SEOF

% 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.

%	$Id: seof.m 2205 2009-10-27 12:18:15Z schloegl $
%	(C) 1997-2005,2007,2008 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

%status=feof(HDR.FILE.FID);  % does not work properly
%if HDR.FILE.POS~=HDR.AS.startrec+HDR.AS.numrec;

if strmatch(HDR.TYPE,{'CTF','RDF','EEG','AVG','SIGIF'}),
	%status=feof(EDF.FILE.FID);  % does not work properly
	%if EDF.FILE.POS~=EDF.AS.startrec+EDF.AS.numrec;
        status = (HDR.FILE.POS >= HDR.NRec);

elseif strmatch(HDR.TYPE,{'RG64','LABVIEW','Nicolet'}),
	status = (HDR.FILE.POS >= (HDR.AS.endpos-HDR.HeadLen));

elseif strmatch(HDR.TYPE,{'ACQ','AINF','BDF','BKR','BrainVision','CNT','CTF','EDF','ET-MEG','GDF','MIT','SMA','CFWB','DEMG','EEProbe-CNT','EEProbe-AVR','MFER','alpha','native','SCP','BCI2000','TMS32','WG1','Sigma'}),
	status = (HDR.FILE.POS >= HDR.SPR*HDR.NRec);

elseif strmatch(HDR.TYPE,{'EGI'}),
        if HDR.FLAG.TRIGGERED,
	        status = (HDR.FILE.POS >= HDR.NRec);
        else
                status = (HDR.FILE.POS >= HDR.SPR);
        end;

elseif strmatch(HDR.TYPE,{'FIF'}),
        [buf, status] = rawdata('next');
        status = strcmp(status,'eof');

else
	status=feof(HDR.FILE.FID);
end;


text:
function msns_l_d
%
%  Deletes the global data that has been generated by 'msns_l_i'.
%
%  Calling sequence:
%
%    msns_l_d
%
%
%  LYAPACK 1.0 (Thilo Penzl, May 1999)

global LP_NU

clear global LP_NU


text:
function Mo = spm_mesh_calc(Mi,Mo,f,varargin)
% Evaluate a function on a mesh's data
% FORMAT Mo = spm_mesh_calc(Mi,Mo,f,opts)
% Mi   - input filenames (char array or cellstr)
%        or cell array of gifti objects or patch structures
% Mo   - output filename
%        if empty, a gifti object is returned and not saved on disk
% f    - MATLAB expression to be evaluated (string or function handle)
%        (e.g., f = '(s1.*s2).^2' or f = @(s1,s2) (s1.*s2).^2)
% opts - optional list of pairs of property names and values
%        dmtx - read images into data matrix X [default: false]
%__________________________________________________________________________
% Copyright (C) 2015-2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_calc.m 7577 2019-04-24 08:59:56Z guillaume $

%-Check input arguments
%--------------------------------------------------------------------------
if ~isempty(Mo)
    Mo = spm_file(Mo,'ext','.gii');
    Mo = spm_file(Mo,'cpath');
end
if numel(varargin) == 1 && isstruct(varargin{1})
    varargin = reshape([fieldnames(varargin{1})';struct2cell(varargin{1})'],1,[]);
end
if mod(numel(varargin),2)
    error('Incorrect number of input arguments.');
end
opts.dmtx = false;
for i=1:2:numel(varargin)
    switch lower(varargin{i})
        case 'dmtx'
            opts.dmtx = logical(varargin{i+1});
        otherwise
            error('Unknown option %s.',varargin{i});
    end
end

%-Load data
%--------------------------------------------------------------------------
if ischar(Mi), Mi = cellstr(Mi); end
for i=1:numel(Mi)
    g = gifti(Mi{i});
    if ~isfield(g,'cdata')
        error('File %s does not contain data.',Mi{i});
    end
    D = full(g.cdata);
    if i==1
        nv = size(D);
    else
        if ~isequal(size(D),nv)
            error('Data dimension mismatch.');
        end
    end
    D = reshape(D,1,[]);
    if opts.dmtx
        X(i,:) = D;
    else
        eval(['s',num2str(i),'=D;']);
    end
end

%-Evaluate function
%--------------------------------------------------------------------------
if ischar(f)
    try
        eval(['S = ' f ';']);
    catch
        l = lasterror;
        error('%s\nCan''t evaluate "%s".',l.message,f);
    end
elseif isa(f,'function_handle')
    try
        if opts.dmtx
            S = feval(f,X);
        else
            list = sprintf('s%d,',1:numel(Mi)); list = list(1:end-1);
            eval(['S = feval(f,' list ');']);
        end
    catch
        l = lasterror;
        error('%s\nCan''t evaluate "%s".',l.message,func2str(f));
    end
else
    error('Unknown function input.');
end

%-Return or save output
%--------------------------------------------------------------------------
gs = gifti(reshape(S,nv));
g  = gifti(Mi{1});
if isfield(g,'vertices') && isfield(g,'faces')
    gs.vertices  = g.vertices;
    gs.faces     = g.faces;
elseif ~isempty(g.private.metadata)
    metadata     = g.private.metadata;
    name         = {metadata.name};
    if any(ismember(name,'SurfaceID'))
        metadata = metadata(ismember(name,'SurfaceID'));
        gs.private.metadata(1) = metadata;
    end
end
if isempty(Mo), Mo = gs; else save(gs,Mo,'ExternalFileBinary'); end


text:
function calc_bootfitrms_a2(a,time,timef,bootloops,maepi)
    % function calc_bootfitrms_a2(a,time,timef,bootloops,maepi);
    % --------------------------------------------------
    % Plots Ncum observed vs. Ncum modeled for specified time windows
    %
    % Input variables:
    % a         : earthquake catalog
    % time      : learning period fo fit Omori parameters
    % timef     : forecast period
    % bootloops : Number of bootstraps
    % maepi     : mainshock

    % S.Neukomm/ S.Wiemer / J.Woessner
    % last update: 04.08.03

report_this_filefun(mfilename('fullpath'));
    % Surpress warnings from fmincon
    warning off;

    %[m_main, main] = max(a.Magnitude);
    date_matlab = datenum(a.Date.Year,a.Date.Month,a.Date.Day,a.Date.Hour,a.Date.Minute,zeros(size(a,1),1));
    date_main = datenum(floor(maepi(3)),maepi(4),maepi(5),maepi(8),maepi(9),0);
    time_aftershock = date_matlab-date_main;

    % Select biggest aftershock earliest in time, but more than 1 day after mainshock
    fDay = 1;
    ft_c=fDay/365; % Time not considered to find biggest aftershock
    vSel = (a.Date > maepi(:,3)+ft_c & a.Date<= maepi(:,3)+time/365);
    mCat = a.subset(vSel);
    vSel = mCat(:,6) == max(mCat(:,6));
    vBigAf = mCat(vSel,:)
    if length(mCat(:,1)) > 1
        vSel = vBigAf(:,3) == min(vBigAf(:,3));
        vBigAf = vBigAf(vSel,:);
    end

    date_biga = datenum(floor(vBigAf(3)),vBigAf(4),vBigAf(5),vBigAf(8),vBigAf(9),0);
    fT1 = date_biga - date_main; % Time of big aftershock

    % Aftershock times
    l = time_aftershock(:) > 0;
    tas = time_aftershock(l);
    eqcatalogue = a.subset(l);

    l = tas <= time;
    time_as=tas(l);

    % Calculate uncertainty and mean values of p,c,and k
    %[pval, pstd, cval, cstd, kval, kstd, loopout] = brutebootloglike_a2(tas(l),bootloops,fT1)
    % Calculate p,c,k for dataset
    nMod = 1;
    %nMod = 2
    if nMod == 1

        [pval1, pval2, cval1, cval2, kval1, kval2] = bruteforcerms_a2(time_as,fT1,nMod)

        pval1 = round(100*pval1)/100;
        cval1 = round(100*cval1)/100;
        kval1 = round(10*kval1)/10;
        kval2 = round(10*kval2)/10;

    else
        [pval1, pval2, cval1, cval2, kval1, kval2] = bruteforcerms_a2(time_as,fT1,nMod)
        pval1 = round(100*pval1)/100;
        pval2 = round(100*pval2)/100;
        cval1 = round(100*cval1)/100;
        cval2 = round(100*cval2)/100;
        kval1 = round(10*kval1)/10;
        kval2 = round(10*kval2)/10;
    end

    if (isnan(pval1) == 0 & isnan(pval2) == 0)

        figure_w_normalized_uicontrolunits('Numbertitle','off','Name','Forecast aftershock occurence')
        %     pval1 = 0.9;
        %     pval2 = 0.9;
        %     cval1 = 0.2;
        %     cval2 = 0.2;
        %     kval1 = 1000;
        %     kval2 = 2000;
        %loopout = [loopout , loopout(:,1)*0];

        %     % Times up to the forecast time
        %     lf = tas <= time+timef ;
        %     time_asf= [tas(lf) ];
        %     time_asf=sort(time_asf);
        %

        %
        %     % Compute the confidence limits
        %     for j = 1:length(loopout(:,1));
        %
        %         cumnr = (1:length(time_asf))'; cumnr_model = [];
        %         pvalb = loopout(j,1);
        %         cvalb = loopout(j,2);
        %         kvalb = loopout(j,3);
        %         for i=1:length(time_asf)
        %             if pval ~= 1
        %                 cm = kvalb/(pvalb-1)*(cvalb^(1-pvalb)-(time_asf(i)+cvalb)^(1-pvalb));
        %             else
        %                 cm = kvalb*log(time_asf(i)/cvalb+1);
        %             end
        %             cumnr_model = [cumnr_model; cm];
        %         end
        %         plot(time_asf,cumnr_model,'color',[0.8 0.8 0.8]);
        %         loopout(j,4) = max(cumnr_model);
        %         hold on
        %         %drawnow
        %     end
        %     % 2nd moment of bootstrap number of forecasted number of events
        %     fStdBst = calc_StdDev(loopout(:,4));
        %
        %     % now plot the forecast ...
        %     cumnrf = (1:length(time_asf))'; cumnr_modelf = [];
        %     for i=1:length(time_asf)
        %         if pval ~= 1
        %             cm = kval/(pval-1)*(cval^(1-pval)-(time_asf(i)+cval)^(1-pval));
        %         else
        %             cm = kval*log(time_asf(i)/cval+1);
        %         end
        %         cumnr_modelf = [cumnr_modelf; cm];
        %     end
        %     pf1 =  plot(time_asf,cumnr_modelf,'g-.','Linewidth',2);
        %     hold on
        %     pf2 =  plot(time_asf,cumnrf, 'b-','Linewidth',2);
        %
        % plot the best fit
        % Cumulative number of pbserved events
        cumnr = (1:length(tas))';
        cumnr_model = [];

        for i=1:length(tas)
            if tas(i) <= fT1
                if pval1 ~= 1
                    cm = kval1/(pval1-1)*(cval1^(1-pval1)-(tas(i)+cval1)^(1-pval1));
                else
                    cm = kval1*log(tas(i)/cval1+1);
                end
                cumnr_model = [cumnr_model; cm];
            else
                if (pval1 ~= 1 & pval2 ~= 1)
                    cm = kval1/(pval1-1)*(cval1^(1-pval1)-(tas(i)+cval1)^(1-pval1))+ kval2/(pval2-1)*(cval2^(1-pval2)-(tas(i)-fT1+cval2)^(1-pval2));
                else
                    cm = kval1*log(tas(i)/cval1+1) + kval2*log((tas(i)-fT1)/cval2+1);
                end
                cumnr_model = [cumnr_model; cm];
            end
        end
        tas=sort(tas);
        cumnr_model=sort(cumnr_model);
        p1 = plot(tas,cumnr_model,'r','Linewidth',2);
        hold on

        p2 = plot(tas,cumnr,'b','Linewidth',2,'Linestyle','--');
        %     string=['p = ' num2str(pval) '+-' num2str(pstd) '; c = ' num2str(cval) '+-' num2str(cstd) '; k = ' num2str(kval) '+-' num2str(kstd)];
        %     title(string)
        %
        %     % Plot observed events in forecast period from endpoint of modeled events in learning period
        %     vSel = time_asf >= max(time_as);
        %     vCumnr_forecast = cumnrf(vSel,:);
        %     vTime_forecast = time_asf(vSel,:);
        %     % Difference of modelled and observed number of events at time_as
        %     fDiff_timeas = cumnr_modelf(length(time_as))-cumnrf(length(time_as));
        %     vCumnr_forecast = vCumnr_forecast+fDiff_timeas;
        %     pf3 = plot(vTime_forecast, vCumnr_forecast,'m-.','Linewidth',2);
        %
        %
        % %     [Y, in] = sort(loopout(:,4));
        % %     %n5Conf = round(0.05*length(loopout))
        % %     %Y = find(in == n5Conf);
        % %     %Y = find(in == 5);
        % %     loops = loopout(in,:);
        % %     n5Conf = floor(min((find(round(loops(:,4))==round(prctile(loops(:,4),5))))));
        % %     %n5Conf = round(0.05*length(loops));
        % %     pvalb = loops(n5Conf,1);
        % %     cvalb = loops(n5Conf,2);
        % %     kvalb = loops(n5Conf,3);
        % %
        % %     cumnr = (1:length(time_as))'; cumnr_model = [];
        % %     for i=1:length(time_as)
        % %         if pval ~= 1
        % %             cm = kvalb/(pvalb-1)*(cvalb^(1-pvalb)-(time_as(i)+cvalb)^(1-pvalb));
        % %         else
        % %             cm = kvalb*log(time_as(i)/cvalb+1);
        % %         end
        % %         cumnr_model = [cumnr_model; cm];
        % %     end
        % %
        % %     %plot(time_as,cumnr_model,'k--','Linewidth',1);
        % %
        % %     %n95Conf = round(0.95*length(loopout));
        % %     %Y = find(in == n95Conf);
        % %     %loops = loopout(in,:);
        % %     n95Conf = ceil(max(find(round(loops(:,4))==round(prctile(loops(:,4),95)))));
        % %     %n95Conf = round(0.95*length(loops));
        % %     pvalb = loops(n95Conf,1);
        % %     cvalb = loops(n95Conf,2);
        % %     kvalb = loops(n95Conf,3);
        % %     cumnr = (1:length(time_as))'; cumnr_model = [];
        % %     for i=1:length(time_as)
        % %         if pval ~= 1
        % %             cm = kvalb/(pvalb-1)*(cvalb^(1-pvalb)-(time_as(i)+cvalb)^(1-pvalb));
        % %         else
        % %             cm = kvalb*log(time_as(i)/cvalb+1);
        % %         end
        % %         cumnr_model = [cumnr_model; cm];
        % %     end
        % %     %pc = plot(time_as,cumnr_model,'k--','Linewidth',1);
        %
        %     xlabel('Time [days]')
        %     ylabel('Cumulative number of aftershocks')
        %     xlim([0 max(time_asf)]);
        %
        %
        %     % calculate uncertainty sigma in forecasted number of aftershocks by
        %     % error propagation law
        %     time1 = time+timef;
        %     if pval == 1
        %         pv = 1-10^(-6);
        %     else
        %         pv = pval;
        %     end
        %     mpm1 = 1-pv;
        %     t1c = time1+cval;
        %     t0c = time+cval;
        %     sigma = (((-t1c^mpm1+t0c^mpm1)/(pv-1)*kstd)^2+...
        %         (kval/(pv-1)*(-t1c^mpm1*mpm1/t1c+t0c^mpm1*mpm1/t0c)*cstd)^2+...
        %         (kval/(pv-1)*(t1c^mpm1*log(t1c)+t1c^mpm1/(pv-1)-t0c^mpm1*log(t0c)-t0c^mpm1/(pv-1))*pstd)^2)^0.5;
        %     % Plot standard deviation error propagation law
        %     ps1=errorbar(max(time_asf),max(cumnr_modelf),sigma,sigma);
        %     set(ps1,'Linewidth',4,'Color',[0 1 0])
        %     % Plot standard deviation from bootstrap
        %     ps2=errorbar(max(time_asf),max(cumnr_modelf),fStdBst,fStdBst);
        %     set(ps2,'Linewidth',2,'Color',[1 0 0])
        %
        %     legend([p2 p1 pf1 pf3 min(ps1) min(ps2)],'data','model','forecast','observed','\sigma (Epl)','\sigma (Bst)',0)
        %
        %     % Title
        %     % Rate change from error propagation law
        %     % Find amount of events in forecast period for modeled data
        %     nummod = max(cumnr_modelf)-cumnr_modelf(length(time_as));
        %     % Find amount of  events in forecast period for observed data
        %     l = time_asf <=time+timef & time_asf > time;
        %     numreal = sum(l); % observed number of aftershocks
        %     fRc_Flaw = (numreal-nummod)/sigma;
        %     fRc_Bst = (numreal-nummod)/fStdBst;
        %     string=['p = ' num2str(pval) '+-' num2str(pstd) '; c = ' num2str(cval) '+-' num2str(cstd) '; k = ' num2str(kval) '+-' num2str(kstd)];
        %     title(string)
        %     string=['\sigma(Epl) = ' num2str(sigma) '; \sigma(Bst) = ' num2str(fStdBst)];
        %     text(0.25,100,string)
        %     % Set line for learning period
        %     yy = get(gca,'ylim');
        %     plot([max(time_as) max(time_as)],[0 yy(2)],'k-.')
    else
        disp('no result')
    end


text:
%%
%   Author: Yash Bansod
%
% GitHub: <https://github.com/YashBansod>

%% Clear the environment and the command line
clear;
close all;
clc;

%% Initialize the maps, color maps, start points and destination points
map_size = [20, 20];
input_map = false(map_size);              % Create an Input Map

input_map (6:7, 8:9) = 1;           % Add an obstacle
input_map (6:9, 12) = 1;            % Add another obstacle
input_map (12, 5:10) = 1;           % Add another obstacle
input_map (13:14, 13:16) = 1;       % Add another obstacle
input_map (16:18, 6:8) = 1;         % Add another obstacle
input_map (4:6, 17:18) = 1;         % Add another obstacle

start_coords = [3, 6];              % Save the location of start coordinate
dest_coords = [18, 15];             % Save location of destination coordinate

drawMapEveryTime = true;            % To see how nodes expand on the grid

cmap = [1   1   1;                  % Create a color map
    0   0   0;
    1   0   0;
    0   0   1;
    0   1   0;
    1   1   0;
    0.5 0.5 0.5];

colormap(cmap);                     % Sets the colormap for the current figure
[nrows, ncols] = size(input_map);   % Save the size of the input_map

map = zeros(nrows,ncols);           % Create map to save the states of each grid cell
map(~input_map) = 1;                % Mark free cells on map
map(input_map) = 2;                 % Mark obstacle cells on map

start_node = sub2ind(size(map), start_coords(1), start_coords(2));      % Generate linear indices of start node
dest_node = sub2ind(size(map), dest_coords(1), dest_coords(2));         % Generate linear indices of dest node

map(start_node) = 5;                % Mark start node on map
map(dest_node) = 6;                 % Mark destination node on map

distanceFromStart = Inf(nrows,ncols);   % Initialize distance from start array to inifinity
distanceFromEnd = Inf(nrows,ncols);     % Initialize distance from end array to inifinity

parent = zeros(nrows, ncols);           % Create a map for holding parent's index for each grid cell

distanceFromStart(start_node) = 0;  % distance of start node from start is zero
distanceFromEnd(dest_node) = 0;     % distance of end node from end is zero

% Update the values of all grid pixels for distance from end
[X, Y] = meshgrid(1:ncols, 1:nrows);
xd = dest_coords(1);
yd = dest_coords(2);
distanceFromEnd = abs(X - yd) + abs(Y - xd);    % Manhattan Distance

image([0.5, map_size(1)-0.5], [0.5, map_size(2)-0.5], map);
ax = gca;
ax.YTick = 0:1:map_size(1);
ax.XTick = 0:1:map_size(2);
grid on;                        % Display grid lines
% drawnow limitrate nocallbacks;  % Update figure
drawnow;

%% Process the map to update the parent information and distance from start
while true                              % Create an infinite loop
    map(start_node) = 5;                % Mark start node on map
    map(dest_node) = 6;                 % Mark destination node on map

    if (drawMapEveryTime)
        image([0.5, map_size(1)-0.5], [0.5, map_size(2)-0.5], map);
        ax = gca;
        ax.YTick = 0:1:map_size(1);
        ax.XTick = 0:1:map_size(2);
        grid on;                        % Display grid lines
        % drawnow limitrate nocallbacks;  % Update figure
        drawnow;
    end

    % Find the node with the minimum heuristic distance
    heuristicDist = distanceFromStart + distanceFromEnd;
    [min_dist, current] = min(heuristicDist(:));

    % Compute row, column coordinates of current node from linear index
    [i, j] = ind2sub(size(heuristicDist), current);

    % Create an exit condition for the infinite loop to end
    if ((current == dest_node) || isinf(min_dist))
        break
    end

    % Update distance value of element right of current element
    if (i+1 <= nrows && distanceFromStart(i+1, j) > distanceFromStart(i,j) + 1)
        if (parent(i+1, j) == 0 && input_map(i+1,j)~=1 && parent(current)~= sub2ind(size(map), i+1, j))
            distanceFromStart(i+1, j) = distanceFromStart(i,j) + 1;
            map(sub2ind(size(map), i+1, j)) = 4;    % Mark the neighbour of current as processing
            parent(i+1, j)= current;
        end
    end

    % Update distance value of element left of current element
    if (i-1 >= 1 && distanceFromStart(i-1, j) > distanceFromStart(i,j) + 1)
        if (parent(i-1, j) == 0 && input_map(i-1,j)~=1 && parent(current)~= sub2ind(size(map), i-1, j))
            distanceFromStart(i-1, j) = distanceFromStart(i,j) + 1;
            map(sub2ind(size(map), i-1, j)) = 4;    % Mark the neighbour of current as processing
            parent(i-1, j)= current;
        end
    end

    % Update distance value of element top of current element
    if (j-1 >= 1 && distanceFromStart(i, j-1) > distanceFromStart(i,j) + 1)
        if (parent(i, j-1) == 0 && input_map(i,j-1)~=1 && parent(current)~= sub2ind(size(map), i, j-1))
            distanceFromStart(i, j-1) = distanceFromStart(i,j) + 1;
            map(sub2ind(size(map), i, j-1)) = 4;    % Mark the neighbour of current as processing
            parent(i, j-1)= current;
        end
    end

    % Update distance value of element bottom of current element
    if (j+1 <= ncols && distanceFromStart(i, j+1) > distanceFromStart(i,j) + 1)
        if (parent(i, j+1) == 0 && input_map(i,j+1)~=1 && parent(current)~= sub2ind(size(map), i, j+1))
            distanceFromStart(i, j+1) = distanceFromStart(i,j) + 1;
            map(sub2ind(size(map), i, j+1)) = 4;    % Mark the neighbour of current as processing
            parent(i, j+1)= current;
        end
    end

    distanceFromStart(current) = -log(0);   % change the distance of current from start as infinity
    map(current) = 3;                       % mark the current point as processed
end

%% Construct route from start to dest by following the parent links
if (isinf(distanceFromStart(dest_node)))
    route = [];    % if distance to destination node is infinity
else
    route = [dest_node];                               % else backtrace the route from destination node
    while (parent(route(1)) ~= 0)                       % check front of route for start node
        route = [parent(route(1)), route];              % add parent of current node to front of route
    end

    for k = 2:length(route) - 1         % To visualize the map and the path
        map(route(k)) = 7;
        pause(0.001);                   % Pause the code for a while
        image([0.5, map_size(1)-0.5], [0.5, map_size(2)-0.5], map);
        ax = gca;
        ax.YTick = 0:1:map_size(1);
        ax.XTick = 0:1:map_size(2);
        grid on;                        % Display grid lines
        % drawnow limitrate nocallbacks;  % Update figure
        drawnow;
    end
end

% Add legends to the colormapped image
hold on;
for K = 1:7
    hidden_h(K) = surf(zeros(2, 2), 'edgecolor', 'none', ...
        'facecolor', cmap(K, :));
end
hold off

uistack(hidden_h, 'bottom');
legend(hidden_h, {'free space', 'obstacles', 'closed nodes', ...
    'open nodes', 'start node', 'goal node', 'shortest path'} )

author: YashBansod

repo: Robotics-Planning-Dynamics-and-Control


text:
function value = pyramid_unit_volume_3d ( )

%*****************************************************************************80
%
%% PYRAMID_UNIT_VOLUME_3D returns the volume of a unit pyramid in 3D.
%
%  Integration region:
%
%    - ( 1 - Z ) <= X <= 1 - Z
%    - ( 1 - Z ) <= Y <= 1 - Z
%              0 <= Z <= 1.
%
%  Discussion:
%
%    A pyramid with square base can be regarded as the upper half of a
%    3D octahedron.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, the volume of the pyramid.
%
  value = 4.0 / 3.0;

  return
end


text:
function [res]=convert_mc(ys,k)

% function [res]=convert_mc(ys,[k])
%
%  convert_mc  - converts multi-class labels from examples x classes matrix
%                of {+1,-1}s to examples x 1 matrix of integers, and back
%                can provide optional k to say how many classes
%                there should be (incase some are missing)
%
%  e.g:
%  convert_mc([1 2 2]') returns [1 -1; -1 1 ; -1 1] and vice-versa
%  if you pass a data object instead, it converts the Y component
%  in the same way and returns the data object with converted Ys.

  if nargin==1 k=0; end;

  d=[];
  if isa(ys,'data')
    d=ys;
    ys=d.Y;
  end

  expand=1;
  if size(ys,2)==1 & min(min(ys))>0
    res=-ones(size(ys,1),max(max(ys),k));
    for i=1:size(ys,1) res(i,ys(i))=1; end;
  else
    expand=0;
    if size(ys,2)==1
      res=ys; res(res==-1)=2;
    else
      for i=1:size(ys,1) res(i)=find(ys(i,:)==1); end;
      res=res';
    end
  end

  if ~isempty(d)
    d.Y=res;
    res=d;
  end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore