text:
function [w,fun,time,iter] = gistLogistic(X,y,lambda,theta,varargin)

% Generalized Iterative Shrinkage and Thresholding (GIST) with Logsitic Regression loss
%
% Non-convex optimization problem:
%
% min_w L(w) + \sum_i r_i(w)
%
% ================================ loss function ==========================
%
% L(w) = 1/n \sum_j log(1 + exp(-y_j*x_j'*w)) (n: number of samples)
%
% ================================ regularizer ============================
%
%  regtype = 1: Capped L1 regularizer (CapL1) (default)
%            r_i(w) = lambda*\min(|w_i|,theta), (theta > 0, lambda >= 0)
%
%  regtype = 2: Log Sum Penalty (LSP)
%            r_i(w) = lambda*\sum_i log(1 + |w_i|/theta), (theta > 0, lambda >= 0)
%
%  regtype = 3: Smoothly Clipped Absolute Deviation (SCAD)
%            r_i(w) = lambda*|w_i|, if |w_i|<=lambda
%            r_i(w) = (-w_i^2 + 2*theta*lambda*|w_i| - lambda^2)/(2(theta - 1)), if lambda<=|w_i|<=theta*lambda
%            r_i(w) = 0.5*(theta + 1)*lambda^2, if |w_i| > theta*lambda, (theta > 2, lambda >= 0)
%
%  regtype = 4: Minimax Concave Penalty (MCP)
%            r_i(w) = lambda*|w_i| - 0.5*w_i^2/theta, if |w_i|<=theta*lambda
%            r_i(w) = 0.5*theta*lambda^2, if |w_i| > theta*lambda, (theta >
%            0, lambda >= 0)
%
% ============================ Input ======================================
%
% X: data matrix with each row as a sample
%
% y: label vector (+1 or -1)
%
% lambda: regularization parameter
%
% theta: theresholding parameter
%
% ======================= varargin: optional settings  ====================
%
% 'regtype': nonconvex regularization type
%          1: CapL1 (default)
%          2: LSP
%          3: SCAD
%          4: MCP
%
% 'stopcriterion': stopping criterion
%                1: relative difference of objective functions
%                   is less than tol (default)
%                0: relative difference of iterative weights is less
%                   than tol
%
% 'startingpoint': starting point (default: zero vector)
%
% 'tolerance': stopping tolerance (default: 1e-5)
%
% 'maxiteration': number of maximum iteration (default: 1000)
%
% 'tinitialization': initialization of t (default: 1)
%
% 'tmin': tmin parameter (default: 1e-20)
%
% 'tmax': tmax parameter (default: 1e-20)
%
% 'eta': eta factor (default: 2)
%
% 'sigma': parameter in the line search (default: 1e-5)
%
% 'nonmonotone': nonmonotone steps in the line search (default: 5)
%
% 'stopnum': number of satisfying stopping criterion (default: 3)
%
% 'maxinneriter': number of maximum inner iteration (line search) (default: 20)
%
% ============================= Output ====================================
%
% w: output weight vector
%
% fun: a vector including all function values at each iteration
%
% time: a vector including all CPU times at each iteration
%
% iter: the number of iterative steps
%
% =========================================================================
%
% Copyright (C) 2012-2013 Pinghua Gong
%
% For any problem, please contact Pinghua Gong via pinghuag@gmail.com
%
% Last modified on March 14, 2013.
%
% Related papers:
%
% [1] Pinghua Gong, Changshui Zhang, Zhaosong Lu, Jianhua Huang, Jieping Ye,
%     A General Iterative Shrinkage and Thresholding Algorithm for Non-convex
%     Regularized Optimization Problems. ICML 2013.
%
% ========================================================================

if nargin < 4
    error('Too few input parameters!');
end

if theta <= 0 || lambda < 0
    error('\theta must be positive and \lambda must be nonneagtive!');
end

% Parse the optional inputs.
if (mod(length(varargin), 2) ~= 0 ),
    error(['Optional Parameters passed to the function ''' mfilename ''' must be passed in pairs!']);
end

% default parameter settings
regtype = 1;
[n,d] = size(X);
w0 = zeros(d,1);
stopcriterion = 1;
tol = 1e-5;
maxiter = 1000;
M = 5;
sigma = 1e-5;

t = 1;
tmin = 1e-20;
tmax = 1e20;

% % t, tmin and tmax are adaptively estimated
% Linfnorm = full(max(sum(abs(X),2)))/n; L1norm = full(max(sum(abs(X),1)))/n; Lmaxnorm = full(max(max(abs(X))))/n;
% t = max([Linfnorm*Linfnorm/d,L1norm*L1norm/n,Lmaxnorm]);
% tmin = min([Linfnorm*Linfnorm/d,L1norm*L1norm/n,Lmaxnorm]);
% tmax = min([Linfnorm*L1norm,n*d*Lmaxnorm,n*Linfnorm*Linfnorm,d*L1norm*L1norm])/(1-sigma);

eta = 2;
stopnum = 3;
maxinneriter = 20;

% Optional parameter settings
parameterCount = length(varargin)/2;

for parameterIndex = 1:parameterCount,
    parameterName = varargin{parameterIndex*2 - 1};
    parameterValue = varargin{parameterIndex*2};
    switch lower(parameterName)
        case 'regtype'
            regtype = parameterValue;
            if regtype == 3 && theta <= 2
                error('\theta must be greater than 2!');
            end
        case 'startingpoint'
            w0 = parameterValue;
        case 'stopcriterion'
            stopcriterion = parameterValue;
        case 'tolerance'
            tol = parameterValue;
        case 'maxiteration'
            maxiter = parameterValue;
        case 'nonmonotone'
            M = parameterValue;
        case 'tinitialization'
            t = parameterValue;
        case 'tmin'
            tmin = parameterValue;
        case 'tmax'
            tmax =  parameterValue;
        case 'sigma'
            sigma =  parameterValue;
        case 'eta'
            eta = parameterValue;
        case 'stopnum'
            stopnum = parameterValue;
        case 'maxinneriter'
            maxinneriter = parameterValue;
        otherwise
            error(['The parameter ''' parameterName ''' is not recognized by the function ''' mfilename '''!']);
    end
end

w = w0;
fun = zeros(maxiter+1,1); time = fun;
logist = zeros(n,1);

% Initial function value
Z = sparse(1:n,1:n,y,n,n)*X;
Zw = -Z*w; posind = (Zw > 0);
logist(posind) = 1 + exp(-Zw(posind));
logist(~posind) = 1 + exp(Zw(~posind));

temp = logist;
temp(posind) = 1./logist(posind);
temp(~posind) = (logist(~posind)-1)./logist(~posind);
grad =  -Z'*temp/n;

fun(1) = (sum(log(logist(~posind))) + sum(Zw(posind) + log(logist(posind))))/n + funRegC(w,d,lambda,theta,regtype);
time(1) = 0;

count = 0;
for iter = 1:maxiter
    tic;

    w_old = w;
    grad_old = grad;
    t = min(max(t,tmin),tmax);

    % line search
    for inneriter = 1:maxinneriter
        w = proximalRegC(w_old - grad_old/t, d, lambda/t, theta,regtype);
        dw = w - w_old;
        Zw = -Z*w; posind = (Zw > 0);
        logist(posind) = 1 + exp(-Zw(posind));
        logist(~posind) = 1 + exp(Zw(~posind));
        fun(iter+1) = (sum(log(logist(~posind))) + sum(Zw(posind) + log(logist(posind))))/n + funRegC(w,d,lambda,theta,regtype);
        if fun(iter+1) <= max(fun(max(iter-M+1,1): iter)) - 0.5*sigma*t*norm(dw)^2
            break;
        else
            t = t*eta;
        end
    end
    time(iter+1) = time(iter) + toc;

    % stopping condition
    if stopcriterion
        relativediff = abs(fun(iter) - fun(iter+1))/fun(iter+1);
    else
        relativediff = norm(w - w_old)/norm(w);
    end
    if relativediff < tol
        count = count + 1;
    else
        count = 0;
    end
    if count >= stopnum
        break;
    end

    temp = logist;
    temp(posind) = 1./logist(posind);
    temp(~posind) = (logist(~posind)-1)./logist(~posind);
    grad =  -Z'*temp/n;

    % BB rule
    st = w - w_old;
    rt = grad - grad_old;
    if norm(st)/d < 1e-12 || norm(rt)/d < 1e-12
        break;
    end
    t = st'*rt/(st'*st);

end
fun = fun(1: min(maxiter,iter)+1);
time = time(1: min(maxiter,iter)+1);

author: thomaskuestner


text:
function varargout = pnes_master(fcn, x0, opt)
%PNES_MASTER  Parameterized Nonlinear Equation Solver wrapper function.
%   [X, F, EXITFLAG, OUTPUT, JAC] = PNES_MASTER(FCN, X0, OPT)
%   [X, F, EXITFLAG, OUTPUT, JAC] = PNES_MASTER(PROBLEM)
%   A common wrapper function for numerical continuation methods for
%   solving parameterized nonlinear equations. Traces the solutions of
%   a parameterized nonlinear equation f(x) = 0, beginning from a starting
%   point x0, where f(x) has dimension n and x has dimension n+1.
%
%   In the current implementation, the last element of x is taken to be
%   the parameter lambda, where lambda = 0 corresponds to the base solution.
%
%   Inputs:
%       FCN : handle to function that evaluates the function f(x) to
%           be solved and its Jacobian, J(x). Calling syntax for this
%           function is:
%               f = FCN(x)
%               [f, J] = FCN(x)
%           For a parameterized function, f is n x 1, x is (n+1) x 1,
%           and J is the n x (n+1) matrix of partial derivatives of
%           f (rows) w.r.t. x (cols).
%       X0 : starting value, x0, of vector x ((n+1) x 1)
%       OPT : optional options structure with the following fields,
%           all of which are also optional (default values shown in
%           parentheses)
%           alg ('DEFAULT') : determines which solver to use
%               'DEFAULT' : automatic, currently there is only one
%               solver implementation, a predictor/corrector method
%           verbose (0) - controls level of progress output displayed
%               0 = no progress output
%               1-5 = increasing levels of progress output
%           nleqs_opt - options struct for NLEQS_MASTER used for
%               corrector stage (see NLEQS_MASTER for details), default
%               sets nleqs_opt.verbose to 0, otherwise to 2 if OPT.verbose > 4
%           solve_base (1) : 0/1 flag that determines whether or not to
%               run a corrector stage for initial solution point, x0
%           parameterization (3) - choice of parameterization
%               1 - natural
%               2 - arc len
%               3 - pseudo arc len
%           stop_at ('NOSE') - determines stopping criterion
%               'NOSE'     - stop when limit or nose point is reached
%               'FULL'     - trace full continuation curve
%               <lam_stop> - stop upon reaching specified target lambda value
%           max_it (2000) - maximum number of continuation steps
%           step (0.05) - continuation step size
%           adapt_step (0) - toggle adaptive step size feature
%               0 - adaptive step size disabled
%               1 - adaptive step size enabled
%           adapt_step_damping (0.7) - damping factor for adaptive step sizing
%           adapt_step_tol (1e-3) - tolerance for adaptive step sizing
%           adapt_step_ws (1) - scale factor for default initial step size
%               when warm-starting with adaptive step size enabled
%           step_min (1e-4) - minimum allowed step size
%           step_max (0.2) - maximum allowed step size
%           default_event_tol (1e-3) - default tolerance for event functions
%           target_lam_tol (0) - tolerance for target lambda detection, 0 means
%               use the value of default_event_tol
%           nose_tol (0) - tolerance for nose point detection, 0 means use
%               the value of default_event_tol
%           events (<empty>) - cell array of specs for user-defined event
%               functions, passed as MY_EVENTS arg to PNE_REGISTER_EVENTS
%               (see PNE_REGISTER_EVENTS for details).
%           callbacks (<empty>) - cell array of specs for user-defined callback
%               functions, to be passed as MY_CBACKS arg to
%               PNE_REGISTER_CALLBACKS (see PNE_REGISTER_CALLBACKS for details).
%           output_fcn (<empty>) - handle to custom output function, called by
%               PNE_CALLBACK_DEFAULT
%           plot - options for plotting of continuation curve by
%               PNE_CALLBACK_DEFAULT
%               .level (0) - control plotting of continuation curve
%                   0 - do not plot continuation curve
%                   1 - plot when completed
%                   2 - plot incrementally at each continuation step
%                   3 - same as 2, with 'pause' at each continuation step
%               .idx (<empty>) - index of quantity to plot, passed to yfcn()
%               .idx_default (<empty>) - function to provide default value
%                   for idx, if none provided
%               .xname ('lam') - name of output field holding values that
%                   determine horizontal coordinates of plot
%               .yname ('x') - name of output field holding values that
%                   determine vertical coordinates of plot
%               .xfcn (<empty>) - handle to function that maps a value from
%                   the field of the OUTPUT indicated by value of plot.xname
%                   to a horizontal coordinate for plotting
%               .yfcn (<empty>) - handle to function that maps a value from
%                   the field of the OUTPUT indicated by value of plot.yname
%                   and an index to be applied to that value into a vertical
%                   coordinate for plotting
%               .xlabel ('\lambda') - label for horizontal axis
%               .ylabel ('Variable Value') - label for vertical axis
%               .title ('Value of Variable %d') - plot title used for plot of
%                   single variable, can use %d as placeholder for var index
%               .title2 ('Value of Multiple Variables') - plot title used for
%                   plot of multiple variables
%               .legend ('Variable %d') - legend label, %d can be used as
%                   placeholder for variable index
%           warmstart (<empty>) - struct containing warm-start state, see
%               warmstart field in OUTPUT below for details of expected
%               fields
%       PROBLEM : The inputs can alternatively be supplied in a single
%           PROBLEM struct with fields corresponding to the input arguments
%           described above: fcn, x0, opt
%
%   Outputs (all optional, except X):
%       X : solution vector x
%       F : final function value, f(x)
%       EXITFLAG : exit flag
%           1 = succeeded
%           0 = failed
%       OUTPUT : output struct with the following fields:
%           corrector - output return value from NLEQS_MASTER from final
%               corrector run (see NLEQS_MASTER for details)
%           iterations - N, total number of continuation steps performed
%           events - struct array of size NE of events detected with fields:
%               k - continuation step at which event was located
%               name - name of detected event
%               idx - index(es) of critical elements in corresponding event
%                   function
%               msg - descriptive text detailing the event
%           done_msg - message describing cause of continuation termination
%           steps - (N+1) row vector of stepsizes taken at each continuation
%               step
%           lam_hat - (N+1) row vector of lambda values from prediction steps
%           lam - (N+1) row vector of lambda values from correction steps
%           max_lam - maximum value of parameter lambda (from OUTPUT.lam)
%           warmstart - optional output with information needed for
%               warm-starting an updated continuation problem, with fields:
%               cont_steps - current value of continuation step counter
%               direction - +1 or -1, for tracing of curve in same or
%                   opposite direction, respectively
%               dir_from_jac_eigs - 0/1 flag to indicate whether to use
%                   the sign of the smallest eigenvalue of the Jacobian to
%                   determine the initial direction
%               x - current solution vector
%               z - current tangent vector
%               xp - previous step solution vector
%               zp - previous step tangent vector
%               parm - function handle for current parameterization function
%               default_parm - function handle for default parameterization fcn
%               default_step - default step size
%               events - current event log, same as OUTPUT.events
%               cbs - struct containing user state information for callbacks
%                   see PNES_CALLBACK_DEFAULT for more details
%           (others) - depends on OPT.output_fcn, by default (i.e. with no
%               explicitly provided output function) includes fields:
%                   x_hat - NX x (N+1) matrix of solution values from
%                       prediction steps
%                   x - NX x (N+1) matrix of solution values from correction
%                       steps
%       JAC : final Jacobian matrix, J(x)
%
%   Calling syntax options:
%       [x, f, exitflag, output, jac] = pnes_master(fcn, x0);
%       [x, f, exitflag, output, jac] = pnes_master(fcn, x0, opt);
%       x = pnes_master(problem);
%               where problem is a struct with fields: fcn, x0, opt
%               where opt is optional
%       x = pnes_master(...);
%       [x, f] = pnes_master(...);
%       [x, f, exitflag] = pnes_master(...);
%       [x, f, exitflag, output] = pnes_master(...);
%       [x, f, exitflag, output, jac] = pnes_master(...);
%
%   Example: (based on https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/)
%       function [f, J] = f1p(x)
%           f = [  x(1)   + x(2) + 6*x(3) - 1;
%                 -x(1)^2 + x(2)          + 5   ];
%           if nargout > 1
%               J = [1 1 6; -2*x(1) 1 0];
%           end
%       end
%       problem = struct( ...
%           'fcn',  @(x)f1p(x), ...
%           'x0',   [-1; 0; 0], ...
%           'opt',  struct('verbose', 2, 'adapt_step', 1, 'step_max', 10) ...
%       );
%       [x, f, exitflag, output, jac] = pnes_master(problem);
%
%   See also PNE_CALLBACK_DEFAULT, PNE_REGISTER_CALLBACKS, PNE_REGISTER_EVENTS

%   MP-Opt-Model
%   Copyright (c) 2013-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell,
%   Shrirang Abhyankar, Argonne National Laboratory,
%   and Alexander Flueck, IIT
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%%----- input argument handling  -----
%% gather inputs
if nargin == 1 && isstruct(fcn) %% problem struct
    p = fcn;
    fcn = p.fcn;
    x0 = p.x0;
    if isfield(p, 'opt'),   opt = p.opt;    else,   opt = struct(); end
else                            %% individual args
    if nargin < 3
        opt = struct();
    end
end

%% default options
if isfield(opt, 'verbose') && opt.verbose > 4
    nleqs_opt_verbose = 2;
else
    nleqs_opt_verbose = 0;
end
dopts = struct( ...
    'alg',              'DEFAULT', ...  %% algorithm
    'verbose',          0, ...
    'nleqs_opt',        struct('verbose', nleqs_opt_verbose), ...
    'solve_base',       1, ...          %% run corrector for initial point
    'parameterization', 3, ...          %% 1 - natural, 2 - arc len, 3 - pseudo arc len
    'stop_at',          'NOSE', ...     %% 'NOSE', 'FULL', <lam_stop>
    'max_it',           2000, ...       %% maximum number of continuation steps
    'step',             0.05, ...       %% continuation step size
    'step_min',         1e-4, ...       %% minimum allowed step size
    'step_max',         0.2, ...        %% maximum allowed step size
    'adapt_step',       0, ...          %% 0/1 toggle, adaptive step size
    'adapt_step_ws',    1, ...          %% warm start inital step size scale factor
    'adapt_step_damping', 0.7, ...      %% adaptive step sizing damping factor
    'adapt_step_tol',   1e-3, ...       %% adaptive step sizing tolerance
    'default_event_tol',1e-3, ...       %% default event function tolerance
    'target_lam_tol',   0, ...          %% event function tolerance for TARGET_LAM event
    'nose_tol',         0, ...          %% event function tolerance for NOSE event
    'events',           {{}}, ...       %% user-defined event detection functions
    'callbacks',        {{}}, ...       %% user-defined callback functions
    'output_fcn',       [], ...         %% custom output fcn, default callback
    'warmstart',        [], ...         %% struct containing warm-start state
    'plot',             struct( ...     %% used by pne_callback_default() for plotting
        'level',        0, ...          %% 0 - no plot, 1 - final, 2 - steps, 3 - steps w/pause
        'idx',          [], ...         %% index of quantity to plot, passed to yfcn()
        'idx_default',  [], ...         %% fcn to provide default value for idx, if none provided
        'xfcn',         [], ...         %% fcn to compute x-coord from data
        'yfcn',         [], ...         %% fcn to compute x-coord from data, idx
        'title',        'Value of Variable %d', ... %% plot title for single var plot
        'title2',       'Value of Multiple Variables', ...  %% plot title for multiple var plot
        'xname',        'lam', ...      %% name of output field holding x vals
        'yname',        'x', ...        %% name of output field holding y vals
        'xlabel',       '\lambda', ...              %% horizontal axis label
        'ylabel',       'Variable Value', ...       %% vertical axis label
        'legend',       'Variable %d' ...           %% legend label
    ) ...
);
opt = nested_struct_copy(dopts, opt);
%% use opt.default_event_tol for NOSE and TARGET_LAM events, unless specified
if opt.target_lam_tol == 0
    opt.target_lam_tol = opt.default_event_tol;
end
if opt.nose_tol == 0
    opt.nose_tol = opt.default_event_tol;
end
if opt.max_it == 0      %% zero means use the default
    opt.max_it == dopts.max_it;
end

%% initialize
warmstarted = ~isempty(opt.warmstart);
s = struct( ...         %% container struct for various variables, flags
    'done',     0, ...      %% flag indicating continuation has terminated
    'done_msg', '', ...     %% termination message
    'warmstart',[], ...     %% warm start state to return when done (to pass
                            ...%% to subsequent warm-started call to PNES_MASTER)
    'rollback', 0, ...      %% flag to indicate a step must be rolled back
    'events',    [], ...    %% struct array for detected events
    'results',  []  );      %% results struct

%% register event and callback functions
switch opt.stop_at
    case 'NOSE'
        my_events = {{'NOSE', @pne_event_nose, opt.nose_tol}, opt.events{:}};
        my_cbacks = {{@pne_callback_nose, 51}, opt.callbacks{:}};
    case 'FULL'
        my_events = {{'TARGET_LAM', @pne_event_target_lam, opt.target_lam_tol}, opt.events{:}};
        my_cbacks = {{@pne_callback_target_lam, 50}, opt.callbacks{:}};
    otherwise   %% numeric, stop at target lam or nose point, whichever is 1st
        my_events = {{'TARGET_LAM', @pne_event_target_lam, opt.target_lam_tol}, ...
                     {'NOSE', @pne_event_nose, opt.nose_tol}, ...
                        opt.events{:}};
        my_cbacks = {{@pne_callback_target_lam, 50}, ...
                     {@pne_callback_nose, 51}, ...
                        opt.callbacks{:}};
end
my_cbacks{end+1} = {@pne_callback_default, 0};
reg_ev = pne_register_events(my_events, opt);   %% registered event functions
reg_cb = pne_register_callbacks(my_cbacks);     %% registered callback functions
nef = length(reg_ev);   %% number of registered event functions
ncb = length(reg_cb);   %% number of registered callback functions

%% initialize continuation step counter
if warmstarted
    cont_steps = opt.warmstart.cont_steps + 1;
    if opt.verbose
        fprintf('... CONTINUATION RESUMED\n');
    end
else
    cont_steps = 0;
    if opt.verbose
        v = mpomver('all');
        fprintf('\nMP-Opt-Model Version %s, %s', v.Version, v.Date);
        fprintf(' -- Predictor/Corrector Continuation Method\n');
    end
end

%% solve corrector step for base point
if opt.solve_base && ~warmstarted
    cfcn = @(xx)pne_corrector_fcn(xx, fcn, @pne_pfcn_natural, x0, 0, []);
    [x, f, exitflag, out] = nleqs_master(cfcn, x0, opt.nleqs_opt);
    if exitflag
        if opt.verbose > 1
            fprintf('step %3d  :                          lambda = %6.3f, %2d corrector steps\n', cont_steps, x0(end), out.iterations);
        end
    else
        s.done = 1;
        s.done_msg = sprintf('base solution did not converge in %d iterations', out.iterations);
        if opt.verbose
            fprintf('%s\n', s.done_msg);
        end
    end
else
    x = x0;     %% ignored for warmstart, overwritten by warmstart.cx
end

%% initialize numerical continuation
if ~s.done
    locating = 0;   %% flag to indicate that an event interval was detected,
                    %% but the event has not yet been located
    rb_cnt_ef = 0;  %% counter for rollback steps triggered by event function intervals
    rb_cnt_cb = 0;  %% counter for rollback steps triggered directly by callbacks

    if warmstarted
        manual_direction_switch = 0;    %% set to 1 to manually prompt for
                                        %% direction change upon warmstart
        ws = opt.warmstart;
        x = ws.x;           %% starting value solution vector
        z = ws.z;           %% starting value of tangent vector
        step = 0;           %% re-solve current point
        parm = ws.parm;
        direction = ws.direction;
        default_parm = ws.default_parm;
        default_step = ws.default_step;
        cbs = ws.cbs;
        event_log = ws.events;
        if manual_direction_switch
            %% decide whether to switch directions
            reply = input('Switch directions? Y/N [N]:','s');
            if strcmp(upper(reply), 'Y')
                direction = -direction;
            end
        elseif isfield(ws, 'dir_from_jac_eigs') && ws.dir_from_jac_eigs
            %% attempt to determine direction from smalles Jacobian eigenvalue
            [~, J] = fcn(x);
            eigs_opt.tol = 1e-3;
            eigs_opt.maxit = 2*length(x);
            direction = sign(z(end) * ...
                        min(real(eigs(J(:,1:end-1), 1, 'sr', eigs_opt))));
        end

        if opt.adapt_step   %% hey, maybe slow down, things might have changed
            default_step = default_step * opt.adapt_step_ws;
        end
    else
        %% initialize parameterization function
        switch opt.parameterization
            case 1
                parm = @pne_pfcn_natural;           %% NAT
            case 2
                parm = @pne_pfcn_arc_len;           %% ARC
            case 3
                parm = @pne_pfcn_pseudo_arc_len;    %% PAL
            otherwise
                error('pnes_master: OPT.parameterization (= %d) must be 1, 2, or 3', opt.parameterization);
        end

        %% finish initializing tangent vector
        direction = 1;
        z0 = zeros(length(x), 1); z0(end) = direction;  %% +ve lambda direction
        z = pne_tangent(x, x, z0, fcn, parm, direction);

        step = opt.step;
        default_step = step;
        default_parm = parm;
        cbs = [];
        event_log = [];
    end

    %% initialize state struct for current continuation step
    cx = struct( ...        %% current state
        'x_hat',        x, ...      %% predicted solution value
        'x',            x, ...      %% corrected solution value
        'z',            z, ...      %% normalized tangent vector
        'default_step', default_step, ...   %% default step size
        'default_parm', default_parm, ...   %% default parameterization
        'this_step', [], ...        %% step size for this step only
        'this_parm', [], ...        %% parameterization for this step only
        'step', step, ...           %% current step size
        'parm', parm, ...           %% current parameterization
        'events', event_log, ...    %% event log
        'cbs', cbs, ...             %% user-defined callback state
        'efv', [] ...               %% event function values
    );

    %% initialize event function values
    cx.efv = cell(nef, 1);
    for k = 1:nef
        cx.efv{k} = reg_ev(k).fcn(cx, opt);
    end

    if warmstarted  %% no need to initialize callbacks
        %% initialize state for previous continuation step
        px = cx;
        px.x = ws.xp;   %% use warm start value for solution value
        px.z = ws.zp;   %% use warm start value for tangent
    else
        %% invoke callbacks - "initialize" context
        for k = 1:ncb
            [nx, cx, s] = reg_cb(k).fcn(cont_steps, cx, cx, cx, s, opt);
        end
        cont_steps = cont_steps + 1;

        %% check for case with base and target the same
        %% evaluate function at lambda = 0 (base)
        if opt.solve_base
            fb = f(1:end-1);
        else
            fb = fcn(x);
            exitflag = 1;
        end

        %% evaluate function at lambda = 1 (target)
        xt = x;
        xt(end) = 1;
        ft = fcn(xt);
        if norm(fb - ft, Inf) < 1e-12
            s.done = 1;
            s.done_msg = 'base and target functions are identical';
        end

        %% initialize state for previous continuation step
        px = cx;
    end
end

%%-----  run numerical continuation  -----
while ~s.done
    %% initialize next candidate with current state
    nx = cx;

    %% predictor step
    nx.x_hat = cx.x + cx.step * cx.z;

    %% corrector step
    cfcn = @(xx)pne_corrector_fcn(xx, fcn, cx.parm, cx.x, cx.step, cx.z);
    [nx.x, f, exitflag, out] = nleqs_master(cfcn, nx.x_hat, opt.nleqs_opt);
    if ~exitflag        %% corrector failed
        s.done = 1;
        s.done_msg = sprintf('Corrector did not converge in %d iterations.', out.iterations);
        if opt.verbose
            fprintf('step %3d  : %s stepsize = %-9.3g lambda = %6.3f  corrector did not converge in %d iterations\n', cont_steps, pne_ptag(cx.parm), cx.step, nx.x(end), out.iterations);
        end
        cont_steps = max(cont_steps - 1, 1);    %% go back to last step, but not to 0
        break;
    end

    %% compute new tangent direction, based on current or prev state: tx
    if nx.step == 0     %% if this is a re-do step, cx and nx are the same
        tx = px;            %% so use px as the previous state
    else                %% otherwise
        tx = cx;            %% use cx as the previous state
    end
    nx.z = pne_tangent(nx.x, tx.x, tx.z, fcn, nx.parm, direction);
    direction = 1;      %% continue in same direction

    %% detect events
    for k = 1:nef
        nx.efv{k} = reg_ev(k).fcn(nx, opt); %% update event function values
    end
    [s.rollback, s.events, nx.efv] = ...
        pne_detect_events(reg_ev, nx.efv, cx.efv, nx.step);

    %% adjust step-size to locate event function zero, if necessary
    if s.rollback               %% current step overshot
        %% roll back & initialize next step size based on rollback and previous
        rx = nx;                    %% save state we're rolling back from
        rx_evnts = s.events;        %% and critical event info
        cx.this_step = s.events.step_scale * rx.step;
        cx.this_parm = rx.parm;     %% keep same parameterization as last step
        locating = 1;               %% enter "locating" mode (or stay in it)
        rb_cnt_ef = rb_cnt_ef + 1;  %% increment rollback counter for ef intervals
        if rb_cnt_ef > 26
            s.done = 1;
            s.done_msg = sprintf('Could not locate %s event!', s.events.name);
        end
        if opt.verbose > 3
            loc_msg = sprintf('OVERSHOOT  : f = [%g, <<%g>>], step <-- %.4g', ...
                        cx.efv{s.events.eidx}(s.events.idx(1)), ...
                        rx.efv{s.events.eidx}(s.events.idx(1)), cx.this_step);
        end
    elseif locating
        if s.events(1).zero      %% found the zero!
            %% step size will be reset to previously used default step size
            locating = 0;           %% exit "locating" mode
            rb_cnt_ef = 0;          %% reset rollback counter for ef intervals
            if opt.verbose > 3
                loc_msg = sprintf('ZERO!      : f = %g, step <-- %.4g', ...
                    nx.efv{rx_evnts.eidx}(rx_evnts.idx(1)), nx.default_step);
            end
        else                    %% prev rollback undershot
            %% initialize next step size based on critical event function
            %% values from prev rollback step and current step
            rx_efv = rx.efv{rx_evnts.eidx}(rx_evnts.idx(1));
            cx_efv = nx.efv{rx_evnts.eidx}(rx_evnts.idx(1));
            step_scale = cx_efv / (cx_efv - rx_efv);
            nx.this_step = step_scale * (rx.step - nx.step);
            rb_cnt_ef = 0;          %% reset rollback counter for ef intervals
            if opt.verbose > 3
                loc_msg = sprintf('UNDERSHOOT : f [<<%g>>, %g], step <-- %.4g', ...
                    cx_efv, rx_efv, nx.this_step);
            end
        end
    else                    %% normal step, not locating anything
        loc_msg = '';
    end

    %% invoke callbacks - "iterations" context
    rb = s.rollback;
    for k = 1:ncb
        [nx, cx, s] = reg_cb(k).fcn(cont_steps, nx, cx, px, s, opt);
    end
    if ~rb && s.rollback    %% rollback triggered by callback (vs event function interval)
        rb_cnt_cb = rb_cnt_cb + 1;  %% increment rollback counter for callbacks
        if rb_cnt_cb > 26
            s.done = 1;
            s.done_msg = 'Too many rollback steps triggered by callbacks!';
        end
    else
        rb_cnt_cb = 0;              %% reset rollback counter for callbacks
    end

    %% print iteration information
    if opt.verbose > 1
        %% set label for rollback step counter
        if rb_cnt_ef
            sub_step = char('a' + rb_cnt_ef - 1);
        elseif rb_cnt_cb
            sub_step = char('A' + rb_cnt_cb - 1);
        else
            sub_step = ' ';
        end

        fprintf('step %3d%s : %s stepsize = %-9.3g lambda = %6.3f', cont_steps, sub_step, pne_ptag(cx.parm), cx.step, nx.x(end));
        if opt.verbose < 5
            fprintf('  %2d corrector steps', out.iterations);
        end
        if s.rollback
            fprintf(' ^ ROLLBACK\n');
        else
            fprintf('\n');
        end
        if opt.verbose > 3 && ~isempty(loc_msg)
            fprintf('    LOCATING -- %s\n', loc_msg);
        end
    end

    %% log events
    for k = 1:length(s.events)
        if s.events(k).log
            e = struct( 'k', cont_steps, ...
                        'name', s.events(k).name, ...
                        'idx', s.events(k).idx, ...
                        'msg', s.events(k).msg   );
            if isempty(nx.events)
                nx.events = e;
            else
                nx.events(end+1) = e;
            end
        end
        if (opt.verbose > 2 && s.events(k).log) || ...
                (opt.verbose > 3 && s.events(k).eidx)
            fprintf('    %s\n', s.events(k).msg);
        end
    end

    %% adapt stepsize if requested and not terminating, locating a zero
    %% or warm starting
    if opt.adapt_step && ~s.done && ~locating && ~s.events(1).zero && nx.step ~= 0
        pred_error = norm(nx.x - nx.x_hat, Inf);

        %% new nominal step size is current size * tol/err, but we reduce
        %% the change from the current size by a damping factor and limit
        %% increases to a factor of 2
        step_scale = min(2, 1 + opt.adapt_step_damping * ...
                        (opt.adapt_step_tol/pred_error - 1));
        nx.default_step = nx.step * step_scale;

        %% limit step-size
        if nx.default_step > opt.step_max
            nx.default_step = opt.step_max;
        end
        if nx.default_step < opt.step_min
            nx.default_step = opt.step_min;
        end
    end

    %% if this is a normal step
    if ~s.rollback
        px = cx;    %% save current state before update
        cx = nx;    %% update current state to next candidate
        if ~s.done
            if cont_steps >= opt.max_it
                s.done = 1;
                s.done_msg = sprintf('Reached maximun number of continuation steps (opt.max_it = %d)', opt.max_it);
            else
                cont_steps = cont_steps + 1;
            end
        end
    end

    %% set step size and parameterization, from one-time or defaults
    if isempty(cx.this_step)
        cx.step = cx.default_step;
    else
        cx.step = cx.this_step;
        cx.this_step = [];      %% disable for next time
    end
    if isempty(cx.this_parm)
        cx.parm = cx.default_parm;
    else
        cx.parm = cx.this_parm;
        cx.this_parm = [];      %% disable for next time
    end
end     %% while ~s.done

%% invoke callbacks - "final" context
s.results = struct();   %% initialize results struct
for k = 1:ncb
    [nx, cx, s] = reg_cb(k).fcn(-cont_steps, nx, cx, px, s, opt);
end
output = s.results;
output.done_msg = s.done_msg;
output.events = cx.events;  %% copy eventlog to results
output.corrector = out;     %% output from last corrector run

%% prepare to exit
if isempty(s.warmstart)
    if opt.verbose
        fprintf('CONTINUATION TERMINATION: %s\n', s.done_msg);
    end
else
    %% save warmstart values
    ws = s.warmstart;
    ws.cont_steps = cont_steps;
    ws.direction = direction;

    %% from state at current step
    ws.x = cx.x;            %% state from current step
    ws.z = cx.z;            %% tangent vector from current step
    ws.parm = cx.parm;
    ws.default_parm = cx.default_parm;
    ws.default_step = cx.default_step;
    ws.events = cx.events;
    ws.cbs = cx.cbs;

    %% from state at previous step
    ws.xp = px.x;           %% state from previous step
    ws.zp = px.z;           %% tangent vector from previous step

    output.warmstart = ws;

    if opt.verbose
        fprintf('%s : CONTINUATION SUSPENDED ...\n', s.done_msg);
    end
end

%% output arguments
if nargout > 4
    [f, J] = fcn(cx.x);
elseif nargout > 1
    f = fcn(cx.x);
end
varargout{1} = cx.x;
if nargout > 1
    varargout{2} = f;
    if nargout > 2
        varargout{3} = exitflag;
        if nargout > 3
            varargout{4} = output;
            if nargout > 4
                varargout{5} = J;
            end
        end
    end
end

%%-----  pne_corrector_fcn  -----
%% fcn(x) combined with parameterization constraint
function [fp, dfp] = pne_corrector_fcn(x, fcn, parm, cx_x, step, z)
if nargout < 2
    fp = [ fcn(x); parm(x, cx_x, step, z) ];
else
    [f, df] = fcn(x);
    [p, dp] = parm(x, cx_x, step, z);
    fp = [f; p];
    dfp = [df; dp];
end

%%-----  pne_tangent  -----
%% find normalized tangent vector
function z = pne_tangent(x, xp, zp, fcn, parm, direction)
[f, df] = fcn(x);
[p, dp] = parm(x, xp, 0, zp);
rhs = [ zeros(length(f), 1); direction ];
z = [df; dp] \ rhs;
z = z / norm(z);    %% normalize it

%%-----  pne_ptag  -----
%% return 3-letter string to indicate parameterization scheme
%% NAT - natural, ARC - arc length, PAL - pseudo arc length
function ptag = pne_ptag(parm)
switch func2str(parm)
    case 'pne_pfcn_natural'
        ptag = 'NAT';
    case 'pne_pfcn_arc_len'
        ptag = 'ARC';
    case 'pne_pfcn_pseudo_arc_len'
        ptag = 'PAL';
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%% TOMLAB PROPT Parameter Estimation Problems
% http://tomdyn.com/parameter_estimation_dynamic_systems.html
clc
clear

%% http://tomdyn.com/examples/catalyticCracking.html
clc
ode = @(t,z,p) [-(p(1)+p(3))*z(1).^2
                p(1)*z(1).^2-p(2)*z(2)];
z0 = [1;0]; %ic

% Various constants and expressions
y1meas = [1.0;0.8105;0.6208;0.5258;0.4345;0.3903;...
    0.3342;0.3034;0.2735;0.2405;0.2283;0.2071;0.1669;...
    0.153;0.1339;0.1265;0.12;0.099;0.087;0.077;0.069];
y2meas = [0;0.2;0.2886;0.301;0.3215;0.3123;0.2716;...
    0.2551;0.2258;0.1959;0.1789;0.1457;0.1198;0.0909...
    ;0.0719;0.0561;0.046;0.028;0.019;0.014;0.010];
tmeas = [0;0.025;0.05;0.075;0.1;0.125;...
    0.15;0.175;0.2;0.225;0.25;0.3;0.35;0.4;...
    0.45;0.5;0.55;0.65;0.75;0.85;0.95];

%Build OPTI Object
theta0 = [10;8;1]; %inital parameter guess
opts = optiset('display','iter');
Opt = opti('ode',ode,'data',tmeas,[y1meas' y2meas'],'theta0',theta0,'z0',z0,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% http://tomdyn.com/examples/isometrizationAlpha.html
clc
ode = @(t,z,p) [-(p(1)+p(2))*z(1)
                p(1)*z(1)
                p(2)*z(1)-(p(3)+p(4))*z(3)+p(5)*z(5)
                p(3)*z(3)
                p(4)*z(3)-p(5)*z(5)];
z0 = [100;0;0;0;0]; %ic

y1meas = [88.35; 76.4; 65.1; 50.4; 37.5; 25.9; 14.0; 4.5];
y2meas = [7.3; 15.6; 23.1; 32.9; 42.7; 49.1; 57.4; 63.1];
y3meas = [2.3; 4.5; 5.3; 6.0; 6.0; 5.9; 5.1; 3.8];
y4meas = [0.4; 0.7; 1.1; 1.5; 1.9; 2.2; 2.6; 2.9];
y5meas = [1.75; 2.8; 5.8; 9.3; 12.0; 17.0; 21.0; 25.7];
tmeas  = [1230; 3060; 4920; 7800; 10680; 15030; 22620; 36420];

%Build OPTI Object
theta0 = [0;0;0;0;0]; %inital parameter guess
opts = optiset('display','iter','dynamicOpts',optidynset('sensitivity','cs','initialT',0));
Opt = opti('ode',ode,'data',tmeas,[y1meas' y2meas' y3meas' y4meas' y5meas'],'theta0',theta0,'z0',z0,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% http://tomdyn.com/examples/marinePopulation.html
clc
ode = @(t,z,p) [[0; p(9:end)].*[0; z(1:7)] - (p(1:8)+[p(9:end);0]).*z];

ymeas = [20000 17000 10000 15000 12000 9000 7000 3000
    12445 15411 13040 13338 13484 8426 6615 4022
     7705 13074 14623 11976 12453 9272 6891 5020
     4664  8579 12434 12603 11738 9710 6821 5722
     2977  7053 11219 11340 13665 8534 6242 5695
     1769  5054 10065 11232 12112 9600 6647 7034
      943  3907  9473 10334 11115 8826 6842 7348
      581  2624  7421 10297 12427 8747 7199 7684
      355  1744  5369  7748 10057 8698 6542 7410
      223  1272  4713  6869  9564 8766 6810 6961
      137   821  3451  6050  8671 8291 6827 7525
       87   577  2649  5454  8430 7411 6423 8388
       49   337  2058  4115  7435 7627 6268 7189
       32   228  1440  3790  6474 6658 5859 7467
       17   168  1178  3087  6524 5880 5562 7144
       11    99   919  2596  5360 5762 4480 7256
        7    65   647  1873  4556 5058 4944 7538
        4    44   509  1571  4009 4527 4233 6649
        2    27   345  1227  3677 4229 3805 6378
        1    20   231   934  3197 3695 3159 6454
        1    12   198   707  2562 3163 3232 5566];
tmeas  = 0:0.5:10;

z0 = ymeas(1,:); %ic

%Build OPTI Object
theta0 = [zeros(8,1);zeros(7,1)]; %inital parameter guess
opts = optiset('display','iter','dynamicOpts',optidynset('sensitivity','cs'));
Opt = opti('ode',ode,'data',tmeas,ymeas,'theta0',theta0,'z0',z0,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% http://tomdyn.com/examples/methanolToHydrocarbons.html
clc
ode = @(t,z,p) [-(2*p(2)-(p(1)*z(2))./((p(2)+p(5))*z(1)+z(2))+p(3)+p(4)).*z(1)
                (p(1)*z(1).*(p(2)*z(1)-z(2)))./((p(2)+p(5))*z(1)+z(2))+p(3)*z(1)
                (p(1)*z(1).*(z(2)+p(5)*z(1)))./((p(2)+p(5))*z(1)+z(2))+p(4)*z(1)];

y1meas = [0.7085;0.5971;0.5537;0.3684;0.1712;...
    0.1198;0.0747;0.0529;0.0415;0.0261;0.0208;...
    0.0085;0.0053;0.0019;0.0018];
y2meas = [0.1621;0.1855;0.1989;0.2845;0.3491;...
    0.3098;0.3576;0.3347;0.3388;0.3557;0.3483;...
    0.3836;0.3611;0.3609;0.3485];
y3meas = [0.0811;0.0965;0.1198;0.1535;0.2097;...
    0.2628;0.2467;0.2884;0.2757;0.3167;0.2954;...
    0.295;0.2937;0.2831;0.2846];
tmeas = [0.05;0.065;0.08;0.123;0.233;0.273;...
    0.354;0.397;0.418;0.502;0.553;...
    0.681;0.75;0.916;0.937];

z0 = [1;0;0]; %ic
lb = ones(5,1)*sqrt(eps);
ub = 10*ones(5,1);

%Build OPTI Object
theta0 = [1;1;1;1;1]; %inital parameter guess
opts = optiset('display','iter','solver','nl2sol','dynamicOpts',optidynset('initialT',0,'sensitivity','cs'));
Opt = opti('ode',ode,'bounds',lb,ub,'data',tmeas,[y1meas y2meas y3meas],'theta0',theta0,'z0',z0,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% http://tomdyn.com/examples/parameterEstimation.html
clc
ode = @(t,z,p) [p(1)*z(2); 1-2*z(2)-z(1)];

y1meas = [0.264;0.594;0.801;0.959];
tmeas = [1;2;3;5];

z0 = [NaN;NaN]; %ic

lb = -1.5*ones(3,1); lb(1)=1; %not actually interested in p1, but needed for OPTI
ub = 1.5*ones(3,1); ub(1)=1;

%Build OPTI Object
theta0 = [1;0;0]; %inital parameter guess
opts = optiset('display','iter','solver','nl2sol','dynamicOpts',optidynset('stateIndex',1,'initialT',0,'sensitivity','cs'));
Opt = opti('ode',ode,'bounds',lb,ub,'data',tmeas,y1meas,'theta0',theta0,'z0',z0,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% FM Comparison (with SHOT descriptors) - First run preprocess_artist_model.m
clear all; close all; clc

addpath(genpath('./'))
addpath(genpath('./../Tools/'))

model_0 = load('./tf_artist/model_0.mat');
model_1 = load('./tf_artist/model_1.mat');
A = model_0.model_evecs_trans*model_0.model_shot;
B = model_1.model_evecs_trans*model_1.model_shot;
C = mldivide(A',B')'; %B = C*A
P = model_1.model_evecs*C*model_0.model_evecs_trans;
P = normc(P);
[~, matches_0_1] = max(P,[],1);

mesh_0 = load('./artist_models/model_0_remesh'); mesh_0 = mesh_0.part;
mesh_1 = load('./artist_models/model_1_remesh'); mesh_1 = mesh_1.model;
D_0 = load('./tf_artist/model_0_dist.mat'); D_0 = D_0.D;
D_1 = load('./tf_artist/model_1_dist.mat'); D_1 = D_1.D;

colors = create_colormap(mesh_1,mesh_1);
figure;
subplot(1,2,1); colormap(colors);
plot_scalar_map(mesh_1,[1: size(mesh_1.VERT,1)]');freeze_colors;title('Target');

subplot(1,2,2); colormap(colors(matches_0_1,:));
plot_scalar_map(mesh_0,[1: size(mesh_0.VERT,1)]');freeze_colors;title('Source');

save('matches_axiomatic_SHOT_FM.mat','matches_0_1');

%% PMF Refinement
[~, matches_filter] = ...
    my_mfilter(mesh_0.n, mesh_1.n, matches_0_1 , [], 150000/400, D_0, D_1, 0); %corr is not used and set to [];

figure;
subplot(1,2,1); colormap(colors);
plot_scalar_map(mesh_1,[1: size(mesh_1.VERT,1)]');freeze_colors;title('Target');

subplot(1,2,2); colormap(colors(matches_filter,:));
plot_scalar_map(mesh_0,[1: size(mesh_0.VERT,1)]');freeze_colors;title('Source');

save('matches_axiomatic_SHOT_FM_PMF.mat','matches_filter');

author: OshriHalimi


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%

J = (1 / m) * sum( -y'*log(sigmoid(X*theta)) - (1-y)'*log( 1 - sigmoid(X*theta)) );

grad = (1 / m) * sum( X .* repmat((sigmoid(X*theta) - y), 1, size(X,2)) );

% =============================================================

end

author: Borye


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
%GRADIENTDESCENTMULTI Performs gradient descent to learn theta
%   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
n = size(X , 2);
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCostMulti) and gradient here.
    %

	H = X * theta;
	T = zeros(n , 1);
	for i = 1 : m,
		T = T + (H(i) - y(i)) * X(i,:)';
	end

	theta = theta - (alpha * T) / m;

    % ============================================================

    % Save the cost J in every iteration
    J_history(iter) = computeCostMulti(X, y, theta);

end

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtHmxCompressorBox.m                         |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Compare compressor for separated boxes        |
%|  `---'  |                                                              |
%+========================================================================+

clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Taille de la boite
edg = 1;

% Precision ou nombre de points de quadrature polynomiale
tol = 1e-3;

% Regular cube
n       = 10;
x       = 0:1/(n-1):1;
[x,y,z] = meshgrid(x,x,x);

% Recepteurs
Nx     = 1000;
% X      = edg*rand(Nx,3);
X      = edg*[x(:),y(:),z(:)];

% Emmeteurs
Ny     = 1000;
% Y      = edg*rand(Ny,3);
Y      = edg*[x(:),y(:),z(:)];
Y(:,1) = 2*edg + Y(:,1);

% Potentiel aleatoire aux emmeteurs
V  = (-1-1i) + (2+2i)*rand(Ny,1);

% Representation graphique
figure(1)
plot3(X(:,1),X(:,2),X(:,3),'*b',Y(:,1),Y(:,2),Y(:,3),'r*')
grid on
axis equal

% Noyau de green Helmholtz
k     = 5;
rxy   = @(X,Y) sqrt( (X(:,1)-Y(:,1)).^2 + (X(:,2)-Y(:,2)).^2 + (X(:,3)-Y(:,3)).^2 );
green = @(X,Y) exp(1i*k*rxy(X,Y))./rxy(X,Y);

% Gridding
[I,J] = ndgrid(1:Nx,1:Ny);

% Noyau de green
Gxy = green(X(I,:),Y(J,:));
Gxy = reshape(Gxy,Nx,Ny);
% Gxy = [Gxy zeros(Nx,Ny) ; ...
%     zeros(Nx,Ny) Gxy];
% Gxy = [Gxy zeros(Nx,2*Ny) ; ...
%     zeros(Nx,Ny) Gxy zeros(Nx,Ny) ; ...
%     zeros(Nx,2*Ny) Gxy];
% figure
% imagesc(Gxy)

% Calcul du rang
tic
rk = rank(Gxy,tol);
toc
rk

% Compression ACA, pivotage total
tic
[A,B] = hmxACA(Gxy,tol);
toc
size(A)
norm(A*B-Gxy,'fro')./norm(Gxy,'fro')

% Compression ACA, pivotage partiel
tic
[A,B] = hmxACA(X,Y,@(X,Y) green(X,Y),tol);
toc
size(A)
norm(A*B-Gxy,'fro')./norm(Gxy,'fro')

% % Recompression ACA, pivotage total
% tic
% [Ap,Bp] = hmxACA(A,B,tol);
% toc
% size(Ap)
% norm(Ap*Bp-Gxy,'fro')./norm(Gxy,'fro')

% Recompression RSVD
tic
[Ap,Bp,flag] = hmxRSVD(A,B,tol);
toc
size(Ap)
norm(Ap*Bp-Gxy,'fro')./norm(Gxy,'fro')

% Recompression QRSVD
tic
[Ap,Bp] = hmxQRSVD(A,B,tol);
toc
size(Ap)
norm(Ap*Bp-Gxy,'fro')./norm(Gxy,'fro')

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function rf=v_lpcaa2rf(aa)
%V_LPCAA2RF LPC: Convert vocal tract areas to reflection coefficients RF=(AA)

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_lpcaa2rf.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p2]=size(aa);
rf = (aa(:,2:p2)-aa(:,1:p2-1))./(aa(:,2:p2)+aa(:,1:p2-1));

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function twRatio = computeTrueThrustToWeight(bodyInfo, totalThrust, totalMass, altitude)
    %thrust must be in Newtons and mass must be in kg!

    gSlAccel = (bodyInfo.gm / ((bodyInfo.radius+altitude)^2))*1000; %m/s^2
    totalSlWeight = totalMass*gSlAccel; %kg*m/s^2 = N

    twRatio = totalThrust/totalSlWeight;
end


text:
function [H]=hessianScalar(varargin)

% function [H]=hessianScalar(U,v,cellOpt)
% ------------------------------------------------------------------------
%
% This function computes the Hessian matrix of the scalar function U for
% each point in U. U may be a vector, a 2D matrix or a 3D matrix. The
% vector v denotes the points spacing between the data entries in U. If v
% is not supplied the spacing is assumed to be homogeneous and unity. If
% the input is n dimensional array consisting of m entries then the output
% is a matrix (if cellOpt==0) the size of mx(n^2) (whereby the colum
% entries define the entries in a Hessian matrix and row entries relate to
% elements in the input array). If cellOpt==1 then the output is reformed
% into a cell array that matches the size of the input aray. Each cell
% entry then contains the nxn Hessian matrix.
%
% See also: |gradient|,|hessian|,|jacobian|,|cellEig|
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2014/10/10 Updated
%------------------------------------------------------------------------

%% Parse input

switch nargin
    case 1
        U=varargin{1};
        v=[];
        cellOpt=0;
    case 2
        U=varargin{1};
        v=varargin{2};
        cellOpt=0;
    case 3
        U=varargin{1};
        v=varargin{2};
        cellOpt=varargin{3};
end
%%

nDims=ndims(U);
if isvector(U)
    nDims=1;
end

if isempty(v)
    v=ones(1,nDims); %Use default
end

if nDims<=3
    switch nDims
        case 1
            %Compute first order derivative
            [dUdx] = gradient(U,v);

            %Compute second order derivatives
            [dUdxdx] = gradient(dUdx,v);

            H_mat=dUdxdx(:);
        case 2
            %Compute first order derivative
            [dUdx,dUdy] = gradient(U,v(1),v(2));

            %Compute second order derivatives
            [dUdxdx,dUdxdy] = gradient(dUdx,v(1),v(2));
            [dUdydx,dUdydy] = gradient(dUdy,v(1),v(2));

            H_mat=[dUdxdx(:) dUdxdy(:)...
                dUdydx(:) dUdydy(:)];
        case 3
            %Compute first order derivative
            [dUdx,dUdy,dUdz] = gradient(U,v(1),v(2),v(3));

            %Compute second order derivatives
            [dUdxdx,dUdxdy,dUdxdz] = gradient(dUdx,v(1),v(2),v(3));
            [dUdydx,dUdydy,dUdydz] = gradient(dUdy,v(1),v(2),v(3));
            [dUdzdx,dUdzdy,dUdzdz] = gradient(dUdz,v(1),v(2),v(3));

            H_mat=[dUdxdx(:) dUdxdy(:) dUdxdz(:)...
                   dUdydx(:) dUdydy(:) dUdydz(:)...
                   dUdzdx(:) dUdzdy(:) dUdzdz(:)];
    end

    if cellOpt==1
        H_mat=reshape(H_mat',nDims,nDims,size(H_mat,1));

        %Cell array of Hessian matrices
        H=reshape(mat2cell(H_mat,nDims,nDims,ones(size(H_mat,3),1)),size(U));
    else
        H=H_mat;
    end
else
    warning('Function only defined for 1D, 2D or 3D scalar valued arrays');
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef RoiPoolingFreeform < dagnn.Layer
    % This layer has to be used AFTER the ROIPoolingFreeform.
    % It applies a mask to the roi-pooled activations, taking either fg, bg
    % or both of the image.
    %
    % inputs are: rois, masks, blobMasks
    % outputs are: rois
    %
    % Copyright by Holger Caesar, 2015

    properties
        combineFgBox = false;
    end

    properties (Transient)
        mask
    end

    methods
        function outputs = forward(obj, inputs, params) %#ok<INUSD>
            assert(numel(inputs) == 3);
            [outputs{1}, obj.mask] = roiPooling_freeform_forward(inputs{1}, inputs{2}, inputs{3}, obj.combineFgBox);
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs) %#ok<INUSL>
            assert(numel(derOutputs) == 1);
            derInputs{1} = roiPooling_freeform_backward(derOutputs{1}, obj.combineFgBox);
            derInputs{2} = [];
            derInputs{3} = [];
            derParams = {} ;
        end

        function obj = RoiPoolingFreeform(varargin)
            obj.load(varargin) ;
        end
    end
end

author: nightrome

repo: matconvnet-calvin


text:
% POP_LOADBCI - import a BCI2000 ascii file into EEGLAB
%
% Usage:
%   >> OUTEEG = pop_loadbci( filename, srate );
%
% Inputs:
%   filename       - file name
%   srate          - sampling rate
%
% Outputs:
%   OUTEEG         - EEGLAB data structure
%
% Author: Arnaud Delorme, CNL / Salk Institute, 9 July 2002
%
% See also: EEGLAB

% Copyright (C) 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [EEG, command] = pop_loadbci(filename, srate);

    EEG = [];
    command = '';

    if nargin < 1
        % ask user
        [filename, filepath] = uigetfile('*.*', 'Choose a BCI file -- pop_loadbci');
        drawnow;
        if filename == 0 return; end
        filename = [filepath filename];
        promptstr    = { 'Sampling rate' };
        inistr       = { '256' };
        result       = inputdlg2( promptstr, 'Import BCI2000 data -- pop_loadbci()', 1,  inistr, 'pop_loadbci');
        if length(result) == 0 return; end
        srate   = eval( result{1} );
    end

    % import data
    % -----------
    EEG = eeg_emptyset;
    fprintf('Pop_loadbci: importing BCI file...\n');
    try
        % try to read as matlab
        % ---------------------
        bci = load( filename, '-mat');
        allfields = fieldnames(bci);
        allfields = setdiff_bc(allfields, 'signal');
        for index = 1:size(bci.signal,2)
            chanlabels{index} = [ 'C' int2str(index) ];
        end
        for index = 1:length(allfields)
            bci.signal(:,end+1) = getfield(bci, allfields{index});
        end
        EEG.chanlocs = struct('labels', { chanlabels{:} allfields{:} });
        EEG.data     = bci.signal';
        EEG.nbchan   = size(EEG.data, 1);
        EEG.pnts     = size(EEG.data, 2);
        EEG.trials   = 1;
        EEG.srate    = srate;
        EEG.comments = [ 'Original file: ' filename ];
        EEG = eeg_checkset(EEG);
        return;
    catch
        % get file names
        % --------------
        fields = loadtxt(filename, 'nlines', 1, 'verbose', 'off');
        if length(fields) > 300
            error('Not a BCI ASCII file');
        end

        % read data
        % ---------
        fid = fopen(filename, 'r');
        allcollumns = fgetl(fid);
        tmpdata = fscanf(fid, '%d', Inf);
        tmpdata = reshape(tmpdata, length(fields), length(tmpdata)/length(fields));

        EEG.data = tmpdata;
        EEG.chanlocs = struct('labels', fields);
        EEG.nbchan = size(EEG.data, 1);
        EEG.pnts   = size(EEG.data, 2);
        EEG.trials = 1;
        EEG.srate  = srate;
        EEG = eeg_checkset(EEG);
        return;

        % data channel range
        % ------------------
        indices = strmatch('ch', fields);

        bci = [];
        for index = setdiff_bc(1:length(fields), indices)
            bci = setfield(bci, fields{index}, tmpdata(index,:));
        end
        bci.signal = tmpdata(indices,:);
    end

    % ask for which event to import
    % -----------------------------
    geom = {[0.7 0.7 0.7]};
    uilist = { { 'style' 'text' 'string' 'State name' 'fontweight' 'bold' } ...
               { 'style' 'text' 'string' '    Import' 'fontweight' 'bold'  } ...
               { 'style' 'text' 'string' 'Type of'  'fontweight' 'bold' } };
    allfields = setdiff_bc(fieldnames(bci), 'signal');
    latencyfields = { '-----' };
    for index = 1:length(allfields)
        if ~isempty(findstr( lower(allfields{index}), 'time'))
            latencyfields{end+1} = allfields{index};
        end
    end
    for index = 1:length(allfields)
        geom   = { geom{:} [1.3 0.3 0.3 0.3 1] };
        uilist{end+1} = { 'style' 'text' 'string' allfields{index} };
        if ~isempty(findstr( lower(allfields{index}), 'time'))
            uilist{end+1} = { 'style' 'checkbox' 'value' 0 };
            uilist{end+1} = { };
            uilist{end+1} = { };
            uilist{end+1} = { };
        else
            uilist{end+1} = { 'style' 'checkbox' };
            uilist{end+1} = { };
            uilist{end+1} = { };
            uilist{end+1} = { 'style' 'listbox' 'string' strvcat(latencyfields) };
        end
    end
    geom   = { geom{:} [1] [0.08 1] };
    uilist{end+1} = { };
    uilist{end+1} = { 'style' 'checkbox' 'value' 0 };
    uilist{end+1} = { 'style' 'text' 'string' 'Attempt to adjust event latencies using sourcetime?' };

    result = inputgui( geom, uilist, 'pophelp(''pop_loadbci'')', 'Import BCI2000 data files - pop_loadbci()');
    if isempty(result), return; end

    % convert results to command line input
    % -------------------------------------
    listimport = {};
    count = 1;
    for index = 1:length(allfields)
        if ~isempty(findstr( lower(allfields{index}), 'time'))
            if result{count}, listimport{end+1} = 'event'; listimport{end+1} = { allfields{index} }; end
            count = count+1;
        else
            if result{count}
                if result{count+1} ~= 1
                    listimport{end+1} = 'event'; listimport{end+1} = { allfields{index}  allfields{result{count+1}-1} };
                else
                    listimport{end+1} = 'event'; listimport{end+1} = { allfields{index} };
                end
            end
            count = count+2;
        end
    end
    if result{end}, adjust = 1; else adjust = 0; end

    % decode command line input
    % -------------------------
    count = 1;
    for index = 2:2:length(listimport)
        tmpindmatch = strmatch(listimport{index}{1}, allfields, 'exact');
        if ~isempty(tmpindmatch), indeximport(count) = tmpindmatch;
        else error(['State ''' listimport{index}{1} ''' not found']);
        end
        if length( listimport{index} ) > 1
            tmpindmatch = strmatch(listimport{index}{2}, allfields, 'exact');
            if ~isempty(tmpindmatch), corresp(count) = tmpindmatch;
            else error(['State ''' listimport{index}{2} ''' not found']);
            end
        else
            corresp(count) = 0;
        end
        count = count+1;
    end

    % find block size
    % ---------------
    tmpevent = find( diff(getfield(bci, 'SourceTime')) ~= 0);
    diffevent = tmpevent(2:end)-tmpevent(1:end-1);
    blocksize = unique_bc(diffevent);
    if length(blocksize) > 1, error('Error in determining block size');
    else                      fprintf('Blocksize: %d\n', blocksize);
    end

    % find types
    % ----------
    tmpcorresp      = find(corresp);
    indexcorresp    = corresp(tmpcorresp);
    indexcorrespval = indeximport(tmpcorresp);
    if length(tmpcorresp) ~= length(intersect(corresp, indeximport))
        disp('Warning: correspondence problem, some information will be lost');
    end

    % remove type from latency array
    % ------------------------------
    indeximport(tmpcorresp) = [];
    if adjust
        fprintf('Latency of event will be adjusted\n');
    else fprintf('WARNING: Latency of event will not be adjusted (latency uncertainty %2.1f ms)\n', ...
                 blocksize/srate*1000);
    end

	% process events
	% --------------
	fprintf('Pop_loadbci: importing events...\n');
	counte = 1; % event counter
	events(10000).latency = 0;
    indexsource =  strmatch('sourcetime', lower( allfields ), 'exact' );
    sourcetime  = getfield(bci, allfields{ indexsource });
	for index = 1:length(indeximport)
        tmpdata  = getfield(bci, allfields{indeximport(index)});
		tmpevent = find( diff(tmpdata) > 0);
		tmpevent = tmpevent+1;
        tmpcorresp = find(indexcorresp == indeximport(index));

		for tmpi = tmpevent'
            if ~isempty(tmpcorresp)
                events(counte).type = allfields{ indexcorrespval(tmpcorresp) };
            else
                events(counte).type = allfields{ indeximport(index) };
            end
            if adjust
                baselatency = sourcetime(tmpi); % note that this is the first bin a block
                realtmpi    = tmpi+blocksize;   % jump to the end of the block+1
                if curlatency < baselatency, curlatency = curlatency+65536; end; % in ms
                events(counte).latency = realtmpi+(curlatency-baselatency)/1000*srate;
                % (curlatency-baselatency)/1000*srate
                % there is still a potentially large error between baselatency <-> realtmpi
            else
                events(counte).latency = tmpi+(blocksize-1)/2;
            end
            counte = counte+1;
		end
	end;

	EEG.data = bci.signal';
	EEG.nbchan = size(EEG.data, 1);
	EEG.pnts   = size(EEG.data, 2);
	EEG.trials = 1;
	EEG.srate  = srate;
	EEG = eeg_checkset(EEG);
	EEG.event = events(1:counte-1);
	EEG = eeg_checkset( EEG, 'eventconsistency' );
	EEG = pop_editeventvals( EEG, 'sort', { 'latency', [0] } );

% $$$
% $$$
% $$$
% $$$
% $$$
% $$$
% $$$
% $$$
% $$$
% $$$
% $$$
% $$$
% $$$
% $$$
% $$$ dsffd
% $$$
% $$$ 	% find electrode indices
% $$$ 	% ----------------------
% $$$ 	indices = [];
% $$$ 	for index = 1:length(colnames)
% $$$ 		if strcmp( colnames{index}(1:2), 'ch')
% $$$ 			indices = [indices index ];
% $$$ 		end
% $$$ 	end
% $$$
% $$$     EEG.data = tmpdata(indices,:);
% $$$     EEG.nbchan = size(EEG.data, 1);
% $$$     EEG.srate  = srate;
% $$$     try
% $$$         eventindices = setdiff_bc(1:length(colnames), indices);
% $$$         ISIind = eventindices(3 + 9);
% $$$         eventindices(3 + [ 1 2 3 4 7 8 9 10 11 12]) = [];
% $$$         eventindices(1:3) = []; % suppress these event
% $$$
% $$$         % add the trial number
% $$$         % --------------------
% $$$         tmptrial = find( diff(tmpdata(ISIind, :)) ~= 0);
% $$$         tmptrial = tmptrial+1;
% $$$
% $$$         % process events
% $$$         % --------------
% $$$         fprintf('Pop_loadbci: importing events...\n');
% $$$         counte = 1; % event counter
% $$$         events(10000).latency = 0;
% $$$         for index = eventindices
% $$$             counttrial = 1;
% $$$             tmpevent = find( diff(tmpdata(index, :)) ~= 0);
% $$$             tmpevent = tmpevent+1;
% $$$             for tmpi = tmpevent
% $$$                 if tmpdata(index, tmpi)
% $$$                     events(counte).type    = [ colnames{index} int2str(tmpdata(index, tmpi)) ];
% $$$                     events(counte).latency = tmpi;
% $$$                     %events(counte).value   = tmpdata(index, tmpi);
% $$$                     %while tmpi > tmptrial(counttrial) & counttrial < length(tmptrial)
% $$$                     %	counttrial = counttrial+1;
% $$$                     %end
% $$$                     %events(counte).trial = counttrial;
% $$$                     counte = counte+1;
% $$$                     %if mod(counte, 100) == 0, fprintf('%d ', counte); end
% $$$                 end
% $$$             end
% $$$         end
% $$$
% $$$         % add up or down events
% $$$         % ---------------------
% $$$         EEG = eeg_checkset(EEG);
% $$$         EEG.event = events(1:counte-1);
% $$$         EEG = pop_editeventvals( EEG, 'sort', { 'latency', [0] } );
% $$$         for index=1:length(EEG.event)
% $$$             if strcmp(EEG.event(index).type(1:6), 'Target')
% $$$                 targetcode = str2num(EEG.event(index).type(end));
% $$$                 if targetcode == 1
% $$$                     EEG.event(index).type = 'toptarget';
% $$$                 else
% $$$                     EEG.event(index).type = 'bottomtarget';
% $$$                 end
% $$$             else
% $$$                 if strcmp(EEG.event(index).type(1:6), 'Result')
% $$$                     resultcode = str2num(EEG.event(index).type(end));
% $$$                     if resultcode == 1
% $$$                         EEG.event(index).type = 'topresp';
% $$$                     else
% $$$                         EEG.event(index).type = 'bottomresp';
% $$$                     end
% $$$                     EEG.event(end+1).latency = EEG.event(index).latency;
% $$$                     if (resultcode == targetcode)
% $$$                         EEG.event(end).type = 'correct';
% $$$                     else
% $$$                         EEG.event(end).type = 'miss';
% $$$                     end
% $$$                 end
% $$$             end
% $$$         end
% $$$         EEG = pop_editeventvals( EEG, 'sort', { 'latency', [0] } );
% $$$         %EEG.data = tmpdata([72 73 75],:);
% $$$     catch, disp('Failed to import data events');
% $$$     end
% $$$
% $$$ command = sprintf('EEG = pop_loadbci(''%s'', %f);',filename, srate);
% $$$ return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
classdef InertialAeroAnglesPolySteeringModel < AbstractAnglePolySteeringModel
    %InertialAeroAnglesPolySteeringModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        bankModel(1,1) PolynominalModel = PolynominalModel(0,0,0,0);
        aoAModel(1,1) PolynominalModel = PolynominalModel(0,0,0,0);
        slipModel(1,1) PolynominalModel = PolynominalModel(0,0,0,0);

        bankContinuity(1,1) logical = true;
        aoAContinuity(1,1) logical = true;
        slipContinuity(1,1) logical = true;
    end

    methods
        function R_body_2_inertial = getBody2InertialDcmAtTime(obj, ut, rVect, vVect, bodyInfo)
            bankAng = obj.bankModel.getValueAtTime(ut);
            angOfAttack = obj.aoAModel.getValueAtTime(ut);
            angOfSideslip = obj.slipModel.getValueAtTime(ut);

            baseFrame = bodyInfo.getBodyCenteredInertialFrame();
            [~, ~, ~, R_body_2_inertial] = computeInertialBodyAxesFromFrameAeroAngles(ut, rVect, vVect, bodyInfo, bankAng, angOfAttack, angOfSideslip, baseFrame);
        end

        function [angleModel, continuity] = getAngleNModel(obj, n)
            angleModel = PolynominalModel.empty(1,0);

            switch n
                case 1
                    angleModel = obj.bankModel;
                    continuity = obj.bankContinuity;
                case 2
                    angleModel = obj.aoAModel;
                    continuity = obj.aoAContinuity;
                case 3
                    angleModel = obj.slipModel;
                    continuity = obj.slipContinuity;
            end
        end

        function t0 = getT0(obj)
            t0 = obj.slipModel.t0;
        end

        function setT0(obj, newT0)
            obj.bankModel.t0 = newT0;
            obj.aoAModel.t0 = newT0;
            obj.slipModel.t0 = newT0;
        end

        function setConstTerms(obj, bankConst, aoaConst, slipConst)
            obj.bankModel.constTerm = bankConst;
            obj.aoAModel.constTerm = aoaConst;
            obj.slipModel.constTerm = slipConst;
        end

        function setLinearTerms(obj, bank, aoa, slip)
            obj.bankModel.linearTerm = bank;
            obj.aoAModel.linearTerm = aoa;
            obj.slipModel.linearTerm = slip;
        end

        function setAccelTerms(obj, bank, aoa, slip)
            obj.bankModel.accelTerm = bank;
            obj.aoAModel.accelTerm = aoa;
            obj.slipModel.accelTerm = slip;
        end

        function setTimeOffsets(obj, timeOffset)
            obj.bankModel.tOffset = timeOffset;
            obj.aoAModel.tOffset = timeOffset;
            obj.slipModel.tOffset = timeOffset;
        end

        function [angle1Cont, angle2Cont, angle3Cont] = getContinuityTerms(obj)
            angle1Cont = obj.bankContinuity;
            angle2Cont = obj.aoAContinuity;
            angle3Cont = obj.slipContinuity;
        end

        function setContinuityTerms(obj, angle1Cont, angle2Cont, angle3Cont)
            obj.bankContinuity = angle1Cont;
            obj.aoAContinuity = angle2Cont;
            obj.slipContinuity = angle3Cont;
        end

        function setConstsFromDcmAndContinuitySettings(obj, dcm, ut, rVect, vVect, bodyInfo)
            if(obj.bankContinuity || obj.aoAContinuity || obj.slipContinuity)
                [bankAng,angOfAttack,angOfSideslip] = computeInertialAeroAnglesFromBodyAxes(ut, rVect, vVect, bodyInfo, dcm(:,1), dcm(:,2), dcm(:,3));

                if(obj.bankContinuity)
                    obj.bankModel.constTerm = bankAng;
                end

                if(obj.aoAContinuity)
                    obj.aoAModel.constTerm = angOfAttack;
                end

                if(obj.slipContinuity)
                    obj.slipModel.constTerm = angOfSideslip;
                end
            end
        end

        function setInitialAttitudeFromState(obj, stateLogEntry, tOffsetDelta)
            t0 = stateLogEntry.time;
            obj.setT0(t0);

            obj.bankModel.tOffset = obj.bankModel.tOffset + tOffsetDelta;
            obj.aoAModel.tOffset = obj.aoAModel.tOffset + tOffsetDelta;
            obj.slipModel.tOffset = obj.slipModel.tOffset + tOffsetDelta;
        end

        function [angle1Name, angle2Name, angle3Name] = getAngleNames(~)
            angle1Name = 'Bank Angle';
            angle2Name = 'Angle of Attack';
            angle3Name = 'Side Slip Angle';
        end

        function newSteeringModel = deepCopy(obj)
            newSteeringModel = InertialAeroAnglesPolySteeringModel(obj.bankModel.deepCopy(), obj.aoAModel.deepCopy(), obj.slipModel.deepCopy());
            newSteeringModel.bankContinuity = obj.bankContinuity;
            newSteeringModel.aoAContinuity = obj.aoAContinuity;
            newSteeringModel.slipContinuity = obj.slipContinuity;
        end

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

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

        function [addActionTf, steeringModel] = openEditSteeringModelUI(obj, lv, useContinuity)
            output = AppDesignerGUIOutput({false, obj});
            lvd_EditActionSetSteeringModelGUI_App(obj, lv, useContinuity, output);
            addActionTf = output.output{1};
            steeringModel = output.output{2};
        end
    end

    methods(Access=private)
        function obj = InertialAeroAnglesPolySteeringModel(bankModel, aoAModel, slipModel)
            obj.bankModel = bankModel;
            obj.aoAModel = aoAModel;
            obj.slipModel = slipModel;
        end
    end

    methods(Static)
        function model = getDefaultSteeringModel()
            bankModel = PolynominalModel(0,0,0,0);
            aoAModel = PolynominalModel(0,0,0,0);
            slipModel = PolynominalModel(0,0,0,0);

            model = InertialAeroAnglesPolySteeringModel(bankModel, aoAModel, slipModel);
        end

        function typeStr = getTypeNameStr()
            typeStr = SteeringModelEnum.InertialAeroAnglesPoly.nameStr;
        end
    end
end


text:
function y = testfunctions(varargin)
% TESTFUNCTIONS   Integration test functions of A. Genz
%    Y = TESTFUNCTIONS(X1, X2, ..., XD, TYPE, C, W) Evaluates test
%    function TYPE at the point (X1,...,XN). C and W are arrays of
%    constants defining the test function (see below). The test
%    functions are defined according to A. Genz: A package for
%    testing multiple integration subroutines, in Numerical
%    Integration, P. Keast and G. Fairweather (Eds.), D. Riedel,
%    pp. 337-340, 1987. The functions are defined on [0 1]^d.
%
% TYPE = 'oscillatory'   | TYPE = 1:
% f(x) = cos(2*pi*w_1 + sum_{i=1}^d ( c_i * x_i ) )
%
% TYPE = 'product peak'  | TYPE = 2:
% f(x) = prod_{i=1}^d ( c_i^{-2} + (x_i - w_i)^2 )^{-1}
%
% TYPE = 'corner peak'   | TYPE = 3:
% f(x) = ( 1 + sum_{i=1}^d (c_i * x_i) )^(-(d+1))
%
% TYPE = 'gaussian'      | TYPE = 4:
% f(x) = exp( - sum_{i=1}^d c_i^2 * (x_i - w_i)^2 )
%
% TYPE = 'continuous'    | TYPE = 5:
% f(x) = exp( - sum_{i=1}^d c_i * abs(x_i - w_i) )
%
% TYPE = 'discontinuous' | TYPE = 6:
% f(x) = { 0                            , if x_1>w_1 or x_2>w_2,
%        { exp( sum_{i=1}^d c_i * x_i ) , otherwise
%
% With the parameters c = (c_1, ..., c_d) and w = (w_1, ..., w_d).
% d denotes the dimension of the function.
%
% Examples:
%    testfunctions(0.5, 0.5, 'product peak', [2, 5.25], [0.2, 0.7])
%    testfunctions(0.5, 0.5, 2, [2, 5.25], [0.2, 0.7])
%
%    x = linspace(0,1,20);
%    [X,Y] = meshgrid(x,x);
%    surf(X,Y, ...
%      testfunctions(X, Y, 'product peak', [2, 5.25], [0.2 0.7]));

% Author : Andreas Klimke, Universitaet Stuttgart
% Version: 1.0
% Date   : August 2, 2003

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

type = varargin{end-2};
c = varargin{end-1};
w = varargin{end};

if isa(type, 'char')
	fnames = {'oscillatory', 'product peak', 'corner peak', 'gaussian', ...
						'continuous', 'discontinuous'};
	for k = 1:length(fnames)
		if strcmp(fnames{k},type)
			type = k;
			break
		end
	end
end

d = length(c);

switch type
 case 1  % oscillatory
	temp = 2*pi*w(1);
	for i = 1:d
		temp = temp + c(i).*varargin{i};
	end
	y = cos(temp);

 case 2  % product peak
	temp = 1;
	for i = 1:d
		temp = temp .* (c(i)^(-2)+(varargin{i}-w(i)).^2);
	end
	y = 1./temp;

 case 3  % corner peak
	temp = 1;
	for i = 1:d
		temp = temp + c(i).*varargin{i};
	end
	y = temp .^ (-(d+1));

 case 4  % gaussian
	temp = 0;
	for i = 1:d
		temp = temp + c(i)^2 .* (varargin{i} - w(i)).^2;
	end
	y = exp(-temp);

 case 5  % continuous
	temp = 0;
	for i = 1:d
		temp = temp + c(i) .* abs(varargin{i} - w(i));
	end
	y = exp(-temp);

 case 6  % discontinuous
	temp = 0;
	if d >= 2
		mask = varargin{1} > w(1) | varargin{2} > w(2);
	else
		mask = varargin{1} > w(1);
	end
	for i = 1:d
		temp = temp + (c(i) .* varargin{i});
	end
	y = exp(temp) .* (~mask);
end

author: johannesgerer

repo: jburkardt-m


text:
clear
clc

pro = pro_Create();

% pro = pro_AddInput(pro, @(N)pdf_Uniform(N, [-pi pi]), 'X1');
% pro = pro_AddInput(pro, @(N)pdf_Uniform(N, [-pi pi]), 'X2');
% pro = pro_AddInput(pro, @(N)pdf_Uniform(N, [-pi pi]), 'X3');

pro = pro_AddInput(pro, @()pdf_Sobol([-pi pi]), 'X1');
pro = pro_AddInput(pro, @()pdf_Sobol([-pi pi]), 'X2');
pro = pro_AddInput(pro, @()pdf_Sobol([-pi pi]), 'X3');

pro = pro_SetModel(pro, @(x)TestModel2(x), 'model');

pro.N = 10000;

pro = GSA_Init(pro);
[S eS pro] = GSA_GetSy(pro, {1});
[Stot eStot pro] = GSA_GetTotalSy(pro, {1});

author: Sable

repo: mcbench-benchmarks


text:
classdef nncycle < nntest
  methods (Test)
    function testCycle(test)
      % diamond topology
      a = Input() ;
      b = sqrt(a) ;
      c = abs(a) ;
      d = b + c ;
      e = exp(d) ;

      % introduce cycles in the DAG
      try
        b.inputs{1} = e ;
      catch err
        assert(strcmp(err.identifier, 'MatConvNet:CycleCheckFailed'), ...
          'Failed to catch cycle in DAG.') ;
      end

      try
        d.inputs{2} = d ;
      catch err
        assert(strcmp(err.identifier, 'MatConvNet:CycleCheckFailed'), ...
          'Failed to catch cycle in DAG.') ;
      end
    end
  end
end

repo: MFT


text:
%FEATSELB Trainable mapping for backward feature selection
%
%   [W,R] = FEATSELB(A,CRIT,K,T)
%   [W,R] = A*FEATSELB([],CRIT,K,T)
%   [W,R] = A*FEATSELB(CRIT,K,T)
%   [W,R] = FEATSELB(A,CRIT,K,N)
%   [W,R] = A*FEATSELB([],CRIT,K,N)
%   [W,R] = A*FEATSELB(CRIT,K,N)
%
% INPUT
%   A     Dataset
%   CRIT  String name of the criterion or untrained mapping
%         (optional; default: 'NN', i.e. 1-Nearest Neighbor error)
%   K     Number of features to select
%         (optional; default: return optimally ordered set of all features)
%   T     Tuning set (optional)
%   N     Number of cross-validations
%
% OUTPUT
%   W     Output feature selection mapping
%   R     Matrix with step-by-step results of the selection
%
% DESCRIPTION
% Backward selection of K features using the dataset A. CRIT sets the
% criterion used by the feature evaluation routine FEATEVAL. If the
% dataset T is given, it is used as test set for FEATEVAL. Alternatvely a
% a number of cross-validation N may be supplied. For K = 0, the optimal
% feature set (corresponding to the maximum value of FEATEVAL) is returned.
% The result W can be used for selecting features by B*W. In this case,
% features are ranked optimally.
% The selected features are stored in W.DATA and can be found by +W.
% In R, the search is reported step by step as:
%
% 	R(:,1) : number of features
% 	R(:,2) : criterion value
% 	R(:,3) : added / deleted feature
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, FEATEVAL, FEATSELLR, FEATSEL,
% FEATSELO, FEATSELF, FEATSELI, FEATSELP, FEATSELM

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: featselb.m,v 1.6 2008/07/03 09:08:43 duin Exp $

function [w,r] = featselb(varargin)

  varargin = shiftargin(varargin,{'char','prmapping'});
  argin = setdefaults(varargin,[],'NN',0,[],[]);
  if mapping_task(argin,'definition')
    w = define_mapping(argin,'untrained','Backward FeatSel');
    return
  end

  [a,crit,ksel,t,fid] = deal(argin{:});
	[w,r] = featsellr(a,crit,ksel,0,1,t,fid);
  %DXD This is a patch: when the number of features has to be
  %optimized, and all features seem useful, when the list of
  %features is not reshuffled to reflect the relative importance of
  %the features:
  % (Obviously, this should be fixed in featsellr, but I don't
  % understand what is happening in there)
  dim = size(a,2);
  if (ksel==0) & (length(getdata(w))==dim)
    rr = -r(:,3); rr(1) = [];
    rr = [setdiff((1:dim)',rr) rr(end:-1:1)'];
    w = setdata(w,rr);
  end

return

author: marianux

repo: ecg-kit

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


text:
function x = tuple_next_fast ( m, n, rank )

%*****************************************************************************80
%
%% TUPLE_NEXT_FAST computes the next element of a tuple space, "fast".
%
%  Discussion:
%
%    The elements are N vectors.  Each entry is constrained to lie
%    between 1 and M.  The elements are produced one at a time.
%    The first element is
%      (1,1,...,1)
%    and the last element is
%      (M,M,...,M)
%    Intermediate elements are produced in lexicographic order.
%
%    This code was written as a possibly faster version of TUPLE_NEXT.
%
%  Example:
%
%    N = 2,
%    M = 3
%
%    INPUT        OUTPUT
%    -------      -------
%    Rank          X
%    ----          ----
%   -1            -1 -1
%    0             1  1
%    1             1  2
%    2             1  3
%    3             2  1
%    4             2  2
%    5             2  3
%    6             3  1
%    7             3  2
%    8             3  3
%    9             1  1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the maximum entry.
%
%    Input, integer N, the number of components.
%
%    Input, integer RANK, indicates the rank of the tuples.
%    However, if RANK < 0, this indicates that this is the first
%    call for the given values of (M,N).  Initialization is done,
%    and X is set to a dummy value.
%
%    Output, integer X(N), the next tuple of the given rank.
%
  global tuple_next_fast_BASE

  if ( rank < 0 )

    tuple_next_fast_BASE(n) = 1;
    for i = n-1 : -1 : 1
      tuple_next_fast_BASE(i) = tuple_next_fast_BASE(i+1) * m;
    end

    x(1:n) = -1;

  else

    x(1:n) = mod ( floor ( rank ./ tuple_next_fast_BASE(1:n) ), m ) + 1;

  end

  return
end


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


text:
%File Description:
%  Unknown parameters in magnetometer error model are calculated using the
%  lm optimization algorithm.
close all
clc
clear
load MagRaw.mat
%magnetometer calibration parameters in the Pixhawk autopilot which can
 %be obtained from QGC
CAL_MAG_SCALE = [1, 1, 1]';
CAL_MAG_OFF = [0.064, 0.014, -0.053]';
MagRaw = (mag + CAL_MAG_OFF)./CAL_MAG_SCALE;  % raw data

m = length(MagRaw);
MagSum = 0;
for k = 1 : m
    MagSum = MagSum + norm(MagRaw(:, k));
end
MagAver = MagSum/m;  % estimated magnetic field strength
Vdata = MagRaw/MagAver;  % normalization

y_dat = ones(m, 1);
p0 = [1 1 1 0 0 0]';
p_init = [1 1 1 0.01 0.01 0.01]';  % initial value of unknown parameters

y_raw = calFunc(Vdata, p0);  %magnetic field strength measured by magnetometer
y_raw = y_raw(:);
r_raw = y_dat - y_raw;
p_fit = lm('calFunc', p_init, Vdata, y_dat, 0.001);
y_lm = calFunc(Vdata, p_fit);  % magnetic field strength measured by calibrated magnetometer
y_lm = y_lm(:);
r_lm = y_dat - y_lm;
y_px4 = calFunc(mag/MagAver, p0);  % magnetic field strength measured by calibrated magnetometer of PX4
y_px4 = y_px4(:);
r_px4 = y_dat - y_px4;

kx = p_fit(1);
ky = p_fit(2);
kz = p_fit(3);
bx = p_fit(4);
by = p_fit(5);
bz = p_fit(6);

Km = [kx 0 0;0 ky 0;0 0 kz]
bm = [bx by bz]'

figure
bar([r_raw'*r_raw, r_lm'*r_lm, r_px4'*r_px4])
grid on;
set(gca, 'XTickLabel', {'raw','lm','px4'});
ylabel('Index');

t=1:m;
figure
title('Magnetometer Calibration')
plot(t, r_raw, '-.', t, r_lm, '-')
legend('Uncalibrated','Calibrated-LM')
xlabel('Measurement sampling number')
ylabel('Residual')

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function out = saveobj(obj)
%SAVEOBJ (overloaded)
warning('YALMIP objects cannot be saved in binary format. You will run into troubles if you try to load this file later. You should clear all YALMIP objects first, or avoid having any YALMIP objects in the list of variables which you save.');
out = [];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef PlaneStressTransformerFactory < handle

    properties (Access = private)
        PSTransformer
        tensor
    end

    methods (Access = public)

        function psTransformer = create(obj,Tensor)
            obj.init(Tensor)
            obj.createPlaneStressTransformer()
            psTransformer = obj.getPlaneStressTransformer();
      end
    end

    methods (Access = private)

        function init(obj,tensor)
            obj.tensor = tensor;
        end

        function createPlaneStressTransformer(obj)
            t = obj.tensor;
            ps = [];
            if obj.isVoigt()
                if obj.isStiffnessTensor()
                    if obj.isInvertible()
                        ps = PST4VoigtFourthOrderTensorNumerically(t.getValue());
                    else
                        ps = PST4VoigtFourthOrderTensorSymbolically(t.getValue());
                    end

                elseif obj.isSecondOrder()
                    ps = PST4VoigtSecondOrderTensor(t);
                end

            elseif obj.isTensor()
                if obj.isSecondOrder()
                    ps =  PlaneStressTransformerForSecondOrderTensor(t);
                end

            end

            if isempty(ps)
                error('Not admitted object to make it Plane Stress')
            else
                obj.PSTransformer = ps;
            end

        end

        function PS = getPlaneStressTransformer(obj)
            PS = obj.PSTransformer;
        end

        function itIs = isVoigt(obj)
            itIs = strcmp(obj.tensor.getRepresentation(),'voigt');
        end

        function itIs = isTensor(obj)
            itIs = strcmp(obj.tensor.getRepresentation(),'tensor');
        end

        function ItIs = isSecondOrder(obj)
            ItIs = strcmp(obj.tensor.getOrder(),'second');
        end

         function itIs = isStiffnessTensor(obj)
            itIs = strcmp(obj.tensor.getFieldName,'stiffness');
        end

        function itIs = isInvertible(obj)
            t = obj.tensor.getValue();
            if obj.isSymbolic(t)
                itIs = isAlways(det(t) == 0);
            else
                itIs = min(abs(eig(t))) > 1e-13;
            end
        end

    end

    methods (Access = private, Static)
        function itIs = isSymbolic(a)
            itIs = isa(a,'sym');
        end
    end

end


text:
classdef BT9 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP with bias feature

%------------------------------- Reference --------------------------------
% H. Li, Q. Zhang, and J. Deng, Biased multiobjective optimization and
% decomposition algorithm, IEEE Transactions on Cybernetics, 2017, 47(1):
% 52-66.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 3;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            [N,D] = size(X);
            I1    = 3 : 3 : D;
            I2    = 4 : 3 : D;
            I3    = 5 : 3 : D;
            Y     = X - sin(repmat(1:D,N,1)*pi/2/D);
            PopObj(:,1) = cos(0.5*X(:,1)*pi).*cos(0.5*X(:,2)*pi) + sum(Y(:,I1).^2+(1-exp(-Y(:,I1).^2/1e-9))/5,2);
            PopObj(:,2) = cos(0.5*X(:,1)*pi).*sin(0.5*X(:,2)*pi) + sum(Y(:,I2).^2+(1-exp(-Y(:,I2).^2/1e-9))/5,2);
            PopObj(:,3) = sin(0.5*X(:,1)*pi)                     + sum(Y(:,I3).^2+(1-exp(-Y(:,I3).^2/1e-9))/5,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,3);
            R = R./repmat(sqrt(sum(R.^2,2)),1,3);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            a = linspace(0,pi/2,10)';
            R = {sin(a)*cos(a'),sin(a)*sin(a'),cos(a)*ones(size(a'))};
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO