text:
function [Y, optinf] = cbpdn(D, S, lambda, opt)

% cbpdn -- Convolutional Basis Pursuit DeNoising
%
%         argmin_{x_m} (1/2)||\sum_m d_m * x_m - s||_2^2 +
%                           lambda \sum_m ||x_m||_1
%
%         The solution is computed using an ADMM approach (see
%         boyd-2010-distributed) with efficient solution of the main
%         linear systems (see wohlberg-2016-efficient).
%
% Usage:
%       [Y, optinf] = cbpdn(D, S, lambda, opt);
%
% Input:
%       D           Dictionary filter set (3D array)
%       S           Input image
%       lambda      Regularization parameter
%       opt         Algorithm parameters structure
%
% Output:
%       Y           Dictionary coefficient map set (3D array)
%       optinf      Details of optimisation
%
%
% Options structure fields:
%   Verbose          Flag determining whether iteration status is displayed.
%                    Fields are iteration number, functional value,
%                    data fidelity term, l1 regularisation term, and
%                    primal and dual residuals (see Sec. 3.3 of
%                    boyd-2010-distributed). The value of rho is also
%                    displayed if options request that it is automatically
%                    adjusted.
%   MaxMainIter      Maximum main iterations
%   AbsStopTol       Absolute convergence tolerance (see Sec. 3.3.1 of
%                    boyd-2010-distributed)
%   RelStopTol       Relative convergence tolerance (see Sec. 3.3.1 of
%                    boyd-2010-distributed)
%   L1Weight         Weighting array for coefficients in l1 norm of X
%   Y0               Initial value for Y
%   U0               Initial value for U
%   rho              ADMM penalty parameter
%   AutoRho          Flag determining whether rho is automatically updated
%                    (see Sec. 3.4.1 of boyd-2010-distributed)
%   AutoRhoPeriod    Iteration period on which rho is updated
%   RhoRsdlRatio     Primal/dual residual ratio in rho update test
%   RhoScaling       Multiplier applied to rho when updated
%   AutoRhoScaling   Flag determining whether RhoScaling value is
%                    adaptively determined (see wohlberg-2015-adaptive). If
%                    enabled, RhoScaling specifies a maximum allowed
%                    multiplier instead of a fixed multiplier.
%   RhoRsdlTarget    Residual ratio targeted by auto rho update policy.
%   StdResiduals     Flag determining whether standard residual definitions
%                    (see Sec 3.3 of boyd-2010-distributed) are used instead
%                    of normalised residuals (see wohlberg-2015-adaptive)
%   RelaxParam       Relaxation parameter (see Sec. 3.4.3 of
%                    boyd-2010-distributed)
%   NonNegCoef       Flag indicating whether solution should be forced to
%                    be non-negative
%   NoBndryCross     Flag indicating whether all solution coefficients
%                    corresponding to filters crossing the image boundary
%                    should be forced to zero.
%   AuxVarObj        Flag determining whether objective function is computed
%                    using the auxiliary (split) variable
%   HighMemSolve     Use more memory for a slightly faster solution
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2015-12-28
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if nargin < 4,
  opt = [];
end
checkopt(opt, defaultopts([]));
opt = defaultopts(opt);

% Set up status display for verbose operation
hstr = 'Itn   Fnc       DFid      l1        r         s      ';
sfms = '%4d %9.2e %9.2e %9.2e %9.2e %9.2e';
nsep = 54;
if opt.AutoRho,
  hstr = [hstr '   rho  '];
  sfms = [sfms ' %9.2e'];
  nsep = nsep + 10;
end
if opt.Verbose && opt.MaxMainIter > 0,
  disp(hstr);
  disp(char('-' * ones(1,nsep)));
end

% Collapsing of trailing singleton dimensions greatly complicates
% handling of both SMV and MMV cases. The simplest approach would be
% if S could always be reshaped to 4d, with dimensions consisting of
% image rows, image cols, a single dimensional placeholder for number
% of filters, and number of measurements, but in the single
% measurement case the third dimension is collapsed so that the array
% is only 3d.
if size(S,3) > 1 && size(S,4) == 1,
  xsz = [size(S,1) size(S,2) size(D,3) size(S,3)];
  % Insert singleton 3rd dimension (for number of filters) so that
  % 4th dimension is number of images in input s volume
  S = reshape(S, [size(S,1) size(S,2) 1 size(S,3)]);
else
  xsz = [size(S,1) size(S,2) size(D,3) size(S,4)];
end

% Start timer
tstart = tic;

% Compute filters in DFT domain
Df = fft2(D, size(S,1), size(S,2));
% Convolve-sum and its Hermitian transpose
Dop = @(x) sum(bsxfun(@times, Df, x), 3);
DHop = @(x) bsxfun(@times, conj(Df), x);
% Compute signal in DFT domain
Sf = fft2(S);
% S convolved with all filters in DFT domain
DSf = DHop(Sf);

% Default lambda is 1/10 times the lambda value beyond which the
% solution is a zero vector
if nargin < 3 | isempty(lambda),
  b = ifft2(DHop(Sf), 'symmetric');
  lambda = 0.1*max(vec(abs(b)));
end

% Set up algorithm parameters and initialise variables
rho = opt.rho;
if isempty(rho), rho = 50*lambda+1; end;
if isempty(opt.RhoRsdlTarget),
  if opt.StdResiduals,
    opt.RhoRsdlTarget = 1;
  else
    opt.RhoRsdlTarget = 1 + (18.3).^(log10(lambda) + 1);
  end
end
if opt.HighMemSolve,
  C = bsxfun(@rdivide, Df, sum(Df.*conj(Df), 3) + rho);
else
  C = [];
end
Nx = prod(xsz);
optinf = struct('itstat', [], 'opt', opt);
r = Inf;
s = Inf;
epri = 0;
edua = 0;

% Initialise main working variables
X = [];
if isempty(opt.Y0),
  Y = zeros(xsz, class(S));
else
  Y = opt.Y0;
end
Yprv = Y;
if isempty(opt.U0),
  if isempty(opt.Y0),
    U = zeros(xsz, class(S));
  else
    U = (lambda/rho)*sign(Y);
  end
else
  U = opt.U0;
end

% Main loop
k = 1;
while k <= opt.MaxMainIter && (r > epri | s > edua),
%while k <= opt.MaxMainIter ,
  % Solve X subproblem
  Xf = solvedbi_sm(Df, rho, DSf + rho*fft2(Y - U), C);
  X = ifft2(Xf, 'symmetric');

  % See pg. 21 of boyd-2010-distributed
  if opt.RelaxParam == 1,
    Xr = X;
  else
    Xr = opt.RelaxParam*X + (1-opt.RelaxParam)*Y;
  end

  % Solve Y subproblem
  Y = shrink(Xr + U, (lambda/rho)*opt.L1Weight);
  if opt.NonNegCoef,
    Y(Y < 0) = 0;
  end
  if opt.NoBndryCross,
    Y((end-size(D,1)+2):end,:,:,:) = 0;
    Y(:,(end-size(D,2)+2):end,:,:) = 0;
  end

  % Update dual variable
  U = U + Xr - Y;

  % Compute data fidelity term in Fourier domain (note normalisation)
  if opt.AuxVarObj,
    Yf = fft2(Y); % This represents unnecessary computational cost
    Jdf = sum(vec(abs(sum(bsxfun(@times,Df,Yf),3)-Sf).^2))/(2*xsz(1)*xsz(2));
    Jl1 = sum(abs(vec(bsxfun(@times, opt.L1Weight, Y))));
  else
    Jdf = sum(vec(abs(sum(bsxfun(@times,Df,Xf),3)-Sf).^2))/(2*xsz(1)*xsz(2));
    Jl1 = sum(abs(vec(bsxfun(@times, opt.L1Weight, X))));
  end
  Jfn = Jdf + lambda*Jl1;

  nX = norm(X(:)); nY = norm(Y(:)); nU = norm(U(:));
  if opt.StdResiduals,
    % See pp. 19-20 of boyd-2010-distributed
    r = norm(vec(X - Y));
    s = norm(vec(rho*(Yprv - Y)));
    epri = sqrt(Nx)*opt.AbsStopTol+max(nX,nY)*opt.RelStopTol;
    edua = sqrt(Nx)*opt.AbsStopTol+rho*nU*opt.RelStopTol;
  else
    % See wohlberg-2015-adaptive
    r = norm(vec(X - Y))/max(nX,nY);
    s = norm(vec(Yprv - Y))/nU;
    epri = sqrt(Nx)*opt.AbsStopTol/max(nX,nY)+opt.RelStopTol;
    edua = sqrt(Nx)*opt.AbsStopTol/(rho*nU)+opt.RelStopTol;
  end

  % Record and display iteration details
  tk = toc(tstart);
  optinf.itstat = [optinf.itstat; [k Jfn Jdf Jl1 r s epri edua rho tk]];
  if opt.Verbose,
    if opt.AutoRho,
      disp(sprintf(sfms, k, Jfn, Jdf, Jl1, r, s, rho));
    else
      disp(sprintf(sfms, k, Jfn, Jdf, Jl1, r, s));
    end
  end

  % See wohlberg-2015-adaptive and pp. 20-21 of boyd-2010-distributed
  if opt.AutoRho,
    if k ~= 1 && mod(k, opt.AutoRhoPeriod) == 0,
      if opt.AutoRhoScaling,
        rhomlt = sqrt(r/(s*opt.RhoRsdlTarget));
        if rhomlt < 1, rhomlt = 1/rhomlt; end
        if rhomlt > opt.RhoScaling, rhomlt = opt.RhoScaling; end
      else
        rhomlt = opt.RhoScaling;
      end
      rsf = 1;
      if r > opt.RhoRsdlTarget*opt.RhoRsdlRatio*s, rsf = rhomlt; end
      if s > (opt.RhoRsdlRatio/opt.RhoRsdlTarget)*r, rsf = 1/rhomlt; end
      rho = rsf*rho;
      U = U/rsf;
      if opt.HighMemSolve && rsf ~= 1,
        C = bsxfun(@rdivide, Df, sum(Df.*conj(Df), 3) + rho);
      end
    end
  end

  Yprv = Y;
  k = k + 1;

end

% Record run time and working variables
optinf.runtime = toc(tstart);
optinf.X = X;
optinf.Xf = Xf;
optinf.Y = Y;
optinf.U = U;
optinf.lambda = lambda;
optinf.rho = rho;

% End status display for verbose operation
if opt.Verbose && opt.MaxMainIter > 0,
  disp(char('-' * ones(1,nsep)));
end

return

function u = vec(v)

  u = v(:);

return

function u = shrink(v, lambda)

  if isscalar(lambda),
    u = sign(v).*max(0, abs(v) - lambda);
  else
    u = sign(v).*max(0, bsxfun(@minus, abs(v), lambda));
  end

return

function opt = defaultopts(opt)

  if ~isfield(opt,'Verbose'),
    opt.Verbose = 0;
  end
  if ~isfield(opt,'MaxMainIter'),
    opt.MaxMainIter = 1000;
  end
  if ~isfield(opt,'AbsStopTol'),
    opt.AbsStopTol = 0;
  end
  if ~isfield(opt,'RelStopTol'),
    opt.RelStopTol = 1e-4;
  end
  if ~isfield(opt,'L1Weight'),
    opt.L1Weight = 1;
  end
  if ~isfield(opt,'Y0'),
    opt.Y0 = [];
  end
  if ~isfield(opt,'U0'),
    opt.U0 = [];
  end
  if ~isfield(opt,'rho'),
    opt.rho = [];
  end
  if ~isfield(opt,'AutoRho'),
    opt.AutoRho = 1;
  end
  if ~isfield(opt,'AutoRhoPeriod'),
    opt.AutoRhoPeriod = 1;
  end
  if ~isfield(opt,'RhoRsdlRatio'),
    opt.RhoRsdlRatio = 1.2;
  end
  if ~isfield(opt,'RhoScaling'),
    opt.RhoScaling = 100;
  end
  if ~isfield(opt,'AutoRhoScaling'),
    opt.AutoRhoScaling = 1;
  end
  if ~isfield(opt,'RhoRsdlTarget'),
    opt.RhoRsdlTarget = [];
  end
  if ~isfield(opt,'StdResiduals'),
    opt.StdResiduals = 0;
  end
  if ~isfield(opt,'RelaxParam'),
    opt.RelaxParam = 1.8;
  end
  if ~isfield(opt,'NonNegCoef'),
    opt.NonNegCoef = 0;
  end
  if ~isfield(opt,'NoBndryCross'),
    opt.NoBndryCross = 0;
  end
  if ~isfield(opt,'AuxVarObj'),
    opt.AuxVarObj = 0;
  end
  if ~isfield(opt,'HighMemSolve'),
    opt.HighMemSolve = 0;
  end

return

author: thfylsty


text:
function [N,X,sp] = histogramDistributionPlot(varargin)
% HISTOGRAM generates a histogram using the "optimal" number of bins
%
% If called with no output argument, histogram plots into the current axes
%
% SYNOPSIS [N,X,sp] = histogram(data,factor,normalize)
%          [...] = histogram(data,'smooth')
%          [...] = histogram(axesHandle,...)
%
% INPUT    data: vector of input data
%          factor: (opt) factor by which the bin-widths are multiplied
%                   if 'smooth' (or 's'), a smooth histogram will be formed.
%                   (requires the spline toolbox). For an alternative
%                   approach to a smooth histogram, see ksdensity.m
%                   if 'discrete' (or 'd'), the data is assumed to be a discrete
%                   collection of values. Note that if every data point is,
%                   on average, repeated at least 3 times, histogram will
%                   consider it a discrete distribution automatically.
%                   if 'continuous' (or 'c'), histogram is not automatically
%                   checking for discreteness.
%          normalize : if 1 (default), integral of histogram equals number
%                       data points. If 0, height of bins equals counts.
%                       This option is exclusive to non-"smooth" histograms
%          axesHandle: (opt) if given, histogram will be plotted into these
%                       axes, even if output arguments are requested
%
% OUTPUT   N   : number of points per bin (value of spline)
%          X   : center position of bins (sorted input data)
%          sp  : definition of the smooth spline
%
% REMARKS: The smooth histogram is formed by calculating the cumulative
%           histogram, fitting it with a smoothening spline and then taking
%           the analytical derivative. If the number of data points is
%           markedly above 1000, the spline is fitting the curve too
%           locally, so that the derivative can have huge peaks. Therefore,
%           only 1000-1999 points are used for estimation.
%           Note that the integral of the spline is almost exactly the
%           total number of data points. For a standard histogram, the sum
%           of the hights of the bins (but not their integral) equals the
%           total number of data points. Therefore, the counts might seem
%           off.
%
%           WARNING: If there are multiples of the minimum value, the
%           smooth histogram might get very steep at the beginning and
%           produce an unwanted peak. In such a case, remove the
%           multiple small values first (for example, using isApproxEqual)
%
%
% c: 2/05 jonas
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% test input
if nargin < 1
    error('not enough input arguments for histogram')
end

% check for axes handle
if length(varargin{1}) == 1 && ishandle(varargin{1});
    axesHandle = varargin{1};
    varargin(1) = [];
else
    % ensure compatibility to when axesHandle was given as last input
    if nargin == 3 && ishandle(varargin{end}) && varargin{end} ~= 0
        axesHandle = varargin{end};
        varargin(end) = [];
    else
        axesHandle = 0;
    end
end

% assign data
numArgIn = length(varargin);
data = varargin{1};
data = data(:);

% check for non-finite data points
data(~isfinite(data)) = [];

% check for "factor"
if numArgIn < 2 || isempty(varargin{2})
    factor = 1;
else
    factor = varargin{2};
end
if ischar(factor)
    switch factor
        case {'smooth','s'}
        factor = -1;
        case {'discrete','d'}
            factor = -2;
        case {'continuous','c'}
            factor = -3;
    otherwise
        error('The only string inputs permitted for histogram.m are ''smooth'',''discrete'', or ''continuous''')
    end
else
    % check for normalize, but do so only if there is no "smooth". Note
    % that numArgIn is not necessarily equal to nargin
    if numArgIn < 3 || isempty(varargin{3})
        normalize = true;
    else
        normalize = varargin{3};
    end
end

% doPlot is set to 1 for now. We change it to 0 below if necessary.
doPlot = 1;

nData = length(data);
% check whether we do a standard or a smooth histogram
if factor ~= -1
    % check for discrete distribution
    [xx,nn] = countEntries(data);
    % consider the distribution discrete if there are, on average, 3
    % entries per bin
    nBins = length(xx);
    if factor == -2 || (factor ~= -3 && nBins*3 < nData)
        % discrete distribution.
        nn = nn';
        xx = xx';
    else
        % not a discrete distribution
        if nData < 20
            warning('HISTOGRAM:notEnoughDataPoints','Less than 20 data points!')
            nBins = ceil(nData/4);
        else

            % create bins with the optimal bin width
            % W = 2*(IQD)*N^(-1/3)
            interQuartileDist = diff(prctile(data,[25,75]));
            binLength = 2*interQuartileDist*length(data)^(-1/3)*factor;

            % number of bins: divide data range by binLength
            nBins = round((max(data)-min(data))/binLength);

            if ~isfinite(nBins)
                nBins = length(unique(data));
            end

        end

        % histogram
        [nn,xx] = hist(data,nBins);
        % adjust the height of the histogram
        if normalize
            Z = trapz(xx,nn);
            nn = nn * nData/Z;
        end

    end
    if nargout > 0
        N = nn;
        X = xx;
        doPlot = axesHandle;
    end
    if doPlot
        if axesHandle
            bar(axesHandle,xx,nn,1);
        else
            bar(xx,nn,1);
        end
    end

else
    % make cdf, smooth with spline, then take the derivative of the spline

    % cdf
    xData = sort(data);
    yData = 1:nData;

    % when using too many data points, the spline fits very locally, and
    % the derivatives can still be huge. Good results can be obtained with
    % 500-1000 points. Use 1000 for now
    step = max(floor(nData/1000),1);
    xData2 = xData(1:step:end);
    yData2 = yData(1:step:end);

    % spline. Use strong smoothing
    cdfSpline = csaps(xData2,yData2,1./(1+mean(diff(xData2))^3/0.0006));

    % pdf is the derivative of the cdf
    pdfSpline = fnder(cdfSpline);

    % histogram
    if nargout > 0
        xDataU = unique(xData);
        N = fnval(pdfSpline,xDataU);
        X = xDataU;
        % adjust the height of the histogram
        Z = trapz(X,N);
        N = N * nData/Z;
        sp = pdfSpline;
        % set doPlot. If there is an axesHandle, we will plot
        doPlot = axesHandle;
    end
    % check if we have to plot. If we assigned an output, there will only
    % be plotting if there is an axesHandle.
    if doPlot
        if axesHandle
            plot(axesHandle,xData,fnval(pdfSpline,xData));
        else
            plot(xData,fnval(pdfSpline,xData));
        end
    end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function KEGG = createUniversalReactionModel(KEGGFilename, KEGGBlackList)
% Creates the `U` matrix using the universal data from the KEGG database
%
% USAGE:
%
%    KEGG = createUniversalReactionModel(KEGGFilename, KEGGBlackList)
%
% INPUTS:
%    KEGGFilename:     downloaded from KEGG database (ie. 'reaction.lst')
%    KEGGBlackList:    KEGG reactions not to use
%
% OUTPUT:
%    KEGG:             `U` Matrix
%
% .. Author: - 11-10-07 IT

if ~exist('KEGGFilename','var') || isempty(KEGGFilename)
    KEGGFilename='11-20-08-KEGG-reaction.lst';
end
if ~exist('KEGGBlackList','var') || isempty(KEGGBlackList)
    KEGGBlackList = {};
end

KEGGReactionList = importdata(KEGGFilename);
KEGG = createModel;
cnt=1;
cnti=1;

HTABLE = java.util.Hashtable; % hashes Kegg.mets
showprogress(0,'KEGG reaction list ...');
for i = 1: length(KEGGReactionList)
    clear Rxn rxnFormulas;
    [Rxn, rxnFormulas] = strtok(KEGGReactionList(i),':');
     %continue if reaction is not in KEGGBlacklist

    if isempty(strmatch(Rxn, KEGGBlackList, 'exact'))%length(strmatch(Rxn,KEGGBlackList,'exact'))==0

        KEGG.rxns(cnti,1)=Rxn;

        %reformats syntax of reaction
        rxnFormulas= regexprep(rxnFormulas,': ','');
        rxnFormulas= regexprep(rxnFormulas,'\+ C','\+ 1 C');
        rxnFormulas= regexprep(rxnFormulas,' \+','[c] \+');
        rxnFormulas= regexprep(rxnFormulas,'=> C','=> 1 C');
        rxnFormulas= regexprep(rxnFormulas,' <','[c] <');
        rxnFormulas= regexprep(rxnFormulas,'^(C)','1 C');
        rxnFormulas= regexprep(rxnFormulas,' \[c]','[c]');

        rxnFormulas= regexprep(rxnFormulas,'\+ G','\+ 1 G');
        rxnFormulas= regexprep(rxnFormulas,'=> G','=> 1 G');
        rxnFormulas= regexprep(rxnFormulas,'^(G)','1 G');
        rxnFormulas= regexprep(rxnFormulas,'^(n) ','2 ');
        rxnFormulas= regexprep(rxnFormulas,'\+ n ','\+ 2 ');
        rxnFormulas= regexprep(rxnFormulas,'\> n ','\> 2 ');
        rxnFormulas= regexprep(rxnFormulas,'\ n-1 ','\ 1 ');
        rxnFormulas= regexprep(rxnFormulas,'^(n-1) ','1 ');
        rxnFormulas= regexprep(rxnFormulas,'\+ 2n ','\+ 2 ');
        rxnFormulas= regexprep(rxnFormulas,'\+ 4n ','\+ 4 ');
        rxnFormulas= regexprep(rxnFormulas,'\+ (n+1) ','\+ 3 ');
        rxnFormulas= regexprep(rxnFormulas,'\=> (n+1) ','\=> 3 ');
        rxnFormulas= regexprep(rxnFormulas,' 3C',' 3 C');
        rxnFormulas= regexprep(rxnFormulas,' 2C',' 2 C');
        rxnFormulas= regexprep(rxnFormulas,' 4C',' 4 C');

        rxnFormulas = strcat(rxnFormulas,'[c]');
        rxnFormulas= regexprep(rxnFormulas,'<=>','<==>');
        rxnFormulas= regexprep(rxnFormulas,'\=>>','=>');

        KEGG.rxnFormulas(cnti,1)=rxnFormulas;
        cnti=cnti+1;
        %compounds is a list of each of metabolites involved in the
        %reaction that has a KEGGID starting with 'C'. b1 and c1 are just
        %so that the output has the correct number of variable
        [compounds, b1, c1] = regexp(char(rxnFormulas),'C\w+\[c]','match','start','end');
        for j=1:length(compounds)
            if (~isempty(compounds(j)))
                %condition1 = length(strmatch(compounds(j),KEGG.mets))==0
                condition2 = isempty(HTABLE.get(compounds{j}));
                %if condition1 ~= condition2, pause; end
                if (condition2)
                    HTABLE.put(compounds{j}, cnt);
                    KEGG.mets(cnt,1)=compounds(j);
                    cnt=cnt+1;
                end
            end
        end
        clear  compounds

        %compounds is a list of each of metabolites involved in the
        %reaction that has a KEGGID starting with 'G'
        [compounds, b1, c1] = regexp(char(rxnFormulas),'G\w+\[c]','match','start','end');
        for j=1:length(compounds)
            if (~isempty(compounds(j)))
                %condition1 = length(strmatch(compounds(j),KEGG.mets))==0
                condition2 = isempty(HTABLE.get(compounds{j}));
                %if condition1 ~= condition2, pause; end
                if (condition2)
                    HTABLE.put(compounds{j}, cnt);
                    KEGG.mets(cnt,1)=compounds(j);
                    cnt=cnt+1;
                end
            end
        end
    end
    showprogress(i/length(KEGGReactionList))
end

%Allocate a nMets x nRxns matrix with approximately 4.5 metabolites per
%reaction. (Recon2 has 4.23 per reaction)
KEGG.S=spalloc(length(KEGG.mets), length(KEGG.rxns), floor(4.5 * length(KEGG.rxns)));
%Initialize required fields properly.
KEGG.csense = repmat('E',size(KEGG.mets));
KEGG.b = zeros(size(KEGG.mets));
[KEGG] = addReactionGEM(KEGG,KEGG.rxns,KEGG.rxns,KEGG.rxnFormulas,ones(length(KEGG.rxns),1),-10000*ones(length(KEGG.rxns),1),10000*ones(length(KEGG.rxns),1),1);
%a=length(KEGG.mets);
%b=length(KEGG.rxns);
%KEGG.S(a+1:end,:)=[];
%KEGG.S(:,b+1:end) = [];

author: opencobra

repo: cobratoolbox


text:
% EX_MAXWELL_SRC_RING: solve Maxwell source problem in one quarter of a ring.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_ring.txt';

% Type of boundary conditions
problem_data.nmnn_sides   = [];
problem_data.drchlt_sides = [1 2 3 4];

% Physical parameters
problem_data.c_stiff = @(x, y) ones(size(x));
problem_data.c_mass  = @(x, y) ones(size(x));

% Source and boundary terms
problem_data.f = @(x, y) cat(1, ...
             reshape (2*sin(y), [1, size(x)]), ...
             reshape (2*sin(x), [1, size(x)]));
problem_data.g = @(x, y, ind) test_maxwell_ring_g_nmnn (x, y, ind);
problem_data.h = @(x, y, ind) cat(1, ...
                    reshape (sin(y), [1, size(x)]), ...
                    reshape (sin(x), [1, size(x)]));

% Exact solution (optional)
problem_data.uex     = @(x, y) cat(1, ...
                    reshape (sin(y), [1, size(x)]), ...
                    reshape (sin(x), [1, size(x)]));
problem_data.curluex = @(x, y) cos(x) - cos(y);

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [3 3];     % Degree of the bsplines
method_data.regularity = [2 2];     % Regularity of the splines
method_data.nsub       = [9 9];     % Number of subdivisions
method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, u] = solve_maxwell_src (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) EXPORT TO PARAVIEW
output_file = 'maxwell_ring_Deg3_Reg2_Sub9';

vtk_pts = {linspace(0, 1, 20), linspace(0, 1, 20)};
fprintf ('The result is saved in the file %s \n \n', output_file);
sp_to_vtk (u, space, geometry, vtk_pts, output_file, 'u')

% 4.2) Plot in Matlab. Comparison with the exact solution
[eu, F] = sp_eval (u, space, geometry, vtk_pts);
[X, Y]  = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)));
eu2     = problem_data.uex (X, Y);

subplot(1,2,1),
quiver (X, Y, squeeze(eu(1,:,:)), squeeze(eu(2,:,:)))
axis equal tight
title('Computed solution')
subplot(1,2,2)
quiver (X, Y, squeeze(eu2(1,:,:)), squeeze(eu2(2,:,:)))
axis equal tight
title('Exact solution')

[error_hcurl, error_l2] = ...
    sp_hcurl_error (space, msh, u, problem_data.uex, problem_data.curluex)

%!demo
%! ex_maxwell_src_ring

%!test
%! problem_data.geo_name = 'geo_ring.txt';
%! problem_data.nmnn_sides   = [];
%! problem_data.drchlt_sides = [1 2 3 4];
%! problem_data.c_stiff = @(x, y) ones(size(x));
%! problem_data.c_mass  = @(x, y) ones(size(x));
%! problem_data.f = @(x, y) cat(1, ...
%!              reshape (2*sin(y), [1, size(x)]), ...
%!              reshape (2*sin(x), [1, size(x)]));
%! problem_data.g = @(x, y, ind) test_maxwell_ring_g_nmnn (x, y, ind);
%! problem_data.h = @(x, y, ind) cat(1, ...
%!                     reshape (sin(y), [1, size(x)]), ...
%!                     reshape (sin(x), [1, size(x)]));
%! problem_data.uex     = @(x, y) cat(1, ...
%!                     reshape (sin(y), [1, size(x)]), ...
%!                     reshape (sin(x), [1, size(x)]));
%! problem_data.curluex = @(x, y) cos(x) - cos(y);
%! method_data.degree     = [3 3];     % Degree of the bsplines
%! method_data.regularity = [2 2];     % Regularity of the splines
%! method_data.nsub       = [9 9];     % Number of subdivisions
%! method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule
%! [geometry, msh, space, u] = solve_maxwell_src (problem_data, method_data);
%! [error_hcurl, error_l2] = sp_hcurl_error (space, msh, u, problem_data.uex, problem_data.curluex);
%! assert (msh.nel, 81)
%! assert (space.ndof, 264)
%! assert (error_l2, 3.85922857182716e-04, 1e-16)
%! assert (error_hcurl, 5.72616592482421e-04, 1e-16)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function G = gradient8(DEM,unit,varargin)

%GRADIENT8 8-connected neighborhood gradient of a digital elevation model
%
% Syntax
%
%     G = gradient8(DEM)
%     G = gradient8(DEM,unit)
%     G = gradient8(DEM,unit,pn,pv,...)
%
% Description
%
%     gradient8 returns the numerical steepest downward gradient and aspect
%     of a digital elevation model using an 8-connected neighborhood.
%
% Input
%
%     DEM       digital elevation model (class: GRIDobj)
%     unit      'tan' --> tangent (default)
%               'rad' --> radian
%               'deg' --> degree
%               'sin' --> sine
%               'per' --> percent
%
%     Parameter name value/pairs (pn,pv,...)
%
%     'useblockproc'    true or {false}: use block processing
%                       (see function blockproc)
%     'useparallel'     true or {false}: use parallel computing toolbox
%     'blocksize'       blocksize for blockproc (default: 5000)
%
% Output
%
%     G         gradient (class: GRIDobj)
%
% Example
%
%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
%     G = gradient8(DEM,'degree');
%     subplot(2,1,1)
%     imagesc(DEM)
%     subplot(2,1,2)
%     imagesc(G)
%
%
% See also: GRIDobj, GRIDobj/CURVATURE, GRIDobj/ASPECT, GRIDobj/arcslope
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 18. August, 2017

if nargin == 1;
    unit = 'tangent';
else
    unit = validatestring(unit,{'tangent' 'degree' 'radian' 'percent' 'sine'},'gradient8','unit',2);
end

p = inputParser;
p.FunctionName = 'GRIDobj/gradient8';
addParamValue(p,'useblockproc',false,@(x) isscalar(x));
addParamValue(p,'blocksize',5000,@(x) isscalar(x));
addParamValue(p,'useparallel',false,@(x) isscalar(x));
parse(p,varargin{:});

% create a copy of the DEM instance
G = DEM;
c = class(DEM.Z);
switch c
    case 'double'
        G.Z = double.empty(0,0);
    otherwise
        G.Z = single.empty(0,0);
        c   = 'single';
end

% I found Large matrix support using blockproc inefficient for gradient8.
% Matrix dimensions have thus been increased to an out-of-range value to
% avoid calling blockproc.
% Large matrix support. Break calculations in chunks using blockproc.

if p.Results.useblockproc
    blksiz = bestblk(size(DEM.Z),p.Results.blocksize);
    c   = class(DEM.Z);

    switch c
        case {'double', 'single'}
            padval = inf;
        case 'logical'
            padval = true;
        otherwise
            padval = intmax(c);
    end
    cs  = G.cellsize;
    fun = @(x) steepestgradient(x,cs,c);
    G.Z = blockproc(DEM.Z,blksiz,fun,...
           'BorderSize',[1 1],...
           'Padmethod',padval,...
           'UseParallel',p.Results.useparallel);
else
    G.Z = steepestgradient(DEM.Z,G.cellsize,c);
end

G.name = 'gradient';
G.zunit = unit;

switch unit
    case 'tangent'
        % do nothing
    case 'degree'
        G.Z = atand(G.Z);
    case 'radian'
        G.Z = atan(G.Z);
    case 'sine'
        G.Z = sin(atan(G.Z));
    case 'percent'
        G.Z = G.Z*100;
end
end

function G = steepestgradient(z,cellsize,c)

if isstruct(z);
    z = z.data;
end

% check for nans;
I = isnan(z);
flagnan = any(I(:));
if flagnan
    z(I) = inf;
end

NEIGH = false(3);
% calculate along orthogonal neighbors
NEIGH(2,:) = true;
NEIGH(:,2) = true;

switch c
    case 'double'
        G = (z-imerode(z,NEIGH))/cellsize;
    case 'single'
        G = single(z-imerode(z,NEIGH))/cellsize;
end

% calculate along diagonal neighbors
NEIGH(:,:) = false;
NEIGH([1 5 9 3 7]) = true;

switch c
    case 'double'
        G = max(G,(z-imerode(z,NEIGH))/norm([cellsize cellsize]));
    case 'single'
        G = max(G,single(z-imerode(z,NEIGH))/single(norm([cellsize cellsize])));
end

if flagnan
    G(I) = nan;
end
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function [g1, g2, g3] = sdlfmaXsdlfmKernGradientBlock(lfmKern1, lfmKern2, ...
    t1, t2, kyy, kyv, kvy, kvv, i, j, generalConst, generalConstGrad, ...
    covGrad)

% SDLFMAXSDLFMKERNGRADIENTBLOCK Gradients of the parameters in block i,j
% FORMAT
% DESC computes the kernel parameters gradients for the SDLFM kernel
% function in the block specified at indeces i,j. It assumes the
% computation for functions that system 1 describes an acceleration and
% system 2 describes a position.
% ARG lfmKern1 : structure containing parameters for the system 1
% ARG lfmKern2 : structure containing parameters for the system 2
% ARG t1 : times at which the system 1 is evaluated
% ARG t2 : times at which the system 2 is evaluated
% ARG kyy : covariance for the initial conditions between position 1 and
% position 2 at block i,j
% ARG kyv : covariance for the initial conditions between position 1 and
% velocity 2 at block i,j
% ARG kvy : covariance for the initial conditions between velocity 1 and
% position 2 at block i,j
% ARG kvv : covariance for the initial conditions between velocity 1 and
% velocity 2 at block i,j
% ARG i : interval to be evaluated for system 1
% ARG j : interval to be evaluated for system 2
% ARG generalConstant : constants evaluated with sdlfmKernComputeConstant.m
% ARG generalConstGrad : derivatives of the constants computed with
% sdlfmKernGradientConstant.m
% ARG covGrad : partial derivatives of the objective function wrt portion
% of the kernel matrix in block i,j
% RETURN g1 : gradients of parameters for the system 1
% RETURN g2 : gradients of parameters for the system 2
% RETURN g3 : gradients of switching points
%
% COPYRIGHT : Mauricio A. Alvarez, 2010.

% KERN

if nargin<11
    j = i;
    generalConst = [];
end

g3 = [];

% Compute derivatives of the mean terms with respect to the parameters

[g1Mean, g2Mean, gsp1Mean, gsp2Mean] = sdlfmKernGradientMean(lfmKern1(1), ...
    lfmKern2(1), t1, t2, kyy, kyv, kvy, kvv, covGrad, {'sdlfma','sdlfm'}, ...
    {'sdlfmj', 'sdlfmv'});

if i==j
    [g1, g2, g3t] = sdlfmaXsdlfmKernGradientBlockIEJ(lfmKern1, lfmKern2, t1, ...
        t2, covGrad, g1Mean, g2Mean, gsp1Mean, gsp2Mean, {'lfma', 'lfm'}, ...
        {'lfmj', 'lfm'}, {'lfma', 'lfmv'});
    g3(i) = g3t;
else
    if i>j
        [g1, g2, g3] = sdlfmXsdlfmKernGradientBlockIGJ(lfmKern1, lfmKern2, ...
            t1, t2, i, j, generalConst, generalConstGrad, covGrad, g1Mean, ...
            g2Mean, gsp1Mean, gsp2Mean, 'sdlfma', 'sdlfmj', 'lfmXlfm', ...
            'lfmvXlfm', {'lfmvXlfm', 'lfmvXlfm'}, {'lfmaXlfm', 'lfmvXlfmv'});
    else
        [g1, g2, g3] = sdlfmaXsdlfmKernGradientBlockILJ(lfmKern1, lfmKern2, ...
            t1, t2, i, j, generalConst, generalConstGrad, covGrad, g1Mean, ...
            g2Mean, gsp1Mean, gsp2Mean, 'sdlfm', 'sdlfmv', 'lfmaXlfm', ...
            'lfmaXlfmv',{'lfmjXlfm', 'lfmaXlfmv'}, {'lfmjXlfmv', 'lfmaXlfma'});
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function demtrain(action);
%DEMTRAIN Demonstrate training of MLP network.
%
%	Description
%	DEMTRAIN brings up a simple GUI to show the training of an MLP
%	network on classification and regression problems.  The user should
%	load in a dataset (which should be in Netlab format: see  DATREAD),
%	select the output activation function, the  number of cycles and
%	hidden units and then train the network. The scaled conjugate
%	gradient algorithm is used. A graph shows the evolution of the error:
%	the value is shown  MAX(CEIL(ITERATIONS / 50), 5) cycles.
%
%	Once the network is trained, it is saved to the file MLPTRAIN.NET.
%	The results can then be viewed as a confusion matrix (for
%	classification problems) or a plot of output versus target (for
%	regression problems).
%
%	See also
%	CONFMAT, DATREAD, MLP, NETOPT, SCG
%

%	Copyright (c) Ian T Nabney (1996-2001)

% If run without parameters, initialise gui.
if nargin<1,
  action='initialise';
end;

% Global variable to reference GUI figure
global DEMTRAIN_FIG
% Global array to reference sub-figures for results plots
global DEMTRAIN_RES_FIGS
global NUM_DEMTRAIN_RES_FIGS

if strcmp(action,'initialise'),

  file = '';
  path = '.';

  % Create FIGURE
  fig = figure( ...
	'Name', 'Netlab Demo', ...
	'NumberTitle', 'off', ...
	'Menubar', 'none', ...
	'Color', [0.7529 0.7529 0.7529], ...
	'Visible', 'on');
  % Initialise the globals
  DEMTRAIN_FIG = fig;
  DEMTRAIN_RES_FIGS = 0;
  NUM_DEMTRAIN_RES_FIGS = 0;

  % Create GROUP for buttons
  uicontrol(fig, ...
	'Style', 'frame', ...
	'Units', 'normalized', ...
	'Position', [0.03 0.08 0.94 0.22], ...
	'BackgroundColor', [0.5 0.5 0.5]);

  % Create MAIN axis
  hMain = axes( ...
	'Units', 'normalized', ...
	'Position', [0.10 0.5 0.80 0.40], ...
	'XColor', [0 0 0], ...
	'YColor', [0 0 0], ...
	'Visible', 'on');

  % Create static text for FILENAME and PATH
  hFilename = uicontrol(fig, ...
	'Style', 'text', ...
	'Units', 'normalized', ...
	'BackgroundColor', [0.7529 0.7529 0.7529], ...
	'Position', [0.05 0.32 0.90 0.05], ...
	'HorizontalAlignment', 'center', ...
	'String', 'Please load data file.', ...
	'Visible', 'on');
  hPath = uicontrol(fig, ...
	'Style', 'text', ...
	'Units', 'normalized', ...
	'BackgroundColor', [0.7529 0.7529 0.7529], ...
	'Position', [0.05 0.37 0.90 0.05], ...
	'HorizontalAlignment', 'center', ...
	'String', '', ...
	'Visible', 'on');

  % Create NO OF HIDDEN UNITS slider and text
  hSliderText = uicontrol(fig, ...
	'Style', 'text', ...
	'BackgroundColor', [0.5 0.5 0.5], ...
	'Units', 'normalized', ...
	'Position', [0.27 0.12 0.17 0.04], ...
	'HorizontalAlignment', 'right', ...
	'String', 'Hidden Units: 5');
  hSlider = uicontrol(fig, ...
	'Style', 'slider', ...
	'Units', 'normalized', ...
	'Position', [0.45 0.12 0.26 0.04], ...
	'String', 'Slider', ...
	'Min', 1, 'Max', 25, ...
	'Value', 5, ...
	'Callback', 'demtrain slider_moved');

  % Create ITERATIONS slider and text
  hIterationsText = uicontrol(fig, ...
	'Style', 'text', ...
	'BackgroundColor', [0.5 0.5 0.5], ...
	'Units', 'normalized', ...
	'Position', [0.27 0.21 0.17 0.04], ...
	'HorizontalAlignment', 'right', ...
	'String', 'Iterations: 50');
  hIterations = uicontrol(fig, ...
	'Style', 'slider', ...
	'Units', 'normalized', ...
	'Position', [0.45 0.21 0.26 0.04], ...
	'String', 'Slider', ...
	'Min', 10, 'Max', 500, ...
	'Value', 50, ...
	'Callback', 'demtrain iterations_moved');

  % Create ACTIVATION FUNCTION popup and text
  uicontrol(fig, ...
	'Style', 'text', ...
	'BackgroundColor', [0.5 0.5 0.5], ...
	'Units', 'normalized', ...
	'Position', [0.05 0.20 0.20 0.04], ...
	'HorizontalAlignment', 'center', ...
	'String', 'Activation Function:');
  hPopup = uicontrol(fig, ...
	'Style', 'popup', ...
	'Units', 'normalized', ...
	'Position' , [0.05 0.10 0.20 0.08], ...
	'String', 'Linear|Logistic|Softmax', ...
	'Callback', '');

  % Create MENU
  hMenu1 = uimenu('Label', 'Load Data file...', 'Callback', '');
  uimenu(hMenu1, 'Label', 'Select training data file', ...
	'Callback', 'demtrain get_ip_file');
  hMenu2 = uimenu('Label', 'Show Results...', 'Callback', '');
  uimenu(hMenu2, 'Label', 'Show classification results', ...
	'Callback', 'demtrain classify');
  uimenu(hMenu2, 'Label', 'Show regression results', ...
	'Callback', 'demtrain predict');

  % Create START button
  hStart = uicontrol(fig, ...
	'Units', 'normalized', ...
	'Position' , [0.75 0.2 0.20 0.08], ...
	'String', 'Start Training', ...
	'Enable', 'off',...
	'Callback', 'demtrain start');

  % Create CLOSE button
  uicontrol(fig, ...
	'Units', 'normalized', ...
	'Position' , [0.75 0.1 0.20 0.08], ...
	'String', 'Close', ...
	'Callback', 'demtrain close');

  % Save handles of important UI objects
  hndlList = [hSlider hSliderText hFilename hPath hPopup ...
      hIterations hIterationsText hStart];
  set(fig, 'UserData', hndlList);
  % Hide window from command line
  set(fig, 'HandleVisibility', 'callback');

elseif strcmp(action, 'slider_moved'),

  % Slider has been moved.

  hndlList = get(gcf, 'UserData');
  hSlider = hndlList(1);
  hSliderText = hndlList(2);

  val = get(hSlider, 'Value');
  if rem(val, 1) < 0.5,  % Force up and down arrows to work!
	val = ceil(val);
  else
	val = floor(val);
  end;
  set(hSlider, 'Value', val);
  set(hSliderText, 'String', ['Hidden Units: ' int2str(val)]);

elseif strcmp(action, 'iterations_moved'),

  % Slider has been moved.

  hndlList = get(gcf, 'UserData');
  hSlider = hndlList(6);
  hSliderText = hndlList(7);

  val = get(hSlider, 'Value');
  set(hSliderText, 'String', ['Iterations: ' int2str(val)]);

elseif strcmp(action, 'get_ip_file'),

  % Get data file button pressed.

  hndlList = get(gcf, 'UserData');

  [file, path] = uigetfile('*.dat', 'Get Data File', 50, 50);

  if strcmp(file, '') | file == 0,
    set(hndlList(3), 'String', 'No data file loaded.');
    set(hndlList(4), 'String', '');
  else
    set(hndlList(3), 'String', file);
    set(hndlList(4), 'String', path);
  end;

  % Enable training button
  set(hndlList(8), 'Enable', 'on');

  set(gcf, 'UserData', hndlList);

elseif strcmp(action, 'start'),

  % Start training

  % Get handles of and values from UI objects
  hndlList = get(gcf, 'UserData');
  hSlider = hndlList(1); % 				No of hidden units
  hIterations = hndlList(6);
  iterations = get(hIterations, 'Value');

  hFilename = hndlList(3);	% 			Data file name
  filename = get(hFilename, 'String');

  hPath = hndlList(4);	% 				Data file path
  path = get(hPath, 'String');

  hPopup = hndlList(5);		% 			Activation function
  if get(hPopup, 'Value') == 1,
	act_fn = 'linear';
  elseif get(hPopup, 'Value') == 2,
	act_fn = 'logistic';
  else
	act_fn = 'softmax';
  end;
  nhidden = get(hSlider, 'Value');

  % Check data file exists
  if fopen([path '/' filename]) == -1,
	errordlg('Training data file has not been selected.', 'Error');
  else
	% Load data file
	[x,t,nin,nout,ndata] = datread([path filename]);

	% Call MLPTRAIN function repeatedly, while drawing training graph.
	figure(DEMTRAIN_FIG);
	hold on;

	title('Training - please wait.');

	% Create net and find initial error
	net = mlp(size(x, 2), nhidden, size(t, 2), act_fn);
	% Initialise network with inverse variance of 10
	net = mlpinit(net, 10);
	error = mlperr(net, x, t);
	% Work out reporting step: should be sufficiently big to let training
	% algorithm have a chance
	step = max(ceil(iterations / 50), 5);

	% Refresh and rescale axis.
	cla;
	max = error;
	min = max/10;
	set(gca, 'YScale', 'log');
	ylabel('log Error');
	xlabel('No. iterations');
	axis([0 iterations min max+1]);
	iold = 0;
	errold = error;
	% Plot circle to show error of last iteration
	% Setting erase mode to none prevents screen flashing during
	% training
	plot(0, error, 'ro', 'EraseMode', 'none');
	hold on
	drawnow; % Force redraw
	for i = step-1:step:iterations,
	  [net, error] = mlptrain(net, x, t, step);
	  % Plot line from last point to new point.
	  line([iold i], [errold error], 'Color', 'r', 'EraseMode', 'none');
	  iold = i;
	  errold = error;

	  % If new point off scale, redraw axes.
	  if error > max,
	    max = error;
	    axis([0 iterations min max+1]);
	  end;
	  if error < min
	    min = error/10;
	    axis([0 iterations min max+1]);
	  end
	  % Plot circle to show error of last iteration
	  plot(i, error, 'ro', 'EraseMode', 'none');
	  drawnow; % Force redraw
	end;
	save mlptrain.net net
	zoom on;

	title(['Training complete. Final error=', num2str(error)]);

  end;

elseif strcmp(action, 'close'),

  % Close all the figures we have created
  close(DEMTRAIN_FIG);
  for n = 1:NUM_DEMTRAIN_RES_FIGS
    if ishandle(DEMTRAIN_RES_FIGS(n))
      close(DEMTRAIN_RES_FIGS(n));
    end
  end

elseif strcmp(action, 'classify'),

  if fopen('mlptrain.net') == -1,
	errordlg('You have not yet trained the network.', 'Error');
  else

	hndlList = get(gcf, 'UserData');
	filename = get(hndlList(3), 'String');
	path = get(hndlList(4), 'String');
	[x,t,nin,nout,ndata] = datread([path filename]);
	load mlptrain.net net -mat
	y = mlpfwd(net, x);

	% Save results figure so that it can be closed later
	NUM_DEMTRAIN_RES_FIGS = NUM_DEMTRAIN_RES_FIGS + 1;
	DEMTRAIN_RES_FIGS(NUM_DEMTRAIN_RES_FIGS)=conffig(y,t);

  end;

elseif strcmp(action, 'predict'),

  if fopen('mlptrain.net') == -1,
	errordlg('You have not yet trained the network.', 'Error');
  else

	hndlList = get(gcf, 'UserData');
	filename = get(hndlList(3), 'String');
	path = get(hndlList(4), 'String');
	[x,t,nin,nout,ndata] = datread([path filename]);
	load mlptrain.net net -mat
	y = mlpfwd(net, x);

	for i = 1:size(y,2),
	  % Save results figure so that it can be closed later
	  NUM_DEMTRAIN_RES_FIGS = NUM_DEMTRAIN_RES_FIGS + 1;
	  DEMTRAIN_RES_FIGS(NUM_DEMTRAIN_RES_FIGS) = figure;
	  hold on;
	  title(['Output no ' num2str(i)]);
	  plot([0 1], [0 1], 'r:');
	  plot(y(:,i),t(:,i), 'o');
	  hold off;
	end;
  end;

end;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function determ = householder_determinant ( n, x )

%*****************************************************************************80
%
%% HOUSEHOLDER_DETERMINANT returns the determinant of a HOUSEHOLDER matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 May 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X(N), the vector that defines the
%    Householder matrix.
%
%    Output, real DETERM, the determinant.
%
  determ = -1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
%
% This file is part of HRVAS
%
% HRVAS 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.
%
% HRVAS 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 HRVAS.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function output=poincareHRV(ibi)
%poincareHRV(ibi) - calculates poincare HRV
%
%Inputs:    ibi = 2dim array containing [t (s),ibi (s)]
%
%Outputs:   output is a structure containg HRV.

    %check inputs
    ibi(:,2)=ibi(:,2).*1000; %convert ibi to ms
    %assumes ibi units are seconds

%     if abs(range(ibi(:,2)))<50 %assume ibi units are seconds
%             ibi(:,2)=ibi(:,2).*1000; %convert ibi to ms
%     end
%     if abs(range(diff(ibi(:,1))))>50 %assume time unites are ms
%         ibi(:,1)=ibi(:,1)./1000; %convert time to s
%     end

    sd=diff(ibi(:,2)); %successive differences
    rr=ibi(:,2);
    SD1=sqrt( 0.5*std(sd)^2 );
    SD2=sqrt( 2*(std(rr)^2) - (0.5*std(sd)^2) );

    %format decimal places
    output.SD1=round(SD1*10)/10; %ms
    output.SD2=round(SD2*10)/10; %ms

end

author: jramshur

repo: HRVAS

save_path: github-repos/MATLAB/jramshur-HRVAS


text:
function [ grid_level, grid_order, grid_point ] = cc_levels_minmax ( ...
  dim_num, level_min, level_max, grid_num, point_num )

%*****************************************************************************80
%
%% CC_LEVELS_MINMAX computes CC grids with LEVEL_MIN <= LEVEL <= LEVEL_MAX.
%
%  Discussion:
%
%    The CC grids are required to have an order that is 2**LEVEL + 1.
%
%    The necessary dimensions of GRID_LEVEL, GRID_ORDER and GRID_POINT can be
%    determined by calling CC_LEVELS_MINMAX first.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MIN, LEVEL_MAX, the minimum and maximum values of
%    LEVEL.
%
%    Input, integer GRID_NUM, the number of Clenshaw Curtis
%    grids whose LEVEL value is between LEVEL_MIN and LEVEL_MAX.
%
%    Input, integer POINT_NUM, the total number of points in the grids.
%
%    Output, integer GRID_LEVEL(DIM_NUM,GRID_NUM), contains, for each
%    grid, the level of the Clenshaw-Curtis rule in each dimension.
%
%    Output, integer GRID_ORDER(DIM_NUM,GRID_NUM), contains, for each
%    grid, the order of the Clenshaw-Curtis rule in each dimension.
%
%    Output, real GRID_POINT(DIM_NUM,POINT_NUM), contains
%    a list of all the abscissas of all the rules, listed one grid at
%    a time.  If a point occurs in several grids, it will be listed
%    several times.
%

%
%  Outer loop generates LEVEL's from LEVEL_MIN to LEVEL_MAX.
%
  point_num = 0;
  grid_num = 0;

  for level = level_min : level_max
%
%  Middle loop generates next partition that adds up to LEVEL.
%
    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Inner (hidden) loop generates all CC points corresponding to given grid.
%
      order_1d = cc_level_to_order ( dim_num, level_1d );

      order_nd = prod ( order_1d(1:dim_num) );

      grid_point(1:dim_num,point_num+1:point_num+order_nd) = cc_grid ( ...
        dim_num, order_1d, order_nd );

      point_num = point_num + order_nd;

      grid_num = grid_num + 1;
      grid_level(1:dim_num,grid_num) = level_1d(1:dim_num);
      grid_order(1:dim_num,grid_num) = order_1d(1:dim_num);

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
function title = p09_title ( )

%*****************************************************************************80
%
%% P09_TITLE returns the title of problem 09.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Sunspot data, 1700-1960.';

  return
end


text:
clear
% define the root name of database
root = '../data_preparation/prepared_data/';

% which scales we're doing
sigma = 1;
num_samples = 2e6;

scales = [0.25,0.35,0.5,1.0];
frontalView = 1;

profileViewInds = [2];

version = 'wild';
ratio_neg = 10;
norm = 1;

data_loc = 'wild_';
rng(0);

% where to save generated patches for external training
% (e.g. for training CEN patches)
patches_loc = './patches/';
patch_folder = [patches_loc version '/'];

for s=scales
    Save_all_patches(root, frontalView, profileViewInds,...
        s, sigma, version, patch_folder, 'ratio_neg', ratio_neg,...
        'num_samples', num_samples, 'data_loc', data_loc,...
        'normalisation_size', 19);
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function y = gt(X,Y)
%GT (overloaded)

if isa(X,'blkvar')
    X = sdpvar(X);
end

if isa(Y,'blkvar')
    Y = sdpvar(Y);
end

try
    y = constraint(X,'>',Y);
catch
    error(lasterr)
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function actdistr()
    % activate deviation input editors in random parameters window
    %
    % This code is responsable for the activation (enable, on) of the std. deviation input
    % editors in the random parameters input window. It is called from randomcat.m,
    % Francesco Pacchiani 6/2000
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    if distr == 4

        set(input9, 'enable', 'on');
        set(input10, 'enable', 'on');
        set(input11, 'enable', 'on');
        set(tx9, 'color', 'k');
        set(tx10, 'color', 'k');
        set(tx11, 'color', 'k');

    elseif distr == 6

        set(input3, 'enable', 'off');
        set(input4, 'enable', 'off');
        set(input5, 'enable', 'off');
        set(input6, 'enable', 'off');
        set(input7, 'enable', 'off');
        set(input8, 'enable', 'off');
        set(input9, 'enable', 'off');
        set(input10, 'enable', 'off');
        set(input11, 'enable', 'off');
        set(tx3, 'color', 'w');
        set(tx4, 'color', 'w');
        set(tx5, 'color', 'w');
        set(tx6, 'color', 'w');
        set(tx7, 'color', 'w');
        set(tx8, 'color', 'w');
        set(tx9, 'color', 'w');
        set(tx10, 'color', 'w');
        set(tx11, 'color', 'w');

    elseif distr == 1  ||  distr == 2  ||  distr == 3  ||  distr == 5

        set(input3, 'enable', 'on');
        set(input4, 'enable', 'on');
        set(input5, 'enable', 'on');
        set(input6, 'enable', 'on');
        set(input7, 'enable', 'on');
        set(input8, 'enable', 'on');
        set(input9, 'enable', 'off');
        set(input10, 'enable', 'off');
        set(input11, 'enable', 'off');
        set(tx3, 'color', 'k');
        set(tx4, 'color', 'k');
        set(tx5, 'color', 'k');
        set(tx6, 'color', 'k');
        set(tx7, 'color', 'k');
        set(tx8, 'color', 'k');
        set(tx9, 'color', 'w');
        set(tx10, 'color', 'w');
        set(tx11, 'color', 'w');

    end

end

author: CelsoReyes


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%
classdef Cost < handle
  properties
    value
    userdata_p
  end

  methods

    function self = Cost(userdata)
      self.value = 0;
      self.userdata_p = userdata;
    end

    function r = userdata(self)
      r = self.userdata_p;
    end

    function add(self,val)
      % add(self,val)
      self.value = self.value + ocl.Variable.getValueAsColumn(val);
    end
  end
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function plotPredictionResults_LGG(pathResults,nameOutcome,fSetNames,metric,maxOrder,pathFig)
% -------------------------------------------------------------------------
% function plotPredictionResults_LGG(pathResults,nameOutcome,fSetNames,metrics,maxOrder)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function plots prediction performance estimation results for all the
% different feature set types entered as inputs in the LGG study.
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathResults: Full path to the 'RESULTS' folder where prediction results
%                 are saved.
%                 --> Ex: '/myProject/WORKSPACE/LOGISTIC_REGRESSION/RESULTS'
% 2. nameOutcome: String specifying the name of the outcome being displayed
%                 --> Ex: 'progression'
% 3. fSetNames: Cell of strings specifying the name of the type of feature
%               set analyzed.
%               --> Ex: {'T1W_T2W','T1W_T2F','T1CE_T2W','T1CE_T2F'}
% 4. metric: String specifying the metric to display.
%            --> 'AUC632'
% 5. maxOrder: Integer specifying the maximal multivariable model order.
%              --> Ex: 10
% 6. pathFig: (optional).  Full path to where figure is saved without
%             displaying it. Put '' for displaying the figure and not
%             saving it to 'pathFig' (default).
%             --> Ex: ''
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  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/>.
% -------------------------------------------------------------------------

if nargin < 6
    pathFig = '';
end

startpath = pwd;
cd(pathResults)

signs = {'-r',':b','--g','-.y'};
nFSET = numel(fSetNames);

maxOrderChosen = maxOrder;

if isempty(pathFig)
    figure
else
    h = figure('visible','off');
end
for i = 1:nFSET
    fSET = fSetNames{i};
    results = load(['RESULTS_',fSET,'_',nameOutcome]); results = struct2cell(results); results = results{1};
    nOrders = numel(fieldnames(results));
    if nOrders < maxOrderChosen
        maxOrder = nOrders;
    else
        maxOrder = maxOrderChosen;
    end
    val = zeros(maxOrder,1);
    val_SE = zeros(maxOrder,1);
    for j = 1:maxOrder
        orderName = ['Order',num2str(j)];
        val(j,1) = results.(orderName).(metric);
        val_SE(j,1) = results.(orderName).(['SE_',metric]);
    end
    errorbar(1:maxOrder,val(:,1),val_SE(:,1),signs{i},'LineWidth',3,'MarkerFaceColor',signs{i}(end),'MarkerSize',6)
    hold on
end
set(gca,'FontSize',20)
xlabel('Model Order','FontSize',24)
ylabel('Prediction performance','FontSize',24)
ind = strfind(metric,'632');
if ~isempty(ind)
    metric = [metric,'+'];
end
metric = [metric(1:ind-1),'_{',metric(ind:end),'}'];
ind = strfind(nameOutcome,'Death');
if ~isempty(ind)
    nameOutcome(ind:ind+4) = [];
    nameOutcome = ['Survival',nameOutcome];
end
titleName = [nameOutcome,' -- ',metric];
title(titleName,'FontSize',30,'FontWeight','bold')
legend(fSetNames,'Location','SouthEast')
axis([0 maxOrderChosen+1 0.5 1])
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10])

if ~isempty(pathFig)
    cd(pathFig)
    saveas(h,titleName,'fig')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function ApplyQiB()
%A rountine to execute the QIB on a computer cluster
% Written By: Issam El Naqa    Date: 10/3/03
% Revised by:                  Date:
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

% read a treatment plan in "mat" format
[filename, pathname] = uigetfile('*.mat', 'Load a treatment plan');
load([pathname,filename]);
D=1; % downsample factor (optional)
indexS = planC{end};
scaleX=planC{indexS.scan}.scanInfo(1).grid1Units*D;
scaleY=planC{indexS.scan}.scanInfo(1).grid2Units*D;
scaleZ=planC{indexS.scan}.scanInfo(1).sliceThickness;
% get the CT file
inputCT=planC{indexS.scan}.scanArray;
% apply down sampling and cropping
vecx=145:D:394; vecy=126:D:375; vecz=40:99;
inputCTc=inputCT(vecx,vecy,vecz);  clear inputCT
CTdim=size(inputCTc);
% add margins of 1 cm to the targets
Margin=1;
TargetStructNum=[17,18];
inputMask1 = getMask3D(TargetStructNum(1),planC);
inputMask2 = getMask3D(TargetStructNum(2),planC);
inputMask=inputMask1 | inputMask2;
inputMasks=inputMask(vecx,vecy,vecz); % just cropping
clear planC inputMask1 inputMask2 inputMask
b=genball(7,Margin/2,[scaleX,scaleY,scaleZ]); % plan specific
% treat the targets separately!
inputMaskb=convn(inputMasks,b,'same');
einputMask=uint8(inputMaskb>0);
ind=regexp(filename,'.mat');
maskfile=filename(1:ind-1);
%save([pathname,maskfile,'_targetmask.mat'],'inputMasks','einputMask','CTdim','scaleX','scaleY','scaleZ');
gridX_distance100=1; gridY_distance100=1;
Energy=6; % plan specific and is inversely proportional to the squared uncertainity
Nbeams=9; % plan specific
Source_Distance=100;
[x,y]=circle([0,0], Source_Distance,Nbeams);
SCRCM=[x(1:Nbeams);y(1:Nbeams);zeros(1,Nbeams)];
% not clear about conventions?!
for N=1:Nbeams
    beam_num=N
    % need to be worked out?!
    [targetBeamGridCoor, numberPB(N),OSC]=mybeamgridcoordinatesVMC(einputMask,scaleX,scaleY,scaleZ,gridX_distance100,gridY_distance100,SCRCM(:,N)); % this is wrong!
    for i=1:numberPB(N)
        [DoseF, downCT, vCT, depthCor] = get_PB_3D(Energy,gridX_distance100,gridX_distance100,Source_Distance,angle,shift,planC);
        %temporal saving! need to be changed to other format?!
        str=[filename,'_beam',num2str(N),'_pb',num2str(i)];
        fid=fopen(str,'w');
        fwrite(fid,CTdim(1),'int32'); fwrite(fid,CTdim(2),'int32');  fwrite(fid,CTdim(3),'int32');
        fwrite(fid,downCT,'float32'); fwrite(fid,vCT,'float32');  fwrite(fid,depthCor,'float32');
        fwrite(fid,DoseF,'float32');
    end
end

function b=genball(N,r,Scale)
% make a ball of dimension NxNxN and radius r, scaled by Scale
Lc=floor(N/2);
vec=-Lc:Lc;
[x,y,z]=meshgrid(vec*Scale(1),vec*Scale(2),vec*Scale(3));
r2=x.^2+y.^2+z.^2;
b=zeros(N,N,N);
b(find(r2<r^2))=1;
return

function [x,y]=circle(c, r, nsides)
% a poly circle routine with center c, radius r, and nsides
nsides = round(nsides);  % make sure it is an integer
a = [0:2*pi/nsides:2*pi];
x=r*cos(a)+c(1);
y=r*sin(a)+c(2);
%line(x,y); uncomment for plotting
return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [err1, err2] = llGrad(p1,p2,type,tolGrad,tolEval)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% [err1,err2] = llGrad(p1, p2, type [, tolGrad ,tolEval])
%
%    Compute the gradient with respect to the means of p1,p2 of the
%       average log-likelihood of p1 evaluated at the locations of p2
%       (p1=p2 => gradient of a resubstitution estimate of negentropy)
%    type: what to take gradient with respect to.
%          0 = means, 1 = bw (variance), 2=weights
%    tolGrad is an acceptable error tolerance on the gradient values
%    tolEval is a *percent* tolerance on the evaluation done as a
%       subroutine in the calculation; it should be tight enough to allow
%       tolGrad to be achieved but beyond that, larger => faster.
%
% Specifically, using the resubstition estimate
%    \hat H = \sum_j v_j \log \sum_i w_i K(x_i - y_j)
% The outputs are
%    err1(:,i) = \sum_j v_j 1/p(y_j) dp(y_j)/dx_i = E_y[ 1/p dp/dx_i ]
%    err2(:,j) = v_j 1/p(y_j) p'(y_j)
%
% Note: [err1,err2] = llGrad(p1, type [, tolGrad ,tolEval]) computes the leave-
%   one-out resub. estimate, similar to llGrad(p1,p1,type...)
%
% Some useful estimates which can be made using this function:
%
% ENTROPY GRADIENT
%    [err1, err2] = llGrad(p,p,1e-3,1e-3);
%    p -= (err1 + err2)    moves p in the direction of increasing entropy
%
% KL-DIVERGENCE GRADIENT
%    [errXX1, errXX2] = llGrad(p1,p1,1e-3,1e-3);
%    [errXYY, errXYX] = llGrad(p2,p1,1e-3,1e-3);
%    err1 = (errXX1 + errXX2 - errXYX);
%    err2 = (-errXYY);
%    p1 += err1, p2 += err2  moves p1,p2 in the direction of increasing KLDiv
%
% See also:  klGrad, miGrad, entropyGrad, adjustPoints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txt

%#mex
error('MEX-file kde/llGrad not found -- please recompile if necessary');

% Incorrect:  Mean shift corresponds to an estimate of 1/p(y) p'(y) where
%    p'(y) is computed using the Epan. kernel while p(y) is computed using
%    a uniform kernel with the same support.  This does *not* correspond to
%    the derivative of log-likelihood for any KDE (that I know of)
%
% MEAN-SHIFT
%    For equal weight uniform bandwidth Epan. kernel densities,
%      the mean shift algorithm (Comaniciu '99) is given by
%    [err1, err2] = llGrad(p,p,1e-3,1e-3);
%    p += .5*getNpts(p)*getBW(p,1)^2 * err2
%

author: ShapeNet

repo: RenderForCNN


text:
function sgmga_vcn_tests ( )

%*****************************************************************************80
%
%% SGMGA_VCN_TESTS tests SGMGA_VCN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2010
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../sandia_rules' );

  timestamp ( );

  test_num = 12;

  dim_num_array(1:test_num,1) = [ 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4 ]';
  level_max_array(1:test_num,1) = [ 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 2, 3 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_VCN_TESTS\n' );
  fprintf ( 1, '  calls SGMGA_VCN_TEST.\n' );
%
%  Isotropic examples.
%
  for test = 1 : test_num

    dim_num = dim_num_array(test);
    importance(1:dim_num,1) = 1.0;
    level_weight = sgmga_importance_to_aniso ( dim_num, importance );
    level_max = level_max_array(test);
    q_min = level_max - sum ( level_weight(1:dim_num) );
    q_max = level_max;

    sgmga_vcn_test ( dim_num, importance, level_weight, q_min, q_max );

  end
%
%  Zero weight example.
%
  dim_num = 3;
  importance = [ 1.0, 0.0, 1.0]';
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2
  q_min = level_max - sum ( level_weight(1:dim_num) );
  q_max = level_max;

  sgmga_vcn_test ( dim_num, importance, level_weight, q_min, q_max );
%
%  Anisotropic examples.
%
  for test = 1 : test_num

    dim_num = dim_num_array(test);
    importance(1:dim_num,1) = ( 1 : dim_num );
    level_weight = sgmga_importance_to_aniso ( dim_num, importance );
    level_max = level_max_array(test);
    q_min = level_max - sum ( level_weight(1:dim_num) );
    q_max = level_max;

    sgmga_vcn_test ( dim_num, importance, level_weight, q_min, q_max );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_VCN_TESTS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  rmpath ( '../sandia_rules' );

  return
end


text:
function [favar,favar_hd_estimates]=favar_hdestimates(favar,hd_estimates,n,IRFt,endo,strctident,L_g)

%% prelim
% number of identified shocks & create labels for the contributions
if IRFt==1||IRFt==2||IRFt==3
    identified=n; % fully identified
    labels=endo;
elseif IRFt==4 || IRFt==6 %if the model is identified by sign restrictions or sign restrictions (+ IV)
    identified=size(strctident.signreslabels_shocks,1); % count the labels provided in the sign res sheet (+ IV)
    labels=strctident.signreslabels_shocks; % signreslabels
elseif IRFt==5
    identified=1; % one IV shock
    labels{identified,1}=strcat('IV Shock (',strctident.Instrument,')'); % one label for the IV shock
end

% if favar.onestep==1
%
% else
%     L=favar.L(favar.plotX_index,:);
% end

L=L_g;

% if the number of identified shocks is smaller than the number endo
% variables (not fully identified) and we have blocks, (and sum the
% shock contributions) then we have to compute blocks_index_shocks, to
% know to which blocks the shocks are assigned to, this is identified
% analogue to the Factors by checking for 'Blockidentifier'. [favar.bnames{jj,1} '.']&& identified<n
% endo_exfactors=endo(favar.variablestrings_exfactors);
if favar.blocks==1
    for jj=1:favar.nbnames
        for ii=1:size(labels,1)
            %for ll=1:favar.bnumpc{jj,1}
            %blocks_indexlogical{jj,1}(ll,ii)=strcmp(endo{ii,1},favar.factorlabels_blocks{jj,ll});
            blocks_indexlogical_shocks{jj,1}(ii,1)=contains(labels{ii,1},[favar.bnames{jj,1} '.'])==1;%((ll,ii))
            %end
        end
        %blocks_indexlogical{jj,1}=sum(blocks_indexlogical{jj,1},1);
        favar.blocks_index_shocks{jj,1}=find(blocks_indexlogical_shocks{jj,1}==1);
    end
%     nbnames_shocks=size(favar.blocks_index_shocks,1)
end

% favar.HD.sumShockcontributions=0;

    for jj=1:favar.npltX %for selected information variables
        for ii=1:n %for variables
            for ll=1:n%size(hd_estimates,1) %for shock contributions
                favar_hd_estimates{ll,ii,jj}(1,:)=hd_estimates{ll,ii}(2,:)*L(jj,ii);%summed over variables {ll,ii,jj} %%% summed over contributions {ii,ll,jj}
            end
        end

        for ll=n+1:size(hd_estimates,1) %for the rest of the contributions
            for ii=1:n %for variables
                favar_hd_estimates{ll,ii,jj}(1,:)=hd_estimates{ll,ii}(2,:)*L(jj,ii);
            end
        end
    end

author: european-central-bank

repo: BEAR-toolbox


text:
function value = year_length_months_english ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_MONTHS_ENGLISH returns the number of months in an English year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of months in the year.
%
  value = 12;

  return
end


text:
function indx = multigrid_index_cfn ( dim_num, order_1d, order_nd )

%*****************************************************************************80
%
%% MULTIGRID_INDEX_CFN indexes a sparse grid based on CFN 1D rules.
%
%  Discussion:
%
%    The sparse grid is presumed to have been created from products
%    of CLOSED FULLY NESTED 1D quadrature rules.
%
%    CFN rules include Clenshaw Curtis rules.
%
%    For dimension DIM, the second index of INDX may vary from
%    0 to ORDER_1D(DIM)-1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension of the points.
%
%    Input, integer ORDER_1D(DIM_NUM), the order of the
%    rule in each dimension.
%
%    Input, integer ORDER_ND, the product of the entries of ORDER_1D.
%
%    Output, integer INDX(DIM_NUM,ORDER_ND), the indices of the points in
%    the grid.  The second dimension of this array is equal to the
%    product of the entries of ORDER_1D.
%
  p = 0;
  indx = zeros(dim_num,order_nd);

  a = [];
  more = 0;

  while ( 1 )

    [ a, more ] = vec_colex_next2 ( dim_num, order_1d, a, more );

    if ( ~more )
      break
    end

    p = p + 1;

    indx(1:dim_num,p) = a(1:dim_num)';

  end

  return
end


text:
function dt = datatypes
% Datatype
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: datatypes.m 1143 2008-02-07 19:33:33Z spm $

persistent dtype
if isempty(dtype),
    t = true;
    f = false;
    table = {...
        0   ,'UNKNOWN'   ,'uint8'   ,@uint8  ,1,1  ,t,t,f
        1   ,'BINARY'    ,'uint1'   ,@logical,1,1/8,t,t,f
        256 ,'INT8'      ,'int8'    ,@int8   ,1,1  ,t,f,t
        2   ,'UINT8'     ,'uint8'   ,@uint8  ,1,1  ,t,t,t
        4   ,'INT16'     ,'int16'   ,@int16  ,1,2  ,t,f,t
        512 ,'UINT16'    ,'uint16'  ,@uint16 ,1,2  ,t,t,t
        8   ,'INT32'     ,'int32'   ,@int32  ,1,4  ,t,f,t
        768 ,'UINT32'    ,'uint32'  ,@uint32 ,1,4  ,t,t,t
        1024,'INT64'     ,'int64'   ,@int64  ,1,8  ,t,f,f
        1280,'UINT64'    ,'uint64'  ,@uint64 ,1,8  ,t,t,f
        16  ,'FLOAT32'   ,'float32' ,@single ,1,4  ,f,f,t
        64  ,'FLOAT64'   ,'double'  ,@double ,1,8  ,f,f,t
        1536,'FLOAT128'  ,'float128',@crash  ,1,16 ,f,f,f
        32  ,'COMPLEX64' ,'float32' ,@single ,2,4  ,f,f,f
        1792,'COMPLEX128','double'  ,@double ,2,8  ,f,f,f
        2048,'COMPLEX256','float128',@crash  ,2,16 ,f,f,f
        128 ,'RGB24'     ,'uint8'   ,@uint8  ,3,1  ,t,t,f};
    dtype = struct(...
        'code'     ,table(:,1),...
        'label'    ,table(:,2),...
        'prec'     ,table(:,3),...
        'conv'     ,table(:,4),...
        'nelem'    ,table(:,5),...
        'size'     ,table(:,6),...
        'isint'    ,table(:,7),...
        'unsigned' ,table(:,8),...
        'min',-Inf,'max',Inf',...
        'supported',table(:,9));
    for i=1:length(dtype),
        if dtype(i).isint
            if dtype(i).unsigned
                dtype(i).min =  0;
                dtype(i).max =  2^(8*dtype(i).size)-1;
            else
                dtype(i).min = -2^(8*dtype(i).size-1);
                dtype(i).max =  2^(8*dtype(i).size-1)-1;
            end;
        end;
    end;
end;

dt = dtype;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function theta = initializeParameters(hiddenSize, visibleSize)

%% Initialize parameters randomly based on layer sizes.
r  = sqrt(6) / sqrt(hiddenSize+visibleSize+1);   % we'll choose weights uniformly from the interval [-r, r]
W1 = rand(hiddenSize, visibleSize) * 2 * r - r;
W2 = rand(visibleSize, hiddenSize) * 2 * r - r;

b1 = zeros(hiddenSize, 1);
b2 = zeros(visibleSize, 1);

% Convert weights and bias gradients to the vector form.
% This step will "unroll" (flatten and concatenate together) all
% your parameters into a vector, which can then be used with minFunc.
theta = [W1(:) ; W2(:) ; b1(:) ; b2(:)];

end

repo: MachineLearning


text:
%Function to integrate to calculate the area of a figure
function partieaire = functionviewfactorarea(t)

%Global variables: pt1 to pt2 are the points that define the segment to integrate, normale is the unit normal vector.
global pt1;
global pt2;
global normale;

%Parametric equations of the segment that join pt1 to pt2
x=pt1(1)+(pt2(1)-pt1(1)).*t;
y=pt1(2)+(pt2(2)-pt1(2)).*t;
z=pt1(3)+(pt2(3)-pt1(3)).*t;

%Function to integrate to calculate the area of a figure
partieaire=normale(2).*z.*(pt2(1)-pt1(1))+normale(3).*x.*(pt2(2)-pt1(2))+normale(1).*y.*(pt2(3)-pt1(3));

author: Sable

repo: mcbench-benchmarks


text:
function [x,g,xn,gg] = kmeanhar(d,k,l,e,x0)
%KMEANS Vector quantisation using K-harmonic means algorithm [X,G,XN,GG]=(D,K,L,E,X0)
%
%  Inputs:
%
%    D(N,P)  contains N data vectors of dimension P
%    K       is number of centres required
%    L       integer portion is max loop count, fractional portion
%            gives stopping threshold as fractional reduction in performance criterion
%    E       is exponent in the cost function. Significantly faster if this is an even integer. [default 4]
%    X0(K,P) are the initial centres (optional)
%            Alternatively, X0 can be a character determining the initialization method:
%                'f'    Initialize with K randomly selected data points [default]
%                'p'    Initialize with centroids and variances of random partitions
%
%  Outputs:
%
%    X(K,P)  is output row vectors
%    G       is the final performance criterion value (normalized by N)
%    XN      nearest centre for each input point
%    GG(L+1) value of performance criterion before each iteration and at end
%
% The k-harmonic means algorithm selects K cluster centres to minimize
%                           sum_n(K/sum_k((d_n-x_k)^-e))
% where sum_n is over the N inputs points d_n and sum_k is over the K cluster centres x_k.
%
% It is often a good idea to scale the input data so that it has equal variance in each
% dimension before calling KMEANHAR so that approximately equal weight is given
% to each dimension in the distance calculation.

%  [1] Bin Zhang, "Generalized K-Harmonic Means - Boosting in Unsupervised Learning",
%      Hewlett-Packartd Labs, Technical Report HPL-2000-137, 2000 [Zhang2000]
%      http://www.hpl.hp.com/techreports/2000/HPL-2000-137.pdf

%  Bugs:
%      (1) Could use nested blocking to allow very large data arrays
%      (2) Could then allow incremental calling with partial data arrays (but messy)

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: kmeanhar.m,v 1.6 2008/06/02 07:20:52 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% sort out the input arguments

if nargin<5
    x0='f';
    if nargin<4
        e=[];
        if nargin<3
            l=[];
        end
    end
end
if isempty(e)
    e=4;  % default value
end
if isempty(l)
    l=50+1e-3; % default value
end
sd=5;       % number of times we must be below threshold

% split into chunks if there are lots of data points

memsize=voicebox('memsize');
[n,p] = size(d);
nb=min(n,max(1,floor(memsize/(8*p*k))));    % block size for testing data points
nl=ceil(n/nb);                  % number of blocks

% initialize if X0 argument is not supplied

if ischar(x0)
    if k<n
        if any(x0=='p')                  % Initialize using a random partition
            ix=ceil(rand(1,n)*k);       % allocate to random clusters
            ix(rnsubset(k,n))=1:k;      % but force at least one point per cluster
            x=zeros(k,p);
            for i=1:k
                x(i,:)=mean(d(ix==i,:),1);
            end
        else                                % Forgy initialization: choose k random points [default]
            x=d(rnsubset(k,n),:);         % sample k centres without replacement
        end
    else
        x=d(mod((1:k)-1,n)+1,:);    % just include all points several times
    end
else
    x=x0;
end
eh=e/2;
th=l-floor(l);
l=floor(l)+(nargout>1);   % extra loop needed to calculate final performance value
if l<=0
    l=100;      % max number of iterations ever
end
if th==0
    th=-1;      % prevent any stopping if l has no fractional part
end
gg=zeros(l+1,1);
im=repmat(1:k,1,nb); im=im(:);

% index arrays for replication

wk=ones(k,1);
wp=ones(1,p);
% wn=ones(1,n);
%
% % Main calculation loop
%
% We have the following relationships to [1] where i and k index
% the data values and cluster centres respectively:
%
%   This program     [Zhang2000]                            Equation
%
%     d(i,:)            x_i                                 input data
%     x(k,:)            m_k                                 cluster centres
%     py(k,i)           (d_ik)^2
%     dm(i)'            d_i,min^2
%     pr(k,i)           (d_i,min/d_ik)^2
%     pe(k,i)           (d_i,min/d_ik)^p                    (7.6)
%     qik(k,i)          q_ik                                (7.2)
%     qk(k)             q_k                                 (7.3)
%     qik(k,i)./qk(k)   p_ik                                (7.4)
%     se(i)'            d_i,min^p * sumk(d_ik^-p)
%     xf(i)'            d_i,min^-2 / sumk(d_ik^-p)
%     xg(i)'            d_i,min^-(p+2) / sumk(d_ik^-p)^2

ss=sd+1;        % one extra loop at the start
g=0;                % dummy initial value of g
xn=zeros(n,1);
for j=1:l

    g1=g;                           % save old performance
    x1=x;                           % save old centres
    % first do partial chunk

    jx=n-(nl-1)*nb;
    ii=1:jx;
    kx=repmat(ii,k,1);
    km=repmat(1:k,1,jx);
    py=reshape(sum((d(kx(:),:)-x(km(:),:)).^2,2),k,jx);
    [dm,xn(ii)]=min(py,[],1);                 % min value in each column gives nearest centre
    dmk=dm(wk,:);                   % expand into a matrix
    dq=py>dmk;                      % update only these values
    pr=ones(k,jx);                   % leaving others at 1
    pr(dq)=dmk(dq)./py(dq);            % ratio of min(py)./py
    pe=pr.^eh;
    se=sum(pe,1);
    xf=dm.^(eh-1)./se;
    g=xf*dm.';                     % performance criterion (divided by k)
    xg=xf./se;
    qik=xg(wk,:).*pe.*pr;           % qik(k,i) is equal to q_ik in [Zhang2000]
    qk=sum(qik,2);
    xs=qik*d(ii,:);
    ix=jx+1;
    for il=2:nl
        jx=jx+nb;        % increment upper limit
        ii=ix:jx;
        kx=ii(wk,:);
        py=reshape(sum((d(kx(:),:)-x(im,:)).^2,2),k,nb);
        [dm,xn(ii)]=min(py,[],1);                 % min value in each column gives nearest centre
        dmk=dm(wk,:);                   % expand into a matrix
        dq=py>dmk;                      % update only these values
        pr=ones(k,nb);                   % leaving others at 1
        pr(dq)=dmk(dq)./py(dq);            % ratio of min(py)./py
        pe=pr.^eh;
        se=sum(pe,1);
        xf=dm.^(eh-1)./se;
        g=g+xf*dm.';                     % performance criterion (divided by k)
        xg=xf./se;
        qik=xg(wk,:).*pe.*pr;           % qik(k,i) is equal to q_ik in [Zhang2000]
        qk=qk+sum(qik,2);
        xs=xs+qik*d(ii,:);
        ix=jx+1;
    end
    gg(j)=g;
    x=xs./qk(:,wp);
    if g1-g<=th*g1
        ss=ss-1;
        if ~ss break; end  %  stop if improvement < threshold for sd consecutive iterations
    else
        ss=sd;
    end
end
gg=gg(1:j)*k/n;                       % scale and trim the performance criterion vector
g=g(end);
% gg' % *** DEBUIG ***
if nargout>1
    x=x1;                               % go back to the previous x values if G and/or XN value is output
end

author: decouples