text:
% Script to run a suite of test problems
% We have six test problems wherein the local Matlab client (possibly) uses
% the ithaca cluster with the Parallel Computing Toolbox

% The various test examples are in folders within the current directory

%%
clear all

local = pwd; addpath(local);

%% md_wrap
% md_program runs a 'molecular dynamics' simulation and writes some results
% to the the Command Window

 try
   fprintf(1,'\n \n Run the parallel MD example \n');
   cd('./john_md');
   nw = 4;
   results = md_wrap(nw);
% Show that we have retrieved results to the local workspace
   fprintf(1,'\n Final potential energy  %14f \n ',results{1}.pe(end))
   cd(local)
   clear all
   fprintf(1,'Results have been cleared \n')
 catch MD
     disp('Trouble with MD_WRAP')
     cd(local)
     MD %#ok
     MD.message
 end

%% systolic_ithaca
% systolic_ithaca creates and submits a job with a systolic array to filter
% an audio file. The job runs on ithaca.arc.vt.edu and requires a parallel
% configiguration file ithaca_2009b.
% The original and the filtered files are played on the local machine
% (assuming it has speakers).

 local = pwd;
 try
   fprintf(1,'\n \n Run the systolic array example \n');
   cd ./lyrics
   systolic_ithaca
   cd(local)
   clear all
 catch SYSTOLIC
     disp('Trouble with systolic_ithaca')
     cd(local)
     SYSTOLIC %#ok
 end

%% dist_knap
% dist_knap  creates a job with four tasks. Each task loops through
% a range of possible weight combinations to find those that sum to a given
% value

 local = pwd;
 try
   fprintf(1,'\n \n Run the knapsack example \n');
   cd ./knap
   dist_knap
   cd(local)
   clear all
 catch KNAP
     disp('Trouble with dist_knap')
     cd(local)
     KNAP %#ok
 end

%% color_segmentation
% colorSegmentationSol2  creates a matlabpool and uses an spmd construct
% to extract parts of an image with a given color. This version creates
% several images on the cluster - these are a Composite variable on the
% local client. After running the filter we write a jpg file on the client.

 local = pwd;
 try
   fprintf(1,'\n \n Run the color segmentation example \n');
   cd ./color_segmentation
   filteredImage=colorSegmentationSol2('ithaca');
   cd(local)
% Display some information about the result
   fprintf(1,'\n  The binary image files are referenced as a Composite variable on the client \n')
   fprintf(1,' We retrieve one of these (four) images to a jpeg file on the client \n');
   fprintf(1,' This could take awhile \n \n ');
%
    jj = 3; % pick one
%  for jj=1:length(filteredImage)
    f_Image = filteredImage{jj};
    disp('size of the image array  ');  size(f_Image)
    imwrite( f_Image, ['fab_' int2str(jj) '.jpg'], 'jpg')
%  end
   fprintf(1,'Results have been written to the local directory \n');
   fprintf(1,'Now close the matlabpool \n');
   matlabpool close
   clear all
 catch COLOR
     disp('Trouble with color_segmentation')
     if matlabpool('size') > 0
        matlabpool close force ithaca_2009b
     end
     cd(local)
     COLOR %#ok
     COLOR.message
 end

%% birthday
% brunbirthday  creates a batch job that uses a matlabpool and a parfor
% construct. The job runs a number of realizations (Monte Carlo) to
% generate an estimate for the probability that at least two people in a
% group (or 30) have the same birthday.

 local = pwd;
 try
   fprintf(1,'\n \n Run the birthday example \n');
   cd ./birthday
   nw = 4;
   prob = brunbirthday(nw)  %#ok
   cd(local)
   fprintf(1,'The brunbirthday script destroys the batch job \n');
   clear all
 catch BIRTH
     disp('Trouble with birthday')
     if matlabpool('size') > 0
        matlabpool close force ithaca_2009b
     end
     cd(local)
     BIRTH %#ok
     BIRTH.message
 end


text:
function [sMap, distNeuronsToData, sTrain] = rsom_batchtrain(sMap, Distances, varargin)

%RSOM_BATCHTRAIN Use a batch algorithm to train the Self-Organizing Map on
%relational data.
%
% [sM,distToData,sT] = som_batchtrain(sM, D, [argID, value, ...])
%
%  sM = som_batchtrain(sM, D);
%  [sM, distToData] = som_batchtrain(sM, D, 'epochs', 50, 'lambda0', 2.5);
%
% Input and output arguments:
%   sM         (struct) rsom map struct, the trained and updated map is
%                       returned
%   D          (matrix) dissimilarity data for training, size nData x nData
%   [argID,    (string) See below.
%    value]    (varies)
%
%   distToData (matrix) distance  from the neurons to data,
%                       size nNeurons x nData
%   sT         (struct) learning parameters used during the training
%
% Here are the valid argument IDs and corresponding values.
%   'lambda0'     (scalar) initial decay constant
%   'lambdaFin'   (scalar) final decay value
%   'epochs'      (scalar) number of epochs
%
% For more help, try 'type rsom_batchtrain' or check out online documentation.
% See also RSOM_RANDINIT, RSOM_LININIT.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% rsom_batchtrain
%
% PURPOSE
%
% Trains a Relational Self-Organizing Map on dissimilarity data using the
% batch algorithm.
%
% SYNTAX
%
%  sM = rsom_batchtrain(sM,D);
%  sM = rsom_batchtrain(...,'argID',value,...);
%
% DESCRIPTION
%
% Trains the given and already initialized RSOM with the batch algorithm.
% Training data are expected to be given in the form of dissimilarity data
% D. If D is generated from euclidean data, i.e.
% (D)_ij = ||x_i - x_j||^2, training the relational SOM is equivalent to
% the standard SOM and the neuron position in the input space can be
% acquired by sMap.cCodebook * x, where x are training data.
%
% REFERENCES
%
%   Barbara Hammer, Alexander Hasenfuss: Topographic Mapping of Large
%   Dissimilarity Data Sets. Neural Computation 22(9): 2229-2284 (2010)
%
% REQUIRED INPUT ARGUMENTS
%
%   sM         (struct) initialized rsom map struct
%   D          (matrix) dissimilarity data for training, size nData x nData
%
% OPTIONAL INPUT ARGUMENTS
%
%  argID (string) Argument identifier string (see below).
%  value (varies) Value for the argument (see below).
%
%  The optional arguments can be given as 'argID',value -pairs.
%  The valid IDs and corresponding values are listed below.
%
%  Below is the list of valid arguments:
%   'lambda0'     (scalar) initial decay constant
%   'lambdaFin'   (scalar) final decay value
%   'epochs'      (scalar) number of training epochs
%
% OUTPUT ARGUMENTS
%
%   sM         (struct) the trained rsom map struct. The current training
%                       is added to the training history (sM.trainhist).
%   distToData (matrix)
%   sT         (struct) train struct; information of the accomplished
%                       training
%
% EXAMPLES
%
%  mData = size(D,1);
%  sM = rsom_randinit(nData, [20 3]); % create a 20 x 3 grid
%  sM = rsom_batchtrain(sM,D);
%
% SEE ALSO
%
%   rsom_randinit   Initialize a RSOM randomly

% Contributed to SOM Toolbox vs2, December 7th, 2012 by Alexander Schulz
% Copyright (c) Alexander Hasenfuss
% http://www.cis.hut.fi/projects/somtoolbox/

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

%% Prase input
p = inputParser;
p.addRequired('sMap', @isstruct);
p.addRequired('Distances', @isnumeric);

% Get size
m = size( Distances, 1 );
GridDistances = som_unit_dists(sMap.topol).^2;
n = size( GridDistances, 1 );

p.addParamValue('lambda0', n/2, @isnumeric);
p.addParamValue('lambdaFin', 0.01, @isnumeric);
p.addParamValue('epochs', 100, @(x) ~(x-floor(x)) & (length(x) == 1));

p.parse(sMap, Distances, varargin{:});
lambda0   = p.Results.lambda0;
lambdaFin = p.Results.lambdaFin;
epochs    = p.Results.epochs;

%% Initialize
lambda = lambda0 * (lambdaFin/lambda0).^((0:(epochs-1))/(epochs-1));

QuantizationErrors = zeros(1,epochs);

cNeurons = sMap.cCodebook;

sTrain = som_train_struct('algorithm','batch');
sTrain = som_set(sTrain, 'neigh', sMap.neigh, 'trainlen', epochs, ...
    'radius_ini', sqrt(lambda0/2), 'radius_fin', sqrt(lambdaFin/2));

%% Perform optimization
for i=1:epochs,

  % Determine quadratic distances between neurons and data
   distNeuronsToData = determine_distance_relational_neurons_to_data( Distances, cNeurons );

  % Determine winner w_I*(x_j) for every datapoint (w* := [w_I*(x_1), ..., w_I*(x_m)])
   [err, w_star] = min( distNeuronsToData, [], 1 );

  % Determine k~_ij (distances to winner neurons)
   KK = GridDistances(:,w_star);

  % Neighbourhood Function
   HH = exp(-(KK)/lambda(i));

  % Update rule
   cNeurons = HH ./ (sum(HH,2) * ones(1,m));

  % Determine quantization error
   QuantizationErrors(i) = determine_qerror( sqrt( distNeuronsToData ) );

end

% Determine distance between neurons and data
distNeuronsToData = sqrt( determine_distance_relational_neurons_to_data( Distances, cNeurons ) );
sMap.cCodebook = cNeurons;

sTrain = som_set(sTrain,'time',datestr(now,0));
tl = length(sMap.trainhist);
sMap.trainhist(tl+1) = sTrain;

function distNeuronsToData = determine_distance_relational_neurons_to_data( Distances, cNeurons )

%
% Determines dissimilarities between relational neurons and data on which they are defined
%
% Distances ... (m x m)-Matrix of dissimilarities between data points
% cNeurons ... Coefficient (n x m)-Matrix of relational neurons
%
% Use squared distances to compare with other methods on euclidean distances !
%
% A.Hasenfuss (c) 2007-2008 (Optimized Version)
%

[n,m] = size( cNeurons );

% Determine distances between neurons and data
  temp1 = cNeurons * Distances;

  temp2 = zeros(n,1);
  for i = 1:n
    temp2(i,1) = 0.5 * (temp1(i,:) * cNeurons(i,:)');
  end

  distNeuronsToData = temp1 - temp2*ones(1,m);

function QuantizationError = determine_qerror( distNeuronsToData )

%
% Determines the quantization error for median neurons
%
% A.Hasenfuss (c) 2006
%

[WinnerDist, WinnerIdx] = min( distNeuronsToData, [], 1  );
QuantizationError = sum( WinnerDist.^2 )/2;

author: ilarinieminen

repo: SOM-Toolbox


text:
% Test file for chebtech/alias.m

function pass = test_alias(pref)

if ( nargin < 1 )
    pref = chebtech.techPref();
end

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    F = @sin;
    f = testclass.make(@(x) F(x), [], pref);

    k = 11;
    g = f;
    g.coeffs = testclass.alias(f.coeffs, k);
    x = testclass.chebpts(k);
    values = testclass.coeffs2vals(g.coeffs);
    pass(n, 1) = norm(values - F(x), inf) < 10*vscale(g)*eps;

    g.coeffs = testclass.alias(f.coeffs, 2);
    exact_values = sin(testclass.chebpts(2));
    values = testclass.coeffs2vals(g.coeffs);
    pass(n, 2) = size(g,1) == 2 && norm(values - exact_values, inf) < ...
        1e1*vscale(g)*eps;

    F = @sin;
    f = testclass.make(@(x) [F(x), -F(x)], [], pref);
    k = 11;
    g = f;
    g.coeffs = testclass.alias(f.coeffs, k);
    x = testclass.chebpts(k);
    values = testclass.coeffs2vals(g.coeffs);
    pass(n, 3) = norm(values - [F(x), -F(x)], inf) < 10*max(vscale(g)*eps);

    g = f;
    g.coeffs = testclass.alias(f.coeffs, 2);
    y = testclass.chebpts(2);
    exact_values = [sin(y) -sin(y)];
    values = testclass.coeffs2vals(g.coeffs);
    pass(n, 4) = size(g,1) == 2 && ...
        norm(values - exact_values, inf) < 10*max(vscale(g)*eps);

    F = @(x) sin(1000*x);
    f = testclass.make(@(x) [F(x), -F(x)], [], pref);
    k = 32;
    g = f;
    g.coeffs = testclass.alias(f.coeffs, k);
    x = testclass.chebpts(k);
    values = testclass.coeffs2vals(g.coeffs);
    pass(n, 5) = norm(values - [F(x), -F(x)], inf) < 1e3*max(vscale(g)*eps);

    k = 100;
    g = f;
    g.coeffs = testclass.alias(f.coeffs, k);
    x = testclass.chebpts(k);
    values = testclass.coeffs2vals(g.coeffs);
    pass(n, 6) = norm(values - [F(x), -F(x)], inf) < 1e4*max(vscale(g)*eps);

    g = f;
    g.coeffs = testclass.alias(f.coeffs, 2);
    y = testclass.chebpts(2);
    exact_values = [sin(1000*y) -sin(1000*y)];
    values = testclass.coeffs2vals(g.coeffs);
    pass(n, 7) = size(g, 1) == 2 && ...
        norm(values - exact_values, inf) < 1e3*eps;

    F = @(x) cos(1000*x);
    f = testclass.make(@(x) [F(x), -F(x)], [], pref);

    g = f;
    g.coeffs = testclass.alias(f.coeffs, 1);
    values = testclass.coeffs2vals(g.coeffs);
    pass(n, 8) = size(g, 1) == 1 && norm(values - [1, -1], inf) < ...
        1e2*max(vscale(g)*eps);

    g = f;
    g.coeffs = testclass.alias(f.coeffs, 2);
    y = testclass.chebpts(2);
    exact_values = [cos(1000*y) -cos(1000*y)];
    values = testclass.coeffs2vals(g.coeffs);
    pass(n, 9) = length(g) == 2 && ...
        norm(values - exact_values, inf) < 1e3*max(vscale(g)*eps);

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [weight]=best_expert(data, weight_o)
% expert function for the market function
% Input
%     data    -- market historical prices
%     weight_o -- portfolio weight for last trading day
%
% Output
%     weight  -- weight column vector
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 National Renewable Energy Laboratory and National
% Technology & Engineering Solutions of Sandia, LLC (NTESS).
% Under the terms of Contract DE-NA0003525 with NTESS,
% the U.S. Government retains certain rights in this software.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function blkStruct = slblocks

Browser.Library='PTOSim_Lib';
Browser.Name='PTO-Sim';
Browser.Choice=1; % Option to generate the required repository info in memory.

blkStruct.Browser=Browser;

author: WEC-Sim

repo: WEC-Sim

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


text:
function filename = utl_whichfile(identifier)
% Find the .m file which defines the given function identifier.
% Filename = utl_whichfile(Identifier)
%
% This is a special-purpose alternative to the builtin 'which', for use in deployed situations.
%
% In these cases, 'which' returns the filename of only an encrypted version of the defining
% function, so it cannot be used to access the m code. This function assumes that the source code
% for each function ships with the release, and returns the file name for a given identifier.
%
% In:
%   Identifier : identifier of some function
%
% Out:
%   Filename : .m file that contains the code which implements the
%              identified function
%
% Notes:
%   This function does not have all the facilities of which(), and assumes that there is one unique
%   source file for every looked up identifier.
%
%   If the global cell string tracking.paths.toolboxes exists, its contents will be used as root
%   directories; otherwise, env_translatepath('functions:/') will be used.
%
% See also:
%   which
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-04

persistent filenames; % table of m file names, indexed by identifier

try
    % look up the record for this file name
    filename = filenames.(identifier);
catch %#ok<CTCH>
    if isempty(filenames)
        % file table is empty (first-time run): build it recursively from root directories
        global tracking; %#ok<TLEV>
        if isfield(tracking,'paths') && isfield(tracking.paths,'toolboxes')
            root_dirs = tracking.paths.toolboxes;
        else
            root_dirs = {env_translatepath('functions:/')};
        end
        for d=1:length(root_dirs)
            populate_records(root_dirs{d}); end
    else
        error('No file was found that corresponds to the identifier %s. Make sure it is in the searched directories.',hlp_tostring(identifier));
    end
    try
        filename = filenames.(identifier);
    catch %#ok<CTCH>
        error('No file was found that corresponds to the identifier %s. Make sure it is in the searched directories.',hlp_tostring(identifier));
    end
end

    % recursively populate file table from an initial directory.
    function populate_records(path)
        if path(end) == filesep
            path(end) = []; end
        for info = dir(path)'
            if info.isdir
                % found a directory: recurse
                if ~any(strcmp(info.name,{'.','..'}))
                    populate_records([path filesep info.name]); end
            else
                % found a file? add it to the table if an .m file
                if length(info.name) > 2 && strcmp(info.name(end-1:end),'.m')
                    try
                        filenames.(info.name(1:end-2)) = [path filesep info.name];
                    catch
                        disp_once('utl_whichfile: found .m file that is not a function name: %s (ignoring)',[path filesep info.name]);
                    end
                end
            end
        end
    end

end

repo: REST


text:
function plotRelaxedFBA(sol, model,tol)

if sol.stat ~= 1
    disp('relaxedFBA did not complete successfully, nothing to display');
    disp('relaxedFBA problem infeasible, check relaxOption fields');
    return
end

if ~exist('tol','var')
    %Set the tolerance to distinguish between zero and non-zero flux
    feasTol = getCobraSolverParams('LP', 'feasTol');
    tol=feasTol*100;
end

%useful for numerical debugging
if 0
    tol=0;
end

optTol = getCobraSolverParams('LP', 'optTol');
sol.p(abs(sol.p)<optTol*10)=0;
sol.q(abs(sol.q)<optTol*10)=0;

if norm(model.S*sol.v-model.b)>optTol*10
    warning(['relaxedFBA relaxed steady state constraints too much, norm(S*v-b)= ' num2str(norm(model.S*sol.v-model.b))])
end

[v,r,p,q] = deal(sol.v,sol.r,sol.p,sol.q);

p(p<tol) = 0;%lower bound relaxation
q(q<tol) = 0;%upper bound relaxation
r(abs(r)<tol) = 0;%steady state constraint relaxation

%Summarise the proposed relaxation sol
printFlag=0;
lineChangeFlag=0;

[nMet,nRxn]=size(model.S);

%totals
fprintf('%u%s\n',nnz(r),' steady state constraints relaxed');
fprintf('%u%s\n',nnz(p & q & model.SIntRxnBool),' internal lower and upper bounds relaxed');
fprintf('%u%s\n',nnz(p & q & ~model.SIntRxnBool),' external lower and upper bounds relaxed');
fprintf('%u%s\n',nnz(p | q & ~model.SIntRxnBool),' external lower or upper bounds relaxed');
maxUB = max(max(model.ub),-min(model.lb));
minLB = min(-max(model.ub),min(model.lb));
intRxnFiniteBound = ((model.ub < maxUB) & (model.lb > minLB));
fprintf('%u%s\n',nnz(p & intRxnFiniteBound),' finite lower bounds relaxed');
exRxn00 = ((model.ub == 0) & (model.lb == 0));
fprintf('%u%s\n',nnz(p & exRxn00),' lower bounds relaxed on fixed reactions (lb=ub=0)');
fprintf('\n');

bool=r~=0;
if any(bool)
    fprintf('%u%s\n',nnz(r),' steady state constraints relaxed');
    fprintf('%-20s%12s%12s\n','mets{n}','r','dxdt')
    for m=1:nMet
        if bool(m)
            fprintf('%-20s%12.4g%12.4g\n',...
                model.mets{m},sol.r(m),model.b(m));
        end
    end
else
    fprintf('%s\n','No steady state constraints relaxed');
end
fprintf('\n');

bool=p~=0 & q==0 & model.SIntRxnBool;
if any(bool)
    fprintf('%u%s\n',nnz(bool),' internal reaction lower bounds relaxed...');
    fprintf('%-20s%12.4s%12.4s%12.4s%12.4s%12.4s%12.8s\n','rxns','p','lb','v','ub','q','formula');
    for n=1:nRxn
        if bool(n)
            rxnAbbrList=model.rxns(n);
            formulas = printRxnFormula(model, rxnAbbrList, printFlag, lineChangeFlag);
            fprintf('%-20s%12.4g%12.4g%12.4g%12.4g%12.4g%6g%-100s\n',...
            model.rxns{n},sol.p(n),model.lb(n),sol.v(n),model.ub(n),sol.q(n),formulas{1});
        end
    end
else
    fprintf('%s\n','No internal reaction lower bounds relaxed');
end
fprintf('\n');

bool=p==0 & q~=0 & model.SIntRxnBool;
if any(bool)
    fprintf('%u%s\n',nnz(bool),' internal reaction lower bounds relaxed...');
    fprintf('%-20s%12.4s%12.4s%12.4s%12.4s%12.4s%12.8s\n','rxns','p','lb','v','ub','q','formula');
    for n=1:nRxn
        if bool(n)
            rxnAbbrList=model.rxns(n);
            formulas = printRxnFormula(model, rxnAbbrList, printFlag, lineChangeFlag);
            fprintf('%-20s%12.4g%12.4g%12.4g%12.4g%12.4g%6g%-100s\n',...
            model.rxns{n},sol.p(n),model.lb(n),sol.v(n),model.ub(n),sol.q(n),formulas{1});
        end
    end
else
    fprintf('%s\n','No internal reaction lower bounds relaxed');
end
fprintf('\n');

bool=p~=0 & q==0 & ~model.SIntRxnBool;
if any(bool)
    fprintf('%u%s\n',nnz(bool),' external reaction lower bounds relaxed...');
    fprintf('%-20s%12.4s%12.4s%12.4s%12.4s%12.4s%12.8s\n','rxns','p','lb','v','ub','q','formula');
     for n=1:nRxn
        if bool(n)
            rxnAbbrList=model.rxns(n);
            formulas = printRxnFormula(model, rxnAbbrList, printFlag, lineChangeFlag);
            fprintf('%-20s%12.4g%12.4g%12.4g%12.4g%12.4g%6g%-100s\n',...
            model.rxns{n},sol.p(n),model.lb(n),sol.v(n),model.ub(n),sol.q(n),formulas{1});
        end
    end
else
    fprintf('%s\n','No external reaction lower bounds relaxed');
end
fprintf('\n');

bool=p==0 & q~=0 & ~model.SIntRxnBool;
if any(bool)
    fprintf('%u%s\n',nnz(bool),' external reaction lower bounds relaxed...');
    fprintf('%-20s%12.4s%12.4s%12.4s%12.4s%12.4s%12.8s\n','rxns','p','lb','v','ub','q','formula');
     for n=1:nRxn
        if bool(n)
            rxnAbbrList=model.rxns(n);
            formulas = printRxnFormula(model, rxnAbbrList, printFlag, lineChangeFlag);
            fprintf('%-20s%12.4g%12.4g%12.4g%12.4g%12.4g%6g%-100s\n',...
            model.rxns{n},sol.p(n),model.lb(n),sol.v(n),model.ub(n),sol.q(n),formulas{1});
        end
    end
else
    fprintf('%s\n','No external reaction lower bounds relaxed');
end
fprintf('\n');

author: opencobra

repo: cobratoolbox


text:
%  Figure 3.28      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%script to generate Fig. 3.28
%  Example 3.25

clf;
u=-1;  % magnitude of impulsive elevator input (Deg)
num=30*[1 -6];
den=[1 4 13 0];
t=0:.02:5;
y=impulse(u*num,den,t);
zero=[0 0];
tzero=[0 5];
axis([0 5 -2 16])
plot(t,y,'-',tzero,zero,'-'),grid
title('Fig. 3.28  Impulse response of aircraft altitude')
xlabel('Time (sec)')
ylabel('Altitude (ft)')

author: Sable

repo: mcbench-benchmarks


text:
function [ cobb ] = genOBBfrom2dOffset( pobb, rp )
%GENOBBFROM2DOFFSET generate box with absolute position from offset format
%
% INPUT: pobb - the parent box with absolute position.
%               rp - relative position (offsets) of child box under its
%               parent's local frame, including:
%                    (1) relative orientation - 1 bit
%                    (2) absolute sizes (under its local frame) - 2 bits
%                    (3) offsets - 2 bits
%                    (4) 4 classes (the relative layout, which two edges are close) - 4 bits
%                    (5) attachments (4 possible cases) - 4 bits
%                    (6) alignment (whether they are aligned or not) - 1 bit
%
% OUTPUT: cobb - the child box with absolute position

% relpos = [nfront(1);csize(1);csize(3);offsets(1);offsets(2);classvec;alignvec];

csize = [rp(2);0;rp(3)];
offsets = [rp(4);rp(5)];
classvec = rp(6:9);
[~,classindex] = max(classvec);
attachvec = rp(10:13);
[~,attachindex] = max(attachvec);
alignment = rp(14);

if(abs(alignment)>0.5)
    if(abs(rp(1))<0.5)
        rotation = 0;
    else
        rotation = 1;
    end
else
    rotation = rp(1);
end
nfront =  [rotation;0;sqrt(1-rotation*rotation)];

% use the alignment to adjust the offsets
if(mod(attachindex-1,2)==1)
    offset(1) = 0;
end
if(floor((attachindex-1)/2)==1)
    offset(2) = 0;
end

% parent box info
pcent = pobb(1:3);
pfront = pobb(4:6);
pup = pobb(7:9);
psize = pobb(10:12);
paxes = cross(pfront,pup);
paxes = paxes/norm(paxes,2);
transform = genTransMat(pfront,pup);
transform = inv(transform);

% child box info
cfront = transform*nfront;
cfront = abs(cfront);% assume all elements in orientation vec are positive
cup = [0;1;0];
caxes = cross(cfront,cup);
caxes = caxes/norm(caxes,2);

% parent edges info
pminx = -psize(1)/2;
pmaxx = psize(1)/2;
pminy = -psize(3)/2;
pmaxy = psize(3)/2;

if(mod(classindex-1,2)==0)
    cmaxx = pmaxx - offsets(1);
    cminx = cmaxx - csize(1);
else
    cminx = pminx + offsets(1);
    cmaxx = cminx + csize(1);
end
if(floor((classindex-1)/2)==0)
    cmaxy = pmaxy - offsets(2);
    cminy = cmaxy - csize(3);
else
    cminy = pminy + offsets(2);
    cmaxy = cminy + csize(3);
end

ncent = [(cmaxx+cminx)/2;0;(cmaxy+cminy)/2];
ccent = transform*ncent + pcent;
cobb = [ccent;cfront;cup;csize];

end

repo: GRAINS


text:
function [rVect, vVect] = getStateGMm(bodyInfo, time, celBodyData)

% Here we are looking for Kepler elements of the body for a given time with
% G(M+m) factor
%     sma = bodyInfo.sma;
%     ecc = bodyInfo.ecc;
%     inc = AngleZero2Pi(deg2rad(bodyInfo.inc));
%     raan = AngleZero2Pi(deg2rad(bodyInfo.raan));
%     argp = AngleZero2Pi(deg2rad(bodyInfo.arg));
%     M0 = deg2rad(bodyInfo.mean);
%
%     n = computeMeanMotion(sma,bodyInfo.gm+parentbodyInfo.gm);
%     deltaT = time - bodyInfo.epoch;
%     M = M0 + n.*deltaT;
%     tru = computeTrueAnomFromMean(M, ecc);
    parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);

    [rVectB, vVectB] = getStateAtTime(bodyInfo, time, getParentGM(bodyInfo, celBodyData));
    [sma, ecc, inc, raan, arg, tru] = getKeplerFromState(rVectB, vVectB, getParentGM(bodyInfo, celBodyData));

% Here we are calculating position and velocity for these parameters, but
% with GM factor
    [rVect,vVect]=getStatefromKepler(sma, ecc, inc, raan, arg, tru, parentBodyInfo.gm, true);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout = mygui(varargin)
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @mygui_OpeningFcn, ...
                   'gui_OutputFcn',  @mygui_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end

function mygui_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
guidata(hObject, handles);
% changeIcon
jFrame=get(handles.figure1,'javaframe');
jIcon=javax.swing.ImageIcon('icon.gif');
jFrame.setFigureIcon(jIcon);
% changeUserData
set(handles.pushbutton_run,'UserData',0);

function varargout = mygui_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

function pushbutton_run_Callback(hObject, eventdata, handles)
if get(handles.pushbutton_run,'UserData')==1, return; end
set(handles.pushbutton_run,'UserData',1);
for n=1:40
    if get(handles.pushbutton_run,'UserData')==0, break; end
    cla; text(0.5,0.5,num2str(n)); pause(0.25);
end

function pushbutton_stop_Callback(hObject, eventdata, handles)
set(handles.pushbutton_run,'UserData',0);


text:
% Test file for bndfun/innerProduct.m

function pass = test_innerProduct(pref)

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

% Set a tolerance.  (pref.chebfuneps doesn't matter here.)
tol = 10*eps;

% Set a domain
dom = [-2 7];

% Fixed arbitrary numbers to use as multiplicative constants.
alpha = -0.194758928283640 + 0.075474485412665i;
beta = -0.526634844879922 - 0.685484380523668i;

%%
% Spot-check a few known results.

f = bndfun(@(x) sin(2*pi*x), struct('domain', dom), pref);
g = bndfun(@(x) cos(2*pi*x), struct('domain', dom), pref);
pass(1) = abs(innerProduct(f, g)) < ...
    10*eps;

g = bndfun(@(x) cos(4*pi*x), struct('domain', dom), pref);
pass(2) = abs(innerProduct(f, g)) < ...
    10*eps;

f = bndfun(@(x) exp(x), struct('domain', dom), pref);
g = bndfun(@(x) exp(-x), struct('domain', dom), pref);
pass(3) = abs(innerProduct(f, g) - 9) < 1e2*max(get(f, 'vscale'), ...
    get(g, 'vscale'))*eps;

g = bndfun(@(x) sin(x), struct('domain', dom), pref);
exact = exp(7)*(sin(7) - cos(7))/2 - exp(-2)*(sin(-2) - cos(-2))/2;
pass(4) = abs(innerProduct(f, g) - exact) < max(get(f, 'vscale'), ...
    get(g, 'vscale'))*10*eps;

%%
% Check a few known properties.

f = bndfun(@(x) exp(1i*x) - 1, struct('domain', dom), pref);
g = bndfun(@(x) 1./(1 + 1i*x.^2), struct('domain', dom), pref);
h = bndfun(@(x) sinh(x*exp(pi*1i/6)), struct('domain', dom), pref);

ip1 = innerProduct(alpha*f, beta*g);
ip2 = conj(alpha)*beta*innerProduct(f, g);
pass(5) = abs(ip1 - ip2) < 10*tol;

ip1 = innerProduct(g, h);
ip2 = innerProduct(h, g);
pass(6) = abs(ip1 - conj(ip2)) < tol;

ip1 = innerProduct(f + g, h);
ip2 = innerProduct(f, h) + innerProduct(g, h);
pass(7) = abs(ip1 - ip2) < ...
    max((get(f, 'vscale') + get(g, 'vscale')), get(h, 'vscale'))*tol;

ip1 = innerProduct(f, g + h);
ip2 = innerProduct(f, g) + innerProduct(f, h);
pass(8) = abs(ip1 - ip2) < ...
    max((get(g, 'vscale') + get(h, 'vscale')), get(f, 'vscale'))*tol;

nf2 = innerProduct(f, f);
ng2 = innerProduct(g, g);
nh2 = innerProduct(h, h);
n2vals = [nf2 ; ng2 ; nh2];
pass(9) = isreal(n2vals) && all(n2vals >= 0);

%%
% Check operation for array-valued bndfun objects.

f = bndfun(@(x) [sin(x) cos(x)], struct('domain', dom), pref);
g = bndfun(@(x) [exp(x) 1./(1 + x.^2) airy(x)], struct('domain', dom), pref);
ip = innerProduct(f, g);
exact = [-53.1070904269318222 0.0025548835039100  -0.4683303433821355;
         773.70343924989359096771 1.3148120368924471 0.6450791915572742];
pass(10) = norm(ip(:) - exact(:), inf) < 10*max([eps, ...
    eps])*max([get(f, 'vscale') get(g, 'vscale')]);

%%
% Check error conditition
% Can't take the inner product of a bndfun and a non-bndfun.
try
    ip = innerProduct(f, 2); %#ok<NASGU>
    pass(11) = false;
catch ME
    pass(11) = strcmp(ME.identifier, ...
        'CHEBFUN:BNDFUN:innerProduct:input');
end

%% Test on singular function:

pow1 = -0.3;
pow2 = -0.5;
op1 = @(x) (x - dom(2)).^pow1.*sin(x);
op2 = @(x) (x - dom(2)).^pow2.*cos(3*x);
pref.blowup = true;
data.domain = dom;
data.exponents = [0 pow1];
f = bndfun(op1, data, pref);
data.domain = dom;
data.exponents = [0 pow2];
g = bndfun(op2, data, pref);
I = innerProduct(f,g);
I_exact = -0.65182492763883119+0.47357853074362785i;
err = abs(I - I_exact);
tol = 5e2*eps*abs(I_exact);
pass(12) = err < tol;

end


text:
function out = brighterremap(x)
% BRIGHTERREMAP Brighter set of parameters for density remap
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

    out = uint8( amplitudetodensity(x,60,40) );
end
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function a = r82row_permute ( n, a, p )

%*****************************************************************************80
%
%% R82ROW_PERMUTE permutes an R82ROW in place.
%
%  Discussion:
%
%    An R82ROW is a (2,N) array of R8's.
%
%    This routine permutes an array of real "objects", but the same
%    logic can be used to permute an array of objects of any arithmetic
%    type, or an array of objects of any complexity.  The only temporary
%    storage required is enough to store a single object.  The number
%    of data movements made is N + the number of cycles of order 2 or more,
%    which is never more than N + N/2.
%
%  Example:
%
%    Input:
%
%      N = 5
%      P = (   2,    4,    5,    1,    3 )
%      A = ( 1.0,  2.0,  3.0,  4.0,  5.0 )
%          (11.0, 22.0, 33.0, 44.0, 55.0 )
%
%    Output:
%
%      A    = (  2.0,  4.0,  5.0,  1.0,  3.0 )
%             ( 22.0, 44.0, 55.0, 11.0, 33.0 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of objects.
%
%    Input, real A(2,N), the array to be permuted.
%
%    Input, integer P(N), the permutation.  P(I) = J means
%    that the I-th element of the output array should be the J-th
%    element of the input array.  P must be a legal permutation
%    of the integers from 1 to N, otherwise the algorithm will
%    fail catastrophically.
%
%    Output, real A(2,N), the permuted array.
%

%
%  Search for the next element of the permutation that has not been used.
%
  for istart = 1 : n

    if ( p(istart) < 0 )

      continue

    elseif ( p(istart) == istart )

      p(istart) = - p(istart);
      continue

    else

      a_temp(1:2) = a(1:2,istart);
      iget = istart;
%
%  Copy the new value into the vacated entry.
%
      while ( 1 )

        iput = iget;
        iget = p(iget);

        p(iput) = - p(iput);

        if ( iget < 1 || n < iget )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'R82ROW_PERMUTE - Fatal error!\n' );
          fprintf ( 1, '  A permutation index is out of range.\n' );
          fprintf ( 1, '  P(%d) = %d\n', iput, iget );
          error ( 'R82ROW_PERMUTE - Fatal error!' );
        end

        if ( iget == istart )
          a(1:2,iput) = a_temp(1:2)';
          break
        end

        a(1:2,iput) = a(1:2,iget);

      end

    end

  end
%
%  Restore the signs of the entries.
%
% p(1:n) = -p(1:n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h = plotTable(table, varargin);
% h = plotTable(table, ['property', val]): display a text
% table in a set of axes.
%
% Can pass in pairs of properties and values for the text.
% Returns an array of handles to the text.
%
% Note this is very primitive.
%
% 08/04 ras.
h = [];

nrows = size(table,1);
ncols = size(table,2);

% ensure table is a cell-of-strings
for row = 1:nrows
    for col = 1:ncols
        if ~ischar(table{row,col})
            table{row,col} = num2str(table{row,col});
        end
    end
end

% count the text width of each cell, get column
% width
for col = 1:ncols
    for row = 1:nrows
        sz(row,col) = length(table(row,col));
    end

    width(1,col) = max(sz(:,col))+4;
end
xpos = cumsum(width);

% set axes
cla;
AX = [1 sum(width)+3 1 nrows+3];
axis(AX); axis off
set(gca,'XTick',[],'YTick',[],'Box','on');

% plot the text
for row = 1:nrows
    for col = 1:ncols
        X = xpos(col)-2;
        Y = nrows-row+2;
        if row==1
            h(end+1) = text(X,Y,table{row,col},'FontSize',10);
        elseif col==1
            h(end+1) = text(X,Y,table{row,col},'FontSize',8,'FontWeight','bold');
        else
            h(end+1) = text(X,Y,table{row,col},'FontSize',8);
        end

    end
end

% apply optional text properties
if length(varargin)>1
    for j = 2:2:length(varargin)
        set(h, varargin{j-1}, varargin{j});
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% std_rmalldatafields - remove all data fields from STUDY (before saving
%                       it for instance.
%
% Usage:
%   STUDY = std_rmalldatafields(STUDY, type);
%
% Input:
%   STUDY   - EEGLAB study structure
%
% Optional input:
%   type    - ['chan'|'clust'|'both'] remove from changrp channel location
%             structure, cluster structure or both. Default is 'both'.
%
% Ouput:
%   STUDY   - updated EEGLAB study structure
%
% Author: Arnaud Delorme, CERCO/CNRS, UCSD, 2010-

% Copyright (C) Arnaud Delorme, arno@salk.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

function STUDY = std_rmalldatafields(STUDY, chanorcomp);

if nargin < 1
    help std_rmalldatafields;
    return;
end;
if nargin < 2
    chanorcomp = 'both';
end;

fields = { 'erpdata'  'erptimes'  'erpdatatrials'  'erptrialinfo' ...
               'specdata' 'specfreqs' 'specdatatrials' 'spectrialinfo' ...
               'erspdata' 'erspbase'  'erspfreqs'      'ersptimes' 'erspdatatrials' 'ersptrialinfo'  ...
               'topo' 'topox' 'topoy' 'topoall' 'topopol' ...
               'itcdata'  'itcfreqs'  'itctimes'       'itcdatatrials' 'itctrialinfo'  ...
               'erpimdata' 'erpimtrials' 'erpimevents' 'erpimtimes' 'erspsubjinds' 'itcsubjinds' 'dipoles' ...
               'data' 'datatimes' 'datasortvals' 'datacontinds' 'centroid' };
for ff = 1:length(fields)
    if strcmpi(chanorcomp, 'data') || strcmpi(chanorcomp, 'both')
        if isfield(STUDY.changrp, fields{ff})
            STUDY.changrp = rmfield(STUDY.changrp, fields{ff} );
        end;
    end;
    if strcmpi(chanorcomp, 'clust') || strcmpi(chanorcomp, 'both')
        if isfield(STUDY.cluster, fields{ff})
            STUDY.cluster = rmfield(STUDY.cluster, fields{ff} );
        end;
    end;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [parameters, stderrors, LLF, ht, resids, summary] = garchfind(data, models, distributions, ar, ma, p, q, options)
%{
[parameters, stderrors, LLF, ht, resids, summary] = garch(data, model, distr, ar, ma, x, p, q, y,startingvalues, options)
-----------------------------------------------------------------------
 PURPOSE:
 Finds the combination of models and distributions that better fits the
 data based on a set of criteria (i.e.  largest log likelihood value
 and the smallest AIC and BIC criteria).
----------------------------------------------------------------------
 USAGE:
 garchfind(data, models, distributions, ar, ma, p, q)

 INPUTS:
 data:          (T x 1) vector of data
 models:        a character vector with models to be estimated
 distributions: a character vector with distributions
 ar:            length of AR
 ma:            length of MA
 p:             length of ARCH
 q:             length of GARCH
 options:       set of options

 OUTPUTS based on the estimation of the best model:
 parameters:   a vector of parameters a0, a1, b0, b1, b2, ...
 stderrors:    standard errors estimated by the inverse Hessian (fmincon)
 LFF:          the value of the Log-Likelihood Function
 ht:           vector of conditional variances
 resids:       vector of residuals
 summary:      summary of results including:
                -model specification, distribution and statistics
                -optimization options
                -t-statistics
                -robust standard errors: (HESSIAN^-1)*cov(scores)*(HESSIAN^-1)
                -scores: numerical scores for M testing
-----------------------------------------------------------------------
Author:
 Alexandros Gabrielsen, a.gabrielsen@city.ac.uk
 Date: 08/2011
-----------------------------------------------------------------------
%}
if nargin == 0
    error('Data, GARCH Models, Distribution, AR, MA, ARCH, GARCH, Options')
end

if size(data,2) > 1
   error('Data vector should be a column vector')
end

if (length(ar) > 1) | (length(ma) > 1) | ar < 0 | ma < 0
    error('AR and MA should be positive scalars')
end

if (length(p) > 1) | (length(q) > 1) | p < 0 | q < 0
    error('P and Q should be positive scalars')
end

Criteria=[];
index=1;
for i = 1:size(models,1)
    for j = 1:size(distributions,1)
        for a=1:ar
            for b=1:ma
                for c=1:p
                    for d=1:q
                         fprintf('Progress: ARMA(%1.0f,%1.0f,%1.0f)-%s(%1.0f,%1.0f,%1.0f) - %s distribution\n', a, b, 0, strcat(models(i,:)), c, d, 0, strcat(distributions(j,:)))
                         eval(['[parameters, stderrors, LLF_temp] = garch(data,''' strcat(models(i,:)),''',''', strcat(distributions(j,:)), ''',a,b,0,c,d,0,[],options);']);
                         Criteria(index,:) = [LLF_temp, -2*LLF_temp+2*size(parameters,1), -2*LLF_temp+size(parameters,1)*log(size(data,1)), i, j, a, b, c,d];
                         index=index+1;
                    end
                end
            end
        end
    end
end

% This estimates a vector of results as long as the criteria are satisfied
holder=sortrows(Criteria,1);
z=size(holder,1)+1;
fprintf('\n The top models are: \n')
for i = 1:size(holder,1)
    if  find(holder == max(holder(1:z-i,1))) == find(holder == min(holder(1:z-i,2))) - size(holder(1:z-i,1),1) & find(holder == max(holder(1:z-i,1))) == find(holder == min(holder(1:z-i,3))) - 2*size(holder(1:z-i,1),1)
        fprintf('ARMA(%1.0f,%1.0f,%1.0f)-%s(%1.0f,%1.0f,%1.0f) - %s distribution\n', holder(z-i,6), holder(z-i,7), 0, strcat(models(holder(z-i,4),:)), holder(z-i,8), holder(z-i,9), 0, strcat(distributions(holder(z-i,5),:)))
    end
end

% Finds the best model so as to estimate the model
if find(Criteria == max(Criteria(:,1))) == find(Criteria == min(Criteria(:,2))) - size(Criteria,1) & find(Criteria == max(Criteria(:,1))) == find(Criteria == min(Criteria(:,3))) - 2*size(Criteria,1)
[parameters, stderrors, LLF, ht, resids, summary] = garch(data, strcat(models(Criteria(find(Criteria == max(Criteria(:,1))),4),:)), ...
    strcat(distributions(Criteria(find(Criteria == max(Criteria(:,1))),5),:)), ...
    Criteria(find(Criteria == max(Criteria(:,1))),6),...
    Criteria(find(Criteria == max(Criteria(:,1))),7),0, ...
    Criteria(find(Criteria == max(Criteria(:,1))),8), ...
    Criteria(find(Criteria == max(Criteria(:,1))),9),0,[],options);
end
end


text:
function r = scale(a,f)

%tstoolbox/@achse/scale
%   Syntax:
%     * r = scale(a,f)
%
%   Scale achses delta by factor f.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

if isa(f, 'double')
	r = a;
	r.delta = f * a.delta;
	r.values = f * a.values;
else
	error('need scalar factor to scale achse')
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
classdef CopeakLoss < dagnn.Loss
    properties
        eps = 10^-5;
        NumInputs = [];
        c = [];
        x = [];
    end

    methods
        function outputs = forward(obj, inputs, params)
            if length(inputs) == 2
                obj.NumInputs = [length(inputs{1}) length(inputs{2})];
                obj.c = cat(1, gpuArray.ones(obj.NumInputs(1), 1), -gpuArray.ones(obj.NumInputs(2), 1));
                obj.x = cat(1, inputs{1}, inputs{2});
                a = -obj.c.*obj.x ;
                b = max(0, a) ;
                t = b + log(exp(-b) + exp(a-b)) ;
                outputs{1} = mean(t(1:obj.NumInputs(1))) + mean(t(obj.NumInputs(1)+1:end));
                NumData = length(inputs{1}) + length(inputs{2});
            else
                obj.NumInputs = length(inputs{1});
                obj.c = gpuArray.ones(obj.NumInputs(1), 1);
                obj.x = inputs{1};
                a = -obj.c.*obj.x ;
                b = max(0, a) ;
                t = b + log(exp(-b) + exp(a-b)) ;
                outputs{1} = mean(t);
                NumData = length(inputs{1});
            end
            n = obj.numAveraged;
            m = n + NumData;
            obj.average = (n * obj.average + double(gather(outputs{1} * NumData))) / m ;
            obj.numAveraged = m;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            y = - derOutputs{1} .* obj.c ./ (1 + exp(obj.c.*obj.x)) ;
            derInputs{1} = y(1:obj.NumInputs(1)) / obj.NumInputs(1);
            if length(inputs) == 2
                derInputs{2} = y(obj.NumInputs(1)+1:end) / obj.NumInputs(2);
            end
            derParams = {};
        end

        function reset(obj)
            obj.average = 0 ;
            obj.numAveraged = 0 ;
        end
        function obj = CopeakLoss(varargin)
            obj.load(varargin) ;
        end
    end
end

author: KuangJuiHsu


text:
function dt = samplePeriod

% SAMPLEPERIOD  Sample period.
%   SAMPLEPERIOD returns the time elapsed sinc last time called.
%
%   See also TIC, TOC.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

dt = toc;
tic;

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [nScales, scale_step, scaleFactors, scale_filter, params] = init_scale_filter(params)

% Initialize the scale filter parameters. Uses the fDSST scale filter.

init_target_sz = params.init_sz(:)';

nScales = params.number_of_scales_filter;
scale_step = params.scale_step_filter;

scale_sigma = params.number_of_interp_scales * params.scale_sigma_factor;

scale_exp = (-floor((nScales-1)/2):ceil((nScales-1)/2)) * params.number_of_interp_scales/nScales;
scale_exp_shift = circshift(scale_exp, [0 -floor((nScales-1)/2)]);

interp_scale_exp = -floor((params.number_of_interp_scales-1)/2):ceil((params.number_of_interp_scales-1)/2);
interp_scale_exp_shift = circshift(interp_scale_exp, [0 -floor((params.number_of_interp_scales-1)/2)]);

scale_filter.scaleSizeFactors = scale_step .^ scale_exp;
scale_filter.interpScaleFactors = scale_step .^ interp_scale_exp_shift;

ys = exp(-0.5 * (scale_exp_shift.^2) /scale_sigma^2);
scale_filter.yf = single(fft(ys));
scale_filter.window = single(hann(size(ys,2)))';

%make sure the scale model is not to large, to save computation time
if params.scale_model_factor^2 * prod(init_target_sz) > params.scale_model_max_area
    params.scale_model_factor = sqrt(params.scale_model_max_area/prod(init_target_sz));
end

%set the scale model size
params.scale_model_sz = max(floor(init_target_sz * params.scale_model_factor), [8 8]);

scale_filter.max_scale_dim = strcmp(params.s_num_compressed_dim,'MAX');
if scale_filter.max_scale_dim
    params.s_num_compressed_dim = length(scale_filter.scaleSizeFactors);
end

% Scale factor for the translation filter
scaleFactors = 1;

author: martin-danelljan

repo: ECO


text:
% -------------------------------------------------------------------------------------------------------------------------
function y = max_score_err(x, y_gt)
% x is [m1, m2, 1, b]
% numel(y_gt) is b
% The dimensions m1 and m2 are odd numbers.
%
%   Luca Bertinetto, Jack Valmadre, Joao F. Henriques, 2016
% -------------------------------------------------------------------------------------------------------------------------

[m1, m2, k, b] = size(x);
assert(mod(m1, 2) == 1);
assert(mod(m2, 2) == 1);
assert(k == 1);

assert(numel(y_gt) == b);

y_gt = reshape(y_gt, [1, 1, 1, b]);
pos = y_gt > 0;
neg = y_gt < 0;

x = gather(x);

% % Express error e as a linear function of scores x.
% h_center = center_mask(m1, m2, b);
% h_max = max_mask(x);
% h = bsxfun(@times, pos, h_center) - bsxfun(@times, neg, h_max);
% % e is positive if classified correctly, negative if incorrectly.
% e = sum(sum(h .* x, 1), 2);

x_center = center_score(x);
x_max = max_score(x);
e = zeros(b, 1);
e(pos) = x_center(pos);
e(neg) = -x_max(neg);

y = sum(e <= 0);

end

function v = center_score(x)
    [m1, m2, ~, b] = size(x);
    c1 = (m1+1) / 2;
    c2 = (m2+1) / 2;
    v = x(c1, c2, :, :);
end

function h = center_mask(m1, m2, b)
% This should satisfy x .* center_mask(m1, m2, b) == center_score(x).
    c1 = (m1+1) / 2;
    c2 = (m2+1) / 2;
    h = zeros(m1, m2, 1, b, 'single');
    h(c1, c2, :, :) = 1;
end

function v = max_score(x)
    [m1, m2, ~, b] = size(x);
    v = max(max(x, [], 1), [], 2);
end

function h = max_mask(x)
% This should satisfy x .* max_mask(x) == max_score(x).
    [m1, m2, ~, b] = size(x);
    x = reshape(x, [m1*m2, b]);
    [~, u] = max(x);
    assert(numel(u) == b);
    h = zeros(m1*m2, b, 'single');
    h(sub2ind(size(x), u, 1:8)) = 1;
    h = reshape(h, [m1, m2, 1, b]);
end

author: bertinetto

repo: siamese-fc


text:
% transform radians to degrees

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

deg=rad*180/pi;

repo: ARTE


text:
% pop_importegimat() - import EGI Matlab segmented file
%
% Usage:
%   >> EEG = pop_importegimat(filename, srate, latpoint0);
%
% Inputs:
%   filename    - Matlab file name
%   srate       - sampling rate
%   latpoint0   - latency in sample ms of stimulus presentation.
%                 When data files are exported using Netstation, the user specify
%                 a time range (-100 ms to 500 ms for instance). In this
%                 case, the latency of the stimulus is 100 (ms). Default is 0 (ms)
%
% Output:
%   EEG        - EEGLAB dataset structure
%
% Authors: Arnaud Delorme, CERCO/UCSD, Jan 2010

% Copyright (C) Arnaud Delorme, arno@salk.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

function [EEG com] = pop_importegimat(filename, srate, latpoint0, dataField);

    EEG = [];
    com = '';
    if nargin < 3, latpoint0 = 0; end;
    if nargin < 4, dataField = 'Session'; end;
    if nargin < 1
        % ask user
        [filename, filepath] = uigetfile('*.mat', 'Choose a Matlab file from Netstation -- pop_importegimat()');
        if filename == 0 return; end;
        filename = fullfile(filepath, filename);
        tmpdata = load('-mat', filename);
        fieldValues = fieldnames(tmpdata);
        sessionPos = strmatch('Session', fieldValues);
        posFieldData = 1;
        if ~isempty(sessionPos), posFieldData = sessionPos; end;

        if ~isfield(tmpdata, 'samplingRate'), srate = 250; else srate = tmpdata.samplingRate; end;
        % epoch data files only
        promptstr    = { { 'style' 'text'       'string' 'Sampling rate (Hz)' } ...
                         { 'style' 'edit'       'string'  int2str(srate) } ...
                         { 'style' 'text'       'string' 'Sample latency for stimulus (ms)' } ...
                         { 'style' 'edit'       'string' '0' } ...
                         { 'style' 'text'       'string' 'Field containing data' } ...
                         { 'style' 'popupmenu'  'string' fieldValues 'value' posFieldData } ...
                         };
        geometry = { [1 1] [1 1] [1 1] };
        result       = inputgui( 'geometry', geometry, 'uilist', promptstr, ...
                                 'helpcom', 'pophelp(''pop_importegimat'')', ...
                                 'title', 'Import a Matlab file from Netstation -- pop_importegimat()');

        if length(result) == 0 return; end;
        srate = str2num(result{1});
        latpoint0 = str2num(result{2});
        dataField = fieldValues{result{3}};
        if isempty(latpoint0), latpoint0 = 0; end;
    end;

    EEG = eeg_emptyset;
    fprintf('Reading EGI Matlab file %s\n', filename);
    tmpdata = load('-mat', filename);
    if isfield(tmpdata, 'samplingRate') % continuous file
        srate = tmpdata.samplingRate;
    end;

    fieldValues = fieldnames(tmpdata);
    if all(cellfun(@(x)isempty(findstr(x, 'Segment')), fieldValues))
        EEG.srate = srate;
        indData = strmatch(dataField, fieldValues);
        EEG.data = tmpdata.(fieldValues{indData(1)});
        EEG = eeg_checkset(EEG);
        EEG = readegilocs(EEG);

        com = sprintf('EEG = pop_importegimat(''%s'');', filename);
    else
        % get data types
        % --------------
        allfields = fieldnames(tmpdata);
        for index = 1:length(allfields)
            allfields{index} = allfields{index}(1:findstr(allfields{index}, 'Segment')-2);
        end;
        datatypes = unique_bc(allfields);
        datatypes(cellfun(@isempty, datatypes)) = [];

        % read all data
        % -------------
        counttrial = 1;
        EEG.srate  = srate;
        latency = (latpoint0/1000)*EEG.srate+1;
        for index = 1:length(datatypes)
            tindex = 1;
            for tindex = 1:length(allfields)
                if isfield(tmpdata, sprintf('%s_Segment%d', datatypes{index}, tindex))
                    datatrial = getfield(tmpdata, sprintf('%s_Segment%d', datatypes{index}, tindex));
                    if counttrial == 1
                        EEG.pnts = size(datatrial,2);
                        EEG.data = repmat(single(0), [size(datatrial,1), size(datatrial,2), 1000]);
                    end;
                    EEG.data(:,:,counttrial) = datatrial;

                    EEG.event(counttrial).type    = datatypes{index};
                    EEG.event(counttrial).latency = latency;
                    EEG.event(counttrial).epoch   = counttrial;

                    counttrial = counttrial+1;
                    latency = latency + EEG.pnts;
                end;
            end;
        end;
        fprintf('%d trials read\n', counttrial-1);
        EEG.data(:,:,counttrial:end) = [];

        EEG.setname  = filename(1:end-4);
        EEG.nbchan   = size(EEG.data,1);
        EEG.trials   = counttrial-1;
        if latpoint0 ~= 1
            EEG.xmin     = -latpoint0/1000;
        end;
        EEG = eeg_checkset(EEG);

        % channel location
        % ----------------
        if all(EEG.data(end,1:10) == 0)
            disp('Deleting empty data reference channel (reference channel location is retained)');
            EEG.data(end,:)   = [];
            EEG.nbchan        = size(EEG.data,1);
            EEG = eeg_checkset(EEG);
        end;
        EEG = readegilocs(EEG);

        com = sprintf('EEG = pop_importegimat(''%s'', %3.2f, %3.2f, %d);', filename, srate, latpoint0, dataField);
    end;


text:
function line_monte_carlo_test ( )

%*****************************************************************************80
%
%% LINE_MONTE_CARLO_TEST tests the LINE_MONTE_CARLO library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_MONTE_CARLO_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the LINE_MONTE_CARLO library.\n' );

  line_monte_carlo_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_MONTE_CARLO_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
2     % problem
4     % grid parameter
2     % discretisation
0.25  % Stokes stabilization parameter

%% Data file for test problem S2


text:
function saveExpression(obj, export_path, varargin)
    % save the symbolic expressions of the constraints matrices to MX
    % binary files
    %
    % Parameters:
    %  export_path: the path to export the file @type char
    %  varargin: variable input parameters @type varargin
    %   ForceExport: force the export @type logical

    cellfun(@(x)save_funcs(x,export_path, varargin{:}), obj.ActualFuncs, ...
        'UniformOutput',false);
    cellfun(@(x)save_funcs(x,export_path, varargin{:}), obj.DesiredFuncs, ...
        'UniformOutput',false);

    %     if ~isempty(obj.OutputFuncs)
    %         cellfun(@(x)save_funcs(x,export_path, varargin{:}), obj.OutputFuncs, ...
    %             'UniformOutput',false);
    %     end

    if ~isempty(obj.tau_)
        cellfun(@(x)save_funcs(x,export_path, varargin{:}), obj.PhaseFuncs, ...
            'UniformOutput',false);
    end

    function save_funcs(x, export_path, varargin)
        if ~isempty(x), save(x,export_path, varargin{:}); end
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
classdef GrEA < ALGORITHM
% <many> <real/integer/label/binary/permutation>
% Grid-based evolutionary algorithm
% div --- --- The number of divisions in each objective

%------------------------------- Reference --------------------------------
% S. Yang, M. Li, X. Liu, and J. Zheng, A grid-based evolutionary algorithm
% for many-objective optimization, IEEE Transactions on Evolutionary
% Computation, 2013, 17(5): 721-736.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            Div = [0 45 15 10 9 9 8 8 10 12];
            div = Algorithm.ParameterSet(Div(min(Problem.M,10)));

            %% Generate random population
            Population = Problem.Initialization();

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = MatingSelection(Population.objs,div);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                Population = EnvironmentalSelection([Population,Offspring],Problem.N,div);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ a, seed ] = r8mat_orth_uniform ( n, seed )

%*****************************************************************************80
%
%% R8MAT_ORTH_UNIFORM returns a random orthogonal matrix.
%
%  Properties:
%
%    The inverse of A is equal to A'.
%
%    A * A'  = A' * A = I.
%
%    Columns and rows of A have unit Euclidean norm.
%
%    Distinct pairs of columns of A are orthogonal.
%
%    Distinct pairs of rows of A are orthogonal.
%
%    The L2 vector norm of A*x = the L2 vector norm of x for any vector x.
%
%    The L2 matrix norm of A*B = the L2 matrix norm of B for any matrix B.
%
%    The determinant of A is +1 or -1.
%
%    All the eigenvalues of A have modulus 1.
%
%    All singular values of A are 1.
%
%    All entries of A are between -1 and 1.
%
%  Discussion:
%
%    Thanks to Eugene Petrov, B I Stepanov Institute of Physics,
%    National Academy of Sciences of Belarus, for convincingly
%    pointing out the severe deficiencies of an earlier version of
%    this routine.
%
%    Essentially, the computation involves saving the Q factor of the
%    QR factorization of a matrix whose entries are normally distributed.
%    However, it is only necessary to generate this matrix a column at
%    a time, since it can be shown that when it comes time to annihilate
%    the subdiagonal elements of column K, these (transformed) elements of
%    column K are still normally distributed random values.  Hence, there
%    is no need to generate them at the beginning of the process and
%    transform them K-1 times.
%
%    For computational efficiency, the individual Householder transformations
%    could be saved, as recommended in the reference, instead of being
%    accumulated into an explicit matrix format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Pete Stewart,
%    Efficient Generation of Random Orthogonal Matrices With an Application
%    to Condition Estimators,
%    SIAM Journal on Numerical Analysis,
%    Volume 17, Number 3, June 1980, pages 403-409.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A(N,N), the orthogonal matrix.
%
%    Output, integer SEED, a seed for the random number generator.
%
  a = zeros ( n, n );
%
%  Start with A = the identity matrix.
%
  for i = 1 : n
    for j = 1 : n
      if ( i == j )
        a(i,j) = 1.0;
      else
        a(i,j) = 0.0;
      end
    end
  end
%
%  Now behave as though we were computing the QR factorization of
%  some other random matrix.  Generate the N elements of the first column,
%  compute the Householder matrix H1 that annihilates the subdiagonal elements,
%  and set A := A * H1' = A * H.
%
%  On the second step, generate the lower N-1 elements of the second column,
%  compute the Householder matrix H2 that annihilates them,
%  and set A := A * H2' = A * H2 = H1 * H2.
%
%  On the N-1 step, generate the lower 2 elements of column N-1,
%  compute the Householder matrix HN-1 that annihilates them, and
%  and set A := A * H(N-1)' = A * H(N-1) = H1 * H2 * ... * H(N-1).
%  This is our random orthogonal matrix.
%
  for j = 1 : n-1
%
%  Set the vector that represents the J-th column to be annihilated.
%
    x(1:j-1) = 0.0;

    for i = j : n
      [ x(i), seed ] = r8_normal_01 ( seed );
    end
%
%  Compute the vector V that defines a Householder transformation matrix
%  H(V) that annihilates the subdiagonal elements of X.
%
    v = r8vec_house_column ( n, x, j );
%
%  Postmultiply the matrix A by H'(V) = H(V).
%
    a = r8mat_house_axh ( n, a, v );

  end

  return
end


text:
function [data, opt] = HogRawPixelNormExtractor(im, tmpl, opt)

[dataHog] = HogExtractor(im, tmpl, opt);
temp = opt;
temp.FeatureExtractor.tmplsize = temp.FeatureExtractor.tmplsize / 2;
[dataRaw] = RawPixelExtractor(im, tmpl, temp);

data.feat = [dataHog.feat; dataRaw.feat];
dataNorm = sqrt(sum(data.feat .* data.feat));
data.feat = bsxfun(@rdivide, data.feat, dataNorm);
data.tmpl = dataHog.tmpl;

author: flyers

repo: drone-tracking


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function plotGeoverErrorNoRot(decentr_state, fig_num, label)

%% Build struct with rel poses and ground truth
rel_poses = {};
for robot_i = 1:numel(decentr_state)
    r_i_data = decentr_state{robot_i};
    for query_i = 1:numel(r_i_data.place_matches)
        match = r_i_data.place_matches{query_i};
        if (numel(match) > 0)
            rel_pose = [];
            rel_pose.Sim_M_Q = match.Sim_M_Q;
            rel_pose.est_dist = norm(match.Sim_M_Q(1:3, 4));
            gt_p_W_M = ...
                decentr_state{match.robot_i}.gt_T_W_C{...
                match.frame_i}(1:3, 4);
            gt_p_W_Q = r_i_data.gt_T_W_C{query_i}(1:3, 4);
            rel_pose.gt_dist = ...
                norm(gt_p_W_M - gt_p_W_Q);
            rel_poses = [rel_poses; {rel_pose}];
        end
    end
end

%% Eval errors: position
figure(fig_num);
pos_error = cellfun(@(x) norm(x.est_dist - x.gt_dist), rel_poses);
subplot(1, 2, 1);
plot(pos_error);
title([label ': position errors']);
xlabel('rel pose #');
ylabel('position error [m]');

subplot(1, 2, 2);
plot(sort(pos_error), 1:numel(pos_error));
title([label ': position errors'])
xlabel('position error [m]');
ylabel('num rel poses with lower error');
end

author: uzh-rpg


text:
function [ps,params] = prmtrzd_sigmoid(w,input_data)
%
% Algorithm: ps = sigmoid( w0+w1'*input_data ), where
%    w = [w1;w0]; w0 is scalar; and w1 is vector
%
%
% Inputs:
%  w = [w1;w0]; w0 is scalar; and w1 is vector.
%     Use w=[] to let output ps be an MV2DF function handle.
%     If w is a function handle to an MV2DF then ps is the function handle
%     to the composition of w and this function.
%
%  input_data: D-by-T matrix
%
%
% Outputs:
%   ps: function handle (if w=[], or w is handle), or numeric T-by-1
%   params.get_w0(ssat): returns w0 for optimization initialization,
%                        0<ssat<1 is required average sigmoid output.
%   params.tail: is tail of parameter w, which is not consumed by this
%                function.

if nargin==0
    test_this();
    return;
end

m = size(input_data,1);
wsz = m+1;
[whead,wtail] = splitvec_fh(wsz,w);
params.get_w0 = @(ssat) init_w0(ssat);
params.tail = wtail;

y = fusion_mv2df(whead,input_data);
ps = sigmoid_mv2df(y);

    function w0 = init_w0(ssat)
    w0 = zeros(wsz,1);
    offs = logit(ssat);
    w0(end) = offs;
    end

end

function test_this()

D = 3;
K = 5;
data = randn(D,K);

ssat = 0.99;
[sys,params] = prmtrzd_sigmoid([],data);

w0 = params.get_w0(ssat);
test_MV2DF(sys,w0);

ps = prmtrzd_sigmoid(w0,data),

end

author: nesl


text:
function [currentSol, allObjValues, allSolutions] = optimizeCbModelNLP(model, varargin)
% Optimizes constraint-based model using a non-linear objective
%
% USAGE:
%
%    [currentSol, allObjValues, allSolutions] = optimizeCbModelNLP(model, varargin)
%
% INPUT:
%    model:            COBRA model structure
%
% Optional inputs are parameter/value pairs
%
% OPTIONAL INPUT:
%    objFunction:      Name of the non-linear matlab function to be optimized (the
%                      corresponding m-file must be in the current matlab path)
%                      The function receives two arguments, the current flux
%                      vector, and the `NLPProblem` structure.
%    initFunction:     Name of the matlab function used to generate random initial
%                      starting points. The function will be supplied with two
%                      arguments: the model and a cell array of input arguments
%                      (specified in the initArgs parameter)
%    osenseStr:        Optimization direction ('max' or 'min'), this will
%                      override any mention in the model.
%    nOpt:             Number of independent optimization runs performed
%    objArgs:          Cell array of arguments that are supplied to the
%                      objective function as `objArguments` in the NLPProblem
%                      structure (i.e. the second element, will have a field
%                      `objArguments`.)
%    initArgs:         Cell array of arguments to the 'initFunction', will be
%                      provided as second input Argument to the `initFunction`
%    solveroptions:    A Struct with options for the solver used. This is
%                      specific to the solver in question, but the fields should
%                      relate to options accepted by the solver.
%
% OUTPUTS:
%    currentSol:       Solution structure
%    allObjValues:     Array of objective value of each iteration
%    allSolutions:     Array of flux distribution of each iteration
%
% .. Authors:
%       - Markus Herrgard 8/24/07
%       - Modified for new options in solveCobraNLP by Daniel Zielinski 3/19/10
%       - Changed the function to use parameter/value pairs, Thomas Pfau 07/22/17

global CBT_NLP_SOLVER
defaultosenseStr = 'max';
if isfield(model,'osenseStr')
       defaultosenseStr = model.osenseStr;
end
defaultObjFunction = 'NLPobjPerFlux';
defaultInitFunction = 'randomObjFBASol';

p = inputParser;
addRequired(p,'model',@isstruct)
addParamValue(p,'objFunction',defaultObjFunction,@ischar)
addParamValue(p,'initFunction',defaultInitFunction,@ischar)
addParamValue(p,'osenseStr',defaultosenseStr,@(x) strcmp(x,'min') | strcmp(x,'max'));
addParamValue(p,'nOpt',100,@(x) rem(x,1) == 0);
addParamValue(p,'objArgs',{},@iscell)
addParamValue(p,'initArgs',[],@iscell)
addParamValue(p,'solverOptions',[],@isstruct)

parse(p,model,varargin{:});

[osenseStr,objFunction,initFunction,nOpt,objArgs,initArgs,solverOptions] = ...
    deal( p.Results.osenseStr,p.Results.objFunction,p.Results.initFunction, ...
    p.Results.nOpt,p.Results.objArgs,p.Results.initArgs, p.Results.solverOptions);
%To assure, that we don't have any complications, we replace the model
%osenseStr during the processing.
model.osenseStr = osenseStr;

if strcmp(osenseStr,'max')
    osense = -1;
else
    osense = 1;
end
%Set some defaults for fmincon:

if isnumeric(solverOptions)
    if strcmp(CBT_NLP_SOLVER,'matlab')
        if verLessThan('matlab','9')
                solverOptions = struct('TolFun',1e-20);
           else
        solverOptions = struct('FunctionTolerance', 1e-20,'OptimalityTolerance',1e-20);
        end
    else
        solverOptions = struct();
    end

end
%Set default values, if they are not set.
if ~isfield(solverOptions,'iterationLimit')
    solverOptions.iterationLimit = 100000;
end

if ~isfield(solverOptions,'printLevel')
    solverOptions.printLevel = 3; %3 prints every iteration.  1 does a summary.  0 = silent.
end

%The same as above for the objective Function
if strcmp(initFunction,defaultInitFunction) && isnumeric(initArgs)
    solOpt = optimizeCbModel(model,osenseStr);
    %Minimum fraction of the objective function to select start points from
    %is 50% of the maximal objective
    initArgs = {osenseStr,0.5, 0.5 * solOpt.f};
else
    if isnumeric(initArgs)
        initArgs = {};
    end
end

[nMets,nRxns] = size(model.S);
NLPproblem = buildLPproblemFromModel(model);

NLPproblem.objFunction = objFunction;
NLPproblem.objArguments = objArgs;

% Current best solution
currentSol.f = osense*inf;
allObjValues = zeros(nOpt,1);
allSolutions = zeros(nRxns,nOpt);

NLPproblem.user.model = model; %pass the model into the problem for access by the nonlinear objective function

for i = 1:nOpt
    x0 = feval(initFunction,model,initArgs);
    NLPproblem.x0 = x0; %x0 now a cell within the NLP problem structure
    solNLP = solveCobraNLP(NLPproblem,solverOptions); %New function call
    %solNLP = solveCobraNLP(NLPproblem,[],objArgs); Old Code
    fprintf('%d\t%f\n',i,osense*solNLP.obj);
    allObjValues(i) = osense*solNLP.obj;
    allSolutions(:,i) = solNLP.full;
    if osense*solNLP.obj > currentSol.f
       currentSol.f = osense*solNLP.obj;
       currentSol.x = solNLP.full;
       currentSol.stat = solNLP.stat;
    end
end

author: opencobra


text:
function show(fname,opts,leng_str,output)
%SHOW   Display contents of NetCDF file
%
%   Syntax:
%      SHOW(FILE,OPT,LENG_STR,OUTPUT)
%
%   Inputs:
%      FILE       NetCDF file [ <none> ]
%      OPT        Show variables long_name, units, file  dimensions
%                 and attributes [ <1 1 0 0> ]
%      LENG_STR   Maximum length of variable name and long_name
%                 string [ 25 ]
%      OUTPUT     Write (fid or filename)/show results [ {1} <fid> <filename> ]
%
%   Requirement:
%      The NetCDF toolbox available at
%      http://crusty.er.usgs.gov/~cdenham/MexCDF/nc4ml5.html
%
%   Comment:
%      If scale or offset are present for a var, then an '*' is written
%      before it, first '*' is for scale and second is for offset
%
%   Examples:
%      show
%      show('file.nc')
%      show('file.nc',[1 1 1 1],50,'file_contents.txt')
%
%   MMA 2001, martinho@fis.ua.pt
%
%   See also USE

%   Department of Physics
%   University of Aveiro, Portugal

%   26-01-2005 - Small correction to be used also by r14
%   07-02-2005 - Added OPT(3:4)

openf=0;
if nargin < 4
  fid=1;
  fprintf(1,'\n');
else
  if isstr(output)
    fid=fopen(output,'w');
    openf=1;
  else
    fid=output;
  end
end

if nargin < 3
  smax=25;
else
  smax=leng_str;
end

if nargin < 2
  show_lname = 1;
  show_units = 1;
  show_att   = 0;
  show_dim   = 0;
else
  show_lname=opts(1);
  eval('show_units = opts(2);','show_units = 0;');
  eval('show_att   = opts(3);','show_att   = 0;');
  eval('show_dim   = opts(4);','show_dim   = 0;');
end

if nargin == 0
  nc=netcdf('nowrite');
else
  nc=netcdf(fname,'nowrite');
end

if isempty(nc)
  return
end

fname=name(nc);
vv=var(nc);
nvars=length(vv);

% ---------------------------------------------------------- lengths :
warning off
l_Name      = length(name(vv{1}));
l_Long_name = length(nc{name(vv{1})}.long_name(:));
l_maxSize   = length(num2str(max(size(vv{1}))));
l_Units     = length(nc{name(vv{1})}.units(:));

for i=2:nvars
  l_Name_      = length(name(vv{i}));
  l_Long_name_ = length(nc{name(vv{i})}.long_name(:));
  l_maxSize_   = length(num2str(max(size(vv{i}))));
  l_Units_     = length(nc{name(vv{i})}.units(:));

  l_Name      = max(l_Name,l_Name_);
  l_Long_name = max(l_Long_name_,l_Long_name);
  l_maxSize   = max(l_maxSize_,l_maxSize);
  l_Units     = max(l_Units_,l_Units);
end
warning on

LN  = min(smax,l_Name);
LLN = min(smax,l_Long_name);
LU  = min(smax,l_Units);

% deal with size:
l_Size = 0;
for i=1:nvars
  Size      = size(vv{i});
  sformat  = [];
  for n=1:length(Size)
    sformat = [sformat,'  %',num2str(l_maxSize),'d '];
  end
  size_ = sprintf(sformat,Size);
  l_Size = max(l_Size,length(size_));
end
LS  = min(smax,l_Size);

% --------------------------------------------------------------------

fprintf(fid,'-------------------------------------------------------------------\n');
fprintf(fid,'# Contents of the NetCDF file\n');
fprintf(fid,'  %s\n',fname);
fprintf(fid,'%s\n','');

% --------------------------------------------------------------------
% Atributes
% --------------------------------------------------------------------
if show_att
  fprintf(fid,':: Atributes\n');

  strn = '';
  strl = '';
  a=att(nc);
  for i=1:length(a)
    b=a{i};
    attrib.name{i}  = name(b);
    attrib.value{i} = b(:);

    strn = strvcat(strn,attrib.name{i});
    % deal with numeric values:
    if isnumeric(attrib.value{i})
      if any(size(attrib.value{i}) == 1) & length(attrib.value{i}) == 1
        attrib.value{i} = ['<numeric: ',num2str(attrib.value{i}),'>'];
      else
        attrib.value{i} = ['<numeric: size = ',num2str(size(attrib.value{i})),'>'];
      end
    % deal with long strings:
    else
      if length(attrib.value{i}) > 50;
        attrib.value{i} = [attrib.value{i}(1:49),'+'];
      end
    end
    % deal with empty values:
    if isempty(attrib.value{i})
      attrib.value{i} = '<empty>';
    end
    strl = strvcat(strl,num2str(attrib.value{i}));
  end
  maxn = size(strn,2);
  maxl = size(strl,2);
  for i=1:length(attrib.name)
    format = ['   %',num2str(maxn),'s     %s\n'];
    fprintf(fid,format,strn(i,:),rtrim(strl(i,:)));
  end

  %  n_fileatt(fname);
  fprintf(fid,'\n');
end

% --------------------------------------------------------------------
% Dimensions
% --------------------------------------------------------------------
if show_dim
  fprintf(fid,':: Dimensions\n');

  strn = '';
  strl = '';
  d=dim(nc);
  for i=1:length(d)
    thedim.name{i}   = name(d{i}); % NetCDF_Dimension
    thedim.length{i} = d{i}(:);    % itsLength

    strn = strvcat(strn,thedim.name{i});
    strl = strvcat(strl,num2str(thedim.length{i}));
  end
  maxn = size(strn,2);
  maxl = size(strl,2);
  for i=1:length(thedim.name)
    format = ['   %',num2str(maxn),'s   %',num2str(maxl),'s\n'];
    fprintf(fid,format,strn(i,:),rtrim(strl(i,:)));
  end

  %  n_filedim(fname);
  fprintf(fid,'\n');
end

% --------------------------------------------------------------------
% Variables
% --------------------------------------------------------------------
fprintf(fid,':: Variables\n\n');

warning off
for i=1:nvars
  % search for scale and offset:
  is_scale=0;  str1='-';
  is_offset=0; str2='-';
  a=att(vv{i});
  for n=1:length(a)
    b=a{n};
    att_name=name(b);
    if isequal(att_name,'scale_factor');
      is_scale=1;
      str1='*';
    end
    if isequal(att_name,'add_offset');
      is_offset=1;
      str2='*';
    end
  end

  Name      = name(vv{i});
  Long_name = nc{name(vv{i})}.long_name(:);
  Size      = size(vv{i});
  Units     = nc{name(vv{i})}.units(:);

  sname = repmat(' ',1,LN);
  lname = repmat(' ',1,LLN);
  size_ = repmat(' ',1,LS);
  units = repmat(' ',1,LU);

  lN  = min(length(Name),LN);
  lLN = min(length(Long_name),LLN);
  lU  = min(length(Units),LU);

  sname(1:lN)  = Name(1:lN);
  lname(1:lLN) = Long_name(1:lLN);
  units(1:lU)  = Units(1:lU);

  % deal with size:
  sformat  = [];
  for n=1:length(Size)
    sformat = [sformat,'  %',num2str(l_maxSize),'d '];
  end
  size__ = sprintf(sformat,Size);
  lS  = min(length(size__),LS);
  size_(1:lS) = size__(1:lS);

  % add + char at the end of strings:
  if length(Name) > lN,       sname(end) = '+'; end
  if length(Long_name) > lLN, lname(end) = '+'; end
  if length(size__) > lS,     size_(end) = '+'; end
  if length(Units) > lU,      units(end) = '+'; end

  if show_lname & show_units
    if i==1, fprintf(fid,'%s\n\n','[scale(*) offset(*) - Var - Long_Name - Size - Units]'); end
    str = sprintf('%s%s %s  %s  %s  %s',str1,str2,sname,lname,size_,units);
  elseif show_lname
    if i==1, fprintf(fid,'%s\n\n','[scale(*) offset(*) - Var - Long_Name - Size]');         end
    str = sprintf('%s%s %s  %s  %s',str1,str2,sname,lname,size_);
  elseif show_units
    if i==1, fprintf(fid,'%s\n\n','[scale(*) offset(*) - Var -  Size - Units]');            end
    str = sprintf('%s%s %s  %s  %s',str1,str2,sname,size_,units);
  else
    if i==1, fprintf(fid,'%s\n\n','[scale(*) offset(*) - Var - Size]');                     end
    str = sprintf('%s%s %s %s',str1,str2,sname,size_);
  end
  fprintf(fid,'%s\n',rtrim(str));
end
warning on

nc=close(nc);
fprintf(fid,'\n-------------------------------------------------------------------\n');
if openf
  fid=fclose(fid);
else
  fprintf(1,'\n');
end

function theString = rtrim(theString,charList)
if nargin < 2
  charList = ' ';
end
n=length(theString);
while 1
  if ~any(theString(n) == charList), break, end
  n=n-1;
  if n==0, break, end
end
theString = theString(1:n);

author: Sable


text:
%%
% Compile all HTML files.

a = dir('test_*.m');
opts.outputDir = 'html';

for i=1:length(a)
    name = a(i).name(1:end-2);
	publish(name,opts);
end

author: gpeyre

repo: matlab-toolboxes


text:
function [ boolout ] = isnisar( filename )
%ISNISAR Is NASA-ISRO SAR Mission SLC dataset in HFD5 file format
%
% Written by: Wade Schwartzkopf, NGA/R
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

fid = fopen(filename,'r','b');
boolout = isequal(fread( fid, 8, 'uint8' )',[137 72 68 70 13 10 26 10]); % HDF5 signature
fclose(fid);
% Now that we know it is HDF5, is it NISAR SLC HDF5?
if boolout
    try % Error if root attribute 'mission_name' does not exist
        boolout=strcmpi(h5readatt(filename,'/','mission_name'), 'NISAR') && ...
            strcmpi(deblank(h5read(filename,'/science/LSAR/identification/productType')), 'SLC');
    catch % h5readatt replaced hdf5read at some point in MATLAB history
        try
            boolout=strcmpi(hdf5read(filename,'/','mission_name'), 'NISAR') && ...
                strcmpi(deblank(hdf5read(filename,'/science/LSAR/identification/productType')), 'SLC');
        catch % Is there a more direct way to check for attribute existence first?
            boolout=false;
        end
    end
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [M] = spm_ADEM_M_set(M)
% Set indices and perform checks on hierarchical action models
% FORMAT [M] = spm_ADEM_M_set(M)
%
% for each level (i); required fields
%
%   M(i).g  = y(t)  = g(x,v,a,P)    {inline function, string or m-file}
%   M(i).f  = dx/dt = f(x,v,a,P)    {inline function, string or m-file}
%
% and
%
%   M(i).m  = number of inputs v(i + 1);
%   M(i).n  = number of states x(i);
%   M(i).l  = number of output v(i);
%   M(i).k  = number of action a(i);
%
% or
%
%   M(i).x  = hidden states;
%   M(i).v  = causal states;
%   M(i).a  = action states;
%
% for each level (i); optional fields
%
%   M(i).pE = prior expectation of p model-parameters
%   M(i).V  = precision (input noise)
%   M(i).W  = precision (state noise)
%   M(i).U  = precision (action)
%
%
% sets fields, checks internal consistency of model specification and sets
% estimation parameters.  If (V,W) are not specified infinite precision is
% assumed.
%--------------------------------------------------------------------------
%
%   M(1).E.s;     = smoothness (s.d. in time bins)
%   M(1).E.d;     = embedding order q(v)  (i.e., number of derivatives)
%   M(1).E.n;     = embedding order q(x)
%
% If the highest level involves any dynamic or static transformation
% of its inputs a further level is added with flat priors
%__________________________________________________________________________
% Copyright (C) 2008-2014 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_ADEM_M_set.m 5962 2014-04-17 12:47:43Z spm $

% order
%--------------------------------------------------------------------------
g      = length(M);

% set missing fields
%==========================================================================

% check for specification of hidden states
%--------------------------------------------------------------------------
if isfield(M,'f') && ~isfield(M,'n') && ~isfield(M,'x')
    msgbox('please specify hidden states or their number')
    error(' ')
end

% check supra-ordinate level and add one (with flat priors) if necessary
%--------------------------------------------------------------------------
try
    fcnchk(M(g).g);
    g      = g + 1;
    M(g).l = M(g - 1).m;
end
M(g).m = 0;
M(g).n = 0;

% default fields for static models (hidden states)
%--------------------------------------------------------------------------
if ~isfield(M,'f')
    [M.f] = deal(inline('sparse(0,1)','x','v','a','P'));
    [M.x] = deal(sparse(0,1));
    [M.n] = deal(0);
end
for i  = 1:g
    try
        fcnchk(M(i).f);
    catch
        M(i).f = inline('sparse(0,1)','x','v','a','P');
        M(i).x = sparse(0,1);
    end
end

% consistency and format check on states, parameters and functions
%==========================================================================

% prior expectation of parameters M.pE
%--------------------------------------------------------------------------
try
    M.pE;
catch
    % Assume fixed parameters
    %----------------------------------------------------------------------
    for i = 1:g
        M(i).pE = sparse(0,0);
    end
end

% get inputs
%--------------------------------------------------------------------------
try
    v  = M(g).v;
catch
    v  = sparse(0,1);
end
if isempty(v)
    try
        v = sparse(M(g - 1).m,1);
    end
end
if isempty(v)
    try
        v = sparse(M(g).l,1);
    end
end
M(g).l    = length(spm_vec(v));
M(g).v    = v;

% ensure action is specified
%--------------------------------------------------------------------------
for i = 1:g
    try
        a  = M(i).a;
    catch
        a  = sparse(0,1);
    end
    if isempty(a)
        try
            a = sparse(M(i).k,1);
        end
    end
    M(i).k    = length(spm_vec(a));
    M(i).a    = a;
end

% check functions
%--------------------------------------------------------------------------
for i = (g - 1):-1:1
    try
        x = M(i).x;
    catch
        x = sparse(M(i).n,1);
    end
    if isempty(x) && M(i).n
        x = sparse(M(i).n,1);
    end

    % check f(x,v,P)
    %----------------------------------------------------------------------
    try
        M(i).f = fcnchk(M(i).f);
        if nargin(M(i).f) ~= 4
            M(i).f = inline(char(M(i).f),'x','v','a','P');
        end
    end
    try
        f      = feval(M(i).f,x,v,a,M(i).pE);
        if length(spm_vec(x)) ~= length(spm_vec(f))
            errordlg(sprintf('please check nargout: G(%i).f(x,v,a,P)',i));
        end
    catch
        errordlg(sprintf('evaluation failure: G(%i).f(x,v,a,P)',i))
    end

    % check g(x,v,P)
    %----------------------------------------------------------------------
    try
        M(i).g = fcnchk(M(i).g);
        if nargin(M(i).g) ~= 4
            M(i).g = inline(char(M(i).g),'x','v','a','P');
        end
    end
    try
        M(i).m = length(spm_vec(v));
        v      = feval(M(i).g,x,v,a,M(i).pE);
        a      = M(i).a;
        M(i).k = length(spm_vec(a));
        M(i).l = length(spm_vec(v));
        M(i).n = length(spm_vec(x));

        M(i).a = a;
        M(i).v = v;
        M(i).x = x;

    catch
        errordlg(sprintf('evaluation failure: G(%i).g(x,v,a,P)',i))
    end
end

% remove empty levels
%--------------------------------------------------------------------------
try
    g  = find(~spm_vec(M.m),1);
    M  = M(1:g);
catch
    errordlg('please specify number of variables')
end

% number of x (hidden states)
%--------------------------------------------------------------------------
nx     = sum(spm_vec(M.n));

% check precisions
%==========================================================================
try, M.V;  catch, M(1).V  = []; end
try, M.W;  catch, M(1).W  = []; end

% check precisions
%--------------------------------------------------------------------------
for i = 1:g

    % check V and assume unit precision if improperly specified
    %----------------------------------------------------------------------
    if isvector(M(i).V), M(i).V = diag(M(i).V); end
    if length(M(i).V) ~= M(i).l
        try
            M(i).V = speye(M(i).l,M(i).l)*M(i).V(1);
        catch
            M(i).V = speye(M(i).l,M(i).l)*exp(32);
        end
    end

    % check W and assume unit precision if improperly specified
    %----------------------------------------------------------------------
    if isvector(M(i).W), M(i).W = diag(M(i).W); end
    if length(M(i).W) ~= M(i).n
        try
            M(i).W = speye(M(i).n,M(i).n)*M(i).W(1);
        catch
            M(i).W = speye(M(i).n,M(i).n)*exp(32);
        end
    end
end

% check restiction of precision for action (gain)
%==========================================================================
try, M(1).U;  catch, M(1).U = []; end

if isvector(M(1).U), M(1).U = diag(M(1).U); end
if length(M(1).U) ~= M(1).l
    try
        M(1).U = speye(M(1).l,M(1).l)*M(1).U(1);
    catch
        M(1).U = speye(M(1).l,M(1).l)*exp(2);
    end
end

% estimation parameters M(1).E.s, n,...
%==========================================================================
% E.s;                               % smoothness (seconds)
% E.dt;                              % time step
% E.d;                               % approximation order of q(x,v)
% E.n;                               % order of embedding (n >= d)

% temporal smoothness - s.d. of kernel
%--------------------------------------------------------------------------
try M(1).E.s;  catch, if nx, M(1).E.s = 1/2; else M(1).E.s = 0; end, end

% time step
%--------------------------------------------------------------------------
try M(1).E.dt; catch, M(1).E.dt = 1; end

% embedding orders
%--------------------------------------------------------------------------
try M(1).E.d;  catch, if nx, M(1).E.d = 2;  else M(1).E.d = 0;  end, end
try M(1).E.n;  catch, if nx, M(1).E.n = 6;  else M(1).E.n = 0;  end, end

M(1).E.d = min(M(1).E.d,M(1).E.n);

% checks on smoothness hyperparameter
%==========================================================================
for i = 1:g

    try, M(i).sv; catch,   M(i).sv = M(1).E.s; end
    try, M(i).sw; catch,   M(i).sw = M(1).E.s; end

    if ~isscalar(M(i).sv), M(i).sv = M(1).E.s; end
    if ~isscalar(M(i).sw), M(i).sw = M(1).E.s; end
end

author: spm


text:
%VAR Dataset overload
%
%   [V,U] = VAR(A,W)
%
% Computes variance V and mean U in a single run for consistency with datafile overload.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit