text:
function [y, Sigma_y, beta] = ml_gmr(Priors, Mu, Sigma, x, in, out)
%
% This function performs Gaussian Mixture Regression (GMR), using the
% parameters of a Gaussian Mixture Model (GMM). Given partial input data,
% the algorithm computes the expected distribution for the resulting
% dimensions. By providing temporal values as inputs, it thus outputs a
% smooth generalized version of the data encoded in GMM, and associated
% constraints expressed by covariance matrices.
%
% Inputs -----------------------------------------------------------------
%   o Priors:  1 x K array representing the prior probabilities of the K GMM
%              components.
%   o Mu:      D x K array representing the centers of the K GMM components.
%   o Sigma:   D x D x K array representing the covariance matrices of the
%              K GMM components.
%   o x:       P x N array representing N datapoints of P dimensions.
%   o in:      1 x P array representing the dimensions to consider as
%              inputs.
%   o out:     1 x Q array representing the dimensions to consider as
%              outputs (D=P+Q).
% Outputs ----------------------------------------------------------------
%   o y:       Q x N array representing the retrieved N datapoints of
%              Q dimensions, i.e. expected means.
%   o Sigma_y: Q x Q x N array representing the N expected covariance
%              matrices retrieved.
%
% Copyright (c) 2006 Sylvain Calinon, LASA Lab, EPFL, CH-1015 Lausanne,
%               Switzerland, http://lasa.epfl.ch
%
% The program is free for non-commercial academic use.
% Please contact the authors if you are interested in using the
% software for commercial purposes. The software must not be modified or
% distributed without prior permission of the authors.
% Please acknowledge the authors in any academic publications that have
% made use of this code or part of it. Please use this BibTex reference:
%
% @article{Calinon06SMC,
%   title="On Learning, Representing and Generalizing a Task in a Humanoid
%     Robot",
%   author="S. Calinon and F. Guenter and A. Billard",
%   journal="IEEE Transactions on Systems, Man and Cybernetics, Part B.
%     Special issue on robot learning by observation, demonstration and
%     imitation",
%   year="2006",
%   volume="36",
%   number="5"
% }

nbData = size(x,2);
nbVar = size(Mu,1);
nbStates = size(Sigma,3);

%% Fast matrix computation (see the commented code for a version involving
%% one-by-one computation, which is easier to understand).
%%
%% Compute the influence of each GMM component, given input x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:nbStates
  Pxi(:,i) = Priors(i).*ml_gaussPDF(x, Mu(in,i), Sigma(in,in,i));
end
beta = Pxi./repmat(sum(Pxi,2)+realmin,1,nbStates);

% ind = find(sum(beta,2) == 0);
% if ~isempty(ind)
%     for i=1:nbStates
%         tmp = x(:,ind)' - repmat(Mu(in,i)',length(ind),1);
%         score(i,:) = Priors(i)*(sum((tmp/Sigma(in,in,i)).*tmp, 2));
%     end
%     [i i]=min(score);
%     beta=beta';
%     beta(nbStates.*(ind-1)+i')=1;
%     beta=beta';
% end
%% Compute expected means y, given input x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:nbStates
  y_tmp(:,:,j) = repmat(Mu(out,j),1,nbData) + Sigma(out,in,j)/(Sigma(in,in,j)) * (x-repmat(Mu(in,j),1,nbData));
end
beta_tmp = reshape(beta,[1 size(beta)]);
y_tmp2 = repmat(beta_tmp,[length(out) 1 1]) .* y_tmp;
y = sum(y_tmp2,3);
%% Compute expected covariance matrices Sigma_y, given input x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This seems to be the computation of the marginal covariance matrix SigmaY
% if nargout > 1
%     for j=1:nbStates
%         Sigma_y_tmp(:,:,1,j) = Sigma(out,out,j) - (Sigma(out,in,j)/(Sigma(in,in,j))*Sigma(in,out,j));
%     end
%     beta_tmp = reshape(beta,[1 1 size(beta)]);
%     Sigma_y_tmp2 = repmat(beta_tmp.*beta_tmp, [length(out) length(out) 1 1]) .* repmat(Sigma_y_tmp,[1 1 nbData 1]);
%     Sigma_y = sum(Sigma_y_tmp2,4);
% end

% %%%%%% Compute expected output distribution, given input x^i %%%%%%
if nargout > 1
    M = nbData; K = nbStates;
    Sigma_y = zeros(length(out), length(out), M);
    for i=1:M
        %%%% Eq.11: Compute expected covariance matrices, given input x^i %%%%
        for k=1:K
            % Full conditional variance equations (Nadia's way)
            mu_k_tmp     = Mu(out,k) + Sigma(out,in,k)*inv(Sigma(in,in,k)) * (x(:,i)-Mu(in,k));
            Sigmak_y_tmp = Sigma(out,out,k) - (Sigma(out,in,k)*inv(Sigma(in,in,k))*Sigma(in,out,k));
            Sigma_y(:,:,i) = Sigma_y(:,:,i) + beta(i,k) .* (mu_k_tmp'*mu_k_tmp + Sigmak_y_tmp);
        end

        % Full conditional variance equations (Nadia's way)
        Sigma_y(:,:,i) = Sigma_y(:,:,i) - (y(:,i)'*y(:,i));
    end
end

% %% Slow one-by-one computation (better suited to understand the algorithm)
% %%
% %% Compute the influence of each GMM component, given input x
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for i=1:nbStates
%   Pxi(:,i) = gaussPDF(x, Mu(in,i), Sigma(in,in,i));
% end
% beta = (Pxi./repmat(sum(Pxi,2)+realmin,1,nbStates))';
% %% Compute expected output distribution, given input x
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% y = zeros(length(out), nbData);
% Sigma_y = zeros(length(out), length(out), nbData);
% for i=1:nbData
%   % Compute expected means y, given input x
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   for j=1:nbStates
%     yj_tmp = Mu(out,j) + Sigma(out,in,j)*inv(Sigma(in,in,j)) * (x(:,i)-Mu(in,j));
%     y(:,i) = y(:,i) + beta(j,i).*yj_tmp;
%   end
%   % Compute expected covariance matrices Sigma_y, given input x
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   for j=1:nbStates
%     Sigmaj_y_tmp = Sigma(out,out,j) - (Sigma(out,in,j)*inv(Sigma(in,in,j))*Sigma(in,out,j));
%     Sigma_y(:,:,i) = Sigma_y(:,:,i) + beta(j,i)^2.* Sigmaj_y_tmp;
%   end
% end

author: epfl-lasa


text:
%DISTM Compute square Euclidean distance matrix
%
%   D = DISTM(A,B)
%
% INPUT
%   A,B   Datasets or matrices; B is optional, default B = A
%
% OUTPUT
%   D     Square Euclidean distance dataset or matrix
%
% DESCRIPTION
% Computation of the square Euclidean distance matrix D between two
% sets A and B. If A has M objects and B has N objects, then D is
% [M x N]. If A and B are datasets, then D is a dataset as well with
% the labels defined by the labels of A and the feature labels defined
% by the labels of B.
%
% Unlabeled objects in B are neglected, unless B is entirely unlabeled.
%
% If A is not a dataset, but a matrix of doubles then D is also not a
% dataset, but a set of doubles.
%
% NOTE
% DISTM(A,B) is equivalent to A*PROXM(B,'d',2)).
%
% SEE ALSO
% DATASETS, PROXM

% Copyright: R.P.W. Duin, r.p.w.duin@prtools.org
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: distm.m,v 1.9 2006/03/07 16:14:41 duin Exp $

function D = distm_new(A,B,DirectionalFeatures)
	prtrace(mfilename);
	if nargin < 2
		B = A;
	end
	B = cdats(B,1);
	[ma,ka] = size(A);
	[mb,kb] = size(B);

	if (ka ~= kb)
		error('Feature sizes should be equal')
    end

    if( isempty(DirectionalFeatures) )
        % The order of operations below is good for the accuracy.
        D = ones(ma,1)*sum(B'.*B',1);
        D = D + sum(A'.*A',1)'*ones(1,mb);
        D = D - 2 .* (+A)*(+B)';
    else

        %simulo los dos wraps posibles
        B1 = B(:,DirectionalFeatures) + 2*pi;
        B2 = B(:,DirectionalFeatures) - 2*pi;

        %mido las distancias con y sin los 2 wraps
        D = ones(ma,1)*sum(B'.*B',1);
        D = D + sum(A'.*A',1)'*ones(1,mb);
        D = D - 2 .* (+A)*(+B)';

        D1 = ones(ma,1)*sum(B1'.*B1',1);
        D1 = D1 + sum(A'.*A',1)'*ones(1,mb);
        D1 = D1 - 2 .* (+A)*(+B1)';

        D2 = ones(ma,1)*sum(B2'.*B2',1);
        D2 = D2 + sum(A'.*A',1)'*ones(1,mb);
        D2 = D2 - 2 .* (+A)*(+B2)';

        %me quedo con la menor distancia.
        D = reshape(min([ D(:) D1(:) D2(:) ],[],2), ma, ma);

    end

	J = find(D<0);                  % Check for a numerical inaccuracy.
	D(J) = zeros(size(J));          % D should be nonnegative.

	if ((nargin < 2) & (ma == mb)) % take care of symmetric distance matrix
		D = (D + D')/2;
		D([1:ma+1:ma*ma]) = zeros(1,ma);
	end

	if isa(A,'prdataset')   % set object and feature labels
		if isa(B,'prdataset')
			D = setdata(A,D,getlab(B));
		else
			D = setdata(A,D);
		end
	end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% difftopo - compute and plot component decomposition for the difference ERP
%            between two EEG datasets. Plots into the current axis (gca);
%            plot into a new empty figure as below.
% Usage:
%          >> figure; difftopo(ALLEEG,eeg1,eeg2,interval);
% Inputs:
%         ALLEEG - array of leaded EEG datasets
%         eeg1   - index of subtrahend (+) dataset
%         eeg2   - index of minuend (-) dataset
%         interval - [minms maxms] latency interval in which to find
%                  and plot largest contributing components {default: whole epoch}
%         limits - [stms endms] latenc plotting limits (in ms) {default|0: whole epoch}
%         subcomps - array of component indices to remove from data
%                  (e.g. eye movement components) {default|0: none}
%
% Outputs: none
%
% Author: Scott Makeig, SCCN/INC/UCSD 3/28/05

% Copyright (C) Scott Makeig, SCCN/INC/UCSD 3/28/05
%
% 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
function difftopo(ALLEEG,eeg1,eeg2,interval,limits,subcomps);

if nargin < 3
   help difftopo
   return
end
if eeg1 < 1 | eeg1 > length(ALLEEG)
   help difftopo
   return
end
if eeg2 < 1 | eeg2 > length(ALLEEG)
   help difftopo
   return
end
if eeg1 == eeg2
   help difftopo
   return
end
if ndims(ALLEEG(eeg1).data) ~=3  | ndims(ALLEEG(eeg2).data) ~=3
  error('EEG datasets must be epoched data');
end

if nargin < 4
   interval = [ALLEEG(eeg1).xmin*999.9 ALLEEG(eeg1).xmax*999.9];
end
if nargin < 5
   limits = 0; % [ALLEEG(eeg1).xmin*1000 ALLEEG(eeg1).xmax*1000];
end
if nargin < 6
   subcomps = [];
end

if ~isfield(ALLEEG(eeg1),'icaweights')
  error('EEG datasets must have icaweights');
end
if length(interval) ~= 2
   help difftopo
   return
end

if ALLEEG(eeg1).pnts ~= ALLEEG(eeg1).pnts
  error('EEG datasets must have the same number of frames per epoch');
end

set(gcf,'Name','difftopo()');
diff = mean(ALLEEG(eeg1).data,3) - mean(ALLEEG(eeg2).data,3);

plottitle = [ ALLEEG(eeg1).setname ' - ' ALLEEG(eeg2).setname];
envtopo(diff,ALLEEG(eeg1).icaweights*ALLEEG(eeg1).icasphere,...
       'chanlocs',ALLEEG(eeg1).chanlocs, ...
       'timerange',[ALLEEG(eeg1).xmin*1000 ALLEEG(eeg1).xmax*1000],...
       'limits',limits,...
       'limcontrib',interval,...
       'title',plottitle, ...
       'subcomps',subcomps);

author: PatternRecognition

repo: OpenBMI


text:
function [grayscale_sequence] = sequence_grayscale(sequence)
% sequence_grayscale Returns grayscale sequence
%
% This sequence converter returns a grayscale version of an input sequence.
%
% Cache notice: The results of this function are cached in the workspace cache directory.
%
% Input:
% - sequence (structure): A valid sequence structure.
%
% Output:
% - grayscale_sequence (structure): A sequence descriptor of a converted sequence.

cache_directory = fullfile(get_global_variable('directory'), 'cache', 'grayscale', sequence.name);

mkpath(cache_directory);

cache_groundtruth = fullfile(cache_directory, 'groundtruth.txt');

sequence_groundtruth = fullfile(sequence.directory, 'groundtruth.txt');

if file_newer_than(cache_groundtruth, sequence_groundtruth)
    grayscale_sequence = sequence_create(cache_directory, 'name', sequence.name);
    grayscale_sequence.tags.names = sequence.tags.names;
    grayscale_sequence.tags.data = sequence.tags.data;
    grayscale_sequence.properties = sequence.properties;
    return;
end;

print_debug('Generating cached grayscale sequence ''%s''...', sequence.name);

for i = 1:sequence.length

    color_image = imread(sequence_get_image(sequence, i));

    gray_image = repmat(rgb2gray(color_image), [1 1 3]);

    imwrite(gray_image, fullfile(cache_directory, sprintf('%08d.jpg', i)));

end;

write_trajectory(cache_groundtruth, sequence.groundtruth);

grayscale_sequence = sequence_create(cache_directory, 'name', sequence.name);

grayscale_sequence.tags.names = sequence.tags.names;
grayscale_sequence.tags.data = sequence.tags.data;

author: votchallenge

repo: toolkit-legacy


text:
function ws = wavelet_wavescale(w)
	ws = w.wavScale;

author: thomaskuestner


text:
function linplus_test295 ( )

%*****************************************************************************80
%
%% TEST295 tests R8GE_DILU.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  ncol = 3;
  nrow = 3;
  n = nrow * ncol;
  m = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST295\n' );
  fprintf ( 1, '  For a matrix in general storage,\n' );
  fprintf ( 1, '  R8GE_DILU returns the DILU factors.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =    %d\n', m );
  fprintf ( 1, '  Matrix columns N = %d\n', n );

  for i = 1 : nrow * ncol
    for j = 1 : nrow * ncol

      if ( i == j )
        a(i,j) = 4.0E+00;
      elseif ( i == j + 1 | i == j - 1 | i == j + nrow | i == j - nrow )
        a(i,j) = -1.0E+00;
      else
        a(i,j) = 0.0E+00;
      end

    end
  end

  r8ge_print ( m, n, a, '  Matrix A:' );
%
%  Compute the incomplete LU factorization.
%
  d = r8ge_dilu ( m, n, a );

  r8vec_print ( m, d, '  DILU factor:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ y, m ] = month_carry_bahai ( y, m )

%*****************************************************************************80
%
%% MONTH_CARRY_BAHAI carries a year of months on the Bahai calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, the YM date.
%
  months = year_length_months_bahai ( y );

  while ( 1 )

    if ( m <= months )
      break;
    end if

    m = m - months;
    y = y + 1;

  end

  return
end


text:
function [im] = nonmax(im,theta)
% function [im] = nonmax(im,theta)
%
% Perform non-max suppression on im orthogonal to theta.  Theta can be
% a matrix providing a different theta for each pixel or a scalar
% proving the same theta for every pixel.
%
% David R. Martin <dmartin@eecs.berkeley.edu>
% March 2003

if numel(theta)==1,
  theta = theta .* ones(size(im));
end

% Do non-max suppression orthogonal to theta.
theta = mod(theta+pi/2,pi);

% The following diagram depicts the 8 cases for non-max suppression.
% Theta is valued in [0,pi), measured clockwise from the positive x
% axis.  The 'o' marks the pixel of interest, and the eight
% neighboring pixels are marked with '.'.  The orientation is divided
% into 8 45-degree blocks.  Within each block, we interpolate the
% image value between the two neighboring pixels.
%
%        .66.77.
%        5\ | /8
%        5 \|/ 8
%        .--o--.-----> x-axis
%        4 /|\ 1
%        4/ | \1
%        .33.22.
%           |
%           |
%           v
%         y-axis
%
% In the code below, d is always the distance from A, so the distance
% to B is (1-d).  A and B are the two neighboring pixels of interest
% in each of the 8 cases.  Note that the clockwise ordering of A and B
% changes from case to case in order to make it easier to compute d.

% Determine which pixels belong to which cases.
mask15 = ( theta>=0 & theta<pi/4 );
mask26 = ( theta>=pi/4 & theta<pi/2 );
mask37 = ( theta>=pi/2 & theta<pi*3/4 );
mask48 = ( theta>=pi*3/4 & theta<pi );

mask = ones(size(im));
[h,w] = size(im);
[ix,iy] = meshgrid(1:w,1:h);

% case 1
idx = find( mask15 & ix<w & iy<h);
idxA = idx + h;
idxB = idx + h + 1;
d = tan(theta(idx));
imI = im(idxA).*(1-d) + im(idxB).*d;
mask(idx(find(im(idx)<imI))) = 0;

% case 5
idx = find( mask15 & ix>1 & iy>1);
idxA = idx - h;
idxB = idx - h - 1;
d = tan(theta(idx));
imI = im(idxA).*(1-d) + im(idxB).*d;
mask(idx(find(im(idx)<imI))) = 0;

% case 2
idx = find( mask26 & ix<w & iy<h );
idxA = idx + 1;
idxB = idx + h + 1;
d = tan(pi/2-theta(idx));
imI = im(idxA).*(1-d) + im(idxB).*d;
mask(idx(find(im(idx)<imI))) = 0;

% case 6
idx = find( mask26 & ix>1 & iy>1 );
idxA = idx - 1;
idxB = idx - h - 1;
d = tan(pi/2-theta(idx));
imI = im(idxA).*(1-d) + im(idxB).*d;
mask(idx(find(im(idx)<imI))) = 0;

% case 3
idx = find( mask37 & ix>1 & iy<h );
idxA = idx + 1;
idxB = idx - h + 1;
d = tan(theta(idx)-pi/2);
imI = im(idxA).*(1-d) + im(idxB).*d;
mask(idx(find(im(idx)<imI))) = 0;

% case 7
idx = find( mask37 & ix<w & iy>1 );
idxA = idx - 1;
idxB = idx + h - 1;
d = tan(theta(idx)-pi/2);
imI = im(idxA).*(1-d) + im(idxB).*d;
mask(idx(find(im(idx)<imI))) = 0;

% case 4
idx = find( mask48 & ix>1 & iy<h );
idxA = idx - h;
idxB = idx - h + 1;
d = tan(pi-theta(idx));
imI = im(idxA).*(1-d) + im(idxB).*d;
mask(idx(find(im(idx)<imI))) = 0;

% case 8
idx = find( mask48 & ix<w & iy>1 );
idxA = idx + h;
idxB = idx + h - 1;
d = tan(pi-theta(idx));
imI = im(idxA).*(1-d) + im(idxB).*d;
mask(idx(find(im(idx)<imI))) = 0;

% apply mask
im = im .* mask;

author: Cloud-CV

repo: object-proposals


text:
function runCommissioned(energy, fieldSize, FFA, FFB, FFDistance, nhist, r, c, s, batch)
% 'runBinEnergy'
% JC 1/29/06
% How to use
%> runBinEnergy 18 20 4 100000000
%
% Written by JC. Dec 2005
% Run DPM for the different energy bins, for optimize the photon spectrum later.
%"input" is the .MAT file have DPMIN, CTscan, and indV
%DPMIN - DPM input structure
%CTscan - DPM scan structure
%indV = sub2ind(getUniformizedSize(planC), r, c, s); see DPMInfluence.m
%line 72
%DPMIN_CTscan_indV DMPIN, CTScan, indV
%"number" is how many bins we're going to divide the whole energy range
% Typical value for "number" is 20
%nhis - number of histories
%
% JC Oct 18, 2006  use target to flattening filter distance as 12.5cm instead of 9.5cm
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

setappdata(0, 'usenativesystemdialogs', false)
currentDir = cd;

[FileName,path] = uigetfile('*.mat','Select MAT file containing DPMIN, CTscan, & indV, p');

if path == 0
    errordlg('File Should exist');
    error('File Should exist');
end

cd(path);
load (FileName, 'DPMIN', 'CTscan', 'indV')
cd(currentDir);

%% Get the commissioned parameters from a .mAT file
[FileName,path] = uigetfile('*.mat','Select MAT file containing p ener a enerFF aFF');

if path == 0
    errordlg('File Should exist');
    error('File Should exist');
end

cd(path);
load (FileName, 'p', 'ener', 'a', 'enerFF', 'aFF')
cd(currentDir);

if (class(energy)== 'char')
    energy = str2num(energy);nhist = str2num(nhist); batch = str2num(batch);
    r = str2num(r); c = str2num(c); s = str2num(s); fieldSize = str2num(fieldSize);
end

[FileName,path] = uigetfile('*.*','Select dpm.dll to execute');

if path == 0
    errordlg('File Should exist');
    error('File Should exist');
end

cd(path)
% need ./local to run dpm. if no. crush.
% so creat it
mkdir local

rand('state',sum(100*clock));

% DPMIN.FlatFilterA = -0.21300000000000;
% DPMIN.FlatFilterB= -0.10000000000000;
% DPMIN.FlatFilterDist = 12.50000000000000;
DPMIN.FlatFilterA = str2num(FFA);
DPMIN.FlatFilterB= str2num(FFB);
DPMIN.FlatFilterDist = str2num(FFDistance);

% Add input ParticleType as 0 for photon or -1 for electron
DPMIN.ParticleType = 0;
DPMIN.SourceEnergy = energy;
% Jan 12, 2008
% Now "nhist" is defined as photons/cm^2 beamlet.
% Turn on Softening flag, for source model.
DPMIN.NumParticles = nhist*fieldSize*fieldSize;
DPMIN.Softening = 1;

DPMIN.IsoDistance = 100;
DPMIN.IncreaseFluence = 0;
DPMIN.OnlyHorn = 0;
DPMIN.UseFlatFilter = 0;
DPMIN.UsePhotSpectrum = 1;
DPMIN.OutputError = 0;
DPMIN.RandSeeds = [rand*1000 rand*10000];
DPMIN.PhotSpectrum =[ener', a'];

if (DPMIN.ParticleType == 0)
    DPMIN.Prefix = 'pre4phot'
elseif (DPMIN.ParticleType == -1)
    DPMIN.Prefix = 'pre4elec'
else
    error('DPMIN.ParticleType has to be 0 or -1')
end

% generratePB.
beamletSize = min(fieldSize,1);
generateOpenfieldPB(fieldSize, fieldSize, beamletSize);
filename = ['PB_', num2str(fieldSize), 'x', num2str(fieldSize), '_PBsize', num2str(beamletSize), 'cm.mat']
disp ('load PB info from above filename.')
load (filename, 'xPosV', 'yPosV', 'beamlet_delta_x', 'beamlet_delta_y', 'w_field');

% Load the Horn Coefficients.
Energy = '6MV10x10MDA.spectrum';
inds = max(strfind(Energy, '.'));
% filename:
fid = fopen([Energy(1:inds-1), '_CosHorn', Energy(inds:end)]);
% Here, the file extension ".spectrum" is only for the purpose of being
% consistent with other source model parameters files.
if (fid == -1),
    disp(['Can not Open CosHorn file',[Energy(1:inds-1), '_CosHorn', Energy(inds:end)]]);
end
tline = fgetl(fid);
tline = fgetl(fid);
CosHorn = fscanf(fid,'%g',[2 inf]); % It has two rows now.
CosHorn = CosHorn';
cosOffAxisAngle = DPMIN.IsoDistance ./ ...
    (sqrt((xPosV.^2+ yPosV.^2 + DPMIN.IsoDistance.^2)));
Horn = interp1(CosHorn(:,1), CosHorn(:,2), cosOffAxisAngle, 'linear');

DPMIN

doseV_PM_total = zeros(size(indV'));
doseV_EF_total = zeros(size(indV'));

for i = 1 : length(xPosV)

    source = DPMIN.SourceXYZ;

    % scale 40cm/100cm back from the iso center.
    % 0.5*0.4 = (half field size on the plus/minus sides; 0.5;
    %           (distance from source-to-jaw: 40cm / isodistance 100cm = 0.4)
    %  == 0.2
    % DPMIN.ClmtCorner1 = [source(1)-0.2*fieldSize    source(2)+40    source(3)+0.2*fieldSize];
    % DPMIN.ClmtCorner2 = [source(1)+0.2*fieldSize    source(2)+40    source(3)+0.2*fieldSize];
    % DPMIN.ClmtCorner3 = [source(1)+0.2*fieldSize    source(2)+40    source(3)-0.2*fieldSize];
    %load (filename, 'xPosV', 'yPosV', 'beamlet_delta_x', 'beamlet_delta_y', 'w_field');

    DPMIN.ClmtCorner1 = [source(1)-0.2*beamlet_delta_x(i)    source(2)+40    source(3)+0.2*beamlet_delta_y(i)];
    DPMIN.ClmtCorner2 = [source(1)+0.2*beamlet_delta_x(i)    source(2)+40    source(3)+0.2*beamlet_delta_y(i)];
    DPMIN.ClmtCorner3 = [source(1)+0.2*beamlet_delta_x(i)    source(2)+40    source(3)-0.2*beamlet_delta_y(i)];

    tic; doseV = dpm(DPMIN, CTscan); toc
    % Add the Horn coefficients
    doseV_PM_total = doseV_PM_total + (1+Horn(i))*doseV;

    %% Run only Flat Filter
    DPMIN.OnlyHorn = 0;
    DPMIN.UseFlatFilter = 1;
    %% Test, use 5 times particiles for EF.
    DPMIN.NumParticles = DPMIN.NumParticles*p(6)*5;
    DPMIN.PhotSpectrum = [enerFF', aFF'];
    DPMIN.RandSeeds = [rand*1000 rand*10000];
    tic; doseV = dpm(DPMIN, CTscan); toc
    doseV_EF_total = doseV_PM_total + doseV;

end

% At the end of the run, save the results.

cd(currentDir);
dose3D = zeros(r,c,s);
dose3D(indV) = doseV_PM_total;
filename = ['dose3D_PM_', num2str(batch),'.mat'];
save(filename, 'dose3D');

dose3Dsum = sum(a)*dose3D;

dose3D = zeros(r,c,s);
dose3D(indV) = doseV_EF_total;
filename =  ['dose3D_EF_',  num2str(batch), '.mat'];
save (filename, 'dose3D');
dose3Dsum = dose3Dsum + dose3D*sum(aFF);          % sum(aFF) == sum(a)*p(6);

filename =  ['dose3Dsum_',  num2str(batch), '.mat'];
save(filename,'dose3Dsum');

return;

% No Longer valide. Do not use.
% % % %% Run OnlyHorn
% % % DPMIN.NumParticles = DPMIN.NumParticles*p(7); % scale down the number of particles by the relative fluence weight of Horn.
% % % DPMIN.OnlyHorn = 1;
% % % DPMIN.RandSeeds = [rand*1000 rand*10000];
% % %
% % % cd(path);
% % % tic; doseV = dpm(DPMIN, CTscan); toc
% % % cd(currentDir);
% % %
% % % dose3D = zeros(r,c,s);
% % % dose3D(indV) = doseV;
% % % filename =  ['dose3D_Horn', num2str(batch), '.mat'];
% % % save (filename, 'dose3D');
% % % dose3Dsum = dose3Dsum + dose3D*sum(a)*p(7);

% Run electron
% Since for 6MV, the relative electron contribution is only about 0.1%, so
% it's negligible.
%  DPMIN.
%   DPMIN.Prefix = 'pre4elec';
%     DPMIN.UseFlatFilter = 0
%
%     cd(path);
%     tic; doseV = dpm(DPMIN, CTscan); toc
%     cd(currentDir);
%
%     dose3D = zeros(r,c,s);
%     dose3D(indV) = doseV;
%     filename = ['dose3D_elec.mat'];
%     save(filename,'dose3D');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function  bzip2Matlab

% This function prompts the user to select a file via gui method
%   If the file ends with *.bz2 suffix
 %      function assumes the file is compressed and will uncompress it
 %      using bzip2
%   Else if the file ends with any other suffix
 %      function attempts to compress it using bzip2
 %  The objective is easier handlng of compressed files in Matlab
%  LM AIB 11/24/2002

% copyright (c) 2001-2006, Washington University in St. Louis.
% Permission is granted to use or modify only for non-commercial,
% non-treatment-decision applications, and further only if this header is
% not removed from any file. No warranty is expressed or implied for any
% use whatever: use at your own risk.  Users can request use of CERR for
% institutional review board-approved protocols.  Commercial users can
% request a license.  Contact Joe Deasy for more information
% (radonc.wustl.edu@jdeasy, reversed).

[fname, pathname] = uigetfile('*.*', ...
    'Select a file you wish to compress or decompress');
tic
oldDir = pwd;
if isdeployed
    pathStr = fullfile(getCERRPath,'bin','Compression');
else
    pathStr = fullfile(getCERRPath,'Compression');
end
cd(pathname);
l = length(fname);
fmat = ''; outstr = '';
fmat = (strcat('"',oldDir, '\', fname,'"'));
if (~strcmpi(midstring(fname,l-2,l),'bz2')) % compress file
        cd(pathStr);
        outstr = ['bzip2-102-x86-win32.exe -vv9 ', fmat];
        system(outstr);
        cd(oldDir);
elseif (strcmpi(midstring(fname,l-2,l),'bz2')) % uncompress file
        cd(pathStr);
        outstr = ['bzip2-102-x86-win32.exe -dvv ', fmat];
        system(outstr);
        cd(oldDir);
else
    error('Incorrect filename chosen. Exiting bzip2Matlab.')
end

toc

return


text:
function plotprmm(catalogObject)
    %CATALOG.PLOTPRMM Plot the peak rate and maximum magnitude of a
    %set of events
    figure
    symsize = get_symsize(catalogObject);
    t=catalogObject.gettimerange();
    days = t(2) - t(1);
    if all(isnan(catalogObject.mag))
        warning('No magnitude data to plot');
    else

        % plot magnitudes
        subplot(2,1,1), scatter(catalogObject.otime, catalogObject.mag, symsize);
        %stem(catalogObject.otime, catalogObject.mag);
        set(gca, 'XLim', [floor(t(1)) ceil(t(2))]);
        datetick('x');
        xlabel('Date');
        ylabel('Magnitude');
        grid on;

        % put 'MM' label by max mag event
        [mm, mmi] = max(catalogObject.mag);
        text(catalogObject.otime(mmi), catalogObject.mag(mmi), 'MM','color','r');
        disp(sprintf('MM=%.1f occurs at %.1f%% of time series',mm,100*(catalogObject.otime(mmi) - t(1))/days));

        % plot event rate in 100 equal bins

        binsize = days/100;
        erobj = catalogObject.eventrate('binsize',binsize);
        %erobj = catalogObject.eventrate();
        subplot(2,1,2),plot(erobj.time, erobj.counts);
        set(gca, 'XLim', [floor(t(1)) ceil(t(2))]);
        datetick('x');
        xlabel('Date');
        ylabel('Event Rate');
        grid on;

        % put 'PR' label by peak-rate
        [pr, pri] = max(erobj.counts);
        text(erobj.time(pri), erobj.counts(pri), 'PR','color','r');
        disp(sprintf('PR=%d occurs at %.1f%% of time series',pr,100*(erobj.time(pri) - erobj.snum)/(erobj.enum-erobj.snum)));
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function disp(obj)
% display MoxUnitTestSuite object
%
% disp(obj)
%
% Inputs:
%   obj             MoxUnitTestSuite object
%
% NNO 2015

    disp(str(obj));

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [registered]=nonrigidICP(target,source,Ft,Fs,iterations)

% INPUT
% -target: vertices of target mesh; n*3 array of xyz coordinates
% -source: vertices of source mesh; n*3 array of xyz coordinates
% -Ft: faces of target mesh; n*3 array
% -Fs: faces of source mesh; n*3 array
% -iterations: number of iterations; usually > 100
%
% OUTPUT
% -registered: registered source vertices on target mesh. Faces are not affected and remain the same is before the registration (Fs).

%EXAMPLE

% load EXAMPLE
% [registered]=nonrigidICP(target,source,Ft,Fs,200);

tic
clf
%initial allignment and scaling
[error1,source,transform]=rigidICP(target,source,0);

%plot of the meshes
h=trisurf(Fs,source(:,1),source(:,2),source(:,3),0.3,'Edgecolor','none');
hold
light
lighting phong;
set(gca, 'visible', 'off')
set(gcf,'Color',[1 1 0.88])
view(90,90)
set(gca,'DataAspectRatio',[1 1 1],'PlotBoxAspectRatio',[1 1 1]);
trisurf(Ft,target(:,1),target(:,2),target(:,3),0.75,'Edgecolor','none');
alpha(0.6)

for i =1:iterations
    % find index of max deviating point
    [distancemax,I,error,Reallignedsource]=ICPmanu_allign2(target,source);
     distancemax=distancemax
     % define measure for decreasing surface area to be transformed
     areafactor=0.5+24*i/iterations;

    for j=1:4
        areafactor=areafactor*j;
        r=areafactor/25;

        if r>1
            r=1;
        end
        %decrease sample size for large surfaces

        [nfs,nvs] = reducepatch(Fs,source,r);
        [nft,nvt] = reducepatch(Ft,target,r);
        %define index of max deviating point on reduced surface
        Inieuw=knnsearch(nvs,source(I,:));
        % define order of neighbouring vertices
        [distancemap]=surfacemap(nvs,nfs,Inieuw);
        d2=distancemap(:,3);
        %use distance to define specific surface area of interest
        d2=d2.^areafactor;
        [error,Reallignedsource,transform]=ICPmanu2weigthed(nvt,nvs,d2);
        [distancemap]=surfacemap(source,Fs,I);
        d2=distancemap(:,3);
        d2=d2.^areafactor;
        correction=1-d2;
        Reallignedsource=transform.b*source*transform.T+repmat(transform.c(1,1:3),length(source(:,1)),1);
        source=horzcat(source(:,1).*correction,source(:,2).*correction,source(:,3).*correction)+horzcat(Reallignedsource(:,1).*d2,Reallignedsource(:,2).*d2,Reallignedsource(:,3).*d2);
        [error,source,transform]=rigidICP(target,source,1);
            delete(h)
            h=trisurf(Fs,source(:,1),source(:,2),source(:,3),d2,'Edgecolor','none');
            pause (0.1)
    end
end
registered=source;
toc

author: Sable

repo: mcbench-benchmarks


text:
function [transV, yawRotV, heightV] = visual_odometry_down(rawImg)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     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/>.

    % The simple visual odometry with scanline intensity profile algorithm.
    % the input is a raw image
    % the output including horizontal translational velocity, rotational velocity,
    % vertical translational velocity (vertical)

    %% start to set up the visual odometry

    % define the veriable for drawing sub images used by estimating translational,
    % rotational, and pitch velocity
    global SUB_TRANS_IMG;
    global SUB_YAW_ROT_IMG;
    global SUB_HEIGHT_V_IMG;

    % definition the Y (vertical) range of images for odomentry, including image for
    % translational velocity, image for rotational velocity, and image for
    % height change velocity
    global ODO_IMG_HEIGHT_V_Y_RANGE;
    global ODO_IMG_YAW_ROT_Y_RANGE;

    % definition of X (horizontal) range of the image for odometry
    global ODO_IMG_YAW_ROT_X_RANGE;
    global ODO_IMG_HEIGHT_V_X_RANGE;

    global ODO_IMG_TRANS_Y_RANGE;
    global ODO_IMG_TRANS_X_RANGE;

    % define the size of resized images for odo
    global ODO_IMG_YAW_ROT_RESIZE_RANGE;
    global ODO_IMG_HEIGHT_V_RESIZE_RANGE;
    global ODO_IMG_TRANS_RESIZE_RANGE;

    % define the scale of translational velocity, rotational velocity, and pitch velocity
    global ODO_TRANS_V_SCALE;

    global ODO_YAW_ROT_V_SCALE;
    global ODO_HEIGHT_V_SCALE;

    % difine the maximum threshold of translational velocity, rotational velocity and pitch velocity
    global MAX_TRANS_V_THRESHOLD;
    global MAX_YAW_ROT_V_THRESHOLD;
    global MAX_HEIGHT_V_THRESHOLD;

    % define the veriable for visual odometry shift match in vertical and
    % horizontal
    global ODO_SHIFT_MATCH_VERT;
    global ODO_SHIFT_MATCH_HORI;

    % define the degree of the field of view in horizontal and vertical.
    % Field of View (FOV), Degree (DEG) the horizontal, vertical, and diagonal
    % degrees for all FOVs
    global FOV_HORI_DEGREE;
    global FOV_VERT_DEGREE;

    % x_sums, sum each column of intensity in the current image and previous
    % image perspectively.
    % form one dimensional vector, 1*N array
    global PREV_TRANS_V_IMG_X_SUMS;
    global PREV_YAW_ROT_V_IMG_X_SUMS;
    global PREV_HEIGHT_V_IMG_Y_SUMS;

    % define the previous velocity for keeping stable speed
    global PREV_TRANS_V;
    global PREV_YAW_ROT_V;
    global PREV_HEIGHT_V;

    global DEGREE_TO_RADIAN;

    %%% End up for setting up the visual odometry

    global OFFSET_YAW_ROT;
    global OFFSET_HEIGHT_V;

    %% start to compute the horizontal rotational velocity (yaw)

    % get the sub_image for rotational velocity from raw image with range constrait
    subRawImg = rawImg(ODO_IMG_YAW_ROT_Y_RANGE, ODO_IMG_YAW_ROT_X_RANGE);
    subRawImg = imresize(subRawImg, ODO_IMG_YAW_ROT_RESIZE_RANGE);
    horiDegPerPixel = FOV_HORI_DEGREE / size(subRawImg, 2);

    SUB_YAW_ROT_IMG = subRawImg;
    SUB_TRANS_IMG = subRawImg;

%     % get the size of template image after resized
%     ySizeODOImg = ODO_IMG_YAW_ROT_RESIZE_RANGE(1);
%     xSizeODOImg = ODO_IMG_YAW_ROT_RESIZE_RANGE(2);
%     ySizeNormImg = ySizeODOImg;
%
%     PATCH_SIZE_Y_K = 5;
%     PATCH_SIZE_X_K = 5;
%
%     % define a temp variable for patch normalization
%     % extent the dimension of raw image for patch normalization (extODOImg, extension sub image of vtResizedImg)
%     extODOImg = zeros(ySizeODOImg + PATCH_SIZE_Y_K - 1, xSizeODOImg + PATCH_SIZE_X_K - 1);
%     extODOImg(fix((PATCH_SIZE_Y_K + 1 )/2) : fix((PATCH_SIZE_Y_K + 1 )/2) ...
%         + ySizeNormImg - 1 ,  fix((PATCH_SIZE_X_K + 1 )/2) : fix((PATCH_SIZE_X_K + 1 )/2) ...
%         + xSizeODOImg - 1 ) = subRawImg;
%
%     %% patch normalisation is applied to compensate for changes in lighting condition
%     for v = 1: ySizeNormImg
%         for u = 1 : xSizeODOImg
%             % get patch image
%             patchImg = extODOImg(v : v + PATCH_SIZE_Y_K - 1, u : u + PATCH_SIZE_X_K -1);
%
%             % Find the average of the matrix patch image
%             meanPatchImg = mean2(patchImg);
%
%             % Find the standard deviation of the matrix patch image
%             stdPatchIMG = std2(patchImg);
%
%             % Normalize the sub raw image
%             % normODOImg(v,u) = 11 * 11 * (vtResizedImg(v,u) - meanPatchImg ) / stdPatchIMG ;
%             % normODOImg(v,u) =  PATCH_SIZE_Y_K * PATCH_SIZE_X_K * (vtResizedImg(v,u) - meanPatchImg ) / stdPatchIMG ;
%             normODOImg(v,u) =  (subRawImg(v,u) - meanPatchImg ) / stdPatchIMG/ 255;
%
%         end
%     end
%
%     SUB_YAW_ROT_IMG = normODOImg;
%     SUB_TRANS_IMG = normODOImg;

    % get the x_sum of average sum intensity values in every column of image
    imgXSums = sum(subRawImg);
    avgIntensity = sum(imgXSums) / size(imgXSums, 2);
    imgXSums = imgXSums / avgIntensity;

    % compare the current image with the previous image
    % get the minimum offset and minimum difference of intensity between two images
    [minOffsetYawRot, minDiffIntensityRot] = compare_segments(imgXSums, PREV_YAW_ROT_V_IMG_X_SUMS, ODO_SHIFT_MATCH_HORI, size(imgXSums, 2));

    OFFSET_YAW_ROT = minOffsetYawRot;
    yawRotV = ODO_YAW_ROT_V_SCALE * minOffsetYawRot * horiDegPerPixel;  % in deg

    if abs(yawRotV) > MAX_YAW_ROT_V_THRESHOLD
        yawRotV = PREV_YAW_ROT_V;
    else
        PREV_YAW_ROT_V = yawRotV;
    end

    PREV_YAW_ROT_V_IMG_X_SUMS = imgXSums;
    PREV_TRANS_V_IMG_X_SUMS = imgXSums;
    %%% end up to compute the translational velocity

    %% start to compute total translational velocity

    % principle
    % speeds are estimates based on the rate of image change.
    % the speed measure v is obtained from the filtered average absolute
    % intensity difference between consecutive scanline intensity profiles at
    % the best match for rotation with best offest in yaw and pitch shift

    transV = minDiffIntensityRot * ODO_TRANS_V_SCALE;

    if transV > MAX_TRANS_V_THRESHOLD
       transV = PREV_TRANS_V;
    else
       PREV_TRANS_V = transV;
    end

    %% start to compute the height change velocity

    % get the sub_image for pitch velocity from raw image with range constrait
    subRawImg = rawImg(ODO_IMG_HEIGHT_V_Y_RANGE, ODO_IMG_HEIGHT_V_X_RANGE);
    subRawImg = imresize(subRawImg, ODO_IMG_HEIGHT_V_RESIZE_RANGE);
    vertDegPerPixel = FOV_VERT_DEGREE / size(subRawImg, 1);

    if minOffsetYawRot > 0
        subRawImg = subRawImg(:, minOffsetYawRot + 1 : end);
    else
        subRawImg = subRawImg(:, 1 : end -(-minOffsetYawRot));
    end

    SUB_HEIGHT_V_IMG = subRawImg;

    imageYSums = sum(subRawImg,2);
    avgIntensity = sum(imageYSums) / size(imageYSums, 1);
    imageYSums = imageYSums / avgIntensity;

    [minOffsetHeightV, minDiffIntensityHeight] = compare_segments(imageYSums, PREV_HEIGHT_V_IMG_Y_SUMS, ODO_SHIFT_MATCH_VERT, size(imageYSums, 1));

    if minOffsetHeightV < 0
        minDiffIntensityHeight = - minDiffIntensityHeight;
    end

    OFFSET_HEIGHT_V = minOffsetHeightV;

    % covert the perceptual speed into a physical speed with an empirically
    % determined constant TRANSLATIONAL_VELOCITY_SCALE

%     if abs(minOffsetYawRot) < 2
%         heightV = ODO_HEIGHT_V_SCALE * minDiffIntensityHeight;
%     else
%         heightV = 0;
%     end

%     if abs(minOffsetYawRot) < 2
%         heightV = minOffsetHeightV * 0.1;
%     else
%         heightV = 0;
%     end
%
%     if abs(minOffsetHeightV) < 2
%         heightV = 0;
%     end
    if minOffsetHeightV < 0
        heightV = ODO_HEIGHT_V_SCALE * minDiffIntensityHeight;
    elseif PREV_HEIGHT_V < 0
        heightV = PREV_HEIGHT_V;
    else
        heightV = 0;
    end

    if abs(heightV) > MAX_HEIGHT_V_THRESHOLD
        heightV = PREV_HEIGHT_V;
    else
        PREV_HEIGHT_V = heightV;
    end

%   if abs(minOffsetHeightV) > 1
%         heightV = ODO_HEIGHT_V_SCALE * minDiffIntensityHeight;
%
%     else
%         heightV = 0;
%   end
%
%     if abs(minOffsetHeightV) > 1 && abs(minOffsetYawRot) < 2
%         heightV = ODO_HEIGHT_V_SCALE * minDiffIntensityHeight;
%
%     else
%         heightV = 0;
%     end

    % to detect excessively large translational velocity
    % the threshold Vmax ensured that spuriously high image differences were
    % not used. Large image differences could be caused by sudden illumination
    % changes such as when travelling uphill facing directly into the sun.

    % define a maximum velocity threshold according to the average motion speed

%     if abs(minOffsetHeightV) >3
%         transV = 0;
%     end

    PREV_HEIGHT_V_IMG_Y_SUMS = imageYSums;
    %%% end up to compute the pitch velocity

end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Goal:
% Input:
%       inObjectName - File name of an input object
%       resampleFactor - Resampling factor to down- or up-sample a binary
%       object
%       saveBIM - Control flag to determine whether binary volumetric data should be
%       saved as a file for the later use. Its name will be the same as the
%       inObjectName with a different postfix('_bim.mat')
%
% Output:
%       bim - a 3D matrix of binary voxels.
%       mins - lowest coordinate values in x, y, and z axis (1 x 3)
%
% Created by Sungeun Kim (10-09-08)

function [bim, mins] = conv2bim(inObjectName, resampleFactor, saveBIM)

% Create place holders for output
bim = [];
mins = [];

if ~ischar(inObjectName)
    disp('inObjectName should a string of input filename or list of input filenames');
    return;
elseif ~exist(inObjectName, 'file')
    disp('Input object file does not exist');
    return;
end

% Parse input object filename
[path,name,ext] = fileparts(inObjectName);

if strcmp(lower(ext), '.txt')
    token = 'TRACE';
elseif strcmp(lower(ext), '.hdr') | strcmp(lower(ext), '.nii')
    token = 'NIfTI';
else
    disp('This file format can not be processed');
    return;
end

switch (token)
    case 'TRACE'
%        [bim, mins] = zroi2bim(inObjectName, resampleFactor);  %zroi2bim
%        is incomplete and need to be further modified (10-10-08)
    case 'NIfTI'
        [bim, mins] = NIfTI2bim(inObjectName, resampleFactor);
end

% Save binary volumetric object to a new file
if (saveBIM)
    new_name = [name '_bim'];
    save(fullfile(path, new_name), 'bim', 'mins');
end

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%% Generating the leading mode-n vectors
% The leading mode-n vectors are those vectors that span the subspace of
% the mode-n fibers. In other words, the left singular vectors of the
% n-mode matricization of X.
%% Using nvecs to calculate the leading mode-n vectors
% The |nvecs| command efficient computes the leading n-mode vectors.
rand('state',0);
X = sptenrand([4,3,2],6) %<-- A sparse tensor
%%
nvecs(X,1,2) %<-- The 2 leading mode-1 vectors
%%
nvecs(X,1,3) % <-- The 3 leading mode-1 vectors
%%
nvecs(full(X),1,3) %<-- The same thing for a dense tensor
%%
X = ktensor({rand(3,2),rand(3,2),rand(2,2)}) %<-- A random ktensor
%%
nvecs(X,2,1) %<-- The 1 leading mode-2 vector
%%
nvecs(full(X),2,1) %<-- Same thing for a dense tensor
%%
X = ttensor(tenrand([2,2,2,2]),{rand(3,2),rand(3,2),rand(2,2),rand(2,2)}); %<-- A random ttensor
%%
nvecs(X,4,2) %<-- The 1 leading mode-2 vector
%%
nvecs(full(X),4,2) %<-- Same thing for a dense tensor
%% Using nvecs for the HOSVD
X = tenrand([4 3 2]) %<-- Generate data
%%
U1 = nvecs(X,1,4); %<-- Mode 1
U2 = nvecs(X,2,3); %<-- Mode 2
U3 = nvecs(X,3,2); %<-- Mode 3
S = ttm(X,{pinv(U1),pinv(U2),pinv(U3)}); %<-- Core
Y = ttensor(S,{U1,U2,U3}) %<-- HOSVD of X
%%
norm(full(Y) - X) %<-- Reproduces the same result.

%%
U1 = nvecs(X,1,2); %<-- Mode 1
U2 = nvecs(X,2,2); %<-- Mode 2
U3 = nvecs(X,3,2); %<-- Mode 3
S = ttm(X,{pinv(U1),pinv(U2),pinv(U3)}); %<-- Core
Y = ttensor(S,{U1,U2,U3}) %<-- Rank-(2,2,2) HOSVD approximation of X
%%
100*(1-norm(full(Y)-X)/norm(X)) %<-- Percentage explained by approximation

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [energy] = getSpecOrbitEnergyFromSma(sma, gmu)
%getSpecOrbitEnergyFromSma Summary of this function goes here
%   Detailed explanation goes here

    energy = -gmu/(2*sma);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function planC = updateStructureMatrices(planC, editStructNum, sliceNumsV)
%"updateStructureMatrices"
%    Updates uniformized data for changes in structure contours.  editStructNum
%    is the number of the structure that has changed.  sliceNumsV is an optional
%    parameter specifying which uniformized slices should be updated.  If
%    it is not included, all relevant slices are re-uniformized.
%
%WARNING: sliceNumsV input parameter is in UNIFORMIZED slices.
%WARNING: Uniformization is performed using raster segments.  Raster
%         segments must be re-generated first!
%
% Created by   : V H Clark, 22 Aug 02.
% Last Modified: JOD, 14 Oct 02.
%              : JRA, 3 Oct 03.
%              : JRA, 10 Nov 04.
%
%Usage:
%   planC = updateStructureMatrices(planC, editStructNum, sliceNumsV)

	indexS = planC{end};

    %Figure out what relative structNum this is for the specified scanNum;
    [scanNum, relStructNum] = getStructureAssociatedScan(editStructNum, planC);

	%Get background data on scan, uniformized data assoc. with structure.
	scanInfo            = planC{indexS.scan}(scanNum).scanInfo;
	sizeArray           = getUniformScanSize(planC{indexS.scan}(scanNum));
    [xV, yV, zV]        = getUniformScanXYZVals(planC{indexS.scan}(scanNum));
    [xVs, yVs, zVs]     = getScanXYZVals(planC{indexS.scan}(scanNum));
	CTOriginalZValues   = [scanInfo(:).zValue];
	CTSliceThickness    = [scanInfo(:).sliceThickness];
	CTdeltaX            = abs(xV(2) - xV(1));
	CTdeltaY            = abs(yV(2) - yV(1));

	%Retrieve the structure's raster segments.
	[allSegmentsM, planC, isError] = getRasterSegments(editStructNum, planC);

    reusableZerosM = repmat(logical(0),[sizeArray(1), sizeArray(2)]);

    %Initialize structureArray and bits array if they dont exist.
    if length(planC{indexS.structureArray}) < scanNum;
        tmp = uint16([0 0 0]);
        tmp(1,:) = [];
        planC{indexS.structureArray}(scanNum).indicesArray = tmp;
        tmp = uint8(0);
        tmp(1) = [];
        planC{indexS.structureArray}(scanNum).bitsArray = tmp;
    end

	%If we have exceeded the max for current datatype, increase it.
	switch class(planC{indexS.structureArray}(scanNum).bitsArray)
        case 'uint8'
            if relStructNum > 8
                planC{indexS.structureArray}(scanNum).bitsArray = uint16(planC{indexS.structureArray}(scanNum).bitsArray);
            end
        case 'uint16'
            if relStructNum > 16
                planC{indexS.structureArray}(scanNum).bitsArray = uint32(planC{indexS.structureArray}(scanNum).bitsArray);
            end
        case 'uint32'
            if relStructNum > 32
                planC{indexS.structureArray}(scanNum).bitsArray = double(planC{indexS.structureArray}(scanNum).bitsArray);
            end
        case 'double'
            %Can up to uint64 here if 52 is not enough in the future.
	end

    %Clear old data for this structure.  If sliceNumsV is specified clear
    %only those slices, otherwise clear all slices.
    if exist('sliceNumsV')
        if MLVersion < 6.5
            %ML6.1 ismember cannot handle uint16.  Must double in this case.
            tf = ismember(double(planC{indexS.structureArray}(scanNum).indicesArray(:,3)), sliceNumsV);
        else
            tf = ismember(planC{indexS.structureArray}(scanNum).indicesArray(:,3), sliceNumsV);
        end
        planC{indexS.structureArray}(scanNum).bitsArray(tf) = bitset(planC{indexS.structureArray}(scanNum).bitsArray(tf), relStructNum, 0);
    else
        planC{indexS.structureArray}(scanNum).bitsArray = bitset(planC{indexS.structureArray}(scanNum).bitsArray, relStructNum, 0);
    end

	%Now we need to generate the indicesM entries for the modified structure
	k = 0;
	if ~isempty(allSegmentsM)
      wb = waitbar(0, ['Saving updated data for ' planC{indexS.structures}(editStructNum).structureName '...']);

      %If no predefined slice numbers were passed in...
      if ~exist('sliceNumsV')
          %Only use slice numbers bookending the raster segments' Zs
          minRasterZ = min(allSegmentsM(:,1));
          maxRasterZ = max(allSegmentsM(:,1));

          %Add a one slice buffer to cover all relevant uni. slices.
          minRasterZ = zVs(max(max(find(zVs <= minRasterZ)) - 1, 1));
          maxRasterZ = zVs(min(min(find(zVs >= maxRasterZ)) + 1, length(zVs)));

          %1 and length(zV) are included to prevent boundary errors.
          minUniSlice = max([find(zV < minRasterZ) 1]);
          maxUniSlice = min([find(zV > maxRasterZ) length(zV)]);

          sliceNumsV = minUniSlice:maxUniSlice;
      end

      %Preallocate the mask array.
      structM = repmat(reusableZerosM, [1 1 length(sliceNumsV)]);

      %Build the mask array, then add the mask to the uniform data.
      for zSliceUniformValue = zV(sliceNumsV)
        k = k+1;
        waitbar(k/length(sliceNumsV), wb);
        structSlc = findStructureMatrixForOneZSlice(allSegmentsM, zSliceUniformValue, CTOriginalZValues, CTSliceThickness, CTdeltaX, CTdeltaY, reusableZerosM);
        structM(:,:,k) = logical(structSlc);
      end
      planC = uniformAdd(structM, scanNum, sliceNumsV, relStructNum, planC);
	end

	if exist('wb')
      close(wb);
	end
	return

function planC = uniformAdd(mask3M, scanNum, sliceNumV, structNum, planC)
%"uniformAdd"
%   Add a single structure's slice(s) to the Uniformized Data for
%   structNum.

	indexS = planC{end};

	siz = getUniformScanSize(planC{indexS.scan}(scanNum));

	indicesM = planC{indexS.structureArray}(scanNum).indicesArray;
	bitsM    = planC{indexS.structureArray}(scanNum).bitsArray;

	%Get indices on this slice.

    if MLVersion < 6.5
        %ismember cannot handle uint16s in ML6.1.
    	ind = find(ismember(double(indicesM(:,3)), sliceNumV));
    else
        ind = find(ismember(indicesM(:,3), sliceNumV));
    end

	[r,c,s] = find3d(mask3M);
	s = sliceNumV(s);
	myDataPoints = [r',c',s'];

	%Convert from rcs to index, makes intersect and setdiff faster.
    %Do this manually since sub2ind has major overhead.
    multiM = cumprod(siz);
    mDP = double(r) + (double(c)-1)*multiM(1) + (double(s)-1)*multiM(2);
	iM  = double(indicesM(ind,1)) + (double(indicesM(ind,2))-1)*multiM(1) + (double(indicesM(ind,3))-1)*multiM(2);

	%Find voxels in both sets.
	[c, ia, ib] = intersect(mDP, iM);

	%Find points not already in Uniformized data.
	bitsM(ind) = bitset(bitsM(ind), structNum, 0);
	bitsM(ind(ib)) = bitset(bitsM(ind(ib)), structNum, 1);
	[c, i] = setdiff(mDP, c);

	newValueIndicies = i;
	numNewIndicies = length(i);

	%Add these new points.
	indicesM = [indicesM;myDataPoints(newValueIndicies, :)];

	a = 0;
	%set the bits of the new points--since they are new only this struct must
	%be flipped.
	a = bitset(a, structNum, 1);
	bitsM(length(bitsM)+1:length(indicesM)) = a;

	planC{indexS.structureArray}(scanNum).indicesArray = indicesM;
	planC{indexS.structureArray}(scanNum).bitsArray = bitsM;


text:
function C = subsref (F,ij)
%SUBSREF A(i,j) or (i,j)th entry of inv(A) if F is inverted.
% Otherwise, explicit entries in the inverse are computed.  This method
% also extracts the contents of F (A, L, U, Q, R, p, q, is_inverse, and
% kind).
%
% Example
%   F = factorize(A)
%   F(1,2)              % same as A(1,2)
%   F.L                 % the L factor of the factorization of A
%   S = inverse(A)
%   S(1,2)              % the (1,2) entry of inv(A), but only computes
%                       % the 2nd column of inv(A) via backslash.
%
% See also factorize.

% Copyright 2009, Timothy A. Davis, University of Florida

switch (ij(1).type)

    case '.'

        % F.U usage, for example: extract one of the matrices from F
        if (length (ij) > 2)
            error ('Improper index matrix reference.') ;
        end

        switch ij(1).subs
            case 'A'
                C = F.A ;
            case 'L'
                C = F.L ;
            case 'U'
                C = F.U ;
            case 'Q'
                C = F.Q ;
            case 'R'
                C = F.R ;
            case 'p'
                C = F.p ;
            case 'q'
                C = F.q ;
            case 'is_inverse'
                C = F.is_inverse ;
            case 'kind'
                C = F.kind ;
            otherwise
                error ('Reference to non-existent field ''%s''.', ...
                    ij(1).subs) ;
        end

        % F.U(2,3) usage, return U(2,3)
        if (length (ij) > 1)
            C = subsref (C, ij (2)) ;
        end

    case '()'

        % F(2,3) usage, return A(2,3) or the (2,3) entry of inv(A).
        if (length (ij) > 1)
            error ('Improper index matrix reference.') ;
        end
        A = F.A ;
        if (F.is_inverse)
            % The caller is requesting explicit entries of the inverse.
            if (length (ij.subs) ~= 2)
                error ('Linear indexing of inverse not supported.') ;
            end
            [m n] = size (A) ;
            ilen = length (ij.subs {1}) ;
            if (strcmp (ij.subs {1}, ':'))
                ilen = n ;
            end
            jlen = length (ij.subs {2}) ;
            if (strcmp (ij.subs {2}, ':'))
                jlen = m ;
            end
            j = ij ;
            j.subs {1} = ':' ;
            i = ij ;
            i.subs {2} = ':' ;
            if (jlen <= ilen)
                % For F(i,j), compute cols S(:,j) of the inverse S=inv(A)
                if (issparse (A))
                    I = speye (m) ;
                else
                    I = eye (m) ;
                end
                C = subsref (mldivide (F, subsref (I,j), 0), i) ;
            else
                % For F(i,j), compute rows S(i,:) of the inverse S=inv(A)
                if (issparse (A))
                    I = speye (n) ;
                else
                    I = eye (n) ;
                end
                C = subsref (mrdivide (subsref (I,i), F, 0), j) ;
            end
        else
            % F is not inverted, so just return A(i,j)
            C = subsref (A, ij) ;
        end

    case '{}'

        error ('Cell contents reference from a non-cell array object.') ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% ONLY FOR DEBUGGING PURPOSE (Javier)
%
%

function dispcond(cond,text,storecode, storesign, storetime, storeflag, storewrite)
% Despliega texto en el command window si cond=1
% Pensado para ocupar mensajes durante debbugin.
% Author: Javier
if nargin==2
        if cond
                disp(text)
        end
elseif nargin==7
        if cond
                disp(text)
                storecode %#ok<NOPRT>
                storesign %#ok<NOPRT>
                storetime %#ok<NOPRT>
                storeflag %#ok<NOPRT>
                storewrite %#ok<NOPRT>
        end
else
        return
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function scr = read_mat(infilename)
% Creates a Scores object from the information in a mat file.
% Inputs:
%   infilename: A mat file containing scores.
% Outputs:
%   scr: A Scores object encoding the scores in the input mat
%     file.

assert(nargin==1)
assert(isa(infilename,'char'))

load(infilename,'scr');
scr = Scores(scr);

assert(scr.validate())

author: nesl


text:
function a = rutis1 ( )

%*****************************************************************************80
%
%% RUTIS1 returns the RUTIS1 matrix.
%
%  Example:
%
%    6 4 4 1
%    4 6 1 4
%    4 1 6 4
%    1 4 4 6
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A has constant row sums.
%
%    Because it has a constant row sum of 15,
%    A has an eigenvalue of 15, and
%    a (right) eigenvector of ( 1, 1, 1, 1 ).
%
%    A has constant column sums.
%
%    Because it has a constant column sum of 15,
%    A has an eigenvalue of 15, and
%    a (left) eigenvector of ( 1, 1, 1, ..., 1 ).
%
%    A has a repeated eigenvalue.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joan Westlake,
%    A Handbook of Numerical Matrix Inversion and Solution of
%    Linear Equations,
%    John Wiley, 1968,
%    ISBN13: 978-0471936756,
%    LC: QA263.W47.
%
%  Parameters:
%
%    Output, real A(4,4), the matrix.
%

%
%  Note that the matrix entries are listed by row.
%
  a(1:4,1:4) = [ ...
    6.0,  4.0,  4.0,  1.0; ...
    4.0,  6.0,  1.0,  4.0; ...
    4.0,  1.0,  6.0,  4.0; ...
    1.0,  4.0,  4.0,  6.0 ];

  return
end


text:
function W=twoMatDiag(C1,C2,algorithm)
%%TWOMATDIAG Given two real or complex Hermitian matrices, the first of
%            which must be positive definite, find a matrix W that
%            diagonalizes both of them. Specifically, W*C1*W'=I and
%            W*C2*W'=a diagonal matrix.
%
%INPUTS: C1 An NXN real or complex positive definite Hermitian matrix.
%        C2 An NXN real or complex Hermitian matrix. This matrix does not
%           have to be positive definite.
% algorithm An optional parameter specifying the algorithm to use. Possible
%           values are:
%           0 (The default if omitted or an empty matrix is passed) Use the
%             singular value decomposition (SVD)-based algorithm given in
%             [1].
%           1 Use the eigenvalue-based decomposition given in [2]. There is
%             generally not a need to use the algorithmic variant. It is
%             not implemented in the C++/mex version of this function.
%
%OUTPUTS: W An NXN matrix such that W*C1*W'=identity matrix and W*C2*W'=a
%           diagonal matrix.
%
%Joint diagonalization of a pair of matrices arises in the fusion problem
%discussed in [1] and [2], among other applications. If both matrices are not
%positive definite, then the jointMatDiagFrob function can produce the
%desired W diagonalization matrix.
%
%EXAMPLE 1:
%Here, we diagonalize two real, positive definite matrices:
% C1=[87, 25, 18, 31;
%     25, 63, 20, 17;
%     18, 20, 65, 29;
%     31, 17, 29, 65];
% C2=[57,  7, 12, 17;
%      7, 47, 17, 22;
%     12, 17, 37, 27;
%     17, 22, 27, 27];
% %C1 is positive definite. C2 is positive semi-definite.
% W=twoMatDiag(C1,C2);
% offDiagErr1=W*C1*W'-diag(diag(W*C1*W'))
% offDiagErr2=W*C2*W'-diag(diag(W*C2*W'))
%One will see that the off-diagonal errors are on the order of 1e-16, which
%is around what one would expect with finite precision limitiations.
%
%EXAMPLE 2:
%This is the diagonalization of two complex Hermitian matrices, the first
%of which is positive definite and the second of which has some negative
%eigenvalues.
%C1 is positive definite.
% C1=[  9+  0*1i,   -65+  0*1i,  -11-153*1i,  -91-173*1i;
%     -65+  0*1i,    83+  0*1i,   54- 38*1i,   31+ 28*1i;
%     -11+153*1i,    54+ 38*1i,  130+  0*1i,   16- 47*1i;
%     -91+173*1i,    31- 28*1i,   16+ 47*1i,   22+  0*1i]+215*eye(4);
% %C2 has both positive and negative eigenvalues.
% C2=[-16+  0*1i,  -32- 56*1i,  -12-128*1i,   16+114*1i;
%     -32+ 56*1i,   79+  0*1i,  -87- 67*1i,  -48+ 51*1i;
%     -12+128*1i,  -87+ 67*1i,   76+  0*1i,   -7- 96*1i;
%      16-114*1i,  -48- 51*1i,   -7+ 96*1i, -147+  0*1i];
% W=twoMatDiag(C1,C2);
% offDiagErr=W*C1*W'-diag(diag(W*C1*W'))
% offDiagErr=W*C2*W'-diag(diag(W*C2*W'))
%One will see that the errors are on the order of 1e-14 or less, which is
%around what one would expect with finite precision errors.
%
%REFERENCES:
%[1] J. Nygårds, V. Deleskog, and G. Hendeby, "Safe fusion compared to
%    established distributed fusion methods," in IEEE International
%    Conference on Multisensor Fusion and Integration for Intelligent
%    Systems, Baden-Baden, Germany, 19-21 Sep. 2016, pp. 265-271.
%[2] M. Reinhardt, B. Noack, and U. D. Hanebeck, "Closed-form optimization
%    of covariance intersection for low-dimensional matrices," in
%    Proceedings of the 15th International Conference on Information
%    Fusion, Singapore, 9-12 Jun. 2012, pp. 1891-1896.
%
%February 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(algorithm))
    algorithm=0;
end

switch(algorithm)
    case 0
        %Use the SVD-based algorithm.
        %Note that C1=U1*D1*U1' for a valid covariance matrix.
        [U1,D1,~]=svd(C1);%Equation 6
        d1=diag(D1);

        D1Root=diag(1./sqrt(d1));

        temp=U1*D1Root;
        [U2,~,~]=svd(temp'*C2*temp);%Equation 7 in [1].

        W=U2'*D1Root*U1';%Equation 8a in [1].
        %Note that W*C1*W'=eye(xDim,xDim);
    case 1
        %Use the eigenvalue-based algorithm.
        [V1,E1]=eig(C1);
        T1=inv(V1*diag(sqrt(diag(E1))));
        C2p=T1*C2*T1';
        [V2p,~]=eig(C2p);

        %The transformation matrix in Equation 2 of [2].
        W=V2p'*diag(1./sqrt(diag(E1)))*V1';
    otherwise
        error('Unknown algorithm specified.')
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function Kc = knCenter(kn, X, X1, X2)
% Centerize the data in the kernel space
% Input:
%   kn: kernel function
%   X: d x n data matrix of which the center in the kernel space is computed
%   X1, X2: d x n1 and d x n2 data matrix. the kernel k(x1,x2) is computed
%       where the origin of the kernel space is the center of phi(X)
% Ouput:
%   Kc: n1 x n2 kernel matrix between X1 and X2 in kernel space centered by
%       center of phi(X)
% Written by Mo Chen (sth4nth@gmail.com).
K = kn(X,X);
mK = mean(K);
mmK = mean(mK);
if nargin == 2     % compute the pairwise centerized version of the kernel of X. eq knCenter(kn,X,X,X)
    Kc = K+mmK-bsxfun(@plus,mK',mK);        % Kc = K-M*K-K*M+M*K*M; where M = ones(n,n)/n;
elseif nargin == 3  % compute the norms (k(x,x)) of X1 w.r.t. the center of X as the origin. eq diag(knCenter(kn,X,X1,X1))
    Kc = kn(X1)+mmK-2*mean(kn(X,X1));
elseif nargin == 4  % compute the kernel of X1 and X2 w.r.t. the center of X as the origin
    Kc = kn(X1,X2)+mmK-bsxfun(@plus,mean(kn(X,X1))',mean(kn(X,X2)));
end

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
% getFWHM_uhe(r)
%
%   Get FWHM (mm) of the ultra-high-energy collimator PSF
%   at some distances 'r' (mm) from the collimator
%
% Qiang (Victor) Lin
% Modified by A. Yendiki, 2/5/02

function FWHM = getFWHM_uhe(r)

r(r < 0) = 0;

%  Radius (mm) FWHM (mm?)
%   20          5.6135
%   75         13.8277
%   130        18.2010
%   185        22.2908
%   245        26.9498

fw20  = 5.6135;
fw75  = 13.8277;
fw130 = 18.2010;
fw185 = 22.2908;
fw245 = 26.9498;

FWHM = (fw75 - (75-r) * ((fw75-fw20)/(75-20))) .* (r <= 75) ...
     + (fw130 - (130-r) * ((fw130-fw75)/(130-75))) .* (r <= 130 & r > 75) ...
     + (fw185 - (185-r) * ((fw185-fw130)/(185-130))) .* (r <= 185 & r > 130) ...
     + (fw185 - (185-r) * ((fw245-fw185)/(245-185))) .* (r > 185);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function itc_ds=cosmo_phase_itc(ds,varargin)
% compute phase inter trial coherence
%
% itc_ds=cosmo_phase_itc(ds,varargin)
%
% Inputs:
%   ds                      dataset struct with fields:
%       .samples            PxQ complex matrix for P samples (trials,
%                           observations) and Q features (e.g. combinations
%                           of time points, frequencies and channels)
%       .sa.targets         Px1 array with trial conditions. Each condition
%                           must occur equally often; that is, the
%                           samples must be balanced.
%                           In the typical case of two conditions,
%                           .sa.targets must have exactly two unique
%                           values.
%       .sa.chunks          Px1 array indicating which samples can be
%                           considered to be independent. It is required
%                           that all samples are independent, therefore
%                           all values in .sa.chunks must be different from
%                           each other
%       .fa                 } optional feature attributes
%       .a                  } optional sample attributes
%  'samples_are_unit_length',u  (optional, default=false)
%                           If u==true, then all elements in ds.samples
%                           are assumed to be already of unit length. If
%                           this is indeed true, this can speed up the
%                           computation of the output.
%  'check_dataset',c        (optional, default=true)
%                           if c==false, there is no check for consistency
%                           of the ds input.
%
% Output:
%   itc_ds                  dataset struct with fields:
%       .samples            (N+1)xQ array with inter-trial coherence
%                           measure, where U=unique(ds.sa.targets) and
%                           N=numel(U). The first N rows correspond to the
%                           inter trial coherence for each condition. The
%                           last row is the inter trial coherence for all
%                           samples together.
%       .sa.targets         (N+1)x1 vector containing the values
%                           [U(:);NaN]' with trial conditions
%       .a                  } if present in the input, then the output
%       .fa                 } contains these fields as well
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    defaults=struct();
    defaults.samples_are_unit_length=false;
    defaults.check_dataset=true;

    opt=cosmo_structjoin(defaults,varargin{:});

    check_input(ds,opt);

    samples=ds.samples;
    if opt.samples_are_unit_length
        quick_check_some_samples_being_unit_length(samples);
    else
        % normalize
        samples=samples./abs(samples);
    end

    % split based on .sa.targets
    [idxs,classes]=cosmo_index_unique(ds.sa.targets);
    nclasses=numel(classes);
    nfeatures=size(samples,2);

    % allocate space for output
    itc=zeros(nclasses+1,nfeatures);

    % ITC for each class
    for k=1:nclasses
        samples_k=samples(idxs{k},:);
        itc(k,:)=itc_on_unit_length_elements(samples_k);
    end

    % overall ITC
    itc(nclasses+1,:)=itc_on_unit_length_elements(samples);

    % set output
    itc_ds=set_output(itc,ds,classes);

function itc_ds=set_output(itc,ds,classes)
    % store results
    itc_ds=struct();
    itc_ds.samples=itc;
    itc_ds.sa.targets=[classes(:); NaN];

    % copy .a and .fa fields, if present
    if isfield(ds,'a')
        itc_ds.a=ds.a;

        if isfield(ds.a,'sdim')
            % remove sample dimensions if present
            itc_ds.a=rmfield(itc_ds.a,'sdim');
        end
    end

    if isfield(ds,'fa')
        itc_ds.fa=ds.fa;
    end

function itc=itc_on_unit_length_elements(samples)
    % computes inter-trial coherence for each column seperately
    itc=abs(sum(samples,1)/size(samples,1));

function quick_check_some_samples_being_unit_length(samples)
    % instead of checking all values, only verify for a subset of values.
    % This should prevent most use cases where the user accidentally
    % uses non-normalized data, whereas checking all values would be
    % equivalent to actually computing their length for each of them.
    count_to_check=10;

    % generate random positions to check for unit length
    nelem=numel(samples);
    pos=ceil(rand(1,count_to_check)*nelem);

    samples_subset=samples(pos);
    lengths=abs(samples_subset);

    % safety margin
    delta=10*eps('single');
    if any(lengths+delta<1 | lengths-delta>1)
        error('.samples input is not of unit length');
    end

function check_input(ds,opt)
    % must be a proper dataset
    if opt.check_dataset
        raise_exception=true;
        cosmo_check_dataset(ds,raise_exception);

        % must have targets and chunks
        cosmo_isfield(ds,{'sa.targets','sa.chunks'},raise_exception);
    end

    % all chunks must be unique
    if ~isequal(sort(ds.sa.chunks),unique(ds.sa.chunks))
        error(['All values in .sa.chunks must be different '...
                    'from each other']);
    end

    % trial counts must be balanced
    [idxs,classes]=cosmo_index_unique(ds.sa.targets);
    class_count=cellfun(@numel,idxs);
    unequal_pos=find(class_count~=class_count(1),1);
    if ~isempty(unequal_pos)
        error(['.sa.targets indicates unbalanced targets, with '...
                        '.sa.targets==%d occuring %d times, and '...
                        '.sa.targets==%d occuring %d times.\n'...
                        'To obtain balanced targets, consider '...
                        'using cosmo_balance_dataset.'],...
                    classes(1),class_count(1),...
                    classes(unequal_pos),class_count(unequal_pos));
    end

    % input must be complex
    if isreal(ds.samples)
        error('.samples must be complex');
    end

    v=opt.samples_are_unit_length;
    if ~(islogical(v) ...
            && isscalar(v))
        error('option ''samples_are_unit_length'' must be logical scalar');
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [states] = tapas_sem_multiv_init_states(data, model, inference)
%% Generates the data structure to store the samples.
%
% Input
%
%
% Output
%   data        -- Cell array with nx1 cells, where n = niter
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%

states = cell(inference.niter, 1);

end

author: translationalneuromodeling

repo: tapas


text:
function [ dmap ] = dataDensity( x, y, width, height, limits, fudge )
%DATADENSITY Get a data density image of data
%   x, y - two vectors of equal length giving scatterplot x, y co-ords
%   width, height - dimensions of the data density plot, in pixels
%   limits - [xmin xmax ymin ymax] - defaults to data max/min
%   fudge - the amount of smear, defaults to size of pixel diagonal
%
% By Malcolm McLean
%
    if(nargin == 4)
        limits(1) = min(x);
        limits(2) = max(x);
        limits(3) = min(y);
        limits(4) = max(y);
    end
    deltax = (limits(2) - limits(1)) / width;
    deltay = (limits(4) - limits(3)) / height;
    if(nargin < 6)
        fudge = sqrt(deltax^2 + deltay^2);
    end
    dmap = zeros(height, width);
    for ii = 0: height - 1
        yi = limits(3) + ii * deltay + deltay/2;
        for jj = 0 : width - 1
            xi = limits(1) + jj * deltax + deltax/2;
            dd = 0;
            for kk = 1: length(x)
                dist2 = (x(kk) - xi)^2 + (y(kk) - yi)^2;
                dd = dd + 1 / ( dist2 + fudge);
            end
            dmap(ii+1,jj+1) = dd;
        end
    end

end


text:
function d = Point2LineDistance(P, A, B)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Annibale Panichella

    d = zeros(size(P,1),1);
    for i = 1:size(P,1)
        pa = P(i,:) - A;
        ba = B - A;
        t = dot(pa, ba)/dot(ba, ba);
        d(i,1) = norm(pa - t * ba,2);
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
clear;
clc
close all;

%% EXAMPLE 2.2(a) Conversion of Curvilinear Position to Cartesian ECEF Position

lat = deg2rad(45);
lon = deg2rad(30);
hgt = 1000;

[R_meridian, R_transverse, ~, ~, ~]= ch_earth(lat, lon, hgt);

xyz = ch_lla2ecef(lat, lon, hgt);
xyz

%% EXAMPLE 2.2(b) Conversion of Cartesian ECEF Position to Curvilinear Position
xyz = [3912960.837 2259148.993 4488055.516]';
[lat lon hgt] = ch_ecef2lla(xyz);
lat
lon
hgt

author: yandld


text:
%helpme_it    iterative solvers interactive help
%   IFISS scriptfile: HCE; 23 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

fprintf(' \n');
fprintf('  To use a preconditioned Krylov subspace iterative method\n');
fprintf('  to solve the current discrete system, run the script file: \n');
fprintf('  it_solve\n');
fprintf(' \n');

author: johannesgerer

repo: jburkardt-m


text:
function ebsd = calcParentEBSD(job)
% reconstruct parent EBSD
%
% Syntax
%   job.calcParentEBSD
%
% Input
%  job - @parentGrainReconstructor
%
% Output
%  ebsd - reconstructed parent @EBSD
%

% copy prior EBSD
ebsd = job.ebsdPrior;

% find Ids of grains that are either parent or child phase
isRecData = ebsd.phaseId == job.childPhaseId | ...
  ebsd.phaseId == job.parentPhaseId;

% compute new grainIds
ebsd.grainId(ebsd.grainId>0) = job.mergeId(ebsd.grainId(ebsd.grainId>0));

grainIds = max(1,ebsd.grainId);
grainIds(~isRecData) = 1;

% consider only child pixels that have been reconstructed to parent
% grains
isNowParent = ebsd.phaseId == job.childPhaseId &...
  job.grains.phaseId(grainIds) == job.parentPhaseId;

% maybe there is nothing to do
if nnz(isNowParent) == 0, return; end

% compute parent orientation
[ori,fit] = calcParent(ebsd(isNowParent).orientations,...
  job.grains(ebsd.grainId(isNowParent)).meanOrientation,job.p2c);

% compute variantId
vId = calcVariantId(ori, ebsd(isNowParent).orientations,job.p2c);
ebsd.prop.variantId = NaN(size(ebsd));
ebsd.prop.variantId(isNowParent) = vId;

% adjust parent and child orientations such that the misorientation is
% closest to the given OR job.p2c
ori = ori.project2FundamentalRegion .* inv(variants(job.p2c,vId)) * job.p2c;

% setup parent ebsd
ebsd.prop.fit = nan(size(ebsd));
ebsd(isNowParent).orientations = ori;
ebsd.prop.fit(isNowParent) = fit;

end

author: mtex-toolbox

repo: mtex

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


text:
% HMPD: Estimate the short-term deviation of Phase Distortion
%
% Inputs
%  PD   : [NxM rad] A matrix of Phase Distortion to measure the deviation from.
%         N is the number of frames, M is the order of the PD (either the
%         maximum number of harmonics or the number of bins).
%  nbat : The number of frames to consider in the smoothing window, i.e. the
%         window size.
%
% Outputs
%  PDD  : [NxM rad] The Phase Distortion Deviation (PDD)
%
% Copyright (c) 2013 University of Crete - Computer Science Department(UOC-CSD)/
%                    Foundation for Research and Technology-Hellas - Institute
%                    of Computer Science (FORTH-ICS)
%
% License
%  This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function PDD = hmpd_phase_deviation(PD, nbat)

    winlen = round(nbat/2)*2+1;
    % win = hann(winlen);
    win = ones(winlen,1); % Rectangular window is better than smooth window
                          % It better discriminates voiced/unvoiced segments
    win = win./sum(win);

    % Compute the std in polar coordinates
    % Compute first the center of gravity of the exp(i*angles)
    PDc = filtfilt(win, 1, cos(PD));
    PDs = filtfilt(win, 1, sin(PD));

    z = abs(PDc + 1j*PDs); % For the variance, we need only the magnitude

    PDD = zeros(size(PDc));
    idx = find(z<1); % To avoid neg sqrt or sqrt(-log(0))
    PDD(idx) = sqrt(-2*log(z(idx))); % Fisher's standard-deviation

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
% khobe ama ba hplot mishkel dare c_type kar nemikone Ti Td , Kc ro ehtemalan baid global bedam
function pid
    global C_controller
    hfigure = findall(0,'tag','pid');
    if hfigure
        figure(hfigure)
    else
        h_pid = figure('unit','pixel','pos',[172  412  450  102],'ToolBar','none',...
            'menubar','none','resize','off','CloseRequestFcn','','tag','pid','name','controler',...
            'NumberTitle','off');
        h_panel = uipanel('parent',h_pid,'unit','pixel','pos',[9    14   432    76]);
        h_Ti = uicontrol('parent',h_panel,'unit','pixel','pos',[250.0000   53.0000   77.0000   16.0000],...
            'style','text','string','1/Ti','tag','text2');
        h_Td = uicontrol('parent',h_panel,'unit','pixel','pos',[350.0000   53.0000   77.0000   16.0000],...
            'style','text','string','Td','tag','text3');
        h_P = uicontrol('parent',h_panel,'unit','pixel','pos',[150.0000   53.0000   77.0000   16.0000],...
            'style','text','string','P','tag','text1');
        h_Ti_edit = uicontrol('parent',h_panel,'unit','pixel','pos',[250.0000   32.0000   76.0000   21.0000],...
            'style','edit','string','1','tag','edit2');
        h_Td_edit = uicontrol('parent',h_panel,'unit','pixel','pos',[350.0000   32.0000   76.0000   21.0000],...
            'style','edit','string','1','tag','edit3');
        h_P_edit = uicontrol('parent',h_panel,'unit','pixel','pos',[150.0000   32.0000   76.0000   21.0000],...
            'style','edit','string','1','tag','edit1');
        h_popupmenu = uicontrol('parent',h_panel,'unit','pixel','pos',[13    30    88    21],...
            'style','popupmenu','string', {'P' 'P.I' 'P.I.D'},'BackgroundColor',[1 1 1],...
            'callback',@value_callback,'tag','popup');
        txths = findall(0,'style','text','-not','tag','text1');
        set(txths,'vis','off')
        txthe = findall(0,'style','edit','-not','tag','edit1');
        set(txthe,'vis','off')
        C_controller=1;
    end
    function value_callback(hObject, eventdata, handles)
        global C_controller
        h_popupmenu = findall(0,'tag','popup');
        hpid= get(h_popupmenu,'value');
        switch hpid
            case 1
                txths = findall(0,'style','text','-not','tag','text1');
                set(txths,'vis','off')
                txthe = findall(0,'style','edit','-not','tag','edit1');
                set(txthe,'vis','off')
                C_controller=1;
            case 2
                txths = findall(0,'style','text');
                set(txths,'vis','on')
                txthe = findall(0,'style','edit');
                set(txthe,'vis','on')
                txths = findall(0,'tag','text3');
                set(txths,'vis','off')
                txthe = findall(0,'tag','edit3');
                set(txthe,'vis','off')
                C_controller=2;
            case 3
                txths = findall(0,'style','text');
                set(txths,'vis','on')
                txthe = findall(0,'style','edit');
                set(txthe,'vis','on')
                C_controller=3;
        end

author: Sable


text:
function engine = jtree_2TBN_inf_engine(bnet, varargin)
% JTREE_ONLINE_INF_ENGINE Online Junction tree inference algorithm for DBNs.
% engine = jtree_online_inf_engine(bnet, ...)
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default value in brackets]
%
% clusters - specifies variables that must be grouped in the 1.5 slice DBN
% maximize - 1 means do max-product, 0 means sum-product [0]
%
% The same nodes must be observed in every slice.

ss = length(bnet.intra);
clusters = {};
engine.maximize = 0;

args = varargin;
nargs = length(args);
for i=1:2:length(args)
  switch args{i},
   case 'clusters', clusters = args{i+1};
   case 'maximize', engine.maximize = args{i+1};
   otherwise, error(['unrecognized argument ' args{i}])
  end
end

engine.evidence = [];
engine.node_sizes = [];

%int = compute_interface_nodes(bnet.intra, bnet.inter);
int = [];

if 1
% include nodes with any outgoing arcs
for u=1:ss
  if any(bnet.inter(u,:))
    int = [int u];
  end
end
end

if 0
% include nodes with any incoming  arcs
incoming = [];
for u=1:ss
  if any(bnet.inter(:,u))
    int = [int u];
    incoming = [incoming u];
  end
end
% include nodes which are parents of nodes with incoming
for u=1:ss
  cs = children(bnet.intra, u);
  if ~isempty(cs) & mysubset(cs, incoming)
    int = [int u];
  end
end
int = unique(int);
end % if

int
engine.interface = int;
engine.nonint = mysetdiff(1:ss, int);

onodes = bnet.observed;

% Create a "1.5 slice" jtree, containing the interface nodes of slice 1
% and all the nodes of slice 2
% To keep the node numbering the same, we simply disconnect the non-interface nodes
% from slice 1, and set their size to 1.
% We do this to speed things up, and so that the likelihood is computed correctly - we do not need to do
% this if we just want to compute marginals (i.e., we can include nodes whose potentials will
% be left as all 1s).
intra15 = bnet.intra;
for i=engine.nonint(:)'
  intra15(:,i) = 0;
  intra15(i,:) = 0;
end
dag15 = [intra15      bnet.inter;
	 zeros(ss)    bnet.intra];
ns = bnet.node_sizes(:);
%ns(engine.nonint) = 1; % disconnected nodes get size 1
obs_nodes = [onodes(:) onodes(:)+ss];
bnet15 = mk_bnet(dag15, ns, 'discrete', bnet.dnodes, 'equiv_class', bnet.equiv_class(:), ...
		 'observed', obs_nodes(:));

% use unconstrained elimination,
% but force there to be a clique containing both interfaces
clusters(end+1:end+2) = {int, int+ss};
engine.jtree_engine = jtree_inf_engine(bnet15, 'clusters', clusters, 'root', int+ss);
jtree_engine = struct(engine.jtree_engine); % violate object privacy

engine.in_clq = clq_containing_nodes(engine.jtree_engine, int);
engine.out_clq = clq_containing_nodes(engine.jtree_engine, int+ss);
engine.clq_ass_to_node = jtree_engine.clq_ass_to_node;
engine.root = jtree_engine.root_clq;

% Also create an engine just for slice 1
bnet1 = mk_bnet(bnet.intra1, bnet.node_sizes_slice, 'discrete', myintersect(bnet.dnodes,1:ss), ...
		'equiv_class', bnet.equiv_class(:,1), 'observed', onodes);
for i=1:max(bnet1.equiv_class)
  bnet1.CPD{i} = bnet.CPD{i};
end
engine.jtree_engine1 = jtree_inf_engine(bnet1, 'clusters', {int}, 'root', int);
jtree_engine1 = struct(engine.jtree_engine1); % violate object privacy
engine.int_clq1 = clq_containing_nodes(engine.jtree_engine1, int);
engine.clq_ass_to_node1 = jtree_engine1.clq_ass_to_node;
engine.root1 = jtree_engine1.root_clq;

engine.observed = [onodes onodes+ss];
engine.observed1 = onodes;
engine.pot_type = determine_pot_type(bnet, onodes);
engine.slice_size = bnet.nnodes_per_slice;

engine = class(engine, 'jtree_2TBN_inf_engine', inf_engine(bnet));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%-----------------------------------------------------------------------
% Job configuration created by cfg_util (rev $Rev: 4899 $)
%-----------------------------------------------------------------------
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_files.type = 'cfg_files';
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_files.name = 'Directory';
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_files.tag = 'dir';
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_files.filter = 'dir';
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_files.ufilter = '.*';
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_files.dir = '';
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_files.num = [1 1];
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_files.check = [];
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_files.help = {'New working directory.'};
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_files.def = [];
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.type = 'cfg_exbranch';
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.name = 'Change Directory';
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.tag = 'cfg_cd';
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1) = cfg_dep;
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1).tname = 'Val Item';
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1).tgt_spec = {};
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1).sname = 'Directory (cfg_files)';
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1).src_exbranch = substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1).src_output = substruct('()',{1});
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.prog = @cfg_run_cd;
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.vout = [];
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.check = [];
matlabbatch{2}.menu_cfg{1}.menu_struct{1}.conf_exbranch.help = {'Change working directory.'};

author: spm


text:
function obj = downloadLatestGithubVersion(obj)
	% Biafra Ahanonu
	% Started: 2021.03.25 [22:11:25] (branched from ciatah.m)

	% Blank function
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
%% Copyright (C) 2016-2017 Lagu
%% Copyright (C) 2017, 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym ellipticCK (@var{m})
%% Complementary complete elliptic integral of the first kind.
%%
%% The complete elliptic integral (of the first kind) with the
%% complementary parameter @code{1 - @var{m}} is given by:
%% @example
%% @group
%% syms m
%% ellipticCK (m)
%%   @result{} ans = (sym) K(1 - m)
%% @end group
%% @end example
%%
%% Example:
%% @example
%% @group
%% ellipticCK (sym (1)/4)
%%   @result{} ans = (sym) K(3/4)
%% vpa (ans)
%%   @result{} (sym) 2.1565156474996432354386749988003
%% @end group
%% @end example
%%
%% There are other conventions for the inputs of elliptic integrals,
%% @pxref{@@sym/ellipticF}.
%%
%% @seealso{@@sym/ellipticK}
%% @end defmethod

function y = ellipticCK (m)
  if (nargin > 1)
    print_usage ();
  end

  y = ellipticK (1 - m);

end

%!error ellipticCK (sym (1), 2)

%!assert (double (ellipticCK (sym (1)/2)), 1.8541, 10e-5)
%!assert (double (ellipticCK (sym (101)/10)), 0.812691836806976, -3*eps)
%!assert (isequal (ellipticCK (sym (1)), sym(pi)/2))

repo: octsympy


text:
% function spm_dcm_Granger_Jacobian
% Demo routine for induced responses
%==========================================================================
%
% This routine illustrates the derivation of spectral Granger causal
% measures from the inversion of a simple state-space DCM paramterised in
% explcicity in terms of its Jacobian (null model).
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_Granger_Jacobian.m 5908 2014-03-05 20:31:57Z karl $

% Process (model) specification
%==========================================================================
rng('default')

% number of regions
%--------------------------------------------------------------------------
Nc    = 2;                                       % number of channels
Ns    = 2;                                       % number of sources
ns    = 2*128;                                   % sampling frequency
dt    = 1/ns;                                    % time bins
Hz    = 1:128;                                   % frequency
p     = 16;                                      % autoregression order
options.spatial  = 'LFP';
options.model    = 'CMC';
options.analysis = 'CSD';
M.dipfit.model = options.model;
M.dipfit.type  = options.spatial;
M.dipfit.Nc    = Nc;
M.dipfit.Ns    = Ns;
M.pF.D         = [1 4];

% extrinsic connections (forward an backward)
%--------------------------------------------------------------------------
A{1} = [0 0; 0 0];
A{2} = [0 0; 0 0];
A{3} = [0 0; 0 0];
B    = {};
C    = sparse(2,0);

% get priors
%--------------------------------------------------------------------------
pE    = spm_dcm_neural_priors(A,B,C,options.model);
pE    = spm_L_priors(M.dipfit,pE);
pE    = spm_ssr_priors(pE);
[x,f] = spm_dcm_x_neural(pE,options.model);

% create forward model
%--------------------------------------------------------------------------
M.f   = f;
M.g   = 'spm_gx_erp';
M.x   = x;
M.n   = length(spm_vec(x));
M.pE  = pE;
M.m   = Ns;
M.l   = Nc;
M.Hz  = Hz;
M.Rft = 4;

% specify M.u - endogenous input (fluctuations) and intial states
%--------------------------------------------------------------------------
M.u   = sparse(Ns,1);
M.x   = spm_dcm_neural_x(pE,M);

% get expected CSD (observations)
%==========================================================================

% (log) connectivity parameters (forward connection only)
%--------------------------------------------------------------------------
pE.A{1}(2,1) = 2;
pE.S         = 1/8;

% (log) amplitude of fluctations and noise
%--------------------------------------------------------------------------
pE.a(1,:) = -2;
pE.b(1,:) = -4;
pE.c(1,:) = -4;

% expected cross spectral density
%--------------------------------------------------------------------------
[csd,Hz,mtf] = spm_csd_mtf(pE,M);

% place in data structure
%--------------------------------------------------------------------------
DCM.xY.y  = csd;
DCM.xY.dt = dt;
DCM.xY.Hz = Hz;

% Invert using a null model
%==========================================================================
DCM.options.model   = 'NULL';
DCM.options.spatial = 'LFP';
DCM.options.DATA    = 0;

DCM.M.dipfit.Nc  = Nc;
DCM.M.dipfit.Ns  = Ns;
DCM.M.U  = eye(Nc,Nc);

DCM.A = {ones(Nc,Nc)};
DCM.B = {};
DCM.C = sparse(Nc,0);

% estimate
%--------------------------------------------------------------------------
DCM   = spm_dcm_csd(DCM);

% show results in terms of transfer functions and Granger causality
%==========================================================================
spm_figure('GetWin','Figure 1'); clf

% transfer functions in the absence of measurement noise
%--------------------------------------------------------------------------
Ep    = DCM.Ep;
Ep.b  = Ep.b - 32;              % and suppress non-specific and
Ep.c  = Ep.c - 32;              % specific channel noise

Gu           = spm_csd_mtf_gu(pE,Hz);
[psd Hz dtf] = spm_csd_mtf(Ep,DCM.M);

mtf   = mtf{1};
csd   = csd{1};
psd   = psd{1};

tew   = spm_dtf2gew(mtf,Gu);
ccf   = spm_csd2ccf(psd,Hz,dt);
qew   = spm_ccf2gew(ccf,Hz,dt,p);
ccf   = spm_csd2ccf(csd,Hz,dt);
gew   = spm_ccf2gew(ccf,Hz,dt,p);

spm_spectral_plot(Hz,tew, 'b', 'frequency','density')
spm_spectral_plot(Hz,qew, 'r',  'frequency','density')
spm_spectral_plot(Hz,gew, 'g',  'frequency','density')
legend('Granger causality (true)',...
       'Granger causality (source)',...
       'Granger causality (channel)')

subplot(2,2,3), a = axis;
subplot(2,2,1), axis(a);
subplot(2,2,2), axis(a);
subplot(2,2,4), axis(a);


text:
function [ G ] = gsp_nn_hypergraph( Xin, param )
%GSP_NN_HYPERGRAPH Create a nearest neighbors hypergraph from a point cloud
%   Usage :  G = gsp_nn_hypergraph( Xin );
%            G = gsp_nn_hypergraph( Xin, param );
%
%   Input parameters:
%       Xin         : Input points
%       param       : Structure of optional parameters
%
%   Output parameters:
%       G           : Resulting graph
%
%   Example:::
%
%           P = rand(100,2);
%           G = gsp_nn_hypergraph(P)
%
%   Additional parameters
%   ---------------------
%
%   * *param.use_flann* : [0, 1]              use the FLANN library
%   * *param.center*    : [0, 1]              center the data
%   * *param.rescale*   : [0, 1]              rescale the data (in a 1-ball)
%   * *param.sigma*     : float               the variance of the distance kernel
%   * *param.k*         : int                 number of neighbors for knn
%
%   See also: gsp_nn_graph
%

