text:
function computeAllPredictionRF_batchHN(pathExperiments,expNumber,maxOrder,nBoot,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% function computeAllPrediction_HN(pathWORK,fSetName,outcomes,freedomMat,maxOrder,nBoot,imbalance,nBatch,matlabPATH)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes prediction performance estimation for a given
% feature set type, and for all model orders of all experiments with
% different degrees of freedom. See ref. [1,2] for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathExperiments: Full path to where all experiments need to be
%                     performed.
%                     --> Ex: /myProject/WORKSPACE/RESULTS
% 2. expNumber: Number of the single experiment to perform
%                  --> Ex: 4
% 3. maxOrder: Integer specifying the maximal model order to construct.
%              --> Ex: 10
% 4. nBoot: Number of bootstrap samples to use.
%           --> Ex: 100
% 5. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (see ref.[2]).
%               --> Ex: 'IALR'
% 6. nBatch: Number of parallel batch.
%            --> Ex: 8
% 7. matlabPATH: Full path to the MATLAB excutable on the system.
%                --> 'matlab' if a symbolic link to the matlab executable
%                     was previously created.
% -------------------------------------------------------------------------
% OUTPUTS: Prediction performance results are saved in a folder named 'RESULTS' in the
% corresponding folder of each experiment (e.g. 'Experiment1', 'Experiment2', etc.)
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;

% INITIALIZATON
time = 60; % Number of seconds to wait before checking if parallel computations are done
cd(fullfile(pathExperiments,['Experiment',num2str(expNumber)])), load('training')
pathModels = fullfile(pwd,'MODELS'); mkdir('RESULTS'), cd('RESULTS'), pathResults = pwd;
mkdir('batchLog_Results'), cd('batchLog_Results'), pathBatch = pwd;
nameOutcomes = fieldnames(training); nOutcomes = numel(nameOutcomes);
for o = 1:nOutcomes
    outcomes.(nameOutcomes{o}) = training.(nameOutcomes{o}).outcome;
end
setNames = fieldnames(training.(nameOutcomes{1}).variables);
[param] = batchExperiments(setNames,outcomes,nBatch); nBatch = length(param);

% PRODUCE BATCH COMPUTATIONS
save('workspace','pathModels','pathResults','training','param','maxOrder','nBoot','imbalance','seed'), pause(5);
for i = 1:nBatch
    nameScript = ['batch',num2str(i),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    for j = 1:numel(param{i})
        fprintf(fid,['computeAllPredictionRF_HN(pathModels,pathResults,training,param{',num2str(i),'}{',num2str(j),'},maxOrder,nBoot,imbalance,seed)\n']);
    end
    fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nBatch)
delete('workspace.mat')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function cellInfo = AirSpeed_Bus(varargin)
% AIRSPEED_BUS returns a cell array containing bus object information
%
% Optional Input: 'false' will suppress a call to Simulink.Bus.cellToObject
%                 when the MATLAB file is executed.
% The order of bus element attributes is as follows:
%   ElementName, Dimensions, DataType, SampleTime, Complexity, SamplingMode, DimensionsMode, Min, Max, DocUnits, Description

suppressObject = false;
if nargin == 1 && islogical(varargin{1}) && varargin{1} == false
    suppressObject = true;
elseif nargin > 1
    error('Invalid input argument(s) encountered');
end

cellInfo = { ...
  { ...
    'AirSpeed_Bus', ...
    '', ...
    '', ...
    'Auto', ...
    '-1', {...
{'timestamp', 1, 'uint32', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('ms'), ''}; ...
{'diff_pressure', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('Pa'), ''}; ...
{'temperature', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('degC'), ''}; ...
    } ...
  } ...
}';

if ~suppressObject
    % Create bus objects in the MATLAB base workspace
    Simulink.Bus.cellToObject(cellInfo)
end

author: Firmament-Autopilot

repo: FMT-Model


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

%%
clear; close all; clc;

%% Syntax
% |opt=fmin_fvalfind_ppform(pp,x,y);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
report_this_filefun(mfilename('fullpath'));

b = a;

% Hier vielleicht dx, dy aendern ...
ve = [];dx = 0.1; dy = 0.1;

x0 = min(a.Longitude);
x1 = max(a.Longitude);
y0 = min(a.Latitude);
y1 = max(a.Latitude);

for x = x0:dx:x1
    for y = y0:dy:y1
        ve =    [ ve ; x y ];
    end
end

le = length(ve);
Y0 = zeros(le,1);

% Hier muss die richtige abminderung rein ...
for i = 1:length(b)
    di2 = deg2km((distance(ve(:,2),ve(:,1),repmat(b(i,2),le,1),repmat(b(i,1),le,1))));
    R = di2;
    r = sqrt(R.^2 + 5.57^2);
    M = b(i,6); % wenn mit felhler: + dM  (random mit standart)
    Y = -0.136 + 0.229*(M-6) - 0.778 * log10(r) ;
    Y = 10.^Y;
    c = [Y , Y0];
    mapga = max(c');
    Y0 = mapga';
end

mapga = mapga';

l1 = length(x0:dx:x1);
l2 = length(y0:dy:y1);

re = reshape(mapga,l2,l1);
rey = reshape(ve(:,2),l2,l1);
rex = reshape(ve(:,1),l2,l1);

figure
pcolor(rex,rey,re)
hold on
shading interp
overlay
colorbar

author: CelsoReyes


text:
function featuresS = ngldmToScalarFeatures(s,numVoxels)
% function featuresS = ngldmToScalarFeatures(s,numVoxels)
%
% APA, 03/16/2017

% Coarseness
Ns = sum(s(:));
Nn = size(s,2);
Ng = size(s,1);
lenV = 1:Nn;
levV = 1:Ng;

% Low dependence emphasis
sLdeM = bsxfun(@rdivide,s,lenV.^2);
featuresS.LowDependenceEmphasis = sum(sLdeM(:))/Ns;

% High dependence emphasis
sHdeM = bsxfun(@times,s,lenV.^2);
featuresS.HighDependenceEmphasis = sum(sHdeM(:))/Ns;

% Low grey level count emphasis
sLgceM = bsxfun(@rdivide,s',(1:Ng).^2);
featuresS.LowGrayLevelCountEmphasis = sum(sLgceM(:))/Ns;

% High grey level count emphasis
sHgceM = bsxfun(@times,s',(1:Ng).^2);
featuresS.HighGrayLevelCountEmphasis = sum(sHgceM(:))/Ns;

% Low dependence low grey level emphasis
sLdlgeM = bsxfun(@rdivide,bsxfun(@rdivide,s,lenV.^2)',(1:Ng).^2);
featuresS.LowDependenceLowGrayLevelEmphasis = sum(sLdlgeM(:))/Ns;

% Low dependence high grey level emphasis
sLdhgeM = bsxfun(@times,bsxfun(@rdivide,s,lenV.^2)',(1:Ng).^2);
featuresS.LowDependenceHighGrayLevelEmphasis = sum(sLdhgeM(:))/Ns;

% High dependence low grey level emphasis
sHdlgeM = bsxfun(@rdivide,bsxfun(@times,s,lenV.^2)',(1:Ng).^2);
featuresS.HighDependenceLowGrayLevelEmphasis = sum(sHdlgeM(:))/Ns;

% High dependence high grey level emphasis
sHdhgeM = bsxfun(@times,bsxfun(@times,s,lenV.^2)',(1:Ng).^2);
featuresS.HighDependenceHighGrayLevelEmphasis = sum(sHdhgeM(:))/Ns;

% Grey level non-uniformity
featuresS.GrayLevelNonuniformity = sum(sum(s,2).^2)/Ns;

% Grey level non-uniformity normalised
featuresS.GrayLevelNonuniformityNorm = sum(sum(s,2).^2)/Ns^2;

% Dependence count non-uniformity
featuresS.DependenceCountNonuniformity = sum(sum(s,1).^2)/Ns;

% Dependence count non-uniformity normalised
featuresS.DependenceCountNonuniformityNorm = sum(sum(s,1).^2)/Ns^2;

% Dependence count percentage
featuresS.DependenceCountPercentage = Ns/numVoxels;

% Grey level variance
iPij = bsxfun(@times,s'/sum(s(:)),levV);
mu = sum(iPij(:));
iMinusMuPij = bsxfun(@times,s'/sum(s(:)),(levV-mu).^2);
featuresS.GrayLevelVariance = sum(iMinusMuPij(:));

% Dependence count variance
jPij = bsxfun(@times,s/sum(s(:)),lenV);
mu = sum(jPij(:));
jMinusMuPij = bsxfun(@times,s/sum(s(:)),(lenV-mu).^2);
featuresS.DependenceCountVariance = sum(jMinusMuPij(:));

% Dependence count entropy
p = s(:)/sum(s(:));
featuresS.Entropy = -sum(p .* log2(p+eps));

% Dependence count energy
p = s(:)/sum(s(:));
featuresS.Energy = sum(p .^2);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function result = fit_ML_log_normal( x,hAx )
% fit_ML_normal - Maximum Likelihood fit of the log-normal distribution of i.i.d. samples!.
%                  Given the samples of a log-normal distribution, the PDF parameter is found
%
%    fits data to the probability of the form:
%        p(x) = sqrt(1/(2*pi))/(s*x)*exp(- (log(x-m)^2)/(2*s^2))
%    with parameters: m,s
%
% format:   result = fit_ML_log_normal( x,hAx )
%
% input:    x   - vector, samples with log-normal distribution to be parameterized
%           hAx - handle of an axis, on which the fitted distribution is plotted
%                 if h is given empty, a figure is created.
%
% output:   result  - structure with the fields
%                      m,s          - fitted parameters
%                      CRB_m,CRB_s  - Cram?r-Rao Bound for the estimator value
%                      RMS          - RMS error of the estimation
%                      type         - 'ML'
%

%
% Algorithm
% ===========
%
% We use the ML algorithm to estimate the PDF from the samples.
% The log-normal destribution is given by:
%
%    p(x;m,s) = sqrt(1/(2*pi))/(s*x)*exp(- (log(x-m)^2)/(2*s^2))
%
%    where x are the samples which distribute by the function p(x;m,s)
%            and are assumed to be i.i.d !!!
%
% The ML estimator is given by:
%
%    a         = parameters vector = [m,s]
%    f(Xn,a)   = sqrt(1/(2*pi))/(s*Xn)*exp(- (log(Xn-m)^2)/(2*s^2))
%    L(a)      = f(X,a) = product_by_n( f(Xn,a) )
%              = (2*pi*s^2)^(-N/2) * PI{1/Xn} * exp(-sum((log(Xn)-m)^2)/(2*s^2))
%    log(L(a)) = -N/2*log(2*pi) - N*log(s) - sum(log(Xn)) - sum((log(Xn)-m)^2)/(2*s^2)
%
%    The maximum likelihood point is found by the derivative of log(L(a)) with respect to "a":
%
%    diff(log(L(a)),s^2) = N/(2*s^4) * ( sum( (log(Xn)-m)^2 )/N - s^2 )
%                        = J(s^2) * (s^2_estimation - s^2)
%    diff(log(L(a)),m)   = 2*N/(s^2) * ( sum(log(Xn))/N - m )
%                        = J(m) * (m_estimator - m)
%
%    Therefore, the (efficient) estimators are given by:
%
%               m   = sum( log( Xn ) )/N
%               s^2 = sum(( log(Xn)-m)^2)/N
%
%    The Cram?r-Rao Bounds for these estimator are:
%
%               VAR( m )   = 1/J(m)   = (s^2) / N
%               VAR( s^2 ) = 1/J(s^2) = (2*s^4) / N
%
%    NOTE: the ML estimator does not detect a deviation from the model.
%          therefore, check the RMS value !
%

if (nargin<1)
    error( 'fit_ML_log_normal - insufficient input arguments' );
end

% Estimation
% =============
lnx     = log(x(:));                 % should be column vectors !
N       = length(x);
m       = sum( lnx )/N;
s2      = (lnx-m)'*(lnx-m)/N;
CRB_m   = s2 / N;
CRB_s2  = (2*s2^2) / N;
[n,x_c] = hist( x,100 );
n       = n / sum(n*abs(x_c(2)-x_c(1)));
y       = sqrt(1/(2*pi))./(sqrt(s2)*x_c).*exp(- ((log(x_c)-m).^2)/(2*s2));
RMS     = sqrt( (y-n)*((y-n)')/ (x_c(2)-x_c(1))^2 / (length(x_c)-1) );

% finish summarizing results
% ============================
result = struct( 'm',m,'s2',s2,'CRB_m',CRB_m,'CRB_s2',CRB_s2,'RMS',RMS,'type','ML' );

% plot distribution if asked for
% ===============================
if (nargin>1)
    xspan = linspace(min(x),max(x),100);
    if ishandle( hAx )
        plot_log_normal( xspan,result,hAx,1 );
    else
        figure;
        plot_log_normal( xspan,result,gca,1 );
    end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef OmoriModel < int32
    enumeration
        pck   (1) % 3 free parameters: p, c , k      % Modified Omori law (pck)
        pckk  (2) % 4 free parameters: p, c , k1, k2 % MOL with secondary aftershock (pckk)
        ppckk (3) % 5 free parameters: p1,p2,c,k1,k2 % MOL with secondary aftershock
        ppcckk(4) % 6 free parameters: p1,p2,c1, c2,k1,k2 % MOL with secondary aftershock
    end

    methods(Static)
        function cm = doForecast(nMod, t, p_1, c_1, k_1, t_break, k_2, p_2, c_2) % notice funky order!
            % log likelyhood
            %   where p = pvalues as n x 1
            %   where c = cvalues as n x 1
            %   where k = kvalues as n x 1
            %   where t = amount of time after the main shock, as m x 1 duration
            %
            %   returns an n x m matrix of values
            %
            % each row represents ONE time.
            % each column represents one itteration of p,c, and k
            %
            %  for scalar p,c,k, this would return a single row of values. one for each time.
            %  for scalar t, this would returna  single column of values.
            %
            %  p1t1   p2t1  p3t1     params ->
            %  p1t2   p2t2  p3t2
            %  p1t3   p2t3  p3t3
            %  p1t4   p2t4  p3t4
            %
            %  time
            %    |
            %    v

            assert(isrow(p_1) && isrow(c_1) && isrow(k_1) , 'All p,c,k input values should be rows');
            assert(all(numel(p_1) == [ numel(k_1), numel(c_1)]),'p, c, and k values should have same length');
            assert(isvector(t), 'time should be a vector');
            if isduration(t)
                t=days(t);
            end
            if ~iscolumn(t)
                t=t';
            end

            cm = nan(numel(t),numel(p_1));
            sz=size(cm);
            switch nMod
                case OmoriModel.pck
                    % deal with non infinite solutions first
                    idx = p_1 ~= 1;
                    c=c_1(idx);
                    k=k_1(idx);
                    p=p_1(idx);
                    cm(:,idx) = k ./ (p-1) .* (c .^ (1-p)-(t +c).^ (1-p));

                    c=c_1(~idx);
                    k=k_1(~idx);
                    p=[];
                    cm(:,~idx) = k .* log(t ./ c + 1);

                case {OmoriModel.pckk, OmoriModel.ppckk, OmoriModel.ppcckk}

                    assert(isrow(p_2) && isrow(c_2) && isrow(k_2) , 'All p2,c2,k2 input values should be rows');
                    assert(all(numel(p_1) == [ numel(p_2) numel(k_2), numel(c_2)]),'pn, cn, and kn values should have same length');
                    assert(isscalar(t_break));
                    if isduration(t_break)
                        t_break=days(t_break);
                    end

                    isafter = t >= t_break;
                    cm(~isafter,:) =  OmoriModel.doForecast(OmoriModel.pck, t(~isafter), p_1, c_1, k_1);

                    idx = p_1 ~= 1 & p_2 ~= 1;
                    if any(idx)
                        c1 = c_1(idx);
                        c2 = c_2(idx);
                        p1 = p_1(idx);
                        p2 = p_2(idx);
                        k1 = k_1(idx);
                        k2 = k_2(idx);
                        cm(isafter,idx) = k1./(p1-1) .* (c1.^(1-p1)-(t(isafter)+c1).^(1-p1)) + k2./(p2-1).*(c2.^(1-p2)-(t(isafter)-t_break + c2).^(1-p2));
                    end

                    idx = p_1 ~= 1 & p_2 == 1;
                    if any(idx)
                        c1 = c_1(idx);
                        c2 = c_2(idx);
                        p1 = p_1(idx);
                        p2 = [];
                        k1 = k_1(idx);
                        k2 = k_2(idx);
                        cm(isafter,idx) = k1./(p1-1) .* (c1.^(1-p1)-(t(isafter)+c1).^(1-p1)) + k2 .* log((t(isafter)-t_break)./c2+1);
                    end

                    idx = p_1 == 1 & p_2 ~= 1;
                    if any(idx)
                        c1 = c_1(idx);
                        c2 = c_2(idx);
                        p1 = [];
                        p2 = p_2(idx);
                        k1 = k_1(idx);
                        k2 = k_2(idx);
                        cm(isafter,idx) =k1.*log(t(isafter)./c1+1) + k2./(p2-1).*(c2.^(1-p2)-(t(isafter)-t_break+c2).^(1-p2));
                    end

                    idx = p_1 == 1 & p_2 == 1;
                    if any(idx)
                        c1 = c_1(idx);
                        c2 = c_2(idx);
                        p1 = [];
                        p2 = [];
                        k1 = k_1(idx);
                        k2 = k_2(idx);
                        cm(isafter,idx) = k1.*log(t(isafter)./c1+1) + k2.*log((t(isafter)-t_break)./c2+1);
                    end

            end
            assert(isequal(size(cm),sz)); % make sure sizes are as expected
        end
    end
end


text:
function model = addMissingReactions(SampledModel, completeModel)
% Verify the consistency between the sampled model and the
% complete reference model by checking that all reactions from complete
% model are in the sampled model
%
% USAGE:
%
%    model = addMissingReactions(SampledModel, completeModel)
%
% INPUTS:
%    SampledModel:     Sampled model
%    completeModel:    The complete reference model
%
% OUTPUTS:
%    model:            Consistent sampled model wrt the complete model
%
% .. Authors:
%       - Nathan E. Lewis, May 2010-May 2011
%       - Anne Richelle, May 2017

model = completeModel;

% find all nonmissing reactions and make a new points matrix
NonMissingRxns = completeModel.rxns(ismember(completeModel.rxns,SampledModel.rxns));
model.points = zeros(length(model.rxns),length(SampledModel.points(1,:)));

% Move all sampled model points to the new model, and leave missing reaction points as zeros
SM_Ind = findRxnIDs(SampledModel,NonMissingRxns);
CM_Ind = findRxnIDs(model,NonMissingRxns);
model.points(CM_Ind,:) = SampledModel.points(SM_Ind,:);
end

author: opencobra

repo: cobratoolbox


text:
% LEARN - Increments the specified example into the current SVM solution.
%         Assumes alpha_c = 0 initially.
%
% Syntax: nstatus = learn(indc,rflag)
%
% nstatus: new status for indc
%    indc: index of the example to learn
%   rflag: flag indicating whether or not to check if any reserve vectors
%          become margin vectors during learning
%
% Version 3.22e -- Comments to diehl@alumni.cmu.edu
%

function nstatus = learn(indc,rflag)

% flags for example state
MARGIN    = 1;
ERROR     = 2;
RESERVE   = 3;
UNLEARNED = 4;

% define global variables
global a;             % alpha coefficients
global b;             % bias
global C;             % regularization parameters
global deps;          % jitter factor in kernel matrix
global g;             % partial derivatives of cost function w.r.t. alpha coefficients
global ind;           % structure containing indices of margin, error, reserve and unlearned vectors
global perturbations; % number of perturbations
global Q;             % extended kernel matrix for all vectors
global Rs;            % inverse of extended kernel matrix for margin vectors
global scale;         % kernel scale
global type;          % kernel type
global X;             % matrix of margin, error, reserve and unlearned vectors stored columnwise
global y;             % column vector of class labels (-1/+1) for margin, error, reserve and unlearned vectors

% compute g(indc)
[f_c,K] = svmeval(X(:,indc));
g(indc) = y(indc)*f_c - 1;

% if g(indc) > 0, place this example into the reserve set directly
if (g(indc) >= 0)

   % move the example to the reserve set
   bookkeeping(indc,UNLEARNED,RESERVE);
   nstatus = RESERVE;

   return;
end;

% compute Qcc and Qc if necessary
num_MVs = length(ind{MARGIN});
Qc = cell(3,1);
if (num_MVs == 0)
	if (length(ind{ERROR}) > 0)
   	Qc{ERROR} = (y(ind{ERROR})*y(indc)).*kernel(X(:,ind{ERROR}),X(:,indc),type,scale);
   end;
else
	Qc{MARGIN} = (y(ind{MARGIN})*y(indc)).*K(1:num_MVs);
	if (length(ind{ERROR}) > 0)
   	Qc{ERROR} = (y(ind{ERROR})*y(indc)).*K(num_MVs+1:length(K));
	end;
end;
if (length(ind{RESERVE}) > 0)
   Qc{RESERVE} = (y(ind{RESERVE})*y(indc)).*kernel(X(:,ind{RESERVE}),X(:,indc),type,scale);
end;
Qcc = kernel(X(:,indc),X(:,indc),type,scale) + deps;

converged = 0;
while (~converged)

   perturbations = perturbations + 1;

   if (num_MVs > 0)  % change in alpha_c permitted

      % compute Qc, beta and gamma
      beta = -Rs*[y(indc) ; Qc{MARGIN}];
      gamma = zeros(size(Q,2),1);
      ind_temp = [ind{ERROR} ind{RESERVE} indc];
      gamma(ind_temp) = [Qc{ERROR} ; Qc{RESERVE} ; Qcc] + Q(:,ind_temp)'*beta;

      % check if gamma_c < 0 (kernel matrix is not positive semi-definite)
      if (gamma(indc) < 0)
         error('LEARN: gamma_c < 0');
      end;

   else  % change in alpha_c not permitted since the constraint on the sum of the
         % alphas must be preserved.  only b can change.

      % set beta and gamma
      beta = y(indc);
      gamma = y(indc)*y;

   end;

   % minimum acceptable parameter change (change in alpha_c (num_MVs > 0) or b (num_MVs = 0))
   [min_delta_param,indss,cstatus,nstatus] = min_delta_acb(indc,gamma,beta,1,rflag);

   % update a, b, and g
   if (num_MVs > 0)
      a(indc) = a(indc) + min_delta_param;
      a(ind{MARGIN}) = a(ind{MARGIN}) + beta(2:num_MVs+1)*min_delta_param;
   end;
   b = b + beta(1)*min_delta_param;
   g = g + gamma*min_delta_param;

   % update Qc and perform bookkeeping
   converged = (indss == indc);
   if (converged)
      cstatus = UNLEARNED;
     	Qc{nstatus} = [Qc{nstatus} ; Qcc];
  	else
  		ind_temp = find(ind{cstatus} == indss);
  		Qc{nstatus} = [Qc{nstatus} ; Qc{cstatus}(ind_temp)];
  		Qc{cstatus}(ind_temp) = [];
   end;
   [indco,removed_i] = bookkeeping(indss,cstatus,nstatus);
   if ((nstatus == RESERVE) & (removed_i > 0))
      Qc{nstatus}(removed_i) = [];
   end;

   % set g(ind{MARGIN}) to zero
   g(ind{MARGIN}) = 0;

   % update Rs and Q if necessary
   if (nstatus == MARGIN)

      num_MVs = num_MVs + 1;
      if (num_MVs > 1)
         if (converged)
            gamma = gamma(indss);
         else

            % compute beta and gamma for indss
            beta = -Rs*Q(:,indss);
            gamma = kernel(X(:,indss),X(:,indss),type,scale) + deps + Q(:,indss)'*beta;

         end;
      end;

      % expand Rs and Q
      updateRQ(beta,gamma,indss);

   elseif (cstatus == MARGIN)

      % compress Rs and Q
      num_MVs = num_MVs - 1;
      updateRQ(indco);

   end;

end;

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ML,MU] = munu_m_d
%
%  Deletes the global data that has been generated by 'munu_m_i'.
%
%  Calling sequence:
%
%    munu_m_d
%    [ML,MU] = munu_m_d
%
%  Output:
%
%    ML, MU      LU factors of M
%
%  ML and MU are returned just in case the LU factors of M are needed for some
%  reason after the computation is done, but have been deleted to save memory.
%
%
%  LYAPACK 1.0 (Thilo Penzl, September 1999)

global LP_ML LP_MU LP_N LP_M

if nargout==1
  ML = LP_ML;
  MU = LP_MU;
end

clear global LP_ML LP_MU LP_N LP_M

author: Sable

repo: mcbench-benchmarks


text:
function Archive = UpdateArchive(Archive,Population,popsize)
% Update the archive in MSOPS-II

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

    %% Combine the archive with the population
    Archive = [Archive,Population];
    ArchObj = Archive.objs;
    N       = length(Archive);

    %% Update the archive by weighted min-max metric
    % Calculate the weighted min-max metric between each two solutions
    WMM = CalMetric(ArchObj,ArchObj);
    WMM_diagonal = WMM(logical(eye(N)))';
    % Delete the solutions which do not have the lowest metric value than
    % others according to its own weight vector
    Remain = true(1,N);
    for i = 1 : N
        if Remain(i)
            if WMM(i,i) > min(WMM(Remain,i))
                Remain(i) = false;
            else
                Remain(WMM(i,:)<WMM_diagonal) = false;
            end
        end
    end
    Archive = Archive(Remain);
    % If the archive contains too many solutions, randomly delete some. The
    % original algorithm does not limit the size of archive, so that the
    % size of archive will increase unrestrainedly
    if length(Archive) > 10*popsize
        Archive = Archive(randperm(length(Archive),5*popsize));
    end
end


text:
function gpcf = gpcf_prod(varargin)
%GPCF_PROD  Create a product form covariance function
%
%  Description
%    GPCF = GPCF_PROD('cf', {GPCF_1, GPCF_2, ...})
%    creates a product form covariance function
%          GPCF = GPCF_1 .* GPCF_2 .* ... .* GPCF_N
%
%  See also
%    GP_SET, GPCF_*
%
% Copyright (c) 2009-2010, 2018 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari
% Copyright (c) 2014 Arno Solin and Jukka Koskenranta

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

  ip=inputParser;
  ip.FunctionName = 'GPCF_PROD';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('cf',[], @iscell);
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

  if isempty(gpcf)
    init=true;
    gpcf.type = 'gpcf_prod';
  else
    if ~isfield(gpcf,'type') && ~isequal(gpcf.type,'gpcf_prod')
      error('First argument does not seem to be a valid covariance function structure')
    end
    init=false;
  end

  if init || ~ismember('cf',ip.UsingDefaults)
    % Initialize parameters
    gpcf.cf = {};
    cfs=ip.Results.cf;
    if ~isempty(cfs)
      for i = 1:length(cfs)
        gpcf.cf{i} = cfs{i};
      end
    else
      error('At least one covariance function has to be given in cf');
    end
  end

  if init
    % Set the function handles to the subfunctions
    gpcf.fh.pak = @gpcf_prod_pak;
    gpcf.fh.unpak = @gpcf_prod_unpak;
    gpcf.fh.lp = @gpcf_prod_lp;
    gpcf.fh.lpg = @gpcf_prod_lpg;
    gpcf.fh.cfg = @gpcf_prod_cfg;
    gpcf.fh.cfdg = @gpcf_prod_cfdg;
    gpcf.fh.cfdg2 = @gpcf_prod_cfdg2;
    gpcf.fh.ginput = @gpcf_prod_ginput;
    gpcf.fh.ginput2 = @gpcf_prod_ginput2;
    gpcf.fh.ginput3 = @gpcf_prod_ginput3;
    gpcf.fh.ginput4 = @gpcf_prod_ginput4;
    gpcf.fh.cov = @gpcf_prod_cov;
    gpcf.fh.trcov  = @gpcf_prod_trcov;
    gpcf.fh.trvar  = @gpcf_prod_trvar;
    gpcf.fh.recappend = @gpcf_prod_recappend;
    gpcf.fh.cf2ss = @gpcf_prod_cf2ss;
  end

end

function [w, s, h] = gpcf_prod_pak(gpcf)
%GPCF_PROD_PAK  Combine GP covariance function parameters into one vector
%
%  Description
%    W = GPCF_PROD_PAK(GPCF, W) loops through all the covariance
%    functions and packs their parameters into one vector as
%    described in the respective functions. This is a mandatory
%    subfunction used for example in energy and gradient computations.
%
%  See also
%    GPCF_PROD_UNPAK

  ncf = length(gpcf.cf);
  w = []; s = {}; h=[];

  for i=1:ncf
    cf = gpcf.cf{i};
    [wi, si, hi] = cf.fh.pak(cf);
    w = [w wi];
    s = [s; si];
    h = [h hi];
  end
end

function [gpcf, w] = gpcf_prod_unpak(gpcf, w)
%GPCF_PROD_UNPAK  Sets the covariance function parameters into
%                 the structures
%
%  Description
%    [GPCF, W] = GPCF_PROD_UNPAK(GPCF, W) loops through all the
%    covariance functions and unpacks their parameters from W to
%    each covariance function structure. This is a mandatory
%    subfunction used for example in energy and gradient computations.
%
%  See also
%    GPCF_PROD_PAK
%
  ncf = length(gpcf.cf);

  for i=1:ncf
    cf = gpcf.cf{i};
    [cf, w] = cf.fh.unpak(cf, w);
    gpcf.cf{i} = cf;
  end

end

function lp = gpcf_prod_lp(gpcf)
%GPCF_PROD_LP  Evaluate the log prior of covariance function parameters
%
%  Description
%    LP = GPCF_PROD_LP(GPCF, X, T) takes a covariance function
%    structure GPCF and returns log(p(th)), where th collects the
%    parameters. This is a mandatory subfunction used for example
%    in energy computations.
%
%  See also
%    GPCF_PROD_PAK, GPCF_PROD_UNPAK, GPCF_PROD_LPG, GP_E

  lp = 0;
  ncf = length(gpcf.cf);
  for i=1:ncf
    cf = gpcf.cf{i};
    lp = lp + cf.fh.lp(cf);
  end

end

function lpg = gpcf_prod_lpg(gpcf)
%GPCF_PROD_LPG  Evaluate gradient of the log prior with respect
%               to the parameters.
%
%  Description
%    LPG = GPCF_PROD_LPG(GPCF) takes a covariance function
%    structure GPCF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters. This is a mandatory subfunction
%    used for example in gradient computations.
%
%  See also
%    GPCF_PROD_PAK, GPCF_PROD_UNPAK, GPCF_PROD_LP, GP_G
  lpg = [];
  ncf = length(gpcf.cf);

  % Evaluate the gradients
  for i=1:ncf
    cf = gpcf.cf{i};
    lpg=[lpg cf.fh.lpg(cf)];
  end

end

function DKff = gpcf_prod_cfg(gpcf, x, x2, mask, i1)
%GPCF_PROD_CFG  Evaluate gradient of covariance function
%               with respect to the parameters.
%
%  Description
%    DKff = GPCF_PROD_CFG(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to th (cell array with matrix elements). This is a
%    mandatory subfunction used in gradient computations.
%
%    DKff = GPCF_PROD_CFG(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_PROD_CFG(GPCF, X, [], MASK) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the diagonal of gradients of covariance matrix
%    Kff = k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_PROD_CFG(GPCF, X, X2, [], i) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith
%    hyperparameter. This subfunction is needed when using
%    memory save option in gp_set.
%
%  See also
%    GPCF_PROD_PAK, GPCF_PROD_UNPAK, GPCF_PROD_LP, GP_G

  [n, m] =size(x);
  ncf = length(gpcf.cf);

  DKff = {};

  if nargin==5
    % Use memory save option
    savememory=1;
    i3=0;
    for k=1:ncf
      % Number of hyperparameters for each covariance function
      cf=gpcf.cf{k};
      i3(k)=cf.fh.cfg(cf,[],[],[],0);
    end
    if i1==0
      % Return number of hyperparameters
      DKff=sum(i3);
      return
    end
    % Help indices
    i3=cumsum(i3);
    ind=find(cumsum(i3 >= i1)==1);
    if ind>1
      i1=[ind i1-i3(ind-1)];
    else
      i1=[ind i1];
    end
  else
    savememory=0;
  end

  % Evaluate: DKff{1} = d Kff / d magnSigma2
  %           DKff{2} = d Kff / d lengthScale
  % NOTE! Here we have already taken into account that the parameters are transformed
  % through log() and thus dK/dlog(p) = p * dK/dp

  % evaluate the gradient for training covariance
  if nargin == 2 || (isempty(x2) && isempty(mask))

    % evaluate the individual covariance functions
    for i=1:ncf
      cf = gpcf.cf{i};
      C{i} = cf.fh.trcov(cf, x);
    end

    % Evaluate the gradients
    ind = 1:ncf;
    DKff = {};
    if ~savememory
      i3=1:ncf;
    else
      i3=i1(1);
    end
    for i=i3
      cf = gpcf.cf{i};
      if ~savememory
        DK = cf.fh.cfg(cf, x);
      else
        DK = {cf.fh.cfg(cf,x,[],[],i1(2))};
      end

      CC = 1;
      for kk = ind(ind~=i)
        CC = CC.*C{kk};
      end
      for j = 1:length(DK)
        DKff{end+1} = DK{j}.*CC;
      end
    end

    % Evaluate the gradient of non-symmetric covariance (e.g. K_fu)
  elseif nargin == 3 || isempty(mask)
    if size(x,2) ~= size(x2,2)
      error('gpcf_prod -> _ghyper: The number of columns in x and x2 has to be the same. ')
    end

    % evaluate the individual covariance functions
    for i=1:ncf
      cf = gpcf.cf{i};
      C{i} = cf.fh.cov(cf, x, x2);
    end

    % Evaluate the gradients
    ind = 1:ncf;
    DKff = {};
    if ~savememory
      i3=1:ncf;
    else
      i3=i1(1);
    end
    for i=i3
      cf = gpcf.cf{i};
      if ~savememory
        DK = cf.fh.cfg(cf, x,x2);
      else
        DK = {cf.fh.cfg(cf,x,x2,[],i1(2))};
      end

      CC = 1;
      for kk = ind(ind~=i)
        CC = CC.*C{kk};
      end

      for j = 1:length(DK)
        DKff{end+1} = DK{j}.*CC;
      end
    end

    % Evaluate: DKff{1}    = d mask(Kff,I) / d magnSigma2
    %           DKff{2...} = d mask(Kff,I) / d lengthScale
  elseif nargin == 4 || nargin == 5

    % evaluate the individual covariance functions
    for i=1:ncf
      cf = gpcf.cf{i};
      C{i} = cf.fh.trvar(cf, x);
    end

    % Evaluate the gradients
    ind = 1:ncf;
    DKff = {};
    if ~savememory
      i3=1:ncf;
    else
      i3=i1(1);
    end
    for i=i3
      cf = gpcf.cf{i};
      if ~savememory
        DK = cf.fh.cfg(cf, x, [], 1);
      else
        DK = {cf.fh.cfg(cf, x, [], 1, i1(2))};
      end

      CC = 1;
      for kk = ind(ind~=i)
        CC = CC.*C{kk};
      end

      for j = 1:length(DK)
        DKff{end+1} = DK{j}.*CC;
      end
    end
  end
  if savememory
    DKff=DKff{1};
  end
end

function DKff = gpcf_prod_cfdg(gpcf, x, x2, dims)
%GPCF_SEXP_CFDG  Evaluate gradient of covariance function, of
%                which has been taken partial derivative with
%                respect to x, with respect to parameters.
%
%  Description
%    DKff = GPCF_SEXP_CFDG(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of derivatived covariance matrix
%    dK(df,f)/dhyp = d(d k(X,X)/dx)/dhyp, with respect to the
%    parameters
%
%    Evaluate: DKff{1:m} = d Kff / d magnSigma2
%              DKff{m+1:2m} = d Kff / d lengthScale_m
%    m is the dimension of inputs. If ARD is used, then multiple
%    lengthScales. This subfunction is needed when using derivative
%    observations.
%
%         dims - is a vector of input dimensions with respect to which the
%                derivatives of the covariance function have been calculated
%                [by default dims=1:size(x,2)]
%
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%    GPCF_SEXP_GINPUT

ncf = length(gpcf.cf);

[~, m] =size(x);
if nargin <3 || isempty(x2)
    x2=x;
end
if nargin < 4 || isempty(dims)
    %dims1 = 1:m;
    error('dims1 needs to be given')
end

% check whether there are any covariance functions that depend on dims
% covariate
ind = [];
DKff = {};
for i1 = 1:ncf
    if (~isfield(gpcf.cf{i1}, 'selectedVariables') || any(ismember(gpcf.cf{i1}.selectedVariables,dims)))
        ind(end+1) = i1;
    end
end
if ~isempty(ind)
    % Calculate derivative
    % calculate the individual covariance functions
    Ci = {};
    C = 1;
    for i1=1:ncf
        Ci{i1} = gpcf.cf{i1}.fh.cov(gpcf.cf{i1}, x, x2);
        C = C.*Ci{i1};
    end

    if ncf==2
        if any(ind==1)
            DKdx1 = gpcf.cf{1}.fh.cfdg(gpcf.cf{1}, x, x2, dims);
            DKx1 = gpcf.cf{1}.fh.ginput4(gpcf.cf{1}, x, x2, dims);
        else
            DKdx1{1} = 0;
            DKx1{1} = 0;
        end
        DK1 = gpcf.cf{1}.fh.cfg(gpcf.cf{1}, x, x2);
        if any(ind==2)
            DKdx2 = gpcf.cf{2}.fh.cfdg(gpcf.cf{2}, x, x2, dims);
            DKx2 = gpcf.cf{2}.fh.ginput4(gpcf.cf{2}, x, x2, dims);
        else
            DKdx2{1} = 0;
            DKx2{1} = 0;
        end
        DK2 = gpcf.cf{2}.fh.cfg(gpcf.cf{2}, x, x2);

        for i2 = 1:length(DK1)
            DKff{end+1} = DKdx1{i2}.*C./Ci{1} + DK1{i2}.*DKx2{1};
        end
        for i2 = 1:length(DK2)
            DKff{end+1} = DKdx2{i2}.*C./Ci{2} + DK2{i2}.*DKx1{1};
        end

    else
        error('not implemented yet')
        for i1=1:ncf

            DKdx = gpcf.cf{i1}.fh.cfdg(gpcf.cf{i1}, x, x2, dims);
            DK = gpcf.cf{i1}.fh.cfg(gpcf.cf{i1}, x, x2);

            gpcftmp = gpcf;
            gpcftmp.cf = {gpcftmp.cf(1:ncf~=i1)};
            % continue from here
            %DKx2 =
        end
    end
else
    DK = gpcf.fh.lpg(gpcf);
    for i2 = 1:length(DK)
        DKff{end+1} = 0;
    end
end

end

function DKff = gpcf_prod_cfdg2(gpcf, x, x2, dims1, dims2)
%GPCF_SEXP_CFDG  Evaluate gradient of covariance function, of
%                which has been taken partial derivative with
%                respect to x, with respect to parameters.
%
%  Description
%    DKff = GPCF_SEXP_CFDG(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of derivatived covariance matrix
%    dK(df,f)/dhyp = d(d k(X,X)/dx)/dhyp, with respect to the
%    parameters
%
%    Evaluate: DKff{1:m} = d Kff / d magnSigma2
%              DKff{m+1:2m} = d Kff / d lengthScale_m
%    m is the dimension of inputs. If ARD is used, then multiple
%    lengthScales. This subfunction is needed when using derivative
%    observations.
%
%         dims - is a vector of input dimensions with respect to which the
%                derivatives of the covariance function have been calculated
%                [by default dims=1:size(x,2)]
%
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%    GPCF_SEXP_GINPUT

ncf = length(gpcf.cf);

[~, m] =size(x);
if nargin <3 || isempty(x2)
    x2=x;
end
if nargin < 4 || isempty(dims1)
    %dims1 = 1:m;
    error('dims1 needs to be given')
end
if nargin < 5 || isempty(dims2)
    %dims2 = 1:m;
    error('dims2 needs to be given')
end

% NOTICE. AS OF NOW we assume that dims1 and dims2 are scalars

% check whether there are any covariance functions that depend on dims
% covariate
ind = [];
DKff = {};
for i1 = 1:ncf
    if (~isfield(gpcf.cf{i1}, 'selectedVariables') || any(ismember(gpcf.cf{i1}.selectedVariables,dims1)) || any(ismember(gpcf.cf{i1}.selectedVariables,dims2)))
        ind(end+1) = i1;
    end
end
if ~isempty(ind)
    % Calculate derivative
    % calculate the individual covariance functions
    Ci = {};
    C = 1;
    for i1=1:ncf
        Ci{i1} = gpcf.cf{i1}.fh.cov(gpcf.cf{i1}, x, x2);
        %C = C.*Ci{i1};
    end

    if ncf==2
        if any(ind==1)
            DKdx11 = gpcf.cf{1}.fh.cfdg(gpcf.cf{1}, x, x2, dims1);
            DKdx12 = gpcf.cf{1}.fh.cfdg(gpcf.cf{1}, x2, x, dims2);
            DKdx1_2 = gpcf.cf{1}.fh.cfdg2(gpcf.cf{1}, x, x2, dims1, dims2);
            DKx11 = gpcf.cf{1}.fh.ginput4(gpcf.cf{1}, x, x2, dims1);
            DKx12 = gpcf.cf{1}.fh.ginput4(gpcf.cf{1}, x2, x, dims2);
        else
            DKdx11{1} = 0;
            DKdx12{1} = 0;
            DKdx1_2{1} = 0;
            DKx11{1} = 0;
            DKx12{1} = 0;
        end
        if any(ind==2)
            DKdx21 = gpcf.cf{2}.fh.cfdg(gpcf.cf{2}, x, x2, dims1);
            DKdx22 = gpcf.cf{2}.fh.cfdg(gpcf.cf{2}, x2, x, dims2);
            DKdx2_2 = gpcf.cf{2}.fh.cfdg2(gpcf.cf{2}, x, x2, dims1, dims2);
            DKx21 = gpcf.cf{2}.fh.ginput4(gpcf.cf{2}, x, x2, dims1);
            DKx22 = gpcf.cf{2}.fh.ginput4(gpcf.cf{2}, x2, x, dims2);
        else
            DKdx21{1} = 0;
            DKdx22{1} = 0;
            DKdx2_2{1} = 0;
            DKx21{1} = 0;
            DKx22{1} = 0;
        end

        if dims1==dims2
            if any(ind==1)
                D2Kx1 = gpcf.cf{1}.fh.ginput2(gpcf.cf{1}, x, x2, dims1, dims2);
            else
                D2Kx1{1} = 0;
            end
            if any(ind==2)
                D2Kx2 = gpcf.cf{2}.fh.ginput2(gpcf.cf{2}, x, x2, dims1, dims2);
            else
                D2Kx2{1} = 0;
            end
        else
%             D2Kx1 = gpcf.cf{1}.fh.ginput3(gpcf.cf{1}, x, x2, dims1, dims2);
%             D2Kx2 = gpcf.cf{2}.fh.ginput3(gpcf.cf{2}, x, x2, dims1, dims2);
            if any(ind==1)
                D2Kx1 = gpcf.cf{1}.fh.ginput3(gpcf.cf{1}, x, x2, dims1, dims2);
            else
                D2Kx1{1} = 0;
            end
            if any(ind==2)
                D2Kx2 = gpcf.cf{2}.fh.ginput3(gpcf.cf{2}, x, x2, dims1, dims2);
            else
                D2Kx2{1} = 0;
            end
        end

        DK1 = gpcf.cf{1}.fh.cfg(gpcf.cf{1}, x, x2);
        DK2 = gpcf.cf{2}.fh.cfg(gpcf.cf{2}, x, x2);

        for i2 = 1:length(DK1)
            DKff{end+1} = DKdx1_2{i2}.*Ci{2} + DKdx11{i2}.*DKx22{1}' + DKdx12{i2}'.*DKx21{1} + DK1{i2}.*D2Kx2{1};
        end
        for i2 = 1:length(DK2)
            DKff{end+1} = DKdx2_2{i2}.*Ci{1} + DKdx21{i2}.*DKx12{1}' + DKdx22{i2}'.*DKx11{1} + DK2{i2}.*D2Kx1{1};
        end

    else
        error('not implemented yet')
        for i1=1:ncf

            DKdx = gpcf.cf{i1}.fh.cfdg(gpcf.cf{i1}, x, x2, dims);
            DK = gpcf.cf{i1}.fh.cfg(gpcf.cf{i1}, x, x2);

            gpcftmp = gpcf;
            gpcftmp.cf = {gpcftmp.cf(1:ncf~=i1)};
            % continue from here
            %DKx2 =
        end
    end
else
    DK = gpcf.fh.lpg(gpcf);
    for i2 = 1:length(DK)
        DKff{end+1} = 0;
    end
end

% [n, m] =size(x);
% ncf = length(gpcf.cf);
%
% if nargin <3 || isempty(x2)
%     x2=x;
% end
% if nargin < 4 || isempty(dims1)
%     %dims1 = 1:m;
%     error('dims1 needs to be given')
% end
% if nargin < 5 || isempty(dims2)
%     %dims2 = 1:m;
%     error('dims2 needs to be given')
% end
%
%
% % check whether there are any covariance functions that depend on dims
% % covariate
% ind = [];
% DKff = {};
% for i1 = 1:ncf
%     if (~isfield(gpcf.cf{i1}, 'selectedVariables') || any(ismember(gpcf.cf{i1}.selectedVariables,dims1)) || any(ismember(gpcf.cf{i1}.selectedVariables,dims2)))
%         ind(end+1) = i1;
%     end
% end
% if ~isempty(ind)
%     % Calculate derivative
%
%     % calculate the individual covariance functions
%     Ci = {};
%     C = 1;
%     for i1=1:ncf
%         Ci{i1} = gpcf.cf{i1}.fh.cov(gpcf.cf{i1}, x, x2);
%         C = C.*Ci{i1};
%     end
%
%     for i1=1:ncf
%         if any(i1 == ind)
%             DK = gpcf.cf{i1}.fh.cfdg2(gpcf.cf{i1}, x, x2, dims1, dims2);
%             for i2 = 1:length(DK)
%                 DKff{end+1} = DK{i2}.*C./Ci{i1};
%             end
%         else
%             DK = gpcf.cf{i1}.fh.lpg(gpcf.cf{i1});
%             for i2 = 1:length(DK)
%                 DKff{end+1} = 0;
%             end
%         end
%     end
% else
%     DK = gpcf.fh.lpg(gpcf);
%     for i2 = 1:length(DK)
%         DKff{end+1} = 0;
%     end
% end

end

function DKff = gpcf_prod_ginput(gpcf, x, x2, i1)
%GPCF_PROD_GINPUT  Evaluate gradient of covariance function with
%                  respect to x
%
%  Description
%    DKff = GPCF_PROD_GINPUT(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to X (cell array with matrix elements). This subfunction
%    is needed when computing gradients with respect to inducing
%    inputs in sparse approximations.
%
%    DKff = GPCF_PROD_GINPUT(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_PROD_GINPUT(GPCF, X, X2, i) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith
%    covariate in X (cell array with matrix elements). This
%    subfunction is needed when using memory save option in
%    gp_set.
%
%  See also
%    GPCF_PROD_PAK, GPCF_PROD_UNPAK, GPCF_PROD_LP, GP_G

  [n, m] =size(x);

  if nargin==4
    % Use memory save option
    savememory=1;
    if i1==0
      % Return number of covariates
      if isfield(gpcf,'selectedVariables')
        DKff=length(gpcf.selectedVariables);
      else
        DKff=m;
      end
      return
    end
  else
    savememory=0;
  end
  % evaluate the gradient for training covariance
  if nargin == 2 || isempty(x2)

    ncf = length(gpcf.cf);

    % evaluate the individual covariance functions
    for i=1:ncf
      cf = gpcf.cf{i};
      C{i} = cf.fh.trcov(cf, x);
    end

    % Evaluate the gradients
    ind = 1:ncf;
    if ~savememory
      DKff=cellfun(@(a) zeros(n,n), cell(1,numel(x)), 'UniformOutput', 0);
    else
      DKff=cellfun(@(a) zeros(n,n), cell(1,n), 'UniformOutput', 0);
    end
    for i=1:ncf
      cf = gpcf.cf{i};
      if ~savememory
        DK = cf.fh.ginput(cf, x);
      else
        DK = cf.fh.ginput(cf, x, [], i1);
      end

      CC = 1;
      for kk = ind(ind~=i)
        CC = CC.*C{kk};
      end

      for j = 1:length(DK)
        DKff{j} = DKff{j} + DK{j}.*CC;
      end
    end

    % Evaluate the gradient of non-symmetric covariance (e.g. K_fu)
  elseif nargin == 3 ||  nargin == 4
    if size(x,2) ~= size(x2,2)
      error('gpcf_prod -> _ghyper: The number of columns in x and x2 has to be the same. ')
    end

    ncf = length(gpcf.cf);

    % evaluate the individual covariance functions
    for i=1:ncf
      cf = gpcf.cf{i};
      C{i} = cf.fh.cov(cf, x, x2);
    end

    % Evaluate the gradients
    ind = 1:ncf;
    if ~savememory
      DKff=cellfun(@(a) zeros(n,n), cell(1,numel(x)), 'UniformOutput', 0);
    else
      DKff=cellfun(@(a) zeros(n,n), cell(1,n), 'UniformOutput', 0);
    end
    for i=1:ncf
      cf = gpcf.cf{i};
      if ~savememory
        DK = cf.fh.ginput(cf, x, x2);
      else
        DK = cf.fh.ginput(cf, x, x2, i1);
      end

      CC = 1;
      for kk = ind(ind~=i)
        CC = CC.*C{kk};
      end

      for j = 1:length(DK)
        DKff{j} = DKff{j} + DK{j}.*CC;
      end
    end
  end

end

function [DKff, DKff1, DKff2]  = gpcf_prod_ginput2(gpcf, x, x2, dims, takeOnlyDiag)
%GPCF_PROD_GINPUT2  Evaluate gradient of covariance function with
%                   respect to both input variables x and x2 (in
%                   same dimension).
%
%  Description
%    DKff = GPCF_PROD_GINPUT2(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of twice derivatived covariance
%    matrix K(df,df) = dk(X1,X2)/dX1dX2 (cell array with matrix
%    elements). Input variable's dimensions are expected to be
%    same. The function returns also DKff1 and DKff2 which are
%    parts of DKff and needed with CFDG2. DKff = DKff1 -
%    DKff2. This subfunction is needed when using derivative
%    observations.
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%    GPCF_PROD_GINPUT, GPCF_PROD_GINPUT2, GPCF_PROD_CFDG2
if nargin < 5
    takeOnlyDiag=[];
end

ncf = length(gpcf.cf);

% check whether there are any covariance functions that depend on dims
% covariate
ind = [];
for i1 = 1:ncf
    if (~isfield(gpcf.cf{i1}, 'selectedVariables') || any(ismember(gpcf.cf{i1}.selectedVariables,dims)))
        ind(end+1) = i1;
    end
end
% Calculate derivative
if ~isempty(ind)
    % evaluate the individual covariance functions
    C = 1;
    K = 0;
    for i1=1:ncf
        if ~isempty(takeOnlyDiag) && strcmp(takeOnlyDiag, 'takeOnlyDiag')
            Ci = gpcf.cf{i1}.fh.trvar(gpcf.cf{i1}, x);
        else
            Ci = gpcf.cf{i1}.fh.cov(gpcf.cf{i1}, x, x2);
        end
        C = C.*Ci;
        if any(i1 == ind)
            Dt = gpcf.cf{i1}.fh.ginput2(gpcf.cf{i1}, x, x2, dims, takeOnlyDiag);
            if length(Dt)>1
                error('gpcf_prod->ginput2: something is wrong')
            end
            K = K+Dt{1}./Ci;
        end
    end
    DKff = {C.*K};
else
    DKff = {0};
end

end

function DKff = gpcf_prod_ginput3(gpcf, x, x2, dims1, dims2)
%GPCF_PROD_GINPUT3  Evaluate gradient of covariance function with
%                   respect to both input variables x and x2 (in
%                   different dimensions).
%
%  Description
%    DKff = GPCF_PROD_GINPUT3(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of twice derivatived covariance
%    matrix K(df,df) = dk(X1,X2)/dX1dX2 (cell array with matrix
%    elements). The derivative is calculated in multidimensional
%    problem between input's observation dimensions which are not
%    same. This subfunction is needed when using derivative
%    observations.
%
%    ---- !!note this help text needs to be corrected !! ---
%    DKff is a cell array with the following elements:
%      DKff{1} = dk(X1,X2)/dX1_1dX2_2
%      DKff{2} = dk(X1,X2)/dX1_1dX2_3
%       ...
%      DKff{m-1} = dk(X1,X2)/dX1_1dX2_m
%      DKff{m} = dk(X1,X2)/dX1_2dX2_3
%       ...
%      DKff{m} = dk(X1,X2)/dX1_(m-1)dX2_m
%    where _m denotes the input dimension with respect to which the
%    gradient is calculated.
%     ---- clip ---
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%    GPCF_PROD_GINPUT, GPCF_PROD_GINPUT2, GPCF_PROD_CFDG2

ncf = length(gpcf.cf);

% check whether there are any covariance functions that depend on dims
% covariate
ind = [];
for i1 = 1:ncf
    if (~isfield(gpcf.cf{i1}, 'selectedVariables') || any(ismember(gpcf.cf{i1}.selectedVariables,dims1)) || any(ismember(gpcf.cf{i1}.selectedVariables,dims2)) )
        ind(end+1) = i1;
    end
end
% Calculate derivative
if ~isempty(ind)
    % evaluate the individual covariance functions
    C = 1;
    K = 0;
    for i1=1:ncf
        Ci{i1} = gpcf.cf{i1}.fh.cov(gpcf.cf{i1}, x, x2);
        C = C.*Ci{i1};
    end
    for i1=1:ncf
        if any(i1 == ind)
            % double derivative of i'th covariance function
            Dt = gpcf.cf{i1}.fh.ginput3(gpcf.cf{i1}, x, x2, dims1, dims2);
            if length(Dt)>1
                error('gpcf_prod->ginput3: something is wrong')
            end
            K = K+Dt{1}./Ci{i1};
            % pairwise derivatives of i'th and j'th covariance functions
            Dti1 = gpcf.cf{i1}.fh.ginput4(gpcf.cf{i1}, x, x2, dims1);
            Dti2 = gpcf.cf{i1}.fh.ginput4(gpcf.cf{i1}, x2, x, dims2);
            for j1=i1+1:ncf
                if any(j1 == ind)
                    Dtj1 = gpcf.cf{j1}.fh.ginput4(gpcf.cf{j1}, x, x2, dims1);
                    Dtj2 = gpcf.cf{j1}.fh.ginput4(gpcf.cf{j1}, x2, x, dims2);

                    K = K + (Dti1{1}.*Dtj2{1}' + Dtj1{1}.*Dti2{1}')./Ci{i1}./Ci{j1};
                end
            end
        end
    end
    DKff = {C.*K};
else
    DKff = {0};
end

end

function DKff = gpcf_prod_ginput4(gpcf, x, x2, dims)
%GPCF_PROD_GINPUT4  Evaluate gradient of covariance function with
%                  respect to x. Simplified and faster version of
%                  sexp_ginput, returns full matrices.
%
%  Description
%    DKff = GPCF_PROD_GINPUT4takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to X (whole matrix). This subfunction is needed when
%    using derivative observations.
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%    GPCF_PROD_GINPUT, GPCF_PROD_GINPUT2, GPCF_PROD_CFDG2

ncf = length(gpcf.cf);

% check whether there are any covariance functions that depend on dims
% covariate
ind = [];
for i1 = 1:ncf
    if (~isfield(gpcf.cf{i1}, 'selectedVariables') || any(ismember(gpcf.cf{i1}.selectedVariables,dims)))
        ind(end+1) = i1;
    end
end
% Calculate derivative
if ~isempty(ind)
    % evaluate the individual covariance functions
    C = 1;
    K = 0;
    for i1=1:ncf
        Ci = gpcf.cf{i1}.fh.cov(gpcf.cf{i1}, x, x2);
        C = C.*Ci;
        if any(i1 == ind)
            Dt = gpcf.cf{i1}.fh.ginput4(gpcf.cf{i1}, x, x2, dims);
            if length(Dt)>1
                error('gpcf_prod->ginput4: something is wrong')
            end
            K = K+Dt{1}./Ci;
        end
    end
    DKff = {C.*K};
else
    DKff = {0};
end

end

function C = gpcf_prod_cov(gpcf, x1, x2)
%GP_PROD_COV  Evaluate covariance matrix between two input vectors
%
%  Description
%    C = GP_PROD_COV(GP, TX, X) takes in covariance function of a
%    Gaussian process GP and two matrixes TX and X that contain
%    input vectors to GP. Returns covariance matrix C. Every
%    element ij of C contains covariance between inputs i in TX
%    and j in X. This is a mandatory subfunction used for example in
%    prediction and energy computations.
%
%
%  See also
%    GPCF_PROD_TRCOV, GPCF_PROD_TRVAR, GP_COV, GP_TRCOV

  if isempty(x2)
    x2=x1;
  end
  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  ncf = length(gpcf.cf);

  % evaluate the individual covariance functions
  C = 1;
  for i=1:ncf
    cf = gpcf.cf{i};
    C = C.*cf.fh.cov(cf, x1, x2);
  end
end

function C = gpcf_prod_trcov(gpcf, x)
%GP_PROD_TRCOV     Evaluate training covariance matrix of inputs
%
%  Description
%    C = GP_PROD_TRCOV(GP, TX) takes in covariance function of a
%    Gaussian process GP and matrix TX that contains training
%    input vectors. Returns covariance matrix C. Every element ij
%    of C contains covariance between inputs i and j in TX. This
%    is a mandatory subfunction used for example in prediction and
%    energy computations.
%
%  See also
%    GPCF_PROD_COV, GPCF_PROD_TRVAR, GP_COV, GP_TRCOV
  ncf = length(gpcf.cf);

  % evaluate the individual covariance functions
  C = 1;
  for i=1:ncf
    cf = gpcf.cf{i};
    C = C.*cf.fh.trcov(cf, x);
  end
end

function C = gpcf_prod_trvar(gpcf, x)
% GP_PROD_TRVAR     Evaluate training variance vector
%
%  Description
%    C = GP_PROD_TRVAR(GPCF, TX) takes in covariance function of
%    a Gaussian process GPCF and matrix TX that contains training
%    inputs. Returns variance vector C. Every element i of C
%    contains variance of input i in TX. This is a mandatory
%    subfunction used for example in prediction and energy computations.
%
%  See also
%    GPCF_PROD_COV, GP_COV, GP_TRCOV

  ncf = length(gpcf.cf);

  % evaluate the individual covariance functions
  C = 1;
  for i=1:ncf
    cf = gpcf.cf{i};
    C = C.*cf.fh.trvar(cf, x);
  end
end

function reccf = gpcf_prod_recappend(reccf, ri, gpcf)
%RECAPPEND  Record append
%
%  Description
%    RECCF = GPCF_PROD_RECAPPEND(RECCF, RI, GPCF) takes a
%    covariance function record structure RECCF, record index RI
%    and covariance function structure GPCF with the current MCMC
%    samples of the parameters. Returns RECCF which contains all
%    the old samples and the current samples from GPCF. This
%    subfunction is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC, GP_MC->RECAPPEND

  if nargin == 2
    % Initialize the record
    reccf.type = 'gpcf_prod';

    % Initialize parameters
    ncf = length(ri.cf);
    for i=1:ncf
      cf = ri.cf{i};
      reccf.cf{i} = cf.fh.recappend([], ri.cf{i});
    end

    % Set the function handles
    reccf.fh.pak = @gpcf_prod_pak;
    reccf.fh.unpak = @gpcf_prod_unpak;
    reccf.fh.lp = @gpcf_prod_lp;
    reccf.fh.lpg = @gpcf_prod_lpg;
    reccf.fh.cfg = @gpcf_prod_cfg;
    reccf.fh.cov = @gpcf_prod_cov;
    reccf.fh.trcov  = @gpcf_prod_trcov;
    reccf.fh.trvar  = @gpcf_prod_trvar;
    reccf.fh.recappend = @gpcf_prod_recappend;
  else
    % Append to the record

    % Loop over all of the covariance functions
    ncf = length(gpcf.cf);
    for i=1:ncf
      cf = gpcf.cf{i};
      reccf.cf{i} = cf.fh.recappend(reccf.cf{i}, ri, cf);
    end
  end
end

function [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = gpcf_prod_cf2ss(gpcf,x)
%GPCF_PROD_CF2SS Convert the covariance function to state space form
%
%  Description
%    Convert the sum of two covariance functions to the corresponding
%    sum of two state space models. Details on how this is done can
%    be found in the reference:
%
%  References
%    Arno Solin and Simo Sarkka (2014). Explicit link between periodic
%    covariance functions and state space models. Accepted for
%    publication in Proceedings of the Seventeenth International
%    Conference on Artifcial Intelligence and Statistics (AISTATS 2014).
%

  % Check arguments
  if nargin < 2, x = []; end

  % Vector of function handles of conversion functions
  % from covariance functions to state space
  cf2ssvect = cell(length(gpcf.cf),1);
  for k = 1:length(gpcf.cf)

      % Initial function handle
      fh = @(y) gpcf.cf{k}.fh.cf2ss(gpcf.cf{k},x);

      % Deal with the periodic (deterministic) covariance function
      % as a special case
      if isequal(gpcf.cf{k}.type,'gpcf_periodic') && gpcf.cf{k}.decay == 0
        cf2ssvect{k} = @(y) cf2ss_periodicprod(fh);
      else
        cf2ssvect{k} = fh;
      end

  end

  % Return model matrices, derivatives and parameter information
  [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = ...
      cf_prod_to_ss(cf2ssvect);

end

function [F,L,Pinf1,H,Pinf,dF,dPinf,dPinf1,params] = cf2ss_periodicprod(fh)
% CF2SS_PERIODICPROD - Computes Qc values for periodic covariance production

 [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = fh();
 Pinf1 = Pinf; dPinf1 = dPinf;
end

author: gpstuff-dev

repo: gpstuff

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


text:
function Hp = essential_flat(H)
    %Reshape a [3x6xk] matrix to a [3x3x2k] matrix
    Hp = reshape(H,3,3,[]);
end

author: NicolasBoumal

repo: manopt


text:
function test_sameobject
% Result should be 1 in both cases below.

a = rand(4);
b = a;
if sameobject(a,b) ~= 1
  error('failed');
end
if helper(a,a) ~= 1
  error('failed');
end
disp('Test passed.')

function x = helper(a,b)

x = sameobject(a,b);

author: andrewssobral

repo: lrslibrary


text:
function [ gFile, qFile, vFile, wFile, pFile, rFile ] = tsgMakeFilenames( sGridName )
%
% [ gFile, qFile, vFile, wFile, pFile, rFile ] = tsgMakeFilenames( sGridName )
%
% given the paths from tsgGetPaths() it returns the filenames of all files associated with the name
%

[ sFiles, sTasGrid ] = tsgGetPaths();

gFile = [sFiles,'gFile_',sGridName]; % the filename to store the grid

qFile = [sFiles,'qFile_',sGridName]; % file with the quadrature points

vFile = [sFiles,'vFile_',sGridName]; % file with values to be loaded

wFile = [sFiles,'wFile_',sGridName]; % file with anisotropic weights

pFile = [sFiles,'pFile_',sGridName]; % file with points to evaluate the surrogate

rFile = [sFiles,'rFile_',sGridName]; % file with the result of the evaluation

end

author: johannesgerer

repo: jburkardt-m


text:
function planC = dicomDirToPlanC(dirPath,optS,mergeScansFlag)
% function planC = dicomDirToPlanC(dirPath,optS,mergeScansFlag)
%
% APA, 9/23/2021

% Read options if not passed
if exist('optS','var') && ~isempty(optS) && isstruct(optS) % optS is filename of json file
    %passed optS is already a structure. Use as is.
elseif exist('optS','var') && ~isempty(optS) && ischar(optS)  % optS is filename of json file
    optS = opts4Exe(optS);
else  % default CERROptions.json from CERR distribution
    pathStr = getCERRPath;
    optName = [pathStr 'CERROptions.json'];
    optS = opts4Exe(optName);
end

if ~exist('mergeScansFlag','var')
    mergeScansFlag = 'No';
end

planC = {};
dcmdirS = [];
patientNum = 1;
recursiveFlag = true;
excludePixelDataFlag = true;
patient = scandir_mldcm(dirPath, excludePixelDataFlag, recursiveFlag);
if ~isempty(patient)
    for j = 1:length(patient.PATIENT)
        dcmdirS.(['patient_' num2str(patientNum)]) = patient.PATIENT(j);
        patientNum = patientNum + 1;
    end
end

if isempty(dcmdirS)
    return;
end

patNameC = fieldnames(dcmdirS);

combinedDcmdirS = struct('STUDY',dcmdirS.(patNameC{1}).STUDY,'info',dcmdirS.(patNameC{1}).info);
count = 0;
newCombinedDcmdirS = struct('STUDY','');
for studyCount = 1:length(combinedDcmdirS.STUDY)
    for seriesCount = 1:length(combinedDcmdirS.STUDY(studyCount).SERIES)
        count = count + 1;
        newCombinedDcmdirS.STUDY.SERIES(count) = combinedDcmdirS.STUDY(studyCount).SERIES(seriesCount);
    end
end
combinedDcmdirS = newCombinedDcmdirS;
for i = 2:length(patNameC)
    for  k = 1:length(dcmdirS.(patNameC{i}).STUDY)
        for j = 1:length(dcmdirS.(patNameC{i}).STUDY(k).SERIES)
            combinedDcmdirS.STUDY.SERIES(end+1) = dcmdirS.(patNameC{i}).STUDY(k).SERIES(j);
        end
    end
end
% Pass the java dicom structures to function to create CERR plan
planC = dcmdir2planC(combinedDcmdirS,mergeScansFlag,optS);


text:
%--- help for statespace/estimate ---
%
% ESTIMATE Maximum likelihood parameter estimation of state-space models
%
%  Syntax:
%
%    [EstMdl,estParams,EstParamCov,logL,Output] = estimate(Mdl,Y,params0)
%    [EstMdl,estParams,EstParamCov,logL,Output] = estimate(Mdl,Y,params0,
%                                                          name,value,...)
%
%  Description:
%
%    For observation vector y(t) and state vector x(t), estimate parameters
%    of the following general state-space model by maximum likelihood:
%
%    State equation:       x(t) = A(t) * x(t-1) + B(t) * u(t)
%    Observation equation: y(t) = C(t) * x(t)   + D(t) * e(t)
%
%    where u(t) and e(t) are uncorrelated, unit-variance white noise vector
%    processes. The length of x(t), y(t), u(t), and e(t) is m, n, k, and h,
%    respectively.
%
%    For models created explicitly by specifying coefficients A, B, C, and D,
%    unknown parameters to estimate are identified by the presence of NaNs in
%    these model coefficients. Unknown parameters identified by the presence
%    of NaNs in the mean vector (Mean0) and covariance matrix (Cov0) of initial
%    states x(0) are optionally estimated as well.
%
%    For models created implicitly by specifying a parameter mapping function
%    ParamMap, the mapping function is responsible for managing the presence
%    and placement of unknown parameters. In the implicit approach, the mapping
%    function alone defines the model, and is particularly convenient for
%    estimating complex models and for imposing certain parameter constraints.
%    Moreover, in more general settings in which the initial states are also
%    determined by unknown parameters, ParamMap may include additional output
%    arguments; refer to the SSM/DSSM constructor for more details.
%
%  Input Arguments:
%
%    Mdl - A state-space model with unknown parameters to estimate, created
%      by the SSM/DSSM constructor.
%
%    Y - Observed response data to which the model is fit. For time-invariant
%      models in which the length of each observation vector (n) is the same,
%      Y is a T-by-n matrix. For time-varying models in which the length of
%      the observation vector changes, Y is a T-by-1 cell array in which
%      each element contains a time-varying n-element vector of observations,
%      y(t), associated with the corresponding period. The last observation
%      is the most recent.
%
%    params0 - A vector containing the initial values of unknown
%      parameters associated with model coefficients A, B, C, and D, and
%      optionally the mean vector (Mean0) and covariance matrix (Cov0) of
%      initial states x(0), estimated by maximum likelihood. For models created
%      explicitly, parameters mapped to NaN values are found by a column-wise
%      search of A, followed by B, then C, then D, and finally Mean0 and Cov0.
%      For models created implicitly, the parameter function ParamMap is
%      solely responsible for mapping the initial parameter vector into model
%      coefficients A, B, C, and D, as well as additional information
%      regarding initial states and types if necessary.
%
%  Optional Input Name/Value Pairs:
%
%    'Univariate'  Logical value indicating whether to use the univariate
%                  treatment of a multivariate series. The default is false.
%
%    'SquareRoot'  Logical value indicating whether to use the square-root
%                  filter. The default is false.
%                  SSM only. No effects on DSSM.
%
%    'Tolerance'   A small, non-negative variance tolerance that controls
%                  whether an observed series is ignored if its forecast
%                  uncertainty falls below this threshold. Setting tolerance
%                  to a small number, say 1e-15, may help overcome numerical
%                  problems of the filter. The default is 0.
%
%    'Predictors'  T-by-d matrix of common predictor variables used to
%                  include a regression component in the observation equation.
%                  Observations at time t are deflated such that
%
%                  [y(t) - z(t)*b] = C * x(t) + D * e(t)
%
%                  where z(t) is a vector of predictor variables and b is
%                  the regression coefficient vector (see below). The default
%                  is an empty matrix (no regression component)
%
%    'Beta0'       d-by-n matrix of initial values of regression coefficients
%                  associated with predictors (see above). If the model contains
%                  a regression component, coefficients are estimated along
%                  with other unknown parameters in A, B, C, D, Mean0, and
%                  Cov0; the default initial values are obtained by ordinary
%                  least squares (OLS) by regressing Y on the explanatory
%                  variables.
%
%    'SwitchTime'  a positive integer that specifies the date after which the
%                  diffuse filter switches to the standard filter.
%                  The default is the earliest date when the smoothed
%                  initial states have a full-rank covariance matrix.
%                  DSSM only. No effects on SSM.
%
%    'CovMethod'   String or character vector indicating the method for
%                  computing the asymptotic parameter error covariance
%                  matrix of estimated parameters. Values are:
%
%                  VALUE             METHOD
%
%                  o 'hessian'       Negative inverted Hessian matrix
%                  o 'opg'           Outer product of gradients (default)
%                  o 'sandwich'      Both Hessian and outer product of gradients
%
%    'Options'     Optimization options created with OPTIMOPTIONS. If
%                  specified, default optimization parameters are replaced
%                  by those in options. The default OPTIMOPTIONS object depends
%                  on the optimization function used to estimate parameters.
%                  For constrained optimization, FMINCON is called and the
%                  algorithm used is interior point; for unconstrained
%                  optimization, FMINUNC is called and the algorithm is
%                  quasi-Newton. See documentation for OPTIMOPTIONS, FMINCON,
%                  and FMINUNC for details.
%
%    'Display'    String vector or cell vector of character vectors
%                 indicating what information to display in the command
%                 window. Values are:
%
%                 VALUE           DISPLAY
%
%                 o 'off'         No display to the command window.
%
%                 o 'params'      Display maximum likelihood parameter
%                                 estimates, standard errors, and t statistics.
%                                 This is the default.
%
%                 o 'iter'        Display iterative optimization information.
%
%                 o 'diagnostics' Display optimization diagnostics.
%
%                 o 'full'        Display 'params', 'iter', and 'diagnostics'.
%
%    The following optional name/value pairs are related to constrained
%    optimization performed by FMINCON (see FMINCON for details):
%
%    'Aineq'  Linear inequality matrix, such that for solution vector X
%             Aineq*X <= bineq. The number of rows is determined by the
%             number of constraints, and the number of columns by the number
%             of estimated parameters. Columns are ordered by estimated
%             parameters in A, B, C, D, Mean0, Cov0, and finally regression
%             coefficients for models with a regression component.
%
%    'bineq'  Linear inequality column vector, such that for solution vector
%             X Aineq*X <= bineq.
%
%    'Aeq'    Linear equality matrix, such that for solution vector X
%             Aeq*X = beq. The number of rows is determined by the
%             number of constraints, and the number of columns by the number
%             of estimated parameters. Columns are ordered by estimated
%             parameters in A, B, C, D, Mean0, Cov0, and finally regression
%             coefficients for models with a regression component.
%
%    'beq'    Linear equality column vector, such that for solution vector X
%             Aeq*X = beq.
%
%    'lb'     Lower bounds column vector on estimated parameters. Vector
%             elements are ordered by estimated parameters in A, B, C, D,
%             Mean0, Cov0, and finally regression coefficients for models with
%             a regression component.
%
%    'ub'     Upper bounds column vector on estimated parameters. Vector
%             elements are ordered by estimated parameters in A, B, C, D,
%             Mean0, Cov0, and finally regression coefficients for models with
%             a regression component.
%
%  Output Arguments:
%
%    EstMdl - A fitted state-space model with estimated model coefficients.
%      Provided the optimization converged successfully, the model is explicit
%      in its coefficients A, B, C, D, Mean0, and Cov0 regardless of whether
%      the input model (Mdl) was created explicitly or implicitly.
%
%    estParams - Vector of estimated parameter. The elements are ordered such
%      that any estimated parameters in A appear first, then B, then C, then
%      D, and finally Mean0 and Cov0. Additionally, if the observation equation
%      includes a regression component, then estimates of any regression
%      coefficients appear last.
%
%    EstParamCov - Variance-covariance matrix of estimated parameters. The
%      rows/columns are ordered such that variances/covariances of any
%      estimated parameters in A appear first, then B, then C, then D, and
%      finally Mean0 and Cov0. Additionally, if the observation equation
%      includes a regression component, then variances/covariances of any
%      estimated regression coefficients appear last.
%
%    logL - Log-likelihood of the observations.
%
%    Output - Output structure with the following fields:
%
%       o ExitFlag - Optimization exit flag that describes the exit condition.
%           See FMINCON or FMINUNC for additional details.
%
%       o Options - Optimization options, created with OPTIMOPTIONS, and used
%           by the optimizer.
%
%  Notes:
%
%    o Missing observations in Y are indicated by NaNs.
%
%    o Although creating a model explicitly by directly specifying parameters
%      (A, B, C, D, etc.) with NaN placeholders to indicate parameters to
%      estimate is more convenient than specifying a user-defined mapping
%      function ParamMap, the utility of an explicit approach is limited in
%      that each estimated parameter affects and is uniquely associated with
%      a single element of a coefficient matrix.
%
%    o The option of univariate treatment requires a diagonal D(t)*D(t)'.
%
%    o The state-space model itself does not store regression data or
%      coefficients. Instead, a regression component is used to deflate the
%      observations, such that the deflated data is given by y(t) - z(t)*b.
%
%    o Regression coefficients in the observation equation are estimated along
%      with other unknown parameters, and so parameter constraints may be
%      placed on regression coefficients as well other parameters by specifying
%      appropriate entries in 'Aineq', 'bineq', etc.
%
%    o If observations are multivariate, then the same regressors apply to all.
%
%    o If the state equation requires predictors, or each individual component
%      of the observation equation requires a set of distinct predictors, try
%      one of the following methods:
%
%        o Expand the states by a constant one.
%        o Expand the states by predictors.
%        o Return 8 output arguments from the user-supplied function
%          ParamMap, the last of which is the deflated observation data.
%
%    o Regression components in the observation equation are allowed only for
%      time-invariant observations, in which the input observation series Y
%      is of constant length and specified as a T-by-n matrix. Y specified as
%      a T-by-1 cell array indicates a time-varying observation series whose
%      length may change over time. If the length of each observation y(t)
%      changes, then it is unclear which regression coefficients are needed
%      to deflate a given observation.
%
%    o If no constraints are specified, FMINUNC is used; otherwise, FMINCON
%      is used for constrained optimization. Therefore, when specifying the
%      input Options, the input should be consistent with the solver (see
%      OPTIMOPTIONS, FMINCON, and FMINUNC for details).
%
%      Whenever possible, it is recommended to avoid equality/inequality
%      constraints by reparameterizing the model. For example, various
%      parameters may be exponentiated using EXP to ensure positivity.
%
%    o The observations in the starting periods are treated as if they were
%      presample data for exact initialization of the diffuse initial
%      states. We define the likelihood function as the joint density
%      of the observations after the algorithm is switched to
%      the standard Kalman filter. The variable 'SwitchTime' determines when
%      the standard Kalman filter starts.
%
%  See also SSM, DSSM, FILTER, SMOOTH, FORECAST, SIMULATE, SIMSMOOTH.
%
%    Other functions named estimate
%
%       abstvar/estimate         garch/estimate
%       arima/estimate           generic/estimate
%       conjugateblm/estimate    gjr/estimate
%       customblm/estimate       regARIMA/estimate
%       diffuseblm/estimate      semiconjugateblm/estimate
%       dsge/estimate            ssm/estimate
%       dssm/estimate            varm/estimate
%       egarch/estimate          vecm/estimate
%       empiricalblm/estimate
%

author: jmaih


text:
function fnamefullpath = findfile(rootdirexamples, fname)

fnamefullpath = [rootdirexamples, fname];
if ~exist(fnamefullpath, 'file')
    msg = sprintf('%s does not exist. You first have to generate it by running demo_snirf\n', fnamefullpath);
    MenuBox(msg, 'OK');
    fnamefullpath = '';
    return;
end

author: BUNPC


text:
function vr = videoReader(url, varargin)
% videoReader class constructor
%   Creates a object that reads video streams.  We use a plugin
%   architecture in the backend to do the actual reading.  For example,
%   on Windows, DirectShow will typically be used and on Linux, the
%   ffmpeg library is often used.
%
%   vr = videoReader(url)
%     Opens the given video file for reading using the default plugin.
%     On Windows, 'DirectShow' is used by default and on Linux,
%     'ffmpegPopen2' is used by default.  For most plugins, the url will
%     really be a filename.
%
%   vr = videoReader(url, ..., 'plugin',pluginName, ...)
%   vr = videoReader(url,pluginName, ...)
%     Opens the file using the manually specified plugin implementation.
%     Available plugins include:
%
%     'DirectShow': preferred method on Windows
%       - Only available on Windows
%       - See INSTALL.dshow.txt for installation instructions
%       - Can load virtually any video file that can be played in
%         Microsoft's Windows Media Player.  Note that the correct codec
%         must be installed to read a file.  For example, to read
%         tests/numbers.3ivx.avi, the user must have installed an MPEG4
%         codec such as 3ivx (www.3ivx.com), DivX (www.divx.com), or XviD
%         (www.xvid.org).
%       - The URL parameter should be a filename.
%         - As a convenience, all forward slashes ('/') are automatically
%           converted to backslashes ('\')
%
%     'ffmpegPopen2': safe method on Linux
%       - Only supported on GNU/Linux (might work on BSD systems too like Mac
%         OS X, but this is untested)
%       - See INSTALL.ffmpeg.txt for installation instructions
%       - Creates a separate server process to communicate with the
%         ffmpeg libraries.
%         - Works when the system's version of GCC is very different from
%           the one that MathWorks used to compile Matlab.
%         - Isolates ffmpeg errors so they typically cannot crash
%           Matlab.
%         - May allow for more flexible distribution terms for your code
%           when it uses videoIO (ffmpeg may be compiled with either
%           the LGPL or GPL license).
%
%     'ffmpegDirect': low-overhead method on Linux
%       - same as ffmpegPopen2, but the ffmpeg libraries are loaded
%         directly by the MEX file.
%         - May not work if MathWorks' and your version of GCC are
%           incompatible.
%         - Slightly faster than ffmpegPopen2 since there is no
%           interprocess communication overhead.
%
%   vr = videoReader(url, ..., param,arg,...)
%     Allows the user to pass extra configuration arguments to plugin.
%     Currently no plugin arguments are supported right now.  In the
%     future, we may allow the user to do things like have DirectShow
%     automatically convert to grayscale, or give options to trade off
%     speed with seeking precision.
%
% Once you have created a videoReader object, you must next call NEXT,
% SEEK, or STEP at least once so that it will read some frame from disk.
% *After* calling one of these methods, you may call GETFRAME as many
% times as you would like and it will decode and return the current frame
% (without advancing to a different frame).  GETINFO may be called at any
% time (even before NEXT, SEEK, or STEP).  It returns basic information
% about the video stream.  Once you are done using the videoReader, make
% sure you call CLOSE so that any system resources allocated by the plugin
% may be released.  Here's a simple example of how you might use
% videoReader:
%
%   % take us to the videoReader directory since we know there's a video
%   % file there.
%   chdir(fileparts(which('videoReaderWrapper.cpp')));
%
%   % Construct a videoReader object
%   vr = videoReader('tests/numbers.uncompressed.avi');
%
%   % Do some processing on the video and display the results
%   avgIntensity = [];
%   i = 1;
%   figure;
%   while (next(vr))
%     img = getframe(vr);
%     avgIntensity(i) = mean(img(:));
%     subplot(121); imshow(img);        title('current frame');
%     subplot(122); plot(avgIntensity); title('avg. intensity vs. frame');
%     drawnow;      pause(0.1);         i = i+1;
%   end
%   vr = close(vr);
%
% SEE ALSO:
%   buildVideoMex
%   videoReader/close
%   videoReader/getframe
%   videoReader/getinfo
%   videoReader/getnext
%   videoReader/next
%   videoReader/seek
%   videoReader/step
%   videoWriter
%
%Copyright (c) 2006 Gerald Dalley
%See "MIT.txt" in the installation directory for licensing details (especially
%when using this library on GNU/Linux).

if (mod(length(varargin),2) == 0)
  plugin     = defaultVideoIOPlugin;
  pluginArgs = varargin;
else
  plugin     = varargin{1};
  pluginArgs = {varargin{2:end}};
end
[plugin,pluginArgs] = parsePlugin(plugin, pluginArgs);

vr = struct('plugin',mexName(plugin), 'handle',int32(-1));
vr = class(vr, 'videoReader');
[pathstr, name, ext, versn] = fileparts(url);
strArgs = cell(size(pluginArgs));
for i=1:numel(pluginArgs), strArgs{i} = num2str(pluginArgs{i}); end
vr.handle = feval(vr.plugin, 'open', vr.handle, ...
  fullfile(pathstr,[name ext versn]), strArgs{:});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function n = mexName(plugin)
n = ['videoReader_' plugin];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [plugin,pluginArgs] = parsePlugin(plugin, pluginArgs)
if (length(pluginArgs) > 0)
  [pluginSpecified,idx] = ismember('plugin', {pluginArgs{1:2:end}});
  if pluginSpecified
    plugin = pluginArgs{idx*2};
    pluginArgs = { pluginArgs{1:idx*2-2}, pluginArgs{idx*2+1:end} };
  end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [pos, tlight] = littim (tjd, idbody, pose, tlite)

% this function computes the position of a solar system body,
% as antedated for light-time.

%      tjd    = tdb julian date of observation (in)

%      idbody = id number of body, used in calls to solsys (in)

%      pose   = position vector of observer (or the geocenter),
%               with respect to origin at solar system barycenter,
%               referred to icrs axes, components in au (in)

%      tlite  = first approximation to light-time, in days (in)
%               (can be set to 0.0d0 if unknown)

%      pos    = position vector of body, with respect to origin at
%               observer (or the geocenter), referred to icrs axes,
%               components in au (out)

%      tlight = final light-time, in days (out)

% ported from NOVAS 3.1

%%%%%%%%%%%%%%%%%%%%%%%

% set light-time convergence tolerance

tol = 1.0d-9;

t0 = 0.0d0;

t1 = tjd - t0;

t2 = t1 - tlite;

% iterate to obtain correct light-time (usually converges rapidly)

for iter = 1:10

    [pos1, vel1, ierr] = solsys (t2, str2num(idbody), 0);

    [pos, tlight] = geocen (pos1, pose);

    if (ierr ~= 0)

        fprintf ('\nplace: cannot obtain coordinates of object at jd %16.8f', t0 + t2);

        return

    end

    t3 = t1 - tlight;

    if (abs(t3 - t2) > tol)

        t2 = t3;

    else

        break

    end

end

author: goGPS-Project


text:
function v = reshape(v,varargin)
% overloads reshape

v.x = reshape(v.x,varargin{:});
v.y = reshape(v.y,varargin{:});
v.z = reshape(v.z,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function conv2conv_f()
    global config mem;
    curr_layer_idx = config.misc.current_layer;
    mem.layer_inputs{curr_layer_idx} = mem.activations{curr_layer_idx-1};
end

author: jimmy-ren


text:
function [v]=arsim(w,A,C,n,ndisc)
%ARSIM	Simulation of AR process.
%
%  v=ARSIM(w,A,C,n) simulates n time steps of the AR(p) process
%
%     v(k,:)' = w' + A1*v(k-1,:)' +...+ Ap*v(k-p,:)' + eta(k,:)',
%
%  where A=[A1 ... Ap] is the coefficient matrix, and w is a vector of
%  intercept terms that is included to allow for a nonzero mean of the
%  process. The vectors eta(k,:) are independent Gaussian noise
%  vectors with mean zero and covariance matrix C.
%
%  The p vectors of initial values for the simulation are taken to
%  be equal to the mean value of the process. (The process mean is
%  calculated from the parameters A and w.) To avoid spin-up effects,
%  the first 10^3 time steps are discarded. Alternatively,
%  ARSIM(w,A,C,n,ndisc) discards the first ndisc time steps.

%  Modified 13-Oct-00
%  Author: Tapio Schneider
%          tapio@gps.caltech.edu

  m       = size(C,1);                  % dimension of state vectors
  p       = size(A,2)/m;                % order of process

  if (p ~= round(p))
    error('Bad arguments.');
  end

  if (length(w) ~= m | min(size(w)) ~= 1)
    error('Dimensions of arguments are mutually incompatible.')
  end
  w       = w(:)';                      % force w to be row vector

  % Check whether specified model is stable
  A1 	  = [A; eye((p-1)*m) zeros((p-1)*m,m)];
  lambda  = eig(A1);
  if any(abs(lambda) > 1)
    warning('The specified AR model is unstable.')
  end

  % Discard the first ndisc time steps; if ndisc is not given as input
  % argument, use default
  if (nargin < 5)
    ndisc = 10^3;
  end

  % Compute Cholesky factor of covariance matrix C
  [R, err]= chol(C);                    % R is upper triangular
  if err ~= 0
    error('Covariance matrix not positive definite.')
  end

  % Get ndisc+n independent Gaussian pseudo-random vectors with
  % covariance matrix C=R'*R
  randvec = randn([ndisc+n,m])*R;

  % Add intercept vector to random vectors
  randvec = randvec + ones(ndisc+n,1)*w;

  % Get transpose of system matrix A (use transpose in simulation because
  % we want to obtain the states as row vectors)
  AT      = A';

  % Take the p initial values of the simulation to equal the process mean,
  % which is calculated from the parameters A and w
  if any(w)
    %  Process has nonzero mean    mval = inv(B)*w'    where
    %             B = eye(m) - A1 -... - Ap;
    %  Assemble B
    B 	 = eye(m);
    for j=1:p
      B = B - A(:, (j-1)*m+1:j*m);
    end
    %  Get mean value of process
    mval = w / B';

    %  The optimal forecast of the next state given the p previous
    %  states is stored in the vector x. The vector x is initialized
    %  with the process mean.
    x    = ones(p,1)*mval;
  else
    %  Process has zero mean
    x    = zeros(p,m);
  end

  % Initialize state vectors
  u      = [x; zeros(ndisc+n,m)];

  % Simulate n+ndisc observations. In order to be able to make use of
  % Matlab's vectorization capabilities, the cases p=1 and p>1 must be
  % treated separately.
  if p==1
    for k=2:ndisc+n+1;
      x(1,:) = u(k-1,:)*AT;
      u(k,:) = x + randvec(k-1,:);
    end
  else
    for k=p+1:ndisc+n+p;
      for j=1:p;
	x(j,:) = u(k-j,:)*AT((j-1)*m+1:j*m,:);
      end
      u(k,:) = sum(x)+randvec(k-p,:);
    end
  end

  % return only the last n simulated state vectors
  v = u(ndisc+p+1:ndisc+n+p,:);


text:
function bohach2_test ( )

%*****************************************************************************80
%
%% BOHACH2_TEST works with the Bohachevsky function #2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOHACH2_TEST:\n' );
  fprintf ( 1, '  Test COORDINATE_SEARCH with the Bohachevsky function #2.\n' );
  n = 2;

  x = [ 0.6, 1.3 ];
  r8vec_print ( n, x, '  Initial point X0:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', bohach2 ( x ) );

  flag = 0;
  x = coordinate_search ( x, @bohach2, flag );
  r8vec_print ( n, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g\n', bohach2 ( x ) );
%
%  Demonstrate correct minimizer.
%
  x = [ 0.0, 0.0 ];
  r8vec_print ( n, x, '  Correct minimizer X*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', bohach2 ( x ) );

  return
end


text:
function val = get_params(CPD, name)
% GET_PARAMS Get the parameters (fields) for a gaussian_CPD object
% val = get_params(CPD, name)
%
% The following fields can be accessed
%
% mean       - mu(:,i) is the mean given Q=i
% cov        - Sigma(:,:,i) is the covariance given Q=i
% weights    - W(:,:,i) is the regression matrix given Q=i
%
% e.g., mean = get_params(CPD, 'mean')

switch name
 case 'mean',      val = CPD.mean;
 case 'cov',       val = CPD.cov;
 case 'weights',   val = CPD.weights;
 otherwise,
  error(['invalid argument name ' name]);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function [bb,conf] = tldExample(opt)

global tld; % holds results and temporal variables

% INITIALIZATION ----------------------------------------------------------

opt.source = tldInitSource(opt.source); % select data source, camera/directory

figure(2); set(2,'KeyPressFcn', @handleKey); % open figure for display of results
finish = 0; function handleKey(~,~), finish = 1; end % by pressing any key, the process will exit

while 1
    source = tldInitFirstFrame(tld,opt.source,opt.model.min_win); % get initial bounding box, return 'empty' if bounding box is too small
    if ~isempty(source), opt.source = source; break; end % check size
end

tld = tldInit(opt,[]); % train initial detector and initialize the 'tld' structure
tld = tldDisplay(0,tld); % initialize display

% RUN-TIME ----------------------------------------------------------------

for i = 2:length(tld.source.idx) % for every frame

    tld = tldProcessFrame(tld,i); % process frame i
    tldDisplay(1,tld,i); % display results on frame i

    if finish % finish if any key was pressed
        if tld.source.camera
            stoppreview(tld.source.vid);
            closepreview(tld.source.vid);
             close(1);
        end
        close(2);
        bb = tld.bb; conf = tld.conf; % return results
        return;
    end

    if tld.plot.save == 1
        img = getframe;
        imwrite(img.cdata,[tld.output num2str(i,'%05d') '.png']);
    end

end

bb = tld.bb; conf = tld.conf; % return results

end

repo: OpenTLD


text:
% std_findsameica() - find groups of datasets with identical ICA decomposiotions
%
% Usage:
%        >> clusters = std_findsameica(ALLEEG);
% Inputs:
%   ALLEEG  - a vector of loaded EEG dataset structures of all sets in the STUDY set.
%
% Outputs:
%   cluster - cell array of groups of datasets
%
% Authors:  Arnaud Delorme, SCCN, INC, UCSD, July 2009-

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, arno@sccn.ucsd.edu
%
% 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

% Coding notes: Useful information on functions and global variables used.

function cluster = std_findsameica(ALLEEG);

cluster = { [1] };
for index = 2:length(ALLEEG)

    found = 0;
    for c = 1:length(cluster)
        if all(size(ALLEEG(cluster{c}(1)).icaweights) == size(ALLEEG(index).icaweights))
            %if isequal(ALLEEG(cluster{c}(1)).icaweights, ALLEEG(index).icaweights)
            if sum(sum(abs(ALLEEG(cluster{c}(1)).icaweights-ALLEEG(index).icaweights))) < 2e-5
                cluster{c}(end+1) = index;
                found = 1;
                break;
            end;
        end;
    end;
    if ~found
        cluster{end+1} = index;
    end;
end;

author: PatternRecognition

repo: OpenBMI


text:
function test_sep
%TEST_SEP test cs_sep, and compare with Gilbert's meshpart vtxsep
% (requires MESHPART).
%
% Example:
%   test_sep
%
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

clear functions

index = UFget ;
[ignore f] = sort (max (index.nrows, index.ncols)) ;

clf

for k = 1:length(f)

    i = f (k) ;
    Prob = UFget (i) ;
    disp (Prob) ;
    A = spones (Prob.A) ;
    [m n] = size (A) ;
    if (m ~= n)
        A = A'*A ;
    end

    A = A|A' ;

    p = symrcm (A) ;

    n = size (A,1) ;
    n2 = fix (n/2) ;
    a = p (1:n2) ;
    b = p ((n2+1):n) ;

    clf

    subplot (2,3,1) ; spy (A) ;
    subplot (2,3,2) ; spy (A (p,p)) ;

    hold on
    plot ([.5 n2+.5 n2+.5 .5 .5], [.5 .5 n2+.5 n2+.5 .5], 'r', 'LineWidth', 2) ;
    hold off

    subplot (2,3,3) ; spy (A (a,b)) ; title ('edge sep') ;
    subplot (2,3,6) ; cs_dmspy (A (a,b)) ; title ('node sep') ;

    [s as bs] = vtxsep (A,a,b) ;                                        %#ok
    [s2 a2 b2] = cs_sep (A,a,b) ;

    p2 = [a2 b2 s2] ;
    B = A (p2,p2) ;
    subplot (2,3,5) ; spy (B) ;
    hold on

    px = [s2 a2 b2] ;
    if (any (sort (px) ~= 1:n))
        px      %#ok
        n       %#ok
        error ('!') ;
    end

    na = length (a2) ;
    nb = length (b2) ;
    ns = length (s2) ;                                                  %#ok

    nab = na + nb ;

    plot ([.5 na+.5 na+.5 .5 .5], [.5 .5 na+.5 na+.5 .5], 'r', 'LineWidth', 2) ;

    plot ([na nab nab na na]+0.5, [na na nab nab na]+0.5, 'r', 'LineWidth', 2) ;

    plot ([.5 nab+.5 nab+.5 .5 .5], [.5 .5 nab+.5 nab+.5 .5], 'g', 'LineWidth', 1) ;

    hold off

    nz1 = nnz (A (a2,b2)) ;
    if (nz1 ~= 0)
        nz1     %#ok
        error ('!') ;
    end

    nz2 = nnz (A (a2,b2)) ;
    if (nz2 ~= 0)
        nz2     %#ok
        error ('!') ;
    end

    if (length (s) ~= length (s2))
        fprintf ('lengths differ: %d %d\n', length (s), length (s2)) ;
    end

    drawnow
    % pause

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef dbPitts < dbBase

    methods

        function db= dbPitts(fullSize, whichSet)
            % fullSize is: true or false
            % whichSet is one of: train, val, test

            assert( ismember(whichSet, {'train', 'val', 'test'}) );

            if fullSize
                sizeSuffix= '250k';
            else
                sizeSuffix= '30k';
            end

            db.name= sprintf('pitts%s_%s', sizeSuffix, whichSet);

            paths= localPaths();
            dbRoot= paths.dsetRootPitts;
            db.dbPath= [dbRoot, 'images/'];
            db.qPath= [dbRoot, 'queries/'];

            db.dbLoad();
        end

    end

end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad