text:
function [fFactorHi, fStdHi, fResHi, fPerHi, fFactorLow, fStdLow, fResLow, fPerLow, fMRateFit] = calc_ratefac(mCat1, mCat2, fPeriod1, fPeriod2, fMc1, fMc2)
% [fFactorHi, fStdHi, fResHi, fPerHi, fFactorLow, fStdLow, fResLow, fPerLow, fMRateFit] = calc_ratefac(mCat1, mCat2, fPeriod1, fPeriod2, fMc1, fMc2)
% --------------------------------------------------------------------------------------------------------------------------------------------------
% Determine rate factor between amount of events between two time periods, differentiated for EQ with M<Mc and M>=Mc.
% Mc = min(fMc1, fMc2);
% N(Per2) = fFac * N (Per1)
% fFac is determined as mean of the factors in the magnitude bins, using a binning of 0.1.A factor for
% magnitudes M<Mc and M>=Mc is determined to explore a difference.Magnitude bins with no earthquakes
% are not used.
%
% Incoming variables:
% mCat1 : EQ catalog period 1
% mCat2 : EQ catalog period 2
% fPeriod1 : Time length period 1 /[dec. year]
% fPeriod2 : Time length period 2 /[dec. year]
% fMc1 : Magnitude of completeness period 1
% fMc2 : Magnitude of completeness period 2
%
% Outgoing variable:
% fFactorHi : Mean rate factor for EQ with M>=Mc
% fStdLow   : Standard deviation of fFactorHi
% fResLow   : Residuum for model application
% fPerHi    : Total percentage of change for EQ with M>=Mc
% fFactorLow : Mean rate factor for EQ with M<Mc
% fStdHi   : Standard deviation offFactorLow
% fResHi   : Residuum for model application
% fPerLow  : Total percentage of change for EQ with M<Mc
% fMRateFit: Goodness of fit by modeling second period with rate factor multiplication of first period
%
% J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 09.10.02

% Track of changes:

% Initialization
vMaxMag = [max(mCat1(:,6)) max(mCat2(:,6))];
fMaxMag =  max(vMaxMag);
fMc = min([fMc1 fMc2]);

% Binning
%% EQ with M < Mc
vSel1 = (mCat1(:,6) < fMc);
vSel2 = (mCat2(:,6) < fMc);
mCat1Mc = mCat1(vSel1,:);
mCat2Mc = mCat2(vSel2,:);
[vMag1,vBin1]=histogram(mCat1Mc(:,6),0:0.1:fMaxMag);
[vMag1sort,vBin1sort]=sort(vMag1);
[vMag2,vBin2]=histogram(mCat2Mc(:,6),0:0.1:fMaxMag);
[vMag2sort,vBin2sort]=sort(vMag2);

%% EQ with M >= Mc
mCat1McHi = mCat1(~vSel1,:);
mCat2McHi = mCat2(~vSel2,:);
[vMag1Hi,vBin1Hi]=histogram(mCat1McHi(:,6),0:0.1:fMaxMag);
[vMag1Hisort,vBin1Hisort]=sort(vMag1Hi);
[vMag2Hi,vBin2Hi]=histogram(mCat2McHi(:,6),0:0.1:fMaxMag);
[vMag2Hisort,vBin2Hisort]=sort(vMag2Hi);

%% Overall percentage of change of EQs
fPerLow = length(mCat2Mc(:,1))/length(mCat1Mc(:,1));
fPerHi = length(mCat2McHi(:,1))/length(mCat1McHi(:,1));
fPerLow = 100*fPerLow-100;
fPerHi = 100*fPerHi-100;

% Estimate factor between rates
% Factor for values below Mc of Period 1
vMag1 = vMag1/fPeriod1; % Time normalization
vMag2 = vMag2/fPeriod2; % Time normalization
fRatioLow = vMag2./vMag1;
%vSelRat = 1 - (isnan(fRatioLow) + isinf(fRatioLow));
vSelRatLow = (~isnan(fRatioLow) & ~isinf(fRatioLow));
fFactorLow = mean(fRatioLow(vSelRatLow));
fStdLow = std(fRatioLow(vSelRatLow));
vMagModelow = vMag1*fFactorLow;
fResLow = sqrt(sum((vMag2-(vMagModelow)).^2)/length(vMag2)); % Residuum

% Factor for values above or equal Mc of background
vMag1Hi = vMag1Hi/fPeriod1; % Time normalization
vMag2Hi = vMag2Hi/fPeriod2; % Time normalization
fRatioHi = vMag2Hi./vMag1Hi;
%vSelRatHi = 1 - (isnan(fRatioHi) + isinf(fRatioHi));
vSelRatHi = (~isnan(fRatioHi) & ~isinf(fRatioHi));
fFactorHi = mean(fRatioHi(vSelRatHi));
fStdHi = std(fRatioHi(vSelRatHi));
% vSelStd = (fRatioHi <= (fFactorHi+fStdHi) & fRatioHi >= (fFactorHi-fStdHi));
% fFactorHi = mean(fRatioHi(vSelStd))
vMagModel = vMag1Hi*fFactorHi;
fResHi = sqrt(sum((vMag2Hi-(vMagModel)).^2)/length(vMag2Hi)); % Residuum

% Determine goodness fit of modeled activity or M>=Mc
vEventsum2 = cumsum(vMag2Hi);
vEventsumMod = cumsum(vMagModel);
fMRateFit  = sum(abs(vEventsum2-vEventsumMod))/sum(vEventsum2);
fMRateFit = 100-fMRateFit*100;


text:
function planar = spm_eeg_planarchannelset(data)
% Define the planar gradiometer channel combinations
% FORMAT planar = spm_eeg_planarchannelset(data)
%
% The output cell array contains the horizontal label, vertical label and
% the label after combining the two.
%__________________________________________________________________________
% Copyright (C) 2013-2015 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_planarchannelset.m 7445 2018-10-12 13:24:48Z vladimir $

if isa(data, 'cell') && ~isempty(data) && isa(data{1}, 'char')
    input.label = data;
else
    input = data;
end

senstype = lower(ft_senstype(input));

k = strfind(senstype, '_combined');
if ~isempty(k)
    senstype(k:(k+length('_combined')-1)) = [];
end

try
    planar = ft_senslabel(senstype, 'output', 'planarcombined');
catch
    planar = ft_senslabel([senstype '_planar'], 'output', 'planarcombined');
end

author: spm


text:
function fx1 = p01_fx1 ( x )

%*****************************************************************************80
%
%% P01_FX1 evaluates the derivative of the function for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX1, the first derivative of the function at X.
%
  fx1 = cos ( x ) - 0.5;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fieldName = GetFieldNameFromFilePath(filepath, nStage)
if nargin<2
    nStage = 1;
end

orig_filepath = filepath;
[filepath, uttID] = fileparts(filepath);
fieldName = uttID;
for i=2:nStage
    [filepath, uttID] = fileparts(filepath);
    fieldName = [uttID '_' fieldName];
    if isempty(filepath)
        break;
    end
end
fieldName = ['U_' fieldName];
fieldName = regexprep(fieldName, '-', '_');

end

author: singaxiong

repo: SignalGraph


text:
function c = char(v,varargin)
% convert to char
%
% Flags
%  LATEX -
%

% short summary for long vectors
if length(v) > 4

  c = [size2str(v), ' points'];
  if ~check_option(varargin,'short') && v.resolution<2*pi
    c = [c, ', res.: ',xnum2str(v.resolution * 180/pi),mtexdegchar];
  end

  return
end

% list all elements

if max(abs(v.x)) < 1e-14, v.x = zeros(size(v.x));end
if max(abs(v.y)) < 1e-14, v.y = zeros(size(v.y));end
if max(abs(v.z)) < 1e-14, v.z = zeros(size(v.z));end

c = [];
for i = 1:length(v.x)
  if check_option(varargin,{'LATEX','tex'})
    if v == xvector
      c = [c,' x'];
    elseif v == yvector
        c = [c,' y'];
    elseif v == zvector
      c = [c,' z'];
    else
      iv = vec2int([v.x(i),v.y(i),v.z(i)]);
      if ~isempty(iv)
        c = [c,' ',barchar(iv(1),varargin{:}),...
          barchar(iv(2),varargin{:}),...
          barchar(iv(3),varargin{:})];
      else
        c = [c,' ',num2str([v.x(i),v.y(i),v.z(i)],'(%3.2f,%3.2f,%3.2f)')];
      end
    end
  else
    c = [c,' ',num2str(v.x(i)),',',num2str(v.y(i)),',',num2str(v.z(i))]; %#ok<AGROW>
  end
end

if ~isempty(c), c(1)=[];end
if ~isempty(c) && check_option(varargin,{'LaTeX'}), c = ['$' c '$'];end

function iv = vec2int(v)

% find common divisor
nz = find(abs(v)==max(abs(v)),1,'first');

for i = 1:9
  iv = v / v(nz) * i;
  e(i) = sum(abs(iv-round(iv)));
end

j = find(e<10e-2,1,'first');

if ~isempty(j)
  iv = round(v / abs(v(nz)) * j);
else
  iv = [];
end

function s=barchar(i,varargin)

if i<0 && check_option(varargin,'latex')
  s = ['\bar{',int2str(-i),'}'];
else
  s = int2str(i);
end

author: mtex-toolbox

repo: mtex

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


text:
function K = knPoly(X, Y, o, c)
% Polynomial kernel k(x,y)=(x'y+c)^o
% Input:
%   X: d x nx data matrix
%   Y: d x ny data matrix
%   o: order of polynomial
%   c: constant
% Ouput:
%   K: nx x ny kernel matrix
% Written by Mo Chen (sth4nth@gmail.com).
if nargin < 4
    c = 0;
end

if nargin < 3
    o = 3;
end

if nargin < 2 || isempty(Y)
    K = (dot(X,X,1)+c).^o;            % norm in kernel space
else
    K = (X'*Y+c).^o;
end

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function  [s_idx, seg]    =  Proc_cls_idx( cls_idx )

[idx  s_idx]    =  sort(cls_idx);

idx2   =  idx(1:end-1) - idx(2:end);
seq    =  find(idx2);

seg    =  [0; seq; length(cls_idx)];

author: lbasek

repo: image-denoising-benchmark


text:
function redoxsteady
% Steady state redox zones
%    using MATLAB ode
%
%   $Ekkehard Holzbecher  $Date: 2006/03/31 $
%--------------------------------------------------------------------------
L = 100;                        % length [m]
v = 1;                          % velocity [m/s]
D = 0.2;                        % diffusivity [m*m/s]
lambda = 0.01;                  % organic carbon degradation parameter [1/m]
k1 = [0.1; 1; 0.9];             % 1. Michaelis-Menten parameter
k2 = [0.035; 1; 1];             % 2. Michaelis-Menten parameter [kg/m*m*m]
k3 = [3.5e-3; 1];               % inhibition coefficient [kg/m*m*m]
corg = 1;                       % organic carbon concentration at interface
cin = [4; 3; 0.001];            % interface concentrations [kg/m*m*m]
N = 100;                        % number of nodes

%----------------------execution-------------------------------------------

x = linspace(0,L,N);
solinit = bvpinit (x,[cin; zeros(3,1)]);
sol = bvp4c(@redox,@bcs,solinit,odeset,D,v,lambda,k1,k2,k3,corg,cin);

%---------------------- graphical output ----------------------------------

plot (x,corg*exp(-lambda*x),sol.x,sol.y(1:3,:));
legend ('C_{org}','O_2','NO_2','Mn'); grid;

%----------------------functions------------------------------
function dydx = redox(x,y,D,v,lambda,k1,k2,k3,corg,cin)

c0 = corg*exp(-lambda*x);
monod = k1.*(y(1:3)>0).*y(1:3)./(k2+y(1:3));
monod(3) = k1(3);
inhib = k3./(k3+y(1:2));
dydx = zeros (6,1);
dydx(1) =  y(4);
dydx(4) =  (v*y(4)+c0*monod(1))/D;
dydx(2) =  y(5);
dydx(5) =  (v*y(5)+c0*monod(2)*inhib(1))/D;
dydx(3) =  y(6);
dydx(6) =  (v*y(6)-c0*monod(3)*inhib(1)*inhib(2))/D;

function res = bcs (ya,yb,D,v,lambda,k1,k2,k3,corg,cin)
res = [ya(1:3)-cin; yb(4:6)-zeros(3,1)];

author: Sable

repo: mcbench-benchmarks


text:
function K = kuratowski_subgraph(A,varargin)
% KURATOWSKI_SUBGRAPH Identify a Kuratowski subgraph
%
% K = kuratowski_subgraph(A) is just a wrapper around the
% boyer_myrvold_planarity_test to get the subgraph as the first output
% argument.
%
% K is empty if the graph A is planar.
%
% Example:
%   A = clique_graph([3,3]); % Generate K_3,3
%   K = kuratowski_subgraph(A);
%   isequal(A,K) % K_3,3 is a Kuratowski graph!

% David Gleich
% Copyright, Stanford University, 2008

%% History
%  2007-09-29: Initial coding
%%

[is_planar K] = boyer_myrvold_planarity_test(A,varargin{:});

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function showrate3(N1,err1,k1,opt1,str1,N2,err2,k2,opt2,str2,N3,err3,k3,opt3,str3)
%% SHOWRATE3 rate of two error sequences
%
% showrate3(N1,err1,k1,opt1,str1,N2,err2,k2,opt2,str2,N3,err3,k3,opt3,str3)
% plots the err1 vs N1, err2 vs N2, and err3 vs N3 in the loglog scale.
% Additional input
%
%   - k1, k2, k3: specify the starting indices; see showrate
%   - opt1, opt2, opt3: the line color and style
%   - str1, str2, str3: strings used in legend
%
% Example
%
% showrate2(N,energyErr,1,'r-+','||u-u_h||_A',...
%           N,L2Err,1,'b-+','||u-u_h||');
%
% See also showrate, showresult, showmesh, showsolution
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if (nargin<=2)
    k1 = 1; opt1 = '-*';
end
r1 = showrate(N1,err1,k1,opt1);
hold on
r2 = showrate(N2,err2,k2,opt2);
r3 = showrate(N3,err3,k3,opt3);
h_legend = legend(str1,['C_1N^{' num2str(r1) '}'],...
                  str2,['C_2N^{' num2str(r2) '}'],...
                  str3,['C_3N^{' num2str(r3) '}'],'LOCATION','Best');
set(h_legend,'FontSize',12);

repo: Aero-matlab


text:
function value = c8_nint ( c1 )

%*****************************************************************************80
%
%% C8_NINT returns the nearest complex integer of a C8.
%
%  Discussion:
%
%    A C8 is a complex ( kind = 8 ) value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, complex C1, the value to be NINT'ed.
%
%    Output, complex VALUE, the NINT'ed value.
%
  xc = real ( c1 );
  yc = imag ( c1 );
%
%  Lower left.
%
  x = floor ( real ( c1 ) );
  y = floor ( imag ( c1 ) );
  r = ( x - xc ) .^2 + ( y - yc ) .^2;
  r_min = r;
  x_min = x;
  y_min = y;
%
%  Lower right.
%
  x = floor ( real ( c1 ) ) + 1.0;
  y = floor ( imag ( c1 ) );
  r = ( x - xc ) .^2 + ( y - yc ) .^2 ;
  if ( r < r_min )
    r_min = r;
    x_min = x;
    y_min = y;
  end
%
%  Upper right.
%
  x = floor ( real ( c1 ) ) + 1.0;
  y = floor ( imag ( c1 ) ) + 1.0;
  r = ( x - xc ) .^2 + ( y - yc ) .^ 2;
  if ( r < r_min )
    r_min = r;
    x_min = x;
    y_min = y;
  end
%
%  Upper left.
%
  x = floor ( real ( c1 ) );
  y = floor ( imag ( c1 ) ) + 1.0;
  r = ( x - xc ) .^ 2 + ( y - yc ) .^ 2;
  if ( r < r_min )
    r_min = r;
    x_min = x;
    y_min = y;
  end

  value = x_min + i * y_min;

  return
end


text:
function [ n_data, n, c ] = catalan_values ( n_data )

%*****************************************************************************80
%
%% CATALAN_VALUES returns some values of the Catalan numbers.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Binomial[2*n,n] / ( n + 1 )
%
%  First values:
%
%     C(0)     1
%     C(1)     1
%     C(2)     2
%     C(3)     5
%     C(4)    14
%     C(5)    42
%     C(6)   132
%     C(7)   429
%     C(8)  1430
%     C(9)  4862
%    C(10) 16796
%
%  Formula:
%
%    C(N) = (2*N)! / ( (N+1) * (N!) * (N!) )
%         = 1 / (N+1) * COMB ( 2N, N )
%         = 1 / (2N+1) * COMB ( 2N+1, N+1).
%
%  Recursion:
%
%    C(N) = 2 * (2*N-1) * C(N-1) / (N+1)
%    C(N) = sum ( 1 <= I <= N-1 ) C(I) * C(N-I)
%
%  Discussion:
%
%    The Catalan number C(N) counts:
%
%    1) the number of binary trees on N vertices;
%    2) the number of ordered trees on N+1 vertices;
%    3) the number of full binary trees on 2N+1 vertices;
%    4) the number of well formed sequences of 2N parentheses;
%    5) the number of ways 2N ballots can be counted, in order,
%       with N positive and N negative, so that the running sum
%       is never negative;
%    6) the number of standard tableaus in a 2 by N rectangular Ferrers diagram;
%    7) the number of monotone functions from [1..N} to [1..N} which
%       satisfy f(i) <= i for all i;
%    8) the number of ways to triangulate a polygon with N+2 vertices.
%
%  Example:
%
%    N = 3
%
%    ()()()
%    ()(())
%    (()())
%    (())()
%    ((()))
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the order of the Catalan number.
%
%    Output, integer C, the value of the Catalan number.
%
  n_max = 11;

  c_vec = [ ...
    1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796 ];

  n_vec = [ ...
     0,  1,  2,  3,  4, 5,  6,  7,  8,  9,  10 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    c = 0;
  else
    n = n_vec(n_data);
    c = c_vec(n_data);
  end

  return
end


text:
matlabbatch{1}.spm.tools.physio.save_dir = {'physio_out'};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Siemens';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'siemens_PAV.ecg'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {''};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {''};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = 0.0025;
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 0;
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 20;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 2.41;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 5;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 400;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 11;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sync.nominal = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'ECG';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.filter.no = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.max_heart_rate_bpm = 90;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.respiratory.filter.passband = [0.01 2];
matlabbatch{1}.spm.tools.physio.preproc.respiratory.despike = false;
matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt';
matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat';
matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none';
matlabbatch{1}.spm.tools.physio.model.censor_unreliable_recording_intervals = false;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1;
matlabbatch{1}.spm.tools.physio.model.rvt.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.hrv.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.movement.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.other.no = struct([]);
matlabbatch{1}.spm.tools.physio.verbose.level = 2;
matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = 'physio.fig';
matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;

author: translationalneuromodeling

repo: tapas


text:
function [T, M, S, DISTR, df] = dtiValTestStat(g1, g2, Y, mask)

% Computes voxel-wise statistics for two groups from a data array of
% diffusion tensors in dt6 format.
%
% The test is H0: both groups have the same eigenvalues, with possibly
% different unknown eigenvectors.
%
%   [T, M, S, DISTR, df] = dtiValTestStat(g1, g2, DT_ARRAY, [MASK])
%
% Input:
%   g1, g2      List of indices that correspond to each group out of 1:N
%                   E.g: g1 = 1:7, g2 = 8:14, N = 14
%   DT_ARRAY    Data array of size XxYxZx6xN (or nx6xN), where X, Y, Z are the volume
%                   dimensions and N is the number of subjects.
%                   (n is the number of voxels).
%   MASK        Optional XxYxZ binary array. Values of T, M and S are computed
%                   where mask = 1; in other voxels, T, M and S are set to 0.
%                   Default is entire volume.
%
% Output:
%   T           XxYxZx1 (or nx1) array of test statistics (0 where mask = 0)
%   M           XxYxZx6x2 (or nx6x2) array of mean tensors (0 where mask = 0)
%   S           XxYxZx1 (or nx1) array of variances (0 where mask = 0)
%   DISTR       'chi2' or 'f'
%   df          degrees of freedom of the appropriate distribution
%
% Utilities:    ndfun.m, dtiSplitTensor.m
%
% WARNING: If using Pentium 4, eliminate NaN's from array before running
% (processor bug).
%
% Copyright by Armin Schwartzman, 2005

% HISTORY:
%   2004.08.29 ASH (armins@stanford.edu) wrote it.
%

% Check inputs
if (ndims(Y)==2 | ndims(Y)==3),
    Ind = 1;    % Data in indexed nx6xN format
    Y = shiftdim(Y, -2);
else
    Ind = 0;    % Data in XxYxZx6xN format
end
if (ndims(Y)<4 | ndims(Y)>5),
    error('Wrong input format');
end
if (~exist('mask')),
    mask = ones([size(Y,1) size(Y,2) size(Y,3)]);
end

% Computations
N1 = length(g1);
N2 = length(g2);
N  = N1 + N2;

q = size(Y, 4);
p = max(roots([1/2 1/2 -q]));
Yavg1 = mean(Y(:,:,:,:,g1),5);
Yavg2 = mean(Y(:,:,:,:,g2),5);
Yavg = mean(Y(:,:,:,:,:),5);
M = cat(5, Yavg1, Yavg2);

[V1,L1] = dtiSplitTensor(Yavg1);
[V2,L2] = dtiSplitTensor(Yavg2);
[V,L] = dtiSplitTensor(Yavg);

% Total variance
d1 = Y(:,:,:,:,g1) - repmat(Yavg1,[1 1 1 1 N1]);
d2 = Y(:,:,:,:,g2) - repmat(Yavg2,[1 1 1 1 N2]);
S = sum(d1(:,:,:,1:p,:).^2, 4) + 2*sum(d1(:,:,:,p+1:q,:).^2, 4) + ...
    sum(d2(:,:,:,1:p,:).^2, 4) + 2*sum(d2(:,:,:,p+1:q,:).^2, 4);
S = sum(S, 5)/(q*(N-2));

% F version
DISTR = 'f';
df = [p, q*(N-2)];
T = N1*N2/N^2 * sum((L1 - L2).^2, 4);
T = df(2)/df(1) * T./(S*(N-2)/N) / q;

% Adjust output
if Ind,
    T = shiftdim(T, 2);
    M = shiftdim(M, 2);
    S = shiftdim(S, 2);
    L = shiftdim(L, 2);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [xpos, u] = firstmax(s)

%tstoolbox/@signal/firstmax
%   Syntax:
%     * [xpos, unit] = firstmax(s)
%
%   Give information about first local maximum of scalar signal s.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt
narginchk(1,1);

if ndim(s) ~= 1
	help(mfilename)
	return
end

d = data(s);
N = dlens(s,1);

c1 = d(1:N-2);
c2 = d(2:N-1);
c3 = d(3:N);

m = find((c1<c2) .* (c3<c2));
if isempty(m)
	m = find((c1<=c2) .* (c3<c2));
end
if isempty(m)
	m = find((c1<c2) .* (c3<=c2));
end

if isempty(m)
	error('No local minimum found')
else
	A = getaxis(s, 1);
	xpos = first(A) + m(1)*delta(A);
	u = unit(A);
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% CONDSTAT - accumulate surrogate data for comparing two data conditions
%
% Usage:
%     >> [diffres, accres, res1, res2] = condstat(formula, naccu, alpha, ...
%                                              bootside, condboot, arg1, arg2 ...);
%
% Inputs:
%    formula  - [string or cell array of strings] formula(s) to compute a given measure.
%               Takes arguments 'arg1', 'arg2' ... and X as inputs. e.g.,
%        'sum(arg1(:,:,X),3) ./ sqrt(sum(arg2(:,:,X))) ./ sqrt(sum(arg3(:,:,X)))'
%    naccu    - [integer] number of accumulations of surrogate data. e.g., 200
%    alpha    - [float] significance level (0<alpha<0.5)
%    bootside - ['both'|'upper'] side of the surrogate distribution to
%               consider for significance. This parameter affect the size
%               of the last dimension of accumulation array 'accres' (size
%               is 2 for 'both' and 1 for 'upper').
%    condboot - ['abs'|'angle'|'complex'|''] When comparing two conditions, compare
%               either absolute vales ('abs'), angles ('angles') or complex values
%               ('complex'). Either '' or 'complex' leave the formula unchanged;
%               'abs' takes its norm before subtraction, and 'angle' normalizes
%               each value (to norm 1) before taking the difference.
%    arg1     - [cell_array] of two 1D,2D,or 3D matrices of values to compare.
%               The last dimension of the array is shuffled to accumulate data, the
%               other dimensions must be the same size across matrices.
%               e.g. size(arg1{1})=[100 200 500], size(arg1{2})=[100 200 395]
%    arg2     - same as arg1, note that it is compared only to itself, and has
%               nothing to do with arg1 besides using the same formula, alpha, etc.
% ...argn     - may call n number of argument pairs
%
% Outputs:
%    diffres  - difference array for the actual (non-shuffled) data, if more than one
%               arg pair is called, format is a cell array of matrices.
%    accres  -  [cell array of 3D numerical arrays] for shuffled data, one per formula.
%    res1     - result for first condition
%    res2     - result for second condition
%
% Authors: Arnaud Delorme & Scott Makeig
%          CNL/Salk Institute 1998-2001; SCCN/INC/UCSD, La Jolla, 2002-
%
% See also: TIMEF, CROSSF

% Copyright (C) 2002  Arnaud Delorme, Lars Kai Hansen & Scott Makeig, SCCN/INC/UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [diffres, accres, res1, res2] = condstat(formula, naccu, alpha, bootside, condboot, varargin);

if nargin < 6
	help condstat;
	return;
end

if ~ischar(formula) && ~iscell(formula)
	error('The first argument must be a string formula or cell array of string');
end
if ischar(formula)
	formula = { formula };
end
if ischar(bootside)
	bootside = { bootside };
end
for index = 1:length(bootside)
	if ~strcmpi(bootside{index}, 'both') && ~strcmpi(bootside{index}, 'upper')
		error('Bootside must be either ''both'' or ''upper''');
	end
end
if ischar(condboot)
	condboot = { condboot };
end
for index = 1:length(condboot)
	if isempty(condboot{index}), condboot{index} = 'complex'; end
end

for index = 1:length(varargin)
	if ~iscell(varargin) || length(varargin{index}) ~=2
		error('Except for the first arguments, all other arguments given to the function must be cell arrays of two numerical array');
	end
end

% accumulate coherence images (all arrays [nb_points x timesout x trials])
% ---------------------------
for index=1:length(varargin)
	tmpvar1 = varargin{index}{1};
	tmpvar2 = varargin{index}{2};
	if index == 1   % Shouldn't this be recalculated for arg2, etc.? TF 2007.06.04
		cond1trials = size(tmpvar1,ndims(tmpvar1));
		cond2trials = size(tmpvar2,ndims(tmpvar2));
		for tmpi = 1:length(formula)
			accres{tmpi} = zeros(size(tmpvar1,1), size(tmpvar1,2), naccu);
		end
	end

	if ndims(tmpvar1) == 2
		eval( [ 'arg' int2str(index) '=zeros(size(tmpvar1,1), cond1trials+cond2trials);' ] );
		eval( [ 'arg' int2str(index) '(:,1:cond1trials)=tmpvar1;' ] );
		eval( [ 'arg' int2str(index) '(:,cond1trials+1:end)=tmpvar2;' ] );
	else
		eval( [ 'arg' int2str(index) '=zeros(size(tmpvar1,1), size(tmpvar1,2), cond1trials+cond2trials);' ] );
		eval( [ 'arg' int2str(index) '(:,:,1:cond1trials)=tmpvar1;' ] );
		eval( [ 'arg' int2str(index) '(:,:,cond1trials+1:end)=tmpvar2;' ] );
	end
end

fprintf('Accumulating permutation statistics:');
alltrials = [1:cond1trials+cond2trials];

% processing formulas
% -------------------
formula1 = [];
formula2 = [];
for index = 1:length(formula)
	% updating formula
	% ----------------
	switch lower(condboot{ min(length(condboot), index) })
	 case 'abs', formula{index} = [ 'abs(' formula{index} ')' ];
	 %case 'angle', formula{index} = [ 'exp(j*angle(' formula{index} '))' ];
	 case 'angle', formula{index} = [ 'angle(' formula{index} ')/(2*pi)' ];
	 case 'complex', ;
	 otherwise, condboot, error('condstat argument must be either ''abs'', ''angle'', ''complex'' or empty');
	end

	% computing difference (non-shuffled)
	% -----------------------------------
	X = 1:cond1trials;
	eval( [ 'res1{index} = ' formula{index} ';'] );
	X = cond1trials+1:cond1trials+cond2trials;
	eval( [ 'res2{index} = ' formula{index}  ';'] );
	diffres{index} = res1{index}-res2{index};

	% build formula to execute
	% ------------------------
	arrayname = [ 'accres{' int2str(index) '}' ];
	if ndims(tmpvar1) == 2 % 2 dimensions
		formula1 = [ formula1 arrayname '(:,index) = ' formula{index} ';'];
		formula2 = [ formula2 arrayname '(:,index) = ' arrayname '(:,index) - ' formula{index}  ';'];
	else % 3 dimensions
		formula1 = [ formula1 arrayname '(:,:,index) = ' formula{index}  ';'];
		formula2 = [ formula2 arrayname '(:,:,index) = ' arrayname '(:,:,index) - ' formula{index}  ';'];
	end
end

% accumulating (shuffling)
% -----------------------
for index=1:naccu
	if rem(index,10) == 0,  fprintf(' %d',index); end
	if rem(index,120) == 0, fprintf('\n'); end

	alltrials = shuffle(alltrials);
	X = alltrials(1:cond1trials);
	eval( formula1 );
	X = alltrials(cond1trials+1:end);
	eval( formula2 );
end

% significance level
% ------------------
for index= 1:length(formula)
	accarray = accres{index};

    % size = nb_points*naccu
    % size = nb_points*naccu*times
	if ~isreal(accarray)    % might want to introduce a warning here: a complex
                            % result may not be desirable, and a single complex value
                            % in accarray could turn this from a 2-tail to 1-tail
                            % bootstrap test, leading to false positives. Hard to
                            % think of a meaningful warning though...
                            % TF 2007.06.04
		accarray = sqrt(accarray .* conj(accarray)); % faster than abs()
    end

    % compute bootstrap significance level
    i = round(naccu*alpha);
    switch ndims(accarray)
	 case 3
	     accarray = sort(accarray,3); % always sort on naccu (when 3D, naccu is the second dim)
         if strcmpi(bootside{min(length(bootside), index)}, 'upper')
			 accarray = accarray(:,:,floor(naccu-i/2+1));
	     else
			 accarray = accarray(:,:,[end:-1:1]);
			 accarraytmp(:,:,2) = accarray(:,:,ceil(i/2));
			 accarraytmp(:,:,1) = accarray(:,:,floor(naccu-i/2+1));
			 accarray = accarraytmp;
		 end

	 case 2
	     accarray = sort(accarray,2); % always sort on naccu (when 3D, naccu is the second dim)
         if strcmpi(bootside{min(length(bootside), index)}, 'upper')
			 accarray = accarray(:,floor(naccu-i/2+1));
	     else
			 accarraytmp(:,2) = accarray(:,ceil(i/2));
			 accarraytmp(:,1) = accarray(:,floor(naccu-i/2+1));
			 accarray = accarraytmp;
		 end
	 case 1
	     accarray = sort(accarray,1); % always sort on naccu (when 3D, naccu is the second dim)
         if strcmpi(bootside{min(length(bootside), index)}, 'upper')
			 accarray = accarray(floor(naccu-i/2+1));
	     else
			 accarraytmp(2) = accarray(ceil(i/2));
			 accarraytmp(1) = accarray(floor(naccu-i/2+1));
			 accarray = accarraytmp;
		 end
    end
    accres{index} = accarray;
end

if length(res1) == 1
	res1 = res1{1};
	res2 = res2{1};
	diffres = diffres{1};
	accres = accres{1};
end
fprintf('\n');
return;

% writing a function
% ------------------
% $$$ fid = fopen('tmpfunc.m', 'w');
% $$$ fprintf(fid, 'function [accres] = tmpfunc(alltrials, cond1trials, naccu,');
% $$$ for index=1:length(varargin)
% $$$ 	fprintf(fid, 'arg%d', index);
% $$$ 	if index ~=length(varargin), fprintf(fid,','); end
% $$$ end
% $$$ fprintf(fid, ')\n');
% $$$ commandstr = [ 'for index=1:naccu, ' ]
% $$$ % 			   'if rem(index,10) == 0,  disp(index); end;' ];
% $$$ % 			   'if rem(index,10) == 0,  fprintf('' %d'',index); end;' ...
% $$$ %	           'if rem(index,120) == 0, fprintf(''\n''); end;' ];
% $$$ commandstr = [ 	commandstr 'shuffle(alltrials);' ];
% $$$ commandstr = [ 	commandstr 'X = alltrials(1:cond1trials);' ];
% $$$ if ndims(tmpvar1) == 2 % 2 dimensions
% $$$ 	commandstr = [ 	commandstr 'accres(:,index) = ' formula ';'];
% $$$ 	commandstr = [ 	commandstr 'X = alltrials(cond1trials+1:end);'];
% $$$ 	commandstr = [ 	commandstr 'accres(:,index) = accres(:,index)-' formula ';end;'];
% $$$ else
% $$$ 	commandstr = [ 	commandstr 'res1 = ' formula ';' 10];
% $$$ 	commandstr = [ 	commandstr 'X = alltrials(cond1trials+1:end);' 10];
% $$$ 	commandstr = [ 	commandstr 'res2 = ' formula ';' 10];
% $$$ 	commandstr = [ 	commandstr 'accres(:,:,index) = res1 - res2; end;'];
% $$$ end;
% $$$ fprintf(fid, commandstr);
% $$$ fclose(fid);
% $$$ profile on;
% $$$ accres = tmpfunc(alltrials, cond1trials, naccu, arg1);
% $$$ profile report;
% $$$ profile off;
% $$$ return;

% evaluating a command
% --------------------
% $$$ commandstr = [ 'for index=1:naccu, ' ...
% $$$ 			   'if rem(index,10) == 0,  fprintf('' %d'',index); end;' ...
% $$$ 	           'if rem(index,120) == 0, fprintf(''\n''); end;' ];
% $$$ commandstr = [ 	commandstr 'shuffle(alltrials);' ];
% $$$ commandstr = [ 	commandstr 'X = alltrials(1:cond1trials);' ];
% $$$ if ndims(tmpvar1) == 2 % 2 dimensions
% $$$ 	commandstr = [ 	commandstr 'accres(:,index) = ' formula ';'];
% $$$ 	commandstr = [ 	commandstr 'X = alltrials(cond1trials+1:end);'];
% $$$ 	commandstr = [ 	commandstr 'accres(:,index) = accres(:,index)-' formula ';end;'];
% $$$ else
% $$$ 	commandstr = [ 	commandstr 'accres(:,:,index) = ' formula ';'];
% $$$ 	commandstr = [ 	commandstr 'X = alltrials(cond1trials+1:end);'];
% $$$ 	commandstr = [ 	commandstr 'accres(:,:,index) = accres(:,:,index)-' formula ';end;'];
% $$$ end;
% $$$ eval(commandstr);
% $$$ return;
% $$$

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function m = rvecrep(v,c)

% RVECREP  Row vector replicate
%
%   M = rvecrep(V, C) Replicates a 1xN dimensional row vector V, C times to generate a
%   CxN dimensional matrix M.
%
%   See also
%   CVECREP, REPMAT

%   Copyright (c) Oregon Health & Science University (2006)
%
%   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for
%   academic use only (see included license file) and can be obtained from
%   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the
%   software should contact rebel@csee.ogi.edu for commercial licensing information.
%
%   See LICENSE (which should be part of the main toolkit distribution) for more
%   detail.

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

if isempty(v)

    m = zeros(c,0);

else

    m = v(ones(1,c),:);

end

repo: ekfmonoslam


text:
function inputString = writeLineFromVar(defaultClass, variableName, maskVars, maskViz, classNum, structName)
% This function creates a string that is written to a wecSimInputFile.m
% line based on a class abbreviation, variable name and its mask value
%
% It is possible to write the input file with a long list of fprintf lines
% such as:
%     fprintf(fid,'simu.rampTime = %s; \r\n',maskVars.rampTime);
%     fprintf(fid,'simu.endTime = %s; \r\n',maskVars.endTime);

% However this doesn't allow the maskVariable to be checked against a
% default value or its visibility (its necessity) without an excessive
% number of if statements in writeInputFromBlocks.m. Comparing against
% these parameters before writing reduces the number of extraneous values written to the input file.
%
% Parameters
% ------------
%     defaultClass : WEC-Sim Object
%         Default instance of a WEC-Sim class. Must contain the variableName that is
%         being written to the input file
%
%     variableName : string
%         Name of the variable within defaultClass that is being written to
%         the input file.
%
%     maskVars : struct
%         Structure containing all mask parameter values
%
%     maskViz : string
%          Structure containing all mask parameter visibilities
%
%     classNum : int
%          Number specifying the index of a class being written
%          i.e. 1 for body(1), 4 for pto(4), 7 for constraint(7), ...
%
%     structName : string
%          Structure containing all mask parameter visibilities
%
% Returns
% -------
%     inputString : string
%         string to write to the input file
%
% TODO
% ----
%     Link mask tool tips with input line comments
%

    % Track if variable should be written
    hasStruct = ~isempty(structName); % check if variable is in a class' struct (correspond to a mask tab)
    isVisible = strcmp(maskViz.(variableName),'on'); % check if mask variable is visible

    % Get the default value of the variable within its WEC-Sim class
    if hasStruct
        isDefault = isequal(defaultClass.(structName).(variableName), eval(maskVars.(variableName)));
    else
        isDefault = isequal(defaultClass.(variableName), eval(maskVars.(variableName)));
    end

    % Only write parameters if they are visible (turned on and relevant) and
    % are different from the class default
    if isVisible && ~isDefault
        % Append the class index if necessary. E.g. 'body' --> 'body(1)'
        classAbbrev = inputname(1);
        if ~isempty(classNum)
            classAbbrev = [classAbbrev '(' num2str(classNum) ')'];
        end

        if hasStruct
            % e.g. 'body(1).initial.displacement = [1 1 1]; \r\n'
            inputString = [classAbbrev '.' structName '.' variableName ' = ' maskVars.(variableName) '; \r\n'];
        else
            % e.g. 'simu.stateSpace = 'on'; \r\n'
            inputString = [classAbbrev '.' variableName ' = ' maskVars.(variableName) '; \r\n'];
        end
    else
        % Write nothing to input file
        inputString = '';
    end

end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function kpn_test ( )

%*****************************************************************************80
%
%% KPN_TEST uses the KPN function for 1D quadrature over (-oo,+oo).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 May 2012
%
%  Author:
%
%    John Burkardt
%
  d = 1;
  exact = hermite_integral ( d );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'KPN_TEST:\n' );
  fprintf ( 1, '  Kronrod-Patterson-Hermite quadrature over (-oo,+oo):\n' );
  fprintf ( 1, '  Exact integral is %g\n', exact );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   Level   Nodes    Estimate  Error\n' );
  fprintf ( 1, '\n' );

  for l = 1 : 5

    n = kpn_order ( l );

    nh = floor ( ( n + 1 ) / 2 );

    [ xh, wh ] = kpn ( l );

    if ( mod ( n, 2 ) == 1 )
      x = [ -xh(nh:-1:2); xh ];
      w = [ wh(nh:-1:2); wh ];
    else
      x = [ -xh(nh:-1:1); xh ];
      w = [ wh(nh:-1:1); wh ];
    end

    fx = hermite_integrand ( d, n, x );
    q = w' * fx;
    e = sqrt ( ( q - exact ).^2 ) / exact;

    fprintf( '  %2d     %6d  %10.5g  %10.5g\n', l, n, q, e )

  end

  return
end


text:
function [M q]=modularity_probtune_und_sign(W,qtype,M,p)
%MODULARITY_PROBTUNE_UND_SIGN     Optimal community structure, modularity and degeneracy
%
%   Ci     = modularity_probtune_und_sign(W,   [],Ci0,0.01);
%   Ci     = modularity_probtune_und_sign(W,'sta',Ci0,0.01);
%   [Ci Q] = modularity_probtune_und_sign(W,'sta',Ci0,0.01);
%
%   The optimal community structure is a subdivision of the network into
%   nonoverlapping groups of nodes in a way that maximizes the number of
%   within-group edges, and minimizes the number of between-group edges.
%   The modularity is a statistic that quantifies the degree to which the
%   network may be subdivided into such clearly delineated groups.
%   High-modularity degeneracy is the presence of many topologically
%   distinct high-modularity partitions of the network.
%
%   This algorithm is inspired by the Kernighan-Lin fine-tuning algorithm
%   and is designed to probabilistically refine a previously detected
%   community by incorporating random node moves into a finetuning
%   algorithm.
%
%   Input:      W,      undirected (weighted or binary) connection matrix
%                       with positive and negative weights
%
%               qtype,  modularity type (see Rubinov and Sporns, 2011)
%                           'sta',  Q_* (default if qtype is not specified)
%                           'pos',  Q_+
%                           'smp',  Q_simple
%                           'gja',  Q_GJA
%                           'neg',  Q_-
%
%               Ci0,    initial community affiliation vector (optional)
%
%               p,      probability of random node moves
%
%
%   Output:     Ci,     refined community affiliation vector
%               Q,      modularity (qtype dependent)
%
%   Note: Ci and Q may vary from run to run, due to heuristics in the
%   algorithm. Consequently, it may be worth to compare multiple runs.
%
%   References:
%   Rubinov and Sporns (2011) NeuroImage.
%   Sun et al. (2008)  Europhysics Lett 86, 28004.
%
%
%   Mika Rubinov, UNSW, 2011

%   Modification History:
%   Mar 2011: Original

n=length(W);                                                %number of nodes/modules
if ~exist('qtype','var') || isempty(qtype);
    qtype = 'sta';
end
if isempty(M);
    M = 1:n;
else
    [dum dum M] = unique(M(:).');                           %align module indices
end

W0= W.*(W>0);                                               %positive weights matrix
W1=-W.*(W<0);                                               %negative weights matrix
s0=sum(W0(:));                                              %positive sum of weights
s1=sum(W1(:));                                              %negative sum of weights
Knm0=zeros(n,n);                                            %positive node-to-module degree
Knm1=zeros(n,n);                                            %negative node-to-module degree
for m=1:max(M)                                              %loop over modules
    Knm0(:,m)=sum(W0(:,M==m),2);
    Knm1(:,m)=sum(W1(:,M==m),2);
end
Kn0=sum(Knm0,2);                                            %positive node degree
Kn1=sum(Knm1,2);                                            %negative node degree
Km0=sum(Knm0,1);                                            %positive module degree
Km1=sum(Knm1,1);                                            %negative module degree

switch qtype
    case 'smp';  d0 = 1/s0;       d1 = 1/s1;                %dQ = dQ0/s0 - dQ1/s1;
    case 'gja';  d0 = 1/(s0+s1);  d1 = 1/(s0+s1);           %dQ = (dQ0 - dQ1)/(s0+s1);
    case 'sta';  d0 = 1/s0;       d1 = 1/(s0+s1);           %dQ = dQ0/s0 - dQ1/(s0+s1);
    case 'pos';  d0 = 1/s0;       d1 = 0;                   %dQ = dQ0/s0;
    case 'neg';  d0 = 0;          d1 = 1/s1;                %dQ = -dQ1/s1;
    otherwise; error('qtype unknown');
end
if ~s0                                                      %adjust for absent positive weights
    s0=1;
    d0=0;
end
if ~s1                                                      %adjust for absent negative weights
    s1=1;
    d1=0;
end

for u=randperm(n);                                          %loop over all nodes in random order
    ma = M(u);                                              %current module of u
    r=rand<p;
    if r                                                    %with probability p
        mb = ceil(rand*n);                                  %choose random new module
    else
        dQ0 = (Knm0(u,:)+W0(u,u)-Knm0(u,ma)) - Kn0(u).*(Km0+Kn0(u)-Km0(ma))/s0;     %positive dQ
        dQ1 = (Knm1(u,:)+W1(u,u)-Knm1(u,ma)) - Kn1(u).*(Km1+Kn1(u)-Km1(ma))/s1;     %negative dQ
        dQ = d0*dQ0 - d1*dQ1;                               %rescaled changes in modularity
        dQ(ma) = 0;                                         %no changes for same module
        [max_dQ mb] = max(dQ);                              %maximal increase in modularity and corresponding module
    end
    if r || (max_dQ>1e-10)                                  %if maximal increase is positive (equiv. dQ(mb)>dQ(ma))
        M(u) = mb;                                          %reassign module

        Knm0(:,mb)=Knm0(:,mb)+W0(:,u);
        Knm1(:,mb)=Knm1(:,mb)+W1(:,u);
        Knm0(:,ma)=Knm0(:,ma)-W0(:,u);
        Knm1(:,ma)=Knm1(:,ma)-W1(:,u);
        Km0(mb)=Km0(mb)+Kn0(u);
        Km1(mb)=Km1(mb)+Kn1(u);
        Km0(ma)=Km0(ma)-Kn0(u);
        Km1(ma)=Km1(ma)-Kn1(u);
    end
end

[dum dum M]=unique(M(:).');                                 %realign module indices
if nargout==2                                               %compute modularity
    m = M(ones(1,n),:);
    Q0 = (W0-(Kn0*Kn0.')/s0).*(m==m.');
    Q1 = (W1-(Kn1*Kn1.')/s1).*(m==m.');
    q = d0*sum(Q0(:)) - d1*sum(Q1(:));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function rate_na_lt = awgn_na_lt(snr_db,nn,error)

P=10^(snr_db/10);

rate_na_lt=zeros(1,length(nn));

%find the minimum blocklength required for the long-term power constraint to be beneficial
C=log(1+P);
V = 1-1/(1+P)^2;
n_min = ( (1+P)/P*sqrt(2*3.1415926*V)*(1-error)*exp(qfuncinv(error)^2/2) + qfuncinv(error)/(1+P)^2/sqrt(V))^2;

n_index=find(nn <= n_min);

% for every n<=n_min, the normal approx. for the st power constraint coincides
if length(n_index)>0
    rat_na_lt(n_index) = (log(1+P) - sqrt(1-(1+P)^(-2))./sqrt(nn(n_index)) + 0.5*log(nn(n_index))./nn(n_index))/log(2);
end

R_st = C - sqrt(V./nn)*qfuncinv(error);
C_lt=log(1+P/(1-error));

index_large_n = find(nn>n_min);

for ii=index_large_n
    n=nn(ii);
    for R_star = R_st(ii):0.0001:C_lt
        omega = P: 0.00001:P*1.2; %note that larger endpoing for omega (P*1.2) gives a better estimate for rate_na;

        q_arg=sqrt(n)*(log(1+omega)-R_star)./sqrt(1-(1+omega).^(-2)) ;

        error_temp = min( 1-P./omega + P.*qfunc( q_arg)./omega);
        if error_temp> error
            rate_na_lt(ii)= (R_star +0.5*log(n)./n)./log(2) ;
            break
        end
    end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function lattice_print ( dim_num, m, z, title )

%*****************************************************************************80
%
%% LATTICE_PRINT prints the points in a lattice rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994,
%    ISBN: 0198534728,
%    LC: QA311.S56
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer M, the number of points to use.
%
%    Input, integer Z(DIM_NUM), the generator vector.
%
%    Input, string TITLE, an optional title.
%
  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  fprintf ( 1, '\n' );

  for i = 0 : m - 1
    y(1:dim_num) = mod ( i * z(1:dim_num), m );
    fprintf ( 1, '%4d    ', i + 1 );
    for dim = 1 : dim_num
      fprintf ( 1, '%4d', y(dim) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%% RECOMPUTES THE REPROJECTION ERROR %%%%%%%%%%%%%%%%%%%%%%%%

check_active_images;

% Reproject the patterns on the images, and compute the pixel errors:

ex = []; % Global error vector
x = []; % Detected corners on the image plane
y = []; % Reprojected points

if ~exist('alpha_c'),
   alpha_c = 0;
end;

for kk = 1:n_ima,

   eval(['omckk = omc_' num2str(kk) ';']);
   eval(['Tckk = Tc_' num2str(kk) ';']);

   if active_images(kk) & (~isnan(omckk(1,1))),

      %Rkk = rodrigues(omckk);

      eval(['y_' num2str(kk) '  = project_points2(X_' num2str(kk) ',omckk,Tckk,fc,cc,kc,alpha_c);']);

      eval(['ex_' num2str(kk) ' = x_' num2str(kk) ' - y_' num2str(kk) ';']);

      eval(['x_kk = x_' num2str(kk) ';']);

      eval(['ex = [ex ex_' num2str(kk) '];']);
      eval(['x = [x x_' num2str(kk) '];']);
      eval(['y = [y y_' num2str(kk) '];']);

   else

      %	eval(['y_' num2str(kk) '  = NaN*ones(2,1);']);

      % If inactivated image, the error does not make sense:
      eval(['ex_' num2str(kk) ' = NaN*ones(2,1);']);

   end;

end;

err_std = std(ex')';


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: September 9th, 2016
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Make_Your_Plates_and_Input_Files()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  64;       % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
dx = Lx/Nx;      % Grid spatial resolution

%
% Immersed Structure Geometric / Dynamic Parameters %
%
ds = 0.5*dx;              % Lagrangian Pt. Spacing (2x resolution of Eulerian grid)
struct_name = 'corals';   % Name for .vertex, .spring, etc files. (must match what's in 'input2d')

% Call function to construct geometry
[xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,Nx,Lx);

% Plot Geometry to test
plot(xLag,yLag,'r-'); hold on;
plot(xLag,yLag,'*'); hold on;
xlabel('x'); ylabel('y');
axis([0 Lx 0 Lx]);

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .spring file!
% k_Spring = 2.5e4;                    % Spring stiffness (does not need to be equal for all springs)
% ds_Plates = dist;                    % Spring resting length (does not need to be equal for all springs)
% print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Plates,ds,struct_name);

% Prints .beam file!
% k_Beam = 0.5;                      % Beam Stiffness (does not need to be equal for all beams)
% C = compute_Curvatures(xLag,yLag)  % Computes curvature of initial configuration
%print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name);

% Prints .target file!
k_Target = 1e7;
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called 'struct_name'.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag); % Total # of Lag. Pts

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called 'struct_name'.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Plates,ds,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N-2 + N/2 );    % Print # of springs

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N/2
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds);
            elseif s > N/2
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds);
            end
    end

    % SPRINGS ACROSS PLATES
    for s=1:N/2
        fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+N/2, k_Spring, ds_Plates);
    end

    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called 'struct_name'.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called 'struct_name'.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C(s) );
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C(s) );
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes "curvature" of starting configuration
%
% NOTE: not curvature in the traditional geometric sense, in the 'discrete'
% sense through cross product.
%
% NOTE: assumes a CLOSED structure
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C = compute_Curvatures(xLag,yLag)

N = length(xLag);
C = zeros( N );

%Note: needs to be done same order as you print .beam file!
for i=1:N

   % Pts Xp -> Xq -> Xr (same as beam force calc.)

   if ( (i > 1) && (i < N) )

        Xp = xLag(i-1); Xq = xLag(i); Xr = xLag(i+1);
        Yp = yLag(i-1); Yq = yLag(i); Yr = yLag(i+1);

   elseif (i==1)

        Xp = xLag(N); Xq = xLag(i); Xr = xLag(i+1);
        Yp = yLag(N); Yq = yLag(i); Yr = yLag(i+1);

   elseif (i==N)

        Xp = xLag(N-1); Xq = xLag(N); Xr = xLag(1);
        Yp = yLag(N-1); Yq = yLag(N); Yr = yLag(1);

   end

   C(i) = (Xr-Xq)*(Yq-Yp) - (Yr-Yq)*(Xq-Xp); %Cross product btwn vectors

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,Nx,L)

% ds: Lagrangian pt. spacing
% Nx: Eulerian grid resolution
% L:  Length of computational domain

ratio = 1024 / (2*Nx);

% Gives LEFT Coral
[x1,y1,N] = please_Give_Single_Polyp_Geometry(ratio);
x1 = x1 + L/10;

% Makes RIGHT Coral
x2 = x1 + L/5;
y2 = y1;

% TESTS GEOMETRY BY PLOTTING
plot(x1,y1,'k*'); hold on;
plot(x1(1),y1(1),'g*'); hold on;
plot(x1(end/2),y1(end/2),'b*'); hold on;
plot(x1(end/2+1),y1(end/2+1),'g*'); hold on;
plot(x1(end),y1(end),'b*'); hold on;
%
plot(x2,y2,'r*'); hold on;
plot(x2(1),y2(1),'g*'); hold on;
plot(x2(end/2),y2(end/2),'b*'); hold on;
plot(x2(end/2+1),y2(end/2+1),'g*'); hold on;
plot(x2(end),y2(end),'b*'); hold on;

xRef = [x1 x2]; yRef = [y1 y2];
ang = pi/4;

% Store Values for Centers of Rotation
xL_1 = xRef(1);       yL_1 = yRef(1);           % Left side of Left Pair
xR_1 = xRef(end/4+1); yR_1 = yRef(end/4+1);     % Right side of Left Pair

xL_2 = xRef(end/2+1);   yL_2 = yRef(end/2+1);   % Left side of Right Pair
xR_2 = xRef(3*end/4+1); yR_2 = yRef(3*end/4+1); % Right side of Right Pair

% -> Rotate Geometry <- %
% LEFT PAIR %
[xR_Ref_1,yR_Ref_1] = rotate_Geometry(ang,xR_1,yR_1,xRef(end/4+1:end/2),yRef(end/4+1:end/2) );
[xL_Ref_1,yL_Ref_1] = rotate_Geometry(-ang,xL_1,yL_1,xRef(1:end/4),yRef(1:end/4) );
% RIGHT PAIR %
[xR_Ref_2,yR_Ref_2] = rotate_Geometry(ang,xR_2,yR_2,xRef(3*end/4+1:end),yRef(3*end/4+1:end) );
[xL_Ref_2,yL_Ref_2] = rotate_Geometry(-ang,xL_2,yL_2,xRef(end/2+1:3*end/4),yRef(end/2+1:3*end/4) );

xLag = [xL_Ref_1 xR_Ref_1 xL_Ref_2 xR_Ref_2];
yLag = [yL_Ref_1 yR_Ref_1 yL_Ref_2 yR_Ref_2];

plot(xLag,yLag,'ro'); hold on;
axis([0 1 0 1]);

% Combine into ONE Vector
%xLag = [xLag_L xLag_R];
%yLag = [yLag yLag];

% Plot the Geometry
% plot(xLag,yLag,'*'); hold on;
% axis([0 L 0 L]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: rotate geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,y] = rotate_Geometry(ang,xC,yC,xRef,yRef)

len = length(xRef);
x = zeros(1,len); y=x;

xRef = xRef - xC;
yRef = yRef - yC;

for i=1:len
   x(i) =  xRef(i)*cos(ang) - yRef(i)*sin(ang);
   y(i) =  xRef(i)*sin(ang) + yRef(i)*cos(ang);
end

x = x + xC;
y = y + yC;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the vertices for one coral polyp
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,y,N] = please_Give_Single_Polyp_Geometry(ratio)

xoffset = 0.3; % To put geometry into QUADRANT-1
yoffset = 0.5; % To put geometry into QUADRANT-1

% Note: (1) ds for stored data is 0.6/(2*1024)
%       (2) 'ratio' is comparing 1024:desired resolution

% Get LEFT side geometry
struct_name1 = 'coral2d_left_1024';
[~,x1,y1] = read_Vertex_Points(struct_name1);

% Reverse order so first point is on bottom
x1 = x1(end:-1:1); y1 = y1(end:-1:1);

% Get RIGHT side geometry
struct_name2 = 'coral2d_right_1024';
[~,x2,y2] = read_Vertex_Points(struct_name2);

% Reverse order so first point is on bottom
x2 = x2(end:-1:1); y2 = y2(end:-1:1);

% Put Geometry Together for One Polyp
xAux = [x1; x2]; yAux = [y1; y2];
xAux = xAux+xoffset; yAux = yAux+yoffset;

% Pull out correct resolution
x=xAux(1:ratio:end)';
y=yAux(1:ratio:end)';
N = length(x)/2;

%plot(xAux,yAux,'*'); hold on;
%plot(x,y,'ro'); hold on;

%NOTE: (1) N here is the # of pts. on one arm (NOT ENTIRE POLYP)!!
%      (2) The geometry here is listed for 1024x1024 meshes -> will need
%          to take every 8 or 16 pts to render geometry usable for MATLAB
%          code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of vertex pts and all the vertex pts from the
%           .vertex file.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [N,xLag,yLag] = read_Vertex_Points(struct_name)

filename = [struct_name '.vertex'];  %Name of file to read in
fileID = fopen(filename);

% Read in the file, use 'CollectOutput' to gather all similar data together
% and 'CommentStyle' to to end and be able to skip lines in file.
C = textscan(fileID,'%f %f','CollectOutput',1);

fclose(fileID);     %Close the data file.

vertices = C{1};    %Stores all read in data in vertices (N+1,2) array

N = vertices(1,1);  % # of Lagrangian Pts
xLag = zeros(N,1);  % Initialize storage for Lagrangian Pts.
yLag = xLag;        % Initialize storage for Lagrangian Pts.

for i=1:N
   xLag(i,1) = vertices(i+1,1); %Stores x-values of Lagrangian Mesh
   yLag(i,1) = vertices(i+1,2); %Stores y-values of Lagrangian Mesh

end

author: nickabattista


text:
%SerialLink.coriolis Coriolis matrix
%
% C = R.coriolis(Q, QD) is the Coriolis/centripetal matrix (NxN) for
% the robot in configuration Q and velocity QD, where N is the number of
% joints.  The product C*QD is the vector of joint force/torque due to velocity
% coupling.  The diagonal elements are due to centripetal effects and the
% off-diagonal elements are due to Coriolis effects.  This matrix is also
% known as the velocity coupling matrix, since it describes the disturbance forces
% on any joint due to velocity of all other joints.
%
% If Q and QD are matrices (KxN), each row is interpretted as a joint state
% vector, and the result (NxNxK) is a 3d-matrix where each plane corresponds
% to a row of Q and QD.
%
% C = R.coriolis( QQD) as above but the matrix QQD (1x2N) is [Q QD].
%
% Notes::
% - Joint viscous friction is also a joint force proportional to velocity but it is
%   eliminated in the computation of this value.
% - Computationally slow, involves N^2/2 invocations of RNE.
%
% See also SerialLink.rne.

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

function C = coriolis(robot, q, qd)

    n = robot.n;

    if nargin == 2
        % coriolis( [q qd] )
        if numcols(q) ~= 2*n
            error('RTB:coriolis:badarg', 'arg must have %d columns', 2*n);
        end
        qd = q(:,n+1:end);
        q = q(:,1:n);
    else
        if numcols(q) ~= n
            error('RTB:coriolis:badarg', 'Cq must have %d columns', n);
        end
        if numcols(qd) ~= n
            error('RTB:coriolis:badarg', 'qd must have %d columns', n);
        end
    end

    % we need to create a clone robot with no friciton, since friction
    % is also proportional to joint velocity
    robot2 = robot.nofriction('all');

    if numrows(q) > 1
        if numrows(q) ~= numrows(qd)
            error('RTB:coriolis:badarg', 'for trajectory q and qd must have same number of rows');
        end
        C = [];
        for i=1:numrows(q)
            C = cat(3, C, robot2.coriolis(q(i,:), qd(i,:)));
        end
        return
    end

    N = robot2.n;

    if isa(q, 'sym')
        C(N,N) = sym();
        Csq(N,N) = sym();
    else

        C = zeros(N,N);
        Csq = zeros(N,N);
    end

    % find the torques that depend on a single finite joint speed,
    % these are due to the squared (centripetal) terms
    %
    %  set QD = [1 0 0 ...] then resulting torque is due to qd_1^2
    for j=1:N
        QD = zeros(1,N);
        QD(j) = 1;
        tau = robot2.rne(q, QD, zeros(size(q)), [0 0 0]');
        Csq(:,j) = Csq(:,j) + tau.';
    end

    % find the torques that depend on a pair of finite joint speeds,
    % these are due to the product (Coridolis) terms
    %  set QD = [1 1 0 ...] then resulting torque is due to
    %    qd_1 qd_2 + qd_1^2 + qd_2^2
    for j=1:N
        for k=j+1:N
            % find a product term  qd_j * qd_k
            QD = zeros(1,N);
            QD(j) = 1;
            QD(k) = 1;
            tau = robot2.rne(q, QD, zeros(size(q)), [0 0 0]');
            C(:,k) = C(:,k) + (tau.' - Csq(:,k) - Csq(:,j)) * qd(j)/2;
            C(:,j) = C(:,j) + (tau.' - Csq(:,k) - Csq(:,j)) * qd(k)/2;

        end
    end

    C = C + Csq * diag(qd);

    if isa(q, 'sym')
        C = simplify(C);
    end

author: Allopart

repo: rbpf-gmapping


text:
function p=dir_pdf(x,a)
%DIR_PDF   Probability density function of uniform Dirichlet
%          distribution
%
%       Description:
%       P = DIR_PDF(X, A) returns the pdf of Dirichlet distribution
%       with A at X
%
% Copyright (c) 2000 Aki Vehtari

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

p=exp(gammaln(sum(a))-sum(gammaln)+sum(log(x).^(a-1)));

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
% gettext() - This function prints a dialog box on screen and waits for
%             the user to enter a string. There is a cancel button which
%             returns a value of [].
% Usage:
%   >> out = gettext(label1,label2,...,label7);
%
% Author: Colin Humphries, CNL / Salk Institute, La Jolla, 1997

% Copyright (C) 1997 Colin Humphries, Salk Institute
%
% 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

% 01-25-02 reformated help & license -ad

function out = gettext(label1,label2,label3,label4,label5,label6,label7);

global is_text_entered
is_text_entered = 0;

for i = 1:nargin
   eval(['leng(i) = length(label',num2str(i),');'])
end
figurelength = 15*1.0*max(leng);

if figurelength < 240
   figurelength = 290;
end

figureheight = 80+nargin*18;
% Set figure

figure('Position',[302 396 figurelength figureheight],'color',[.6 .6 .6],'NumberTitle','off')
f1 = gcf;

ax = axes('Units','Pixels','Position',[0 0 figurelength figureheight],'Visible','off','XLim',[0 figurelength],'YLim',[0 figureheight]);

for i = 1:nargin

   eval(['text(figurelength/2,',num2str(figureheight-(i-1)*18-20),',label',num2str(i),',''color'',''k'',''FontSize'',16,''HorizontalAlignment'',''center'')'])

end
% Set up uicontrols

TIMESTRING = ['global is_text_entered;is_text_entered = 1;clear is_text_entered'];

u = uicontrol('Style','Edit','Units','Pixels','Position',[15 20 figurelength-95 25],'HorizontalAlignment','left','Callback',TIMESTRING);

TIMESTRING = ['global is_text_entered;is_text_entered = 2;clear is_text_entered'];
v = uicontrol('Style','Pushbutton','Units','Pixels','Position',[figurelength-70 20 55 25],'String','Cancel','Callback',TIMESTRING);

while(is_text_entered == 0)
   drawnow
end

if is_text_entered == 1
   out = get(u,'string');
else
   out = [];
end

clear is_text_entered
delete(f1)

author: PatternRecognition

repo: OpenBMI


text:
function i4vec_sorted_unique_test (  )

%*****************************************************************************80
%
%% I4VEC_SORTED_UNIQUE_TEST tests I4VEC_SORTED_UNIQUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 20;
  b = 0;
  c = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_SORTED_UNIQUE_TEST\n' );
  fprintf ( 1, '  I4VEC_SORTED_UNIQUE finds unique entries in a sorted array.\n' );

  seed = 123456789;

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  a = i4vec_sort_heap_a ( n, a );

  i4vec_print ( n, a, '  Input vector:' );

  [ unique_num, a_unique ] = i4vec_sorted_unique ( n, a );

  i4vec_print ( unique_num, a_unique, '  Unique entries:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [fileNameSourceArray, fileNameTargetArray, ...
    fileNameSourceMovedArray, fileNameTargetMovedArray] = tapas_uniqc_move_with_hdr(...
    fileNameSourceArray, fileNameTargetArray)
% Moves given files; for .nii (nifti-) files, also moves .mat-header,
% if existing; for given .img (analyze) files, also moves .hdr-header
%
%   [fileNameSourceArray, fileNameTargetArray] = tapas_uniqc_move_with_hdr(...
%    fileNameSourceArray, fileNameTargetArray)
%
% IN
%   fileNameSourceArray   cell of filenames to be moved (source)
%   fileNameTargetArray   cell of filenames to be moved to (targets)
% OUT
%   fileNameSourceArray   cell of filenames that were tried to be moved
%                   (includes   .mat files corresponding to .nii and
%                               .hdr files corresponding to .img)
%   fileNameSourceMovedArray   cell of source filenames that were actually moved
%                   (includes   .mat files corresponding to .nii and
%                               .hdr files corresponding to .img)
%   fileNameTargetMovedArray   cell of target filenames that were actually moved
%                   (includes   .mat files corresponding to .nii and
%                               .hdr files corresponding to .img)
%
% EXAMPLE
%   tapas_uniqc_move_with_hdr('from.nii', 'to.nii')
%
%   See also tapas_uniqc_delete_with_hdr copy_with_mat

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-08
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

doMove = true;
[fileNameSourceArray, fileNameTargetArray, ...
    fileNameSourceMovedArray, fileNameTargetMovedArray] = tapas_uniqc_copy_with_hdr(...
    fileNameSourceArray, fileNameTargetArray, doMove);


text:
function setcustompath(option_path)

% SETCUSTOMPATH is a helper function that updates the path to contain all the
% custom directories, while keeping the peer or qsub path untouched.

% Copyright (C) 2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are for remembering the path on subsequent calls with the same input arguments
persistent previous_argin previous_path

% the following pertains to path settings, which does not apply in case
% this code is included in a compiled MATLAB application
if ~isdeployed

  if isequal(previous_argin, option_path) && isequal(previous_path, path)
    % no reason to change the path again
    return
  end

  if ~isempty(option_path)
    f = mfilename('fullpath'); % this is .../peer/private/getcustompath.m
    f = fileparts(f);          % this is .../peer/private
    f = fileparts(f);          % this is .../peer
    [p, d] = getcustompath;
    % ensure that the default path is maintained
    path(d);
    % add the custom path
    path(option_path, path);
    % ensure that the peer functions are maintained on the path
    path(f, path);
  end

  % remember the current settings for the next call
  previous_argin = option_path;
  previous_path  = path;

end


text:
function modelPlane = sr_extract_plane(imgPath, imgName, opt)

% SR_EXTRACT_PLANE:
%
% Extract plane model from an image
%
% The code is adapted and modified from the paper
%
% Jia-Bin Huang, Sing Bing Kang, Narendra Ahuja, Johannes Kopf
% Image Completion using Planar Structure Guidance
% ACM Transactions on Graphics (Proceedings of SIGGRAPH 2014).
%
% Output: modelPlane
% The model plane has the following data structure
%
% modelPlane
% modelPlane.vpData:                        Vanishing point data
% modelPlane.numPlane                       Number of planes
% modelPlane.plane{indPlane}.vLine          The vanishing line of the plane
% modelPlane.plane{indPlane}.imgPlaneProb;  The planar location density
% modelPlane.plane{indPlane}.sourceVP       Which two VPs form the plane
% modelPlane.plane{indPlane}.rotPar(vpInd)  Rotation parameters for aligning
%                                           two sets of lines with the x-axis                                           x-axis
% modelPlane.plane{indPlane}.postProb       Posteior probability of the plane

% =========================================================================
% Vanishing point detection
% =========================================================================

vpFilePath = 'cache\vpdetection';
%vpFilePath = [fileparts(fileparts(mfilename('fullpath'))), '\cache\vpdetection'];
vpFileName  = [imgName(1:end-4), '-vanishingpoints.txt'];

recomputeFlag = 1;
if(~exist(fullfile(vpFilePath, 'text', vpFileName), 'file') || recomputeFlag)
    vpExeFile = 'source\EdgeDetectTest.exe';
    vpDetectCMD = [vpExeFile, ' -indir ', imgPath, ' -infile ', imgName, ' -outdir ', vpFilePath];
    system(vpDetectCMD);
end
% Read vanishing point data
vpData = sr_read_vpdata(fullfile(vpFilePath, 'text', vpFileName));

img = imread(fullfile(imgPath, imgName));

% =========================================================================
% Plane localization
% =========================================================================

modelPlane = sr_detect_plane_from_vp(vpData, img, opt);
% figure(1), imshow(modelPlane.postProb(:,:,1:3));

visVPFlag = 0;
if(visVPFlag)
    [vpVis, planeVis] = vis_vp(img, vpData, modelPlane.postProb(:,:,1:3));
end
end

function [vpVis, planeVis] = vis_vp(img, vpData, postProb)

img = im2double(img);

[imgH, imgW, nCh] = size(img);
vpVis = zeros(imgH, imgW, nCh);
planeVis = zeros(imgH, imgW, nCh);
whiteImg = ones(imgH, imgW, nCh);
% blackImg = zeros()
alphaW = 0.5;

imgW = img*(1 - alphaW) + whiteImg*(alphaW);

% Plot VP
h1 = figure(1);
imshow(imgW); hold on;
lineWidth = 4;
for vpInd = 1: vpData.numVP
    vpCurr = vpData.vp{vpInd};
    for lineInd = 1: vpCurr.numLines
        lineCur = vpCurr.lines(lineInd,:);
        lineCur = lineCur + 1;
        if(vpInd == 1)
            plot(lineCur([1,3]), lineCur([2,4]), 'r', 'LineWidth', lineWidth);
        elseif(vpInd == 2)
            plot(lineCur([1,3]), lineCur([2,4]), 'g', 'LineWidth', lineWidth);
        elseif(vpInd == 3)
            plot(lineCur([1,3]), lineCur([2,4]), 'b', 'LineWidth', lineWidth);
        end
    end
end
hold off;
print(h1, '-dpng', fullfile('paper\planar_struct_SR_CVPR2015\figures\plane', 'vpdetection.png'));

% Plot posterior
alphaP = 0.9;
imgP = (1 - alphaP)*img + alphaP*postProb;
figure(2); imshow(imgP);
imwrite(imgP, fullfile('paper\planar_struct_SR_CVPR2015\figures\plane', 'planedetection.png'));

end
function vpData = sr_read_vpdata(fileName)

% SC_READ_VPDATA: read the data from vanishing point detection algorithm
% Input:
%   - fileName: the txt file containing the pre-computed vanishing point
%   detection code
% Output:
%   - vpData
%   The data structure of vpData
%       - vpData.numVP: number of detected vanishing points
%       - vpData.vp{i}.pos: the vanishing point position in the homogenous coordiante
%       - vpData.vp{i}.score: the score of the vanishing point
%       - vpData.vp{i}.numLines: number of lines supporting the vanishing point
%       - vpData.vp{i}.lines{j}.p1: (x1, y1): starting position
%       - vpData.vp{i}.lines{j}.p2: (x2, y2): ending position
%       - vpData.vp{i}.lines{j}.length: length of the line segment

vpData = [];

% Read data
fid = fopen(fileName);

%% Parse VP positions
temp = fscanf(fid, '%s ', [1 5]);
numVP = 0;
readVPFlag = 1;
VP = [];
while(readVPFlag)
    numVP = numVP + 1;
    vpCurr = fscanf(fid, '%g %g %g %g %g', [5 1]);
    if(~isempty(vpCurr))
        VP(:,numVP) = vpCurr;
    else
        temp = fscanf(fid, '%s ', [1 6]);
        readVPFlag = 0;
    end
end
VP = VP';

vpData.numVP = size(VP, 1);

% Save VP position data
for i = 1: vpData.numVP
    vpData.vp{i}.pos = VP(i, 1:3);
    vpData.vp{i}.score = VP(i, 4);
    vpData.vp{i}.numLines = VP(i, 5);
end

%% Parse each set of line segments for the corresponding VP
for i = 1: vpData.numVP
    numLine = fscanf(fid, '%d ', [1 1]);
    lines = fscanf(fid, '%g %g %g %g %g', [5 numLine]);
    vpData.vp{i}.lines = lines';
end

fclose(fid);
end

function modelPlane = sr_detect_plane_from_vp(vpData, img, opt)

% SC_DETECT_PLANE_FROM_VP: simple plane detection algorithm
% Input:
%     - vpData: vanishing point data
%     - img: input image
%     - mask: hole mask
% Output:
%     - modelPlane

%%

modelPlane = [];

% === Setting up ===
[imgH, imgW, ch] = size(img);
HfilterX = fspecial('gaussian', [1, opt.filterSize], opt.filterSigma);
HfilterY = HfilterX';
% fspecial('gaussian', opt.filterSize, opt.filterSigma);

img = im2double(img);

% === Supporting lines spatial support estimation ===
shapeInserter = vision.ShapeInserter('Shape', 'Lines','BorderColor', 'White');

for i = 1: vpData.numVP
    % The support lines
    imgLines = zeros(imgH, imgW);
    imgLines = step(shapeInserter, imgLines, int16(round(vpData.vp{i}.lines(:,1:4))));
    % Spatial density estimation via blurring
    imgLinesPosMap = imgLines;
    for k = 1:opt.numFilterIter
        imgLinesPosMap = imfilter(imgLinesPosMap, HfilterX, 'conv', 'replicate');
    end
    for k = 1:opt.numFilterIter
        imgLinesPosMap = imfilter(imgLinesPosMap, HfilterY, 'conv', 'replicate');
    end

    % Save results
    modelPlane.vp{i}.imgLines = imgLines;
    modelPlane.vp{i}.imgLinesPosMap = imgLinesPosMap;
end

% === Estimate plane support and plane parameters ===
numPlane = (vpData.numVP)*(vpData.numVP-1)/2;
% Initialize plane data
modelPlane.plane = cell(numPlane, 1);

indPlane = 1;
% A pair of vanishing points forms a plane hypothesis
for i = 1: vpData.numVP - 1
    for j = i+1: vpData.numVP
        % Compute the vanishing line
        modelPlane.plane{indPlane}.vLine = vLineFromTwoVP(vpData.vp{i}.pos, vpData.vp{j}.pos);
        % Element-wise product of two support line density
        modelPlane.plane{indPlane}.imgPlaneProb = modelPlane.vp{i}.imgLinesPosMap.*modelPlane.vp{j}.imgLinesPosMap; % Product of two probability maps

        %         modelPlane.plane{indPlane}.imgPlaneProb(mask) = 1e-10;
        modelPlane.plane{indPlane}.sourceVP = [i, j];

        indPlane = indPlane + 1;
    end
end

% === Compute rectified rotation parameters ===

for i = 1: numPlane
    for vpInd = 1: 2

        linesCurr = vpData.vp{modelPlane.plane{i}.sourceVP(vpInd)}.lines;
        invalidLineInd = linesCurr(:,5) == 0;
        linesCurr = linesCurr(~invalidLineInd,:);
        numLines = size(linesCurr, 1);

        vLineCurr = modelPlane.plane{i}.vLine;

        % Rectified homography
        H = eye(3);
        H(3,:) = vLineCurr;

    end
end

% === Add a fronto-parallel plane ===

modelPlane.plane{indPlane}.vLine = [0 0 1];
modelPlane.plane{indPlane}.imgPlaneProb = opt.fpPlaneProb*ones(imgH, imgW);
modelPlane.plane{indPlane}.score = sum(modelPlane.plane{indPlane}.imgPlaneProb(:));

numPlane = numPlane + 1;

modelPlane.numPlane = numPlane;

% === Compute posterior probability ===

planeProb = zeros(imgH, imgW, numPlane);
for i = 1 : numPlane
    planeProb(:,:,i) = modelPlane.plane{i}.imgPlaneProb;
end
planeProbSum = sum(planeProb, 3);
planeProb = bsxfun(@rdivide, planeProb, planeProbSum);
planeProb = planeProb + 0.1; % blur the posterior map
planeProb = bsxfun(@rdivide, planeProb, sum(planeProb, 3));

modelPlane.postProb = planeProb;

end

function vLine = vLineFromTwoVP(vp1, vp2)

A = cat(1, vp1, vp2);

[U S V] = svd(A, 0);
vLine = V(:,end);
vLine = vLine/vLine(3); % [h7, h8, 1]

end

author: thomas-koehler

repo: SupER


text:
function versionnumber = SFS_version()
%SFS_VERSION Return the version number of the SFS Toolbox
%
%   Usage: versionnumber = SFS_version;
%
%   SFS_VERSION returns the version number of the used Sound Field Synthesis
%   Toolbox.
%
%   See also: SFS_start, SFS_config

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

%% ===== Configuration ===================================================
% Version number
versionnumber = '2.5.0';

author: sfstoolbox

repo: sfs-matlab


text:
function ranks=cosmo_tiedrank(data, dim)
% Compute ranks for the input along the specified dimension
%
% ranks=cosmo_tiedrank(data[, dim])
%
% Inputs:
%   data                        numeric N-dimensional array
%   dim                         optional dimension along which the ranks
%                               are computed (default: 1)
%
% Output:
%   ranks                       numeric N-dimensional array with the same
%                               size as the input containing the rank of
%                               each vector along the dim-th dimension.
%                               Equal values have the same rank, which is
%                               the average of the rank the values would
%                               have if they differed by a minimal amount.
%                               NaN values in the input result in a NaN
%                               values in the output at the corresponding
%                               locations.
%                               If dim is greater than the number of
%                               dimensions in data, then all values in rank
%                               are one (or NaN of the corresponding value
%                               in data is NaN).
%
% Examples:
%     cosmo_tiedrank([1 2 2],2)
%     %|| [ 1 2.5 2.5]
%
%     cosmo_tiedrank([NaN 2 2;3 NaN 4],1)
%     %|| [ NaN     1     1;
%     %||     1   NaN     2];
%
%     cosmo_tiedrank([NaN 2 2;3 NaN 4],2)
%     %|| [ NaN   1.5   1.5;
%     %||     1   NaN     2];
%
%     cosmo_tiedrank([2 4 3 3 3 3 5 5 5],2)
%     %|| [ 1.0 6.0 3.5 3.5 3.5 3.5 8.0 8.0 8.0 ]
%
% Notes:
% - Unlike the Matlab builtin function 'tiedrank' (part of the statistics
%   toolbox), the meaning of the second argument is the dimension along
%   which the ranks are computed.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    if nargin<2
        dim=1;
    end

    check_inputs(data, dim);
    orig_size=size(data);

    if numel(orig_size)<dim
        % if the input data does not have enough data, then the output
        % consists of an array with only ones (or NaNs, if present)
        ranks=singleton_ranks(data);
        return;
    end

    [values,idx]=sort(data,dim);

    data_is_vector=numel(orig_size)<=2 && orig_size(3-dim)==1;
    if data_is_vector
        ranks=vector_tied_rank(values(:), idx(:));
        if orig_size(1)==1
            % transpose to turn it back into a row vector
            ranks=ranks';
        end
        return
    end

    % make the dim-th dimension the first dimension
    values_sh=shiftdim(values,dim-1);
    idx_sh=shiftdim(idx,dim-1);
    sh_size=size(values_sh);

    count_along_dim=size(values_sh,1);

    % reshape into a matrix
    values_mat=reshape(values_sh,count_along_dim,[]);
    idx_mat=reshape(idx_sh,count_along_dim,[]);

    % space for output
    ranks_mat=zeros(size(idx_mat));

    % compute for each column vector
    n_col=size(ranks_mat,2);
    for k=1:n_col
        ranks_mat(:,k)=vector_tied_rank(values_mat(:,k),idx_mat(:,k));
    end

    % put back in shape after shiftdim
    ranks_sh=reshape(ranks_mat,sh_size);

    % undo shiftdim
    unshift_count=numel(orig_size)-dim+1;
    ranks=reshape(shiftdim(ranks_sh,unshift_count),orig_size);

function ranks=vector_tied_rank(sorted_values, sort_idx)
% sorted_values and sort_idx are the output from 'sort'
% it is assumes that sorted_values is a vector
    n_values=numel(sorted_values);
    nan_msk=isnan(sorted_values);
    nan_count=sum(nan_msk);
    non_nan_count=numel(sorted_values)-nan_count;

    % first set ranks for values without ties
    ranks=sort_idx+NaN;
    ranks(sort_idx(1:non_nan_count))=1:non_nan_count;

    % now deal with ties
    tie_msk=sorted_values(2:end)==sorted_values(1:(end-1));
    tie_idx=find(tie_msk);
    tie_count=numel(tie_idx);

    k=0;
    while k<tie_count
        k=k+1;

        tie_start=tie_idx(k);
        tie_end=tie_start+1;

        while tie_end<n_values ...
                && sorted_values(tie_end)==sorted_values(tie_end+1)
            tie_end=tie_end+1;
            k=k+1;
        end

        tie_value=(tie_start+tie_end)/2;
        pos=tie_start+(0:(tie_end-tie_start));
        ranks(sort_idx(pos))=tie_value;
    end

function ranks=singleton_ranks(data)
    % all ranks are either NaN or 1
    ranks=ones(size(data));
    ranks(isnan(data))=NaN;

function check_inputs(data, dim)
    if ~isnumeric(data)
        error('First input must be numeric')
    end

    if ~(isnumeric(dim) ...
            && isscalar(dim) ...
            && round(dim)==dim ...
            && dim>0)
        error('Second argument must be numeric integer');
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function M = convectionTermSpherical1D(u)
% This function uses the central difference scheme to discretize a 1D
% convection term in the form \grad (u \phi) where u is a face vactor
% It is for a cylindrical coordinate in the r direction
%
% SYNOPSIS:
%   M = convectionTermCylindrical1D(u)
%
% PARAMETERS:
%	u   - FaceVariable
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nr = u.domain.dims(1);
G = 1:Nr+2;
DXe = u.domain.cellsize.x(3:end);
DXw = u.domain.cellsize.x(1:end-2);
DXp = u.domain.cellsize.x(2:end-1);
% rp = u.domain.cellcenters.x;
rf = u.domain.facecenters.x;

% define the vectors to stores the sparse matrix data
iix = zeros(3*(Nr+2),1);
jjx = zeros(3*(Nr+2),1);
sx = zeros(3*(Nr+2),1);

% reassign the east, west, north, and south velocity vectors for the
% code readability
ue = u.xvalue(2:Nr+1).*DXp./(DXp+DXe).*rf(2:Nr+1).^2./(1/3*(rf(2:Nr+1).^3-rf(1:Nr).^3));
uw = u.xvalue(1:Nr).*DXp./(DXp+DXw).*rf(1:Nr).^2./(1/3*(rf(2:Nr+1).^3-rf(1:Nr).^3));

% calculate the coefficients for the internal cells
AE = reshape(ue,Nr,1);
AW = reshape(-uw,Nr,1);
APx = reshape((ue.*DXe-uw.*DXw)./DXp,Nr,1);

% build the sparse matrix based on the numbering system
rowx_index = reshape(G(2:Nr+1),Nr,1); % main diagonal x
iix(1:3*Nr) = repmat(rowx_index,3,1);
jjx(1:3*Nr) = [reshape(G(1:Nr),Nr,1); ...
		reshape(G(2:Nr+1),Nr,1); reshape(G(3:Nr+2),Nr,1)];
sx(1:3*Nr) = [AW; APx; AE];

% build the sparse matrix
kx = 3*Nr;
M = sparse(iix(1:kx), jjx(1:kx), sx(1:kx), Nr+2, Nr+2);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function jed2 = jed_to_next_noon ( jed1 )

%*****************************************************************************80
%
%% JED_TO_NEXT_NOON converts a JED to the JED of the next noon.
%
%  Discussion:
%
%    This is primarily to make a fair test of the weekday routines,
%    which have trouble when the JED is at midnight.
%
%    Note that noon corresponds to an integral JED value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real JED1, the Julian Ephemeris Date.
%
%    Output, real JED2, the Julian Ephemeris Date
%    of the next noon.
%
  jed2 = round ( jed1 );
%
%  The integer part of JED1 is one of the two integers that
%  bracket JED1.  If it's the smaller one (which it should
%  be as long as JED1 is positive), make it the bigger one.
%
%  This correctly leaves undisturbed cases where JED1 is
%  already an integer, and where JED1 is negative (which
%  is not a case we expect to occur often).
%
  if ( jed2 < jed1 )
    jed2 = jed2 + 1.0;
  end

  return
end


text:
classdef prtRegressPlsda < prtRegress
    % prtRegressPlsda  Partial least squares discriminant regression
    %
    %    REGRESS = prtRegressPlsda returns a Partial least squares
    %    discriminant regressor
    %
    %    REGRESS = prtRegressPlsda(PROPERTY1, VALUE1, ...) constructs a
    %    prtRegressPlsda object REGRESS with properties as specified by
    %    PROPERTY/VALUE pairs.
    %
    %    A prtRegressPlsda object inherits all properties from the abstract
    %    Regress prtRegress. In addition is has the following properties:
    %
    %    nComponents  -  The number of components
    %    Bpls         -  The regression weights, estimated during training
    %    xMeans       -  The xMeans, estimated during training
    %    yMeans       -  The yMeana, estimated during training
    %
    %    For information on the partial least squares discriminant
    %    algorithm, please refer to the following URL:
    %
    %    http://en.wikipedia.org/wiki/Partial_least_squares_regression
    %
    %    A prtRegressPlsda object inherits the TRAIN, RUN, CROSSVALIDATE and
    %    KFOLDS methods from prtAction. It also inherits the PLOT method
    %    from prtRegress.
    %
    %    Example:
    %
    %   TestDataSet = prtDataGenNoisyLine;         % Create some test and
    %   TrainingDataSet = prtDataGenNoisyLine;     % training data
    %   regress = prtRegressPlsda;                % Create a regressor
    %   regress = regress.train(TrainingDataSet); % Train
    %   regressed = run(regress, TestDataSet);    % Test
    %   regress.plot;
    %
    %    See also prtRegress, prtRegressLslr, prtRegressRvm

    properties (SetAccess=private)
        name = 'Partial Least Squares Discriminant' % Partial Least Squares Discriminant
        nameAbbreviation = 'PLSDA' % PLSDA
        isNativeMary = true;  % True
    end

    properties
        % w is a DataSet.nDimensions x 1 vector of projection weights
        % learned during Fld.train(DataSet)
        nComponents = 2;
    end

    properties (SetAccess=protected)
        Bpls     % The prediction weights
        loadings % T
        xFactors % P
        yFactors % Q
        yMeansFactor % Factor to be added into regression output (accounts for X means and yMeans);
    end

    methods

        function Obj = prtRegressPlsda(varargin)
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end

        function Obj = set.nComponents(Obj,val)
            if ~prtUtilIsPositiveInteger(val)
                error('prt:prtRegressPlsda:nComponents','nComponents must be a positive integer');
            end
            Obj.nComponents = val;
        end

    end

    methods (Access=protected, Hidden = true)

        function Obj = trainAction(Obj,DataSet)

            X = DataSet.getObservations;

            Y = DataSet.getY;

            maxComps = min(size(X));
            if Obj.nComponents > maxComps;
                Obj.nComponents = maxComps;
            end

            xMeans = mean(X,1);
            yMeans = mean(Y,1);
            X = bsxfun(@minus, X, xMeans);
            Y = bsxfun(@minus, Y, yMeans);

            [Obj.Bpls, R, Obj.xFactors, Obj.yFactors, Obj.loadings, U] = prtUtilSimpls(X,Y,Obj.nComponents);  %#ok<ASGLU,NASGU>

            Obj.yMeansFactor = yMeans - xMeans*Obj.Bpls;

        end

        function DataSet = runAction(Obj,DataSet)
            yOut = bsxfun(@plus,DataSet.getObservations*Obj.Bpls, Obj.yMeansFactor);
            DataSet = DataSet.setObservations(yOut);
        end

        function xOut = runActionFast(Obj,xIn,ds) %#ok<INUSD>
           xOut = bsxfun(@plus,xIn*Obj.Bpls, Obj.yMeansFactor);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [bestEpsilon bestF1] = selectThreshold(yval, pval)
%SELECTTHRESHOLD Find the best threshold (epsilon) to use for selecting
%outliers
%   [bestEpsilon bestF1] = SELECTTHRESHOLD(yval, pval) finds the best
%   threshold to use for selecting outliers based on the results from a
%   validation set (pval) and the ground truth (yval).
%

bestEpsilon = 0;
bestF1 = 0;
F1 = 0;

stepsize = (max(pval) - min(pval)) / 1000;
for epsilon = min(pval):stepsize:max(pval)

    % ====================== YOUR CODE HERE ======================
    % Instructions: Compute the F1 score of choosing epsilon as the
    %               threshold and place the value in F1. The code at the
    %               end of the loop will compare the F1 score for this
    %               choice of epsilon and set it to be the best epsilon if
    %               it is better than the current choice of epsilon.
    %
    % Note: You can use predictions = (pval < epsilon) to get a binary vector
    %       of 0's and 1's of the outlier predictions

    predictions = (pval < epsilon);

    fp = sum((predictions == 1) & (yval == 0));
    fn = sum((predictions == 0) & (yval == 1));
    tp = sum((predictions == 1) & (yval == 1));

    prec = tp / (tp + fp);
    rec = tp / (tp + fn);

    F1 = 2 * prec * rec / (prec + rec);

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

    if F1 > bestF1
       bestF1 = F1;
       bestEpsilon = epsilon;
    end
end

end

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 9 Assignments/Anomaly Detection and Recommender Systems/mlclass-ex8/selectThreshold.m


text:
function point_merge_test05 ( m, n, n_unique, tol, seed )

%*****************************************************************************80
%
%% POINT_MERGE_TEST05 times uniqueness indexing with a tolerance.
%
%  Discussion:
%
%    POINT_RADIAL_TOL_UNIQUE_COUNT uses an algorithm that should be,
%      in general, O(N);
%    POINT_TOL_UNIQUE_COUNT uses an O(N^2) algorithm.
%
%    For this test, we just want to make sure the algorithms agree
%    in the counting.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POINT_MERGE_TEST05\n' );
  fprintf ( 1, '  We time the computations in TEST04, calling\n' );
  fprintf ( 1, '  POINT_RADIAL_TOL_UNIQUE_COUNT, (with random center)\n' );
  fprintf ( 1, '  POINT_TOL_UNIQUE_COUNT, (with zero tolerance)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  M = %d\n', m );
  fprintf ( 1, '  N = %d\n', n );
  fprintf ( 1, '  TOL = %f\n', tol );
  fprintf ( 1, '  SEED = %d\n', seed );

  [ a, seed ] = r8col_duplicates ( m, n, n_unique, seed );
%
%  The form of the tolerance test means that if two vectors are initially
%  equal, they remain "tolerably equal" after the addition of random
%  perturbation vectors whose 2-norm is no greater than TOL/2.
%
  for j = 1 : n
    [ r, seed ] = r8vec_uniform_01 ( m, seed );
    r(1:m) = r(1:m) / sqrt ( sum ( r(1:m).^2 ) );
    a(1:m,j) = a(1:m,j) + 0.5 * tol * r(1:m,1);
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  N_UNIQUE =                      %d\n', n_unique );

  tic;
  [ unique_num, undx, xdnu, seed ] = point_radial_tol_unique_index ( m, n, a, ...
    tol, seed );
  wtime = toc;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  POINT_RADIAL_TOL_UNIQUE_INDEX\n' );
  fprintf ( 1, '  Unique_num = %d\n', unique_num );
  fprintf ( 1, '  Time = %f\n', wtime );

  tic;
  [ unique_num, xdnu ] = point_tol_unique_index ( m, n, a, tol );
  wtime = toc;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  POINT_TOL_UNIQUE_INDEX\n' );
  fprintf ( 1, '  Unique_num = %d\n', unique_num );
  fprintf ( 1, '  Time = %f\n', wtime );

  return
end


text:
classdef Downloader < adi.utils.libad9361

    properties (Constant)
       PossibleDependencies = {'libad9361'};
    end

    methods

        function download(obj,dep)
            if nargin == 0
                fprintf('Specific dependency name on input. Options are:\n%s',obj.PossibleDependencies);
            end
            switch dep
                case 'libad9361'
                    download_libad9361(obj);
                otherwise
                    error('Unknown dependency %s',dep);
            end
        end

    end
end

author: analogdevicesinc