% Author: Nathanael Perraudin
% Date: 21 October 2015
% Testing: test_rmse

%   * *param.type*      : ['knn', 'radius']   the type of graph (default 'knn')
%   * *param.epsilon*   : float               the radius for the range search
%   * *param.use_l1*    : [0, 1]              use the l1 distance

    if nargin < 2
    % Define parameters
        param = {};
    end

    %Parameters
%     if ~isfield(param, 'type'), param.type = 'knn'; end
    if ~isfield(param, 'use_flann'), param.use_flann = 0; end
    if ~isfield(param, 'center'), param.center = 0; end
    if ~isfield(param, 'rescale'), param.rescale = 1; end
    if ~isfield(param, 'k'), param.k = 10; end

    param.type = 'knn';
%     if ~isfield(param, 'epsilon'), param.epsilon = 0.01; end
%     if ~isfield(param, 'use_l1'), param.use_l1 = 0; end
%     if ~isfield(param, 'target_degree'), param.target_degree = 0; end;
    paramnn = param;
%     paramnn.k = param.k +1;
    [indx, ~, dist] = gsp_nn_distanz(Xin',Xin',paramnn);

%     switch param.type
%         case 'knn'
%             if param.use_l1
%                 if ~isfield(param, 'sigma'), param.sigma = mean(dist); end
%             else
                if ~isfield(param, 'sigma'), param.sigma = mean(dist)^2; end
%             end
%         case 'radius'
%             if param.use_l1
%                 if ~isfield(param, 'sigma'), param.sigma = epsilon/2; end
%             else
%                 if ~isfield(param, 'sigma'), param.sigma = epsilon.^2/2; end
%             end
%         otherwise
%             error('Unknown graph type')
%     end

    w = exp(-dist.^2/param.sigma);
    G.N = size(Xin,1);
    G.Ne = G.N;
    G.W = sparse(G.N,G.Ne);
    G.E = cell(G.Ne,1);
    k = param.k;
    for ii = 1:G.Ne
        edge = indx((1:k)+(ii-1)*k);
        G.E{ii} = edge;
        % Here we use H for HW...
        G.W(edge,ii) = sqrt(sum(w(edge)));
    end
    G.hypergraph = 1;
    G.directed = 0;

    %Fill in the graph structure
    G.coords = Xin;

    G.type = 'Nearest neighboors hypergraph';
    G.lap_type = 'normalized';
    G.sigma = param.sigma;
    G = gsp_graph_default_parameters(G);
end

repo: gspbox