text:
function [D PD] = allfitdist(data,sortby,varargin)
%ALLFITDIST Fit all valid parametric probability distributions to data.
%   [D PD] = ALLFITDIST(DATA) fits all valid parametric probability
%   distributions to the data in vector DATA, and returns a struct D of
%   fitted distributions and parameters and a struct of objects PD
%   representing the fitted distributions. PD is an object in a class
%   derived from the ProbDist class.
%
%   [...] = ALLFITDIST(DATA,SORTBY) returns the struct of valid distributions
%   sorted by the parameter SORTBY
%        NLogL - Negative of the log likelihood
%        BIC - Bayesian information criterion (default)
%        AIC - Akaike information criterion
%        AICc - AIC with a correction for finite sample sizes
%
%   [...] = ALLFITDIST(...,'DISCRETE') specifies it is a discrete
%   distribution and does not attempt to fit a continuous distribution
%   to the data
%
%   [...] = ALLFITDIST(...,'PDF') or (...,'CDF') plots either the PDF or CDF
%   of a subset of the fitted distribution. The distributions are plotted in
%   order of fit, according to SORTBY.
%
%   List of distributions it will try to fit
%     Continuous (default)
%       Beta
%       Birnbaum-Saunders
%       Exponential
%       Extreme value
%       Gamma
%       Generalized extreme value
%       Generalized Pareto
%       Inverse Gaussian
%       Logistic
%       Log-logistic
%       Lognormal
%       Nakagami
%       Normal
%       Rayleigh
%       Rician
%       t location-scale
%       Weibull
%
%     Discrete ('DISCRETE')
%       Binomial
%       Negative binomial
%       Poisson
%
%   Optional inputs:
%   [...] = ALLFITDIST(...,'n',N,...)
%   For the 'binomial' distribution only:
%      'n'            A positive integer specifying the N parameter (number
%                     of trials).  Not allowed for other distributions. If
%                     'n' is not given it is estimate by Method of Moments.
%                     If the estimated 'n' is negative then the maximum
%                     value of data will be used as the estimated value.
%   [...] = ALLFITDIST(...,'theta',THETA,...)
%   For the 'generalized pareto' distribution only:
%      'theta'        The value of the THETA (threshold) parameter for
%                     the generalized Pareto distribution. Not allowed for
%                     other distributions. If 'theta' is not given it is
%                     estimated by the minimum value of the data.
%
%   Note: ALLFITDIST does not handle nonparametric kernel-smoothing,
%   use FITDIST directly instead.
%
%
%   EXAMPLE 1
%     Given random data from an unknown continuous distribution, find the
%     best distribution which fits that data, and plot the PDFs to compare
%     graphically.
%        data = normrnd(5,3,1e4,1);         %Assumed from unknown distribution
%        [D PD] = allfitdist(data,'PDF');   %Compute and plot results
%        D(1)                               %Show output from best fit
%
%   EXAMPLE 2
%     Given random data from a discrete unknown distribution, with frequency
%     data, find the best discrete distribution which would fit that data,
%     sorted by 'NLogL', and plot the PDFs to compare graphically.
%        data = nbinrnd(20,.3,1e4,1);
%        values=unique(data); freq=histc(data,values);
%        [D PD] = allfitdist(values,'NLogL','frequency',freq,'PDF','DISCRETE');
%        PD{1}
%
%  EXAMPLE 3
%     Although the Geometric Distribution is not listed, it is a special
%     case of fitting the more general Negative Binomial Distribution. The
%     parameter 'r' should be close to 1. Show by example.
%        data=geornd(.7,1e4,1); %Random from Geometric
%        [D PD]= allfitdist(data,'PDF','DISCRETE');
%        PD{1}
%
%  EXAMPLE 4
%     Compare the resulting distributions under two different assumptions
%     of discrete data. The first, that it is known to be derived from a
%     Binomial Distribution with known 'n'. The second, that it may be
%     Binomial but 'n' is unknown and should be estimated. Note the second
%     scenario may not yield a Binomial Distribution as the best fit, if
%     'n' is estimated incorrectly. (Best to run example a couple times
%     to see effect)
%        data = binornd(10,.3,1e2,1);
%        [D1 PD1] = allfitdist(data,'n',10,'DISCRETE','PDF'); %Force binomial
%        [D2 PD2] = allfitdist(data,'DISCRETE','PDF');       %May be binomial
%        PD1{1}, PD2{1}                             %Compare distributions
%

%    Mike Sheppard
%    Last Modified: 17-Feb-2012

%% Check Inputs
if nargin == 0
    data = 10.^((normrnd(2,10,1e4,1))/10);
    sortby='BIC';
    varargin={'CDF'};
end
if nargin==1
    sortby='BIC';
end
sortbyname={'NLogL','BIC','AIC','AICc'};
if ~any(ismember(lower(sortby),lower(sortbyname)))
    oldvar=sortby; %May be 'PDF' or 'CDF' or other commands
    if isempty(varargin)
        varargin={oldvar};
    else
        varargin=[oldvar varargin];
    end
    sortby='BIC';
end
if nargin < 2, sortby='BIC'; end
distname={'beta', 'birnbaumsaunders', 'exponential', ...
    'extreme value', 'gamma', 'generalized extreme value', ...
    'generalized pareto', 'inversegaussian', 'logistic', 'loglogistic', ...
    'lognormal', 'nakagami', 'normal', ...
    'rayleigh', 'rician', 'tlocationscale', 'weibull'};
if ~any(strcmpi(sortby,sortbyname))
    error('allfitdist:SortBy','Sorting must be either NLogL, BIC, AIC, or AICc');
end
%Input may be mixed of numeric and strings, find only strings
vin=varargin;
strs=find(cellfun(@(vs)ischar(vs),vin));
vin(strs)=lower(vin(strs));
%Next check to see if 'PDF' or 'CDF' is listed
numplots=sum(ismember(vin(strs),{'pdf' 'cdf'}));
if numplots>=2
    error('ALLFITDIST:PlotType','Either PDF or CDF must be given');
end
if numplots==1
    plotind=true; %plot indicator
    indxpdf=ismember(vin(strs),'pdf');
    plotpdf=any(indxpdf);
    indxcdf=ismember(vin(strs),'cdf');
    vin(strs(indxpdf|indxcdf))=[]; %Delete 'PDF' and 'CDF' in vin
else
    plotind=false;
end
%Check to see if discrete
strs=find(cellfun(@(vs)ischar(vs),vin));
indxdis=ismember(vin(strs),'discrete');
discind=false;
if any(indxdis)
    discind=true;
    distname={'binomial', 'negative binomial', 'poisson'};
    vin(strs(indxdis))=[]; %Delete 'DISCRETE' in vin
end
strs=find(cellfun(@(vs)ischar(vs),vin));
n=numel(data); %Number of data points
data = data(:);
D=[];
%Check for NaN's to delete
deldatanan=isnan(data);
%Check to see if frequency is given
indxf=ismember(vin(strs),'frequency');
if any(indxf)
    freq=vin{1+strs((indxf))}; freq=freq(:);
    if numel(freq)~=numel(data)
        error('ALLFITDIST:PlotType','Matrix dimensions must agree');
    end
    delfnan=isnan(freq);
    data(deldatanan|delfnan)=[]; freq(deldatanan|delfnan)=[];
    %Save back into vin
    vin{1+strs((indxf))}=freq;
else
    data(deldatanan)=[];
end

%% Run through all distributions in FITDIST function
warning('off','all'); %Turn off all future warnings
for indx=1:length(distname)
    try
        dname=distname{indx};
        switch dname
            case 'binomial'
                PD=fitbinocase(data,vin,strs); %Special case
            case 'generalized pareto'
                PD=fitgpcase(data,vin,strs); %Special case
            otherwise
                %Built-in distribution using FITDIST
                PD = fitdist(data,dname,vin{:});
        end

        NLL=PD.NLogL; % -Log(L)
        %If NLL is non-finite number, produce error to ignore distribution
        if ~isfinite(NLL)
            error('non-finite NLL');
        end
        num=length(D)+1;
        PDs(num) = {PD}; %#ok<*AGROW>
        k=numel(PD.Params); %Number of parameters
        D(num).DistName=PD.DistName;
        D(num).NLogL=NLL;
        D(num).BIC=-2*(-NLL)+k*log(n);
        D(num).AIC=-2*(-NLL)+2*k;
        D(num).AICc=(D(num).AIC)+((2*k*(k+1))/(n-k-1));
        D(num).ParamNames=PD.ParamNames;
        D(num).ParamDescription=PD.ParamDescription;
        D(num).Params=PD.Params;
        D(num).Paramci=PD.paramci;
        D(num).ParamCov=PD.ParamCov;
        D(num).Support=PD.Support;
    catch err %#ok<NASGU>
        %Ignore distribution
    end
end
warning('on','all'); %Turn back on warnings
if numel(D)==0
    error('ALLFITDIST:NoDist','No distributions were found');
end

%% Sort distributions
indx1=1:length(D); %Identity Map
sortbyindx=find(strcmpi(sortby,sortbyname));
switch sortbyindx
    case 1
        [~,indx1]=sort([D.NLogL]);
    case 2
        [~,indx1]=sort([D.BIC]);
    case 3
        [~,indx1]=sort([D.AIC]);
    case 4
        [~,indx1]=sort([D.AICc]);
end
%Sort
D=D(indx1); PD = PDs(indx1);

%% Plot if requested
if plotind;
    plotfigs(data,D,PD,vin,strs,plotpdf,discind)
end

end

function PD=fitbinocase(data,vin,strs)
%% Special Case for Binomial
% 'n' is estimated if not given
vinbino=vin;
%Check to see if 'n' is given
indxn=any(ismember(vin(strs),'n'));
%Check to see if 'frequency' is given
indxfreq=ismember(vin(strs),'frequency');
if ~indxn
    %Use Method of Moment estimator
    %E[x]=np, V[x]=np(1-p) -> nhat=E/(1-(V/E));
    if isempty(indxfreq)||~any(indxfreq)
        %Raw data
        mnx=mean(data);
        nhat=round(mnx/(1-(var(data)/mnx)));
    else
        %Frequency data
        freq=vin{1+strs(indxfreq)};
        m1=dot(data,freq)/sum(freq);
        m2=dot(data.^2,freq)/sum(freq);
        mnx=m1; vx=m2-(m1^2);
        nhat=round(mnx/(1-(vx/mnx)));
    end
    %If nhat is negative, use maximum value of data
    if nhat<=0, nhat=max(data(:)); end
    vinbino{end+1}='n'; vinbino{end+1}=nhat;
end
PD = fitdist(data,'binomial',vinbino{:});
end

function PD=fitgpcase(data,vin,strs)
%% Special Case for Generalized Pareto
% 'theta' is estimated if not given
vingp=vin;
%Check to see if 'theta' is given
indxtheta=any(ismember(vin(strs),'theta'));
if ~indxtheta
    %Use minimum value for theta, minus small part
    thetahat=min(data(:))-10*eps;
    vingp{end+1}='theta'; vingp{end+1}=thetahat;
end
PD = fitdist(data,'generalized pareto',vingp{:});
end

function plotfigs(data,D,PD,vin,strs,plotpdf,discind)
%Plot functionality for continuous case due to Jonathan Sullivan
%Modified by author for discrete case

%Maximum number of distributions to include
%max_num_dist=Inf;  %All valid distributions
max_num_dist=4;

%Check to see if frequency is given
indxf=ismember(vin(strs),'frequency');
if any(indxf)
    freq=vin{1+strs((indxf))};
end

figure

%% Probability Density / Mass Plot
if plotpdf
    if ~discind
        %Continuous Data
        nbins = max(min(length(data)./10,100),50);
        xi = linspace(min(data),max(data),nbins);
        dx = mean(diff(xi));
        xi2 = linspace(min(data),max(data),nbins*10)';
        fi = histc(data,xi-dx);
        fi = fi./sum(fi)./dx;
        inds = 1:min([max_num_dist,numel(PD)]);
        ys = cellfun(@(PD) pdf(PD,xi2),PD(inds),'UniformOutput',0);
        ys = cat(2,ys{:});
        bar(xi,fi,'FaceColor',[160 188 254]/255,'EdgeColor','k');
        hold on;
        plot(xi2,ys,'LineWidth',1.5)
        legend(['empirical',{D(inds).DistName}],'Location','NE')
        xlabel('Value');
        ylabel('Probability Density');
        title('Probability Density Function');
        grid on
    else
        %Discrete Data
        xi2=min(data):max(data);
        %xi2=unique(x)'; %If only want observed x-values to be shown
        indxf=ismember(vin(strs),'frequency');
        if any(indxf)
            fi=zeros(size(xi2));
            fi((ismember(xi2,data)))=freq; fi=fi'./sum(fi);
        else
            fi=histc(data,xi2); fi=fi./sum(fi);
        end
        inds = 1:min([max_num_dist,numel(PD)]);
        ys = cellfun(@(PD) pdf(PD,xi2),PD(inds),'UniformOutput',0);
        ys=cat(1,ys{:})';
        bar(xi2,[fi ys]);
        legend(['empirical',{D(inds).DistName}],'Location','NE')
        xlabel('Value');
        ylabel('Probability Mass');
        title('Probability Mass Function');
        grid on
    end
else

%Cumulative Distribution
    if ~discind
        %Continuous Data
        [fi xi] = ecdf(data);
        inds = 1:min([max_num_dist,numel(PD)]);
        ys = cellfun(@(PD) cdf(PD,xi),PD(inds),'UniformOutput',0);
        ys = cat(2,ys{:});
        if max(xi)/min(xi) > 1e4; lgx = true; else lgx = false; end
        subplot(2,1,1)
        if lgx
            semilogx(xi,fi,'k',xi,ys)
        else
            plot(xi,fi,'k',xi,ys)
        end
        legend(['empirical',{D(inds).DistName}],'Location','NE')
        xlabel('Value');
        ylabel('Cumulative Probability');
        title('Cumulative Distribution Function');
        grid on
        subplot(2,1,2)
        y = 1.1*bsxfun(@minus,ys,fi);
        if lgx
            semilogx(xi,bsxfun(@minus,ys,fi))
        else
            plot(xi,bsxfun(@minus,ys,fi))
        end
        ybnds = max(abs(y(:)));
        ax = axis;
        axis([ax(1:2) -ybnds ybnds]);
        legend({D(inds).DistName},'Location','NE')
        xlabel('Value');
        ylabel('Error');
        title('CDF Error');
        grid on
    else
        %Discrete Data
        indxf=ismember(vin(strs),'frequency');
        if any(indxf)
            [fi xi] = ecdf(data,'frequency',freq);
        else
            [fi xi] = ecdf(data);
        end
        %Check unique xi, combine fi
        [xi,ign,indx]=unique(xi); %#ok<ASGLU>
        fi=accumarray(indx,fi);
        inds = 1:min([max_num_dist,numel(PD)]);
        ys = cellfun(@(PD) cdf(PD,xi),PD(inds),'UniformOutput',0);
        ys=cat(2,ys{:});
        subplot(2,1,1)
        stairs(xi,[fi ys]);
        legend(['empirical',{D(inds).DistName}],'Location','NE')
        xlabel('Value');
        ylabel('Cumulative Probability');
        title('Cumulative Distribution Function');
        grid on
        subplot(2,1,2)
        y = 1.1*bsxfun(@minus,ys,fi);
        stairs(xi,bsxfun(@minus,ys,fi))
        ybnds = max(abs(y(:)));
        ax = axis;
        axis([ax(1:2) -ybnds ybnds]);
        legend({D(inds).DistName},'Location','NE')
        xlabel('Value');
        ylabel('Error');
        title('CDF Error');
        grid on
    end
end

end

author: Sable

repo: mcbench-benchmarks


text:
function dist = lines_imp_dist_2d ( a1, b1, c1, a2, b2, c2 )

%*****************************************************************************80
%
%% LINES_IMP_DIST_2D determines the distance between two implicit lines in 2D.
%
%  Discussion:
%
%    If the lines intersect, then their distance is zero.
%    If the two lines are parallel, then they have a nonzero distance.
%
%    The implicit form of a line in 2D is:
%
%      A * X + B * Y + C = 0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A1, B1, C1, define the first line.
%    At least one of A1 and B1 must be nonzero.
%
%    Input, real A2, B2, C2, define the second line.
%    At least one of A2 and B2 must be nonzero.
%
%    Output, real DIST, the distance between the two lines.
%
  dim_num = 2;
%
%  Refuse to handle degenerate lines.
%
  if ( a1 == 0.0 & b1 == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LINES_IMP_DIST_2D - Fatal error!\n' );
    fprintf ( 1, '  Line 1 is degenerate.\n' );
    error ( 'LINES_IMP_DIST_2D - Fatal error!' );
  elseif ( a2 == 0.0 & b2 == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LINES_IMP_DIST_2D - Fatal error!\n' );
    fprintf ( 1, '  Line 2 is degenerate.\n' );
    error ( 'LINES_IMP_DIST_2D - Fatal error!' );
  end
%
%  Determine if the lines intersect.
%
  if ( a1 * b2 ~= a2 * b1 )
    dist = 0.0;
    return
  end
%
%  Determine the distance between the parallel lines.
%
  dist = abs ( c2 / sqrt ( a2 * a2 + b2 * b2 ) - c1 / sqrt ( a1 * a1 + b1 * b1 ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pass = test_emptyObjects()
% Does each chebfun3 command deal with empty chebfun3 objects,
% appropriately?

% Check things work for empty chebfun3 objects.
f = chebfun3();
try
    f + f;
    2*f;
    f*2;
    f.^2;
    2.^f;
    f.^f;
    sqrt(f);
    sum(f);
    sum2(f);
    sum3(f);
    integral(f);
    diff(f);
    sin(f);
    cos(f);
    sinh(f);
    f.^f + f;
    tucker(f)
    mean(f);
    max3(f);
    norm(f);
    minandmax3(f);
    permute(f);
    pass = 1;
catch
    pass = 0;
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: SSD versus rotations, linearInter, HNSP level=4
%
%==============================================================================

clear, close all, help(mfilename);

setup2DHNSPData;
imgModel('set','imgModel','linearInter');
level = 4; omega = ML{level}.omega; m = ML{level}.m;
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);
xc    = getCellCenteredGrid(omega,m);
Rc    = imgModel(R,omega,xc);

center = (omega(2:2:end)-omega(1:2:end))'/2;
trafo('set','trafo','rotation2D','c',center);

wc = pi/2*linspace(-1,1,101);  dc = zeros(size(wc));
figure(1); clf;
for j=1:length(wc),
  yc = trafo(wc(j),xc);
  Tc = imgModel(T,omega,yc);
  dc(j) = SSD(Tc,Rc,omega,m);
  viewImage(128+(Tc-Rc)/2,omega,m); drawnow;
  FAIRpause(1/6)
end;
figure; clf; p1 = plot(wc,dc);
%==============================================================================


text:
% Peak of Wigner-Ville Frequency Estimation
%
% Estimates the instantaneous frequency of the input signal by
% extracting the peaks of the Wigner-Ville distribution.
%
%
% Usage:
%
%     ife = wvpe( signal, lag_window_length, time_res [, fft_length] )
%
% Parameters:
%
%    signal
%
%	  Input one dimensional signal to be analysed. An analytic signal
%	  is required for this function, however, if signal is real, a
%	  default analytic transformer routine will be called from this
%	  function before computing tfd.
%
%    lag_window_length
%
%	  This is the data window length and controls the size of
%	  the kernel used for analysis (lag_window_length must be
%	  odd). The kernel used will be defined from -(lag_window_length+1)/2
%	  to +(lag_window_length+1)/2 in both time and lag dimensions.
%
%    time_res
%
%	  The number of time samples to skip between successive slices.
%	  Default vaule is 1.
%
%    fft_length
%
%	  Zero-padding at the FFT stage of the analysis may be specified by
%	  giving an fft_length larger than lag_window_length. If fft_length
%	  is not specified, or is smaller than the lag_window_length, then the
%	  next highest power of two above lag_window_length is used. If
%	  fft_length is not a power of two, the next highest power of two is
%	  used.
%
%     tfd
%
%	  The computed time-frequency distribution. size(tfd) will
%	  return [a, b], where a is the next largest power of two above
%	  fft_length, and b is floor(length(signal)/time_res) - 1.
%
%
%
%  See Also: wvd
% TFSAP 7.0
% Copyright Prof. B. Boashash
% Qatar University, Doha
% email: tfsap.research@gmail.com

author: HeLiangHIT


text:
function [seqName, seqFolder, imgFolder, frameRate, F, imWidth, imHeight, imgExt] ...
    = getSeqInfoFromFile(seq, dataDir)
% construct variables with relevant information about the sequence 'seq'
%

    seqName=char(seq);
    seqFolder= [dataDir,seqName,filesep];

    seqInfoFile = [dataDir,seqName,filesep,'seqinfo.ini'];
    ini = IniConfig();
    ini.ReadFile(seqInfoFile);

    imgFolder = ini.GetValues('Sequence','imDir');
    frameRate = ini.GetValues('Sequence','frameRate');
    F=ini.GetValues('Sequence','seqLength');
    imWidth=ini.GetValues('Sequence','imWidth');
    imHeight=ini.GetValues('Sequence','imHeight');
    imgExt = ini.GetValues('Sequence','imExt');

end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE: This script and pricing function are in progress, and have not been well tested
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Calculating forward value of the call with Antonov's mapping strategy
T      = 1;  %Time to maturity
r      = 0.00;    %Risk-free interest rate
F_0    = 1.1;    %Initial forward value

ModParams.v0     = 0.2;  %Inital volatility
ModParams.beta   = 0.7;  %Exponent
ModParams.alpha  = 0.08;  %Vol-vol
ModParams.rho    = 0; %Correlation

% ModParams.v0     = 0.25;  %Inital volatility
% ModParams.beta   = 0.6;  %Exponent
% ModParams.alpha     = 0.3;  %Vol-vol
% ModParams.rho     = -0.5; %Correlation

Kvec   = F_0*[0.6 0.8 0.90 0.95 0.999 1.05 1.10 1.2 1.4];
call   = 0;

%%% Price Strikes
prices = zeros(length(Kvec),1);
for k=1:length(Kvec)
    K = Kvec(k);
    prices(k) = SABR_European_AntonovApprox(F_0,K,T,call,r,ModParams);
    fprintf('%.8f\n', prices(k));
end


text:
function g = svargplvmPointGradient(x, model, y)
% VARGPLVMPOINTGRADIENT Wrapper function for gradient of a single point.
% FORMAT
% DESC is a wrapper function for the gradient of the log likelihood
% with respect to a point in the latent space. The GP-LVM
% model is one that is assumed to have already been trained.
% ARG x : the position in the latent space that is being optimised.
% ARG model : the trained GP-LVM model that is being optimised.
% ARG y : the position in data space for which the latent point is
% being optimised.
% RETURN g : the gradient of the log likelihood with respect to the
% latent position.
%
% SEEALSO : vargplvmPointLogLikeGradient, vargplvmOptimisePoint
%
% COPYRIGHT Michalis K. Titsias and Neil D. Lawrence, 2009, 2011
% COPYRIGHT Andreas C. Damianou, 2011

% VARGPLVM

if isfield(model, 'dynamics') && ~isempty(model.dynamics)
    if isfield(model.dynamics, 'reoptimise') && model.dynamics.reoptimise  %%% RE-OPT-CODE-NEW
        [vardistx, model] = vargplvmPartExpand(model, x); %%% RE-OPT-CODE-NEW
    else %%% RE-OPT-CODE-NEW
        % this is doing the expand
        x = reshape(x, model.N+size(y,1), model.dynamics.q*2);
        xtrain = x(1:model.N,:);
        xtest = x(model.N+1:end,:);
        model.dynamics.vardist = vardistExpandParam(model.dynamics.vardist, xtrain);
        vardistx = vardistExpandParam(model.vardistx, xtest);
        % end of expand
    end %%% RE-OPT-CODE-NEW
else
    vardistx = model.vardistx;
    vardistx = vardistExpandParam(vardistx, x);
end

g = - svargplvmPointLogLikeGradient(model, vardistx, y);

author: SheffieldML

repo: vargplvm


text:
function [success] = updatePkg(varargin)
	% Checks ciapkg version against GitHub repo and updates `ciapkg` and `+ciapkg` folders.
	% Biafra Ahanonu
	% started: 2020.08.18 [11:16:56]
	% inputs
		%
	% outputs
		%
	% Notes
		% Follows https://semver.org convention.

	% changelog
		% 2021.03.09 [21:34:52] - Option on whether to update package (only 0 for now). Also alert user if behind a version.
		% 2021.06.21 [16:18:07] - Updated to ciatah from calciumImagingAnalysis URLs.
		% 2021.07.06 [11:33:22] - Updated comparison to ensure ignoring of lower level versions after higher level versions are found to differ.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		% Only check once every couple of days to warn user instead of every time.

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% Char cell array: list of directories to remove then update, lead private, _external_programs, and data directories intact.
	options.removeDirs = {...
		'@calciumImagingAnalysis',...
		'@ciatah',...
		'+ciapkg',...
		'ciapkg',...
		'docs',...
		'file_exchange'};
	% [IGNORE] Char: GitHub API URL to VERSION file on CIAPKG repository.
	options.versionURL = {'https://api.github.com/repos/bahanonu/ciatah/contents/+ciapkg/VERSION','https://api.github.com/repos/bahanonu/ciatah/contents/ciapkg/VERSION','https://api.github.com/repos/bahanonu/ciatah/contents/VERSION'};
	% Binary: 1 = pop-up GUI enabled
	options.showGui = 1;
	% Binary: 1 = update code, 0 = only check for update and notify user behind a version.
	options.updatePackage = 0;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		success = 0;

		% Versions follow MAJOR.MINOR.PATCH convention.
		% Get online CIAPKG version
		[onlineVersion, onlineDateTimeStr] = ciapkg.versionOnline;
		% Get local CIAPKG version
		[currentVersion, currentDateTimeStr] = ciapkg.version;
		% currentVersion = 'v3.19.20201021135616'
		% currentVersion = 'v3.20.20201019131033';

		% Compare to see if running an old version.
		verCompare = subfxn_verCompare(currentVersion,onlineVersion);

		if isnan(verCompare)==1
			verInfoStr = 'Cannot obtain online version, check internet or update manually.';
			verInfoStr2 = verInfoStr;
		elseif verCompare==1
			verInfoStr = 'Running most up-to-date version!';
			verInfoStr2 = 'Running most up-to-date version!';
		elseif verCompare<1
			verInfoStr = 'Running behind! Initiating update [IGNORE for now].';
			verInfoStr2 = 'Running a version behind, consider updating.';
		elseif verCompare>1
			verInfoStr = 'I do not know how, but you are running a version ahead! [Dev build]';
			verInfoStr2 = 'I do not know how, but you are running a version ahead! [Dev build]';
		end

		if options.updatePackage==0
			if iscell(onlineVersion)
				onlineVersionTmp = onlineVersion{1};
			else
				onlineVersionTmp = onlineVersion;
			end
			% verInfoStr = sprintf('%s\n Local version: %s.\n Online version %s.\n',verInfoStr2,currentVersion,onlineVersionTmp);
			verInfoStr = sprintf('%s\n Local version: %s (%s).\n Online version %s (%s).\n',verInfoStr2,currentVersion,currentDateTimeStr,onlineVersionTmp,onlineDateTimeStr);
			warning(verInfoStr)
			if verCompare<1
				ciapkg.overloaded.msgbox(verInfoStr,'Note on CIAtah version.')
			end
			success = 1;
			return
		end

		disp(verInfoStr)
		if options.showGui==1
			uiwait(ciapkg.overloaded.msgbox(verInfoStr));
		end

		% If the user is behind, ask if they want to update then initiate download and update from online repository
		if verCompare<1
			% Get a list of directories to remove, make sure

			% Confirm with the user

			% rmdir
			% nFolders = length(options.removeDirs);

			% for folderNo = 1:nFolders
			% end
		end
		success = 1;
	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
	function [verCompare,compareVector] = subfxn_verCompare(verId1,verId2)
		fprintf('Comparing\nLocal:  %s.\nOnline: %s.\n',verId1,verId2);

		if isempty(verId1)==1|isempty(verId2)==1
			verCompare = NaN;
			compareVector = [];
			return;
		end

		if iscell(verId1)
			verId1 = verId1{1};
		end
		if iscell(verId2)
			verId2 = verId2{1};
		end
		% Remove version tag
		verId1 = strrep(verId1,'v','');
		verId2 = strrep(verId2,'v','');

		% Compare two versions assuming major.minor.patch.build where in each the numbers always increment in a positive manner.
		verId1 = cellfun(@(x) str2num(x),strsplit(verId1,'.'));
		verId2 = cellfun(@(x) str2num(x),strsplit(verId2,'.'));
		compareVector = zeros([1 max(length(verId1),length(verId2))]);
		nLevels = length(compareVector);

		% Following MAJOR.MINOR.PATCH convention
		% Start at the highest level and only continue comparing while version 2 is greater than or equal to version 1.
		lockOut = 0; % 1 = higher level is an old version, so only calculate version difference
		for verLvl = 1:nLevels
			compareVector(verLvl) = verId1(verLvl) - verId2(verLvl);
			% If higher version is older or newer, lower version is irrelevant so skip.
			if lockOut==1
				compareVector(verLvl) = 0;
			end
			if compareVector(verLvl)<0 && lockOut==0 % This level indicates an old version
				verCompare = 0;
				lockOut = 1;
			elseif compareVector(verLvl)>0 && lockOut==0 % This level indicates this level is a newer version
				lockOut = 1;
			end
		end

		% User is running the same exact version.
		if sum(compareVector)==0
			verCompare = 1;
		end

		% If user is running ahead, should only occur in development builds.
		if sum(compareVector)>0
			verCompare = 2;
		end
		% disp(compareVector)
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function varargout = spatialdecomposition3d(x_D,varargin)

if ~check_option(varargin,'voronoi') && size(x_D,1) > 1 %  check_option(varargin,'unitcell') &&

  unitCell = varargin{1};

  dX = abs(2*unitCell([1 13 17]));

  % grid coordinates

  iX = bsxfun(@minus,x_D,min(x_D));
  iX = 1+round(bsxfun(@rdivide,iX,dX));

  sz = max(iX);  %extent, number of voxels in each direction

  % generate a tetragonal unit cell
  id = s2i(sz,iX);
  id = [id(:) id(:)];

  el = []; er = [];
  for dim=1:3
    nX = iX;    pX = iX;
    nX(:,dim) = nX(:,dim)+1;
    pX(:,dim) = pX(:,dim)-1;

    ex = [s2i(sz,nX) s2i(sz,pX)];

    del = pX(:,dim) < 1 | nX(:,dim) > sz(dim);
    el = [el;id(~del,:)];
    er = [er;ex(~del,:)];
  end

  % adjacency of voxels
  varargout{1} = el(:);
  varargout{2} = er(:);
  clear el er
  varargout{3} = sz;
  varargout{4} = dX;
  varargout{5} = min(x_D);

elseif ~check_option(varargin,'voronoi')

  sz = x_D;

  nd = prod(double(sz));
  nf = prod(double(sz+1));

  ind = varargin{1};
  dz = varargin{2};
  lz = double(varargin{3});

  %  adding boundary voxels
  [bind{1} bf{1}] = constr(sz,[0 0 0],0,nf);
  [bind{2} bf{2}] = constr(sz,[1 0 0],0,nf);
  [bind{3} bf{3}] = constr(sz,[0 0 0],2,nf);
  [bind{4} bf{4}] = constr(sz,[0 1 0],2,nf);
  [bind{5} bf{5}] = constr(sz,[0 0 0],4,nf);
  [bind{6} bf{6}] = constr(sz,[0 0 1],4,nf);

  aind = unique([uint32(ind); vertcat(bind{:})]);
  [ix,iy,iz] = ind2sub(sz,aind);
  clear bind

  ixn = ix+1; iyn = iy+1; izn = iz+1; % next voxel index

  % pointels incident to voxels
  VD = [s2i(sz+1,ix,iy,iz) s2i(sz+1,ixn,iy,iz)  s2i(sz+1,ixn,iyn,iz)  s2i(sz+1,ix,iyn,iz) ....
    s2i(sz+1,ix,iy,izn) s2i(sz+1,ixn,iy,izn) s2i(sz+1,ixn,iyn,izn) s2i(sz+1,ix,iyn,izn)];

  % new voxel coordinates
  vind = unique(VD(:));
  [vx,vy,vz] = ind2sub(sz+1,vind);
  varargout{1}(vind,:) = [double(vx-1)*dz(1)-dz(1)/2+lz(1),...
    double(vy-1)*dz(2)-dz(2)/2+lz(2),...
    double(vz-1)*dz(3)-dz(3)/2+lz(3)];
  clear vx vy vz dx dy dz lx ly lz

  % surfels incident to voxel
  FD = [s2i(sz+1,ix,iy,iz) s2i(sz+1,ixn,iy,iz) ...
    s2i(sz+1,ix,iy,iz)+2*nf s2i(sz+1,ix,iyn,iz)+2*nf ...
    s2i(sz+1,ix,iy,iz)+4*nf s2i(sz+1,ix,iy,izn)+4*nf];

  clear ix iy iz ixn iyn izn

  % pointels incident to facets
  tri = [ ...
    8 5 1 4  %4 1 5 8
    2 3 7 6
    5 6 2 1  %1 2 6 5
    3 4 8 7
    1 2 3 4  %4 3 2 1
    5 6 7 8];

  VF = zeros(nf,4,'uint32');
  VF(FD,:) = reshape(VD(:,tri),[],4);
  varargout{2} = VF;
  %   clear VD VF

  id = repmat(double(aind(:)),1,6);
  % surfels as incidence matrix

  FD = double(FD);
  o = ones(size(FD));
  o(:,1:2:6) = -1;

  FD = sparse(FD,id,o,6*nf,nd);
  varargout{3} = FD;

  bf = double(vertcat(bf{:}));

  [F,D,val] = find(FD(bf,:));
  varargout{4} =  sparse(bf(F),D,val,6*nf,nd);

else
  boundary = get_option(varargin,'boundary','cube');

  % extrapolate dummy coordinates %dirty
  if ~check_option(varargin,'3d')
    switch lower(boundary)
      case 'cube'
        v = get_option(varargin,'dx',ceil(nthroot(size(x_D,1),3)));
        a = min(x_D);
        b = max(x_D);
        d = (b-a)./v;

        dx = linspace(a(1),b(1),v);
        dy = linspace(a(2),b(2),v);
        dz = linspace(a(3),b(3),v);

        [x y] = meshgrid(dx,dy);
        z1 = [x(:) y(:)];
        z2 = z1;
        z1(:,3) = a(3)- d(3);
        z2(:,3) = b(3)+ d(3);

        [x z] = meshgrid(dx,dz);
        y1(:,[1 3]) = [x(:) z(:)];
        y2 = y1;
        y1(:,2) = a(2)- d(2);
        y2(:,2) = b(2)+ d(2);

        [y z] = meshgrid(dy,dz);
        x1(:,2:3) = [y(:) z(:)];
        x2 = x1;
        x1(:,1) = a(1)- d(1);
        x2(:,1) = b(1)+ d(1);

        dummy = [x1; x2; y1; y2; z1; z2];
      case 'sphere'
        dx = (max(x_D(:,1)) - min(x_D(:,1)));
        dy = (max(x_D(:,2)) - min(x_D(:,2)));
        dz = (max(x_D(:,3)) - min(x_D(:,3)));

        [x y z] = sphere;

        dummy = [x(:).*dx+dx/2+min(x_D(:,1)) ...
          y(:).*dy+dy/2+min(x_D(:,2)) ...
          z(:).*dz+dz/2+min(x_D(:,3)) ];

        dummy = unique(dummy,'rows');
      otherwise
        error('Uknown boundary type.')
    end
  else
    dx = (max(xy(:,1)) - min(xy(:,1)));
    dy = (max(xy(:,2)) - min(xy(:,2)));
    dz = (max(xy(:,3)) - min(xy(:,3)));

    [x y z] = sphere;

    dummy = [x(:).*dx+dx/2 ...
      y(:).*dy+dy/2 ...
      z(:).*dz+dz/2 ];

    dummy = unique(dummy,'rows');

  end
  n = size(x_D,1);
  x_D = [x_D; dummy];
  [v c] = voronoin(x_D,{'Q7','Q8','Q5','Q3','Qz'});   %Qf {'Qf'} ,{'Q7'}
  c(end-length(dummy)+1:end) = [];

  FD = cell(numel(c),1); VF = FD; % preallocate
  for k=1:numel(c)
    vertex = c{k}; s = [];
    tri = convhulln(v(vertex,:));
    s(1:size(tri,1),1) = k;

    VF{k} = vertex(tri);
    FD{k} = s;
  end

  VF = vertcat(VF{:});
  FD = vertcat(FD{:});
  [VF b faceID] = unique(sort(VF,2),'rows');
  FD = sparse(faceID(:),FD(:),1,max(faceID(:)),n);

  A = triu(FD'*FD,1);
  [varargout{1},varargout{2}] = find(A|A');

  varargout{3} = v;
  varargout{4} = VF;
  varargout{5} = FD;
end

function [bind,bf] = constr(sz,s,d,n)

dx = [1 sz(1)];
dy = [1 sz(2)];
dz = [1 sz(3)];

switch d
  case 0
    dx(~any(s)+1) = dx(any(s)+1);
  case 2
    dy(~any(s)+1) = dy(any(s)+1);
  case 4
    dz(~any(s)+1) = dz(any(s)+1);
end

[bx by bz] = meshgrid(dx(1):dx(2),dy(1):dy(2),dz(1):dz(2));
bind = s2i(sz,bx(:),by(:),bz(:));

[ix,iy,iz] = ind2sub(sz,bind);
bf = s2i(sz+1,ix+s(1),iy+s(2),iz+s(3))+d*n;

function ndx = s2i(sz,ix,iy,iz)
% faster version of sub2ind
if nargin == 4
  ndx = 1 + (ix-1) + (iy-1)*sz(1) +(iz-1)*sz(1)*sz(2);
elseif nargin == 2 && size(ix,2) == 3
  ndx = 1 + (ix(:,1)-1) + (ix(:,2)-1)*sz(1) +(ix(:,3)-1)*sz(1)*sz(2);
end

author: mtex-toolbox

repo: mtex

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


function results = OTB_DEEP_settings(seq, res_path, bSaveImage, parameters)

close all

s_frames = seq.s_frames;

% Feature specific parameters
hog_params.cell_size = 4;
hog_params.compressed_dim = 10;

% grayscale_params.colorspace='gray';
% grayscale_params.cell_size = 1;

cn_params.tablename = 'CNnorm';
cn_params.useForGray = false;
cn_params.cell_size = 4;
cn_params.compressed_dim = 3;

ic_params.tablename = 'intensityChannelNorm6';
ic_params.useForColor = false;
ic_params.cell_size = 4;
ic_params.compressed_dim = 3;

cnn_params.nn_name = 'imagenet-vgg-m-2048.mat'; % Name of the network
cnn_params.output_layer = [3 14];               % Which layers to use
cnn_params.downsample_factor = [2 1];           % How much to downsample each output layer
cnn_params.compressed_dim = [16 64];            % Compressed dimensionality of each output layer
cnn_params.input_size_mode = 'adaptive';        % How to choose the sample size
cnn_params.input_size_scale = 1;                % Extra scale factor of the input samples to the network (1 is no scaling)

% Which features to include
params.t_features = {
    struct('getFeature',@get_cnn_layers, 'fparams',cnn_params),...
    ...struct('getFeature',@get_colorspace, 'fparams',grayscale_params),...
    struct('getFeature',@get_fhog,'fparams',hog_params),...
    struct('getFeature',@get_table_feature, 'fparams',cn_params),...
    struct('getFeature',@get_table_feature, 'fparams',ic_params),...
};

% Global feature parameters1s
params.t_global.normalize_power = 2;    % Lp normalization with this p
params.t_global.normalize_size = true;  % Also normalize with respect to the spatial size of the feature
params.t_global.normalize_dim = true;   % Also normalize with respect to the dimensionality of the feature

% Image sample parameters
params.search_area_shape = 'square';    % The shape of the samples
params.search_area_scale = 4.0;         % The scaling of the target size to get the search area
params.min_image_sample_size = 200^2;   % Minimum area of image samples
params.max_image_sample_size = 250^2;   % Maximum area of image samples

% Detection parameters
params.refinement_iterations = 1;       % Number of iterations used to refine the resulting position in a frame
params.newton_iterations = 5;           % The number of Newton iterations used for optimizing the detection score
params.clamp_position = false;          % Clamp the target position to be inside the image

% Learning parameters
params.output_sigma_factor = 1/12;		% Label function sigma
params.learning_rate = 0.012;	 	 	% Learning rate
params.nSamples = 50;                   % Maximum number of stored training samples
params.sample_replace_strategy = 'lowest_prior';    % Which sample to replace when the memory is full
params.lt_size = 0;                     % The size of the long-term memory (where all samples have equal weight)
params.train_gap = 5;                   % The number of intermediate frames with no training (0 corresponds to training every frame)
params.skip_after_frame = 1;            % After which frame number the sparse update scheme should start (1 is directly)
params.use_detection_sample = true;     % Use the sample that was extracted at the detection stage also for learning

% Factorized convolution parameters
params.use_projection_matrix = true;    % Use projection matrix, i.e. use the factorized convolution formulation
params.update_projection_matrix = true; % Whether the projection matrix should be optimized or not
params.proj_init_method = 'pca';        % Method for initializing the projection matrix
params.projection_reg = 2e-7;	 	 	% Regularization paremeter of the projection matrix

% Generative sample space model parameters
params.use_sample_merge = true;                 % Use the generative sample space model to merge samples
params.sample_update_criteria = 'Merge';        % Strategy for updating the samples
params.weight_update_criteria = 'WeightedAdd';  % Strategy for updating the distance matrix
params.neglect_higher_frequency = false;        % Neglect hiigher frequency components in the distance comparison for speed

% Conjugate Gradient parameters
params.CG_iter = 5;                     % The number of Conjugate Gradient iterations in each update after the first frame
params.init_CG_iter = 10*20;            % The total number of Conjugate Gradient iterations used in the first frame
params.init_GN_iter = 10;               % The number of Gauss-Newton iterations used in the first frame (only if the projection matrix is updated)
params.CG_use_FR = false;               % Use the Fletcher-Reeves (true) or Polak-Ribiere (false) formula in the Conjugate Gradient
params.CG_standard_alpha = true;        % Use the standard formula for computing the step length in Conjugate Gradient
params.CG_forgetting_rate = 75;	 	 	% Forgetting rate of the last conjugate direction
params.precond_data_param = 0.7;	 	% Weight of the data term in the preconditioner
params.precond_reg_param = 0.1;	 	% Weight of the regularization term in the preconditioner
params.precond_proj_param = 30;	 	 	% Weight of the projection matrix part in the preconditioner

% Regularization window parameters
params.use_reg_window = true;           % Use spatial regularization or not
params.reg_window_min = 1e-4;			% The minimum value of the regularization window
params.reg_window_edge = 10e-3;         % The impact of the spatial regularization
params.reg_window_power = 2;            % The degree of the polynomial to use (e.g. 2 is a quadratic window)
params.reg_sparsity_threshold = 0.12;   % A relative threshold of which DFT coefficients that should be set to zero

% Interpolation parameters
params.interpolation_method = 'bicubic';    % The kind of interpolation kernel
params.interpolation_bicubic_a = -0.75;     % The parameter for the bicubic interpolation kernel
params.interpolation_centering = true;      % Center the kernel at the feature sample
params.interpolation_windowing = false;     % Do additional windowing on the Fourier coefficients of the kernel

% Scale parameters for the translation model
% Only used if: params.use_scale_filter = false
params.number_of_scales = 7;            % Number of scales to run the detector
params.scale_step = 1.01;               % The scale factor

% Scale filter parameters
% Only used if: params.use_scale_filter = true
params.use_scale_filter = false;            % Use the fDSST scale filter or not (for speed)
% params.scale_sigma_factor = 1/16;         % Scale label function sigma
% params.scale_learning_rate = 0.025;       % Scale filter learning rate
% params.number_of_scales_filter = 17;    % Number of scales
% params.number_of_interp_scales = 33;    % Number of interpolated scales
% params.scale_model_factor = 1.0;        % Scaling of the scale model
% params.scale_step_filter = 1.02;        % The scale factor for the scale filter
% params.scale_model_max_area = 32*16;    % Maximume area for the scale sample patch
% params.scale_feature = 'HOG4';          % Features for the scale filter (only HOG4 supported)
% params.s_num_compressed_dim = 'MAX';    % Number of compressed feature dimensions in the scale filter
% params.lambda = 1e-2;                   % Scale filter regularization
% params.do_poly_interp = true;           % Do 2nd order polynomial interpolation to obtain more accurate scale

% Other parameters
params.visualization = 0;               % Visualiza tracking and detection scores
params.debug = 0;                       % Do full debug visualization

% Initialize
params.init_sz = [seq.init_rect(1,4), seq.init_rect(1,3)];
params.init_pos = [seq.init_rect(1,2), seq.init_rect(1,1)] + (params.init_sz - 1)/2;
params.s_frames = s_frames;

% Run tracker
results = tracker(params);


text:
function [x,w,P]=GaussRadau(kk)
N = kk-1; % Compute for the number of points kk

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lgrnodes.m
%
% Computes the Legendre-Gauss-Radau nodes, weights and the LGR Vandermonde
% matrix. The LGR nodes are the zeros of P_N(x)+P_{N+1}(x).
%
% References on LGR nodes and weights:
%   C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Tang, "Spectral Methods
%   in Fluid Dynamics," Section 2.3. Springer-Verlag 1987
%
%   F. B. Hildebrand , "Introduction to Numerical Analysis," Section 8.11
%   Dover 1987
%
% Written by Greg von Winckel - 05/02/2004
% Contact: gregvw@chtm.unm.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Truncation + 1
N1=N+1;

% Use Chebyshev-Gauss-Radau nodes as initial guess for LGR nodes
x=-cos(2*pi*(0:N)/(2*N+1))';

% The Legendre Vandermonde Matrix
P=zeros(N1,N1+1);

% Compute P_(N) using the recursion relation
% Compute its first and second derivatives and
% update x using the Newton-Raphson method.

xold=2;

% Free abscissae
free=2:N1;

while max(abs(x-xold))>eps

    xold=x;

    P(1,:)=(-1).^(0:N1);

    P(free,1)=1;    P(free,2)=x(free);

    for k=2:N1
        P(free,k+1)=( (2*k-1)*x(free).*P(free,k)-(k-1)*P(free,k-1) )/k;
    end

    x(free)=xold(free)-((1-xold(free))/N1).*(P(free,N1)+P(free,N1+1))...
        ./(P(free,N1)-P(free,N1+1));
end

% The Legendre-Gauss-Radau Vandermonde
P=P(1:N1,1:N1);

% Compute the weights
w=zeros(N1,1);
w(1)=2/N1^2;
w(free)=(1-x(free))./(N1*P(free,N1)).^2;

repo: Aero-matlab


text:
function value = r8vec_amax ( n, a )

%*****************************************************************************80
%
%% R8VEC_AMAX returns the maximum absolute value in an R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(N), the array.
%
%    Output, real VALUE, the value of the entry
%    of largest magnitude.
%
  value = max ( abs ( a(1:n) ) );

  return
end


text:
function mean = normal_01_mean ( )

%*****************************************************************************80
%
%% NORMAL_01_MEAN returns the mean of the Normal 01 PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = 0.0;

  return
end


text:
function [both, pos, neg] = visualizePathwayInEpistasis(Nall, radius, pathwayNames)
% This function creates a circular network such that all the nodes are
% arranged in a circle. The nodes represent a pathway, the color of the
% node represents whether the interactions within a pathway are none, positive,
% negative, or mixed. The size of the node represents the number of interactions.
% The edges represent which pathways are sharing which
% type of interaction. The color of the edge represents if the interactions
% are positive negative or mixed. The with of the edge represent the total number
% of interactions between any two given pathways.
%
% USAGE:
%
%     [both,pos,neg]=visualizePathwayInEpistasis(Nall,radius,pathwayNames)
%
% INPUT:
%     Nall:          a structure representing the epistatic interaction networks:
%                    Nall.pos: a square matrix representing number of positive interactions shared by any two pathways.
%                    Nall.neg: a square matrix representing number of negative interactions shared by any two pathways.
%     radius:        radius of the generic circles for clock diagram.
%     pathwayNames:  subSystem names corresponding to matrices in Nall
%                    structure
%
% OUTPUT:
%     both:          contains both positive and negative interactions
%     pos:           contains only positive interactions
%     neg:           contains only negative interactions
%
% NOTE:
%    See figures in following publication:
%    Joshi CJ and Prasad A, 2014, "Epistatic interactions among metabolic genes
%    depend upon environmental conditions", Mol. BioSyst., 10, 2578-2589.
%
% .. Authors:
%      - Chintan Joshi 10/26/2018
%      - Chintan Joshi; made modifications to label the nodes using pathway names, 10/26/2018

if nargin < 2
    radius = 40;
end

pos = Nall.pos;
neg = Nall.neg;
both = Nall.pos + Nall.neg;
indn = []; indp = []; mix_ = [];
ps = find(diag(pos) ~= 0); ns = find(diag(neg) ~= 0); nps = find(diag(both) ~= 0);
if ~isempty(nps) && ~isempty(ps)
    indp = intersect(ps, nps);
end
if ~isempty(nps) && ~isempty(ns)
    indn = intersect(ns, nps);
end
if isempty(indn)
    p_ = indp; n_ = []; mix_ = [];
elseif isempty(indp)
    p_ = []; n_ = indn; mix_ = [];
else
    p_ = setdiff(indp, indn); n_ = setdiff(indn, indp); mix_ = intersect(indp, indn);
end
% keeps only those that have both but not either type of interactions alone
c = 0;
for i = 1:length(both)
    for j = 1:length(both)
        if both(i, j) ~= 0 && pos(i, j) ~= 0 && neg(i, j) ~= 0
            c = c + 1;
            indm1(c, 1) = i;
            indm2(c, 1) = j;
        end
    end
end
% conversion to circular cordinates
theta = linspace(0, 2 * pi, length(pos) + 1);
theta = theta(1:end - 1);
[x, y] = pol2cart(theta, 1);
tx = x; ty = y;
x = x * 0.9; y = y * 0.9;
[indp1, indp2] = ind2sub(size(pos), find(pos(:)));
[indn1, indn2] = ind2sub(size(neg), find(neg(:)));
gc = [0.23 0.44 0.34];
rc = [0.58 0.39 0.39];
yc = [255 204 102] / 255;
subplot(3, 3, 5);
% plot the structure of the clock diagram
plot(x, y, '.', 'MarkerEdgeColor', [150 150 150]./255, 'markersize', 1); hold on
% label the the texts for each point in clock diagram
for i = 1:length(x)
    textAngle = (i - 1) * 360 / length(x);  % calculate the angle for text
    text(tx(i), ty(i), pathwayNames(i), 'Rotation', textAngle);
end
if ~isempty(indp1)
    arrayfun(@(p, q)line([x(p), x(q)], [y(p), y(q)], 'Color', gc, 'LineWidth', pos(p, q)), indp1, indp2);  % plot inter-pathway positive interactions
end
if ~isempty(indn1)
    arrayfun(@(p, q)line([x(p), x(q)], [y(p), y(q)], 'Color', rc, 'LineWidth', neg(p, q)), indn1, indn2);  % plot inter-pathway negative interactions
end
if exist('indm1')
    arrayfun(@(p, q)line([x(p), x(q)], [y(p), y(q)], 'Color', yc, 'LineWidth', both(p, q)), indm1, indm2);  % plot inter pathway both types of interactions
end
% plot intra-pathway interactions
plot(x, y, '.', 'MarkerEdgeColor', [150 150 150]./255, 'markersize', radius);
plot(x(p_), y(p_), '.', 'MarkerEdgeColor', gc, 'markersize', radius + length(p_) * 3);
plot(x(n_), y(n_), '.', 'MarkerEdgeColor', rc, 'markersize', radius + length(n_) * 3);
plot(x(mix_), y(mix_), '.', 'MarkerEdgeColor', yc, 'markersize', radius + length(mix_) * 3);
set(gca, 'DataAspectRatio', [1 1 1]);

axis equal off;
hold off;

author: opencobra

repo: cobratoolbox


text:
function i4vec_index_insert_test ( )

%*****************************************************************************80
%
%% I4VEC_INDEX_INSERT_TEST tests I4VEC_INDEX_INSERT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n_max = 25;
  n = 0;
  x = [];
  indx = [];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_INDEX_INSERT_TEST\n' );
  fprintf ( 1, '  I4VEC_INDEX_INSERT inserts values into an\n' );
  fprintf ( 1, '  index sorted array of integers.\n' );

  xval = 8;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  xval = 7;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  b = 0;
  c = 20;
  seed = 123456789;

  for i = 1 : 20
    [ xval, seed ] = i4_uniform_ab ( b, c, seed );
    fprintf ( 1, '  %6d\n', xval );
    [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );
  end

  xval = 7;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  xval = 8;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indexed list of entries:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(I)  X(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %6d  %6d  %6d\n', i, indx(i), x(i), x(indx(i)) );
  end

  return
end


text:
% Test file for trigtech/max.m

function pass = test_max(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

%%
% Spot-check the extrema for a few functions.
pass(1) = test_spotcheck_max(testclass, @(x) exp(-cos(2*pi*x)), exp(1), pref);
pass(2) = test_spotcheck_max(testclass, @(x) sin(10*pi*x), 1, pref);

pass(3) = test_spotcheck_max(testclass, @(x) exp(sin(pi*x).^100), exp(1), pref);
pass(4) = test_spotcheck_max(testclass, @(x) exp(-sin(pi*x).^100), 1, pref);

% Approx to sign function
pass(5) = test_spotcheck_max(testclass, @(x) 4/pi*(sin(pi*x) + ...
    1/3*sin(3*pi*x) + 1/5*sin(5*pi*x) + 1/7*sin(7*pi*x) + 1/9*sin(9*pi*x)), 1.182328208857607, pref);

%%
% Check operation for array-valued inputs.
fun_op = @(x) [exp(-cos(2*pi*x)) sin(10*pi*x) exp(-sin(pi*(x-0.32)).^100)];
f = testclass.make(fun_op, [], pref);
[y, x] = max(f);
exact_max = [exp(1) 1 1];
fx = [exp(-cos(2*pi*x(1))) sin(10*pi*x(2)) exp(-sin(pi*(x(3)-0.32)).^100)];
pass(6) = (all(abs(y - exact_max) < 100*eps) && ...
           all(abs(fx - exact_max) < 10*eps));

%%
% Test for complex-valued TRIGTECH objects.
pass(7) = test_spotcheck_max(testclass, ...
    @(x) 1 + cos(pi*x) + exp(1i*pi*x), ...
    3, pref);

end

% Spot-check the results for a given function.
function result = test_spotcheck_max(testclass, fun_op, exact_max, pref)

f = testclass.make(fun_op,[], pref);
[y, x] = max(f);
fx = fun_op(x);
result = (all(abs(y - exact_max) < 100*vscale(f)*eps) && ...
          all(abs(fx - exact_max) < 100*vscale(f)*eps));

end


text:
function setRavenSolver(solver)
% setRavenSolver
%   Sets the solver for RAVEN and optionally saves a default value
%	to MATLAB startup.
%
%   solver  string with solver name. The following options are available:
%           glpk    uses RAVEN-provided binaries
%           gurobi  requires a working Gurobi installation
%           cobra   requires a working COBRA Toolbox installation. RAVEN
%                   parses the problem to COBRA and uses the solver that is
%                   set by changeCobraSolver.
%
%   Usage: setRavenSolver(solver)
solver=char(solver);
switch solver
    case 'cobra'
        global CBT_LP_SOLVER;
        if isempty(CBT_LP_SOLVER)
            error('COBRA toolbox is not initialized, run ''initCobraToolbox()''')
        else
            setpref('RAVEN','solver','cobra');
        end
    case 'mosek'
        error('MOSEK support has been discontinued since RAVEN 2.3.0.')
    case {'glpk','gurobi','soplex'}
        setpref('RAVEN','solver',solver)
    otherwise
        error('Invalid solver defined')
end
global RAVENSOLVER
RAVENSOLVER = solver;
end

author: SysBioChalmers

repo: RAVEN


text:
function [gK_uf, gK_ff] = gpsimCandidateCovGrads(model, M)

% GPSIMCANDIDATECOVGRADS Sparse objective function gradients wrt Covariance function.
% FORMAT
% DESC gives the gradients of the log likelihood with respect to the
% components of the posterior covariance.
% ARG model : the model for which the gradients are to be computed.
% ARG M : The training data for which the computation is to be made
% RETURN gK_uf : the gradient of the likelihood with respect to the
% elements of K_uf.
% RETURN gK_ff : the gradient of the likelihood with respect to K_ff
%
% COPYRIGHT : Neil D. Lawrence, 2007
%
% SEEALSO : gpsimCreate, gpsimAddCandidate, gpsimCandidateLogLikeGradient

% SHEFFIELDML
%E = model.candidate.K_uf*M;
E = model.candidate.K_uf*model.candidate.invK*M;
EET = E*E';
AinvE = model.candidate.Ainv*E;
%AinvEET = model.candidate.Ainv*EET;
diagK_fuAinvEMT = sum(model.candidate.K_uf.*(model.candidate.Ainv*E*M'), 1)';
K_fuAinvEMT = model.candidate.K_uf'*model.candidate.Ainv*E*M';
AinvEETAinv = AinvE*AinvE';
K_ufdAinvplusAinvEETAinvK_fu = model.candidate.K_uf'*(model.candidate.Ainv ...
                                                  + AinvEETAinv)*model.candidate.K_uf;
invK_uuK_uf = model.invK*model.candidate.K_uf;
invK_uuK_ufDinv = invK_uuK_uf*model.candidate.invK;
MMT = M*M';
Q = -model.candidate.K + MMT ...
    + K_ufdAinvplusAinvEETAinvK_fu...
    -K_fuAinvEMT - K_fuAinvEMT';
%gK_uu = 0.5*((model.invK ...
%              -model.candidate.Ainv) - AinvEETAinv ...
%             + invK_uuK_ufDinv*Q*invK_uuK_ufDinv');
gK_uf = -invK_uuK_ufDinv*Q*model.candidate.invK ...
        -model.candidate.Ainv*model.candidate.K_uf*model.candidate.invK ...
        -AinvEETAinv*model.candidate.K_uf*model.candidate.invK ...
        +model.candidate.Ainv*E*M'*model.candidate.invK;
gK_ff = 0.5*model.candidate.invK*Q*model.candidate.invK;

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_snirf

% WALLTIME 00:10:00
% MEM 2GB
% DEPENDENCY homer2fieldtrip fieldtrip2homer ft_write_data

p = tempdir;
f1 = fullfile(p, 'data1.snirf');
f2 = fullfile(p, 'data2.snirf');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% test the conversion of a Homer file to snirf

% the Homer nirs file has the 's' channel, which is not present in the snirf representation

% read the original data
cfg = [];
cfg.dataset = dccnpath('/home/common/matlab/fieldtrip/data/test/original/nirs/homer/S1001_run01.nirs');
cfg.channel = {'all', '-s*'};
data1 = ft_preprocessing(cfg);
event1 = ft_read_event(cfg.dataset);

% write the data to disk and convert to snirf
chanindx = find(ismember(data1.hdr.label, data1.label));
ft_write_data(f1, data1.trial{1}, 'header', data1.hdr, 'chanindx', chanindx, 'event', event1)

% read the converted data
cfg = [];
cfg.dataset = f1;
data1a = ft_preprocessing(cfg);
event1a = ft_read_event(cfg.dataset);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% test the conversion of an Artinis file to snirf

% read the original data
cfg = [];
cfg.dataset = dccnpath('/home/common/matlab/fieldtrip/data/test/original/nirs/artinis/Helena/190528_fingertap_L.oxy3');
cd(fileparts(cfg.dataset)); % we should be located in the directory with the optodetemplates.xml
data2 = ft_preprocessing(cfg);
event2 = ft_read_event(cfg.dataset, 'chanindx', -1); % do not parse the ADC channels

% replace the original data of the first sample by channel numbers. This helps to assert
% whether data2a is still the same as the original data (data2), even if the
% channel order is switched.
data2.trial{1}(:,1) = [1:length(data2.label)]';

% write the data to disk and convert to snirf
ft_write_data(f2, data2.trial{1}, 'header', data2.hdr, 'event', event2)

% read the converted data
cfg = [];
cfg.dataset = f2;
data2a = ft_preprocessing(cfg);
event2a = ft_read_event(cfg.dataset);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% clean up

delete(f1)
delete(f2)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% compare all fields

assert(isalmostequal(data1.time,        data1a.time, 'abstol', 1e-8));
assert(isalmostequal(data1.fsample,     data1a.fsample, 'abstol', 1e-8));
assert(isalmostequal(data1.sampleinfo,  data1a.sampleinfo, 'abstol', 1e-8));
assert(isequal(data1.label, data1a.label));
assert(isalmostequal(data1.trial,       data1a.trial, 'abstol', 1e-8));

% the duration of the Homer and snirf events is not the same
assert(isequal({event1.type},   {event1a.type}));
assert(isequal([event1.sample], [event1a.sample]));
assert(isequal([event1.value],  [event1a.value]));

%%

% do not compare the label field for Artinis, it differs in Tx/Rx versus S/D
assert(isalmostequal(data2.time,        data2a.time, 'abstol', 1e-8));
assert(isalmostequal(data2.fsample,     data2a.fsample, 'abstol', 1e-8));
assert(isalmostequal(data2.sampleinfo,  data2a.sampleinfo, 'abstol', 1e-8));
assert(length(intersect(data2.trial{1}(:,1), data2a.trial{1}(:,1)))==length(data2.trial{1}(:,1))); % the first sample should contain 1:nchan (see here above)

% the channel names are not the same, neither for the nirs channels, not for the ADC/AUX channels
% assert(isequal(data2.label, data2a.label));

% the type and value of the Artinis and snirf events are not the same
assert(isequal([event2.sample], [event2a.sample]));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function handles = setupPanels(handles)
%SETUPPANELS    Populate the panels on CHEBGUI

% Call the subfunctions to populate the panels.
handles = chebguiController.setupPanelFigures(handles);
handles = chebguiController.setupPanelInput(handles);
handles = chebguiController.setupPanelType(handles);
handles = chebguiController.setupPanelDiscretization(handles);
handles = chebguiController.setupPanelIVPsolver(handles);

end


text:
function nWritten = export_GSPS_IOD(planC, filenameRoot, filenumber)
%"export_GSPS_IOD"
%   Writes to disk all scans contained in planC by constructing the modules
%   specified in part A.3.3 of PS 3.3 of the 2006 DICOM standard.  The
%   planC MUST have first been run through the function:
%
%   generate_DICOM_UID_Relationships.m.
%
%   filenameRoot is the the base path + filename to have file numbers
%   appended to, ie, "C:\tmp\DCM".  filenumber is the starting filenumber
%   to append, with results such as "C:\tmp\DCM00000", "C:\tmp\DCM00001"
%   etc.
%
% APA, 7/12/2019
%
%Usage:
%   nWritten = export_GSPS_IOD(planC, filenameRoot, filenumber);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS  = planC{end};

%Init number of files written to zero.
nWritten = 0;

    scanNum = 1;
    if length(planC{indexS.scan}) > 1
        destDirPath = fullfile(filenameRoot,['Scan_',num2str(scanNum)]);
        if ~exist(destDirPath,'dir')
            mkdir(destDirPath)
        end
    else
        destDirPath = filenameRoot;
    end

    % Get the associated scanNum
    %scanNum = getStructureAssociatedScan(structNum,planC);

    %Get the scan we are currently working with
    scanS   = planC{indexS.scan}(scanNum);

    if isempty(planC{indexS.GSPS})
        return;
    end

    gspsObjS = planC{indexS.GSPS}(1);

    %Build patient module.
    patientattr      = export_module('patient', 'gsps', gspsObjS);

    %Build general study module.
    studyattr        = export_module('general_study', gspsObjS);

    %Build general series module.
    %seriesattr       = export_module('general_series', gspsObjS);
    seriesattr       = export_module('rt_series', 'gsps', gspsObjS);

    %Build a frame of reference module.
    frameattr        = export_module('frame_of_reference', gspsObjS);

    %Build a general equipment module.
    equipattr        = export_module('general_equipment', 'gsps', gspsObjS);

    %Build a content module
    contentattr        = export_module('content_identification', gspsObjS);

    %Build Graphic Layer module
    %graphicLayerAttr        = export_module('graphic_layer');

    %For slice-specific modules iterate over scaninfo.
    %for i=1:length(planC{indexS.GSPS}) % in future, loop over GSPS objects
    %belonging to a particular scanNum. Currently, we assume all GSPS objects belong
    %only to a single scan.

        gspsS = planC{indexS.GSPS};

        %Create a attr to hold a single slice.
        %CHANGE to dcm4che 3attribute
        attr = javaObject('org.dcm4che3.data.Attributes');

        %Get info for the slice we are handling.
        %scanInfoS = scanS.scanInfo(slcNum);
        %contourS = planC{indexS.structures}(structNum).contour(slc);

        % Build GSPS module from this structure
        gspsattr    =  export_module('gsps', gspsS, scanS);

        SOPattr     = export_module('SOP_common', 'gsps', gspsObjS);

        %Combine all modules into a single attribute.
        attr.addAll(patientattr);
        attr.addAll(studyattr);
        attr.addAll(seriesattr);
        attr.addAll(frameattr);
        attr.addAll(equipattr);
        attr.addAll(contentattr);
        attr.addAll(gspsattr);
        attr.addAll(SOPattr);

        %clear imgattr imgplaneattr imgpixelattr ctimageattr SOPattr

        fileNum  = num2str(filenumber + nWritten);
        filename = fullfile(destDirPath,['GSPS_', repmat('0', [1 5-length(fileNum)]), fileNum]);

        writefile_mldcm(attr, filename);

        nWritten = nWritten + 1;

        clear attr;
    %end

    clear patientattr studyattr seriesattr frameattr equipattr

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR