text:
function [ffun,flag] = limgrad(edge,elen,ffun,dfdx,imax)
%LIMGRAD impose "gradient-limits" on a function defined over
%an undirected graph.
%   [FNEW] = LIMGRAD(EDGE,ELEN,FFUN,DFDX,ITER) computes a
%   "gradient-limited" function FNEW on the undirected graph
%   {EDGE,ELEN}, where EDGE is an NE-by-2 array of edge ind-
%   ices, and ELEN is an NE-by-1 array of edge lengths.
%   Gradients are limited over the graph edges, such that
%
%       ABS(FNEW(N2)-FNEW(N1)) <= ELEN(II) * DFDX,
%
%   where N1=EDGE(II,1) and N2=EDGE(II,2) are the two nodes
%   in the II-TH edge. An iterative algorithm is used, swee-
%   ping over an "active-set" of graph edges until converge-
%   nce is achieved. A maximum of IMAX iterations are done.
%
%   [FNEW,FLAG] = LIMGRAD(...) also returns a boolean FLAG,
%   with FLAG=TRUE denoting convergence.
%
%   See also LIMHFN2

%   Darren Engwirda : 2017 --
%   Email           : engwirda@mit.edu
%   Last updated    : 18/04/2017

%---------------------------------------------- basic checks
    if ( ~isnumeric(edge) || ...
         ~isnumeric(elen) || ...
         ~isnumeric(ffun) || ...
         ~isnumeric(dfdx) || ...
         ~isnumeric(imax) )
        error('limgrad:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end
%---------------------------------------------- basic checks
    if (ndims(edge) ~= +2 || ...
        ndims(elen)  > +2 || ...
        ndims(ffun)  > +2 || ...
        numel(dfdx) ~= +1 || ...
        numel(imax) ~= +1 )
        error('limgrad:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(edge,2) < +2 || ...
        size(elen,2)~= +1 || ...
        size(ffun,2)~= +1 || ...
        size(edge,1)~= size(elen,1) )
        error('limgrad:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    nnod = size(ffun,1) ;

%---------------------------------------------- basic checks
    if (dfdx < +0. || imax < +0)
        error('limgrad:invalidInputArgument', ...
            'Invalid input parameter.');
    end
    if (min(min(edge(:,1:2))) < +1 || ...
            max(max(edge(:,1:2))) > nnod )
        error('limgrad:invalidInputArgument', ...
            'Invalid EDGE input array.') ;
    end

%-- IVEC(NPTR(II,1):NPTR(II,2)) are edges adj. to II-TH node
    nvec = [edge(:,1); edge(:,2)];
    ivec = [(1:size(edge,1))'; ...
            (1:size(edge,1))'] ;

   [nvec,pidx] = sort (nvec) ;
    ivec       = ivec (pidx) ;

    mark = false(nnod,1) ;
    mark(edge(:,1)) = true ;
    mark(edge(:,2)) = true ;

    idxx = find(diff(nvec) > +0) ;

    nptr = zeros(nnod,2) ;
    nptr(:,2) = -1 ;
    nptr(mark,1) = [+1; idxx+1];
    nptr(mark,2) = [idxx; nnod];

%----------------------------- ASET=ITER if node is "active"
    aset = zeros(size(ffun,1),1) ;

%----------------------------- exhaustive 'til all satisfied
    ftol = min(ffun) * sqrt(eps) ;

    for iter = +1 : imax

    %------------------------- find "active" nodes this pass
        aidx = find(aset == iter - 1) ;

        if (isempty(aidx)), break; end

    %------------------------- reorder => better convergence
       [~,idxx] = sort(ffun(aidx)) ;

        aidx = aidx(idxx);

    %------------------------- visit adj. edges and set DFDX
        for ipos = 1 : length(aidx)
            npos = aidx(ipos) ;
            for jpos = nptr(npos,1) ...
                     : nptr(npos,2)

                epos = ivec(jpos,1) ;

                nod1 = edge(epos,1) ;
                nod2 = edge(epos,2) ;

            %----------------- calc. limits about min.-value
                if (ffun(nod1) > ffun(nod2))

                fun1 = ffun(nod2) ...
                     + elen(epos) * dfdx ;

                if (ffun(nod1) > fun1+ftol)
                    ffun(nod1) = fun1;
                    aset(nod1) = iter;
                end

                else

                fun2 = ffun(nod1) ...
                     + elen(epos) * dfdx ;

                if (ffun(nod2) > fun2+ftol)
                    ffun(nod2) = fun2;
                    aset(nod2) = iter;
                end

                end

            end
        end

    end

    flag = (iter < imax) ;

end

author: CHLNDDEV


text:
%                                              x = frameit(lab,frmq,frms,K)
% FRAMEIT
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

function x = frameit(lab,frmq,frms,K)
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%
x = [lab(1:K.l); qframeit(lab,frmq,K); psdframeit(lab,frms,K)];

repo: REST


text:
function varargout=zerotrim(varargin)
%
% For calling details please see v_zerotrim.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_zerotrim.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: zerotrim.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_zerotrim(varargin{:});
else
    v_zerotrim(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function model = svargplvmExpandParam(model, params)

% SVARGPLVMEXPANDPARAM Expand a parameter vector into a shared variational GP-LVM model.
% FORMAT
% DESC takes a SVARGPLVM structure and a vector of parameters, and
% fills the structure with the given parameters. Also performs any
% necessary precomputation for likelihood and gradient
% computations, so can be computationally intensive to call.
% Parallelism can be applied w.r.t the models, so that some speed-up is achieved on
% multi-core machines.

% Parameters must be passed as a vector in the following order (left to right)
% - parameter{size} -
% vardistParams{model.vardist.nParams} % mu, S
%       OR
% [dynamicsVardistParams{dynamics.vardist.nParams} dynamics.kernParams{dynamics.kern.nParams}] % mu_bar, lambda
%    % Followed by:
% private params of 1st model {model.comp{1}.nPrivateParams}
% private params of 2nd model {model.comp{2}.nPrivateParams}
%           ...
% private params of i-th model {model.comp{i}.nPrivateParams}
%
% ARG model : the svargplvm model to update with parameters
% ARG params : parameter vector
% RETURN model : model with updated parameters
%
% SEEALSO : svargplvmCreate, svargplvmExtractParam, modelExpandParam
%
% COPYRIGHT : Andreas C. Damianou, 2011

% VARGPLVM

try
    pool_open = matlabpool('size')>0;
catch e
    pool_open = 0;
end

if pool_open && (isfield(model,'parallel') && model.parallel)
    model = svargplvmExpandParamPar(model, params);
else
    model = svargplvmExpandParamOrig(model, params);
end
%!!! Since all dynamical models have the same "dynamics" structure, we can
%avoid the expandParam for ALL the dynamics structures. In fact, we should
%not store all these dynamics structures at all.

function model = svargplvmExpandParamPar(model, params)
%startVal = 1;
%endVal = model.N*model.q;
%paramsX = params(startVal:endVal);
% model.barmu = reshape(params(startVal:endVal), model.N, model.q);
startValS=1;

if isfield(model, 'dynamics') & ~isempty(model.dynamics)
    endValS = model.dynamics.nParams;
else
    endValS = model.vardist.nParams;
end
sharedParams = params(startValS:endValS);

startVal{1} = endValS+1;
endVal{1} = startVal{1} + model.comp{1}.nPrivateParams-1;
for i = 2:model.numModels
    % model.comp{i}.X = model.X;
    startVal{i} = endVal{i-1}+1;
    endVal{i} = startVal{i} + model.comp{i}.nPrivateParams-1;
end

modelTempComp = model.comp;
parfor i = 1:model.numModels
    % model.comp{i}.X = model.X;
   % startVal = endVal+1;
   % endVal = startVal + model.comp{i}.nPrivateParams-1;
    params_i = [sharedParams params(startVal{i}:endVal{i})];
    modelTempComp{i} = vargplvmExpandParam(modelTempComp{i}, params_i);
    if strcmp(modelTempComp{i}.kern.type, 'rbfardjit') || ~iscell(modelTempComp{i}.kern)
        modelTempComp{i}.kern.comp{1}.inputScales = modelTempComp{i}.kern.inputScales;
    end
    %[par, names] = vargplvmExtractParam(model.comp{i}); %%%%5
    %for j=5:-1:1 %%
    %params_i(end-j) %%%
    %end%%
    %model.comp{i}.kern.comp{1}.inputScales
end
model.comp = modelTempComp;

% with vargplvmExpandParam for the submodels, and given that the dynamics
% parameters (if there are dynamics) are shared, barmu has been turned into
% mu and this is the same for all submodels. The same holds for the
% variational distribution.
model.X = model.comp{1}.X;

if isfield(model, 'dynamics') & ~isempty(model.dynamics)
    model.dynamics = model.comp{1}.dynamics; % = model.comp{i}.dynamics for all submodels m
end

model.vardist = model.comp{1}.vardist;

function model = svargplvmExpandParamOrig(model, params)

%startVal = 1;
%endVal = model.N*model.q;
%paramsX = params(startVal:endVal);
% model.barmu = reshape(params(startVal:endVal), model.N, model.q);

startVal=1;

if isfield(model, 'dynamics') & ~isempty(model.dynamics)
    endVal = model.dynamics.nParams;
else
    endVal = model.vardist.nParams;
end
sharedParams = params(startVal:endVal);

for i = 1:model.numModels
    % model.comp{i}.X = model.X;
    startVal = endVal+1;
    endVal = startVal + model.comp{i}.nPrivateParams-1;
    params_i = [sharedParams params(startVal:endVal)];
    model.comp{i} = vargplvmExpandParam(model.comp{i}, params_i);
    % For compatibility ...
    if strcmp(model.comp{i}.kern.type, 'rbfardjit')
        model.comp{i}.kern.comp{1}.inputScales = model.comp{i}.kern.inputScales;
    end
    %[par, names] = vargplvmExtractParam(model.comp{i}); %%%%5
    %for j=5:-1:1 %%
    %params_i(end-j) %%%
    %end%%
    %model.comp{i}.kern.comp{1}.inputScales
end

% with vargplvmExpandParam for the submodels, and given that the dynamics
% parameters (if there are dynamics) are shared, barmu has been turned into
% mu and this is the same for all submodels. The same holds for the
% variational distribution.
model.X = model.comp{1}.X;

if isfield(model, 'dynamics') & ~isempty(model.dynamics)
    model.dynamics = model.comp{1}.dynamics; % = model.comp{i}.dynamics for all submodels m
end

model.vardist = model.comp{1}.vardist;

author: SheffieldML

repo: vargplvm


text:
function manageMiji(varargin)
	% Controls starting and stopping of Miji in a way to reduce issues
	% Biafra Ahanonu
	% started: 2019.04.08
	% inputs
		%
	% outputs
		%

	ciapkg.io.manageMiji('passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [tissueModel, coreRxn, nonCoreRxn,	zeroExpRxns, pruneTime, cRes] = mCADRE(model, ubiquityScore, confidenceScores, protectedRxns, checkFunctionality, eta, tol)
% Uses the mCADRE algorithm (`Wang et al., 2012`) to extract a context
% specific model using expression data. mCADRE algorithm defines a set of core
% reactions and prunes all other reactions based on their expression,
% connectivity to core and confidence score. The removal of reactions is not
% necessary performed to support the core of defined functionalities. Core
% reactions are only removed if supported by a certain number of
% zero-expression reactions.
%
% USAGE:
%   tissueModel = mCADRE(model, ubiquityScore, confidenceScores)
%
% INPUTS:
%    model:                 input model (COBRA model structure)
%    ubiquityScore:         ubiquity scores corresponding to genes
%                           in gene_id quantify how often a gene is expressed accross samples.
%                           (default values defined in function of threshold_high value)
%    confidenceScores:      literature-based evidence for generic model
%                           (default value = 0)
%
% OPTIONAL INPUTS:
%    protectedRxns:         cell with reactions names that are manually added to
%                           the core reaction set (i.e. {'Biomass_reaction'})
%    checkFunctionality:    Boolean variable that determine if the model should be able
%                           to produce the metabolites associated with the protectedRxns
%
%                           0: don't use functionality check (default value)
%
%                           1: include functionality check
%    eta:                   tradeoff between removing core and zero-expression
%                           reactions (default value: 1/3)
%    tol:                   minimum flux threshold for "expressed" reactions
%                           (default 1e-8)
%
% OUTPUTS:
%    tissueModel:           pruned, context-specific model
%    coreRxn:               core reactions in model
%    nonCoreRxn:            non-core reactions in model
%    zeroExpRxns:           reactions with zero expression (i.e., measured zero, not just
%                           missing from expression data)
%    pruneTime:             total reaction pruning time
%    cRes:                  result of model checks (consistency/function)
%                           - vs. +: reaction `r` removed from generic model or not
%                           1 vs. 2: reaction `r` had zero or non-zero expression evidence
%                           `-x.y` : removal of reaction `r` corresponded with removal of `y` (num.) total
%                           core reactions
%                           `+x.1` vs. `x.0` : precursor production possible after removal of
%                           reaction `r` or not
%                           3: removal of reaction r by itself prevented production of required
%                           metabolites (therefore was not removed)
%
%
% `Wang et al. (2012). Reconstruction of genome-scale metabolic models for
% 126 human tissues using mCADRE. BMC Syst. Biol. 6, 153.`
%
% Authors: - This script is an adapted version of the implementation from
%            https://github.com/jaeddy/mcadre.
%          - Modified and commented by S. Opdam and A. Richelle,May 2017

    if nargin < 4 || isempty(protectedRxns)
        protectedRxns = [];
    end
    if nargin < 5 || isempty(checkFunctionality)
        checkFunctionality=0; %don't use functionality check
    end
    if nargin < 6 || isempty(eta)
        eta = 1/3;% default value for tradeoff
    end
    if nargin < 7 || isempty(tol)
        tol = 1e-8;% default tolerance
    end

    disp('Processing inputs and ranking reactions...')

    [coreRxn, nonCoreRxn, rankNonCore, zeroExpRxns] = rankReactions(model, ubiquityScore, confidenceScores, protectedRxns);

    % Precursor are defined as the metabolites involved in the
    % reactions defined in protected reactions
    if checkFunctionality == 1 && ~isempty(protectedRxns)
        precursorMets = model.mets(any(model.S(:, findRxnIDs(model, protectedRxns)) < 0, 2));

        genericStatus = checkModelFunction(model, precursorMets);

    else
    	precursorMets={};
    	genericStatus=1;
    end

    if genericStatus
    	disp('Generic model passed precursor metabolites test');

        disp('Pruning reactions...')
        t0 = clock;

        [tissueModel, cRes] = pruningModel(model, rankNonCore, coreRxn, zeroExpRxns, precursorMets, eta, tol);
        pruneTime = etime(clock,t0);
    else
    	error('Generic model failed precursor metabolites test')
    end

    tissueModel = removeUnusedGenes(tissueModel);

end

author: opencobra

repo: cobratoolbox


text:
function centroid = polygon_centroid_2 ( n, v )

%*****************************************************************************80
%
%% POLYGON_CENTROID_2 computes the centroid of a polygon.
%
%  Discussion:
%
%    The centroid is the area-weighted sum of the centroids of
%    disjoint triangles that make up the polygon.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%
%    Input, real V(2,N), the coordinates of the vertices.
%
%    Output, real CENTROID(2,1), the coordinates of the centroid.
%
  area = 0.0;
  centroid(1:2,1) = 0.0;

  for i = 1 : n - 2

    area_triangle = triangle_area ( ...
      v(1,i),   v(2,i), ...
      v(1,i+1), v(2,i+1), ...
      v(1,n),   v(2,n) );

    area = area + area_triangle;

    centroid(1:2,1) = centroid(1:2,1) + area_triangle ...
      * ( v(1:2,i) + v(1:2,i+1) + v(1:2,n) ) / 3.0;

  end

  if ( area == 0.0 )
    centroid(1:2,1) = v(1:2,1);
  else
    centroid(1:2,1) = centroid(1:2,1) / area;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function chunk = utl_buffer(chunk,buffer,desired_length)
% Append a chunk to a buffer and trim the buffer to the desired length.
% Buffer = utl_buffer(Chunk,Buffer,DesiredLength)
%
% In:
%   Chunk : a chunk of data (continuous EEGLAB dataset struct)
%
%   Buffer : a buffer of previous data (continuous EEGLAB dataset struct)
%            each non-empty time-series field present in the buffer must have the same number of
%            channels as the corresponding field in the chunk
%
%   DesiredLength : the number of samples of desired output
%
% Out:
%   Buffer : the updated and trimmed buffer
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-11-12

% input validation
utl_check_fields(chunk,{'data','event','xmax','srate'},'chunk','signal');
utl_check_fields(buffer,{'data','event'},'buffer','signal');
if ~isempty(chunk.event)
    if ~isfield(chunk.event,'latency')
        error('The chunk is missing the required .event.latency field.'); end
    latency_numels = cellfun('prodofsize',{chunk.event.latency});
    if any(latency_numels == 0)
        error('The given chunk has one or more events with empty .latency field. This is not permitted.');
    elseif any(latency_numels ~= 1)
        error('The given chunk has one or more events with a .latency value that is not a scalar. This is not permitted.');
    end
end

% concatenate markers if necessary
if (desired_length ~= size(chunk.data,2)) && (~isempty(buffer.event) || ~isempty(chunk.event))
    if size(chunk.data,2) > desired_length
        if ~isempty(chunk.event)
            % shift the latency of the markers based on how many samples we have to cut
            latency = [chunk.event.latency] - (size(chunk.data,2)-desired_length);
            [chunk.event.latency] = arraydeal(latency);
            % remove excess markers
            chunk.event(latency<0.5) = [];
        end
    else
        if ~isempty(chunk.event)
            % shift the latency of the chunk's markers based on how many samples we have to prepend
            [chunk.event.latency] = arraydeal([chunk.event.latency]+min(size(buffer.data,2),(desired_length-size(chunk.data,2))));
        end
        if ~isempty(buffer.event)
            % shift the latency of the buffer's markers based on how many samples we cut from that
            latency = [buffer.event.latency] - max(0,(size(buffer.data,2)+size(chunk.data,2)-desired_length));
            [buffer.event.latency] = arraydeal(latency);
            % concatenate markers
            if isempty(chunk.event)
                chunk.event = buffer.event(latency>=0.5);
            else
                try
                    chunk.event = [buffer.event(latency>=0.5) chunk.event];
                catch e
                    fprintf('WARNING: Trying to concatenate incompatible event structs in the buffer (%s) and the incoming chunk (%s); dropping buffer contents; error message: %s\n',hlp_tostring(fieldnames(buffer.event)),hlp_tostring(fieldnames(chunk.event)),e.message);
                end
            end
        end
    end
end

% for each time-series field in the chunk...
for fld = utl_timeseries_fields(chunk)
    field = fld{1};

    % make sure that we have the requested field in the previous buffer to avoid special cases below
    if ~isfield(buffer,field)
        buffer.(field) = []; end

    % if some other amount of data than what's in the chunk was requested
    if desired_length ~= size(chunk.(field),2)
        if size(chunk.(field),2) < desired_length
            try
                if size(buffer.(field),2) == desired_length
                    % we can do an in-place update without reallocations
                    chunk.(field) = cat(2,buffer.(field)(:,size(chunk.(field),2)+1:end,:,:,:,:,:,:),chunk.(field));
                else
                    % append new samples & cut excess data
                    chunk.(field) = cat(2,buffer.(field),chunk.(field));
                    if size(chunk.(field),2) > desired_length
                        chunk.(field) = chunk.(field)(:,end-desired_length+1:end,:,:,:,:,:,:); end
                end
            catch e
                error('Error trying to concatenate time-series field .%s of buffer (size=%s) and incoming chunk (size=%s); error message: %s.',field,hlp_tostring(size(buffer.(field))),hlp_tostring(size(chunk.field)),e.message);
            end
        else
            % if the chunk is longer than what's requested cut the excess data
            chunk.(field) = chunk.(field)(:,end-desired_length+1:end,:,:,:,:,:,:);
        end
    end
end

% update time-related meta-data
chunk.pnts = size(chunk.data,2);
chunk.xmin = chunk.xmax - (chunk.pnts-1)/chunk.srate;


text:
function ord = plot_vpath (vpath,T,continuous,order_trials,behaviour,cm)

if nargin < 2 || isempty(T), T = size(vpath,1); end
if nargin < 3, continuous = length(T)==1 || any(T(1)~=T); end % show it as continuous data?
if nargin < 4, order_trials = false; end % order the trials?
if nargin < 5, behaviour = []; end % behaviour with respect to which order the trials
if nargin < 6, cm = colormap; end % colormap

Gamma = vpath_to_stc(vpath);
ord = plot_Gamma (Gamma,T,continuous,order_trials,behaviour,cm);

end

author: OHBA-analysis

repo: HMM-MAR


text:
function test_suite=test_moxunit_util_is_message_identifier
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_moxunit_util_is_message_identifier_non_char
    aet=@(varargin)assertExceptionThrown(@()...
                    moxunit_util_is_message_identifier(varargin{:}),'');

    aet(struct());                      % struct
    aet({'id:id2'});                    % cell
    aet(2);                             % numeric
    aet(false);                         % boolean
    aet(true);                          % boolean
    aet([]);                          % boolean

function test_moxunit_util_is_message_identifier_true
    true_strs={'a:b',...
                'a___:b___',...
                'abc:def:ghi:jkl:mno',...
                'a_b:c_123',...
                };

    for k=1:numel(true_strs)
        s=true_strs{k};
        assertTrue(moxunit_util_is_message_identifier(s),s);
    end

function test_test_moxunit_util_is_message_identifier_false
    false_strs={'',...                  % empty string
                'a',...                 % missing colon
                'ab',...                % missing colon
                'ab: cd',...            % contains space
                sprintf('ab:c\td'),...  % contains tabs
                '_a:b',...              % starts with underscore
                'a:_b',...              % starts with underscore
                '1ab:def',...           % has a number
                ':abc:def:',...         % starts with colon
                'abc:def:',...          % ends with colon
                'a-b',...               % no colon
                };

    for k=1:numel(false_strs)
        s=false_strs{k};

        assertFalse(moxunit_util_is_message_identifier(s),s);
    end

function test_test_moxunit_util_is_message_identifier_octave
    % in octave, the '-' is allowed to be in message identifers;
    % in Matlab this is not allowed.
    % For compatibility we allow the '-' to be allowed in message
    % identifiers
    true_strs={'a-b:b:cdef-ghij',...
                'a-b:c',...
                'a:b-c',...
                };

    for k=1:numel(true_strs)
        s=true_strs{k};
        assertTrue(moxunit_util_is_message_identifier(s),s);
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [isPos,reportStr]=istriequal(v1Mat,f1Mat,v2Mat,f2Mat,maxTol)
% ISTRIEQUAL checks if two triangulations are equal. Triangulations
% resulting from permunations of edge directions and vertices in
% faces are considered equal
%
% Input:
%   regular:
%       v1Mat: double[n1Verts,3] - vertices of the first triangualtion
%       f1Mat: double[n1Faces,3] - faces of the first triangulation
%       v2Mat: double[n1Verts,3] - vertices of the second triangulation
%       f2Mat: double[n1Faces,3] - faces of the second triangulation
%       maxTol: double[1,1] - tolerances used for comparing vertex
%          coordinates
%
% Output:
%   isPos: logical[1,1] - specifies if result of comparison is true
%   reportStr: char[1,] - describes a reason of negative result
%
% $Author: Peter Gagarinov, PhD  <pgagarinov@gmail.com> $
% $Copyright: Peter Gagarinov, PhD,
%            Moscow State University,
%            Faculty of Computational Mathematics and Computer Science,
%            System Analysis Department 2011-2016 $
%
n1Verts=size(v1Mat,1);
n2Verts=size(v2Mat,1);
isPos=n1Verts==n2Verts;
%
if isPos
    n1Faces=size(f1Mat,1);
    n2Faces=size(f2Mat,1);
    isPos=n1Faces==n2Faces;
    if isPos
        [v1Mat,~,indF1Vec]=sortrowstol(v1Mat,maxTol);
        [v2Mat,~,indF2Vec]=sortrowstol(v2Mat,maxTol);
        %
        f1Mat=indF1Vec(f1Mat);
        f2Mat=indF2Vec(f2Mat);
        realTol=max(max(abs(v1Mat-v2Mat)));
        isPos=realTol<=maxTol;
        %
        if isPos
            nF1Unique=size(unique(f1Mat,'rows'),1);
            nF2Unique=size(unique(f2Mat,'rows'),1);
            isPos=nF1Unique==nF2Unique;
            if isPos
                isPos=all(isface(v2Mat,f2Mat,f1Mat));
                if isPos
                    reportStr='';
                else
                    reportStr='faces are different';
                end
            else
                reportStr='numbers of unique faces are different';
            end
        else
            reportStr=sprintf(...
                'vertices are different, real tol=%f, exp tol=%f',...
                realTol,maxTol);
        end
    else
        reportStr='number of faces is different';
    end
else
    reportStr='numbers of vertices are different';
end

repo: OpenMAS


text:
classdef UF2 < PROBLEM
% <multi> <real> <large/none>
% Unconstrained benchmark MOP

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, S. Zhao, P. N. Suganthan, W. Liu, and S. Tiwari,
% Multiobjective optimization test instances for the CEC 2009 special
% session and competition, School of CS & EE, University of Essex, Working
% Report CES-487, 2009.
%------------------------------- 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
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = [0,zeros(1,obj.D-1)-1];
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            J1 = 3 : 2 : obj.D;
            J2 = 2 : 2 : obj.D;
            Y       = zeros(size(X));
            X1      = repmat(X(:,1),1,length(J1));
            Y(:,J1) = X(:,J1)-(0.3*X1.^2.*cos(24*pi*X1+4*repmat(J1,size(X,1),1)*pi/obj.D)+0.6*X1).*cos(6*pi*X1+repmat(J1,size(X,1),1)*pi/obj.D);
            X1      = repmat(X(:,1),1,length(J2));
            Y(:,J2) = X(:,J2)-(0.3*X1.^2.*cos(24*pi*X1+4*repmat(J2,size(X,1),1)*pi/obj.D)+0.6*X1).*sin(6*pi*X1+repmat(J2,size(X,1),1)*pi/obj.D);
            PopObj(:,1) = X(:,1)         + 2*mean(Y(:,J1).^2,2);
            PopObj(:,2) = 1-sqrt(X(:,1)) + 2*mean(Y(:,J2).^2,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^0.5;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
        	R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


% Copyright Andrew Binning 2013
% Please feel free to use and modify this code as you see if fit. If you
% use this code in any academic work, please cite
% Andrew Binning, 2013.
% "Underidentified SVAR models: A framework for combining short and long-run restrictions with sign-restrictions,"
% Working Paper 2013/14, Norges Bank.
function C = generateDraw(C,k)
%==========================================================================
% Generates a draw that is consistent with the shock variance/covariance
% matrix. Based on Juan F. Rubio-Ramirez & Daniel F. Waggoner & Tao Zha, 2010.
% "Structural Vector Autoregressions: Theory of Identification and
% Algorithms for Inference," Review of Economic Studies, Oxford University Press, vol. 77(2), pages 665-696.
%
% inputs:
% C = initial impact matrix, usually from the cholesky decomposition of the
% forecast error variance decomposition
% k = number of dependent variables
%
% outputs:
% C = new draw of the short run impact matrix
%==========================================================================

newmatrix = randn(k,k);

[Q,R] = qr(newmatrix);

for ii = 1:k
    if R(ii,ii)<0
        Q(:,ii) = -Q(:,ii);
    end
end

C = C*Q;


text:
function A2 = double_diagonal_blocks(A, row_range, col_range)
% function A = double_diagonal_blocks(A, row_range, col_range)
% Double diagonal blocks of a block matrix A whose row range and col range
% are `row_range` and `col_range`
% required: `numel(row_range) == numel(col_range)
% -----------------------------------------------
% Author: Tiep Vu, thv102@psu.edu, 5/11/2016
%         (http://www.personal.psu.edu/thv102/)
% -----------------------------------------------
	if nargin == 0 % test mode
		row_range = [0 3 5 8];
		col_range = [0 4 7 10];
		A = rand(row_range(end), col_range(end));
		Ain = A;
		subplot(1, 3, 1); imagesc(A); title('input'); colormap jet;
	end
	%% check requirements
	if numel(row_range) ~= numel(col_range)
		error('number of blocks in each dimension of the input matrix must be the same');
	end
	%% MAIN
	mask_diagonal_id = build_diagonal_mask(row_range, col_range);
	A(mask_diagonal_id) = 2*A(mask_diagonal_id);
	%% test mode
	if nargin == 0
		subplot(1, 3, 2); imagesc(A); title('output');
		subplot(1, 3, 3); imagesc(A - Ain); title('difference');
		A2 = [];
	end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function pdf = pdf_discrete_value ( x_num, x, y, v_num, v )

%*****************************************************************************80
%
%% PDF_DISCRETE_VALUE evaluates the PDF of a discrete histogram.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer V_NUM, the number of sample values.
%
%    Input, real V(V_NUM,1), the sample values.
%
%    Input, integer X_NUM, the number of values in the discrete histogram.
%
%    Input, real X(X_NUM,1), the histogram data values.
%
%    Input, real Y(X_NUM,1), the normalized histogram values.
%
%    Output, real PDF(V_NUM,1), the values of the discrete PDF.
%
  pdf = zeros ( v_num, 1 );
%
%  Seek bracket intervals for each sample point.
%
  b(1:v_num,1) = r8vec_bracket6 ( x_num, x, v_num, v );
%
%  Sample points outside the interval are ignored.
%
  i = find ( b ~= -1 );
%
%  Linearly interpolate values within the interval.
%
  l = b(i);
  r = l + 1;
%
% pdf(i) = (  ( x(r) - v(i)        ) .* y(l)   ...
%          +  (        v(i) - x(l) ) .* y(r) ) ...
%          ./ ( x(r)        - x(l) );

  for i = 1 : v_num
    if ( b(i) ~= -1 )
      l = b(i);
      r = l + 1;
      pdf(i) = (  ( x(r) - v(i)        ) .* y(l)   ...
               +  (        v(i) - x(l) ) .* y(r) ) ...
               ./ ( x(r)        - x(l) );
    end
  end

  return
end


text:
clear
        % Memorarea timpului pornirii programului
t0=clock;
		% Prealocarea vectorului y
y=zeros(50001,1);
		% Initializarea unui contor de ciclu
n=0;
		% Executarea ciclului for
for t=0:0.001:50
	n=n+1;
	y(n)=sin(t);
end
		% Calcularea timpului parcurs de la
		% pornirea programului
durata=etime(clock,t0)

author: Sable

repo: mcbench-benchmarks


text:
% This is a test/example script which demonstrates the usage of the
% ProbabilisticDataAssociationX class.
% =========================================================================>

%% Load the ground truth data
load('multiple-robot-tracking.mat');

%% Plot settings
ShowPlots = 1;              % Set to 0 to hide plots
numTrueTracks = 3;
GroundTruthTracks = GroundTruthTracks(1:numTrueTracks);
for i=1:length(GroundTruthStateSequence)
    GroundTruthStateSequence{i} = GroundTruthStateSequence{i}(1:numTrueTracks);
end
%% Model parameter shortcuts
lambdaV = 100; % Expected number of clutter measurements over entire surveillance region
V = 10^2;     % Volume of surveillance region (10x10 2D-grid)
V_bounds = [0 10 0 10]; % [x_min x_max y_min y_max]
P_D = 0.9;    % Probability of detection
timestep_duration = duration(0,0,1);

%% Models
transition_model = ConstantVelocityX('VelocityErrVariance', 0.01^2,...
                                     'NumDims', 2,...
                                     'TimestepDuration', timestep_duration);
% measurement_model = LinearGaussianX('NumMeasDims', 2,...
%                                     'NumStateDims', 4,...
%                                     'MeasurementErrVariance', 0.2^2,...
%                                     'Mapping', [1 3]);
measurement_model = RangeBearing2CartesianX('NumStateDims',4,'MeasurementErrVariance',...
                                            [(pi/180)^2,0.1^2],'Mapping',[1 3]);
% clutter_model = PoissonRateUniformPositionX('ClutterRate',lambdaV,...
%                                             'Limits',[V_bounds(1:2);...
%                                                       V_bounds(3:4)]);
clutter_model = PoissonRateUniformPositionX('ClutterRate',lambdaV,...
                                            'Limits',[0, pi/2;...
                                                      0, 12]);
detection_model = ConstantDetectionProbabilityX('DetectionProbability',P_D);

% Compile the State-Space model
model = StateSpaceModelX(transition_model,measurement_model,'Clutter',clutter_model, 'Detection', detection_model);

%% Generate DataList
meas_simulator = MultiTargetMeasurementSimulatorX('Model',model);

% DataList = meas_simulator.simulate(GroundTruthStateSequence);
N = numel(DataList);

%% Base Filter
obs_covar= measurement_model.covar();
dist = GaussianDistributionX(zeros(4,1), transition_model.covar() + blkdiag(obs_covar(1,1), 0, obs_covar(2,2),0));
PriorState = ParticleStateX(dist,10000);
% PriorState = GaussianStateX(zeros(4,1), transition_model.covar() + blkdiag(obs_covar(1,1), 0, obs_covar(2,2),0));
base_filter = ParticleFilterX('Model', model, 'StatePrior', PriorState);

%% Data Associator
config.ClutterModel = clutter_model;
config.Clusterer = NaiveClustererX();
config.Gater = EllipsoidalGaterX(2,'GateLevel',10)';
config.DetectionModel = detection_model;
assocFilter = JointProbabilisticDataAssocX(config);

%% Metric Generator
ospa = OSPAX('CutOffThreshold',1,'Order',1);

%% Initiate TrackList
NumTracks = numTrueTracks;
TrackList = cell(1,NumTracks);
for i=1:NumTracks
    xPrior = [GroundTruth{1}(1,i); 0; GroundTruth{1}(2,i); 0];
    PPrior = 10*transition_model.covar();
    StatePrior = GaussianStateX(xPrior,PPrior);
    TrackList{i} = TrackX(StatePrior);
    TrackList{i}.addprop('Filter');
    TrackList{i}.Filter = copy(base_filter);
    TrackList{i}.Filter.initialise('Model',model,'StatePrior',StatePrior);
end
assocFilter.TrackList = TrackList;

%% START OF SIMULATION
%  ===================>

% Create figure windows
% Create figure windows
if(ShowPlots)
    img = imread('maze.png');

    % set the range of the axes
    % The image will be stretched to this.
    min_x = 0;
    max_x = 10;
    min_y = 0;
    max_y = 10;

    % make data to plot - just a line.
    x = min_x:max_x;
    y = (6/8)*x;

    figure('units','normalized','outerposition',[0 0 .5 1])
    subplot(2,1,1);
end

ospa_vals= zeros(N,3);
for k=2:N
    fprintf('Iteration = %d/%d\n================>\n',k,N);

    %% Extract DataList at time k
    MeasurementList = DataList(k);
    timestamp_km1 = DataList(k-1).Timestamp;
    timestamp_k = MeasurementList.Timestamp;
    dt = timestamp_k - timestamp_km1;
    transition_model.TimestepDuration = dt;
    fprintf('Timestamp = %s\n================>\n',timestamp_k);

    %% Process JPDAF
    assocFilter.MeasurementList = MeasurementList;
    assocFilter.TrackList = TrackList;
    assocFilter.predictTracks();
    assocFilter.associate();
    assocFilter.updateTracks();

    %% Update target trajectories
    for t = 1: numel(assocFilter.TrackList)
        assocFilter.TrackList{t}.Trajectory(end+1) = assocFilter.TrackList{t}.Filter.StatePosterior;
    end

    %% Evaluate performance metric
    [ospa_vals(k,1), ospa_vals(k,2), ospa_vals(k,3)]= ospa.evaluate(GroundTruthStateSequence{k},assocFilter.TrackList);

     %% Plot update step results
    if(ShowPlots)

        subplot(3,1,1:2);
        ax(1) = gca;
        cla(ax(1));
        imagesc(ax(1), V_bounds(1:2), V_bounds(3:4), flipud(img));
        hold on;
        if(exist('data_plot','var'))
            delete(data_plot);
        end
        data_inv = measurement_model.finv(MeasurementList.Vectors);
        data_plot = plot(ax(1), data_inv(1,:), data_inv(3,:),'k*','MarkerSize', 10);

        % Plot tracks
        for j=1:numel(assocFilter.TrackList)
            means = [assocFilter.TrackList{j}.Trajectory.Mean];
            h2 = plot(ax(1), means(1,:),means(3,:),'-','LineWidth',1);
            h2 = plotgaussellipse(assocFilter.TrackList{j}.Filter.StatePosterior.Mean([1 3]),...
                                  assocFilter.TrackList{j}.Filter.StatePosterior.Covar([1 3],[1 3]),...
                                  'Color','r',...
                                  'Axis',ax(1));
        end
        % set the y-axis back to normal.
        str = sprintf('Target positions (Update)');
        set(ax(1), 'ydir','normal');
        title(str)
        xlabel('X position (m)')
        ylabel('Y position (m)')
        axis(ax(1), V_bounds)

        % Plot metric
        subplot(3,1,3);
        cla;
        plot(1:k,ospa_vals(1:k,1));
        title('OSPA vs Time');
        pause(0.01);

    end
end

figure
subplot(2,2,[1 2]), plot(1:k,ospa_vals(1:k,1));
subplot(2,2,3), plot(1:k,ospa_vals(1:k,2));
subplot(2,2,4), plot(1:k,ospa_vals(1:k,3));

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function varargout = DIAOpenSEGY(varargin)
% DIAOPENSEGY Application M-file for DIAOpenSEGY.fig
%    FIG = DIAOPENSEGY launch DIAOpenSEGY GUI.
%    DIAOPENSEGY('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.0 12-Jun-2002 10:25:58

if ((nargin == 0)|isstruct(varargin{1}))  % LAUNCH GUI

	fig = openfig(mfilename,'reuse');

	% Use system color scheme for figure:
	set(fig,'Color',get(0,'defaultUicontrolBackgroundColor'));

	% Generate a structure of handles to pass to callbacks, and store it.
	handles = guihandles(fig);
	guidata(fig, handles);

    data = guihandles(fig); % initialize it to contain handles
    data.JUMP.jump=5;
    data.JUMP.enable=0;
    data.TRACEHEADER.name='offset';
    data.TRACEHEADER.min=1;
    data.TRACEHEADER.max=1e+9;
    data.TRACEHEADER.enable=0;
    data.TIMERANGE.min=1;
    data.TIMERANGE.max=5;
    data.TIMERANGE.enable=0;

    if nargin==0;
      data.pathname='';
      data.filename='';
      guidata(fig,data);
    else
      try;data.JUMP=varargin{1}.JUMP;end
      try;data.TRACEHEADER=varargin{1}.TRACEHEADER;end
      try;data.TIMERANGE=varargin{1}.TIMERANGE;end
      if isfield(varargin{1},'pathname');data.pathname=varargin{1}.pathname;end;
      if isfield(varargin{1},'filename');data.filename=varargin{1}.filename;end
    end
    guidata(fig,data);

    set(fig,'HandleVisibility','On')
    DIAOpenSEGY('actionUpdateGUI',fig,handles)
    set(fig,'HandleVisibility','CallBack')

    % Wait for callbacks to run and window to be dismissed:
	uiwait(fig);

    if nargout > 0
       data=guidata(fig);
       SEGYFILE.JUMP=data.JUMP;
       SEGYFILE.TRACEHEADER=data.TRACEHEADER;
       SEGYFILE.TIMERANGE=data.TIMERANGE;
       if isfield(data,'SegyHeader')
         SEGYFILE.SegyHeader=data.SegyHeader;
       end
       if isfield(data,'filename'),SEGYFILE.filename=data.filename;end
       if isfield(data,'pathname'),SEGYFILE.pathname=data.pathname;end
       varargout{1} = SEGYFILE;
    end

    close(handles.figure1);
    try;set(handles.figure1,'Visible','off');end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch
		disp(lasterr);
	end

end

%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and
%| sets objects' callback properties to call them through the FEVAL
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the
%| callback type separated by '_', e.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(h, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.

function actionUpdateGUI(h,handles)
data=guidata(h);

% RESET FILEMAME,PATHNAME if FILE DOES NOT EXIST
if isfield(data,'filename')
  if exist([data.pathname,data.filename])~=2
      data=rmfield(data,'filename');
      data=rmfield(data,'pathname');
  end
end

% CHECK THAT A FILENAME I SET, and disbale if not.
if isfield(data,'filename'),
    cbEnable='on';
    set(handles.eFilename,'String',[data.pathname,data.filename]);
else
    cbEnable='off';
    set(handles.eFilename,'String','');
end
set(findobj('style','checkbox'),'enable',cbEnable)

% JUMP
set(handles.cbJump,'Value',data.JUMP.enable);
data.Jump.value=get(handles.cbJump,'Value');
if ((get(handles.cbJump,'Value')==1)&(isfield(data,'filename'))),
    enableJump='on';
    data.JUMP.enable=1;
else,
    enableJump='off';
    data.JUMP.enable=0;
end
set(handles.eJump,'Enable',enableJump)
set(handles.eJump,'String',data.JUMP.jump)

% TraceHeader
set(handles.cbTraceHeader,'Value',data.TRACEHEADER.enable);
data.TraceHeader.value=get(handles.cbTraceHeader,'Value');
if ((get(handles.cbTraceHeader,'Value')==1)&(isfield(data,'filename'))),
    enableTraceHeader='on';
    data.TRACEHEADER.enable=1;
else,
    enableTraceHeader='off';
    data.TRACEHEADER.enable=0;
end
hnames=get(handles.popTraceHeader,'string');
selected=get(handles.popTraceHeader,'value');
for i=1:length(hnames)
    if strcmp(hnames{i},data.TRACEHEADER.name),
        set(handles.popTraceHeader,'Value',i);
    end
end
set(handles.popTraceHeader,'Enable',enableTraceHeader)
set(handles.eTraceHeaderMin,'Enable',enableTraceHeader)
set(handles.eTraceHeaderMin,'String',data.TRACEHEADER.min)
set(handles.eTraceHeaderMax,'Enable',enableTraceHeader)
set(handles.eTraceHeaderMax,'String',data.TRACEHEADER.max)

% TimeRange
set(handles.cbTimeRange,'Value',data.TIMERANGE.enable);
data.TimeRange.value=get(handles.cbTimeRange,'Value');
if ((get(handles.cbTimeRange,'Value')==1)&(isfield(data,'filename'))),
    enableTimeRange='on';
    data.TIMERANGE.enable=1;
else,
    data.TIMERANGE.enable=0;
    enableTimeRange='off';
end
set(handles.eTimeRangeMin,'Enable',enableTimeRange)
set(handles.eTimeRangeMin,'String',data.TIMERANGE.min)
set(handles.eTimeRangeMax,'Enable',enableTimeRange)
set(handles.eTimeRangeMax,'String',data.TIMERANGE.max)

guidata(h,data)

% --------------------------------------------------------------------
function varargout = eFilename_Callback(h, eventdata, handles, varargin)

% --------------------------------------------------------------------
function varargout = pbSelectFile_Callback(h, eventdata, handles, varargin)
data=guidata(h);
[filename, pathname] = uigetfile('*.su;*.SU;*.sgy;*.SGY;*.seg*;*.SEGY', 'Pick a SEGY or SU file');

    if isequal(filename,0)|isequal(pathname,0)
       disp('File not found')
       return
    else
       disp(['File ', pathname, filename, ' found'])
       set(handles.eFilename,'String',[pathname,filename])
       data.pathname=pathname;
       data.filename=filename;
    end
guidata(h,data)
actionUpdateGUI(h,handles)

% --------------------------------------------------------------------
function varargout = cbJump_Callback(h, eventdata, handles, varargin)
data=guidata(h);
data.JUMP.enable=get(h,'Value');
guidata(h,data);
actionUpdateGUI(h,handles)

% --------------------------------------------------------------------
function varargout = eJump_Callback(h, eventdata, handles, varargin)
data=guidata(h);
if ~isempty(str2num(get(h,'String')));data.JUMP.jump=str2num(get(h,'String'));end
guidata(h,data);
actionUpdateGUI(h,handles)

% --------------------------------------------------------------------
function varargout = cbTraceHeader_Callback(h, eventdata, handles, varargin)
data=guidata(h);
data.TRACEHEADER.enable=get(h,'Value');
guidata(h,data);
actionUpdateGUI(h,handles)

% --------------------------------------------------------------------
function varargout = popTraceHeader_Callback(h, eventdata, handles, varargin)
data=guidata(h);
hnames=get(h,'string');
selected=get(h,'value');
data.TRACEHEADER.name=hnames{selected};
guidata(h,data);
actionUpdateGUI(h,handles)

% --------------------------------------------------------------------
function varargout = eTraceHeaderMin_Callback(h, eventdata, handles, varargin)
data=guidata(h);
if ~isempty(str2num(get(h,'String')));data.TRACEHEADER.min=str2num(get(h,'String'));end
if data.TRACEHEADER.min>data.TRACEHEADER.max
    data.TRACEHEADER.min=data.TRACEHEADER.max;
end
guidata(h,data);
actionUpdateGUI(h,handles)

% --------------------------------------------------------------------
function varargout = eTraceHeaderMax_Callback(h, eventdata, handles, varargin)
data=guidata(h);
if ~isempty(str2num(get(h,'String')));data.TRACEHEADER.max=str2num(get(h,'String'));end
if data.TRACEHEADER.max<data.TRACEHEADER.min
    data.TRACEHEADER.max=data.TRACEHEADER.min;
end
guidata(h,data);
actionUpdateGUI(h,handles)

% --------------------------------------------------------------------
function varargout = cbTimeRange_Callback(h, eventdata, handles, varargin)
data=guidata(h);
data.TIMERANGE.enable=get(h,'Value');
guidata(h,data);
actionUpdateGUI(h,handles)

% --------------------------------------------------------------------
function varargout = eTimeRangeMin_Callback(h, eventdata, handles, varargin)
data=guidata(h);
if ~isempty(str2num(get(h,'String')));data.TIMERANGE.min=str2num(get(h,'String'));end
if data.TIMERANGE.min>data.TIMERANGE.max;
    data.TIMERANGE.min=data.TIMERANGE.max;
end
guidata(h,data);
actionUpdateGUI(h,handles)

% --------------------------------------------------------------------
function varargout = eTimeRangeMax_Callback(h, eventdata, handles, varargin)
data=guidata(h);
if ~isempty(str2num(get(h,'String')));data.TIMERANGE.max=str2num(get(h,'String'));end
if data.TIMERANGE.max<data.TIMERANGE.min
    data.TIMERANGE.max=data.TIMERANGE.min;
end
guidata(h,data);
actionUpdateGUI(h,handles)

% --------------------------------------------------------------------
function varargout = pbReadSegy_Callback(h, eventdata, handles, varargin)
uiresume

% --------------------------------------------------------------------
function varargout = pbCancel_Callback(h, eventdata, handles, varargin)
data=guidata(h);
if isfield(data,'filename'), data=rmfield(data,'filename'); end
if isfield(data,'pathname'), data=rmfield(data,'pathname'); end
guidata(h,data);
uiresume

% --------------------------------------------------------------------
function varargout = pbGUIEditSegyHeader_Callback(h, eventdata, handles, varargin)
data=guidata(h);

% GET SEGY HEDER FROM FILE, UNLESS WE ALLREADY GOT IT
if isfield(data,'SegyHeader')==0
  data.SegyHeader=GetSegyHeader([data.pathname,data.filename]);
end

% EDIT THE SEGY HEADER
data.SegyHeader=GUIEditSegyHeader(data.SegyHeader);

guidata(h,data);

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
%{
    Copyright (C) Cambridge Electronic Design Limited 2014
    Author: James Thompson
    Web: www.ced.co.uk email: james@ced.co.uk, softhelp@ced.co.uk

    This file is part of CEDS64ML, a MATLAB interface to the SON64 library.

    CEDS64ML 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.

    CEDS64ML 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 CEDS64ML.  If not, see <http://www.gnu.org/licenses/>.
%}

function [ iMaxChans ] = CEDS64MaxChan( fhand )
%CEDS64MAXCHAN Returns the maximum number of channels that can stored in this
%file
%   [ iMaxChans ] = CEDS64MaxChan( fhand )
%   Inputs
%   fhand - Integer file handle
%   Outputs
%   iMaxChans - The maximum number of channels, otherwise a negative error
%   code.

if (nargin ~= 1)
    iMaxChans = -22;
    return;
end

iMaxChans = calllib('ceds64int', 'S64MaxChans', fhand);
end

author: brainstorm-tools


text:
function [out] = recharge_1(p1,S,Smax,flux)
%recharge_1

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Recharge as scaled fraction of incoming flux
% Constraints:  -
% @(Inputs):    p1   - fraction of flux that is recharge [-]
%               S    - current storage [mm]
%               Smax - maximum contributing storage [mm]
%               flux - incoming flux [mm/d]

out = p1*S/Smax*flux;

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function Population = DecompositionSelection(Global,Population,associate,Cosinemax)
% The decomposition-based method environmental selection

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

% This function is written by Shufen Qin
% E-mail: shufen.qin@stu.tyust.edu.cn

    np = length(Population);
    %% Normalization
    Obj = Population.objs;
    Obj = (Obj-repmat(min(Obj),np,1))./(repmat(max(Obj),np,1)-repmat(min(Obj),np,1));

    %% Select one solution for each reference vector
    list = unique(associate)';
    Next = zeros(length(list),1);
    t = 1;
    for i = list
        current = find(associate == i);
        dist = pdist2(Obj(current,:),zeros(1,Global.M),'Euclidean');
        Fan = Cosinemax(current)./dist;
        [~,best] = max(Fan);
        Next(t)  = current(best);
        t = t +1;
    end
    % Population for next generation
    Population = Population(Next);
end


text:
function center = initComponents_2p(obj, Y, K, debug_on, save_avi)
%% a greedy method for detecting ROIs and initializing CNMF. in each iteration,
% it searches the one with large (peak-median)/noise level and large local
% correlation
%% Input:
%   Y:  d X T matrx, imaging data
%   K:  scalar, maximum number of neurons to be detected.
%   debug_on: options for showing procedure of detecting neurons
%% Output:
%       Ain:  d X K' matrix, estimated spatial component
%       Cin:  K'X T matrix, estimated temporal component
%       bin:  d X nb matrix/vector, spatial components of the background
%       Cin:  nb X T matrix/vector, temporal components of the background
%       center: K' X 2, coordinate of each neuron's center
%       res:  d X T, residual after initializing Ain, Cin, bin, fin

%% Author: Pengcheng Zhou, Carnegie Mellon University.
% the method is an modification of greedyROI method used in Neuron paper of Eftychios
% Pnevmatikakis et.al. https://github.com/epnev/ca_source_extraction/blob/master/utilities/greedyROI2d.m
%% In each iteration of initializing neurons, it searchs the one with maximum
% value of (max-median)/noise * Cn, which selects pixels with high SNR and
% local correlation.

%% parameters
if exist('K', 'var')
    K = 200;
end
Y = obj.reshape(Y, 1);
Y_std = get_noise_fft(Y);
options = obj.options;

if ~exist('debug_on', 'var'); debug_on = false; end
if ~exist('save_avi', 'var'); save_avi=false; end
d1 = options.d1;
d2 = options.d2;
gSig = options.gSig;
gSiz = options.gSiz;
if and(isempty(gSiz), isempty(gSig)); gSig = 3; gSiz = 10; end
if isempty(gSiz); gSiz=3*gSig; end
if isempty(gSig); gSig=gSiz/3; end
if isfield(options, 'min_corr')
    min_corr = options.min_corr;    % minimum local correaltion value to start one neuron
else
    min_corr = 0.3;
end
if isfield(options, 'deconv_flag')
    deconv_flag = options.deconv_flag;
    deconv_options_0= options.deconv_options;
else
    deconv_flag = false;
end
if isfield(options, 'min_pnr')
    min_pnr = options.min_pnr;
else
    min_pnr = 2;
end
nb = options.nb;        % number of the background
pSiz = 1;       % after selecting one pixel, take the mean of square box
%near the pixel as temporal activity. the box size is (2*pSiz+1)
psf = ones(gSig)/(gSig^2);

maxIter = 5;            % iterations for refining results
sz = 4;            %distance of neighbouring pixels for computing local correlation

if ~ismatrix(Y); Y = reshape(Y, d1*d2, []); end;
[~, T] = size(Y);       % number of frames
Ain = zeros(d1*d2, K);  % spatial components
Cin = zeros(K, T);      % temporal components
Cin_raw = zeros(K, T);      % temporal components
Sin = zeros(K, T);      % temporal components
center = zeros(K, 2);   % center of the initialized components
kernel_pars = cell(K,1);    % parameters for the convolution kernels of all neurons

%% compute correlation image and (max-median)/std ratio
ind_frame = round(linspace(1, T, min(T, 1000)));    % select few frames for speed issues
%tmp_noise = randn(d1*d2, length(ind_frame));
C1 = correlation_image(full(Y(:, ind_frame)), sz, d1, d2);
Cb =  zeros(size(C1)); %correlation_image(full(Y(:, ind_frame(1:3:end)))+tmp_noise(:, 1:3:end), [gSiz, gSiz+1], d1, d2);  %backgroung correlatin
Cn = C1-Cb; %here Cb is the background correlation. for 2photon imaging results. It might be useful when the background signal is large
Y_median = median(Y(:, ind_frame), 2);
Y = bsxfun(@minus, Y, Y_median);
% Y_std = sqrt(mean(Y.*Y, 2));

%% find local maximum
k = 0;      %number of found components
min_pixel = floor(gSig^2/2);  % minimum number of peaks to be a neuron
peak_ratio = full(max(Y, [], 2))./Y_std; %(max-median)/std
peak_ratio(isinf(peak_ratio)) = 0;  % avoid constant values
peak_ratio(peak_ratio<min_pnr) = 0;

% save_avi = false;   %save procedures for demo
if debug_on
    figure('position', [100, 100, 800, 650]); %#ok<*UNRCH>
      figure('position', [100, 100, 1200, 800], 'color', [1,1,1]*0.9); %#ok<*UNRCH>
    set(gcf, 'defaultAxesFontSize', 20);
    ax_cn = axes('position', [0.04, 0.5, 0.3, 0.4]);
    ax_cn_varying = axes('position', [0.36, 0.5, 0.3, 0.4]);
    ax_cn_box = axes('position', [0.68, 0.54, 0.24, 0.32]);
    ax_raw = axes('position', [0.05, 0.25, 0.92, 0.2]);
    ax_trace = axes('position', [0.05, 0.01, 0.92, 0.2]);
    axes(ax_cn);
    imagesc(Cn);  hold on;
    axis equal off tight;
    title('correlation image');
    axes(ax_cn_varying);
    imagesc(Cn, [0,1]);
    axis equal off tight; hold on;
    if save_avi
        avi_file = VideoWriter('greedyROI_example.avi');
        avi_file.open();
    end
end

max_thresh = min_pnr * (min_corr);
while k<K
    %% find the pixel with the maximum ratio
    [max_v, ind_p] = max(peak_ratio.*(Cn(:)));
    peak_ratio(ind_p) = 0;  % no longer visit this pixel any more
    if max_v<max_thresh; break; end
    if Cn(ind_p)<min_corr; continue; end % ignore this local maximum due to small local correlation
    if  max_v/(min_corr)< min_pnr;     continue;    end

    [r, c] = ind2sub([d1,d2], ind_p);

    % select its neighbours for computing correlation
    rsub = max(1, -gSiz+r):min(d1, gSiz+r);
    csub = max(1, -gSiz+c):min(d2, gSiz+c);
    [cind, rind] = meshgrid(csub, rsub);
    nr = length(rsub);  %size of the neighboring matrix
    nc = length(csub);
    ind_nhood = sub2ind([d1, d2], rind(:), cind(:));
    Y_box = Y(ind_nhood, :);

    % draw a small area near the peak and extract the mean activities
    r0 = rsub(1); c0 = csub(1);
    rsub = (max(1, -pSiz+r):min(d1, pSiz+r)) - r0+1;
    csub = (max(1, -pSiz+c):min(d2, pSiz+c)) -c0+1;
    [cind, rind] = meshgrid(csub, rsub);
    ind_peak = sub2ind([nr, nc], rind(:), cind(:));
    y0 = mean(Y_box(ind_peak, :), 1);
%     y0(y0<0) = 0;

    % compute the correlation between the peak and its neighbours
    temp = reshape(corr(y0', Y_box'), nr, nc);
    active_pixel = full(temp>min_corr/2);
    l = bwlabel(active_pixel, 8);   % remove disconnected components
    active_pixel(l~=mode(l(ind_peak))) = false;
    tmp_v = sum(active_pixel(:));    %number of pixels with above-threshold correlation
    if debug_on
        axes(ax_cn_varying); cla;
        imagesc(reshape(Cn(:), d1, d2), [0, Cn(ind_p)]);
        title(sprintf('neuron %d', k+1));
        axis equal off tight; hold on;
        plot(c,r, 'om');
        axes(ax_cn_box);
        imagesc(temp, [min_corr, 1]);
        axis equal off tight;
        title('local corr.');
        axes(ax_raw); cla; hold on;
        plot(y0/max(y0)); title('activity in the center');
        axis off tight;
        if ~save_avi; pause; end
    end
    if tmp_v<min_pixel;         continue;  end % neuron is too small

    %% save neuron
    %   nonzero area
    data = Y_box(active_pixel(:), :);
    ind_active = ind_nhood(active_pixel(:));  %indices of active pixels within the whole frame
    peak_ratio(ind_nhood(ind_peak)) = 0;    % the small area near the peak is not able to initialize neuron anymore

    % do a rank-1 matrix factorization in this small area
    [ai, ci_raw] = finetune2d(data, y0, maxIter);
    if norm(ai)==0;        continue;   end
    k = k+1;
    if deconv_flag
        % deconv the temporal trace
        [ci, si, deconv_options] = deconvolveCa(ci_raw, deconv_options_0);  % sn is 1 because i normalized c_raw already
        % save this initialization
        Ain(ind_active, k) = ai;
        Cin(k, :) = ci;
        Sin(k, :) = si;
        Cin_raw(k, :) = ci_raw;
        %                 kernel_pars(k, :) = kernel.pars;
        kernel_pars{k} = reshape(deconv_options.pars, 1, []);
    else
        ci = ci_raw;
        Ain(ind_active, k) = ai;
        Cin(k, :) = ci_raw;
        Cin_raw(k, :) = ci_raw;
    end
    ci = reshape(ci, 1, []);
    Y(ind_active, :) = data-ai*ci;
    center(k, :) = [r, c];

    if debug_on
        axes(ax_cn);
        plot(c, r, '.r');
        axes(ax_cn_varying);
        plot(c,r, 'or');
        axes(ax_cn_box);
        temp = zeros(nr, nc); temp(active_pixel) = ai;
        imagesc(temp);
        axis equal off tight;
        title('spatial component');
        axes(ax_trace);  cla; hold on;
        if deconv_flag
                    plot(ci_raw, 'linewidth', 2); title('temporal component'); axis tight;
            plot(ci, 'r', 'linewidth', 1);  axis tight;
            legend('raw trace', 'denoised trace');
        else
            plot(ci_raw, 'r', 'linewidth', 1); title('temporal component'); axis tight;
        end
        title('temporal component');
        if exist('avi_file', 'var')
            temp = getframe(gcf);
            temp.cdata = imresize(temp.cdata, [800, 640]);
            avi_file.writeVideo(temp);
        elseif ~save_avi
            temp = input('type s to stop the debug mode:  ', 's');
            if strcmpi(temp, 's')
                save_avi = true;
            end
        else
            drawnow();
        end
    end

    if mod(k, 10)==0
        fprintf('%d/%d neurons have been detected\n', k, K);
    end

    if k==K;   break; end

    %% udpate peak_ratio and correlation image
    tmp_old = peak_ratio(ind_active);
    tmp_new = max(Y(ind_active, :), [], 2)./Y_std(ind_active);
    temp = zeros(nr, nc);
    temp(active_pixel) = max(0, tmp_old-tmp_new); % after each iteration, the peak ratio can not be increased
    peak_ratio(ind_nhood) = max(0, peak_ratio(ind_nhood) - reshape(imfilter(temp, psf), [], 1)); % update peak_ratio, results are smoothed
    Cn(ind_nhood) = correlation_image(full(Y(ind_nhood, ind_frame)), sz, nr, nc)-reshape(Cb(ind_nhood), nr, nc);  % update local correlation
end

center = center(1:k, :);
obj.A = sparse(Ain(:, 1:k));
obj.C = Cin(1:k, :);
obj.C_raw = Cin_raw(1:k, :);
obj.S = Sin(1:k, :);
obj.P.kernel_pars = kernel_pars(1:k);

if exist('avi_file', 'var'); avi_file.close(); end
res = bsxfun(@plus, Y, Y_median);

% %% initialize background
% tsub = max(1, round(T/1000));
% [bin, f] = nnmf(max(res(:, 1:tsub:T), 0), nb);
% fin = imresize(f, [nb, T]);
% fin = HALS_temporal(max(res, 0), bin, fin, maxIter);
% obj.b = HALS_spatial(max(res, 0), bin, fin, [], maxIter);
% obj.f = fin;
end

function [ai, ci] = finetune2d(data, ci, nIter)
%do matrix factorization given the model data = ai*ci, where ai>=0
%
%Input:
%   data:   d x T matrix, small patch containing one neuron
%   ci:     initial value for trace
%   nIter  number of coordinate descent steps
%
%Output:
%   ai  M x N matrix, result of the fine-tuned neuron shape
%   ci  1 x T matrix, result of the neuron
%% copied from greedyROI.m

if ~exist('nIter', 'var'), nIter = 1; end
data(data<0)= 0;
%do block coordinate descent
for iter = 1:nIter,
    %update basis
    ai = max(0, (data*ci')/(ci*ci'));
    norm_ai = norm(ai, 2);
    if norm_ai==0; break;     end
    ai = ai/norm_ai;
    ci =  (ai'*data);
    %     ci(ci<0) = 0;
end
% [b, sn] = estimate_baseline_noise(ci);
% ci = ci - b;
% ai = ai*sn;
% ci = ci/sn;
end

author: zhoupc


text:
classdef SideSlipAngleTermCondition < AbstractEventTerminationCondition
    %SideSlipAngleTermCondition Summary of this class goes here
    %   Detailed explanation goes here

    properties
        steeringModel(1,1) AbstractSteeringModel = RollPitchYawPolySteeringModel.getDefaultSteeringModel();
        targetSlipAngle(1,1) double = 0;
        bodyInfo KSPTOT_BodyInfo
    end

    methods
        function obj = SideSlipAngleTermCondition(targetSlipAngle)
            obj.targetSlipAngle = targetSlipAngle;
        end

        function evtTermCondFcnHndl = getEventTermCondFuncHandle(obj)
            evtTermCondFcnHndl = @(t,y) obj.eventTermCond(t,y, obj.targetSlipAngle, obj.steeringModel, obj.bodyInfo);
        end

        function initTermCondition(obj, initialStateLogEntry)
            obj.steeringModel = initialStateLogEntry.steeringModel;
            obj.bodyInfo = initialStateLogEntry.centralBody;
        end

        function name = getName(obj)
            name = sprintf('Side Slip Angle (%.3f deg)', rad2deg(obj.targetSlipAngle));
        end

        function tf = shouldBeReinitOnRestart(obj)
            tf = false;
        end

        function params = getTermCondUiStruct(obj)
            params = struct();

            params.paramName = 'Target Side Slip Angle';
            params.paramUnit = 'deg';
            params.useParam = 'on';
            params.useStages = 'off';
            params.useTanks = 'off';
            params.useEngines = 'off';
            params.useStopwatches = 'off';

            params.value = rad2deg(obj.targetSlipAngle);
            params.refStage = LaunchVehicleStage.empty(1,0);
            params.refTank = LaunchVehicleEngine.empty(1,0);
            params.refEngine = LaunchVehicleEngine.empty(1,0);
            params.refStopwatch = LaunchVehicleStopwatch.empty(1,0);
        end

        function optVar = getNewOptVar(obj)
            optVar = SideSlipAngleTermCondOptimVar(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end
    end

    methods(Static)
        function termCond = getTermCondForParams(paramValue, stage, tank, engine, stopwatch)
            termCond = SideSlipAngleTermCondition((paramValue));
        end
    end

    methods(Static, Access=private)
        function [value,isterminal,direction] = eventTermCond(t,y, targetSlipAngle, steeringModel, bodyInfo)
            ut = t;
            rVect = y(1:3);
            vVect = y(4:6);

            dcm = steeringModel.getBody2InertialDcmAtTime(ut, rVect, vVect, bodyInfo);
            [~,~,angOfSideslip] = computeAeroAnglesFromBodyAxes(ut, rVect, vVect, bodyInfo, dcm(:,1), dcm(:,2), dcm(:,3));

            angOfSideslip = AngleZero2Pi(angOfSideslip);

            value = angOfSideslip - targetSlipAngle;
            isterminal = 1;
            direction = 0;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function R = sw_rcu(a, n_vals, eps)
% Compute the achievable sum rate (at the symmetrical rate point) from the SW RCU bound
% for the joint source distribution: [a 1/3*(1-a); 1/3*(1-a) 1/3*(1-a)].

% Input:
% a: parameter to define the joint source distribution, should be in range [0.25, 1)
% n_vals: a vector containing the blocklengths
% eps: target error probability
% Output:
% R: the achievable rates for the blocklengths specified in n_vals.

    p = [a 1/3*(1-a); 1/3*(1-a) 1/3*(1-a)];
    H = sum(sum(p.*log(1./p)));

    % Compute log binomial coefficients
    log_b = binomial_coeff(max(n_vals));

    R = zeros(length(n_vals),1);

    for i = 1:length(n_vals)
        n = n_vals(i);
        m = 0:n;
        % log_Pr stores all the log joint probabilities in ascending order
        log_Pr = m*log(a)+(n-m)*log(1/3*(1-a));

        % Calculate and store the coefficients that will be used iteratively later
        A_m = {};
        A_kj = {};
        for m = 0:n
            A_tempt = log_b{n+1}+(n:-1:0)*log(3);
            A_m{m+1} = my_logsumexp(A_tempt(m+1:n+1));
            for k = 0:n-m
                A_kj{m+1}(k+1) = my_logsumexp(log_b{m+k+1}(m+1:m+k+1))+(n-m-k)*log(2);
            end
        end

        % [x, y]: initial range for the bisection algorithm
        % For n small or eps small, need to adjust the right boundary of this initial interval.
        x = H;
        y = H + 0.5;
        %%%%%%%%%% Use this part to check if the range is valid for bisection algorithm %%%%%%%%%%
        c = x;
        log_err_m = zeros(1,n+1);
        for m = 0:n
            log_err_k = zeros(1,n-m+1);
            for k = 0:n-m
                log_err_j = zeros(1,n-m-k+1);
                for j = 0:n-m-k
                    log_err_j(j+1) = log_b{n-m-k+1}(j+1) + ...
                        min(0, my_logsumexp([A_kj{m+1}(j+1)-n*c/2 A_kj{m+1}(k+1)-n*c/2 A_m{m+1}-n*c]));
                end
                log_err_k(k+1) = log_b{n-m+1}(k+1) + my_logsumexp(log_err_j);
            end
            log_err_m(m+1) = log_Pr(m+1) + log_b{n+1}(m+1) + my_logsumexp(log_err_k);
        end
        err = exp(my_logsumexp(log_err_m));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % Start of bisection iteration
        iter = 0;
        while abs(x - y) > 0.00001 && iter <= 30
            c = (x+y)/2;
            %%%%%%%%%%%%%%%% body of iteration %%%%%%%%%%%%%%%%
            log_err_m = zeros(1,n+1);
            for m = 0:n
                log_err_k = zeros(1,n-m+1);
                for k = 0:n-m
                    log_err_j = zeros(1,n-m-k+1);
                    for j = 0:n-m-k
                        log_err_j(j+1) = log_b{n-m-k+1}(j+1) + ...
                            min(0, my_logsumexp([A_kj{m+1}(j+1)-n*c/2 A_kj{m+1}(k+1)-n*c/2 A_m{m+1}-n*c]));
                    end
                    log_err_k(k+1) = log_b{n-m+1}(k+1) + my_logsumexp(log_err_j);
                end
                log_err_m(m+1) = log_Pr(m+1) + log_b{n+1}(m+1) + my_logsumexp(log_err_k);
            end
            err = exp(my_logsumexp(log_err_m));
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if err - eps < 0
                y = c;
            else
                x = c;
            end
            iter = iter + 1;
        end
        R(i) = c;
    end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function ndvi = NDVI( red,nir )
%NDVI Calculate Normalized Difference Vegetation Index (NDVI) using NIR and
%     Red bands.
%
% Syntax
%
%     ndvi = NDVI(red,nir)
%
% Description
%
%     This function calculates Normalized Difference Vegetation Index (NDVI)
%     using NIR and Red bands (as following equation). This range is between
%     -1 and 1.
%     NDVI=(NIR-Red)/(NIR+Red).
%
% Input arguments
%
%     red         Red band
%     nir         Near-infrared band
%
% Output arguments
%
%     ndvi        Normalized Difference Vegetation Index
%
%
% Author: Shi Qiu (shi.qiu@ttu.edu)
% Date: 19. October, 2017

    % calculate NDVI
    ndvi=(nir-red)./(nir+red);

    % fix unnormal pixels
    ndvi((nir+red)==0)=0.01;

end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
classdef IspPressureCurveElement < AbstractCurveElement
    %IspPressureCurveElement Summary of this class goes here
    %   Detailed explanation goes here

    properties(Constant)
        minIndepValue = 0;
        maxIndepValue = NaN;

        minDepValue = 1E-6;
        maxDepValue = NaN;

        mandatoryIndepValues = [LaunchVehicleEngine.vacPress, LaunchVehicleEngine.seaLvlPress];
    end

    methods
        function obj = IspPressureCurveElement(atmoPress, engineIsp)
            obj.indepVar = atmoPress;
            obj.depVar = engineIsp;
        end

        function indepVarName = getIndepVarName(obj)
            indepVarName = 'Atmospheric Pressure';
        end

        function indepVarUnit = getIndepVarUnit(obj)
            indepVarUnit = 'kPa';
        end

        function depVarName = getDepVarName(obj)
            depVarName = 'Isp';
        end

        function depVarUnit = getDepVarUnit(obj)
            depVarUnit = 'sec';
        end
    end
end


text:
function demo
% Demo of the M2DP descriptor described in the following paper:
%
% Li He, Xiaolong Wang and Hong Zhang, M2DP: A Novel 3D Point Cloud
% Descriptor and Its Application in Loop Closure Detection, IROS 2016.
%
% Li He, Dept. of Computing Science, University of Alberta
% lhe2@ualberta.ca

%% 0. Please run mex CountVote2D.c in Matlab
% mex CountVote2D.c

%% 1. Read a point cloud
% 000000.bin is from KITTI dataset of sequence 00.
% http://www.cvlibs.net/datasets/kitti/eval_odometry.php
[data, numData] = ReadBinData('000000.bin');

% display the input data
figure(10);grid on;hold on
pcshow(data);title('Input point cloud');

%% 2. M2DP
% desM2DP: descriptor of data; A: 2D signatures of data.
tstart = tic;
[desM2DP, A] = M2DP(data);
telapsed = toc(tstart);
disp(['Processing time of M2DP on ' num2str(numData) ' points: ' num2str(telapsed) ' seconds']);

function [data, numData] = ReadBinData(nameDataFile)

fid = fopen(nameDataFile,'r');
if fid==-1
    disp('Cannot find input file\n');
    data=[];
    numData=[];
    return;
end

% get the size of file
fseek(fid,0,'eof');
fsize = ftell(fid);
% number points = total bytes / 4 per float32 / 4 float32 per row
numData = floor(fsize/16);

% now, read data
fseek(fid,0,'bof');
data = fread(fid,fsize,'float32');
fclose(fid);

% reshape data
data = reshape(data,[4, numData]);
data = data';
% ignoring the last feature, laser intensity
data = data(:,1:3);

author: LiHeUA


text:
% Make list of used BNT parameters in file
%
% This function parses input file or script for entries like `<pattern>name = value`,
% and prints out the list of names. An example of such entry is: p.imageFormat = 'png'.
% It can be used to quickly check what BNT parameters are used in the script.
%
%  USAGE
%   helpers.extractParametersList(filename, pattern)
%   filename    - Path to the file to parse or script name. If script name is provided,
%                 then Matlab's function which is used to find out the corresponding
%                 file name.
%   pattern     - Begin of BNT parameters pattern. If omitted, then 'p.' is used.
%
function extractParametersList(filename, pattern)
    if nargin < 2
        pattern = 'p.';
    end
    existCode = exist(filename, 'file');
    if existCode == 0 || existCode == 2
        filename = which(filename);
    end
    parameters = containers.Map('KeyType', 'char', 'ValueType', 'char');

    fid = data.safefopen(filename, 'r');

    while ~feof(fid)
        str = fgetl(fid);
        patternPos = strfind(str, pattern);
        if ~isempty(patternPos)
            originalPos = patternPos;
            for i = 1:length(originalPos)
                patternPos = originalPos(i);

                wholeWord = regexp(str(patternPos-1:patternPos-1), '\w', 'once');
                if ~isempty(wholeWord)
                    continue;
                end
                commentsPos = strfind(str, '%');
                if ~isempty(commentsPos) && commentsPos(1) < patternPos
                    % this string is commented out, skip it
                    continue;
                end
                parameterName = regexp(str(patternPos+2:end), '\w*', 'match');
                parameterName = sprintf('%s%s', pattern, parameterName{1});

                parameters(parameterName) = parameterName;
            end
        end
    end

    cellfun(@(x) fprintf('%%  %s\n', x), values(parameters));
end

author: kavli-ntnu


text:
% Detecting a small subset of infeasible linear inequalities
% Section 5.8, Boyd & Vandenberghe "Convex Optimization"
% Written for CVX by Almir Mutapcic - 02/18/06
%
% We consider a set of linear inequalities A*x <= b which are
% infeasible. Here A is a matrix in R^(m-by-n) and b belongs
% to R^m. We apply a l1-norm heuristic to find a small subset
% of mutually infeasible inequalities from a larger set of
% infeasible inequalities. The heuristic finds a sparse solution
% to the alternative inequality system.
%
% Original system is A*x <= b and it alternative ineq. system is:
%
%   lambda >= 0,   A'*lambda == 0.   b'*lambda < 0
%
% where lambda in R^m. We apply the l1-norm heuristic:
%
%   minimize   sum( lambda )
%       s.t.   A'*lambda == 0
%              b'*lambda == -1
%              lambda >= 0
%
% Positive lambdas gives us a small subset of inequalities from
% the original set which are mutually inconsistent.

% problem dimensions (m inequalities in n-dimensional space)
m = 150;
n = 10;

% fix random number generator so we can repeat the experiment
seed = 0;
randn('state',seed);

% construct infeasible inequalities
A = randn(m,n);
b = randn(m,1);

fprintf(1, ['Starting with an infeasible set of %d inequalities ' ...
            'in %d variables.\n'],m,n);

% you can verify that the set is infeasible
% cvx_begin
%   variable x(n)
%   A*x <= b;
% cvx_end

% solve the l1-norm heuristic problem applied to the alternative system
cvx_begin
   variables lambda(m)
   minimize( sum( lambda ) )
   subject to
     A'*lambda == 0;
     b'*lambda == -1;
     lambda >= 0;
cvx_end

% report the smaller set of mutually inconsistent inequalities
infeas_set = find( abs(b.*lambda) > sqrt(eps)/n );
disp(' ');
fprintf(1,'Found a smaller set of %d mutually inconsistent inequalities.\n',...
        length(infeas_set));
disp(' ');
disp('A smaller set of mutually inconsistent inequalities are the ones');
disp('with row indices:'), infeas_set'

% check that this set is infeasible
% cvx_begin
%    variable x_infeas(n)
%    A(infeas_set,:)*x_infeas <= b(infeas_set);
% cvx_end

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
%COFICOSTFUNC Collaborative filtering cost function
%   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
%   num_features, lambda) returns the cost and gradient for the
%   collaborative filtering problem.
%

% Unfold the U and W matrices from params
X = reshape(params(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(params(num_movies*num_features+1:end), ...
                num_users, num_features);

% You need to return the following values correctly
J = 0;
X_grad = zeros(size(X));
Theta_grad = zeros(size(Theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost function and gradient for collaborative
%               filtering. Concretely, you should first implement the cost
%               function (without regularization) and make sure it is
%               matches our costs. After that, you should implement the
%               gradient and use the checkCostFunction routine to check
%               that the gradient is correct. Finally, you should implement
%               regularization.
%
% Notes: X - num_movies  x num_features matrix of movie features
%        Theta - num_users  x num_features matrix of user features
%        Y - num_movies x num_users matrix of user ratings of movies
%        R - num_movies x num_users matrix, where R(i, j) = 1 if the
%            i-th movie was rated by the j-th user
%
% You should set the following variables correctly:
%
%        X_grad - num_movies x num_features matrix, containing the
%                 partial derivatives w.r.t. to each element of X
%        Theta_grad - num_users x num_features matrix, containing the
%                     partial derivatives w.r.t. to each element of Theta
%

J = sum(sum(((X * Theta') .* R .- Y) .^ 2)) / 2;
X_grad = ((X * Theta') .* R .- Y) * Theta;
Theta_grad = ((X * Theta') .* R .- Y)' * X;

offset = (lambda / 2) * (sum(sum(Theta .^ 2)) + sum(sum(X .^ 2)));
J = J + offset;
X_grad = X_grad + lambda * X;
Theta_grad = Theta_grad + lambda * Theta;
% =============================================================

grad = [X_grad(:); Theta_grad(:)];

end

author: ecmadao

repo: Coding-Guide


text:
function setL(OCP,L)
    OCP.L = L;
end

author: deng-haoyang

repo: ParNMPC


text:
function y=Tsallis_entro(x,q)
  [M,N]=size(x);
       y=zeros(1,N);
          for l=1:N
          sum1=sum(x(:,l)-(x(:,l)).^q);
          sum2=sum1/(q-1);
          y(1,l)=sum2;
       end


text:
function rule_num = lyness_rule_num ( )

%*****************************************************************************80
%
%% LYNESS_RULE_NUM returns the number of Lyness quadrature rules.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    James Lyness, Dennis Jespersen,
%    Moderate Degree Symmetric Quadrature Rules for the Triangle,
%    Journal of the Institute of Mathematics and its Applications,
%    Volume 15, Number 1, February 1975, pages 19-32.
%
%  Parameters:
%
%    Output, integer RULE_NUM, the number of rules.
%
  rule_num = 21;

  return
end

author: johannesgerer


text:
classdef AssessmentClassificationResultsPlotter < handle

    properties (Access = public, Constant)
        FontSize = 18;
        LineWidth = 2;
        MissedEventSymbolSize = 16;
        RectangleYPosToDataRatios = [1.0, 1.05];
        LabelYPosToRectangleRatios = [0.8,0.95];
        RectangleCurvature = 0.1;
    end

    properties (Access = public)
        yRange = [];
        missedEvents;
    end

    properties (Access = private)
        plotAxes;
        classNames;
        isUpperYPos = false;
    end

    methods (Access = public)
        function obj = AssessmentClassificationResultsPlotter(plotAxes,classNames)
            obj.plotAxes = plotAxes;
            obj.classNames = classNames;
        end

        function plotClassificationResults(obj, detailedClassificationResults, signal)

            nSegments = length(detailedClassificationResults.segments);
            for i = 1 : nSegments
                segment = detailedClassificationResults.segments(i);
                if ~isempty(segment.eventIdx)
                    y = signal(segment.eventIdx);
                else
                    y = [];
                end

                truthClass = detailedClassificationResults.classificationResult.truthClasses(i);
                predictedClass = detailedClassificationResults.classificationResult.predictedClasses(i);

                if truthClass ~= predictedClass || truthClass ~= Labeling.kNullClass
                    obj.plotClassificationResult(segment,y,truthClass,predictedClass);
                end
            end

            if ~isempty(obj.missedEvents)
                obj.plotMissedEventAnnotations(signal);
            end
        end
    end

    methods (Access = private)
        function plotHandle = plotClassificationResult(obj, segment, segmentEventY, truthClass, predictedClass)
            colorStr = obj.getColorForPrediction(truthClass,predictedClass);

            rectangleHeight = (obj.yRange(2) - obj.yRange(1)) * obj.RectangleYPosToDataRatios(obj.isUpperYPos+1);
            yOffset = (rectangleHeight - (obj.yRange(2) - obj.yRange(1))) / 2;
            rectangleWidth = single(segment.endSample - segment.startSample);
            rectanglePosition = [single(segment.startSample), obj.yRange(1) - yOffset, single(rectangleWidth), rectangleHeight];
            rectangleHandle = rectangle('Position',rectanglePosition,'Curvature',[obj.RectangleCurvature obj.RectangleCurvature],'LineWidth',obj.LineWidth, 'EdgeColor',colorStr);

            %plot text
            textStr = obj.getTextForPrediction(truthClass,predictedClass);
            xPos = (double(segment.startSample) + double(segment.endSample)) / 2;
            yPos = double(obj.yRange(2)) * obj.LabelYPosToRectangleRatios(obj.isUpperYPos+1);
            textHandle = text(obj.plotAxes,xPos,yPos, textStr,...
                'FontSize',obj.FontSize,'HorizontalAlignment','center','Color',colorStr);
            set(textHandle, 'Clipping', 'on');
            obj.isUpperYPos = ~obj.isUpperYPos;

            %plot symbol
            symbolHandle = [];
            if ~isempty(segmentEventY)
                x = segment.eventIdx;
                symbolHandle = plot(obj.plotAxes,x,segmentEventY,'*','Color',colorStr);
            end

            plotHandle = AssessmentClassificationResultPlotHandle(rectangleHandle, textHandle, symbolHandle);
        end

        function textStr = getTextForPrediction(obj,truthClass,predictedClass)
            predictedClassStr = Labeling.StringForClass(predictedClass,obj.classNames);

            if truthClass ~= predictedClass
                truthClassStr = Labeling.StringForClass(truthClass,obj.classNames);
                textStr = sprintf('%s\n (truth: %s)',predictedClassStr,truthClassStr);
            else
                textStr = predictedClassStr;
            end
        end

        function colorStr = getColorForPrediction(~,truthClass,predictedClass)
            if truthClass == predictedClass
                colorStr = Constants.kCorrectColor;
            else
                colorStr = Constants.kWrongColor;
            end
        end

        function plotMissedEventAnnotations(obj,signal)
            for i = 1 : length(obj.missedEvents)
                missedEvent = obj.missedEvents(i);
                if missedEvent.label ~= Labeling.kNullClass
                    obj.plotMissedEventAnnotation(missedEvent,signal);
                end
            end
        end

        function plotMissedEventAnnotation(obj,event,signal)
            eventX = event.sample;
            eventY = signal(eventX);
            class = event.label;

            classStr = Helper.StringForClass(class,obj.classNames);

            plot(obj.plotAxes,eventX,eventY,'*','Color',Constants.kMissedColor,'LineWidth',AssessmentClassificationResultsPlotter.MissedEventSymbolSize);
            textHandle = text(obj.plotAxes,double(eventX),double(eventY), classStr,'FontSize',AssessmentClassificationResultsPlotter.FontSize);
            set(textHandle, 'Clipping', 'on');
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function convex_hull ( prefix )

%*****************************************************************************80
%
%% CONVEX_HULL: an algorithm for the convex hull of a set of 2D points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CONVEX_HULL:\n' );
  fprintf ( 1, '  MATLAB version.\n' );

  if ( nargin < 1 )
    prefix = 'random points'
    n = 25;
    xy = rand ( n, 2 );
  else
    node_filename = strcat ( prefix, '_nodes.txt' );
    xy = load ( node_filename );
  end
%
%  Plot the initial data.
%
  figure ( 1 )
  clf
  plot ( xy(:,1), xy(:,2), 'b.', 'MarkerSize', 20 );
  grid on
  axis tight
  axis equal
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  title ( prefix );
%
%  Plot the convex hull as computed by MATLAB.
%
  figure ( 2 )
  clf
  k = convhull ( xy(:,1), xy(:,2) );
  plot ( xy(:,1), xy(:,2), 'b.', 'MarkerSize', 20 );
  hold on
  plot( xy(k,1), xy(k,2), 'r-', 'LineWidth', 2 );
  grid on
  axis tight
  axis equal
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  title ( sprintf ( '%s by MATLAB convhull', prefix ) );
  hold off
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CONVEX_HULL:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  return
end
function timestamp ( )

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

  return
end

repo: jburkardt-m


text:
function [eta, Heta, inner_it, stop_tCG] ...
                 = tCG(problem, x, grad, eta, Delta, options, storedb, key)
% tCG - Truncated (Steihaug-Toint) Conjugate-Gradient method
% minimize <eta,grad> + .5*<eta,Hess(eta)>
% subject to <eta,eta>_[inverse precon] <= Delta^2
%
% See also: trustregions

% This file is part of Manopt: www.manopt.org.
% This code is an adaptation to Manopt of the original GenRTR code:
% RTR - Riemannian Trust-Region
% (c) 2004-2007, P.-A. Absil, C. G. Baker, K. A. Gallivan
% Florida State University
% School of Computational Science
% (http://www.math.fsu.edu/~cbaker/GenRTR/?page=download)
% See accompanying license file.
% The adaptation was executed by Nicolas Boumal.
%
% Change log:
%
%   NB Feb. 12, 2013:
%       We do not project r back to the tangent space anymore: it was not
%       necessary, and as of Manopt 1.0.1, the proj operator does not
%       coincide with this notion anymore.
%
%   NB April 3, 2013:
%       tCG now also returns Heta, the Hessian at x along eta. Additional
%       esthetic modifications.
%
%   NB Dec. 2, 2013:
%       If options.useRand is activated, we now make sure the preconditio-
%       ner is not used, as was originally intended in GenRTR. In time, we
%       may want to investigate whether useRand can be modifed to work well
%       with preconditioning too.
%
%   NB Jan. 9, 2014:
%       Now checking explicitly for model decrease at each iteration. The
%       first iteration is a Cauchy point, which necessarily realizes a
%       decrease of the model cost. If a model increase is witnessed
%       (which is theoretically impossible if a linear operator is used for
%       the Hessian approximation), then we return the previous eta. This
%       ensures we always achieve at least the Cauchy decrease, which
%       should be sufficient for convergence.
%
%   NB Feb. 17, 2015:
%       The previous update was in effect verifying that the current eta
%       performed at least as well as the first eta (the Cauchy step) with
%       respect to the model cost. While this is an acceptable strategy,
%       the documentation (and the original intent) was to ensure a
%       monotonic decrease of the model cost at each new eta. This is now
%       the case, with the added line: "model_value = new_model_value;".
%
%   NB April 3, 2015:
%       Works with the new StoreDB class system.

% All terms involving the trust-region radius will use an inner product
% w.r.t. the preconditioner; this is because the iterates grow in
% length w.r.t. the preconditioner, guaranteeing that we will not
% re-enter the trust-region.
%
% The following recurrences for Prec-based norms and inner
% products come from [CGT2000], pg. 205, first edition.
% Below, P is the preconditioner.
%
% <eta_k,P*delta_k> =
%          beta_k-1 * ( <eta_k-1,P*delta_k-1> + alpha_k-1 |delta_k-1|^2_P )
% |delta_k|^2_P = <r_k,z_k> + beta_k-1^2 |delta_k-1|^2_P
%
% therefore, we need to keep track of
% 1)   |delta_k|^2_P
% 2)   <eta_k,P*delta_k> = <eta_k,delta_k>_P
% 3)   |eta_k  |^2_P
%
% initial values are given by:
%    |delta_0|_P = <r,z>
%    |eta_0|_P   = 0
%    <eta_0,delta_0>_P = 0
% because we take eta_0 = 0 (if useRand = false).
%
% [CGT2000] Conn, Gould and Toint: Trust-region methods, 2000.

inner = problem.M.inner;
lincomb = problem.M.lincomb;

theta = options.theta;
kappa = options.kappa;

if ~options.useRand % and therefore, eta == 0
    Heta = problem.M.zerovec(x);
    r = grad;
    e_Pe = 0;
else % and therefore, no preconditioner
    % eta (presumably) ~= 0 was provided by the caller.
    Heta = getHessian(problem, x, eta, storedb, key);
    r = lincomb(x, 1, grad, 1, Heta);
    e_Pe = inner(x, eta, eta);
end
r_r = inner(x, r, r);
norm_r = sqrt(r_r);
norm_r0 = norm_r;

% Precondition the residual.
if ~options.useRand
    z = getPrecon(problem, x, r, storedb, key);
else
    z = r;
end

% Compute z'*r.
z_r = inner(x, z, r);
d_Pd = z_r;

% Initial search direction.
delta  = lincomb(x, -1, z);
if ~options.useRand % and therefore, eta == 0
    e_Pd = 0;
else % and therefore, no preconditioner
    e_Pd = inner(x, eta, delta);
end

% If the Hessian or a linear Hessian approximation is in use, it is
% theoretically guaranteed that the model value decreases strictly
% with each iteration of tCG. Hence, there is no need to monitor the model
% value. But, when a nonlinear Hessian approximation is used (such as the
% built-in finite-difference approximation for example), the model may
% increase. It is then important to terminate the tCG iterations and return
% the previous (the best-so-far) iterate. The variable below will hold the
% model value.
model_fun = @(eta, Heta) inner(x, eta, grad) + .5*inner(x, eta, Heta);
if ~options.useRand
    model_value = 0;
else
    model_value = model_fun(eta, Heta);
end

% Pre-assume termination because j == end.
stop_tCG = 5;

% Begin inner/tCG loop.
j = 0;
for j = 1 : options.maxinner

    % This call is the computationally expensive step.
    Hdelta = getHessian(problem, x, delta, storedb, key);

    % Compute curvature (often called kappa).
    d_Hd = inner(x, delta, Hdelta);

    % Note that if d_Hd == 0, we will exit at the next "if" anyway.
    alpha = z_r/d_Hd;
    % <neweta,neweta>_P =
    % <eta,eta>_P + 2*alpha*<eta,delta>_P + alpha*alpha*<delta,delta>_P
    e_Pe_new = e_Pe + 2.0*alpha*e_Pd + alpha*alpha*d_Pd;

    if options.debug > 2,
        fprintf('DBG:   (r,r)  : %e\n', r_r);
        fprintf('DBG:   (d,Hd) : %e\n', d_Hd);
        fprintf('DBG:   alpha  : %e\n', alpha);
    end

    % Check against negative curvature and trust-region radius violation.
    % If either condition triggers, we bail out.
    if d_Hd <= 0 || e_Pe_new >= Delta^2,
        % want
        %  ee = <eta,eta>_prec,x
        %  ed = <eta,delta>_prec,x
        %  dd = <delta,delta>_prec,x
        tau = (-e_Pd + sqrt(e_Pd*e_Pd + d_Pd*(Delta^2-e_Pe))) / d_Pd;
        if options.debug > 2,
            fprintf('DBG:     tau  : %e\n', tau);
        end
        eta  = lincomb(x, 1,  eta, tau,  delta);

        % If only a nonlinear Hessian approximation is available, this is
        % only approximately correct, but saves an additional Hessian call.
        Heta = lincomb(x, 1, Heta, tau, Hdelta);

        % Technically, we may want to verify that this new eta is indeed
        % better than the previous eta before returning it (this is always
        % the case if the Hessian approximation is linear, but I am unsure
        % whether it is the case or not for nonlinear approximations.)
        % At any rate, the impact should be limited, so in the interest of
        % code conciseness (if we can still hope for that), we omit this.

        if d_Hd <= 0,
            stop_tCG = 1;     % negative curvature
        else
            stop_tCG = 2;     % exceeded trust region
        end
        break;
    end

    % No negative curvature and eta_prop inside TR: accept it.
    e_Pe = e_Pe_new;
    new_eta  = lincomb(x, 1,  eta, alpha,  delta);

    % If only a nonlinear Hessian approximation is available, this is
    % only approximately correct, but saves an additional Hessian call.
    new_Heta = lincomb(x, 1, Heta, alpha, Hdelta);

    % Verify that the model cost decreased in going from eta to new_eta. If
    % it did not (which can only occur if the Hessian approximation is
    % nonlinear or because of numerical errors), then we return the
    % previous eta (which necessarily is the best reached so far, according
    % to the model cost). Otherwise, we accept the new eta and go on.
    new_model_value = model_fun(new_eta, new_Heta);
    if new_model_value >= model_value
        stop_tCG = 6;
        break;
    end

    eta = new_eta;
    Heta = new_Heta;
    model_value = new_model_value; %% added Feb. 17, 2015

    % Update the residual.
    r = lincomb(x, 1, r, alpha, Hdelta);

    % Compute new norm of r.
    r_r = inner(x, r, r);
    norm_r = sqrt(r_r);

    % Check kappa/theta stopping criterion.
    % Note that it is somewhat arbitrary whether to check this stopping
    % criterion on the r's (the gradients) or on the z's (the
    % preconditioned gradients). [CGT2000], page 206, mentions both as
    % acceptable criteria.
    if j >= options.mininner && norm_r <= norm_r0*min(norm_r0^theta, kappa)
        % Residual is small enough to quit
        if kappa < norm_r0^theta,
            stop_tCG = 3;  % linear convergence
        else
            stop_tCG = 4;  % superlinear convergence
        end
        break;
    end

    % Precondition the residual.
    if ~options.useRand
        z = getPrecon(problem, x, r, storedb, key);
    else
        z = r;
    end

    % Save the old z'*r.
    zold_rold = z_r;
    % Compute new z'*r.
    z_r = inner(x, z, r);

    % Compute new search direction.
    beta = z_r/zold_rold;
    delta = lincomb(x, -1, z, beta, delta);

    % Update new P-norms and P-dots [CGT2000, eq. 7.5.6 & 7.5.7].
    e_Pd = beta*(e_Pd + alpha*d_Pd);
    d_Pd = z_r + beta*beta*d_Pd;

end  % of tCG loop
inner_it = j;

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial