text:
function [pdc, pdcvar, n] = ft_connectivity_pdc(input, varargin)

% FT_CONNECTIVITY_PDC computes partial directed coherence. This function implements
% the metrices described in Baccala et al., Biological Cybernetics 2001, 84(6),
% 463-74. and in Baccala et al., 15th Int.Conf.on DSP 2007, 163-66.
%
% The implemented algorithm has been tested against the implementation in the
% SIFT-toolbox. It yields numerically identical results to what is known there as
% 'nPDC' (for PDC) and 'GPDC' for generalized pdc.
%
% Use as
%   [p, v, n] = ft_connectivity_pdc(h, key1, value1, ...)
%
% The input argument H should be a spectral transfer matrix organized as
%   Nrpt x Nchan x Nchan x Nfreq (x Ntime),
% where Nrpt can be 1.
%
% Additional optional input arguments come as key-value pairs:
%   'hasjack'  = 0 (default) is a boolean specifying whether the input
%                contains leave-one-outs, required for correct variance
%                estimate
%   'feedback' = string, determining verbosity (default = 'none'), see FT_PROGRESS
%   'invfun'   = 'inv' (default) or 'pinv', the function used to invert the
%                transfer matrix to obtain the fourier transform of the
%                MVAR coefficients. Use 'pinv' if the data are
%                poorly-conditioned.
%   'noisecov' = matrix containing the covariance of the residuals of the
%                MVAR model. If this matrix is defined, the function
%                returns the generalized partial directed coherence.
%
% Output arguments:
%   p = partial directed coherence matrix Nchan x Nchan x Nfreq (x Ntime).
%       If multiple observations in the input, the average is returned.
%   v = variance of p across observations.
%   n = number of observations.
%
% Typically, nrpt should be 1 (where the spectral transfer matrix is
% computed across observations. When nrpt>1 and hasjack is true the input
% is assumed to contain the leave-one-out estimates of H, thus a more
% reliable estimate of the relevant quantities.
%
% See also FT_CONNECTIVITYANALYSIS

% Copyright (C) 2009-2017, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

hasjack  = ft_getopt(varargin, 'hasjack', 0);
feedback = ft_getopt(varargin, 'feedback', 'none');
invfun   = ft_getopt(varargin, 'invfun', 'inv');
noisecov = ft_getopt(varargin, 'noisecov');

switch invfun
  case {'inv' 'pinv'}
    invfun = str2func(invfun);
  otherwise
    ft_error('unknown specification of inversion-function for the transfer matrix');
end

% crossterms are described by chan_chan_therest
siz = [size(input) 1];
n   = siz(1);

if ~isempty(noisecov)
  ft_progress('init', feedback, 'computing generalized partial directed coherence...');
  scale = diag(sqrt(1./diag(noisecov))); % get the 1./sqrt(var) of the signals
else
  ft_progress('init', feedback, 'computing partial directed coherence...');
  scale = eye(siz(2));
end

% pre-allocate some variables
outsum = zeros(siz(2:end));
outssq = zeros(siz(2:end));

% the mathematics for pdc is most straightforward using the inverse of the
% transfer function
pdim     = prod(siz(4:end));
tmpinput = reshape(input, [siz(1:3) pdim]);
ft_progress('init', feedback, 'inverting the transfer function...');
for k = 1:n
  ft_progress(k/n, 'inverting the transfer function for replicate %d from %d\n', k, n);
  tmp = reshape(tmpinput(k,:,:,:), [siz(2:3) pdim]);
  for m = 1:pdim
    tmp(:,:,m) = scale*invfun(tmp(:,:,m));
  end
  tmpinput(k,:,:,:) = tmp;
end
ft_progress('close');
input = reshape(tmpinput, siz);

for j = 1:n
  ft_progress(j/n, 'computing metric for replicate %d from %d\n', j, n);
  invh   = reshape(input(j,:,:,:,:), siz(2:end));

  den    = sum(abs(invh).^2,1);
  tmppdc = abs(invh)./sqrt(repmat(den, [siz(2) 1 1 1 1]));

  outsum = outsum + tmppdc;
  outssq = outssq + tmppdc.^2;
end
ft_progress('close');

pdc = outsum./n;

if n>1
  if hasjack
    bias = (n-1).^2;
  else
    bias = 1;
  end
  pdcvar = bias*(outssq - (outsum.^2)./n)./(n - 1);
else
  pdcvar = [];
end

% this is to achieve the same convention for all connectivity metrics:
% row -> column
for k = 1:prod(siz(4:end))
  pdc(:,:,k)    = transpose(pdc(:,:,k));
  if ~isempty(pdcvar)
    pdcvar(:,:,k) = transpose(pdcvar(:,:,k));
  end
end


text:
function fe2dx_r_fast_test ( )

%*****************************************************************************80
%
%% FE2DX_R_FAST_TEST tests the FE2DX_R_FAST code.
%
%  Discussion:
%
%    This function sets all parameter values and initial condition information
%    necessary to execute the "fast" version of the fe2dx_r algorithm.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    28 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Reference:
%
%    Marcus R Garvie, John Burkardt, Jeff Morgan,
%    Simple Finite Element Methods for Approximating Predator-Prey Dynamics
%    in Two Dimensions using MATLAB,
%    Submitted to Bulletin of Mathematical Biology, 2014.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FE2DX_R_FAST_TEST:\n' );
  fprintf ( 1, '  Test the FE2DX_R_FAST function\n' );
  fprintf ( 1, '  which applies Robin boundary conditions as it\n' );
  fprintf ( 1, '  approximates a solution to a predator-prey system.\n' );
%
%  Set the parameters.
%
  alpha = 0.4;
  beta = 2.0;
  gamma = 0.6;
  delta = 1.0;
%
%  Use T=150.0 for normal run.
%  Use T=0.50 for a "quick" run that might take 15 minutes of computing.
%
% T = 150.0;
  T = 0.50;
  delt = 1.0 / 384.0;
  k1 = 0.01;
  k2 = 0.01;

  t = tic;
  fe2dx_r_fast ( alpha, beta, gamma, delta, T, delt, @u0f, @v0f, k1, k2 );
  t = toc ( t );

  fprintf ( 1, '  Execution took %10.2g minutes \n', t / 60.0 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FE2DX_R_FAST_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

function value = u0f ( x, y )

%*****************************************************************************80
%
%% U0F evaluates the initial condition for U.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    26 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location in the region.
%
%    Output, real VALUE, the initial condition for U at (X,Y).
%
  value = 6.0 / 35.0 - 2.0E-07 * ( x - 0.1 * y - 225.0 ) * ( x - 0.1 * y - 675.0 );

  return
end

function value = v0f ( x, y )

%*****************************************************************************80
%
%% V0F evaluates the initial condition for V.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    26 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location in the region.
%
%    Output, real VALUE, the initial condition for V at (X,Y).
%
  value = 116.0 / 245.0 - 3.0E-05 * ( x - 450.0 ) - 1.2E-04 * ( y - 150.0 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef MapNested < containers.Map & handle

    % A nested map container
    %
    % A MapNested object implements nested maps (map of maps).
    %
    % MapN is a handle class.
    % Since MapNested is a subclass from containers.Map, all functions from
    % containers.Map will be inherited (eg. keys, values, ...)
    %
    % Description - basic outline
    % ---------------------------
    %
    % A MapNested object is constructed like this:
    %
    %   M = MapNested();
    %
    % Values are stored using M(key1, key2, ...) = value, for example:
    %
    %   M(1, 'a')     = 'a string value';
    %   M(1, 'b')     = 287.2;
    %   M(2)          = [1 2 3; 4 5 6];
    %
    % Attention: you can not assign cell-arrays in this implementation:
    %   M(2, 'x', pi) = {'a' 'cell' 'array'}; %throws error
    %
    % another possibility is to define the keys as cell-array like:
    %   M({key1, key2, key3}) = value;
    %
    %
    % Values are retrieved using M(key1, key2, ...), for example
    %
    %   v = M(1, 'b');
    %   u = M(2);
    %
    % or with using cell-arrays for the keys
    %   v = M({key1, key2, key3});
    %
    % Set and get - methods
    % -----------------------------
    %
    % for setting and retrieving values there are also two methods
    % implemented, for setting a value:
    %
    %   MapObj = setValueNested(MapObj, {key1, key2, key3, ...}, value);
    %
    %   here the second input parameter has to be a cell-array with the
    %   keys.
    %
    % For retrieving values one can use:
    %   value = getValueNested(MapObj, {key1, key2, key3,...});
    %
    %   here the second input parameter has to be a cell-array with the
    %   keys.
    %
    %
    % Updating and removing entries
    % -----------------------------
    %
    % The value for a given key list is updated using the usual assignment;
    % the previous value is overwritten.
    %
    %   M(pi, 'x') = 1;     % 1 is current value
    %   M(pi, 'x') = 2;     % 2 replaces 1 as the value for this key list
    %
    %   NMapobj('Key') = []; %removes the key
    %   remove(NMapobj, 'Key'); %also removes the key
    %
    % check if is a specific key in the map
    % -------------------------------------
    %   M.isKey(keyList) %keyList is a cell array of the nested keys
    %   M.isKey(Key1, Key2, Key3,..) %Key1, Key2, Key3 are the nested keys
    %   isKey(M, keyList) %keyList is a cell array of the nested keys
    %   isKey(M, Key1, Key2, Key3,..) %Key1, Key2, Key3 are the nested keys
    %
    %
    % Methods and properties
    % ----------------------
    %
    % MapN methods:
    %   MapNested        - constructor for MapNested objects
    %   subsref     - implements value = Mobj(keylist)
    %   subsasgn    - implements M(keylist) = value
    %   setValueNested  - implements Mobj = setValueNested(Mobj, keyList,
    %   value)
    %   getValueNested  - implements value = getValueNested(Mobj, keyList);
    %   isKey - implements Mobj.iskey(keyList);
    %
    % See also: containers.Map

    % (c) 2017, Roland Ritt
    methods
        function obj = MapNested(varargin)
            % constructor, calls Superclass-constructor with varargin;
            obj = obj@containers.Map(varargin{:});

        end

        function obj = setValueNested(obj, keyList, value)
            % method which recursively constructs a map of maps and add
            % a a value at the specified keys
            %
            % Implements the syntax
            %
            %   MapNestedObj = setValueNested(MapNestedObj, keyList,
            %   value);
            %       (keyList is of type 'CellArray')
            %
            % See also: MapNested, MapNested/setValueNested
            if ~iscell(keyList)
                %check if the keyList is a cell array
                error('second input (keylist) is no cellArray');
            end

            KeyType = class(keyList{1});
            %find the class-type of first key
            if length(keyList)==1
                % check if the object is defined
                if isempty(obj)
                    %if the object does not exist...
                    obj = MapNested('KeyType', KeyType,'ValueType', 'any');
                    % generate a new MapNested-object with the specified
                    % key-Type, and default ValueType set to 'any'

                end

                if isempty(value)
                    remove(obj, keyList{1});
                else
                obj = [obj;MapNested(keyList{1},value)];
                end

                return
            else
                % if more than one key
                if obj.isKey(keyList{1})
                    % check if the key is in the map
                    temp = values(obj, keyList(1));
                    temp = temp{1};
                    % retrieve the MapNested-object from the map
                    if ~isa(temp ,'containers.Map')
                        % if the value is not a MapObject, generate a new
                        % one
                        temp = MapNested('KeyType', KeyType,'ValueType', 'any');
                    end

                else
                    % if the key does not exist, generate a new
                    % MapNested-object
                    temp = MapNested('KeyType', KeyType,'ValueType', 'any');

                end
                temp = setValueNested(temp, keyList(2:end), value);
                % set the value in the MapNested-object (recursive call)

                if isempty(obj)
                    %if obj is empty, generate a new MapNested-object
                    obj = MapNested(keyList{1}, temp);
                else
                    %if obj exists, concatenate the new entry
                    obj = [obj;MapNested(keyList{1}, temp)];
                end

            end

        end

        function value = getValueNested(obj, keyList)
            % method for retrieving values from a MapNested object by
            % recursively calls to this method
            %
            % Implements the syntax
            %
            %   value = getValueNested(MapNestedObj, keyList, value);
            %       (keyList is of type 'CellArray')
            %
            % See also: MapNested, MapNested/setValueNested
            if ~iscell(keyList)
                %check if the keyList is a cellArray
                error('second input is no cellArray');
            end

            if ~obj.isKey(keyList{1})
                %check if the first key is in the list, if not, throw an
                %error (maybe should be changed to empty object)
                error(['key ''', keyList{1}, ''' is not a key'] );
            end
            if length(keyList)==1
                % if the keyList is only of length 1, return the value of
                % this list
                value = values(obj, {keyList{1}});
                % call the method from the superclass
                value=value{1}; %retrieve the value
                return
            else
                % if there are more than one keys Left in the keyList
                temp = values(obj, {keyList{1}});
                % retrieve the MapNested object (key1)
                temp = temp{1};
                if ~isa(temp ,'containers.Map')
                    % if the retrieved value is not a Map, the
                    % MapNested object was misused --> throw error message
                    error(['key ''', keyList{2}, ''' is not a key'] );
                end
                value = getValueNested(temp, keyList(2:end));
                % recursively call the getValueNested method.

            end

        end

        function v = subsref(M, S)
            % returns value associated with key list
            %
            % Implements the syntax
            %
            %   value = MapNobj(key1, key2, ...)
            %
            % See also: MapNested, MapNested/subsasgn

            switch S(1).type
                case '.'
                    v = builtin('subsref',M,S);

                case '()'
                    %                     if ~isscalar(S) || ~strcmp(S.type, '()') || length(S)<1
                    %                         error('MapNested:Subsref:LimitedIndexing', ...
                    %                             'Only ''()'' indexing is supported by a MapNested');
                    %                     end

                    try

                        if iscell(S(1).subs{1})
                            temp = S(1).subs{1};
                        else
                            temp = S(1).subs;
                        end
                        v = getValueNested(M, temp);
                    catch me
                        % default is handled in subsrefError for efficiency
                        error('MapNested:Subsref:IndexingError', ...
                            ['Something went wrong in indexing: ',me(1).message] );
                    end

                    if length(S)>1
                        v = subsref(v, S(2:end));
                    end

                case '{}'
                    error('MapNested:Subsasgn:LimitedIndexing', ...
                        '''{}'' indexing is not supported by a MapNested');
                otherwise
                    error('Not a valid indexing expression')
            end

        end

        function M = subsasgn(M, S, v)
            % sets value associated with key list
            %
            % Implements the syntax
            %
            %   MapNobj(key1, key2, ...) = value
            %
            % See also: MapNested, MapNested/subsasgn

            if ~isscalar(S) || ~strcmp(S.type, '()') || length(S)<1
                error('MapNested:Subsasgn:LimitedIndexing', ...
                    'Only ''()'' indexing is supported by a MapNested');
            end

            try

                if iscell(S.subs{1})
                    temp = S.subs{1};
                else
                    temp = S.subs;
                end
                M = setValueNested(M, temp, v);
            catch me
                % default is handled in subsrefError for efficiency
                error('MapNested:Subsasgn:IndexingError', ...
                    'Something went wrong in indexing');
            end
        end

        function bisKey = isKey(M, varargin)

            %% test

            if length(varargin)==1
                keyList = varargin{1};
            else
                keyList = varargin;
            end

            if iscell(keyList)
                if length(keyList)== 1
                    bisKey = isKey@containers.Map(M,keyList{1});
                    return
                else

                    KeyTemp = keyList(1);
                    if M.isKey(KeyTemp)
                        Mtemp = M.getValueNested(KeyTemp);
                        bisKey = isKey(Mtemp, keyList(2:end));
                    else
                        bisKey = false;
                    end
                    return;
                end

            else
                bisKey = isKey@containers.Map(M,varargin{1});
                return;
            end

        end

    end

end

author: anicolson

repo: DeepXi

save_path: github-repos/MATLAB/anicolson-DeepXi


text:
function delay=delayest_3point(u2,u1,method,estimator,parameter);
%delay=delayest_3point(u2,u1,method,estimator,parameter);
%Estimates delay of u2 wrt u1 by interpolating the peak of the cross correlatio
%using a three-point peak interpolation
%method may be 'parabola','Gaussian','modGaussian','cosine'.
%estimator may be 'xcorr','ASDF','AMDF' (AMDF is very slow).
%if using modGaussian, parameter is a bias (always used if >=0, only used
%for xc<0 for parameter < 0)

if nargin<3
    method='parabola';
end

if nargin<4
    estimator='xcorr';
end

if nargin<5
    switch method
        case {'modGaussian','modgaussian'}
            parameter=1;
    end
end

N_p=numel(u2);%number of elements

switch estimator
    case {'xcorr','xc','xcorr_fft'}%cross correlation
        U1=fft(u1);
        U2=fft(u2);
        xc=ifft(U2.*conj(U1));%circular xcorr
        [tmp idx]=max(xc);
    case {'ASDF'}%average squared difference function
         xc = (-2*ifft(fft(u2).*conj(fft(u1))) + sum(u1.^2) + sum(u2.^2))/N_p;%ADSF
        [tmp idx]=min(xc);
    case {'AMDF'}%average magnitude difference function (this is a lot slower for large N_p)
        xc=sum(abs(repmat(u1',1,N_p)-hankel(u2',[u2(end)'; u2(1:(end-1))'])))/N_p;%AMDF
        [tmp idx]=min(xc);
end

R=zeros(1,3);%three points to use for interpolation
R(2)=xc(idx);%peak

%find neighbors, assuming periodic
if idx==1
    R(1)=xc(end);
    R(3)=xc(2);
elseif idx==numel(xc)
    R(1)=xc(end-1);
    R(3)=xc(1);
else
    R(1)=xc(idx-1);
    R(3)=xc(idx+1);
end

switch method
    case {'Gaussian','gaussian'}
        %a=exp(log(R(2))+(log(R(3))-log(R(1)))^2/(16*log(R(2))-8*log(R(1))-8*log(R(3))));%peak value
        %b=(2*log(R(2))-log(R(1))-log(R(3)))/2;%width
        c=(log(R(3))-log(R(1)))/(4*log(R(2))-2*log(R(1))-2*log(R(3)));
    case {'parabola','parabolic'}
    	%A=inv([1 -1 1;0 0 1;1 1 1])*R';
        %c=-A(2)/(2*A(1));
        %a=-1/4*A(2)^2/A(1)+A(3);%peak value
        c=(R(3)-R(1))/(2*(2*R(2)-R(1)-R(3)));
    case {'modGaussian','modgaussian'}
        if parameter>=0%always add bias
            R=R-min(R)+parameter*R(2);
        else%only add bias if R has a negative value
            if any(R<=0)
                R=R-min(R)-parameter*R(2);
            end
        end
        %a=exp(log(R(2))+(log(R(3))-log(R(1)))^2/(16*log(R(2))-8*log(R(1))-8*log(R(3))));
        %b=(2*log(R(2))-log(R(1))-log(R(3)))/2;
        c=(log(R(3))-log(R(1)))/(4*log(R(2))-2*log(R(1))-2*log(R(3)));
    case {'cosine'}
        omega=acos((R(1)+R(3))/(2*R(2)));
        theta=atan((R(1)-R(3))/(2*R(2)*sin(omega)));
        %A=R(2)/cos(theta);
        c=-theta/omega;
    otherwise
        error('unknown method')
end

lag=mod(idx-1+floor(N_p/2),N_p)-floor(N_p/2);%integer part
delay=lag+c;%delay estimate

author: Sable

repo: mcbench-benchmarks


text:
function prediction = RUSBoost (TRAIN,TEST,WeakLearn)
% This function implements the RUSBoost Algorithm. For more details on the
% theoretical description of the algorithm please refer to the following
% paper:
% C. Seiffert, T.M. Khoshgoftaar, J. Van Hulse and A. Napolitano,
% "RUSBoost: A Hybrid Approach to Alleviating Class Imbalance, IEEE
% Transaction on Systems, Man and Cybernetics-Part A: Systems and Human,
% Vol.40(1), January 2010.
% Input: TRAIN = Training data as matrix
%        TEST = Test data as matrix
%        WeakLearn = String to choose algortihm. Choices are
%                    'svm','tree','knn' and 'logistic'.
% Output: prediction = size(TEST,1)x 2 matrix. Col 1 is class labels for
%                      all instances. Col 2 is probability of the instances
%                      being classified as positive class.

javaaddpath('weka.jar');

%% Training RUSBoost
% Total number of instances in the training set
m = size(TRAIN,1);
POS_DATA = TRAIN(TRAIN(:,end)==1,:);
NEG_DATA = TRAIN(TRAIN(:,end)==0,:);
pos_size = size(POS_DATA,1);
neg_size = size(NEG_DATA,1);

% Reorganize TRAIN by putting all the positive and negative exampels
% together, respectively.
TRAIN = [POS_DATA;NEG_DATA];

% Converting training set into Weka compatible format
CSVtoARFF (TRAIN, 'train', 'train');
train_reader = javaObject('java.io.FileReader', 'train.arff');
train = javaObject('weka.core.Instances', train_reader);
train.setClassIndex(train.numAttributes() - 1);

% 65% of NEG_DATA
neg65 = round(pos_size * (65/35));

% Total number of iterations of the boosting method
T = 10;

% W stores the weights of the instances in each row for every iteration of
% boosting. Weights for all the instances are initialized by 1/m for the
% first iteration.
W = zeros(1,m);
for i = 1:m
    W(1,i) = 1/m;
end

% L stores pseudo loss values, H stores hypothesis, B stores (1/beta)
% values that is used as the weight of the % hypothesis while forming the
% final hypothesis. % All of the following are of length <=T and stores
% values for every iteration of the boosting process.
L = [];
H = {};
B = [];

% Loop counter
t = 1;

% Keeps counts of the number of times the same boosting iteration have been
% repeated
count = 0;

% Boosting T iterations
while t <= T

    % LOG MESSAGE
    disp (['Boosting iteration #' int2str(t)]);

    % Resampling NEG_DATA with weights of positive example
    RESAM_NEG = NEG_DATA(randsample(1:neg_size,neg65,false),:);
    RESAMPLED = [POS_DATA;RESAM_NEG];

    % Converting resample training set into Weka compatible format
    CSVtoARFF (RESAMPLED,'resampled','resampled');
    reader = javaObject('java.io.FileReader','resampled.arff');
    resampled = javaObject('weka.core.Instances',reader);
    resampled.setClassIndex(resampled.numAttributes()-1);

    % Training a weak learner. 'pred' is the weak hypothesis. However, the
    % hypothesis function is encoded in 'model'.
    switch WeakLearn
        case 'svm'
            model = javaObject('weka.classifiers.functions.SMO');
        case 'tree'
            model = javaObject('weka.classifiers.trees.J48');
        case 'knn'
            model = javaObject('weka.classifiers.lazy.IBk');
            model.setKNN(5);
        case 'logistic'
            model = javaObject('weka.classifiers.functions.Logistic');
    end
    model.buildClassifier(resampled);

    pred = zeros(m,1);
    for i = 0 : m - 1
        pred(i+1) = model.classifyInstance(train.instance(i));
    end

    % Computing the pseudo loss of hypothesis 'model'
    loss = 0;
    for i = 1:m
        if TRAIN(i,end)==pred(i)
            continue;
        else
            loss = loss + W(t,i);
        end
    end

    % If count exceeds a pre-defined threshold (5 in the current
    % implementation), the loop is broken and rolled back to the state
    % where loss > 0.5 was not encountered.
    if count > 5
       L = L(1:t-1);
       H = H(1:t-1);
       B = B(1:t-1);
       disp ('Too many iterations have loss > 0.5');
       disp ('Aborting boosting...');
       break;
    end

    % If the loss is greater than 1/2, it means that an inverted
    % hypothesis would perform better. In such cases, do not take that
    % hypothesis into consideration and repeat the same iteration. 'count'
    % keeps counts of the number of times the same boosting iteration have
    % been repeated
    if loss > 0.5
        count = count + 1;
        continue;
    else
        count = 1;
    end

    L(t) = loss; % Pseudo-loss at each iteration
    H{t} = model; % Hypothesis function
    beta = loss/(1-loss); % Setting weight update parameter 'beta'.
    B(t) = log(1/beta); % Weight of the hypothesis

    % At the final iteration there is no need to update the weights any
    % further
    if t==T
        break;
    end

    % Updating weight
    for i = 1:m
        if TRAIN(i,end)==pred(i)
            W(t+1,i) = W(t,i)*beta;
        else
            W(t+1,i) = W(t,i);
        end
    end

    % Normalizing the weight for the next iteration
    sum_W = sum(W(t+1,:));
    for i = 1:m
        W(t+1,i) = W(t+1,i)/sum_W;
    end

    % Incrementing loop counter
    t = t + 1;
end

% The final hypothesis is calculated and tested on the test set
% simulteneously.

%% Testing RUSBoost
n = size(TEST,1); % Total number of instances in the test set

CSVtoARFF(TEST,'test','test');
test = 'test.arff';
test_reader = javaObject('java.io.FileReader', test);
test = javaObject('weka.core.Instances', test_reader);
test.setClassIndex(test.numAttributes() - 1);

% Normalizing B
sum_B = sum(B);
for i = 1:size(B,2)
   B(i) = B(i)/sum_B;
end

prediction = zeros(n,2);

for i = 1:n
    % Calculating the total weight of the class labels from all the models
    % produced during boosting
    wt_zero = 0;
    wt_one = 0;
    for j = 1:size(H,2)
       p = H{j}.classifyInstance(test.instance(i-1));
       if p==1
           wt_one = wt_one + B(j);
       else
           wt_zero = wt_zero + B(j);
       end
    end

    if (wt_one > wt_zero)
        prediction(i,:) = [1 wt_one];
    else
        prediction(i,:) = [0 wt_one];
    end
end


text:
% Random directed graph construction
% Note 1: if p is omitted, p=0.5 is default
% Note 2: no self-loops, no double edges
%
% INPUTS:  n - number of nodes
%          p - probability, 0<=p<=1
% Output: adjacency matrix, nxn
%
% GB: last updated, Oct 21 2012

function adj = randomDirectedGraph(n,p)

adj=zeros(n); % initialize adjacency matrix

if nargin==1; p=0.5; end; % default probability

% splitting j = 1:i-1,i+1,n avoids the if statement i==j

for i=1:n

  for j=1:i-1
    if rand<=p; adj(i,j)=1; end;
  end

  for j=i+1:n
    if rand<=p; adj(i,j)=1; end;
  end

end

author: aeolianine

repo: octave-networks-toolbox


text:
function [D, ni] = find_nn(X, k)
%FIND_NN Finds k nearest neigbors for all datapoints in the dataset
%
%	[D, ni] = find_nn(X, k)
%
% Finds the k nearest neighbors for all datapoints in the dataset X.
% In X, rows correspond to the observations and columns to the
% dimensions. The value of k is the number of neighbors that is
% stored. The function returns a sparse distance matrix D, in which
% only the distances to the k nearest neighbors are stored. For
% equal datapoints, the distance is set to a tolerance value.
% The method is relatively slow, but has a memory requirement of O(nk).
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

	if ~exist('k', 'var') || isempty(k)
		k = 12;
    end

    % Perform adaptive neighborhood selection if desired
    if ischar(k)
        [D, max_k] = find_nn_adaptive(X);
        ni = zeros(size(X, 1), max_k);
        for i=1:size(X, 1)
            tmp = find(D(i,:) ~= 0);
            tmp(tmp == i) = [];
            tmp = [tmp(2:end) zeros(1, max_k - length(tmp) + 1)];
            ni(i,:) = tmp;
        end

    % Perform normal neighborhood selection
    else

        % Compute distances in batches
        n = size(X, 1);
        sum_X = sum(X .^ 2, 2);
        batch_size = round(2e7 ./ n);
        D = zeros(n, k);
        ni = zeros(n, k);
        for i=1:batch_size:n
            batch_ind = i:min(i + batch_size - 1, n);
            DD = bsxfun(@plus, sum_X', bsxfun(@plus, sum_X(batch_ind), ...
                                                   -2 * (X(batch_ind,:) * X')));
            [DD, ind] = sort(DD, 2, 'ascend');
            D(batch_ind,:) = sqrt(DD(:,2:k + 1));
            ni(batch_ind,:) = ind(:,2:k + 1);
        end
        D(D == 0) = 1e-9;
        D = sparse(repmat(1:n, [1 k])', ni(:), D(:), n, n);
    end

author: BehnoodRasti


text:
function varargout = drawGraph(varargin)
%DRAWGRAPH Draw a graph, given as a set of vertices and edges
%
%   drawGraph(NODES, EDGES)
%   Draws a graph specified by a set of nodes (array N-by-2 or N-by-3,
%   corresponding to coordinate of each node), and a set of edges (an array
%   NE-by-2, containing for each edge the first and the second node).
%   Default drawing is a red circle for nodes and a blue line for edges.
%
%   drawGraph(NODES, EDGES, FACES)
%   Also draws faces of the graph as patches.
%
%   drawGraph(GRAPH)
%   Passes argument in a srtucture with at least 2 fields named 'nodes' and
%   'edges', and possibly one field 'faces', corresponding to previously
%   described parameters.
%   GRAPH can also be a cell array, whose first element is node array,
%   second element is edges array, and third element, if present, is faces
%   array.
%
%
%   drawGraph(..., SNODES)
%   drawGraph(..., SNODES, SEDGES)
%   drawGraph(..., SNODES, SEDGES, SFACES)
%   Specifies the draw mode for each element, as in the classical 'plot'
%   function. To not display some elements, uses 'none'.
%
%
%   H = drawGraph(...)
%   Returns handle to the set of edges.
%
%   [HN, HE] = drawGraph(...)
%   Returns handle to the set of nodes and to the set of edges.
%
%   [HN, HE, HF] = drawGraph(...)
%   Also returns handle to the set of faces.
%
%   See Also
%   graphs, drawGraphEdges, fillGraphFaces, clipGraph, clipGraphPolygon
%

%   -----
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 17/07/2003.
%

%   HISTORY
%   10/02/2004 : documentation
%   06/04/2004 : change name
%   09/07/2004 : add faces
%   05/08/2004 : correct bug when drawing 2D graph
%   06/08/2004 : small bug for drawing poins (length instead of size(.,1))
%   09/08/2004 : rewrite code (separate 2D and 3D, use of plot instead of
%       line, manage faces if present in 2D and 3D, ...), add style
%       management, various input types, and documentation
%   22/09/2004 : correct bug in drawing faces
%   11/11/2005 : forgot a loop index for faces stored as cells
%   22/05/2009 add more drawing options

%% initialisations

% uses empty arrays by default for edges and faces
e = [];
f = [];

% default styles for nodes, edges, and faces

% nodes are drawn as red circles
sn = {'linestyle', 'none', 'color', 'r', 'marker', 'o'};

% edges are drawn as blue lines
se = {'linestyle', '-', 'color', 'b'};

% faces are cyan, their edges are not drawn
sf = {'EdgeColor', 'none', 'Facecolor', 'c'};

%% Process input arguments

% case of a call without arguments
if nargin == 0
    help drawGraph;
    return;
end

% ---------------------------------------------------------------
% First extract the graph structure

var = varargin{1};
if iscell(var)
    % graph is stored as a cell array: first cell is nodes, second one is
    % edges, and third is faces
    n = var{1};
    if length(var)>1
        e = var{2};
    end
    if length(var)>2
        f = var{3};
    end
    varargin(1) = [];
elseif isstruct(var)
    % graph is stored as a structure, with fields 'nodes', 'edges', and
    % eventually 'faces'.
    n = var.nodes;
    e = var.edges;
    if isfield(var, 'faces')
        f = var.faces;
    end
    varargin(1) = [];
else
    % graph is stored as set of variables: nodes, edges, and eventually
    % faces
    n = varargin{1};
    e = varargin{2};
    varargin(1:2) = [];

    if ~isempty(varargin)
        var = varargin{1};
        if isnumeric(var)
            % faces are stored in a numeric array of indices
            f = var;
            varargin(1) = [];
        elseif iscell(var)
            if ~ischar(var{1})
                % faces are stored in a cell array, each cell containing a
                % row vector of indices
                f = var;
                varargin(1) = [];
            end
        end
    end
end

% extract drawing style

if ~isempty(varargin)
    sn = concatArguments(sn, varargin{1});
end

if length(varargin)>1
    se = concatArguments(se, varargin{2});
end

if length(varargin)>2
    sf = concatArguments(sf, varargin{3});
end

%% Main drawing processing

hold on;

if size(n, 2) == 2
    % Draw a 2 dimensional graph ----------------------

    % Draw faces of the graph ------------
    if ~strcmp(sf{1}, 'none') && ~isempty(f)
        if iscell(f)
            % each face is contained in a cell.
            hf = zeros(size(f));
            for fi = 1:length(f)
                hf(fi) = patch('Faces', f{fi}, 'Vertices', n, sf{:});
            end
        else
            % process faces as an Nf*N array. Nf is the number of faces,
            % and all faces have the same number of vertices (nodes).
            hf = patch('Faces', f, 'Vertices', n, sf{:});
        end
    end

    % Draw 2D Edges ----------------------
    if ~strcmp(se{1}, 'none') && size(e, 1) > 0
        he = plot([n(e(:,1),1) n(e(:,2),1)]', [n(e(:,1),2) n(e(:,2),2)]', se{:});
    end

    % Draw 2D nodes ----------------------
    if ~strcmp(sn{1}, 'none')
        hn = plot(n(:,1), n(:,2), sn{:});
    end

elseif size(n, 2)==3
    % Draw a 3 dimensional graph ----------------------

    % use a zbuffer to avoid display pbms.
    set(gcf, 'renderer', 'zbuffer');

    % Draw 3D Faces ----------------------
    if ~strcmp(sf{1}, 'none')
        if iscell(f)
            % each face is contained in a cell.
            hf = zeros(size(f));
            for fi = 1:length(f)
                hf(fi) = patch('Faces', f{fi}, 'Vertices', n, sf{:});
            end
        else
            % process faces as an Nf*N array. Nf i the number of faces,
            % and all faces have the same number of vertices (nodes).
            hf = patch('Faces', f, 'Vertices', n, sf{:});
        end
    end

    % Draw 3D edges ----------------------
    if ~strcmp(se{1}, 'none') && size(e, 1) > 0
        he = line(...
            [n(e(:,1),1) n(e(:,2),1)]', ...
            [n(e(:,1),2) n(e(:,2),2)]', ...
            [n(e(:,1),3) n(e(:,2),3)]', ...
            se{:});
    end

    % Draw 3D nodes ----------------------
    if ~strcmp(sn{1}, 'none')
        hn = plot3(n(:,1), n(:,2), n(:,3), sn{:});
    end

end

%% Format output arguments

% return handle to edges
if nargout==1
    varargout{1} = he;
end

% return handle to nodes and edges
if nargout==2
    varargout{1} = hn;
    varargout{2} = he;
end

% return handle to nodes, edges and faces
if nargout==3
    varargout{1} = hn;
    varargout{2} = he;
    varargout{3} = hf;
end

end

function res = concatArguments(in1, in2)
% in1 is a cell array already initialized
% in2 is an argument that can be:
%   - empty
%   - the string 'none'
%   - another cell array

if isempty(in2)
    res = in1;
    return;
end

if ischar(in2)
    if strcmp('none', in2)
        res = {'none'};
        return;
    end
end

if iscell(in1)
    res = [in1(:)' in2(:)'];
else
    res = [{in1} in2(:)];
end

end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function varargout=subsref(obj,subs)
% SUBSREF Subscripted reference
% An overloaded function...
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: subsref.m 7439 2018-10-10 17:19:31Z john $

if isempty(subs), return; end

switch subs(1).type
    case '{}'
        error('Cell contents reference from a non-cell array object.');
    case '.'
        varargout = access_fields(obj,subs);
        return;
end

if numel(subs)~=1, error('Expression too complicated'); end

dim  = [size(obj) ones(1,16)];
nd   = find(dim>1,1,'last')-1;
sobj = struct(obj);

if ~numel(subs.subs)
    [subs.subs{1:nd+1}] = deal(':');
elseif length(subs.subs) < nd
    l   = length(subs.subs);
    dim = [dim(1:(l-1)) prod(dim(l:end))];
    if numel(sobj) ~= 1
        error('Can only reshape simple file_array objects.');
    else
        if numel(sobj.scl_slope)>1 || numel(sobj.scl_inter)>1
            error('Can not reshape file_array objects with multiple slopes and intercepts.');
        end
        sobj.dim = dim;
    end
end

di   = ones(16,1);
args = cell(1,length(subs.subs));
for i=1:length(subs.subs)
    if ischar(subs.subs{i})
        if ~strcmp(subs.subs{i},':'), error('This shouldn''t happen....'); end
        if length(subs.subs) == 1
            args{i} = 1:prod(dim);
            k = 0;
            for j=1:length(sobj)
                sobj(j).dim = [prod(sobj(j).dim) 1];
                sobj(j).pos = [k+1 1];
                k           = k + sobj(j).dim(1);
            end
        else
            args{i} = 1:dim(i);
        end
    else
        args{i} = subs.subs{i};
    end
    di(i) = length(args{i});
end

if length(sobj)==1
    t = subfun(sobj,args{:});
else
   %dt  = datatypes;
   %dt  = dt([dt.code]==sobj(1).dtype); % assuming identical datatypes
   %t   = zeros(di',func2str(dt.conv));
    t   = zeros(di','double');
    for j=1:length(sobj)
        ps = [sobj(j).pos ones(1,length(args))];
        dm = [sobj(j).dim ones(1,length(args))];
        for i=1:length(args)
            msk      = find(args{i}>=ps(i) & args{i}<(ps(i)+dm(i)));
            args2{i} = msk;
            args3{i} = double(args{i}(msk))-ps(i)+1;
        end

        t  = subsasgn(t,struct('type','()','subs',{args2}),subfun(sobj(j),args3{:}));
    end
end
varargout = {t};

%==========================================================================
% function t = subfun(sobj,varargin)
%==========================================================================
function t = subfun(sobj,varargin)

%sobj.dim = [sobj.dim ones(1,16)];
try
    args = cell(size(varargin));
    for i=1:length(varargin)
        args{i} = int64(varargin{i});
    end
    t = file2mat(sobj,args{:});
catch
    t = multifile2mat(sobj,varargin{:});
end
if ~isempty(sobj.scl_slope) || ~isempty(sobj.scl_inter)
    slope = 1;
    inter = 0;
    if ~isempty(sobj.scl_slope), slope = sobj.scl_slope; end
    if ~isempty(sobj.scl_inter), inter = sobj.scl_inter; end
    if numel(slope)>1
        slope = resize_scales(slope,sobj.dim,varargin);
        t     = double(t).*slope;
    else
        t     = double(t)*slope;
    end
    if numel(inter)>1
        inter = resize_scales(inter,sobj.dim,varargin);
    end
    t = t + inter;
end

%==========================================================================
% function c = access_fields(obj,subs)
%==========================================================================
function c = access_fields(obj,subs)

sobj = struct(obj);
c    = cell(1,numel(sobj));
for i=1:numel(sobj)
    %obj = class(sobj(i),'file_array');
    obj = sobj(i);
    switch(subs(1).subs)
        case 'fname',      t = fname(obj);
        case 'dtype',      t = dtype(obj);
        case 'offset',     t = offset(obj);
        case 'dim',        t = dim(obj);
        case 'scl_slope',  t = scl_slope(obj);
        case 'scl_inter',  t = scl_inter(obj);
        case 'permission', t = permission(obj);
        otherwise
            error(['Reference to non-existent field "' subs(1).subs '".']);
    end
    if numel(subs)>1
        t = subsref(t,subs(2:end));
    end
    c{i} = t;
end

%==========================================================================
% function val = multifile2mat(sobj,varargin)
%==========================================================================
function val = multifile2mat(sobj,varargin)

% Convert subscripts into linear index
[indx2{1:length(varargin)}] = ndgrid(varargin{:},1);
if numel(sobj.dim) == 1, sobj.dim = [sobj.dim 1]; end
ind = sub2ind(sobj.dim,indx2{:});

% Work out the partition
dt  = datatypes;
dt  = dt([dt.code]==sobj.dtype);
sz  = dt.size;
try
    mem = spm('Memory'); % in bytes, has to be a multiple of 16 (max([dt.size]))
catch
    mem = 200 * 1024 * 1024;
end
s   = ceil(prod(sobj.dim) * sz / mem);

% Assign indices to partitions
[x,y] = ind2sub([mem/sz s],ind(:));
c     = histc(y,1:s);
cc    = [0 reshape(cumsum(c),1,[])];

% Read data in relevant partitions
obj = sobj;
val = zeros(length(x),1,func2str(dt.conv));
for i=reshape(find(c),1,[])
    obj.offset = sobj.offset + mem*(i-1);
    obj.dim = [1 min(mem/sz, prod(sobj.dim)-(i-1)*mem/sz)];
    val(cc(i)+1:cc(i+1)) = file2mat(obj,int64(1),int64(x(y==i)));
end
r   = cellfun('length',varargin);
if numel(r) == 1, r = [r 1]; end
val = reshape(val,r);

author: spm


text:
function combo_test33 ( )

%*****************************************************************************80
%
%% COMBO_TEST33 tests RGF_G_TABLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 July 2011
%
%  Author:
%
%    John Burkardt
%
  m = 6;

  fprintf ( 1, ' \n' );
  fprintf ( 1, 'COMBO_TEST33\n' );
  fprintf ( 1, '  RGF_G_TABLE tabulates generalized restricted\n' );
  fprintf ( 1, '  growth functions.\n' );
  fprintf ( 1, ' \n' );

  offset = 1;

  d = rgf_g_table ( m );

  for i = 0 : m
    for j = 0 : m - i
      fprintf ( 1, '%6d', d(i+offset,j+offset) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function [ xmax, ixmax ] = i4row_max ( m, n, x )

%*****************************************************************************80
%
%% I4ROW_MAX returns the maximums of rows of an I4ROW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, integer X(M,N), the array to be examined.
%
%    Output, integer XMAX(M), the maximums of the rows of X.
%
%    Output, integer IXMAX(M); IXMAX(I) is the column of X in which
%    the maximum for row I occurs.
%
  for i = 1 : m

    ixmax(i) = 1;
    xmax(i) = x(i,1);
    for j = 2 : n
      if ( xmax(i) < x(i,j) )
        ixmax(i) = j;
        xmax(i) = x(i,j);
      end
    end

  end

  return
end


text:
function [result]=sv_timeperdiff(mCatalog, fTimePeriod)
% function [result]=sv_timeperdiff(mCatalog, fTimePeriod)
% --------------------------------------------------------------------------------------
% Function to calculate absolute difference of number of events between
% two time periods of an earthquake catalog
% Author: J. Woessner
% woessner@seismo.ifg.ethz.ch
% last update: 08.07.02
%
% Incoming variables:
% mCatalog           : current earthquake catalog
% params.fTimePeriod : Time period in days
%
% Outgoing variable:
% result.dNdiffsumVal      : total difference of number of events in the two time periods
% result.dNdiffsumYearVal  : total difference of number of events in the two time periods normalized to a year
% result.dNdiffsumMonthVal : total difference of number of events in the two time periods normalized to a month
% result.dNdiff            : Difference of seismicity in 0.1 magnitude bins

% Init variable
result=[];

% Create the catalogs for two time periods
fStartTime = min(mCatalog(:,3));
fEndTime = max(mCatalog(:,3));

[result.mFirstCatalog, result.mSecondCatalog, result.fFirstPeriodExact, result.fSecondPeriodExact, result.fFirstPeriod,...
        result.fSecondPeriod] = ex_SplitCatalog(mCatalog, fSplitTime, 0, 100, 0, 100);

[dEv_val dMags dEv_valsum dEv_valsum_rev,  dMags_rev] =fcumulsum(result.mFirstCatalog);
[dEv_val2 dMags2 dEv_valsum2 dEv_valsum_rev2,  dMags_rev2] =fcumulsum(result.mSecondCatalog);

result.dNdiff = dEv_val2-dEv_val;
result.dNdiffYear=result.dNdiff/365;
result.dNdiffMonth = result.dNdiffYear/12;

result.dNdiffsum = cumsum(result.dNdiff');
result.dNdiffsumVal = result.dNdiffsum(length(result.dNdiffsum));
result.dNdiffsumYearVal = result.dNdiffsumVal/365;
result.dNdiffsumMonthVal = result.dNdiffsumYearVal/12;

result.dMags = dMags;
result.dMags = dMags2;
return

author: CelsoReyes


text:
function [y] = cellcolremove(x, cols)

% [Y] = CELLCOLREMOVE(X, COLS) outputs cell-array Y with the same dimensionality as X
% Each cell in Y has the indexed columns COLS removed from the original corresponding cell in X
%
% X (and Y) should be linear cell-array(s) of matrices for which the number of rows
% should be the same for all cells

nx = size(x);
if ~iscell(x) || length(nx)>2 || all(nx>1),
  error('incorrect input for cellrows');
end

y = cellfun(@colr, x, repmat(mat2cell(cols(:),length(cols),1),nx), 'UniformOutput', 0);

function [y] = colr(x, cols)

if ischar(cols)
  if strcmp(cols(1),'r')
    % remove n columns on the right
    cols = str2double(cols(2:end));
    cols = size(x,2)-((cols-1):-1:0);
  elseif strcmp(cols(1),'l')
    % remove n columns on the left
    cols = str2double(cols(2:end));
    cols = 1:cols;
  else
    % numeric
  end
end

x(:,cols) = [];
y         = x;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [nodes2, edges2] = grRemoveEdges(nodes, edges, rmEdges)
%GRREMOVEEDGES Remove several edges from a graph.
%
%   [NODES2 EDGES2] = grRemoveEdges(NODES, EDGES, EDGES2REMOVE)
%   Remove some edges in the edges list, and return the modified graph.
%   The NODES array is not modified.

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-08-13
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

rmEdges = sort(rmEdges);

% do not change the node list
nodes2 = nodes;

% number of edges
N   = size(edges, 1);
NR  = length(rmEdges);
N2  = N - NR;

% allocate memory for new  edges list
edges2 = zeros(N2, 2);

% case of no edge to remove
if NR == 0
    nodes2 = nodes;
    edges2 = edges;
    return;
end

% process the first edge
edges2(1:rmEdges(1)-1,:) = edges(1:rmEdges(1)-1,:);

% process the classical edges
for i = 2:NR
    %if rmEdges(i)-i < rmEdges(i-1)-i+2
    %    continue;
    %end
    edges2(rmEdges(i-1)-i+2:rmEdges(i)-i, :) = edges(rmEdges(i-1)+1:rmEdges(i)-1, :);
end

% process the last edge
edges2(rmEdges(NR)-NR+1:N2, :) = edges(rmEdges(NR)+1:N, :);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%% brute_force_tune
% Code to test the performance of various tuning parameters
% Works sorta like RANSAC I guess?
% Adam Werries 2016, see Apache 2.0 license.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% pos_sd_min = linspace(0,2,2000);
% num_items = length(pos_sd_min);
% rms_error_filter = Inf*ones(1,num_items);
% max_error_filter = Inf*ones(1,num_items);
% parfor i = 1:num_items
%     fprintf('Iteration: %d, pos_sd_min: %08.7f\n', i, pos_sd_min(i));
%     temp_conf = LC_KF_config;
%     temp_conf.pos_sd_min = pos_sd_min(i);
%     [out_profile,out_IMU_bias_est,out_KF_SD, out_R_matrix, out_Q_matrix, corrections] = ...
%         Loosely_coupled_INS_GNSS(init_cond, filter_time, epoch, lla, gps, imu, temp_conf, est_IMU_bias);
%     xyz = out_profile(:,2:4);
%     if ~any(any(isnan(xyz))) && ~any(any(isinf(xyz)))
%         llh = ecef2lla(xyz);
%         [x,y] = deg2utm(llh(:,1),llh(:,2));
%         x = x-min_x;
%         y = y-min_y;
%         distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;
%         rms_error_filter(i) = rms(distance);
%         max_error_filter(i) = max(distance);
%     end
% end
%
% [minmax, i] = min(max_error_filter);
% fprintf('\nBest max: %08.7f, rms is %08.7f\n', minmax, rms_error_filter(i));
% fprintf('Best iteration for max: %d, pos_sd_min: %08.7f\n', i, pos_sd_min(i));
% [minrms, i] = min(rms_error_filter);
% fprintf('Best rms: %08.7f, max is %08.7f\n', minrms, max_error_filter(i));
% fprintf('Best iteration for rms: %d, pos_sd_min: %08.7f\n', i, pos_sd_min(i));
% figurec;
% subplot(211);
% plot(pos_sd_min, rms_error_filter);
% xlabel('pos_sd_min');
% ylabel('max error');
%
% subplot(212);
% plot(pos_sd_min, max_error_filter);
% xlabel('pos_sd_min');
% ylabel('max error');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pos_sd_max = linspace(0,200,2000);
% num_items = length(pos_sd_max);
% rms_error_filter = Inf*ones(1,num_items);
% max_error_filter = Inf*ones(1,num_items);
% parfor i = 1:num_items
%     fprintf('Iteration: %d, pos_sd_max: %08.7f\n', i, pos_sd_max(i));
%     temp_conf = LC_KF_config;
%     temp_conf.pos_sd_max = pos_sd_max(i);
%     [out_profile,out_IMU_bias_est,out_KF_SD] = Loosely_coupled_INS_GNSS(init_cond, filter_time, epoch, lla, gps, imu, temp_conf, est_IMU_bias);
%     xyz = out_profile(:,2:4);
%     if ~any(any(isnan(xyz))) && ~any(any(isinf(xyz)))
%         llh = ecef2lla(xyz);
%         [x,y] = deg2utm(llh(:,1),llh(:,2));
%         x = x-min_x;
%         y = y-min_y;
%         distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;
%         rms_error_filter(i) = rms(distance);
%         max_error_filter(i) = max(distance);
%     end
% end
%
% [minmax, i] = min(max_error_filter);
% fprintf('\nBest max: %08.7f, rms is %08.7f\n', minmax, rms_error_filter(i));
% fprintf('Best iteration for max: %d, pos_sd_max: %08.7f\n', i, pos_sd_max(i));
% [minrms, i] = min(rms_error_filter);
% fprintf('Best rms: %08.7f, max is %08.7f\n', minrms, max_error_filter(i));
% fprintf('Best iteration for rms: %d, pos_sd_max: %08.7f\n', i, pos_sd_max(i));
% figurec;
% plot(pos_sd_max, max_error_filter); hold on;
% plot(pos_sd_max, rms_error_filter);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% vel_sd_min = linspace(0,30,2000);
% num_items = length(vel_sd_min);
% rms_error_filter = Inf*ones(1,num_items);
% max_error_filter = Inf*ones(1,num_items);
% parfor i = 1:num_items
%     fprintf('Iteration: %d, vel_sd_min: %08.7f\n', i, vel_sd_min(i));
%     temp_conf = LC_KF_config;
%     temp_conf.vel_sd_min = vel_sd_min(i);
%     [out_profile,out_IMU_bias_est,out_KF_SD] = Loosely_coupled_INS_GNSS(init_cond, filter_time, epoch, lla, gps, imu, temp_conf, est_IMU_bias);
%     xyz = out_profile(:,2:4);
%     if ~any(any(isnan(xyz))) && ~any(any(isinf(xyz)))
%         llh = ecef2lla(xyz);
%         [x,y] = deg2utm(llh(:,1),llh(:,2));
%         x = x-min_x;
%         y = y-min_y;
% %         h = -llh(:,3);
%         distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;
%         rms_error_filter(i) = rms(distance);
%         max_error_filter(i) = max(distance);
%     end
% end
%
% [minmax, i] = min(max_error_filter);
% fprintf('\nBest max: %08.7f, rms is %08.7f\n', minmax, rms_error_filter(i));
% fprintf('Best iteration for max: %d, vel_sd_min: %08.7f\n', i, vel_sd_min(i));
% [minrms, i] = min(rms_error_filter);
% fprintf('Best rms: %08.7f, max is %08.7f\n', minrms, max_error_filter(i));
% fprintf('Best iteration for rms: %d, vel_sd_min: %08.7f\n', i, vel_sd_min(i));
% figure;
% plot(vel_sd_min, max_error_filter); hold on;
% plot(vel_sd_min, rms_error_filter);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vel_sd_max = linspace(1,50,200);
num_items = length(vel_sd_max);
rms_error_filter = Inf*ones(1,num_items);
max_error_filter = Inf*ones(1,num_items);
parfor i = 1:num_items
    fprintf('Iteration: %d, vel_sd_max: %08.7f\n', i, vel_sd_max(i));
    temp_conf = LC_KF_config;
    temp_conf.vel_sd_max = vel_sd_max(i);
    [out_profile,out_IMU_bias_est,out_KF_SD, out_R_matrix, out_Q_matrix, corrections] = ...
        Loosely_coupled_INS_GNSS(init_cond, filter_time, epoch, lla, gps, imu, temp_conf, est_IMU_bias);
    xyz = out_profile(:,2:4);
    if ~any(any(isnan(xyz))) && ~any(any(isinf(xyz)))
        llh = ecef2lla(xyz);
        [x,y] = deg2utm(llh(:,1),llh(:,2));
        x = x-min_x;
        y = y-min_y;
%         h = -llh(:,3);
%         distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2 + (ground_truth_full(:,3)-h).^2).^0.5;
        distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;
        rms_error_filter(i) = rms(distance);
        max_error_filter(i) = max(distance);
    end
end

[minmax, i] = min(max_error_filter);
fprintf('\nBest max: %08.7f, rms is %08.7f\n', minmax, rms_error_filter(i));
fprintf('Best iteration for max: %d, vel_sd_max: %08.7f\n', i, vel_sd_max(i));
[minrms, i] = min(rms_error_filter);
fprintf('Best rms: %08.7f, max is %08.7f\n', minrms, max_error_filter(i));
fprintf('Best iteration for rms: %d, vel_sd_max: %08.7f\n', i, vel_sd_max(i));
figurec;
plot(vel_sd_max, max_error_filter); hold on;
plot(vel_sd_max, rms_error_filter);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gps_pos_stddev = linspace(.001,10,2000);
% num_items = length(gps_pos_stddev);
% rms_error_filter = Inf*ones(1,num_items);
% max_error_filter = Inf*ones(1,num_items);
% parfor i = 1:num_items
%     fprintf('Iteration: %d, gps_pos_stddev: %08.7f\n', i, gps_pos_stddev(i));
%     temp_conf = LC_KF_config;
%     temp_conf.gps_pos_stddev = gps_pos_stddev(i);
%     [out_profile,out_IMU_bias_est,out_KF_SD, out_R_matrix, out_Q_matrix, corrections] = ...
%         Loosely_coupled_INS_GNSS(init_cond, filter_time, epoch, lla, gps, imu, temp_conf, est_IMU_bias);
%     xyz = out_profile(:,2:4);
%     if ~any(any(isnan(xyz))) && ~any(any(isinf(xyz)))
%         llh = ecef2lla(xyz);
%         [x,y] = deg2utm(llh(:,1),llh(:,2));
%         x = x-min_x;
%         y = y-min_y;
% %         h = -llh(:,3);
% %         distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2 + (ground_truth_full(:,3)-h).^2).^0.5;
%         distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;
%         rms_error_filter(i) = rms(distance);
%         max_error_filter(i) = max(distance);
%     end
% end
%
% [minmax, i] = min(max_error_filter);
% fprintf('\nBest max: %08.7f, rms is %08.7f\n', minmax, rms_error_filter(i));
% fprintf('Best iteration for max: %d, gps_pos_stddev: %08.7f\n', i, gps_pos_stddev(i));
% [minrms, i] = min(rms_error_filter);
% fprintf('Best rms: %08.7f, max is %08.7f\n', minrms, max_error_filter(i));
% fprintf('Best iteration for rms: %d, gps_pos_stddev: %08.7f\n', i, gps_pos_stddev(i));
% figurec;
% plot(gps_pos_stddev, max_error_filter); hold on;
% plot(gps_pos_stddev, rms_error_filter);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gps_vel_stddev = linspace(.01,100,2000);
% num_items = length(gps_vel_stddev);
% rms_error_filter = Inf*ones(1,num_items);
% max_error_filter = Inf*ones(1,num_items);
% parfor i = 1:num_items
%     fprintf('Iteration: %d, vel_sd_max: %08.7f\n', i, gps_vel_stddev(i));
%     temp_conf = LC_KF_config;
%     temp_conf.gps_vel_stddev = gps_vel_stddev(i);
%     [out_profile,out_IMU_bias_est,out_KF_SD, out_R_matrix, out_Q_matrix, corrections] = ...
%         Loosely_coupled_INS_GNSS(init_cond, filter_time, epoch, lla, gps, imu, temp_conf, est_IMU_bias);
%     xyz = out_profile(:,2:4);
%     if ~any(any(isnan(xyz))) && ~any(any(isinf(xyz)))
%         llh = ecef2lla(xyz);
%         [x,y] = deg2utm(llh(:,1),llh(:,2));
%         x = x-min_x;
%         y = y-min_y;
% %         h = -llh(:,3);
% %         distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2 + (ground_truth_full(:,3)-h).^2).^0.5;
%         distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;        rms_error_filter(i) = rms(distance);
%         max_error_filter(i) = max(distance);
%     end
% end
%
% [minmax, i] = min(max_error_filter);
% fprintf('\nBest max: %08.7f, rms is %08.7f\n', minmax, rms_error_filter(i));
% fprintf('Best iteration for max: %d, gps_vel_stddev: %08.7f\n', i, gps_vel_stddev(i));
% [minrms, i] = min(rms_error_filter);
% fprintf('Best rms: %08.7f, max is %08.7f\n', minrms, max_error_filter(i));
% fprintf('Best iteration for rms: %d, gps_vel_stddev: %08.7f\n', i, gps_vel_stddev(i));
% figurec;
% plot(gps_vel_stddev, max_error_filter); hold on;
% plot(gps_vel_stddev, rms_error_filter);

load handel
sound(y,Fs)

author: awerries

repo: kalman-localization


text:
classdef Optimizer_MMA < Optimizer

    properties (GetAccess = public, SetAccess = protected)
        type = 'MMA'
    end

    properties (Access = private)
        kkttol
        maxoutit
        x
        xold1
        xold2
        xmin
        xmax
        low
        outit = 0;
        outeriter = 0;
        upp
        m
        c
        d
        a0
        a
        n
        f0val
        df0dx
        fval
        dfdx
        upperBound
        lowerBound
        hasFinished
        incrementalScheme
        hasConverged
        historicalVariables
       % targetParameters
        KKTnorm
    end

    methods (Access = public)

        function obj = Optimizer_MMA(cParams)
            obj.initOptimizer(cParams);
            obj.init(cParams);
            obj.outputFunction.monitoring.create(cParams);
%             obj.upperBound = cParams.uncOptimizerSettings.ub;
%             obj.lowerBound = cParams.uncOptimizerSettings.lb;
            obj.maxoutit = 1e4;
        end

       function solveProblem(obj)
            % obj.cost.computeFunctionAndGradient();
            % obj.constraint.computeFunctionAndGradient();
%             obj.printOptimizerVariable();
            obj.hasFinished = false;
            obj.printOptimizerVariable();
            while ~obj.hasFinished
                obj.update();
                obj.updateIterInfo();
                obj.updateMonitoring();
                obj.printOptimizerVariable();
            end
%             obj.printOptimizerVariable();
%             obj.printHistory();
            obj.hasConverged = 0;
%             obj.printHistoryFinalValues();
       end

        function update(obj)
            x = obj.designVariable.value;
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();
            obj.checkInitial(x);
            obj.outit = obj.outit+1;
            obj.outeriter = obj.outeriter+1;
            %%%% The MMA subproblem is solved at the point xval:
            [xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,obj.low,obj.upp] = ...
                obj.mmasub(obj.m,obj.n,obj.outeriter,x,obj.xmin,obj.xmax,obj.xold1,obj.xold2, ...
                obj.f0val,obj.df0dx,obj.fval,obj.dfdx,obj.low,obj.upp,obj.a0,obj.a,obj.c,obj.d);
            %%%% Some vectors are updated:
            obj.xold2 = obj.xold1;
            obj.xold1 = x;
            x = xmma;
            %%%% The user should now calculate function values and gradients
            %%%% of the objective- and constraint functions at xval.
            %%%% The results should be put in f0val, df0dx, fval and dfdx.
            obj.designVariable.update(x);
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();

            [obj.f0val,obj.df0dx,obj.fval,obj.dfdx] = obj.funmma();
            %%%% The residual vector of the KKT conditions is calculated:
            [~,kktnorm] = obj.kktcheck(obj.m,obj.n,xmma,ymma,zmma,lam,xsi,eta,mu,zet,s, ...
                obj.xmin,obj.xmax,obj.df0dx,obj.fval,obj.dfdx,obj.a0,obj.a,obj.c,obj.d);

            obj.historicalVariables.kktnorm = kktnorm;
            obj.dualVariable.value = lam;
            obj.updateConvergenceStatus();
            obj.KKTnorm     = kktnorm;
        end

    end

    methods (Access = private)

        function updateMonitoring(obj)
            s.hasFinished          = obj.hasFinished;
            s.nIter                = obj.nIter;
            s.KKTnorm              = obj.KKTnorm;
            s.outitFrac            = obj.outit/obj.maxoutit;
            obj.outputFunction.monitoring.compute(s);
        end

        function init(obj,cParams)
            obj.outputFunction         = cParams.outputFunction.monitoring;
            obj.upperBound             = cParams.uncOptimizerSettings.ub;
            obj.lowerBound             = cParams.uncOptimizerSettings.lb;
            obj.incrementalScheme      = cParams.incrementalScheme;
            obj.hasConverged           = false;
            obj.constraintCase         = cParams.constraintCase;
            obj.targetParameters       = cParams.targetParameters;
        end

        function [f,df,c,dc] = funmma(obj)
            f  = obj.cost.value;
            df = obj.cost.gradient;
            c  = obj.constraint.value;
            dc = obj.constraint.gradient;
            dc = dc';

            [c,dc] = obj.checkConstraintCase(c,dc);

            %% Re-scale constraints
            % In many applications, the constraints are on the form yi(x) < =  ymaxi
            % The user should then preferably scale the constraints in such a way that 1 < =  ymaxi < =  100 for each i
            % (and not ymaxi = 10^10 for example).

            % kconstr = 100;
            kconstr = 1;
            cconstr = 0;
            c = kconstr*c;
            c(c > 0) = c(c > 0) + cconstr;
            c(c < 0) = c(c < 0) - cconstr;
            % dc = kconstr*dc;

            %% Re-scale objective function
            % The objective function f(x) should preferably be scaled such that
            % 1 < =  f0(x) < =  100 for reasonable values on the variables.
            kfun = 1;
            cfun = 0;
            f = kfun*f + cfun;
            df = kfun*df;
        end

        function checkInitial(obj,x0)
            if isempty(obj.x)
                obj.x = x0;
                obj.xold1 = obj.x;
                obj.xold2 = obj.xold1;
                obj.xmin = obj.lowerBound*ones(length(x0),1);
                obj.xmax = obj.upperBound*ones(length(x0),1);
                % obj.low = obj.xmin;
                obj.low = zeros(length(x0),1);
                % obj.upp = obj.xmax;
                obj.upp = ones(length(x0),1);
                [obj.f0val,obj.df0dx,obj.fval,obj.dfdx] = obj.funmma();
                obj.m = length(obj.fval);
                obj.c = 1000*ones(obj.m,1);
                obj.d = 0*ones(obj.m,1);
                obj.a0 = 1;
                obj.a = 0*ones(obj.m,1);
                obj.n = length(obj.x);
            end
        end

        function [xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,low,upp] = ...
                mmasub(obj,m,n,iter,xval,xmin,xmax,xold1,xold2, ...
                ~,df0dx,fval,dfdx,low,upp,a0,a,c,d)
            epsimin = sqrt(m+n)*10^(-9); %10^(-7);
            raa0 = 1e-5;
            move = 1.0;
            albefa = 0.1;
            asyinit = 0.1 ; %0.2 % 0.5;
            asyincr = 1.1; % 1.2;
            asydecr = 0.65; % 0.7
            eeen = ones(n,1);
            eeem = ones(m,1);
            zeron = zeros(n,1);

            % Calculation of the asymptotes low and upp :
            if iter < 2.5
                low = xval - asyinit*(xmax-xmin);
                upp = xval + asyinit*(xmax-xmin);
            else
                zzz = (xval-xold1).*(xold1-xold2);
                factor = eeen;
                factor(find(zzz > 0)) = asyincr;
                factor(find(zzz < 0)) = asydecr;
                low = xval - factor.*(xold1 - low);
                upp = xval + factor.*(upp - xold1);
                lowmin = xval - 10*(xmax-xmin);
                lowmax = xval - 0.01*(xmax-xmin);
                uppmin = xval + 0.01*(xmax-xmin);
                uppmax = xval + 10*(xmax-xmin);
                low = max(low,lowmin);
                low = min(low,lowmax);
                upp = min(upp,uppmax);
                upp = max(upp,uppmin);
            end

            % Calculation of the bounds alfa and beta :

            zzz1 = low + albefa*(xval-low);
            zzz2 = xval - move*(xmax-xmin);
            zzz  = max(zzz1,zzz2);
            alfa = max(zzz,xmin);
            zzz1 = upp - albefa*(upp-xval);
            zzz2 = xval + move*(xmax-xmin);
            zzz  = min(zzz1,zzz2);
            beta = min(zzz,xmax);

            % Calculations of p0, q0, P, Q and b.

            xmami = xmax-xmin;
            xmamieps = 1e-5*eeen;
            xmami = max(xmami,xmamieps);
            xmamiinv = eeen./xmami;
            ux1 = upp-xval;
            ux2 = ux1.*ux1;
            xl1 = xval-low;
            xl2 = xl1.*xl1;
            uxinv = eeen./ux1;
            xlinv = eeen./xl1;
            %
            p0 = zeron;
            q0 = zeron;
            p0 = max(df0dx,0);
            q0 = max(-df0dx,0);
            %p0(find(df0dx > 0)) = df0dx(find(df0dx > 0));
            %q0(find(df0dx < 0)) = -df0dx(find(df0dx < 0));
            pq0 = 0.001*(p0 + q0) + raa0*xmamiinv;
            p0 = p0 + pq0;
            q0 = q0 + pq0;
            p0 = p0.*ux2;
            q0 = q0.*xl2;
            %
            dfdx = sparse(dfdx);
            P = max(dfdx,0);
            Q = max(-dfdx,0);
            %P(find(dfdx > 0)) = dfdx(find(dfdx > 0));
            %Q(find(dfdx < 0)) = -dfdx(find(dfdx < 0));
            PQ = 0.001*(P + Q) + raa0*(eeem*xmamiinv');
            %PQ = 0.001*(P + Q) + raa0*eeem*xmamiinv;
            P = P + PQ;
            Q = Q + PQ;
            P = P * spdiags(ux2,0,n,n);
            Q = Q * spdiags(xl2,0,n,n);
            b = P*uxinv + Q*xlinv - fval ;
            %
            %%% Solving the subproblem by a primal-dual Newton method
            [xmma,ymma,zmma,lam,xsi,eta,mu,zet,s] = ...
                obj.subsolv(m,n,epsimin,low,upp,alfa,beta,p0,q0,P,Q,a0,a,b,c,d);
        end

        function updateConvergenceStatus(obj)
            kktnorm = obj.historicalVariables.kktnorm;
            has_not_converged = kktnorm > obj.kkttol && obj.outit < obj.maxoutit;
            obj.hasConverged = ~has_not_converged;
        end

        function [c,dc] = checkConstraintCase(obj,c,dc)
            if strcmp(obj.constraintCase,'EQUALITY')
                c  = [c;-c];
                dc = [dc;-dc];
            end
        end

        function updateIterInfo(obj)
            obj.increaseIter();
            obj.updateStatus();
        end

        function increaseIter(obj)
            obj.nIter = obj.nIter + 1;
        end

        function updateStatus(obj)
            obj.hasFinished = obj.hasConverged || obj.hasExceededStepIterations();
        end

        function itHas = hasExceededStepIterations(obj)
            iStep = obj.incrementalScheme.iStep;
            nStep = obj.incrementalScheme.nSteps;
            itHas = obj.nIter >= obj.maxIter*(iStep/nStep);
        end

    end

    methods (Access = private, Static)

        function [xmma,ymma,zmma,lamma,xsimma,etamma,mumma,zetmma,smma] = ...
                subsolv(m,n,epsimin,low,upp,alfa,beta,p0,q0,P,Q,a0,a,b,c,d)
            %
            % This function subsolv solves the MMA subproblem:
            %
            % minimize   SUM[ p0j/(uppj-xj) + q0j/(xj-lowj) ] + a0*z +
            %          + SUM[ ci*yi + 0.5*di*(yi)^2 ],
            %
            % subject to SUM[ pij/(uppj-xj) + qij/(xj-lowj) ] - ai*z - yi < =  bi,
            %            alfaj < =   xj < =   betaj,  yi > =  0,  z > =  0.
            %
            % Input:  m, n, low, upp, alfa, beta, p0, q0, P, Q, a0, a, b, c, d.
            % Output: xmma,ymma,zmma, slack variables and Lagrange multiplers.
            %
            een = ones(n,1);
            eem = ones(m,1);
            epsi = 1;
            epsvecn = epsi*een;
            epsvecm = epsi*eem;
            x = 0.5*(alfa+beta);
            y = eem;
            z = 1;
            lam = eem;
            xsi = een./(x-alfa);
            xsi = max(xsi,een);
            eta = een./(beta-x);
            eta = max(eta,een);
            mu  = max(eem,0.5*c);
            zet = 1;
            s = eem;
            itera = 0;
            while epsi > epsimin
                epsvecn = epsi*een;
                epsvecm = epsi*eem;
                ux1 = upp-x;
                xl1 = x-low;
                ux2 = ux1.*ux1;
                xl2 = xl1.*xl1;
                uxinv1 = een./ux1;
                xlinv1 = een./xl1;
                plam = p0 + P'*lam ;
                qlam = q0 + Q'*lam ;
                gvec = P*uxinv1 + Q*xlinv1;
                dpsidx = plam./ux2 - qlam./xl2 ;
                rex = dpsidx - xsi + eta;
                rey = c + d.*y - mu - lam;
                rez = a0 - zet - a'*lam;
                relam = gvec - a*z - y + s - b;
                rexsi = xsi.*(x-alfa) - epsvecn;
                reeta = eta.*(beta-x) - epsvecn;
                remu = mu.*y - epsvecm;
                rezet = zet*z - epsi;
                res = lam.*s - epsvecm;
                residu1 = [rex' rey' rez]';
                residu2 = [relam' rexsi' reeta' remu' rezet res']';
                residu = [residu1' residu2']';
                residunorm = sqrt(residu'*residu);
                residumax = max(abs(residu));
                ittt = 0;
                while residumax > 0.9*epsi & ittt < 200 % 100
                    ittt = ittt + 1;
                    itera = itera + 1;
                    ux1 = upp-x;
                    xl1 = x-low;
                    ux2 = ux1.*ux1;
                    xl2 = xl1.*xl1;
                    ux3 = ux1.*ux2;
                    xl3 = xl1.*xl2;
                    uxinv1 = een./ux1;
                    xlinv1 = een./xl1;
                    uxinv2 = een./ux2;
                    xlinv2 = een./xl2;
                    plam = p0 + P'*lam ;
                    qlam = q0 + Q'*lam ;
                    gvec = P*uxinv1 + Q*xlinv1;
                    GG = P*spdiags(uxinv2,0,n,n) - Q*spdiags(xlinv2,0,n,n);
                    dpsidx = plam./ux2 - qlam./xl2 ;
                    delx = dpsidx - epsvecn./(x-alfa) + epsvecn./(beta-x);
                    dely = c + d.*y - lam - epsvecm./y;
                    delz = a0 - a'*lam - epsi/z;
                    dellam = gvec - a*z - y - b + epsvecm./lam;
                    diagx = plam./ux3 + qlam./xl3;
                    diagx = 2*diagx + xsi./(x-alfa) + eta./(beta-x);
                    diagxinv = een./diagx;
                    diagy = d + mu./y;
                    diagyinv = eem./diagy;
                    diaglam = s./lam;
                    diaglamyi = diaglam+diagyinv;
                    if m < n
                        blam = dellam + dely./diagy - GG*(delx./diagx);
                        bb = [blam' delz]';
                        Alam = spdiags(diaglamyi,0,m,m) + GG*spdiags(diagxinv,0,n,n)*GG';
                        AA = [Alam     a
                            a'    -zet/z ];
                        solut = AA\bb;
                        dlam = solut(1:m);
                        dz = solut(m+1);
                        dx = -delx./diagx - (GG'*dlam)./diagx;
                    else
                        diaglamyiinv = eem./diaglamyi;
                        dellamyi = dellam + dely./diagy;
                        Axx = spdiags(diagx,0,n,n) + GG'*spdiags(diaglamyiinv,0,m,m)*GG;
                        azz = zet/z + a'*(a./diaglamyi);
                        axz = -GG'*(a./diaglamyi);
                        bx = delx + GG'*(dellamyi./diaglamyi);
                        bz  = delz - a'*(dellamyi./diaglamyi);
                        AA = [Axx   axz
                            axz'  azz ];
                        bb = [-bx' -bz]';
                        solut = AA\bb;
                        dx  = solut(1:n);
                        dz = solut(n+1);
                        dlam = (GG*dx)./diaglamyi - dz*(a./diaglamyi) + dellamyi./diaglamyi;
                    end
                    %
                    dy = -dely./diagy + dlam./diagy;
                    dxsi = -xsi + epsvecn./(x-alfa) - (xsi.*dx)./(x-alfa);
                    deta = -eta + epsvecn./(beta-x) + (eta.*dx)./(beta-x);
                    dmu  = -mu + epsvecm./y - (mu.*dy)./y;
                    dzet = -zet + epsi/z - zet*dz/z;
                    ds   = -s + epsvecm./lam - (s.*dlam)./lam;
                    xx  = [ y'  z  lam'  xsi'  eta'  mu'  zet  s']';
                    dxx = [dy' dz dlam' dxsi' deta' dmu' dzet ds']';
                    %
                    stepxx = -1.01*dxx./xx;
                    stmxx  = max(stepxx);
                    stepalfa = -1.01*dx./(x-alfa);
                    stmalfa = max(stepalfa);
                    stepbeta = 1.01*dx./(beta-x);
                    stmbeta = max(stepbeta);
                    stmalbe  = max(stmalfa,stmbeta);
                    stmalbexx = max(stmalbe,stmxx);
                    stminv = max(stmalbexx,1);
                    steg = 1/stminv;
                    %
                    xold   =   x;
                    yold   =   y;
                    zold   =   z;
                    lamold =  lam;
                    xsiold =  xsi;
                    etaold =  eta;
                    muold  =  mu;
                    zetold =  zet;
                    sold   =   s;
                    %
                    itto = 0;
                    resinew = 2*residunorm;
                    while resinew > residunorm & itto < 50
                        itto = itto+1;
                        x   =   xold + steg*dx;
                        y   =   yold + steg*dy;
                        z   =   zold + steg*dz;
                        lam = lamold + steg*dlam;
                        xsi = xsiold + steg*dxsi;
                        eta = etaold + steg*deta;
                        mu  = muold  + steg*dmu;
                        zet = zetold + steg*dzet;
                        s   =   sold + steg*ds;
                        ux1 = upp-x;
                        xl1 = x-low;
                        ux2 = ux1.*ux1;
                        xl2 = xl1.*xl1;
                        uxinv1 = een./ux1;
                        xlinv1 = een./xl1;
                        plam = p0 + P'*lam ;
                        qlam = q0 + Q'*lam ;
                        gvec = P*uxinv1 + Q*xlinv1;
                        dpsidx = plam./ux2 - qlam./xl2 ;
                        rex = dpsidx - xsi + eta;
                        rey = c + d.*y - mu - lam;
                        rez = a0 - zet - a'*lam;
                        relam = gvec - a*z - y + s - b;
                        rexsi = xsi.*(x-alfa) - epsvecn;
                        reeta = eta.*(beta-x) - epsvecn;
                        remu = mu.*y - epsvecm;
                        rezet = zet*z - epsi;
                        res = lam.*s - epsvecm;
                        residu1 = [rex' rey' rez]';
                        residu2 = [relam' rexsi' reeta' remu' rezet res']';
                        residu = [residu1' residu2']';
                        resinew = sqrt(residu'*residu);
                        steg = steg/2;
                    end
                    residunorm = resinew;
                    residumax = max(abs(residu));
                    steg = 2*steg;
                end
                if ittt > 198
                    epsi;
                    ittt;
                end
                epsi = 0.1*epsi;
            end
            xmma   =   x;
            ymma   =   y;
            zmma   =   z;
            lamma =  lam;
            xsimma =  xsi;
            etamma =  eta;
            mumma  =  mu;
            zetmma =  zet;
            smma   =   s;
        end

        function [residu,residunorm,residumax] = ...
                kktcheck(~,~,x,y,z,lam,xsi,eta,mu,zet,s, ...
                xmin,xmax,df0dx,fval,dfdx,a0,a,c,d)

            rex   = df0dx + dfdx'*lam - xsi + eta;
            rey   = c + d.*y - mu - lam;
            rez   = a0 - zet - a'*lam;
            relam = fval - a*z - y + s;
            rexsi = xsi.*(x-xmin);
            reeta = eta.*(xmax-x);
            remu  = mu.*y;
            rezet = zet*z;
            res   = lam.*s;
            %
            residu1 = [rex' rey' rez]';
            residu2 = [relam' rexsi' reeta' remu' rezet res']';
            residu = [residu1' residu2']';
            residunorm = sqrt(residu'*residu);
            residumax = max(abs(residu));
        end

    end

    methods

        function kkttol = get.kkttol(obj)
            kkttol = obj.targetParameters.optimality_tol;
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% set Matlab path and TOOLBOX_PATH environment variable
addpath(fullfile(pwd, 'matlab'));
setenv('TOOLBOX_PATH', pwd);

author: thomaskuestner


text:
function S = Diffmaps(S,varargin)

%% S = Diffmaps(S)
%% Returns component maps in S.DiffMap
%% S.DiffMap contains an m by n by j matrix with m and n being the x and
%% y pixel values, and j being the jth component map
%%
%% Stacks must be imported, aligned and converted to OD using ProcDir.m and
%% AlignOD.m
%%
%% COMPONENTS:
%% j=1 -- inorganic
%% j=2 -- potassium
%% j=3 -- sp2
%% j=4 -- CO3
%% j=5 -- carboxylic acid
%%
%% Example input
%% S =
%%
%%     eVenergy: [95x1 double]
%%       Xvalue: 4
%%       Yvalue: 4
%%       spectr: [67x67x95 double]  %% this is stack data!
%%     particle: '70202041'
%%
if isempty(varargin)
    spThresh=0.35;
    figsav=0;
elseif length(varargin)==1
    spThresh=varargin{1};
    figsav=0;
elseif length(varargin)==3
    spThresh=varargin{1};
    figsav=1;
    rootdir=varargin{2};
    sample=varargin{3};
end

stack=S.spectr;energy=S.eVenergy;particle=S.particle;

DiffMap=zeros(size(stack,1),size(stack,2),5);

peakregion1=find(energy > min(energy) & energy < 283);  %% Pre edge
if isempty(peakregion1)
    S.DiffMap=[];
    S.LabelMat=[];
    beep
    disp('This isnt the carbon edge');
    return
end

%% map inorganics
Preedge=mean(stack(:,:,peakregion1),3); %% compute average.
Preedge(Preedge==Inf)=0;
DetectThresh=3*mean(mean(std(stack(:,:,peakregion1),0,3)));
DiffMap(:,:,1)=ConstThresh(Preedge,DetectThresh);
DiffMap(:,:,1)=InorgMap1(S,DetectThresh,DiffMap(:,:,1));
Inorg=DiffMap(:,:,1);
DiffMap(DiffMap<.5 | DiffMap==NaN | DiffMap==Inf)=0;
DiffMap(:,:,1)=medfilt2(DiffMap(:,:,1));

postidx=find(energy > 315 & energy < 320);
postedge=mean(stack(:,:,postidx),3);

%% Potassium Map
% IaIdx=find(energy > 294 & energy < 295);  %% define first baseline point
IbIdx=find(energy > 302 & energy < 305);  %% second baseline point
DiffMap(:,:,2)=PeakDetect(S,[294,295],[302,305],[296.1,298.1],[298.6,301],[292,296],3);
DiffMap(:,:,2)=medfilt2(DiffMap(:,:,2));

%% 285 eV sp2 peak Minus post edge Map
ECpeakpos=find(energy > 284.75 & energy < 285.9);
Tmp=PeakDetect(S,[energy(1),energy(1)+1],[282,283],...
    [284.75,285.9],[284.75,285.9],[energy(1),283],6);
da=diff(stack(:,:,peakregion1),1,3);
de=diff(energy(peakregion1));
dade=zeros(size(stack(:,:,peakregion1)));
for i=1:length(de)
    dade(:,:,i)=da(:,:,i)./de(i);
end
avdade=mean(dade,3);
newstack=stack;
dep=diff(energy);
% calculate y intercept of preedge
for i=1:length(peakregion1)
    b(:,:,i)=stack(:,:,i)-avdade.*energy(i);
end
avb=mean(b,3);
% subtract linear preedge obtained above
for i=1:length(energy)
    newstack(:,:,i)=stack(:,:,i)-(avdade.*energy(i)+avb);
end

Tot = [280,310];
CCPeak = [284,287];
totidx=find(energy > Tot(1) & energy < Tot(2));
CCidx=find(energy > CCPeak(1) & energy < CCPeak(2));
Atot=trapz(energy(totidx),newstack(:,:,totidx),3);
ACC=trapz(energy(CCidx),newstack(:,:,CCidx),3);
sp2out=ACC./Atot.*8.6;
sp2out=medfilt2(sp2out.*Tmp);
sp2out(sp2out>1)=1;
sp2out(sp2out<0)=0;  %% sp2out is what is shown in the plot
Tmp=sp2out;

Tmp(Tmp<spThresh)=0; %% threshold map for output into partlabel.
DiffMap(:,:,3)=Tmp;
DiffMap(DiffMap(:,:,3)<0 | DiffMap(:,:,3)==NaN)=0;

% Post Edge Minus Pre Edge Map (total carbon)

% DiffMap(:,:,4)=medfilt2(negfilter(postedge(:,:)-Preedge));
% T.DiffMap=DiffMap;

%% CO3 Map

% DiffMap(:,:,4)=PeakDetect(S,[289.3,289.8],[291.1,292.5],...
%     [290.1,290.6],[290.1,290.6],[291.6,293.5]);
DiffMap(:,:,4)=PeakDetect(S,[289.3,289.8],[291.1,292.5],...
    [290.1,290.6],[290.1,290.6],[energy(1),283],6);
DiffMap(:,:,4)=medfilt2(DiffMap(:,:,4));
%% carboxylic map
CarbxIdx=find(energy > 288 & energy < 289);  %% define first baseline point
ICarbx=mean(stack(:,:,IbIdx),3);
ICarbx(ICarbx<DetectThresh)=0;
DiffMap(:,:,5)=ICarbx-Preedge;
tmp=DiffMap(:,:,5);
tmp(tmp<0 | isnan(tmp))=0;
tmp=medfilt2(tmp);
DiffMap(:,:,5)=tmp;

MatSiz=size(S.spectr);
XSiz=S.Xvalue/MatSiz(1);
YSiz=S.Yvalue/MatSiz(2);
xdat=[0:XSiz:S.Xvalue];
ydat=[0:YSiz:S.Yvalue];

%% Show DiffMaps
figure1=figure('Name',particle,'NumberTitle','off','Position',[1,1,715,869]);
title(sprintf(S.particle));
subplot(3,2,1)
plotIn=DiffMap(:,:,1);
plotIn(plotIn>3)=3;
imagesc(xdat,ydat,plotIn)
axis image
title('Inorganic')
colorbar
xlabel('X (\mum)');
ylabel('Y (\mum)');

subplot(3,2,2)
imagesc(xdat,ydat,DiffMap(:,:,2))
axis image
title('Potassium')
colorbar
xlabel('X (\mum)');
ylabel('Y (\mum)');

subplot(3,2,3)
% imagesc(xdat,ydat,DiffMap(:,:,3))
imagesc(xdat,ydat,sp2out)
axis image
title('Sp^{2}')
colorbar
xlabel('X (\mum)');
ylabel('Y (\mum)');

subplot(3,2,4)
imagesc(xdat,ydat,DiffMap(:,:,5))
axis image
% colormap('gray')
title('Carbox')
colorbar
xlabel('X (\mum)');
ylabel('Y (\mum)');

subplot(3,2,5)
imagesc(xdat,ydat,DiffMap(:,:,4))
axis image
title('CO_{3}')
colorbar
xlabel('X (\mum)');
ylabel('Y (\mum)');

subplot(3,2,6)

PartIm=DiffMap(:,:,5)+Preedge+DiffMap(:,:,3);
LabelMat=raw2mask(S);
title('Particle Map')
colorbar
axis image
xlabel('X (\mum)');
ylabel('Y (\mum)');

if figsav==1
filename=sprintf('%s%s%s%s',rootdir,sample,S.particle,'_f1_map');
saveas(gcf,filename,'tiff');
end

S=S;
S.DiffMap=DiffMap;
S.LabelMat=LabelMat;


text:
function L = ivmLogLikelihood(model, x, y);

% IVMLOGLIKELIHOOD Return the log-likelihood for the IVM.
% FORMAT
% DESC computes the log likelihood of the given data points under
% the IVM and the given noise model.
% ARG model : the model for which the log likelihood is computed.
% ARG x : the input locations.
% ARG y : the target values.
%
% SEEALSO : ivmNegLogLikelihood, ivmPosteriorMeanVar,
% noiseLogLikelihood
%
% COPYRIGHT : Neil D. Lawrence

% IVM

if nargin < 3
  % This implies evaluate for the training data.
  mu = model.mu;
  varsigma = model.varSigma;
  y = model.y;
else
  [mu, varsigma] = ivmPosteriorMeanVar(model, x);
end

L = noiseLogLikelihood(model.noise, mu, varsigma, y);

% check if there is a prior over kernel parameters
if isfield(model.kern, 'priors')
  fhandle = str2func([model.kern.type 'KernExpandParams']);
  params = fhandle(model.kern);
  for i = 1:length(model.kern.priors)
    index = model.kern.priors(i).index;
    L = L + priorLogProb(model.kern.priors(i), params(index));
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function f = cloister(xmin,xmax,ymin,ymax,n)

% CLOISTER  Generates features in a 2D cloister shape.
%   CLOISTER(XMIN,XMAX,YMIN,YMAX,N) generates a 2D cloister in the limits
%   indicated as parameters.
%
%   N is the number of rows and columns; it defaults to N = 9.
%
%   See also THICKCLOISTER.

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

if nargin < 5
    n = 9;
end

% Center of cloister
x0 = (xmin+xmax)/2;
y0 = (ymin+ymax)/2;

% Size of cloister
hsize = xmax-xmin;
vsize = ymax-ymin;
tsize = diag([hsize vsize]);

% Integer ordinates of points
outer = (-(n-3)/2 : (n-3)/2);
inner = (-(n-3)/2 : (n-5)/2);

% Outer north coordinates
No = [outer; (n-1)/2*ones(1,numel(outer))];
% Inner north
Ni = [inner ; (n-3)/2*ones(1,numel(inner))];
% East (rotate 90 degrees the North points)
E = [0 -1;1 0] * [No Ni];
% South and West are negatives of N and E respectively.
points = [No Ni E -No -Ni -E];

% Rescale
f = tsize*points/(n-1);

% Move
f(1,:) = f(1,:) + x0;
f(2,:) = f(2,:) + y0;

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

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function triangulation_test32 ( )

%*****************************************************************************80
%
%% TEST32 tests VORONOI_POLYGON_CENTROID.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2006
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  neighbor_num = 4;
  node_num = 5;

  center = 5;
  centroid_exact = [ 0.5, 0.5 ];
  neighbor_index = [ 1, 2, 3, 4 ];
  node_xy = [ ...
    0.0, 0.0; ...
    1.0, 0.0; ...
    1.0, 1.0; ...
    0.0, 1.0; ...
    0.5, 0.5 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST32\n' );
  fprintf ( 1, '  VORONOI_POLYGON_CENTROID computes the centroid of\n' );
  fprintf ( 1, '  a finite Voronoi polygon.\n' );

  centroid = voronoi_polygon_centroid ( center, neighbor_num, ...
    neighbor_index, node_num, node_xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The computed centroid ( %f  %f )\n', centroid(1:dim_num) );
  fprintf ( 1, '  The correct centroid  ( %f  %f )\n', centroid_exact(1:dim_num) );

  return
end


text:
% SP_EVALUATE_ELEMENT_LIST_PARAM: compute the basis functions, in the parametric domain, in a given list of elements.
%
%     sp = sp_evaluate_element_list_param (space, msh_elems, 'option1', value1, ...)
%
% INPUTS:
%
%    space:   object defining the space of discrete functions (see sp_vector)
%    msh_elems: msh structure containing the information of quadrature or
%               visualization points, for a given list of elements
%               (see msh_cartesian/msh_evaluate_element_list)
%   'option', value: additional optional parameters, currently available options are:
%
%              Name     |   Default value |  Meaning
%           ------------+-----------------+----------------------------------
%            value      |      true       |  compute shape_functions
%            gradient   |      false      |  compute shape_function_gradients
%            divergence |      false      |  compute shape_function_divs
%            curl       |      false      |  compute shape_function_curls
%            hessian    |      false      |  compute shape_function_hessians
%
% OUTPUT:
%
%    sp: struct representing the discrete function space, with the following fields:
%              (see the article for a detailed description)
%
%    FIELD_NAME      (SIZE)                                     DESCRIPTION
%    ncomp           (scalar)                                   number of components of the functions of the space
%    ndof            (scalar)                                   total number of degrees of freedom
%    ndof_dir        (ncomp_param x ndim matrix)                for each component, number of degrees of freedom along each direction
%    nsh_max         (scalar)                                   maximum number of shape functions per element
%    nsh             (1 x msh_col.nel vector)                   actual number of shape functions per each element
%    connectivity    (nsh_max x msh_col.nel vector)             indices of basis functions that do not vanish in each element
%    shape_functions (ncomp_param x msh_col.nqn x nsh_max x msh_col.nel)  basis functions evaluated at each quadrature node in each element
%    shape_function_gradients
%       (ncomp_param x ndim x msh_col.nqn x nsh_max x msh_col.nel) basis function gradients evaluated at each quadrature node in each element
%    shape_function_hessians
%       (ncomp_param x ndim x ndim x msh_col.nqn x nsh_max x msh_col.nel) basis function hessians evaluated at each quadrature node in each element
%    shape_function_divs (msh_col.nqn x nsh_max x msh_col.nel)     basis function divergence evaluated at each quadrature node in each element
%    shape_function_curls
%         2D:  (msh_col.nqn x nsh_max x msh_col.nel)               basis function curl evaluated at each quadrature node in each element
%         3D:  (3 x msh_col.nqn x nsh_max x msh_col.nel)
%
% Copyright (C) 2009, 2010, 2011 Carlo de Falco
% Copyright (C) 2011, 2015, 2019 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function sp = sp_evaluate_element_list_param (space, msh, varargin)

value = true;
gradient = false;
divergence = false;
curl = false;
hessian = false;
if (~isempty (varargin))
  if (~rem (length (varargin), 2) == 0)
    error ('sp_evaluate_element_list_param: options must be passed in the [option, value] format');
  end
  for ii=1:2:length(varargin)-1
    if (strcmpi (varargin {ii}, 'value'))
      value = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'gradient'))
      gradient = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'curl'))
      curl = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'divergence'))
      divergence = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'hessian'))
      hessian = varargin {ii+1};
    else
      warning ('Ignoring unknown option %s', varargin {ii});
    end
  end
end

first_der = gradient || divergence || curl;
for icomp = 1:space.ncomp_param
  sp_col_scalar(icomp) = sp_evaluate_element_list_param (space.scalar_spaces{icomp}, msh, 'value', value, 'gradient', first_der, 'hessian', hessian);
end

ndof_scalar = [sp_col_scalar.ndof];

ndof = sum (ndof_scalar);
nsh  = zeros (1, msh.nel);
connectivity = [];
aux = 0;
for icomp = 1:space.ncomp_param
  ndof_dir(icomp,:) = sp_col_scalar(icomp).ndof_dir;
  nsh = nsh + sp_col_scalar(icomp).nsh(:)';

  inds = find (sp_col_scalar(icomp).connectivity);
  sp_col_scalar(icomp).connectivity(inds) = sp_col_scalar(icomp).connectivity(inds) + aux;
  connectivity = [connectivity; sp_col_scalar(icomp).connectivity];
  aux = aux + ndof_scalar(icomp);
end

sp = struct('nsh_max', space.nsh_max, 'nsh', nsh, 'ndof', ndof,  ...
            'ndof_dir', ndof_dir, 'connectivity', connectivity, ...
            'ncomp', space.ncomp, 'ncomp_param', space.ncomp_param);

if (value)
  sp.shape_functions = zeros (sp.ncomp_param, msh.nqn, sp.nsh_max, msh.nel);
  for icomp = 1:space.ncomp_param
    indices = space.cumsum_nsh(icomp)+(1:sp_col_scalar(icomp).nsh_max);
    sp.shape_functions(icomp,:,indices,:) = sp_col_scalar(icomp).shape_functions;
  end
end

if (gradient || curl || divergence)
  shape_fun_grads = zeros (space.ncomp_param, msh.ndim, msh.nqn, sp.nsh_max, msh.nel);

  for icomp = 1:space.ncomp_param
    indices = space.cumsum_nsh(icomp)+(1:sp_col_scalar(icomp).nsh_max);
    shape_fun_grads(icomp,:,:,indices,:) = sp_col_scalar(icomp).shape_function_gradients;
  end

  if (gradient)
    sp.shape_function_gradients = shape_fun_grads;
  end

  if (divergence)
    sp.shape_function_divs = zeros (msh.nqn, sp.nsh_max, msh.nel);
    for icomp = 1:space.ncomp_param
      sp.shape_function_divs = sp.shape_function_divs + ...
        reshape (shape_fun_grads(icomp,icomp,:,:,:), msh.nqn, sp.nsh_max, msh.nel);
    end
  end

  if (curl)
    if (space.ncomp_param == 2)
      sp.shape_function_curls = reshape (shape_fun_grads(2,1,:,:,:) - ...
	 	       shape_fun_grads(1,2,:,:,:), msh.nqn, sp.nsh_max, msh.nel);
    elseif (space.ncomp_param == 3)
      shape_fun_curls = zeros (space.ncomp, msh.nqn, sp.nsh_max, msh.nel);
      for icomp = 1:space.ncomp
        ind1 = mod(icomp,3) + 1;
        ind2 = mod(ind1, 3) + 1;
        shape_fun_curls(icomp,:,:,:) = reshape (shape_fun_grads(ind2,ind1,:,:,:) - ...
                   shape_fun_grads(ind1,ind2,:,:,:), 1, msh.nqn, sp.nsh_max, msh.nel);
      end
      sp.shape_function_curls = shape_fun_curls;
    end
  end
end

if (hessian)
  sp.shape_function_hessians = zeros (space.ncomp, msh.ndim, msh.ndim, msh.nqn, sp.nsh_max, msh.nel);
  for icomp = 1:space.ncomp_param
    indices = space.cumsum_nsh(icomp)+(1:sp_col_scalar(icomp).nsh_max);
    sp.shape_function_hessians(icomp,:,:,:,indices,:) = sp_col_scalar(icomp).shape_function_hessians;
  end
end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [e1, e2, E1_l, E2_l] = hmgLinEndpoints(l,t1,t2)

% HMGLINENDPOINTS  HMG line endpoints.
%   HMGLINENDPOINTS(L,T) returns the 3D Euclidean endpoint of the
%   HMG line L at abscissa T.
%
%   [E1,E2] = HMGLINENDPOINTS(L,T1,T2) returns two 3D Euclidean endpoints
%   at abscissas T1 and T2.
%
%   [e1,e2,E1_l,E2_l] = HMGLINENDPINTS(...) returns the Jacobians of the
%   endopints wrt L.
%
%   See also HMGLINSEGMENT, HMGLIN2SEG, HMGLIN2IDPPNTS.

%   Copyright 2009 Teresa Vidal.

if nargout <= 2

    % support segment
    s  = hmgLin2seg(l);

    % support points
    p1 = s(1:3);
    p2 = s(4:6);

    % endpoints
    e1 = (1-t1)*p1 + t1*p2;
    e2 = (1-t2)*p1 + t2*p2;

else

    [s, S_l] = hmgLin2seg(l);       % support segment

    p1 = s(1:3);                    % support points
    p2 = s(4:6);

    P1_l = S_l(1:3,:);              % jacobians of sup points
    P2_l = S_l(4:6,:);

    e1 = (1-t1)*p1 + t1*p2;         % endpoints
    e2 = (1-t2)*p1 + t2*p2;

    E1_p1 = 1-t1;                   % jacobians of endpoints
    E1_p2 = t1;
    E2_p1 = 1-t2;
    E2_p2 = t2;

    E1_l = E1_p1*P1_l + E1_p2*P2_l; % chain rule
    E2_l = E2_p1*P1_l + E2_p2*P2_l;

end

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

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

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

%   # END GPL LICENSE