text:
function dat = read_spike6mat_data(filename, varargin)

% read_spike6mat_data() - read Matlab files exported from Spike 6
%
% Usage:
%   >> header = read_spike6mat_data(filename, varargin);
%
% Inputs:
%   filename - [string] file name
%
% Optional inputs:
%   'begsample'      first sample to read
%   'endsample'      last sample to read
%   'chanindx'  -    list with channel indices to read
%   'header'    - FILEIO structure header
%
% Outputs:
%   dat    - data over the specified range
% _______________________________________________________________________
% Copyright (C) 2008 Institute of Neurology, UCL
% Vladimir Litvak

if nargin < 1
    help read_spike6mat_data;
    return;
end;

header    = keyval('header',     varargin);
begsample = keyval('begsample',  varargin);
endsample = keyval('endsample',  varargin);
chanindx  = keyval('chanindx',   varargin);

if isempty(header)
    header = read_spike6mat_header(filename);
end

if isempty(begsample), begsample = 1; end;
if isempty(endsample), endsample = header.nSamples; end;

try
    vars = struct2cell(load(filename));
catch
    error('File not found or wrong format.');
end

if isempty(chanindx)
    chanindx = 1:numel(vars);
end

dat = zeros(length(chanindx), endsample-begsample+1);

for i = 1:length(chanindx)
    dat(i, :) = vars{chanindx(i)}.values(begsample:endsample);
end


text:
classdef TP5 < PROBLEM
% <multi> <real> <large/none> <robust>
% Test problem for robust multi-objective optimization
% delta --- 0.05 --- Maximum disturbance degree
% H     ---   50 --- Number of disturbances

%------------------------------- Reference --------------------------------
% A. Gaspar-Cunha, J. Ferreira, and G. Recio, Evolutionary robustness
% analysis for multi-objective optimization: benchmark problems, Structural
% and Multidisciplinary Optimization, 2014, 49: 771-793.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        delta;      % Maximum disturbance degree
        H;          % Number of disturbances
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            [obj.delta,obj.H] = obj.ParameterSet(0.05,50);
            obj.M = 2;
            if isempty(obj.D); obj.D = 10; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj(:,1) = PopDec(:,1);
            g = 1 + 10*mean(PopDec(:,2:end),2);
            h = sin(4*pi*PopDec(:,1))/15 - PopDec(:,1) + 1;
            PopObj(:,2) = h.*g;
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(~,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = sin(4*pi*R(:,1))/15 - R(:,1) + 1;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case {'Mean_IGD','Mean_HV','Worst_IGD','Worst_HV'}
                    score = feval(metName,Population,obj);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Perturb solutions multiple times
        function PopX = Perturb(obj,PopDec,N)
            if nargin < 3; N = obj.H; end
            Delta = repmat(obj.delta.*(obj.upper-obj.lower),N*size(PopDec,1),1);
            w     = UniformPoint(N,obj.D,'Latin');
            Dec   = 2*Delta.*w(reshape(repmat(1:end,size(PopDec,1),1),1,[]),:) + repmat(PopDec,N,1) - Delta;
            Dec   = obj.CalDec(Dec);
            PopX  = SOLUTION(Dec,obj.CalObj(Dec),obj.CalCon(Dec));
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [W] = winding_number(V,F,O,varargin)
  % WINDING_NUMBER Compute the sum of solid angles of a triangle (tetrahedron)
  % described by points (vectors) V
  %
  % [W] = winding_number(V,F,O)
  % [W] = winding_number(V,F,O,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %  V  #V by 3 list of vertex positions
  %  F  #F by 3 list of triangle indices
  %  O  #O by 3 list of origin positions
  %  Optional inputs:
  %    'Hierarchical'  followed by true or false. Use hierarchical evaluation.
  %      for mex: {true}, for matlab this is not supported
  %    'Fast' followed by whether to use fast-winding-number {false}
  %    'RayCast' followed by true or flase. Use ray cast version of approximate
  %      evaluation: {false}
  %    'TwoDRays' followed by true or false. Use 2d rays only.
  %    'NumRays' followed by the number of rays to cast for each origin
  % Outputs:
  %  W  no by 1 list of winding numbers
  %

  warning('not mex...');
  S = solid_angle(V,F,O);
  W = sum(S,2);
  switch size(F,2)
  case 3
    W = W/(2*pi);
  case 4
    W = W/(4*pi);
  end

end

author: alecjacobson

repo: gptoolbox


text:
clear; close all; clc;

group1=[10 12 22 15 21];
group2=[22 16 28 23] ;

mwwtest(group1,group2)

ns=4; nb=5;

mu_t=ns*(ns+nb+1)/2;
std_t=sqrt(ns*nb*(ns+nb+1)/12);

T1=17.5; T2=27.5;

Z_T=(T2-mu_t)/std_t

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function[varargout] = vmoment(varargin)
%VMOMENT Central moment over non-NaN elements along a specfied dimension.
%
%   Y=VMOMENT(X,N,DIM) finds the Nth central moment of all non-NaN elements
%   of X along dimension DIM.
%
%   [Y,NUM]=VMOMENT(X,N,DIM) also outputs the number of non-NaN data points
%   NUM, which has the same dimension as X.
%
%   [Y1,Y2,...YN]=VMOMENT(X1,X2,...XN,N,DIM) also works.
%
%   VMOMENT(X1,X2,...XN,N,DIM);  with no output arguments overwrites the
%   original input variables.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2001--2020 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
  vmoment_test,return
end

n=varargin{end-1};
ndim=varargin{end};

for i=1:length(varargin)-2
  x=varargin{i};
  m=vmean(x,ndim);
  m=vrep(m,size(x,ndim),ndim);

  %previously I had an abs around (x-m), that was incorrect
  [varargout{i},numi{i}]=vmean((x-m).^n,ndim);
end

for i=length(varargin)-1:nargout
  varargout{i}=numi{i-length(varargin)+2};
end

eval(to_overwrite(nargin-2))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[]=vmoment_test
x1=[1 2 3 nan];
x2=x1;
ans1=2/3;

vmoment(x1,x2,2,2);
reporttest('VMOMENT output overwrite', aresame(x1,ans1) && aresame(x2,ans1))

x1=[1 2 3 nan];
ans2=3;
[y1,y2]=vmoment(x1,2,2);
reporttest('VMOMENT moment & num', aresame(y1,ans1) && aresame(y2,ans2))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
%
% lellippi(phi, k, errtol)
%
% Inputs:
%
%   phi     Input angle vector size 1 or 1xN.
%   k       Input parameter vector size 1 or 1xN.
%   n       Input parameter vector size 1 or 1xN.
%   errtol  Error tolerance for Carlson's algorithms.
%
% Matlab function to compute Legendre's (incomplete) elliptic integral
% Pi(phi, k, n).  Uses a vectorized implementation of Carlson's Duplication Algorithms
% for symmetric elliptic integrals as found in "Computing Elliptic
% Integrals by Duplication," by B. C. Carlson, Numer. Math. 33, 1-16 (1979)
% and also found in ACM TOMS Algorithm 577.  Section 4 in the paper cited
% here describes how to convert between the symmetric elliptic integrals
% and Legendre's elliptic integrals.
%
% Returns NaN's for any argument values outside input range.
%

function f = lellippi(phi, k, n, errtol)
% Argument checking for vectorization:
lphi = length(phi);
lk = length(k);
ln = length(n);
errflag = logical(0);
if ( ~ ((lphi==lk) & (lphi==ln) & (lk==ln)) )
    if ( (lk==1) & (ln==1) )
        kvec = k * ones(1,lphi);
        nvec = n * ones(1,lphi);
        phivec = phi;
    elseif ( (lphi==1) & (ln==1) )
        phivec = phi * ones(1,lk);
        nvec = n * ones(1,lk);
        kvec = k;
    elseif ( (lphi==lk) & (ln==1) )
        nvec = n * ones(1,lphi);
        kvec = k;
        phivec = phi;
    elseif ( (lphi==1) & (lk==1) )
        phivec = phi * ones(1,ln);
        kvec = k * ones(1,lk);
        nvec = n;
    elseif ( (lphi==ln) & (lk==1) )
        kvec = k * ones(1,lphi);
        phivec = phi;
        nvec = n;
    elseif ( (lk==ln) & (lphi==1) )
        phivec = phi * ones(1,lk);
        kvec = k;
        nvec = n
    else
        disp('Incompatible input vector dimensions in lellipf!');
        errflag = logical(1);
    end
else
    phivec = phi;
    kvec = k;
    nvec = n;
end
if (~errflag)
    snphi = sin(phivec);
    csphi = cos(phivec);
    snphi2 = snphi.^2;
    csphi2 = csphi.^2;
    k2 = kvec.^2;
    y = 1.0 - k2.*snphi2;
    p = 1.0 + nvec .* snphi2;
    onesvec = ones(1,length(phivec));
    f = snphi .* rf(csphi2,  y, onesvec, errtol) - ...
        nvec .* snphi .* snphi2 .* rj(csphi2, y, onesvec, p, errtol) / 3.0;
else
    f = NaN;
end

author: Sable

repo: mcbench-benchmarks


text:
imageFile = 'data/einstein.png';
lipschitzConstant = 8;
totalVariationWeigth = 10;
noiseAmplitud = 100;
maxIter = 10;
optimizer = 'AcceleratedForwardBackward';

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function replaceLayer(obj, layerName, newLayerName, newBlock, addInputs, addOutputs, addParams, removeOldInputs)
% replaceLayer(obj, layerName, newLayerName, newBlock, [addInputs], [addOutputs], [addParams], [removeOldInputs])
%
% Takes a DAG and replaces an existing layer with a new one.
% If not specified, the inputs, outputs and params are reused from the
% previous layer.
%
% Copyright by Holger Caesar, 2015

% Default settings
if ~exist('removeOldInputs', 'var') || isempty(removeOldInputs),
    removeOldInputs = false;
end;

% Find the old layer
layerIdx = obj.getLayerIndex(layerName);
layer = obj.layers(layerIdx);

% Keep old settings (or not)
if removeOldInputs
    newInputs = {};
else
    newInputs = layer.inputs;
end
newOutputs = layer.outputs;
newParams  = layer.params;

% Add new settings
if exist('addInputs', 'var') && ~isempty(addInputs)
    newInputs = [newInputs, addInputs];
end
if exist('addOutputs', 'var') && ~isempty(addOutputs)
    newOutputs = [newOutputs, addOutputs];
end
if exist('addParams', 'var') && ~isempty(addParams)
    newParams = [newParams, addParams];
end

% Remove the old layer (must come before add to avoid conflicts)
obj.removeLayer(layerName);

% Add the new layer
obj.addLayer(newLayerName, newBlock, newInputs, newOutputs, newParams);

author: nightrome

repo: matconvnet-calvin


text:
 function ob = fatrix2_setup_dim(ob)
%function ob = fatrix2_setup_dim(ob)
%|
%| Set up idim and odim etc.
%|
%| Copyright 2010-12-17, Jeff Fessler, University of Michigan

persistent warned
if ~isvar('warned') || isempty(warned)
	warned = false;
end

arg = ob.arg; % this shall not be modified

if isempty(ob.imask) && isfield(arg, 'mask')
	ob.imask = arg.mask;
end

% input mask (rhs)

if isempty(ob.idim)
	if ~isempty(ob.imask)
		ob.idim = size(ob.imask);
	elseif isfield(arg, 'idim') && ~isempty(arg.idim)
		ob.idim = arg.idim;
	else
		fail('idim or imask or arg.idim required')
	end
end

if ~isempty(ob.imask)
	if ~islogical(ob.imask)
		fail('imask must be logical')
	end
	fatrix2_check_dim1(ob.imask, ob.idim)
	if isfield(arg, 'mask')
		if ~isequal(ob.imask, arg.mask)
			warn('mask inconsistent!?')
		end
	end
end

% output mask (lhs)

if isempty(ob.odim)
	if ~isempty(ob.omask)
		ob.odim = size(ob.omask);
	elseif isfield(arg, 'odim') && ~isempty(arg.odim)
		ob.odim = arg.odim;
	else
		fail('odim or omask or arg.odim required')
	end
end

if ~isempty(ob.omask)
	if ~islogical(ob.omask)
		fail('omask must be logical')
	end
	fatrix2_check_dim1(ob.omask, ob.odim)
end

ob.size = [prod(ob.odim) prod(ob.idim)];
if ~isempty(ob.imask)
	ob.size(2) = sum(ob.imask(:));
end
if ~isempty(ob.omask)
	ob.size(1) = sum(ob.omask(:));
end

% 1D input fully masked with >1D output is ambiguous
if ~ob.accept1d
	if numel(ob.idim) == 1 && numel(ob.odim) > 1 && ob.size(2) == prod(ob.idim)
%		fail 'use 1d odim for 1d idim case with full mask'
		warn '1d idim with full imask and >1d odim may not work'
	end

	if numel(ob.odim) == 1 && numel(ob.idim) > 1 && ob.size(1) == prod(ob.odim)
        if ~warned
            warn '1d odim full omask, >1d idim: transpose may not work'
            warned = true;
        end
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
clear

Extract_table_results_49
%%
scrsz = get(0,'ScreenSize');
figure1 = figure('Position',[20 50 3*scrsz(3)/4 0.9*scrsz(4)]);

set(figure1,'Units','Inches');
pos = get(figure1,'Position');
set(figure1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])

% Create axes
axes1 = axes('Parent',figure1,'FontSize',40,'FontName','Helvetica');

line_width = 6;
hold on;

[error_x, error_y] = cummErrorCurve(ceclm_error);
plot(error_x, error_y, 'r','DisplayName', 'OpenFace 2.0', 'LineWidth',line_width);

[error_x, error_y] = cummErrorCurve(clnf_error);
plot(error_x, error_y, 'DisplayName', 'OpenFace', 'LineWidth',line_width);

[error_x, error_y] = cummErrorCurve(cfss_error);
plot(error_x, error_y, 'DisplayName', 'CFSS', 'LineWidth',line_width);

[error_x, error_y] = cummErrorCurve(sdm_error);
plot(error_x, error_y, 'DisplayName', 'SDM', 'LineWidth',line_width);

[error_x, error_y] = cummErrorCurve(pocr_error);
plot(error_x, error_y,'DisplayName', 'PO-CR', 'LineWidth',line_width);

[error_x, error_y] = cummErrorCurve(cfan_error);
plot(error_x, error_y,'DisplayName', 'CFAN', 'LineWidth',line_width);

[error_x, error_y] = cummErrorCurve(tcdcn_error);
plot(error_x, error_y,'DisplayName', 'TCDCN', 'LineWidth',line_width);

[error_x, error_y] = cummErrorCurve(error_3ddfa);
plot(error_x, error_y,'DisplayName', '3DDFA', 'LineWidth',line_width);

set(gca,'xtick',[0:0.02:0.10])
xlim([0.02,0.10]);
xlabel('Size normalised MSE','FontName','Helvetica');
ylabel('Proportion of images','FontName','Helvetica');
grid on
% title('Fitting in the wild without outline','FontSize',60,'FontName','Helvetica');

leg = legend('show', 'Location', 'SouthEast');
set(leg,'FontSize',40)

print -dpdf results/Janus-no-outline.pdf
print -dpng results/Janus-no-outline.png

author: TadasBaltrusaitis

repo: OpenFace


text:
function arpack_test02 ( )

%*****************************************************************************80
%
%% ARPACK_TEST02 tests EIGS for the smallest eigenvalues, using a function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ARPACK_TEST02:\n' );
  fprintf ( 1, '  EIGS can compute a few eigenvalues\n' );
  fprintf ( 1, '  of least or greatest magnitude.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  When the matrix A is too large to store explicitly,\n' );
  fprintf ( 1, '  the user may supply a routine which computes A*x,\n' );
  fprintf ( 1, '  or inverse(A)*x.\n' );

  opts.disp = 0;
  opts.issym = 1;

  n = 139;
%
%  In general, n = size ( A, 1 ) where A = delsq ( numgrid('C',15) );
%
  lambda_small = eigs ( @(x)ainv_x(x,'C',15), n, 5, 'SM', opts )

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = laguerre ( n )

%*****************************************************************************80
%
%% LAGUERRE returns the LAGUERRE matrix.
%
%  Example:
%
%    N = 8
%
%      1      .     .      .    .    .    .    .
%      1     -1     .      .    .    .    .    .
%      2     -4     1      .    .    .    .    .   /    2
%      6    -18     9     -1    .    .    .    .   /    6
%     24    -96    72    -16    1    .    .    .   /   24
%    120   -600   600   -200   25   -1    .    .   /  120
%    720  -4320  5400  -2400  450  -36    1    .   /  720
%   5040 -35280 52920 -29400 7350 -882   49   -1   / 5040
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is lower triangular.
%
%    The columns of A alternate in sign.
%
%    A(I,1) = 1
%    A(I,I) = (-1)^(I-1) / (I-1)!.
%
%    LAMBDA(I) = (-1)^(I-1) / (I-1)!.
%
%    det ( A ) = product ( 1 <= I <= N ) (-1)^(I-1) / (I-1)!
%
%    The I-th row contains the coefficients of the Laguerre
%    polynomial of degree I-1.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  if ( n <= 0 )
    a = [];
    return
  end

  a(1,1) = 1.0;

  if ( n == 1 )
    return
  end

  a(2,1) = 1.0;
  a(2,2) = -1.0;

  if ( n == 2 )
    return
  end

  for i = 3 : n
    for j = 1 : n
      if ( j == 1 )
        a(i,j) = (  ( 2 * i - 3 ) * a(i-1,j) ...
                 +  (   - i + 2 ) * a(i-2,j) ) ...
                 /  (     i - 1 );
      else
        a(i,j) = (  ( 2 * i - 3 ) * a(i-1,j) ...
                 -  (         1 ) * a(i-1,j-1) ...
                 +  (   - i + 2 ) * a(i-2,j) ) ...
                 /  (     i - 1 );
      end
    end
  end

  return
end


text:
function r = uniform_dataset ( m, n, seed )

%*****************************************************************************80
%
%% UNIFORM_DATASET generates a uniform dataset and writes it to a file.
%
%  Discussion:
%
%    UNIFORM_DATASET generates a uniform random data and writes it to a file.
%
%  Usage:
%
%    r = uniform_dataset ( m, n, seed )
%
%    where
%
%    * M, the spatial dimension,
%    * N, the number of points to generate,
%    * SEED, the seed, a positive integer.
%    * R is the M by N array created.
%
%    creates an M by N uniform random dataset and writes it to the
%    file "uniform_M_N.txt".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 December 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'UNIFORM_DATASET\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Generate a uniform pseudorandom dataset.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program requests input values from the user:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * M, the spatial dimension,\n' );
  fprintf ( 1, '  * N, the number of points to generate,\n' );
  fprintf ( 1, '  * SEED, a seed for the random number generator.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program generates the data and writes it to the file\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    uniform_M_N.txt\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  where "M" and "N" are the numeric values.\n' );
%
%  Get the spatial dimension.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    m = input ( '  Enter the spatial dimension M: ' );
  else
    m = str2num ( m );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension M = %d\n', m );
%
%  Get the number of points.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    n = input ( '  Enter the number of points N: ' );
  else
    n = str2num ( n );
  end

  fprintf ( 1, '  Number of points N = %d\n', n );
%
%  Get the seed.
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    seed = input ( '  Enter the seed: ' );
  else
    seed = str2num ( seed );
  end

  fprintf ( 1, '  The seed = %d\n', seed );
%
%  Compute the data.
%
  [ r, seed ] = r8mat_uniform_01 ( m, n, seed );
%
%  Write it to a file.
%
  output_filename = ...
    strcat ( 'uniform_', num2str ( m ), '_', num2str ( n ), '.txt' );

  r8mat_write ( output_filename, m, n, r );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data was written to the file "%s".\n', ...
    output_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'UNIFORM_DATASET:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
    return;
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function [ r, seed ] = r8mat_uniform_01 ( m, n, seed )

%*****************************************************************************80
%
%% R8MAT_UNIFORM_01 returns a unit pseudorandom R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Second Edition,
%    Springer, 1987,
%    ISBN: 0387964673,
%    LC: QA76.9.C65.B73.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, December 1986, pages 362-376.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley, 1998,
%    ISBN: 0471134031,
%    LC: T57.62.H37.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, Number 2, 1969, pages 136-143.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, integer SEED, the integer "seed" used to generate
%    the output random number.
%
%    Output, real R(M,N), an array of random values between 0 and 1.
%
%    Output, integer SEED, the updated seed.  This would
%    normally be used as the input seed on the next call.
%
  i4_huge = 2147483647;

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8MAT_UNIFORM_01 - Fatal error!' );
  end

  for j = 1 : n
    for i = 1 : m

      seed = floor ( seed );

      seed = mod ( seed, i4_huge );

      if ( seed < 0 )
        seed = seed + i4_huge;
      end

      k = floor ( seed / 127773 );

      seed = 16807 * ( seed - k * 127773 ) - k * 2836;

      if ( seed < 0 )
        seed = seed + i4_huge;
      end

      r(i,j) = seed * 4.656612875E-10;

    end
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
% The code is to generate the memory base.
% Input: RGB images (img_path)
% Input: Label map (label_path)
% Output: segment, corresponding mask
%% output structure:
% library_img (training data): segmented rgb patches based on connectivity
%
% library_mask: connected segment mask, 1 indicates exist
% library_mask_pole: (optional)Special handling for cityscapes dataset,
% because many objects often separated by thin pole structures (pole,
% light, traffic sign),
% this is not very important just for sightly better results
% for other datasets, can be deleted.

%% Function: Get the corresponding segment based on connectivity analysis.

img_path = '../traindata/RGB256Full/';
label_path = '../traindata/label_refine/';
list = dir([img_path '*.png']);
list = struct2cell(list);
load('cityscapes_colormap.mat');
pixel_matrix = zeros(256,512,3,'uint8');

save_path_library = '../traindata/original_segment/';

mkdir(save_path_library);

img_size_h = 256;
img_size_w = 512;
mapping = load('mapping.mat');
mapping = mapping.mapping;
pole_matrix = zeros(img_size_h,img_size_w,3,'uint8');

%% think structure
% cityscape thin structure effect remove, we find in cityscapes, only
% single buding often split out by the pole, tranffic sign, and traffic
% light, segments with the same class that are separated by these thin
% structures are merged together
pole_matrix(:,:,1) = mapping(6,1);
pole_matrix(:,:,2) = mapping(6,2);
pole_matrix(:,:,3) = mapping(6,3);

sign_matrix = pole_matrix;
sign_matrix(:,:,1) = mapping(7,1);
sign_matrix(:,:,2) = mapping(7,2);
sign_matrix(:,:,3) = mapping(7,3);

light_matrix = pole_matrix;
light_matrix(:,:,1) = mapping(8,1);
light_matrix(:,:,2) = mapping(8,2);
light_matrix(:,:,3) = mapping(8,3);

for c = 1:19

    pixel_matrix = zeros(img_size_h,img_size_w,3,'uint8');
    pixel_matrix(:,:,1) = mapping(c,1);
    pixel_matrix(:,:,2) = mapping(c,2);
    pixel_matrix(:,:,3) = mapping(c,3);
    save_path_library_class = [save_path_library sprintf('%02d',c)];
    mkdir([save_path_library_class '/']);

    for i = size(list,2):size(list,2)
        i
        img = imread([img_path, list{1,i}]);
        img = im2double(img);
        label = imread([label_path list{1,i}]);
        mask = single(label)-single(pixel_matrix);
        mask = sum(mask.^2,3);
        mask(mask~=0) = -1;
        mask = mask + 1;
        mask_original = mask;
        %% sanity check for think structures, comment this if you do not use this
        % comment begin
        if(c~=6 & c~=7 & c~=8)
            mask_pole = sum((single(label)-single(pole_matrix)).^2,3);
            mask_pole(mask_pole~=0) = -1;
            mask_pole = mask_pole + 1;

            mask_light = sum((single(label)-single(light_matrix)).^2,3);
            mask_light(mask_light~=0) = -1;
            mask_light = mask_light + 1;

            mask_sign = sum((single(label)-single(sign_matrix)).^2,3);
            mask_sign(mask_sign~=0) = -1;
            mask_sign = mask_sign + 1;
            small_mask = mask_pole + mask_sign + mask_light;
            mask(mask_pole==1 | mask_light==1 | mask_sign == 1) = 1;
        end
        % comment end

        %% find connected component
        connected_component = bwlabel(mask);
        conn = unique(connected_component);

        library_img = zeros(size(img,1),size(img,2),3,size(conn,1),'uint8');
        library_mask = zeros(size(img,1),size(img,2),size(conn,1),'uint8');
        library_mask_pole = zeros(size(img,1),size(img,2),size(conn,1),'uint8');

        empty_ind = [];
        for j = 1:size(conn,1)
            [r, col] = find(connected_component==conn(j));
          if(c~=6 & c~=7 & c~=8)
            %  t = size(r,1) < 1000
            if((sum(sum(small_mask(sub2ind([img_size_h,img_size_w],r,col))))>=sum(sum(mask_original(sub2ind([img_size_h,img_size_w],r,col))))))
                empty_ind(end+1) = j;
                continue;
            end

          end
            if(mask(r(1),col(1))==1)

                tmp_mask = mask;
                tmp_mask((connected_component~=conn(j))) = 0;

                tmp_img = img.*repmat(tmp_mask,[1,1,3]);
                tmp_img = tmp_img.*repmat(mask_original,[1,1,3]);
                library_img(:,:,:,j) = uint8(tmp_img*255.0);
                library_mask(:,:,j) = uint8(tmp_mask.*mask_original);
                [row_lm, col_lm] = find(library_mask(:,:,j)==1);
                range_row = min(row_lm(:)):max(row_lm(:));
                range_col = min(col_lm(:)):max(col_lm(:));
                library_mask_pole(range_row,range_col,j) = uint8(tmp_mask(range_row,range_col));

            end
        end

        library_img(:,:,:,empty_ind) = [];

        library_mask(:,:,empty_ind) = [];
        library_mask_pole(:,:,empty_ind) = [];
        save([save_path_library_class,'/',list{1,i}(1:end-4),'.mat'],'library_img','library_mask','library_mask_pole');

    end
end

author: xjqicuhk

repo: SIMS

save_path: github-repos/MATLAB/xjqicuhk-SIMS


text:
function obj = coef_J_pure_pnp_func_new(in1,in2,in3,scale)
b1 = in1(:,1);
b2 = in1(:,2);
cx = in3(:,3);
cy = in3(:,4);
fx = in3(:,1);
fy = in3(:,2);
r1 = in2(:,1);
r2 = in2(:,2);
r3 = in2(:,3);
t4 = b1.*r3;
t5 = cx.*r3;
t6 = fx.*r1;
t2 = -t4+t5+t6;
t7 = b2.*r3;
t8 = cy.*r3;
t9 = fy.*r2;
t3 = -t7+t8+t9;
t10 = b1.*r2.*2.0;
t16 = cx.*r2.*2.0;
t11 = t10-t16;
t12 = b2.*r2.*2.0;
t13 = fy.*r3.*2.0;
t22 = cy.*r2.*2.0;
t14 = t12+t13-t22;
t15 = t3.^2;
t17 = b1.*r1.*2.0;
t18 = fx.*r3.*2.0;
t26 = cx.*r1.*2.0;
t19 = t17+t18-t26;
t20 = b2.*r1.*2.0;
t24 = cy.*r1.*2.0;
t21 = t20-t24;
t23 = fx.*r2.*t2.*4.0;
t25 = t21.*t3.*2.0;
t27 = t2.^2;
t28 = t11.*t2.*2.0;
t29 = t4-t5+t6;
t30 = t2.*t29.*2.0;
t31 = t7-t8+t9;
t32 = t3.*t31.*2.0;
t33 = t14.*t3.*2.0;
t34 = fy.*r1.*t14.*4.0;
t35 = -t17+t18+t26;
t36 = fy.*r1.*t3.*4.0;
t37 = fx.*r2.*t19.*4.0;
t38 = t11.^2;
t39 = t21.^2;
t40 = -t12+t13+t22;
t41 = fx.^2;
t42 = r2.^2;
t43 = fy.^2;
t44 = r1.^2;
t45 = fy.*r1.*t21.*4.0;
t46 = b1-cx;
t47 = b2-cy;
t48 = t19.*t2.*2.0;
t49 = fx.*r2.*t11.*4.0;
t50 = fx.*r2.*t29.*4.0;
t51 = fy.*r1.*t31.*4.0;
t52 = t30+t32-t41.*t42.*4.0-t43.*t44.*4.0;
t53 = fx.*r2.*t35.*4.0;
t54 = t29.^2;
t55 = fy.*scale.*t3.*2.0;
t56 = t3.*t47.*2.0;
t57 = fy.*r1.*t40.*4.0;
t58 = t11.*t35.*2.0;
t59 = t21.*t40.*2.0;
t60 = r1.*scale.*t43.*4.0;
t61 = fx.*r2.*t46.*4.0;
t62 = t21.*t31.*2.0;
t63 = t29.*t35.*2.0;
t64 = fy.*scale.*t21.*2.0;
t65 = t21.*t47.*2.0;
t66 = t31.^2;
t67 = fx.*scale.*t2.*2.0;
t68 = t2.*t46.*2.0;
t69 = t11.*t29.*2.0;
t70 = t31.*t40.*2.0;
t71 = t11.*t46.*2.0;
t72 = fx.*scale.*t29.*2.0;
t73 = fy.*scale.*t31.*2.0;
obj = [scale.*(t15+t27),-scale.*(t28+t33),scale.*(t25+t48),-scale.*(t23-fy.*r1.*t3.*4.0),scale.*(t15.*-2.0+t30+t38+t14.^2),scale.*(t23+t36-t11.*t19.*2.0-t14.*t21.*2.0),-scale.*(t25+t34-t2.*t35.*2.0-fx.*r2.*t11.*4.0),scale.*(t27.*-2.0+t32+t39+t19.^2),-scale.*(t28+t37-t3.*t40.*2.0-fy.*r1.*t21.*4.0),-scale.*t52,t67,t55,-scale.*(t56+t68),scale.*(t33-t11.*t29.*2.0),-scale.*(t25+t34+t49-t19.*t29.*2.0),-scale.*(t36+t50+t58-t14.*t21.*2.0),scale.*(t28+t37+t45-t14.*t31.*2.0),scale.*(t38.*2.0-t39.*2.0+t19.*t35.*2.0-t14.*t40.*2.0-t41.*t42.*8.0+t43.*t44.*8.0),-scale.*(t45+t53-t11.*t29.*2.0-t14.*t31.*2.0),fx.*scale.*t11.*-2.0,fy.*scale.*t14.*-2.0,scale.*(t71+t14.*t47.*2.0),-scale.*(t48-t21.*t31.*2.0),scale.*(t23+t51+t59-t11.*t19.*2.0),scale.*(t49+t57-t19.*t29.*2.0-t21.*t31.*2.0),fx.*scale.*t19.*2.0,t64,-scale.*(t65+t19.*t46.*2.0),scale.*(t50-t51),r2.*scale.*t41.*-4.0,t60,scale.*(t61-fy.*r1.*t47.*4.0),scale.*(t15+t54),-scale.*(t36-t50),scale.*(t25+t63),-scale.*t52,-scale.*(t45-t53+t69+t3.*t40.*2.0),scale.*(t32+t39-t54.*2.0+t35.^2),t72,-t55,scale.*(t56-t29.*t46.*2.0),-scale.*(t23-t51),-scale.*(t49-t57+t62+t2.*t35.*2.0),-scale.*(t50+t51+t58+t59),r2.*scale.*t41.*4.0,t60,-scale.*(t61+fy.*r1.*t47.*4.0),scale.*(t62-t63),fx.*scale.*t35.*2.0,-t64,scale.*(t65-t35.*t46.*2.0),scale.*(t27+t66),scale.*(t28+t70),scale.*(t30+t38-t66.*2.0+t40.^2),-t67,t73,scale.*(t68-t31.*t47.*2.0),scale.*(t69-t70),fx.*scale.*t11.*-2.0,fy.*scale.*t40.*2.0,scale.*(t71-t40.*t47.*2.0),scale.*(t54+t66),-t72,-t73,scale.*(t29.*t46.*2.0+t31.*t47.*2.0),scale.*t41,fx.*scale.*t46.*-2.0,scale.*t43,fy.*scale.*t47.*-2.0,scale.*(t46.^2+t47.^2)];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function varargout = get(tpl,action,varargin)
%TEMPLATE/GET Access data stored in a Template object
%  TPL = GET(TPL,ACTION,VARARGIN)
%     ACTION 'var'
%     ACTION 'undefined'

%  Copyright (C) 2003 Guillaume Flandin <Guillaume@artefact.tk>
%  $Revision: 1.0 $Date: 2003/05/05 22:19:51 $

narginchk(2,3);

switch lower(action)
	case 'var'
		narginchk(2,3);
		if nargin == 2
			varargout{1} = tpl.varvals;
		elseif iscellstr(varargin{1})
			varargout{1} = {};
			for i=1:length(varargin{1})
				key = find(ismember(tpl.varkeys,varargin{1}{i}));
				if isempty(key)
					%error('[Template] No such variable name.');
					varargout{1}{end+1} = '';
				else
					varargout{1}{end+1} = tpl.varvals{key};
				end
			end
		elseif ischar(varargin{1})
			varargout{1} = char(get(tpl,'var',cellstr(varargin{1})));
		else
			varargout{1} = '';
		end
	case 'undefined'
    narginchk(3,3);
		tpl = loadtpl(tpl,varargin{1});
		str = get(tpl,'var',varargin{1});
		varargout{1} = {};

		%%%%%%%%%%%%%%%%%%%%%%%% WIH REGEXP ONLY %%%%%%%%%%%%%%%%%%%%
		% [b, e] = regexp(str,'{[^ \t\r\n}]+}');
		% for i=1:length(b)
		% 	if ~any(ismember(tpl.varkeys,str(b(i)+1:e(i)-1)))
		% 		varargout{1}{end+1} = str(b(i)+1:e(i)-1);
		%	end
		% end
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	otherwise
		varargout{1} = finish(get(tpl,'var',action),tpl.unknowns);
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function ilt = euler_inversion_sym(f_s, t, M, P)

% ilt = euler_inversion_sym(f_s, t, [M], [P])
%
% Returns an approximation to the inverse Laplace transform of function
% handle f_s evaluated at each value in t (1xn) using the Euler method as
% summarized in the source below. This is a symbolic implementation capable
% of much greater accuracy than euler_inversion. For this reason, it takes
% an additional argument, P, the number of significant digits required by
% the calculation. In general, P should be about 0.6*M.
%
% f_s: Handle to function of s
% t:   Times at which to evaluate the inverse Laplace transformation of f_s
% M:   Optional number of terms to sum for each t (64 is a good guess);
%      highly oscillatory functions require higher M, but this can grow
%      unstable; see example_inversions.m for an example of stability [64]
% P:   Optional precision of calculation in significant digits [default 32]
%
% Requires the Symbolic Toolbox(TM).
%
% Abate, Joseph, and Ward Whitt. "A Unified Framework for Numerically
% Inverting Laplace Transforms." INFORMS Journal of Computing, vol. 18.4
% (2006): 408-421. Print.
%
% The paper is also online: http://www.columbia.edu/~ww2040/allpapers.html.
%
% Tucker McClure
% Copyright 2012, The MathWorks, Inc.

    % Make sure t is n-by-1.
    if size(t, 1) == 1
        t = t';
    elseif size(t, 2) > 1
        error('Input times, t, must be a vector.');
    end

    % Set M to 64 if user didn't specify an M.
    if nargin < 3, M = 32; end

    % Set P to the greater of the default from digits() or 0.6M.
    if nargin < 4, P = max(floor(0.6*M), digits()); end

    % Vectorized Talbot's algorithm

    % Binominal function
    bnml = @(n, z) factorial(n)/(factorial(z)*factorial(n-z));

    xi = sym([0.5, ones(1, M), zeros(1, M-1), 2^-sym(M)]);
    for k = 1:M-1
        xi(2*M-k + 1) = xi(2*M-k + 2) + 2^-sym(M) * bnml(sym(M), sym(k));
    end
    k = sym(0:2*M); % Iteration index
    beta = vpa(sym(M)*log(sym(10))/3 + 1i*pi*k, P);
    eta  = vpa((1-mod(k, 2)*2) .* xi, P);

    % Make a mesh so we can do this entire calculation across all k for all
    % given times without a single loop (it's faster this way).
    [beta_mesh, t_mesh] = meshgrid(beta, sym(t));
    eta_mesh = meshgrid(eta, t);

    % Finally, calculate the inverse Laplace transform for each given time.
    f_s_evals = arrayfun(f_s, beta_mesh./t_mesh, 'UniformOutput', false);
    f_s_evals = vpa(reshape([f_s_evals{:}], size(beta_mesh)), P);
    ilt = vpa(10^(sym(M)/3)./sym(t).*sum(eta_mesh.*real(f_s_evals), 2), P);

end


text:
% This script calculates full-reference and no-reference quality measures
% on the sequence length datasets.

%% Update MATLAB path.
addpath('../algorithms/SRAlgorithms');
addpath('../algorithms/qualityMeasures');

%% Evaluation settings.

% Path to input data and result directory.
inputDir = '../../data';
datasetConfig = 'sequenceLengthDatasets';
resultDir = '../../results/sequenceLengthDatasets';

% Parameters for this evaluation.
initEvaluationParametersForSequenceLength;

% The quality measures that are calculated.
measure_val = 1:length(qualityMeasures);

% Size of the border in pixels that is cropped from super-resolved image to
% evaluate the quality measures without distortion by boundary artifacts.
border = 3;

%% Start evaluation.
calculateQualityMeasures(inputDir, datasetConfig, resultDir, binningFactors, numberOfFrames, numberOfFrames_val, sliding_val, sr_method, binning_val, scenes_val, compressions(compress_val), measure_val, border);

author: thomas-koehler

repo: SupER


text:
function prob_test062 ( )

%*****************************************************************************80
%
%% TEST062 tests EXPONENTIAL_CDF, EXPONENTIAL_CDF_INV, EXPONENTIAL_PDF;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST062\n' );
  fprintf ( 1, '  For the Exponential CDF:\n' );
  fprintf ( 1, '  EXPONENTIAL_CDF evaluates the CDF.\n' );
  fprintf ( 1, '  EXPONENTIAL_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  EXPONENTIAL_PDF evaluates the PDF.\n' );

  a = 1.0;
  b = 2.0;

  check = exponential_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST062 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF parameter B = %14f\n', b );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed  ] = exponential_sample ( a, b, seed );

    pdf = exponential_pdf ( x, a, b );

    cdf = exponential_cdf ( x, a, b );

    x2 = exponential_cdf_inv ( cdf, a, b );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end


text:
function [meshes, fnames] = meshImportFreesurferSurfaces(subjid, hemi, surfaces)
% Import and save freesurfer surfaces in vistasoft form
% [meshes, fnames] = meshImportFreesurferSurfaces(subjid, [hemi], [surfaces], [vw])
%
% Freesurfer autosegmentation produces at least 4 meshes per subject:
% white, pial, inflated, and sphere. The white mesh is closest to what
% vistasoft BuildMesh produces, and we take this to be the base mesh,
% defining the initVertices for all meshes. We then create and one or more
% meshes for one or both hemispheres, and save them as .mat files in
% vistasoft format.
%
% Inputs:
%   subjid: Either the subjid recognized by freesurfer, or the path to the
%               freesurfer subject directory.
%   hemi:   'l', 'r', or 'b' (left, right, or both). [Default = 'b']
%   surfaces: cell array with one or more of 'white', 'pial', 'sphere',
%               'inflated'
%                [Default = {'white' 'pial' 'sphere' 'inflated'};]
%
% Outputs:
%   meshes: structured array with vistasoft compatible meshes
%   fnames: paths to the created meshes
%
% Example:
%   % Install Ernie freesurfer directory and ernie vistasession
%   forceOverwrite = false;
%   fssubjectsdir = getenv('SUBJECTS_DIR');
%   mrtInstallSampleData('anatomy/freesurfer', 'ernie', fssubjectsdir, forceOverwrite);
%   erniePRF = mrtInstallSampleData('functional', 'erniePRF', [], forceOverwrite);
%   cd(erniePRF);
%   % Create and save meshes
%   [meshes, fnames] = meshImportFreesurferSurfaces('ernie');

if notDefined('subjid'), help(mfilename); error('subjid is required'); end
if notDefined('hemi'), hemi = 'b'; end
if notDefined('surfaces'), surfaces = {'white' 'pial' 'sphere' 'inflated'}; end

% open a hidden view to determine location to store meshes
%vw = initHiddenInplane;

% Check for subject directory
if exist(subjid, 'dir')
    subjdir = subjid;
else
    fsPath = getenv('SUBJECTS_DIR');
    if ~exist(fsPath, 'dir'), error('Freesurfer subjects directory not found'); end
    subjdir = fullfile(fsPath, subjid);
end

if ~exist(subjdir, 'dir')
    error('Subject freesurfer directory %s not found', subjdir);
end

switch lower(hemi(1))
    case 'l', hs.fs = 'lh'; hs.vista = 'Left';
    case 'r', hs.fs = 'rh'; hs.vista = 'Right';
    case 'b'
        hs(1).fs = 'lh'; hs(1).vista = 'Left';
        hs(2).fs = 'rh'; hs(2).vista = 'Right';
    otherwise
        error('Hemi %s not recognized', hemi);
end

nummeshes = length(hs) * length(surfaces);
fnames = cell(1, nummeshes);

for h = 1:length(hs)

    % Read in the white mesh for this hemifield. This is the base mesh for
    % the hemifield. All other meshes will have the same initVertices, the
    % same curvature values (for coloration), and the same faces

    fsSurface  = fullfile(subjdir, 'surf', sprintf('%s.%s', hs(h).fs, 'white'));
    msh0         = fs_meshFromSurface(fsSurface);
    msh0.name    = sprintf('%s_white', hs(h).vista);

    % intialize the structured array meshes for output. This will contain
    % all the mesh structs created
    if h==1, meshes = msh0; end

    % Loop over surfaces for this hemisphere
    for ii = 1:length(surfaces)

        % Create this mesh and store it temporarily in the variable msh1
        fsSurface   = fullfile(subjdir, 'surf', sprintf('%s.%s', hs(h).fs, surfaces{ii}));
        msh1         = fs_meshFromSurface(fsSurface);
        msh1.name    = sprintf('%s_%s', hs(h).vista, surfaces{ii});

        % The mesh we will save is called msh, and combines values from the
        % white mesh (msh0) and the current mesh (msh1)
        msh = msh0;
        msh.vertices = msh1.vertices;
        msh.name = msh1.name;

        % Previously, we looked up the mesh directory stored in the mrVista
        %   view structure.
        % savepth = viewGet(vw, 'MeshDir', hs(h).vista);
        savepth = fullfile('./3DAnatomy', hs(h).vista, '3DMeshes');
        if ~exist(savepth, 'dir'), mkdir(savepth); end
        fname = fullfile(savepth, msh.name);
        save(fname, 'msh')

        % save for output
        meshnum = length(surfaces) * (h-1) + ii;
        meshes(meshnum) = msh;
        fnames{meshnum} = fname;

    end

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function value = normal_01_cdf_inv ( p )

%*****************************************************************************80
%
%% NORMAL_01_CDF_INV inverts the standard normal CDF.
%
%  Discussion:
%
%    The result is accurate to about 1 part in 10**16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 December 2004
%
%  Author:
%
%    Original FORTRAN77 version by Michael Wichura.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Michael Wichura,
%    The Percentage Points of the Normal Distribution,
%    Algorithm AS 241,
%    Applied Statistics,
%    Volume 37, Number 3, pages 477-484, 1988.
%
%  Parameters:
%
%    Input, real P, the value of the cumulative probability
%    densitity function.  0 < P < 1.  If P is not in this range, an "infinite"
%    result is returned.
%
%    Output, real VALUE, the normal deviate value with the
%    property that the probability of a standard normal deviate being
%    less than or equal to the value is P.
%
  a = [ 3.3871328727963666080,      1.3314166789178437745e+2, ...
        1.9715909503065514427e+3,   1.3731693765509461125e+4, ...
        4.5921953931549871457e+4,   6.7265770927008700853e+4, ...
        3.3430575583588128105e+4,   2.5090809287301226727e+3 ];
  b = [ 1.0,                        4.2313330701600911252e+1, ...
        6.8718700749205790830e+2,   5.3941960214247511077e+3, ...
        2.1213794301586595867e+4,   3.9307895800092710610e+4, ...
        2.8729085735721942674e+4,   5.2264952788528545610e+3 ];
  c = [ 1.42343711074968357734,     4.63033784615654529590, ...
        5.76949722146069140550,     3.64784832476320460504, ...
        1.27045825245236838258,     2.41780725177450611770e-1, ...
        2.27238449892691845833e-2,  7.74545014278341407640e-4 ];
  const1 = 0.180625;
  const2 = 1.6;
  d = [ 1.0,                        2.05319162663775882187,    ...
        1.67638483018380384940,     6.89767334985100004550e-1, ...
        1.48103976427480074590e-1,  1.51986665636164571966e-2, ...
        5.47593808499534494600e-4,  1.05075007164441684324e-9 ];
  e = [ 6.65790464350110377720,     5.46378491116411436990,    ...
        1.78482653991729133580,     2.96560571828504891230e-1, ...
        2.65321895265761230930e-2,  1.24266094738807843860e-3, ...
        2.71155556874348757815e-5,  2.01033439929228813265e-7 ];
  f = [ 1.0,                        5.99832206555887937690e-1, ...
        1.36929880922735805310e-1,  1.48753612908506148525e-2, ...
        7.86869131145613259100e-4,  1.84631831751005468180e-5, ...
        1.42151175831644588870e-7,  2.04426310338993978564e-15 ];
  split1 = 0.425;
  split2 = 5.0;

  if ( p <= 0.0 )
    value = -r8_huge ( );
    return
  end

  if ( 1.0 <= p )
    value = r8_huge ( );
    return
  end

  q = p - 0.5;

  if ( abs ( q ) <= split1 )

    r = const1 - q * q;
    value = q * r8poly_value ( 8, a, r ) / r8poly_value ( 8, b, r );

  else

    if ( q < 0.0 )
      r = p;
    else
      r = 1.0 - p;
    end

    if ( r <= 0.0 )

      value = r8_huge ( );

    else

      r = sqrt ( -log ( r ) );

      if ( r <= split2 )

        r = r - const2;
        value = r8poly_value ( 8, c, r ) / r8poly_value ( 8, d, r );

       else

         r = r - split2;
         value = r8poly_value ( 8, e, r ) / r8poly_value ( 8, f, r );

      end

    end

    if ( q < 0.0 )
      value = -value;
    end

  end

  return
end


text:
%> @file classDocumentationExample.m
%> @brief File used to show an example of class description
% ======================================================================
%> @brief Here we have a brief description of the class.
%
%> And here we can put some more detailed informations about the class.
% ======================================================================
classdef  (InferiorClasses = {?class1,?class2}) classDocumentationExample

  properties (Access = protected)
    %> Description of a protected property
    protectedProperty
  end
  properties (Access = public)
    %> Description of a public property
    publicProperty
  end
  properties (Access = private)
    %> Description of a private property
    pivateProperty
  end
  properties (Constant = true)
    %> Description of a constant property
    constantProperty = {'1', '2', 'trois'};
  end
  properties
    %> Description of the first property of the class
    first_property = []
    %> Description of the second property of the class
    second_property = []
    %> Description of the third property of the class
    third_property = [1  2];
  end
  events
    %> Description of first event
    FirstEvent
    %> Description of second event
    SecondEvent
  end
  %> Description of the enumeration.
  enumeration
    %> Description of the first item
    one (1)
    %> Description of the second item
    two (2)
    %> Description of the third item
    three
  end
  methods
    % ======================================================================
    %> @brief Class constructor
    %>
    %> More detailed description of what the constructor does.
    %>
    %> @param param1 Description of first parameter
    %> @param anotherParam Description of the second parametere
    %>
    %> @return instance of the classDocumentationExample class.
    % ======================================================================
    function obj = classDocumentationExample(param1, anotherParam)
    end

    % ======================================================================
    %> @brief Brief description of the exampleMethod1 method
    %>
    %> @param obj instance of the classDocumentationExample class.
    % ======================================================================
    function exampleMethod1(obj)
    end

    % ======================================================================
    %> @brief Brief description of the exampleMethod2 method
    %>
    %> @param obj instance of the classDocumentationExample class.
    %> @retval ret return value of this method
    % ======================================================================
    function ret = exampleMethod2(obj)
    end
  end
  methods (Static=true)
    % ======================================================================
    %> @brief Brief description of the exampleStaticMethod method
    %>
    %> This method is static and public, with an inused (~) argument
    %> @param param1 Description of the parameter
    %> @param param2 Description of the parameter
    %> @retval out return value of this method
    % ======================================================================
    function out = exampleStaticMethod(param1, ~, param2)
    end
  end
  methods (Static, Access=private)
    % ======================================================================
    %> @brief Brief description of the exampleStaticPrivateMethod method
    %>
    %> This method is static and private
    %> @param param1 Description of the parameter
    %> @param param2 Description of the parameter
    %> @retval out return value of this method
    % ======================================================================
    function out = exampleStaticPrivateMethod(param1, param2)
    end
  end
  methods (Access=protected, Static)
    % ======================================================================
    %> @brief Brief description of the exampleStaticProtectedMethod method
    %>
    %> This method is static and protected
    %> @param param1 Description of the parameter
    %> @retval out return value of this method
    % ======================================================================
    function out = exampleStaticProtectedMethod(param1)
    end
  end
  methods (Access=public, Static = true)
    % ======================================================================
    %> @brief Brief description of the exampleStaticPublicMethod method
    %>
    %> This method is static and public
    %> @param param1 Description of the parameter
    %> @retval out return value of this method
    % ======================================================================
    function out = exampleStaticPublicMethod(param1)
    end
  end
  methods (Access=private)
    % ======================================================================
    %> @brief Brief description of the examplePrivateMethod method
    %>
    %> This method is private
    %> @param param1 Description of the parameter
    %> @retval out return value of this method
    % ======================================================================
    function out = examplePrivateMethod(param1)
    end
  end
  methods (Access=protected)
    % ======================================================================
    %> @brief Brief description of the exampleProtectedMethod method
    %>
    %> This method is protected
    %> @param param1 Description of the parameter
    %> @retval out return value of this method
    % ======================================================================
    function out = exampleProtectedMethod(param1)
    end
  end
  methods (Access=public, Static=false)
    % ======================================================================
    %> @brief Brief description of the examplePublicMethod2 method
    %>
    %> This method is public and not static
    %> @param param1 Description of the parameter
    %> @retval out return value of this method
    % ======================================================================
    function out = examplePublicMethod2(param1)
    end
  end
  methods (Access=public, ~Static)
    % ======================================================================
    %> @brief Brief description of the exampleNonStaticPublicMethod3 method
    %>
    %> This method is public and not static
    %> @param param1 Description of the parameter
    %> @retval out return value of this method
    % ======================================================================
    function out = exampleNonStaticPublicMethod3(param1)
    end
  end
  methods (Abstract = true)
    % ======================================================================
    %> @brief Brief description of the exampleAbstractMethod method
    %>
    %> This method is abstract : only the signature of this function is
    %> declared.
    %> @param param1 Description of the first parameter
    %> @param param2 Description of the second parameter
    %> @retval out return value of this method
    % ======================================================================
    out = exampleAbstractMethod(param1, param2);
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
% statcond()  - compare two or more data conditions statistically using
%               standard parametric or nonparametric permutation-based ANOVA
%               (1-way or 2-way) or t-test methods. Parametric testing uses
%               fcdf() from the Matlab Statistical Toolbox.
% Usage:
%          >> [stats, df, pvals, surrog] = statcond( data, 'key','val'... );
%
% Inputs:
%   data       = one-or two-dimensional cell array of data matrices.
%                   For nonparametric, permutation-based testing, the
%                last dimension of the data arrays (which may be of up to
%                4 dimensions) is permuted across conditions, either in
%                a 'paired' fashion (not changing the, e.g., subject or
%                trial order in the last dimension) or in an umpaired
%                fashion (not respecting this order). If the number of
%                elements in the last dimension is not the same across
%                conditions, the 'paired' option is turned 'off'.  Note:
%                All other dimensions MUST be constant across conditions.
%                   For example, consider a (1,3) cell array of matrices
%                of size (100,20,x) each holding a (100,20) time/frequency
%                transform from each of x subjects. Only the last dimension
%                (here x, the number of subjects) may differ across the
%                three conditions.
%                   The test used depends on the size of the data array input.
%                When the data cell array has 2 columns and the data are
%                paired, a paired t-test is performed; when the data are
%                unpaired, an unpaired t-test is performed. If 'data'
%                has only one row (paired or unpaired) and more than 2
%                columns, a one-way ANOVA is performed. If the data cell
%                array contains several rows and columns, and the data is
%                paired, a two-way repeated measure ANOVA is performed.
%                NOTE THAT IF THE DATA is unpaired, EEGLAB will use a
%                balanced 1 or 2 way ANOVA and parametric results might not
%                be meaningful (bootstrap and permstatcondutation should be fine).
%
% Optional inputs:
%   'paired'   = ['on'|'off'] pair the data array {default: 'on' unless
%                the last dimension of data array is of different lengths}.
%                For two independent variables, this input is a cell array,
%                for example { 'on' 'off' } indicating that the first
%                independent variable is paired and the second is not.
%   'method'   = ['perm'|'bootstrap'|'param'] method for computing the p-values:
%                 'param' or 'parametric' = parametric testing (standard ANOVA
%                                           or t-test);
%                 'perm' or 'permutation' = non-parametric testing using
%                                           surrogate data
%                 'bootstrap' = non-parametric bootstrap
%                  made by permuting the input data {default: 'param'}
%   'naccu'    = [integer] Number of surrogate data copies to use in 'perm'
%                 or 'bootstrap' method estimation (see above) {default: 200}.
%   'verbose'  = ['on'|'off'] print info on the command line {default: 'on'}.
%   'variance' = ['homegenous'|'inhomogenous'] this option is exclusively
%                for parametric statistics using unpaired t-test. It allows
%                to compute a more accurate value for the degree of freedom
%                using the formula for inhomogenous variance (see
%                ttest2_cell function). Default is 'inhomegenous'.
%   'surrog'   = surrogate data array (see output).
%   'stats'    = F- or T-value array (see output).
%   'tail'     = ['one'|'two'] run one-tailed (F-test) or two tailed
%                (T-test). This option is only relevant when using the
%                'surrog' input. Otherwise it is ignored.
%   'forceanova' = ['on'|'off'] force the use of ANOVA calculation even
%                for 2x1 designs. Default is 'off'.
%   'alpha'    = [float] p-value threshold value. Allow returning
%                confidence intervals and mask (requires structoutput below).
%   'structoutput' = ['on'|'off'] return an output structure instead of
%                the regular output. Allow to output mask and confidence
%                intervals.
%
% Legacy parameters:
%   'threshold' - now 'alpha'
%   'mode'      - now 'method'
%
% Outputs:
%   stats      = F- or T-value array of the same size as input data without
%                the last dimension. A T value is returned only when the data
%                includes exactly two conditions.
%   df         = degrees of freedom, a (2,1) vector, when F-values are returned
%   pvals      = array of p-values. Same size as input data without the last
%                data dimension. All returned p-values are two-tailed.
%   surrog     = surrogate data array (same size as input data with the last
%                dimension filled with a number ('naccu') of surrogate data sets.
%
% Important note: When a two-way ANOVA is performed, outputs are cell arrays
%                 with three elements: output(1) = row effects;
%                 output(2) = column effects; output(3) = interactions
%                 between rows and columns.
%
% Examples:
%      >> a = { rand(1,10) rand(1,10)+0.5 }; % pseudo 'paired' data vectors
%         [t df pvals] = statcond(a);        % perform paired t-test
%           pvals =
%              5.2807e-04 % standard t-test probability value
%         % Note: for different rand() outputs, results will differ.
%
%         [t df pvals surog] = statcond(a, 'method', 'perm', 'naccu', 2000);
%           pvals =
%              0.0065 % nonparametric t-test using 2000 permuted data sets
%
%         a = { rand(2,11) rand(2,10) rand(2,12)+0.5 }; % pseudo 'unpaired'
%         [F df pvals] = statcond(a); % perform an unpaired ANOVA
%           pvals =
%              0.00025 % p-values for difference between columns
%              0.00002 % for each data row
%
%         a = { rand(3,4,10) rand(3,4,10) rand(3,4,10); ...
%               rand(3,4,10) rand(3,4,10) rand(3,4,10)+0.5 };
%         % pseudo (2,3)-condition data array, each entry containing
%         %                                    ten (3,4) data matrices
%         [F df pvals] = statcond(a);  % perform a paired 2-way ANOVA
%         % Output:
%           pvals{1} % a (3,4) matrix of p-values; effects across rows
%           pvals{2} % a (3,4) matrix of p-values; effects across colums
%           pvals{3} % a (3,4) matrix of p-values; interaction effects
%                                      % across rows and columns
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005-
%         With thanks to Robert Oostenveld for fruitful discussions
%         and advice on this function.
%
% See also: anova1_cell(), anova2_cell(), anova2rm_cell, fcdf()

% perform a paired t-test
% -----------------------
% a = { rand(2,10) rand(2,10) };
% [t df pval] = statcond(a); pval
% [h p t stat] = ttest( a{1}(1,:), a{2}(1,:)); p
% [h p t stat] = ttest( a{1}(2,:), a{2}(2,:)); p
%
% compare significance levels
% --------------------------
% a = { rand(1,10) rand(1,10) };
% [F df pval] = statcond(a, 'method', 'perm', 'naccu', 200); pval
% [h p t stat] = ttest( a{1}(1,:), a{2}(1,:)); p

% Copyright (C) Arnaud Delorme
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [ ori_vals, df, pvals, surrogval ] = statcond( data, varargin );

    if nargin < 1
        help statcond;
        return;
    end;
    try, warning('off', 'MATLAB:divideByZero'); catch, end;

    if exist('finputcheck')
        g = finputcheck( varargin, { 'naccu'      'integer'   [1 Inf]             200;
                                     'method'     'string'    { 'param','parametric','perm','permutation','bootstrap' }  'param';
                                     'mode'       'string'    { }                 '';
                                     'paired'     'string'    { 'on','off' }      'on';
                                     'surrog'     { 'real','cell' }      []       [];
                                     'stats'      { 'real','cell' }      []       [];
                                     'structoutput' 'string'  { 'on','off' }      'off';
                                     'forceanova'   'string'  { 'on','off' }      'off';
                                     'arraycomp'  'string'    { 'on','off' }      'on';
                                     'alpha'      'real'      []                  NaN;
                                     'tail'       'string'    { 'one','both','upper','lower'}    'both';
                                     'variance'   'string'    { 'homogenous','inhomogenous' }    'inhomogenous';
                                     'returnresamplingarray' 'string'    { 'on','off' }      'off';
                                     'verbose'    'string'    { 'on','off' }      'on' }, 'statcond');
        if isstr(g), error(g); end;
    else
        g = struct(varargin{:});
        if ~isfield(g, 'naccu'),     g.naccu = 200; end;
        if ~isfield(g, 'method'),    g.method  = 'param'; end;
        if ~isfield(g, 'paired'),    g.paired = 'on'; end;
        if ~isfield(g, 'surrog'),    g.surrog = []; end;
        if ~isfield(g, 'orivals'),   g.orivals = []; end;
        if ~isfield(g, 'arraycomp'), g.arraycomp = 'on'; end;
        if ~isfield(g, 'verbose'),   g.verbose = 'on'; end;
        if ~isfield(g, 'tail'),      g.tail = 'both'; end;
        if ~isfield(g, 'variance'),  g.variance = 'homogenous'; end;
        if ~isfield(g, 'structoutput'), g.structoutput = 'on'; end;
        if ~isfield(g, 'returnresamplingarray'),   g.returnresamplingarray = 'off'; end;
    end;
    if ~isempty(g.mode), g.method = g.mode; end;

    if strcmpi(g.method, 'parametric'), g.method = 'param'; end;
    if strcmpi(g.method, 'permutation'), g.method = 'perm'; end;
    if strcmpi(g.verbose, 'on'), verb = 1; else verb = 0; end;
    if strcmp(g.method, 'param' ) && exist('fcdf') ~= 2
      myfprintf('on',['statcond(): parametric testing requires fcdf() \n' ...
               '            from the Matlab StatsticaL Toolbox.\n' ...
               '            Running nonparametric permutation tests\n.']);
      g.method = 'perm';
    end
    if size(data,2) == 1, data  = transpose(data); end; % cell array transpose
    g.naccu = round(g.naccu);

    % reshape matrices
    % ----------------
    nd = size(data{1});
    nd = nd(1:end-1);
    for index = 1:prod(size(data))
        data{index} = reshape(data{index}, [prod(nd) size(data{index},myndims(data{index}))]);
    end;

    if ~strcmpi(g.method, 'param') && isempty(g.surrog)
         tmpsize   = size(data{1});
         surrogval = zeros([ tmpsize(1:end-1) g.naccu ], 'single');
    else surrogval = [];
    end;

    % check for NaNs or Inf
    % ---------------------
    for iDat = 1:length(data(:))
        if any(isnan(reshape(data{iDat}, prod(size(data{iDat})),1))) || ...
                any(isinf(reshape(data{iDat}, prod(size(data{iDat})),1)))
            error('Statcond: One of the input array contains NaNs or Infinite values');
        end;
    end;

    % bootstrap flag
    % --------------
    if strcmpi(g.method, 'bootstrap'), bootflag = 1;
    else                               bootflag = 0;
    end;

    if isempty(g.surrog)
        % test if data can be paired
        % --------------------------
        if length(unique(cellfun('size', data, ndims(data{1}) ))) > 1
            g.paired = 'off';
        end;
        if strcmpi(g.paired, 'on')
             pairflag = 1;
        else pairflag = 0;
        end;

        % return resampling array
        % -----------------------
        if strcmpi(g.returnresamplingarray, 'on')
            [ datavals datalen datadims ] = concatdata( data );
            if strcmpi(g.arraycomp, 'on')
                ori_vals = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
            else
                ori_vals = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
            end;
            return;
        end;

        % text output
        % -----------
        myfprintf(verb,'%d x %d, ', size(data,1), size(data,2));
        if strcmpi(g.paired, 'on')
             myfprintf(verb,'paired data, ');
        else myfprintf(verb,'unpaired data, ');
        end;
        if size(data,1) == 1 && size(data,2) == 2
             myfprintf(verb,'computing T values\n');
        else myfprintf(verb,'computing F values\n');
        end;
        if size(data,1) > 1
            if strcmpi(g.paired, 'on')
                 myfprintf(verb,'Using 2-way repeated measure ANOVA\n');
            else myfprintf(verb,'Using balanced 2-way ANOVA (not suitable for parametric testing, only bootstrap)\n');
            end;
        elseif size(data,2) > 2
            if strcmpi(g.paired, 'on')
                 myfprintf(verb,'Using 1-way repeated measure ANOVA\n');
            else myfprintf(verb,'Using balanced 1-way ANOVA (equivalent to Matlab anova1)\n');
            end;
        else
            if strcmpi(g.paired, 'on')
                 myfprintf(verb,'Using paired t-test\n');
            else myfprintf(verb,'Using unpaired t-test\n');
            end;
        end;
        if ~strcmpi(g.method, 'param')
            if bootflag, myfprintf(verb,'Bootstraps (of %d):', g.naccu);
            else         myfprintf(verb,'Permutations (of %d):', g.naccu);
            end;
        end;
    end;

    tail = g.tail;
    if isempty(g.surrog)
        if size(data,1) == 1, % only one row

            if size(data,2) == 2 && strcmpi(g.forceanova, 'off')

                % paired t-test (very fast)
                % -------------
                [ori_vals df] = ttest_cell_select(data, g.paired, g.variance);

                if strcmpi(g.method, 'param')

                    % Check if exist tcd.m file from the Statistics Toolbox (Bug 1352 )
                    if exist('tcdf','file') == 2  & license('test', 'Statistics_Toolbox')
                        pvals = 2*tcdf(-abs(ori_vals), df);
                    else
                        pvals = 2*mytcdf(-abs(ori_vals), df);
                    end

                    pvals = reshape(pvals, size(pvals));
                else
                    if strcmpi(g.arraycomp, 'on')
                        try
                            myfprintf(verb,'...');
                            res = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
                            surrogval = ttest_cell_select( res, g.paired, g.variance);
                        catch,
                           lasterr
                           myfprintf(verb,'\nSuperfast array computation failed because of memory limitation, reverting to standard computation');
                           g.arraycomp = 'off';
                        end;
                    end;
                    if strcmpi(g.arraycomp, 'off')
                        [res precomp] = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
                        for index = 1:g.naccu
                            res = surrogdistrib( {}, 'precomp', precomp);
                            if mod(index, 10) == 0, myfprintf(verb,'%d ', index); end;
                            if mod(index, 100) == 0, myfprintf(verb,'\n'); end;
                            if myndims(res{1}) == 1
                                 surrogval(index)     = ttest_cell_select(res, g.paired, g.variance);
                            else surrogval(:,index)   = ttest_cell_select(res, g.paired, g.variance);
                            end;
                        end;
                    end;
                end;
            else
                % one-way ANOVA (paired) this is equivalent to unpaired t-test
                % -------------
                tail = 'one';
                [ori_vals df] = anova1_cell_select( data, g.paired );
                if strcmpi(g.method, 'param')
                    pvals = 1-fcdf(ori_vals, df(1), df(2));
                else
                    if strcmpi(g.arraycomp, 'on')
                        try
                            myfprintf(verb,'...');
                            res = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
                            surrogval = anova1_cell_select( res, g.paired );
                        catch,
                            myfprintf(verb,'\nSuperfast array computation failed because of memory limitation, reverting to standard computing');
                            g.arraycomp = 'off';
                        end;
                    end;
                    if strcmpi(g.arraycomp, 'off')
                        [res precomp] = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
                        for index = 1:g.naccu
                            if mod(index, 10) == 0, myfprintf(verb,'%d ', index); end;
                            if mod(index, 100) == 0, myfprintf(verb,'\n'); end;

                            res = surrogdistrib( {}, 'precomp', precomp);
                            if myndims(data{1}) == 1
                            	 surrogval(index)     = anova1_cell_select( res, g.paired );
                            else surrogval(:,index)   = anova1_cell_select( res, g.paired );
                            end;
                        end;
                    end;
                end;
            end;
        else
            % two-way ANOVA (paired or unpaired)
            % ----------------------------------
            tail = 'one';
            [ ori_vals{1} ori_vals{2} ori_vals{3} df{1} df{2} df{3} ] = anova2_cell_select( data, g.paired );
            if strcmpi(g.method, 'param')
                pvals{1} = 1-fcdf(ori_vals{1}, df{1}(1), df{1}(2));
                pvals{2} = 1-fcdf(ori_vals{2}, df{2}(1), df{2}(2));
                pvals{3} = 1-fcdf(ori_vals{3}, df{3}(1), df{3}(2));
            else
                surrogval = { surrogval surrogval surrogval };
                dataori   = data;
                if strcmpi(g.arraycomp, 'on')
                    try
                        myfprintf(verb,'...');
                        res = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
                        [ surrogval{1} surrogval{2} surrogval{3} ] = anova2_cell_select( res, g.paired );
                    catch,
                        myfprintf(verb,'\nSuperfast array computation failed because of memory limitation, reverting to standard computing');
                        g.arraycomp = 'off';
                    end;
                end;
                if strcmpi(g.arraycomp, 'off')
                    [res precomp] = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
                    for index = 1:g.naccu
                        if mod(index, 10) == 0, myfprintf(verb,'%d ', index); end;
                        if mod(index, 100) == 0, myfprintf(verb,'\n'); end;

                        res = surrogdistrib( {}, 'precomp', precomp);
                        if myndims(data{1}) == 1
                         	 [ surrogval{1}(index)     surrogval{2}(index)     surrogval{3}(index)     ] = anova2_cell_select( res, g.paired );
                        else [ surrogval{1}(:,index)   surrogval{2}(:,index)   surrogval{3}(:,index)   ] = anova2_cell_select( res, g.paired );
                        end;
                    end;
                end;
            end;
        end;
        myfprintf(verb,'\n');
    else
        surrogval = g.surrog;
        ori_vals  = g.stats;
        df        = [];
    end;

    % compute p-values
    % ----------------
    if ~strcmpi(g.method, 'param')
        if iscell( surrogval )
            pvals{1} = stat_surrogate_pvals(surrogval{1}, ori_vals{1}, tail);
            pvals{2} = stat_surrogate_pvals(surrogval{2}, ori_vals{2}, tail);
            pvals{3} = stat_surrogate_pvals(surrogval{3}, ori_vals{3}, tail);
        else
            pvals = stat_surrogate_pvals(surrogval, ori_vals, tail);
        end;
        try, warning('on', 'MATLAB:divideByZero'); catch, end;
    end;

    [ ori_vals, pvals ] = reshape_results( nd, ori_vals, pvals);
    [ surrogval ]       = reshape_results( [nd g.naccu], surrogval);

    % confidence intervals
    % --------------------
    if ~isnan(g.alpha)
        outputstruct.ci = stat_surrogate_ci(surrogval, g.alpha, tail);
        if strcmpi(g.structoutput, 'off')
            disp('Warning: returning confidence interval requires an output structure');
        end;
        if iscell(pvals)
            for ind = 1:length(pvals)
                outputstruct.mask{ind} = pvals{ind} < g.alpha;
            end;
        else
            outputstruct.mask = pvals < g.alpha;
        end;
    end;

    % create a structure for outputing values
    % ---------------------------------------
    if strcmpi(g.structoutput, 'on')
        outputstruct.method = g.method;
        outputstruct.pval   = pvals;
        outputstruct.df     = df;
        outputstruct.surrog = surrogval;
        if length(data(:)) == 2
             outputstruct.t = ori_vals;
        else outputstruct.f = ori_vals;
        end;
        outputstruct.stat   = ori_vals;
        ori_vals = outputstruct;
    end;

% compute ANOVA 2-way
% -------------------
function [f1 f2 f3 df1 df2 df3] = anova2_cell_select( res, paired);
    if strcmpi(paired,'on')
        [f1 f2 f3 df1 df2 df3] = anova2rm_cell( res );
    else
        [f1 f2 f3 df1 df2 df3] = anova2_cell( res );
    end;

% compute ANOVA 1-way
% -------------------
function [f df] = anova1_cell_select( res, paired);
    if strcmpi(paired,'on')
        [f df] = anova1rm_cell( res );
    else
        [f df] = anova1_cell( res );
    end;

% compute t-test
% -------------------
function [t df] = ttest_cell_select( res, paired, homogenous);
    if strcmpi(paired,'on')
        [t df] = ttest_cell( res{1}, res{2});
    else
        [t df] = ttest2_cell( res{1}, res{2}, homogenous);
    end;

% function to compute the number of dimensions
% --------------------------------------------
function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end;
    end;

% function for verbose messages
% -----------------------------
function myfprintf(verb, varargin)
    if verb
        fprintf(varargin{:});
    end;

% function to replace tcdf
% ------------------------
function p = mytcdf(x,v)

if length(v) == 1,
    v = repmat(v, size(x));
end;

x2 = x.^2;
inds1 = (v < x2);
inds2 = (v >= x2);
if any(inds1(:)), p(inds1) = betainc(v(inds1) ./ (v(inds1) + x2(inds1)), v(inds1)/2, 0.5, 'lower') / 2; end;
if any(inds2(:)), p(inds2) = betainc(x2(inds2) ./ (v(inds2) + x2(inds2)), 0.5, v(inds2)/2, 'upper') / 2; end;
inds = (x > 0);
if any(inds)
    p(inds) = 1 - p(inds);
end;

inds = (v > 1e7);
if any(inds(:)), p(inds) = normcum(x(inds)); end;

p(x == 0) = 0.5;
if isempty(p)
    p = ones(size(x));
else
    p = reshape(p, size(x));
end;
function [p] = normcum(z)
p = 0.5 * erfc(-z ./ sqrt(2));

% reshape results
% ---------------
function varargout = reshape_results(nd, varargin)
    if length(varargin) > 1
        for index = 1:length(varargin)
            varargout{index} = reshape_results(nd, varargin{index});
        end;
    elseif iscell(varargin{1})
        for index = 1:length(varargin{1})
            varargout{1}{index} = reshape_results(nd, varargin{1}{index});
        end;
    else
        if ~isempty(varargin{1})
            if length(nd) == 1, nd = [ nd 1 ]; end;
            varargout{1} = reshape(varargin{1}, nd);
        else varargout{1} = [];
        end;
    end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [DCM] = spm_dcm_fmri_check(P, varargin)
% post-hoc diagnostics for DCM (bilinear or nonlinear) of fMRI data
% FORMAT [DCM] = spm_dcm_fmri_check(DCM)
%   DCM     - DCM structure or its filename
%
% FORMAT [GCM] = spm_dcm_fmri_check(GCM)
%   GCM     - Subjects x Models cell array of DCM structures or filenames
%
% FORMAT [DCM] = spm_dcm_fmri_check(DCM, nograph, GCM)
%   DCM     - DCM structure or its filename
%   nograph - (Optional) if true, disables graphical output
%   GCM     - (Optional) full GCM array from which the DCM in P was sourced
%             for use in graphics
%
% This routine provides some diagnostics to ensure model inversion has
% converged. It plots the predicted and observed responses over all regions
% and provides the coefficient of determination - or percent variance
% explained. This should normally be above 10%. An abnormally low
% coefficient of determination is highlighted in red. Quantitatively, one
% would normally expect to see one or more extrinsic (between source)
% connections with the strength of 1/8 Hz or greater. If all the extrinsic
% posterior expectations are below this value, then this suggests a failure
% of convergence or that the data are very noisy (possibly due to using
% very small regions of interest to summarise regional responses). Finally,
% the posterior correlations among all parameters are shown in terms of a
% correlation matrix. The number of effective parameters estimated is
% reported in terms of the (KL) divergence between the posterior and
% prior densities over parameters. This is divided by the log of the
% number of observations, by appealing to the Bayesian information
% criterion. The divergence corresponds to complexity or Bayesian
% surprise. Normally, one would expect the posterior and prior to diverge
% in a non-trivial fashion.
%
% Posterior densities are shown as bars with 90% confidence intervals in
% pink. An informed model inversion would normally provide posterior
% densities with confidence intervals that are, for some connections,
% displaced from prior expectations (at or around zero).
%
% The following diagnostics are stored in the returned DCM:
%
% DCM.diagnostics(1) - Percent variance explained
% DCM.diagnostics(2) - Largest absolute parameter estimate
% DCM.diagnostics(3) - Effective number of parameters estimated
%
% This routine is compatible with DCM8, DCM10 and DCM12 files.
%__________________________________________________________________________
% Copyright (C) 2012-2013 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_fmri_check.m 7639 2019-07-16 14:48:11Z peter $

%-Prepare inputs
%--------------------------------------------------------------------------
if isempty(varargin)
    nograph = false;
else
    nograph = varargin{1};
end

if length(varargin) < 2
    GCM = [];
else
    GCM = varargin{2};
end

%-Load DCM structure
%--------------------------------------------------------------------------
if ~nargin
    [P, sts] = spm_select(1,'^(D|G)CM.*\.mat$','select DCM or GCM mat');
    if ~sts, DCM = []; return; end
end

if isstruct(P)
    DCM = P;
elseif ischar(P)
    DCM = load(P);
    if isfield(DCM,'DCM')
        DCM = DCM.DCM;
    elseif isfield(DCM,'GCM')
        DCM = DCM.GCM;
    else
        error('Unknown DCM format');
    end
elseif iscell(P)
    DCM = P;
end

%-Handle multiple DCMs (group DCM structure)
%--------------------------------------------------------------------------
if iscell(DCM)

    % Prepare figure
    if ~nograph
        datacursormode off;
        f = spm_figure('GetWin','DCM diagnostics'); spm_clf;
        colormap(jet);
        add_title('Loading...',...
                  'Position',[0,0.5,1,0.2], ...
                  'ForegroundColor',[0.7 0.7 0.7]);
        drawnow;
    end

    % Load if filenames given
    if ischar(DCM{1})
        DCM = spm_dcm_load(DCM);
    end

    % Call spm_dcm_fmri_check recursively to assemble diagnostics
    [stats,DCM] = get_diagnostics(DCM);

    if nograph, return; end

    % Extract explained variance
    data = cellfun(@(x)x(1),stats);

    % Add title
    spm_clf; subplot(10,2,1:2); axis off;
    p = get(gca,'Position');
    add_title('DCM for fMRI Diagnostics','Position',[0,p(2)+p(4)/2,1,0.06]);

    % Plot grid of models
    subplot(10,2,3:2:13);
    im = plot_model_space(data, f);

    % Create diagnostic panel
    subplot(10,2,4:2:14); axis off;
    h = create_diagnostic_panel();

    % Store handles and DCM in image
    h.DCM = DCM;
    set(im,'UserData',h);

    % Enable clicks
    datacursormode on;

    % Plot performance over empirical bayes iterations
    if isfield(DCM{1},'FEB')
        subplot(10,3,[25 28]), bar(DCM{1}.FEB(:) - DCM{1}.FEB(1))
        xlabel('iteration','FontSize',12);
        title('Free energy','FontSize',16)
        axis square

        subplot(10,3,[26 29]), bar(cellfun(@mean,DCM{1}.EEB(:)),'b')
        xlabel('iteration','FontSize',12)
        title('Log precision','FontSize',16)
        axis square

        subplot(10,3,[27 30]), bar(DCM{1}.HEB(:) - DCM{1}.HEB(1),'c')
        xlabel('iteration','FontSize',12)
        title('Posterior uncertainty','FontSize',16)
        axis square
    end

    return;
end

% Assemble diagnostics
%==========================================================================

% coefficient of determination (percent variance explained)
%--------------------------------------------------------------------------

% Check if spectral DCM (DCM for cross spectra)
%--------------------------------------------------------------------------
try
    analysis = DCM.options.analysis;
catch
    analysis = '';
end
if strcmp(analysis,'CSD') && isfield(DCM,'Hc')
    PSS   = sum(sum(sum(abs(DCM.Hc).^2)));
    RSS   = sum(sum(sum(abs(DCM.Rc).^2)));
elseif isfield(DCM,'y')
    PSS   = sum(sum(DCM.y.^2));
    RSS   = sum(sum(DCM.R.^2));
else
    PSS   = NaN;
    RSS   = NaN;
end

D(1)  = 100*PSS/(PSS + RSS);

% largest absolute posterior expectation (extrinsic connections)
%--------------------------------------------------------------------------
try
    A = DCM.Ep.A;
catch
    A = DCM.A;
end

if isfield(DCM.options,'two_state') && DCM.options.two_state
    A = exp(A);
end

D(2)  = max(max(abs(A - diag(diag(A)))));

% complexity and effective number of parameters estimated
%--------------------------------------------------------------------------
qE    = spm_vec(DCM.Ep);
pE    = spm_vec(DCM.M.pE);
qC    = DCM.Cp;
pC    = full(DCM.M.pC);
k     = rank(full(pC));
pC    = pinv(pC);

D(3)  = trace(pC*qC) + (pE - qE)'*pC*(pE - qE) - spm_logdet(qC*pC) - k;
D(3)  = D(3)/log(DCM.v);

D     = full(D);
DCM.diagnostics = D;

if nograph
    return;
end

% Plot summary of inversion
%==========================================================================
spm_figure('GetWin','DCM diagnostics'); clf

% back to GCM link if needed
%--------------------------------------------------------------------------
if ~isempty(GCM)
   uicontrol('Style','PushButton','String','Return to models',...
        'Enable','Inactive','ButtonDownFcn',@backbutton_clicked, ...
        'Units','Normalized','Position',[0.04 0.96 0.2 0.03],...
        'HorizontalAlignment','Left','UserData',GCM);
end

% plot predicted and observed regional responses
%--------------------------------------------------------------------------
subplot(2,1,1);

% Check if spectral DCM (DCM for cross spectra)
%--------------------------------------------------------------------------
if strcmp(analysis,'CSD')

    %-CSD data
    %----------------------------------------------------------------------
    Hz   = DCM.Hz;                 % frequencies
    name = {DCM.xY.name};          % names
    ns   = size(DCM.a,1);          % number of regions
    ns   = min(ns,8);              % bounded number of regions

    c = lines(ns);
    h = [];
    for i = 1:ns
        Hc(:,i) = abs(DCM.Hc(:,i,i));
        Yc(:,i) = abs(DCM.Hc(:,i,i) + DCM.Rc(:,i,i));

        h(i)=plot(Hz,Hc(:,i),'Color',c(i,:)); hold on
        plot(Hz,Yc(:,i),':','Color',c(i,:))
    end
    hold off

    str = sprintf('variance explained %0.0f%%', D(1));
    str = {'Responses and Predictions',str};

    if D(1) > 10
        title(str,'FontSize',16);
    else
        title(str,'FontSize',16,'Color','r');
    end

    xlabel('frequency (Hz)')
    ylabel('abs(CSD)')
    axis square, spm_axis tight
    try
        legend(h,name)
    end
else
    t   = (1:DCM.v)*DCM.Y.dt;
    plot(t,DCM.y,t,DCM.y + DCM.R,':');
    str = sprintf('variance explained %0.0f%%', D(1));
    str = {'Responses and Predictions',str};

    if D(1) > 10
        title(str,'FontSize',16);
    else
        title(str,'FontSize',16,'Color','r');
    end

    xlabel('time {seconds}');
end

% posterior densities over A parameters
%--------------------------------------------------------------------------
try
    i = spm_fieldindices(DCM.Ep,'A');
catch
    i = 1 + (1:DCM.n^2);
end
qE  = spm_vec(DCM.Ep);
qC  = DCM.Cp;

if DCM.options.two_state
    qE = exp(qE);
end

subplot(2,2,3)
spm_plot_ci(qE(i),qC(i,i)), hold on
str = sprintf('largest connection strength %0.2f', D(2));
str = {'Intrinsic and Extrinsic connections',str};
if D(2) > 1/8
    title(str,'FontSize',16);
else
    title(str,'FontSize',16,'Color','r');
end
xlabel('parameters');
axis square

% posterior correlations among all parameters
%--------------------------------------------------------------------------
subplot(2,2,4)
imagesc(spm_cov2corr(DCM.Cp))
title('Posterior Correlations','FontSize',16)
str = sprintf('estimable parameters %0.0f', D(3));
str = {'Posterior Correlations',str};
if D(3) > 1
    title(str,'FontSize',16);
else
    title(str,'FontSize',16,'Color','r');
end
axis square

% =========================================================================
function [stats,GCM] = get_diagnostics(GCM)
% Builds a [s x m] cell array of diagnostic stats for an array of DCMs

stats    = cell(size(GCM));
stats(:) = {zeros(1,3)};
[Ns, Nm] = size(GCM);
counter = 1;
for s = 1:Ns
    for m = 1:Nm
        if isfield(GCM{s,m},'Ep')
            GCM{s,m}.v  = GCM{s,1}.v;
            GCM{s,m}    = spm_dcm_fmri_check(GCM{s,m}, true);
            stats{s,m}  = GCM{s,m}.diagnostics;
        else
            stats{s,m}  = zeros(3,1);
        end
        counter = counter + 1;
    end
end

% =========================================================================
function add_title(str,varargin)
% Draws a title
uicontrol('Style','Text','String',str,'Fontsize',20,...
          'Units','Normalized','HorizontalAlignment','Center',...
          'BackgroundColor','w',varargin{:});

% =========================================================================
function im = plot_model_space(data, f)
% Plots the model space. Returns a graphics handle to the image

% Plot explained variances
im = imagesc(data, [0 100]);

[Ns, Nm] = size(data);
set(gca,'xtick', linspace(0.5,Nm+0.5,Nm+1), ...
        'ytick', linspace(0.5,Ns+.5,Ns+1));
set(gca,'xgrid', 'on', 'ygrid', 'on', ...
        'gridlinestyle', '-', 'xcolor', 'w', 'ycolor', 'w');
ylabel('Subjects','FontSize',14,'Color','k');
xlabel('Models','FontSize',14,'Color','k');
title('Variance explained (%)','FontSize',16);
axis equal

% Create small colorbar
cbar = colorbar;
p    = get(cbar,'Position');
t    = p(2) + p(4);
p(2) = t - 0.2;
p(4) = 0.2;
set(cbar,'Position',p);

% Handle mouseclicks
dcm_obj = datacursormode(f);
set(dcm_obj,'UpdateFcn',@gcm_clicked);

% =========================================================================
function h = create_diagnostic_panel()
% Create a panel to show details / controls for the selected model

% Get axes top
p = get(gca,'Position');
t = p(2)+p(4);

% Add panel
he = 0.15; % Height
h = struct();
h.panel = uipanel('Units','normalized',...
'BackgroundColor',[1 1 1],'Title','Diagnostics',...
'Position',[0.5 t-he 0.4 he],'FontSize',12);

% Add explained variance label
h.expvar = uicontrol('Style','text','Units','normalized',...
          'String','Please select a model, left.',...
          'BackgroundColor',[1 1 1],'Position',[0,0.6,1,0.2],...
          'Parent',h.panel,'HorizontalAlignment','center');

% Add button
h.infobutton = uicontrol('Style','pushbutton','Units','normalized',...
          'String','Diagnostics',...
          'Position',[0.1,0.2,0.8,0.25],...
          'Parent',h.panel,'Callback',@info_clicked,'Enable','Off',...
          'Visible','Off');

% =========================================================================
function txt = gcm_clicked(varargin)
% Handle mouseclick on DCM array

% Unpack handles and explained variance matrix
im     = get(varargin{2},'Target');
h      = get(im,'UserData');
expvar = get(im,'CData');

% Get DCM for this click position
xy  = get(varargin{2},'Position');
m   = xy(1);
s   = xy(2);
GCM = h.DCM;
DCM = GCM{s,m};

% Update DCM info panel text
if isnan(expvar(s,m))
    str = 'Variance explained: unavailable';
    set(h.infobutton,'Enable','off');
else
    str = sprintf('Variance explained: %2.2f%%', expvar(s,m));
    set(h.infobutton,'Enable','on');
end
set(h.expvar, 'String', str);
set(h.panel,  'Title', sprintf('Subject %d Model %d',s,m));

% Update DCM info panel button
button_data = struct();
button_data.DCM = DCM;
button_data.GCM = GCM;
set(h.infobutton,'UserData', button_data);
set(h.infobutton,'Visible', 'on');

% Return tooltip
txt = sprintf('Subject %d model %d',s,m);

% =========================================================================
function info_clicked(varargin)
% Callback for the Diagnostics button being clicked

button_data = get(varargin{1},'UserData');
DCM = button_data.DCM;
GCM = button_data.GCM;

if ~isempty(button_data.DCM)
    dcm_obj = datacursormode(gcf);
    set(dcm_obj,'UpdateFcn',[]);
    datacursormode off;
    spm_dcm_fmri_check(DCM, false, GCM);
end

% =========================================================================
function backbutton_clicked(varargin)
% Callback for the back button, to return to the model space
GCM = get(varargin{1},'UserData');
clf; drawnow();
spm_dcm_fmri_check(GCM);

author: spm


text:
function jac = autoJac(fun,x0,varargin)
%AUTOJAC  Returns a Automatically Differentiated (AD) Jacobian
%
%   jac = autoJac(fun,x0) uses the user supplied Matlab function to
%   generate the Jacobian using automatic differentiation. The user
%   function must be a Matlab function and must not call external code or
%   class / toolbox functions. If your function breaks any of these
%   conditions consider using mklJac instead.
%
%   The underlying AD algorithm is adiff by William McIlhagga and its
%   documentation pdf is provided in the Differentiation folder. See the
%   BSD license below the code.

%   Copyright (C) 2011 Jonathan Currie (IPL)

if(~isa(fun,'function_handle'))
    error('Fun should be a function handle!');
end

[~,jac] = adiffget(fun(adiff(x0),varargin{:}));

end

% Copyright (c) 2010, William McIlhagga
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * 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 OWNER 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.

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function chebyshev_series_test03 ( )

%*****************************************************************************80
%
%% CHEBYSHEV_SERIES_TEST03 considers an odd Chebyshev series for SINH(X).
%
%  Discussion:
%
%    TABLE5ODD contains the odd Chebyshev series coefficients for
%    sinh(x) over -1 <= x <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2014
%
%  Author:
%
%    Manfred Zimmer
%
%  Reference:
%
%    Charles Clenshaw,
%    Mathematical Tables, Volume 5,
%    Chebyshev series for mathematical functions,
%    London, 1962.
%
  nc = 9;

  table5odd = [ ...
    1.13031820798497005442, ...
    0.04433684984866380495, ...
    0.00054292631191394375, ...
    0.00000319843646240199, ...
    0.00000001103677172552, ...
    0.00000000002497956617, ...
    0.00000000000003991263, ...
    0.00000000000000004741, ...
    0.00000000000000000004 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBYSHEV_SERIES_TEST03:\n' );
  fprintf ( 1, '  ODDCHEBSER2 computes an odd Chebyshev series approximation.\n' );
  fprintf ( 1, '  and its first two derivatives.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Errors of an odd Chebyshev series y(x) approximating sinh(x):\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    x        err(y)       err(y'')       err(y")\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 10
    x = i / 10.0;
    [ s, s1, s2 ] = oddchebser2 ( x, table5odd, nc );
    y = sinh ( x );
    y1 = cosh ( x );
    s = s - y;
    s1 = s1 - y1;
    s2 = s2 - y;
    fprintf ( 1, '%7.4f  %14.6g  %14.6g  %14.6g\n', x, s, s1, s2 );
  end

  return
end


text:
function varargout = mtimes(varargin)
%*	   Pointwise multiplication for CHEBFUN2 objects.
%
%   c*F or F*c multiplies a CHEBFUN2 F by a scalar c.
%
%   F*G computes the integral of F(s,y)G(x,s) over s, the continuous
%   analogue of matrix-matrix multiplication.
%
% See also TIMES.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = mtimes@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef AssessmentPerformanceMetricResult < handle
   properties (Access = public)
       algorithmName;
       segmentationPerformanceMetrics;
       featureExtractionPerformanceMetrics;
       classificationPerformanceMetrics;
       predictionPerformanceMetrics;
   end

end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK