text:
function [y, stopCondition, sigma, sigma0, t] ...
    = cons_smacof_pip(dx, y, isFree, bnd, w, con, smacof_opts, scip_opts)
% CONS_SMACOF_PIP  SMACOF algorithm with polynomial constraints (PIP file
% format).
%
% Scaling by MAjorizing a COnvex Function (SMACOF) is an iterative solution
% to the Multidimensional Scaling (MDS) problem (see de Leeuw and Mair [1]
% for an up to date review).
%
% The classic implementation of SMACOF uses an iterative algorithm that
% relies on the Guttman transform update. Our function here implements a
% different approach, where the Guttman transform update is replaced by
% solving a Quadratic Program (QP) (as proposed by Dwyer et al [2]) with
% polynomial constraints.
%
% In this implementation, we use the SCIP binary to solve the constrained
% QP. This binary can be downloaded for Linux, MacOS X and Windows from the
% Zuse Institute Berlin (ZIB) website
%
%   http://scip.zib.de/#download
%
% The objective function is
%
%   min_y 1/2 y' * H * y + f' * y
%
% subject to the constraints and bounds provided by the user.
%
% We use the PIP file format to formulate the constrained QP and pass it to
% SCIP.
%
%   http://polip.zib.de/pipformat.php
%
%
% [Y, STOPCONDITION, SIGMA, SIGMA0, T] = CONS_SMACOF_PIP(D, Y0, ISFREE, BND, [], CON)
%
%   D is an (N, N)-distance matrix, with distances between the points in an
%   N-point configuration. D can be full or sparse. D(i,j)=0 means that
%   vertices i and j are not directly connected.
%
%   Y0 is an initial guess of the solution, given as an (N, P)-matrix,
%   where P is the dimensionality of the output points. Currently, P must
%   be either 2 or 3.
%
%   BND is a cell array with the variable bounds in PIP format, and cannot
%   be empty (otherwise SCIP doesn't return a solution). E.g.
%
%      BND = {'Bounds', ' -1 <= x1 <= 4', ' -1 <= y1 <= 2.5', ...
%             ' -1 <= x2 <= 4', ' -1 <= y2 <= 2.5'};
%
%   CON is a cell array with the problem constraints in PIP format, and
%   cannot be empty (otherwise SCIP doesn't return a solution). E.g.
%
%      CON = {'Subject to', ...
%            ' c1: -0.5 x6 y7 +0.5 x3 y7 +0.5 x7 y6 -0.5 x3 y6 -0.5 x7 y3 +0.5 x6 y3 >= 0.1'};
%
%   Y is the best solution computed by the algorithm within all iterations.
%   Y is a point configuration with the same size as Y0. If the algorithm
%   cannot find any valid solution (e.g. because we set a too short time
%   limit, or because a valid solution doesn't exist, Y will be an array of
%   NaNs).
%
%   STOPCONDITION is a cell array with a string for each stop condition
%   that made the algorithm stop at the last iteration.
%
%   SIGMA is a vector with the stress value at each iteration. Weighted
%   stress is given as
%
%     SIGMA = \sum_{i<j} W_ij (D_ij - DY_ij)^2
%
%   where W is a weight matrix the same size as D. W_ij = 0 means that the
%   distance between points i and j does not affect the stress measure.
%
%   If the algorithm could not find any solution, SIGMA is a single scalar
%   with value = Inf. Note that even if the last SIGMA value is Inf, this
%   doesn't mean that the algorithm hasn't found a valid solution in a
%   previous iteration.
%
%   SIGMA0 is a scalar with the stress value of the initial guess Y0. Note
%   that we have no guarantee that Y0 is a valid solution, so be careful
%   when comparing SIGMA0 to SIGMA. In particular, SIGMA0 may be smaller
%   than any SIGMA, but may not be a valid solution according to the
%   constraints.
%
%   T is a vector with the time between the beginning of the algorithm and
%   each iteration. Units in seconds.
%
% Y = CONS_SMACOF_PIP(..., SMACOF_OPTS, SCIP_OPTS)
%
%   SMACOF_OPTS is a struct with parameters to tweak the SMACOF algorithm.
%
%     'MaxIter': (default = 0) Maximum number of majorization iterations we
%                allow the optimisation algorithm.
%
%     'Epsilon': (default = Inf) The algorithm will stop if
%                (SIGMA(I+1)-SIGMA(I))/SIGMA(I) < OPTS.Epsilon.
%
%     'Display': (default = 'off') Do not display any internal information.
%                'iter': display internal information at every iteration.
%
%     'TolFun':  (default = 1e-12) Termination tolerance of the stress
%                value.
%
%   SCIP_OPTS is a struct with parameters to tweak the SCIP algorithm.
%
%     'scipbin': (defaults = 'scip.linux.x86_64.gnu.opt.spx' (Linux),
%                            'scip.darwin.x86_64.gnu.opt.spx' (Mac),
%                            'scip.mingw.x86_64.intel.opt.spx.exe' (Win))
%                Name of the SCIP binary/executable. This binary should be
%                available in the system path. E.g. place it in
%                gerardus/programs. The binaries/executable can be
%                downloaded from http://scip.zib.de/#download. You may need
%                to rename them so that they fit the expected filenames.
%
%     'display_verblevel': (default 4) verbosity level of output (0: SCIP
%                quiet mode).
%
%     'display_freq': (default 100) frequency for displaying node
%                information lines.
%
%     'limits_absgap': (default 0.0) solving stops, if the absolute
%                gap = |primalbound - dualbound| is below the given value.
%
%     'limits_gap': (default 0.0) solving stops, if the relative
%                gap = |primal - dual|/MIN(|dual|,|primal|) is below the
%                given value.
%
%     'limits_time': (default 1e+20) maximal time in seconds to run.
%
%     'limits_solutions': (default -1) solving stops, if the given number
%                of solutions were found (-1: no limit).
%
%     'lp_threads': (default 0: automatic) number of threads used for
%                solving the LP.
%
%     'numerics_feastol': (default 1e-6) feasibility tolerance for
%                constraints in SCIP. If this value is changed from the
%                default, it needs to be passed to the function that
%                creates the constraints too, e.g.
%                tri_ccqp_smacof_nofold_sph_pip().
%
%
% [1] J de Leeuw, P Mair, "Multidimensional scaling using majorization:
% SMACOF in R", Journal of Statistical Software, 31(3), 2009.
%
% [2] T. Dwyer, Y. Koren, and K. Marriott, "Drawing directed graphs using
% quadratic programming," IEEE Transactions on Visualization and Computer
% Graphics, vol. 12, no. 4, pp. 536-548, 2006.
%
% See also: cmdscale, qcqp_smacof, tri_sphparam.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2014, 2016 University of Oxford
% Version: 0.4.6
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

%% Input arguments

% check arguments
narginchk(6, 8);
nargoutchk(0, 5);

% start clock
tic

% number of points
N = size(dx, 1);

% dimensionality of the points
D = size(y, 2);
if ((D ~= 2) && (D ~= 3))
    error('Only implemented for 2D output')
end

% check inputs
if (N ~= size(dx, 2))
    error('D must be a square matrix')
end
if (N ~= size(y, 1))
    error('Y0 must have the same number of rows as D')
end

% defaults
if (nargin < 3 || isempty(isFree))
    isFree = true(N, 1);
end
Nfree = nnz(isFree);

if (isempty(w))
    % if the user doesn't provide a weight matrix, we simply assign 1 if
    % two vertices are connected, and 0 if not
    w = double(dx ~= 0);
    if (any(diag(w)) ~= 0)
        error('Assertion error: W matrix has diagonal elements that are non-zero')
    end
end
if (any(size(w) ~= [N N]))
    error('W must be a square matrix with the same size as D')
end

% if the user doesn't allow any vertices to be optimized, then we can exit
% the function
if (Nfree == 0)
    stopCondition = 'No free vertices to optimise';
    sigma = [];
    dy = dmatrix_con(dx, y);
    sigma0 = 0.5 * sum(sum(w .* (dx - dy).^2));
    t = [];
    return;
end

% SMACOF_OPTS defaults
if (nargin < 7 || isempty(smacof_opts) || ~isfield(smacof_opts, 'MaxIter'))
    smacof_opts.MaxIter = 100;
end
if (nargin < 7 || isempty(smacof_opts) || ~isfield(smacof_opts, 'Epsilon'))
    smacof_opts.Epsilon = 0;
end
if (nargin < 7 || isempty(smacof_opts) || ~isfield(smacof_opts, 'Display'))
    smacof_opts.Display = 'off';
end
if (nargin < 7 || isempty(smacof_opts) || ~isfield(smacof_opts, 'TolFun'))
    smacof_opts.TolFun = 1e-12;
end

% SCIP_OPTS
if (nargin < 8 || isempty(scip_opts) || ~isfield(scip_opts, 'scipbin'))

    % default name of the SCIP binary depending on the architecture
    if (isunix)
        SCIPBIN = 'scip.linux.x86_64.gnu.opt.spx';
    elseif (ismac)
        SCIPBIN = 'scip.darwin.x86_64.gnu.opt.spx';
    elseif (ispc)
        SCIPBIN = 'scip.mingw.x86_64.intel.opt.spx.exe';
    else
        error('Operating system not recognized: I do not know what is the default name for the SCIP binary')
    end
end

scip_opts_comm = {};
if (nargin >= 8 && ~isempty(scip_opts))

    % SCIP binary
    if (isfield(scip_opts, 'scipbin'))
        SCIPBIN = scip_opts.scipbin;
    end

    % display options
    QUIETFLAG = [];
    OUTPUTREDIR = []; % output redirecton, e.g. "> /dev/null"
    if isfield(scip_opts, 'display_verblevel')
        if (scip_opts.display_verblevel == 0)
            QUIETFLAG = ' -q ';

            % bug workaround: there's a bug in
            % scip-3.1.0.linux.x86_64.gnu.opt.spx, which causes the
            % solution to be written as an empty file when the quiet flag
            % is used. As a workaround, we disable the quiet flag in that
            % case, and instead send the output to /dev/null, as suggested
            % by Stefan Vigerske
            if strcmp(SCIPBIN, 'scip-3.1.0.linux.x86_64.gnu.opt.spx')
                QUIETFLAG = [];
                OUTPUTREDIR = ' > /dev/null';
            end

        end
        scip_opts_comm{end+1} = [' -c "set display verblevel ' num2str(scip_opts.display_verblevel) '"'];
    end

    % frequency for displaying node information lines [100]
    if (isfield(scip_opts, 'display_freq'))
        scip_opts_comm{end+1} = [' -c "set display freq ' num2str(scip_opts.display_freq) '"'];
    end

    % limits options

    % solving stops, if the absolute gap = |primalbound - dualbound| is below the given value [0.0]
    if (isfield(scip_opts, 'limits_absgap'))
        scip_opts_comm{end+1} = [' -c "set limits absgap ' num2str(scip_opts.limits_absgap) '"'];
    end

    % solving stops, if the relative gap = |primal - dual|/MIN(|dual|,|primal|) is below the given value [0.0]
    if (isfield(scip_opts, 'limits_gap'))
        scip_opts_comm{end+1} = [' -c "set limits gap ' num2str(scip_opts.limits_gap) '"'];
    end

    % maximal time in seconds to run [1e+20]
    if (isfield(scip_opts, 'limits_time'))
        scip_opts_comm{end+1} = [' -c "set limits time ' num2str(scip_opts.limits_time) '"'];
    end

    % solving stops, if the given number of solutions were found (-1: no limit) [-1]
    if (isfield(scip_opts, 'limits_solutions'))
        scip_opts_comm{end+1} = [' -c "set limits solutions ' num2str(scip_opts.limits_solutions) '"'];
    end

    % lp options

    % number of threads used for solving the LP (0: automatic)
    if (isfield(scip_opts, 'lp_threads'))
        scip_opts_comm{end+1} = [' -c "set lp advanced threads ' num2str(scip_opts.lp_threads) '"'];
    end

    % numerics options

    % feasibility tolerance for constraints in SCIP [1e-6]
    if (isfield(scip_opts, 'numerics_feastol'))
        scip_opts_comm{end+1} = [' -c "set numerics feastol ' num2str(scip_opts.numerics_feastol) '"'];
    end

end

%% Objective function: 1/2 nu' * H * nu + f' * nu

% pre-compute the weighted Laplacian matrix
V = -w;
V(1:N+1:end) = sum(w, 2);

% quadratic terms of the objective function

% upper triangular matrix terms (we assume symmetric matrix)
dx = tril(dx);

% remove pairs where both vertices are fixed, as those only contribute a
% constant to the objective function and can be ignored in the optimization
dx(~isFree, ~isFree) = 0;

Nterms = nnz(dx);
objfunq = cell(1, Nterms + Nfree);
count = 1;
for idx = find(dx)'

    % decode dx matrix index into [I, J] vertices (we swap I,J so that we
    % get I as smaller index, and J larger index
    [J, I] = ind2sub(size(dx), idx);

    % upper triangular terms

    % 2D outputs
    if (D == 2)

        % xi (free), xj (free)
        if (isFree(I) && isFree(J))
            objfunq{count} = sprintf(...
                '+%.16g x%d x%d + %.16g y%d y%d', ...
                2*full(V(I, J)), I, J, ...
                2*full(V(I, J)), I, J);
        % xi (free), xj (fixed)
        elseif (isFree(I) && ~isFree(J))
            objfunq{count} = sprintf(...
                '+%.16g x%d + %.16g y%d', ...
                2*full(V(I, J)) * y(J, 1), I, ...
                2*full(V(I, J)) * y(J, 2), I);
        % xi (fixed), xj (free)
        elseif (~isFree(I) && isFree(J))
            objfunq{count} = sprintf(...
                '+%.16g x%d + %.16g y%d', ...
                2*full(V(I, J)) * y(I, 1), J, ...
                2*full(V(I, J)) * y(I, 2), J);
        end
        count = count + 1;

    % 3D outputs
    elseif (D == 3)

        % xi (free), xj (free)
        if (isFree(I) && isFree(J))
            objfunq{count} = sprintf(...
                '+%.16g x%d x%d + %.16g y%d y%d + %.16g z%d z%d', ...
                2*full(V(I, J)), I, J, ...
                2*full(V(I, J)), I, J, ...
                2*full(V(I, J)), I, J);
        % xi (free), xj (fixed)
        elseif (isFree(I) && ~isFree(J))
            objfunq{count} = sprintf(...
                '+%.16g x%d + %.16g y%d + %.16g z%d', ...
                2*full(V(I, J)) * y(J, 1), I, ...
                2*full(V(I, J)) * y(J, 2), I, ...
                2*full(V(I, J)) * y(J, 3), I);
        % xi (fixed), xj (free)
        elseif (~isFree(I) && isFree(J))
            objfunq{count} = sprintf(...
                '+%.16g x%d + %.16g y%d + %.16g z%d', ...
                2*full(V(I, J)) * y(I, 1), J, ...
                2*full(V(I, J)) * y(I, 2), J, ...
                2*full(V(I, J)) * y(I, 3), J);
        end
        count = count + 1;

    else
        error('Assertion fail: Output dimension D is not 2 or 3')
    end
end

% main diagonal terms
for I = find(isFree)'

    % 2D outputs
    if (D == 2)

        % only free vertices contribute to the objective function
        objfunq{count} = sprintf(...
            '+%.16g x%d^2 + %.16g y%d^2', ...
            full(V(I, I)), I, ...
            full(V(I, I)), I);

    % 3D outputs
    elseif (D == 3)

        % only free vertices contribute to the objective function
        objfunq{count} = sprintf(...
            '+%.16g x%d^2 + %.16g y%d^2 + %.16g z%d^2', ...
            full(V(I, I)), I, ...
            full(V(I, I)), I, ...
            full(V(I, I)), I);

    else
        error('Assertion fail: Output dimension D is not 2 or 3')
    end
    count = count + 1;

end
objfunq{1} = [' obj: ' objfunq{1}];

% the linear term of the objective function (f) has to be computed at each
% iteration of the QPQC-SMACOF algorithm. Thus, it is not computed here

%% SMACOF algorithm

% file name and path to save PIP model, computed solution and initial guess
% (generate unique names so that it is possible to run several instances of
% this function in parallel)
[~, aux] = fileparts(tempname);
pipfilename = [tempdir 'model-' aux '.pip']; % PIP model
solfilename = [tempdir 'model-' aux '.sol']; % output solution
sol0filename = [tempdir 'model-' aux '.sol0']; % initial guess

% init stopCondition
stopCondition = [];

% Euclidean distances between vertices in the current solution
dy = dmatrix_con(dx, y);

% initial stress
sigma0 = 0.5 * sum(sum(w .* (dx - dy).^2));

% if dx, dy are sparse matrices, sigma0 will be a sparse scalar, and this
% gives an error with fprintf below
sigma0 = full(sigma0);

% vector of stress computed by the algorithm
sigma = zeros(1, smacof_opts.MaxIter);

% display algorithm's evolution
t = zeros(1, smacof_opts.MaxIter); % time past from 0th iteration
if (strcmp(smacof_opts.Display, 'iter'))
    fprintf('Iter\t\tSigma\t\t\tTime (sec)\n')
    fprintf('===================================================\n')
    fprintf('%d\t\t%.4e\t\t%.4e\n', 0, sigma0, 0.0)
end

% auxiliary intermediate result
mwdx = -w .* dx;

% initialize the storage of the best solution found by the algorithm. We
% start with Inf so that the initial guess will always be replaced by any
% valid solution. The reason is that we cannot be sure that the initial
% guess fulfills the constraints, but any valid solution from a later
% iteration does
sigmabest = Inf;
ybest = nan(size(y));

% majorization loop
for I = 1:smacof_opts.MaxIter

    % auxiliary matrix B: non-main-diagonal elements
    B = mwdx ./ dy;
    B(isnan(B)) = 0;

    % auxiliary matrix B: main diagonal elements
    B(1:N+1:end) = -sum(B, 2);

    % the linear term (f) of the quadratic objective function
    % 1/2 nu' * H * nu + f' * nu
    % has to be recomputed at every iteration
    f = -2 * B * y;

    % convert linear term to PIP format
    objfunl = cell(1, Nfree);
    J = find(isFree);
    for idx = 1:length(J)

        % 2D output
        if (D == 2)

            objfunl{idx} = sprintf(...
                '+%.16g x%d +%.16g y%d', ...
                f(J(idx), 1), J(idx), f(J(idx), 2), J(idx));

        % 3D output
        elseif (D == 3)

            objfunl{idx} = sprintf(...
                '+%.16g x%d +%.16g y%d +%.16g z%d', ...
                f(J(idx), 1), J(idx), f(J(idx), 2), J(idx), ...
                f(J(idx), 3), J(idx));

        else
            error('Assertion fail: Output dimension D is not 2 or 3')
        end

    end

    % create PIP file to describe problem
    fid = fopen(pipfilename, 'w');
    if (fid == -1)
        error(['Cannot open file ' pipfilename ' to save PIP model'])
    end
    fprintf(fid, '%s\n%s\n%s\n%s\n%s\n%s\n', ...
        'Minimize', objfunq{:}, objfunl{:}, bnd{:}, con{:}, 'End');
    if (fclose(fid) == -1)
        error(['Cannot close file ' pipfilename ' to save PIP model'])
    end

    % create file for the initial guess
    write_solution(sol0filename, y);

    % solve the quadratic problem
    system([...
        SCIPBIN...
        QUIETFLAG ...
        ' -c "read ' pipfilename '"'...
        strcat(scip_opts_comm{:}) ...
        ' -c "optimize"'...
        ' -c "write solution ' solfilename '"'...
        ' -c "quit"' ...
        OUTPUTREDIR]);

    % read solution
    [aux, status] = read_solution(solfilename, size(y));

    % delete temp files
    if (exist(pipfilename, 'file')), delete(pipfilename); end
    if (exist(solfilename, 'file')), delete(solfilename); end
    if (exist(sol0filename, 'file')), delete(sol0filename); end

    % if SCIP cannot find a valid solution (e.g. we set a too short time
    % limit or it doesn't exist), then we cannot update the current
    % solution. This also means that we need to stop the optimization,
    % because the problem matrices and vectors won't change, and basically
    % we would try to solve the same problem again, not reaching a solution
    % either. Note however that it's possible that we have been finding
    % valid solutions before, and it's only at a later iteration that we
    % get to this "not valid solution" situation. Thus, we cannot assume
    % that the stop condition below these lines means "error".
    if (isempty(aux))
        stopCondition{end+1} = ['SCIP: ' status];
        sigma(I) = NaN;
        t(I) = toc;
        break;
    end

    % we only have to update the positions of the free vertices. The values
    % for fixed vertices in aux are all 0, so they need to be ignored
    y(isFree, :) = aux(isFree, :);

    % check that no two vertices are overlapping. Duplicated vertices cause
    % errors in the SMACOF algorithm because they create dy=0 components,
    % that produce Inf values in B = mwdx ./ dy;. Ideally, we would like
    % to move one of them a bit so that they don't overlap, while
    % maintaining the positive orientation of the triangles they are
    % involved with. However, we don't have the triangulation in this
    % function, and leave that for future work. For the time being, we are
    % going to consider overlapping vertices as a failed solution, and exit
    if (size(unique(y, 'rows'), 1) ~= N)
        stopCondition{end+1} = 'Duplicated vertex/vertices';
        break
    end

    % recompute distances between vertices in the current solution
    dy = dmatrix_con(dx, y);

    % compute stress with the current solution
    sigma(I) = 0.5 * sum(sum(w .* (dx - dy).^2));

    % update best solution
    if (sigma(I) < sigmabest)
        sigmabest = sigma(I);
        ybest = y;
    end

    % display algorithm's evolution
    t(I) = toc;
    if (strcmp(smacof_opts.Display, 'iter'))
        fprintf('%d\t\t%.4e\t\t%.4e\n', I, sigma(I), t(I))
    end

    % check whether the stress is under the tolerance level requested by
    % the user
    if (sigma(I) < smacof_opts.TolFun)
        stopCondition{end+1} = 'TolFun';
    end

    % check whether the improvement in stress is below the user's request,
    % and it's positive (don't stop if the stress gets worse, because
    % stress can go up and down in the optimization)
    if (I > 1)
        if ((sigma(I-1)-sigma(I))/sigma(I-1) < smacof_opts.Epsilon ...
                && (sigma(I-1)-sigma(I))/sigma(I-1) >= 0)
            stopCondition{end+1} = 'Epsilon';
        end
    end

    % stop if any stop condition has been met
    if (~isempty(stopCondition))
        break;
    end

end

% return the best solution the algorithm has found in all iterations
y = ybest;

% check whether the "maximum number of iterations" stop condition has been
% met
if (I == smacof_opts.MaxIter)
    stopCondition{end+1} = 'MaxIter';
end

% prune stress and time vectors if convergence was reached before the
% maximum number of iterations
sigma(I+1:end) = [];
t(I+1:end) = [];

end

% read SCIP solution from a text file created by SCIP
%
% file: path and file name of the solution file.
%
% sz: size of the output matrix with the solution. This parameter makes the
%     code easier to write, and it also allows to detect missing variables
%     in the solution
%
% y: matrix with the solution as a point configuration (each row is a
%    point)
function [y, status] = read_solution(file, sz)

if ((sz(2) ~= 2) && (sz(2) ~= 3))
    error('We only know how to read solutions that are sets of 2D or 3D points')
end

fid = fopen(file, 'r');
if (fid == -1)
    error(['Cannot open file ' file ' to read solution'])
end

% read status of the solution
status = fgetl(fid);
if ((isnumeric(status) && status == -1) ...
        || ~strcmp(status(1:16), 'solution status:'))
    error(['Assertion fail: File with SCIP solution does not start with string ''solution status:''. File ' file])
end
status = status(18:end);

% unless we obtained a valid solution, we don't bother with the rest of the
% file (which should be empty anyway), we exit, and the main function
% should detect the empty solution y, and create a stopCondition with the
% status returned by SCIP
if (~strcmp(status, 'optimal solution found') ...
        && ~strcmp(status, 'solution limit reached'))
    fclose(fid);
    y = [];
    return;
end

% the second line in the file should be the value of the objective function
aux = fgetl(fid);
if (~strcmp(aux(1:16), 'objective value:'))
    error(['Assertion fail: Second line in file with SCIP solution does not start with string ''objective value:''. File ' file])
end

% read contents of the file. Example result:
%
% solution status: optimal solution found
% objective value:                     468.345678663118
% x1                                                 -4 	(obj:0)
% y1                                                  2 	(obj:0)
% x2                                  0.613516669331233 	(obj:0)
% y2                                                 -4 	(obj:0)
% x3                                   1.24777861008035 	(obj:0)
% y3                                  -3.43327058768579 	(obj:0)
% x4                                                 -4 	(obj:0)
% y4                                 -0.804343353251697 	(obj:0)
% x5                                                  2 	(obj:0)
% y5                                                 -4 	(obj:0)
% x6                                  -3.31069797707353 	(obj:0)
% y6                                   1.21192102496956 	(obj:0)
% x7                                     1.643412276563 	(obj:0)
% y7                                  -3.72674560989634 	(obj:0)
% quadobjvar                           468.345678663118 	(obj:1)
c = textscan(fid, '%s%f%s', 'Delimiter', ' ', 'MultipleDelimsAsOne', true);
fclose(fid);

if isempty(c{1})
    error('SCIP did not return any solution')
end

% if SCIP has found a solution, read it from the file into the output
% matrix. Note that if a variable = 0, SCIP will not put it in the output
% file
y = zeros(sz);
for I = 1:length(c{1})-1

    % variable name
    varname = c{1}{I};

    % if list of output variables has finished, we can exit the loop
    if (strcmp(varname, 'quadobjvar'))
        break;
    end

    % vertex index
    idx = str2double(varname(2:end));

    if (c{1}{I}(1) == 'x') % this is an x-coordinate
        y(idx, 1) = c{2}(I);
    elseif (c{1}{I}(1) == 'y') % this is a y-coordinate
        y(idx, 2) = c{2}(I);
    elseif (c{1}{I}(1) == 'z') % this is a z-coordinate
        y(idx, 3) = c{2}(I);
    end

end

end

% write SCIP initial solution to a text file that SCIP can understand
%
% y: matrix with the solution as a point configuration (each row is a
%    point), including the points that are not free points
%
% file: path and file name of the solution file.
function write_solution(file, y)

% size of full solution configuration
sz = size(y);

if ((sz(2) ~= 2) && (sz(2) ~= 3))
    error('We only know how to read solutions that are sets of 2D or 3D points')
end

fid = fopen(file, 'w');
if (fid == -1)
    error(['Cannot open file ' file ' to write solution'])
end

% write solution to file. Example result:
%
% x1                                                 -4 	(obj:0)
% y1                                                  2 	(obj:0)
% x2                                  0.613516669331233 	(obj:0)
% y2                                                 -4 	(obj:0)
% x3                                   1.24777861008035 	(obj:0)
% y3                                  -3.43327058768579 	(obj:0)
% x4                                                 -4 	(obj:0)
% y4                                 -0.804343353251697 	(obj:0)
% x5                                                  2 	(obj:0)
% y5                                                 -4 	(obj:0)
% x6                                  -3.31069797707353 	(obj:0)
% y6                                   1.21192102496956 	(obj:0)
% x7                                     1.643412276563 	(obj:0)
% y7                                  -3.72674560989634 	(obj:0)
for I = 1:sz(1)

    % write x-coordinate
    fprintf(fid, '%s%d\t%.16e\t%s\n', 'x', I, y(I, 1), '(obj:0)');

    % write y-coordinate
    fprintf(fid, '%s%d\t%.16e\t%s\n', 'y', I, y(I, 2), '(obj:0)');

    % write z-coordinate
    if (sz(2) == 3)
        fprintf(fid, '%s%d\t%.16e\t%s\n', 'z', I, y(I, 3), '(obj:0)');
    end

end

% close file
st = fclose(fid);
if (st == -1)
    error(['Cannot close file ' file ' after writing solution'])
end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function sF = power(sF1,sF2)
%
% Syntax
%   sF = sF1.^a
%

sF = power@S2FunHarmonic(sF1,sF2);

% try to preserve symmetry
sym = [];
if isnumeric(sF1)
  sym = sF2.CS;
elseif isnumeric(sF2)
  sym = sF1.CS;
elseif isa(sF1,'S2FunHarmonicSym') && isa(sF2,'S2FunHarmonicSym')
  sym = disjoint(sF1.CS,sF2.CS);
end

if ~isempty(sym), sF = S2FunHarmonicSym(sF.fhat,sym); end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [K] = ku0u0(x, xp, hyp, ubar, ubarp, dt, i)

logsigma = hyp(1);
logtheta = hyp(2);

a1 = hyp(3);
a2 = hyp(4);
a3 = hyp(5);

n_x = size(x,1);
n_xp = size(xp,1);

x = repmat(x,1,n_xp);
xp = repmat(xp',n_x,1);

ubar = repmat(ubar,1,n_xp);
ubarp = repmat(ubarp',n_x,1);

switch i

case 0

K=exp(1).^(logsigma+(-8).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*(x+( ...
  -1).*xp).^2).*(exp(1).^(4.*logtheta).*(exp(1).^(2.*logtheta).*(exp(1).^( ...
  2.*logtheta)+a1.*dt.*exp(1).^logtheta.*(a1.*dt.*ubar.*ubarp+(-1).*(ubar+ ...
  (-1).*ubarp).*(x+(-1).*xp))+(-1).*a1.^2.*dt.^2.*ubar.*ubarp.*(x+(-1).* ...
  xp).^2)+a2.*dt.*exp(1).^logtheta.*((-2).*exp(1).^(2.*logtheta)+exp(1) ...
  .^logtheta.*(3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*(x+(-1).*xp)).*(x+(-1).* ...
  xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^3)+a2.^2.*dt.^2.*( ...
  3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1) ...
  .*xp).^4))+a3.*dt.*exp(1).^(2.*logtheta).*((-2).*a2.*dt.*(15.*exp(1).^( ...
  3.*logtheta)+(-45).*exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1) ...
  .^logtheta.*(x+(-1).*xp).^4+(-1).*(x+(-1).*xp).^6)+exp(1).^logtheta.*( ...
  6.*exp(1).^(3.*logtheta)+(-3).*exp(1).^(2.*logtheta).*(5.*a1.*dt.*(ubar+ ...
  (-1).*ubarp)+4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).* ...
  ubarp).*(x+(-1).*xp).^5+2.*exp(1).^logtheta.*(x+(-1).*xp).^3.*(5.*a1.* ...
  dt.*(ubar+(-1).*ubarp)+x+(-1).*xp)))+a3.^2.*dt.^2.*(105.*exp(1).^(4.* ...
  logtheta)+(-420).*exp(1).^(3.*logtheta).*(x+(-1).*xp).^2+210.*exp(1).^( ...
  2.*logtheta).*(x+(-1).*xp).^4+(-28).*exp(1).^logtheta.*(x+(-1).*xp).^6+( ...
  x+(-1).*xp).^8));

case 1 % logsigma

K=exp(1).^(logsigma+(-8).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*(x+( ...
  -1).*xp).^2).*(exp(1).^(4.*logtheta).*(exp(1).^(2.*logtheta).*(exp(1).^( ...
  2.*logtheta)+a1.*dt.*exp(1).^logtheta.*(a1.*dt.*ubar.*ubarp+(-1).*(ubar+ ...
  (-1).*ubarp).*(x+(-1).*xp))+(-1).*a1.^2.*dt.^2.*ubar.*ubarp.*(x+(-1).* ...
  xp).^2)+a2.*dt.*exp(1).^logtheta.*((-2).*exp(1).^(2.*logtheta)+exp(1) ...
  .^logtheta.*(3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*(x+(-1).*xp)).*(x+(-1).* ...
  xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^3)+a2.^2.*dt.^2.*( ...
  3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1) ...
  .*xp).^4))+a3.*dt.*exp(1).^(2.*logtheta).*((-2).*a2.*dt.*(15.*exp(1).^( ...
  3.*logtheta)+(-45).*exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1) ...
  .^logtheta.*(x+(-1).*xp).^4+(-1).*(x+(-1).*xp).^6)+exp(1).^logtheta.*( ...
  6.*exp(1).^(3.*logtheta)+(-3).*exp(1).^(2.*logtheta).*(5.*a1.*dt.*(ubar+ ...
  (-1).*ubarp)+4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).* ...
  ubarp).*(x+(-1).*xp).^5+2.*exp(1).^logtheta.*(x+(-1).*xp).^3.*(5.*a1.* ...
  dt.*(ubar+(-1).*ubarp)+x+(-1).*xp)))+a3.^2.*dt.^2.*(105.*exp(1).^(4.* ...
  logtheta)+(-420).*exp(1).^(3.*logtheta).*(x+(-1).*xp).^2+210.*exp(1).^( ...
  2.*logtheta).*(x+(-1).*xp).^4+(-28).*exp(1).^logtheta.*(x+(-1).*xp).^6+( ...
  x+(-1).*xp).^8));

case 2 % logtheta

K=exp(1).^(logsigma+(-8).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*(x+( ...
  -1).*xp).^2).*(exp(1).^(5.*logtheta).*(6.*a2.^2.*dt.^2.*(exp(1) ...
  .^logtheta+(-1).*(x+(-1).*xp).^2)+exp(1).^logtheta.*(4.*exp(1).^(2.* ...
  logtheta)+3.*a1.*dt.*exp(1).^logtheta.*(a1.*dt.*ubar.*ubarp+(-1).*(ubar+ ...
  (-1).*ubarp).*(x+(-1).*xp))+(-2).*a1.^2.*dt.^2.*ubar.*ubarp.*(x+(-1).* ...
  xp).^2)+a2.*dt.*((-6).*exp(1).^(2.*logtheta)+2.*exp(1).^logtheta.*(3.* ...
  a1.*dt.*(ubar+(-1).*ubarp)+2.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.* ...
  (ubar+(-1).*ubarp).*(x+(-1).*xp).^3))+4.*exp(1).^(4.*logtheta).*(exp(1) ...
  .^(2.*logtheta).*(exp(1).^(2.*logtheta)+a1.*dt.*exp(1).^logtheta.*(a1.* ...
  dt.*ubar.*ubarp+(-1).*(ubar+(-1).*ubarp).*(x+(-1).*xp))+(-1).*a1.^2.* ...
  dt.^2.*ubar.*ubarp.*(x+(-1).*xp).^2)+a2.*dt.*exp(1).^logtheta.*((-2).* ...
  exp(1).^(2.*logtheta)+exp(1).^logtheta.*(3.*a1.*dt.*(ubar+(-1).*ubarp)+ ...
  2.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1) ...
  .*xp).^3)+a2.^2.*dt.^2.*(3.*exp(1).^(2.*logtheta)+(-6).*exp(1) ...
  .^logtheta.*(x+(-1).*xp).^2+(x+(-1).*xp).^4))+2.*a3.*dt.*exp(1).^(2.* ...
  logtheta).*((-2).*a2.*dt.*(15.*exp(1).^(3.*logtheta)+(-45).*exp(1).^(2.* ...
  logtheta).*(x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+(-1).*xp).^4+(-1).*( ...
  x+(-1).*xp).^6)+exp(1).^logtheta.*(6.*exp(1).^(3.*logtheta)+(-3).*exp(1) ...
  .^(2.*logtheta).*(5.*a1.*dt.*(ubar+(-1).*ubarp)+4.*(x+(-1).*xp)).*(x+( ...
  -1).*xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^5+2.*exp(1) ...
  .^logtheta.*(x+(-1).*xp).^3.*(5.*a1.*dt.*(ubar+(-1).*ubarp)+x+(-1).*xp)) ...
  )+(exp(1).^(4.*logtheta).*(exp(1).^(2.*logtheta).*(exp(1).^(2.*logtheta) ...
  +a1.*dt.*exp(1).^logtheta.*(a1.*dt.*ubar.*ubarp+(-1).*(ubar+(-1).*ubarp) ...
  .*(x+(-1).*xp))+(-1).*a1.^2.*dt.^2.*ubar.*ubarp.*(x+(-1).*xp).^2)+a2.* ...
  dt.*exp(1).^logtheta.*((-2).*exp(1).^(2.*logtheta)+exp(1).^logtheta.*( ...
  3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.* ...
  dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^3)+a2.^2.*dt.^2.*(3.*exp(1).^(2.* ...
  logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1).*xp).^4))+a3.* ...
  dt.*exp(1).^(2.*logtheta).*((-2).*a2.*dt.*(15.*exp(1).^(3.*logtheta)+( ...
  -45).*exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+( ...
  -1).*xp).^4+(-1).*(x+(-1).*xp).^6)+exp(1).^logtheta.*(6.*exp(1).^(3.* ...
  logtheta)+(-3).*exp(1).^(2.*logtheta).*(5.*a1.*dt.*(ubar+(-1).*ubarp)+ ...
  4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1) ...
  .*xp).^5+2.*exp(1).^logtheta.*(x+(-1).*xp).^3.*(5.*a1.*dt.*(ubar+(-1).* ...
  ubarp)+x+(-1).*xp)))+a3.^2.*dt.^2.*(105.*exp(1).^(4.*logtheta)+(-420).* ...
  exp(1).^(3.*logtheta).*(x+(-1).*xp).^2+210.*exp(1).^(2.*logtheta).*(x+( ...
  -1).*xp).^4+(-28).*exp(1).^logtheta.*(x+(-1).*xp).^6+(x+(-1).*xp).^8)).* ...
  ((-8)+(1/2).*exp(1).^((-1).*logtheta).*(x+(-1).*xp).^2)+28.*a3.^2.* ...
  dt.^2.*exp(1).^logtheta.*(15.*exp(1).^(3.*logtheta)+(-45).*exp(1).^(2.* ...
  logtheta).*(x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+(-1).*xp).^4+(-1).*( ...
  x+(-1).*xp).^6)+a3.*dt.*exp(1).^(3.*logtheta).*(6.*exp(1).^(3.*logtheta) ...
  +(-30).*a2.*dt.*(3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+( ...
  -1).*xp).^2+(x+(-1).*xp).^4)+2.*exp(1).^logtheta.*(9.*exp(1).^(2.* ...
  logtheta)+(-3).*exp(1).^logtheta.*(5.*a1.*dt.*(ubar+(-1).*ubarp)+4.*(x+( ...
  -1).*xp)).*(x+(-1).*xp)+(x+(-1).*xp).^3.*(5.*a1.*dt.*(ubar+(-1).*ubarp)+ ...
  x+(-1).*xp))+(-3).*exp(1).^(2.*logtheta).*(5.*a1.*dt.*(ubar+(-1).*ubarp) ...
  +4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+( ...
  -1).*xp).^5+2.*exp(1).^logtheta.*(x+(-1).*xp).^3.*(5.*a1.*dt.*(ubar+(-1) ...
  .*ubarp)+x+(-1).*xp)));

case 3 % a1

K=dt.*exp(1).^(logsigma+(-5).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*( ...
  x+(-1).*xp).^2).*(exp(1).^(2.*logtheta).*(exp(1).^logtheta.*(2.*a1.*dt.* ...
  ubar.*ubarp.*(exp(1).^logtheta+(-1).*(x+(-1).*xp).^2)+(-1).*exp(1) ...
  .^logtheta.*(ubar+(-1).*ubarp).*(x+(-1).*xp))+(-1).*a2.*dt.*(ubar+(-1).* ...
  ubarp).*((-3).*exp(1).^logtheta+(x+(-1).*xp).^2).*(x+(-1).*xp))+a3.*dt.* ...
  (ubar+(-1).*ubarp).*((-15).*exp(1).^(2.*logtheta)+10.*exp(1).^logtheta.* ...
  (x+(-1).*xp).^2+(-1).*(x+(-1).*xp).^4).*(x+(-1).*xp));

case 4 % a2

K=dt.*exp(1).^(logsigma+(-6).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*( ...
  x+(-1).*xp).^2).*(exp(1).^(2.*logtheta).*(exp(1).^logtheta.*((-2).*exp( ...
  1).^(2.*logtheta)+exp(1).^logtheta.*(3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*( ...
  x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).* ...
  xp).^3)+2.*a2.*dt.*(3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+ ...
  (-1).*xp).^2+(x+(-1).*xp).^4))+(-2).*a3.*dt.*(15.*exp(1).^(3.*logtheta)+ ...
  (-45).*exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+( ...
  -1).*xp).^4+(-1).*(x+(-1).*xp).^6));

case 5 % a3

K=dt.*exp(1).^(logsigma+(-8).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*( ...
  x+(-1).*xp).^2).*(exp(1).^(2.*logtheta).*((-2).*a2.*dt.*(15.*exp(1).^( ...
  3.*logtheta)+(-45).*exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1) ...
  .^logtheta.*(x+(-1).*xp).^4+(-1).*(x+(-1).*xp).^6)+exp(1).^logtheta.*( ...
  6.*exp(1).^(3.*logtheta)+(-3).*exp(1).^(2.*logtheta).*(5.*a1.*dt.*(ubar+ ...
  (-1).*ubarp)+4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).* ...
  ubarp).*(x+(-1).*xp).^5+2.*exp(1).^logtheta.*(x+(-1).*xp).^3.*(5.*a1.* ...
  dt.*(ubar+(-1).*ubarp)+x+(-1).*xp)))+2.*a3.*dt.*(105.*exp(1).^(4.* ...
  logtheta)+(-420).*exp(1).^(3.*logtheta).*(x+(-1).*xp).^2+210.*exp(1).^( ...
  2.*logtheta).*(x+(-1).*xp).^4+(-28).*exp(1).^logtheta.*(x+(-1).*xp).^6+( ...
  x+(-1).*xp).^8));

otherwise

        K = zeros(n_x, n_xp);
end

if K == 0

    K = zeros(n_x, n_xp);

end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function c = ldivide(a,b)
%LDIVIDE      Interval elementwise left division a .\ b
%

% written  11/02/05     S.M. Rump
%

  c = b ./ a;

repo: OpenMAS


text:
function [ fa ] = gsp_cfa( T,fs )
%GSP_CFA Create frequency axis
%   Usage: fa = gsp_cfa(N);
%          fa = gsp_cfa(N,fs);
%
%   Input parameters:
%       N   : Number of samples
%       fs  : Sampling frequency (default 1)
%   Ouput parameters:
%       fa  : Frequency axis
%
%   This function create a normalized frequency axis corresponding to the
%   frequency of the function fft.
%

% Author: Nathanael Perraudin

if nargin<2
    fs = 1;
end

if mod(T,2)
    fa = linspace(-0.5+1/(2*T),0.5-1/(2*T),T);
else
    fa = linspace(-0.5,0.5-1/T,T);
end

fa = fa*fs;

fa = ifftshift(fa);
fa = fa(:);

end

repo: gspbox


text:
classdef factorize
%FACTORIZE an object-oriented method for solving linear systems
% and least-squares problems, and for representing operations with the
% inverse of a square matrix or the pseudo-inverse of a rectangular matrix.
%
% F = factorize(A) returns an object F that holds the factorization of a
% non-singular matrix A.  x=F\b then solves a linear system or a
% least-squares problem.  S=inverse(F) or S=inverse(A) returns a factorized
% representation of the inverse of A so that inverse(A)*b is mathematically
% equivalent to inv(A)*b, but the former does not actually compute the
% inverse of A.
%
% Example
%
%   F = factorize(A) ;      % LU, QR, or Cholesky factorization of A
%   x = F\b ;               % solve A*x=b; same as x=A\b
%   S = inverse (F) ;       % S represents the factorization of inv(A)
%   x = S*b ;               % same as x = A\b.
%   S = A-B*inverse(D)*C    % efficiently computes the Schur complement
%   S = A-B*inv(D)*C        % bad method for computing the Schur complement
%   S = inverse(A) ; S(:,1) % compute just the first column of inv(A),
%                           % without computing inv(A)
%
% If A is square, symmetric (Hermitian for the complex case), and has a
% real positive diagonal, then use F=factorize(A,1).  If you know
% otherwise, use F=factorize(A,0).  Using this option improves performance,
% since otherwise this condition must be checked to choose between a
% Cholesky or LU factorization.  The option is ignored if A is rectangular.
%
% For more details, type "help factorize1".
% For a demo type "fdemo" or see the html/ directory.
%
% See also inverse, factorize1, mldivide, mrdivide, inv, pinv, linsolve

% Copyright 2009, Timothy A. Davis, University of Florida

    properties (SetAccess = protected)
        % The factorize object holds a QR, LU, Cholesky factorization:
        A = [ ] ;           % a copy of the input matrix
        L = [ ] ;           % lower-triangular factor for LU and Cholesky
        U = [ ] ;           % upper-triangular factor for LU
        Q = [ ] ;           % Q factor for dense QR
        R = [ ] ;           % R factor for QR
        p = [ ] ;           % sparse row permutation matrix
        q = [ ] ;           % sparse column permutation matrix
        is_inverse = false ; % F represents the factorization of A or inv(A)
        kind = 0 ;          % F is one of 8 kinds of factorizations
    end

    methods

        function F = factorize (A,try_chol)

            % factorize constructor: compute a factorization of A

            if (ndims (A) > 2)
                error ('Matrix must be 2D.') ;
            end
            [m n] = size (A) ;
            F.A = A ;

            if (m > n)

                % QR factorization of A
                if (issparse (A))
                    % Q-less econonmy sparse QR: (A*q)'*(A*q) = R'*R
                    q = sparse (colamd (A), 1:n, 1) ;
                    R = qr (A*q, 0) ;
                    F.q = q ;
                    F.kind = 1 ;
                else
                    % dense economy QR factorization: A = Q*R
                    [Q R] = qr (A,0) ;
                    F.Q = Q ;
                    F.kind = 2 ;
                end
                ok = (nnz (diag (R)) == n) ;
                F.R = R ;

            elseif (m < n)

                % QR factorization of A'
                if (issparse (A))
                    % Q-less economy sparse QR: (p*A)*(p*A)' = R'*R
                    C = A' ;
                    p = sparse (1:m, colamd (C), 1) ;
                    R = qr (C*p', 0) ;
                    F.p = p ;
                    F.kind = 3 ;
                else
                    % dense economy LQ factorization: A' = Q*R
                    [Q R] = qr (A',0) ;
                    F.Q = Q ;
                    F.kind = 4 ;
                end
                ok = (nnz (diag (R)) == m) ;
                F.R = R ;

            else

                % Cholesky or LU factorization of A
                g = 1 ;
                if (nargin == 1)
                    % This is an expensive test, so skip it if the caller
                    % already knows the matrix is a candidate for Cholesky.
                    d = diag (A) ;
                    try_chol = (all (d > 0) && nnz (imag (d)) == 0 && ...
                        nnz (A-A') == 0) ;
                end
                if (try_chol)
                    if (issparse (A))
                        % sparse Cholesky factorization: q'*A*q = L*L'
                        [L g q] = chol (A, 'lower') ;
                    else
                        % dense Cholesky factorization: A = R'*R
                        [R g] = chol (A) ;
                    end
                end
                % do an LU factorization if Cholesky failed or was skipped
                ok = (g == 0) ;
                if (ok)
                    % Cholesky was successful
                    if (issparse (A))
                        F.L = L ;
                        F.q = q ;
                        F.kind = 5 ;
                    else
                        F.R = R ;
                        F.kind = 6 ;
                    end
                else
                    % need an LU factorization
                    if (issparse (A))
                        % sparse LU factorization: p*A*q = L*U
                        [L U p q] = lu (A) ;
                        F.q = q ;
                        F.kind = 7 ;
                    else
                        % dense LU factorization: p*A = L*U
                        [L U p] = lu (A, 'vector') ;
                        p = sparse (1:n, p, 1) ;
                        F.kind = 8 ;
                    end
                    F.L = L ;
                    F.U = U ;
                    F.p = p ;
                    ok = (nnz (diag (U)) == n)  ;
                end
            end

            if (~ok)
                error ('Matrix is rank deficient.') ;
            end
        end
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function obj = eq_hand_eye_func_new(in1,in2,in3,in4,in5)
coef_f0_q_sym1 = in1(:,1);
coef_f0_q_sym2 = in1(:,2);
coef_f0_q_sym3 = in1(:,3);
coef_f0_q_sym4 = in1(:,4);
coef_f0_q_sym5 = in1(:,5);
coef_f0_q_sym6 = in1(:,6);
coef_f0_q_sym7 = in1(:,7);
coef_f0_q_sym8 = in1(:,8);
coef_f0_q_sym9 = in1(:,9);
coef_f0_q_sym10 = in1(:,10);
coef_f0_q_sym11 = in1(:,11);
coef_f1_q_sym1 = in2(:,1);
coef_f0_q_sym12 = in1(:,12);
coef_f1_q_sym2 = in2(:,2);
coef_f0_q_sym13 = in1(:,13);
coef_f1_q_sym3 = in2(:,3);
coef_f0_q_sym14 = in1(:,14);
coef_f1_q_sym4 = in2(:,4);
coef_f0_q_sym15 = in1(:,15);
coef_f1_q_sym5 = in2(:,5);
coef_f0_q_sym16 = in1(:,16);
coef_f1_q_sym6 = in2(:,6);
coef_f0_q_sym17 = in1(:,17);
coef_f1_q_sym7 = in2(:,7);
coef_f0_q_sym18 = in1(:,18);
coef_f1_q_sym8 = in2(:,8);
coef_f0_q_sym19 = in1(:,19);
coef_f1_q_sym9 = in2(:,9);
coef_f0_q_sym20 = in1(:,20);
coef_f0_q_sym21 = in1(:,21);
coef_f2_q_sym1 = in3(:,1);
coef_f0_q_sym22 = in1(:,22);
coef_f2_q_sym2 = in3(:,2);
coef_f0_q_sym23 = in1(:,23);
coef_f2_q_sym3 = in3(:,3);
coef_f0_q_sym24 = in1(:,24);
coef_f2_q_sym4 = in3(:,4);
coef_f2_q_sym5 = in3(:,5);
coef_f2_q_sym6 = in3(:,6);
coef_f2_q_sym7 = in3(:,7);
coef_f2_q_sym8 = in3(:,8);
coef_f2_q_sym9 = in3(:,9);
coef_f3_q_sym1 = in4(:,1);
coef_f3_q_sym2 = in4(:,2);
coef_f3_q_sym3 = in4(:,3);
coef_f3_q_sym4 = in4(:,4);
coef_f3_q_sym5 = in4(:,5);
coef_f3_q_sym6 = in4(:,6);
coef_f3_q_sym7 = in4(:,7);
coef_f3_q_sym8 = in4(:,8);
coef_f3_q_sym9 = in4(:,9);
coef_f1_q_sym10 = in2(:,10);
coef_f1_q_sym11 = in2(:,11);
coef_f1_q_sym12 = in2(:,12);
coef_f1_q_sym13 = in2(:,13);
coef_f1_q_sym14 = in2(:,14);
coef_f1_q_sym15 = in2(:,15);
coef_f1_q_sym16 = in2(:,16);
coef_f1_q_sym17 = in2(:,17);
coef_f1_q_sym18 = in2(:,18);
coef_f1_q_sym19 = in2(:,19);
coef_f1_q_sym20 = in2(:,20);
coef_f1_q_sym21 = in2(:,21);
coef_f1_q_sym22 = in2(:,22);
coef_f1_q_sym23 = in2(:,23);
coef_f1_q_sym24 = in2(:,24);
coef_f2_q_sym10 = in3(:,10);
coef_f2_q_sym11 = in3(:,11);
coef_f2_q_sym12 = in3(:,12);
coef_f2_q_sym13 = in3(:,13);
coef_f2_q_sym14 = in3(:,14);
coef_f2_q_sym15 = in3(:,15);
coef_f2_q_sym16 = in3(:,16);
coef_f2_q_sym17 = in3(:,17);
coef_f2_q_sym18 = in3(:,18);
coef_f2_q_sym19 = in3(:,19);
coef_f2_q_sym20 = in3(:,20);
coef_f2_q_sym21 = in3(:,21);
coef_f2_q_sym22 = in3(:,22);
coef_f2_q_sym23 = in3(:,23);
coef_f2_q_sym24 = in3(:,24);
coef_f3_q_sym10 = in4(:,10);
coef_f3_q_sym11 = in4(:,11);
coef_f3_q_sym12 = in4(:,12);
coef_f3_q_sym13 = in4(:,13);
coef_f3_q_sym14 = in4(:,14);
coef_f3_q_sym15 = in4(:,15);
coef_f3_q_sym16 = in4(:,16);
coef_f3_q_sym17 = in4(:,17);
coef_f3_q_sym18 = in4(:,18);
coef_f3_q_sym19 = in4(:,19);
coef_f3_q_sym20 = in4(:,20);
coef_f3_q_sym21 = in4(:,21);
coef_f3_q_sym22 = in4(:,22);
coef_f3_q_sym23 = in4(:,23);
coef_f3_q_sym24 = in4(:,24);
q0 = in5(1,:);
q1 = in5(2,:);
q2 = in5(3,:);
q3 = in5(4,:);
t2 = q0.^2;
t3 = q0.^3;
t4 = q1.^2;
t6 = q1.^3;
t7 = q2.^2;
t8 = q2.^3;
t9 = q3.^2;
t10 = q3.^3;
t5 = t2.^2;
obj = [-coef_f1_q_sym1.*t5+coef_f0_q_sym18.*t4-coef_f1_q_sym11.*t2+coef_f0_q_sym12.*t4.^2+coef_f0_q_sym11.*q0.*q1+coef_f0_q_sym22.*q1.*q2+coef_f0_q_sym24.*q1.*q3-coef_f1_q_sym18.*q0.*q1-coef_f1_q_sym22.*q0.*q2-coef_f1_q_sym24.*q0.*q3+coef_f0_q_sym1.*q1.*t3+coef_f0_q_sym5.*q0.*t6-coef_f1_q_sym2.*q1.*t3-coef_f1_q_sym3.*q2.*t3-coef_f1_q_sym4.*q3.*t3+coef_f0_q_sym13.*q2.*t6+coef_f0_q_sym14.*q3.*t6+coef_f0_q_sym19.*q1.*t8+coef_f0_q_sym23.*q1.*t10-coef_f1_q_sym12.*q0.*t6-coef_f1_q_sym19.*q0.*t8-coef_f1_q_sym23.*q0.*t10+coef_f0_q_sym2.*t2.*t4-coef_f1_q_sym5.*t2.*t4+coef_f0_q_sym15.*t4.*t7-coef_f1_q_sym8.*t2.*t7+coef_f0_q_sym17.*t4.*t9-coef_f1_q_sym10.*t2.*t9+coef_f0_q_sym3.*q1.*q2.*t2+coef_f0_q_sym4.*q1.*q3.*t2+coef_f0_q_sym6.*q0.*q2.*t4+coef_f0_q_sym7.*q0.*q3.*t4+coef_f0_q_sym8.*q0.*q1.*t7+coef_f0_q_sym10.*q0.*q1.*t9-coef_f1_q_sym6.*q1.*q2.*t2-coef_f1_q_sym7.*q1.*q3.*t2+coef_f0_q_sym16.*q2.*q3.*t4-coef_f1_q_sym9.*q2.*q3.*t2+coef_f0_q_sym20.*q1.*q3.*t7+coef_f0_q_sym21.*q1.*q2.*t9-coef_f1_q_sym13.*q0.*q2.*t4-coef_f1_q_sym14.*q0.*q3.*t4-coef_f1_q_sym15.*q0.*q1.*t7-coef_f1_q_sym17.*q0.*q1.*t9-coef_f1_q_sym20.*q0.*q3.*t7-coef_f1_q_sym21.*q0.*q2.*t9+coef_f0_q_sym9.*q0.*q1.*q2.*q3-coef_f1_q_sym16.*q0.*q1.*q2.*q3;-coef_f2_q_sym1.*t5+coef_f0_q_sym22.*t7-coef_f2_q_sym11.*t2+coef_f0_q_sym19.*t7.^2+coef_f0_q_sym11.*q0.*q2+coef_f0_q_sym18.*q1.*q2+coef_f0_q_sym24.*q2.*q3-coef_f2_q_sym18.*q0.*q1-coef_f2_q_sym22.*q0.*q2-coef_f2_q_sym24.*q0.*q3+coef_f0_q_sym1.*q2.*t3+coef_f0_q_sym8.*q0.*t8+coef_f0_q_sym12.*q2.*t6+coef_f0_q_sym15.*q1.*t8-coef_f2_q_sym2.*q1.*t3-coef_f2_q_sym3.*q2.*t3-coef_f2_q_sym4.*q3.*t3+coef_f0_q_sym20.*q3.*t8+coef_f0_q_sym23.*q2.*t10-coef_f2_q_sym12.*q0.*t6-coef_f2_q_sym19.*q0.*t8-coef_f2_q_sym23.*q0.*t10+coef_f0_q_sym3.*t2.*t7+coef_f0_q_sym13.*t4.*t7-coef_f2_q_sym5.*t2.*t4+coef_f0_q_sym21.*t7.*t9-coef_f2_q_sym8.*t2.*t7-coef_f2_q_sym10.*t2.*t9+coef_f0_q_sym2.*q1.*q2.*t2+coef_f0_q_sym4.*q2.*q3.*t2+coef_f0_q_sym5.*q0.*q2.*t4+coef_f0_q_sym6.*q0.*q1.*t7+coef_f0_q_sym9.*q0.*q3.*t7+coef_f0_q_sym10.*q0.*q2.*t9+coef_f0_q_sym14.*q2.*q3.*t4+coef_f0_q_sym16.*q1.*q3.*t7+coef_f0_q_sym17.*q1.*q2.*t9-coef_f2_q_sym6.*q1.*q2.*t2-coef_f2_q_sym7.*q1.*q3.*t2-coef_f2_q_sym9.*q2.*q3.*t2-coef_f2_q_sym13.*q0.*q2.*t4-coef_f2_q_sym14.*q0.*q3.*t4-coef_f2_q_sym15.*q0.*q1.*t7-coef_f2_q_sym17.*q0.*q1.*t9-coef_f2_q_sym20.*q0.*q3.*t7-coef_f2_q_sym21.*q0.*q2.*t9+coef_f0_q_sym7.*q0.*q1.*q2.*q3-coef_f2_q_sym16.*q0.*q1.*q2.*q3;coef_f0_q_sym24.*t9-coef_f3_q_sym1.*t5-coef_f3_q_sym11.*t2+coef_f0_q_sym23.*t9.^2+coef_f0_q_sym11.*q0.*q3+coef_f0_q_sym18.*q1.*q3+coef_f0_q_sym22.*q2.*q3-coef_f3_q_sym18.*q0.*q1-coef_f3_q_sym22.*q0.*q2-coef_f3_q_sym24.*q0.*q3+coef_f0_q_sym1.*q3.*t3+coef_f0_q_sym10.*q0.*t10+coef_f0_q_sym12.*q3.*t6+coef_f0_q_sym17.*q1.*t10+coef_f0_q_sym19.*q3.*t8+coef_f0_q_sym21.*q2.*t10-coef_f3_q_sym2.*q1.*t3-coef_f3_q_sym3.*q2.*t3-coef_f3_q_sym4.*q3.*t3-coef_f3_q_sym12.*q0.*t6-coef_f3_q_sym19.*q0.*t8-coef_f3_q_sym23.*q0.*t10+coef_f0_q_sym4.*t2.*t9+coef_f0_q_sym14.*t4.*t9+coef_f0_q_sym20.*t7.*t9-coef_f3_q_sym5.*t2.*t4-coef_f3_q_sym8.*t2.*t7-coef_f3_q_sym10.*t2.*t9+coef_f0_q_sym2.*q1.*q3.*t2+coef_f0_q_sym3.*q2.*q3.*t2+coef_f0_q_sym5.*q0.*q3.*t4+coef_f0_q_sym7.*q0.*q1.*t9+coef_f0_q_sym8.*q0.*q3.*t7+coef_f0_q_sym9.*q0.*q2.*t9+coef_f0_q_sym13.*q2.*q3.*t4+coef_f0_q_sym15.*q1.*q3.*t7+coef_f0_q_sym16.*q1.*q2.*t9-coef_f3_q_sym6.*q1.*q2.*t2-coef_f3_q_sym7.*q1.*q3.*t2-coef_f3_q_sym9.*q2.*q3.*t2-coef_f3_q_sym13.*q0.*q2.*t4-coef_f3_q_sym14.*q0.*q3.*t4-coef_f3_q_sym15.*q0.*q1.*t7-coef_f3_q_sym17.*q0.*q1.*t9-coef_f3_q_sym20.*q0.*q3.*t7-coef_f3_q_sym21.*q0.*q2.*t9+coef_f0_q_sym6.*q0.*q1.*q2.*q3-coef_f3_q_sym16.*q0.*q1.*q2.*q3;t2+t4+t7+t9-1.0];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%% FUNCTION Logistic_CMTL
%   Convex-relaxed Clustered Multi-Task Learning with Logistic Loss.
%
%% OBJECTIVE
%   argmin_{W,M,C} { sum_i^t (- sum(log (1./ (1+ exp(-X{i}*W(:, i) - Y{i} .* C(i)))))/length(Y{i}))
%            + rho1 * eta (1+eta) trace(W (eta I + M)^-1 W')
%     subject to: trace (M) = k, M \preceq I, M \in S_+^t, eta = rho2/rho1
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   k: cluster number
%   rho1: clustering penalty controlling parameter (rho1=0 then tasks are not
%   related)
%   rho2: L2 norm regularization on model W (rho2=0 then reduce to ASO but
%   this solver does not support this case)
%
%% OUTPUT
%   W: model: d * t
%   C: model: 1 * t
%   funcVal: function value vector.
%   M: relaxed F*F', where F is the clustering assignment matrix
%
%% LICENSE
%   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/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou, Jianhui Chen and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% RELATED PAPERS
%
%   [1] J. Zhou, J. Chen and J. Ye, Clustered Multi-Task Learning via
%       Alternating Structure Optimization, NIPS 2011.
%
%% RELATED FUNCTIONS
%   Least_CMTL, init_opts

%% Code starts here
function [W, C, funcVal, M] = Logistic_CMTL(X, Y, rho1, rho2, k, opts)

if nargin <5
    error('\n Inputs: X, Y, rho1, rho2 and k should be specified!\n');
end
X = multi_transpose(X);

if nargin <6
    opts = [];
end

if rho2<=0 || rho1<=0
    error('rho1 and rho2 should both greater than zero.');
end

% if exist('mosekopt','file')==0
%     error('Mosek is not found. Please install Mosek first. \n')
% end

% initialize options.
opts=init_opts(opts);

task_num  = length (X);
dimension = size(X{1}, 1);
funcVal = [];

eta = rho2 / rho1;
c = rho1 * eta * (1 + eta);

%initialize a starting point
C0_prep = zeros(1, task_num);
for t_idx = 1: task_num
    m1 = nnz(Y{t_idx} == 1);
    m2 = nnz(Y{t_idx} == -1);
    if ( m1==0 || m2==0 )
        C0_prep(t_idx) = 0;
    else
        C0_prep(t_idx) = log(m1/m2);
    end
end

if opts.init==2
    W0 = zeros(dimension, task_num);
    C0 = zeros(1, task_num);
elseif opts.init== 0
    W0 = randn(dimension, task_num);
    C0 = C0_prep;
else
    if isfield(opts,'W0')
        W0=opts.W0;
        if (nnz(size(W0)-[dimension, task_num]))
            error('\n Check the input .W0');
        end
    else
        W0 = zeros(dimension, task_num);
    end
    if isfield(opts,'C0')
        C0=opts.C0;
    else
        C0=C0_prep;
    end
end

M0 = speye (task_num) * k / task_num;

bFlag=0; % this flag tests whether the gradient step only changes a little

Wz= W0;
Cz= C0;
Mz = M0;

Wz_old = W0;
Cz_old = C0;
Mz_old = M0;

t = 1;
t_old = 0;
iter = 0;
gamma = 1;
gamma_inc = 2;

while iter < opts.maxIter
    alpha = (t_old - 1) /t;

    Ws = (1 + alpha) * Wz - alpha * Wz_old;
    Cs = (1 + alpha) * Cz - alpha * Cz_old;
    Ms = (1 + alpha) * Mz - alpha * Mz_old;

    % compute function value and gradients of the search point
    [gWs, gCs, gMs, Fs ]  = gradVal_eval(Ws, Cs, Ms);

    % the Armijo Goldstein line search scheme
    while true
        %[Wzp l1c_wzp] = l1_projection(Ws - gWs/gamma, 2 * rho2 / gamma);
        %Fzp = funVal_eval  (Wzp, Czp, rho1);

        Wzp = Ws - gWs/gamma;
        Czp = Cs - gCs/gamma;
        [Mzp Mzp_Pz Mzp_DiagSigz ] = singular_projection (Ms - gMs/gamma, k);
        Fzp = funVal_eval (Wzp, Czp, Mzp_Pz, Mzp_DiagSigz);

        delta_Wzp = Wzp - Ws;
        delta_Czp = Czp - Cs;
        delta_Mzp = Mzp - Ms;

        nrm_delta_Wzp = norm(delta_Wzp, 'fro')^2;
        nrm_delta_Czp = norm(delta_Czp, 'fro')^2;
        nrm_delta_Mzp = norm(delta_Mzp, 'fro')^2;

        r_sum = (nrm_delta_Wzp+nrm_delta_Czp+nrm_delta_Mzp)/3;

        Fzp_gamma = Fs + sum(sum(delta_Wzp .* gWs)) ...
            + sum(sum(delta_Czp .* gCs))...
            + sum(sum(delta_Mzp .* gMs)) ...
            + gamma/2 * nrm_delta_Wzp ...
            + gamma/2 * nrm_delta_Mzp ...
            + gamma/2 * nrm_delta_Czp;

        if (r_sum <=eps)
            bFlag=1; % this shows that, the gradient step makes little improvement
            break;
        end

        if (Fzp <= Fzp_gamma)
            break;
        else
            gamma = gamma * gamma_inc;
        end
    end

    Wz_old = Wz;
    Cz_old = Cz;
    Mz_old = Mz;

    Wz = Wzp;
    Cz = Czp;
    Mz = Mzp;

    funcVal = cat(1, funcVal, Fzp);

    if (bFlag)
        % fprintf('\n The program terminates as the gradient step changes the solution very small.');
        break;
    end

    % test stop condition.
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

    iter = iter + 1;
    t_old = t;
    t = 0.5 * (1 + (1+ 4 * t^2)^0.5);

end

W = Wzp;
C = Czp;
M = Mzp;

% private functions

    function [Mzp Mzp_Pz Mzp_DiagSigz ] = singular_projection (Msp, k)
        [EVector EValue] = eig(Msp);
        Pz = real(EVector);  diag_EValue = real(diag(EValue));
        %DiagSigz = SingVal_Projection(diag_EValue, k);
        DiagSigz = bsa_ihb(diag_EValue, ones(size(diag_EValue)), k, ones(size(diag_EValue)));
        Mzp = Pz * diag(DiagSigz) *Pz';
        Mzp_Pz = Pz;
        Mzp_DiagSigz = DiagSigz;
    end

    function [grad_W, grad_C, grad_M, funcVal] = gradVal_eval(W, C, M)
        IM = (eta * speye(task_num) + M);
        invEtaMWt = IM\W';

        grad_W = zeros(dimension, task_num);
        grad_C = zeros(1, task_num);
        lossValVect = zeros (1 , task_num);
        if opts.pFlag
            parfor i = 1:task_num
                [ grad_W(:, i), grad_C(:, i), lossValVect(:, i)] = unit_grad_eval( W(:, i), C(i), X{i}, Y{i});
            end
        else
            for i = 1:task_num
                [ grad_W(:, i), grad_C(:, i), lossValVect(:, i)] = unit_grad_eval( W(:, i), C(i), X{i}, Y{i});
            end
        end
        grad_W = grad_W + 2 * c * invEtaMWt';   %W component
        grad_M = - c * (W' * W / IM /IM );      %M component

        funcVal = sum(lossValVect) + c * trace( W * invEtaMWt);
    end

    function [funcVal] = funVal_eval (W, C, M_Pz, M_DiagSigz)
        invIM = M_Pz * (diag( 1./(eta + M_DiagSigz))) * M_Pz';
        invEtaMWt = invIM * W';

        funcVal = 0;
        if opts.pFlag
            parfor i = 1: task_num
                funcVal = funcVal + unit_funcVal_eval( W(:, i), C(i), X{i}, Y{i});
            end
        else
            for i = 1: task_num
                funcVal = funcVal + unit_funcVal_eval( W(:, i), C(i), X{i}, Y{i});
            end
        end
        funcVal = funcVal + c * trace( W * invEtaMWt);
    end

end

function [ grad_w, grad_c, funcVal ] = unit_grad_eval( w, c, x, y)
%gradient and logistic evaluation for each task
m = length(y);
weight = ones(m, 1)/m;
weighty = weight.* y;
aa = -y.*(x'*w + c);
bb = max( aa, 0);
funcVal = weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb );
pp = 1./ (1+exp(aa));
b = -weighty.*(1-pp);
grad_c = sum(b);
grad_w = x * b;
end
function [ funcVal ] = unit_funcVal_eval( w, c, x, y)
%function value evaluation for each task
m = length(y);
weight = ones(m, 1)/m;
aa = -y.*(x'*w + c);
bb = max( aa, 0);
funcVal = weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb );
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function [ops, stat, Fcell, FcellNeu]      = extractSignalsNoOverlaps(ops, m, stat)

ops.saveNeuropil = getOr(ops, 'saveNeuropil', 0);

Nk = numel(stat);

Ny = numel(ops.yrange);
Nx = numel(ops.xrange);

% make new set of basis functions (larger coverage)
ops.neuropilRange = 10;

S = getNeuropilBasis(ops, Ny, Nx, 'raisedcosyne'); % 'raisedcosyne', 'Fourier'
S = normc(S);

% S = m.S;

S = reshape(S, [], size(S, ndims(S)));
nBasis = size(S,2);

% initialize mask
maskNeu = ones(size(S,1), 1);

stat = getNonOverlapROIs(stat, Ny, Nx);

LtS = zeros(Nk, size(S,2));
for k = 1:Nk
    ix = stat(k).ipix(~stat(k).isoverlap);
    maskNeu(stat(k).ipix)= 0;
    if numel(ix)==0 || sum(~stat(k).isoverlap)==0
        LtS(k,:) = 0;
    else
        LtS(k,:) = stat(k).lam(~stat(k).isoverlap)' * S(ix, :);
    end
end

% add all pixels within X um
if isfield(ops, 'exclFracCell') && ops.exclFracCell>0
    H       = fspecial('disk', round(ops.diameter * ops.exclFracCell));
    maskNeu = reshape(maskNeu, Ny, Nx);
    maskNeu = imfilter(maskNeu, H, 'replicate');
    maskNeu = single(maskNeu(:) > 1-1e-3);
end
%% get signals
S    = bsxfun(@times, S, maskNeu(:));
StS = S' * S;
StS = StS + 1e-2 * eye(size(StS));

nimgbatch = 2000;

ix = 0;
fclose all;
fid = fopen(ops.RegFile, 'r');

tic
F        = zeros(Nk, sum(ops.Nframes), 'single');
Fneu    = zeros(Nk, sum(ops.Nframes), 'single');
if ops.saveNeuropil
    Ntraces = zeros(nBasis, sum(ops.Nframes), 'single');
end
% S    = bsxfun(@times, S, maskNeu(:));

[Ly Lx] = size(ops.mimg1);

mimg1 = ops.mimg1(ops.yrange, ops.xrange);
% find indices of good clusters
while 1
    data = fread(fid,  Ly*Lx*nimgbatch, '*int16');
    if isempty(data)
       break;
    end
    data = reshape(data, Ly, Lx, []);
    data = data(ops.yrange, ops.xrange, :);
    data = single(data);
    NT   = size(data,3);

    % process the data
    data = bsxfun(@minus, data, mimg1);
    data = my_conv2(data, ops.sig, [1 2]);
    data = bsxfun(@rdivide, data, m.sdmov);
    data = single(reshape(data, [], NT));

    %
    Ftemp = zeros(Nk, NT, 'single');
    for k = 1:Nk
       ipix = stat(k).ipix(~stat(k).isoverlap)';
       if ~isempty(ipix)
           Ftemp(k,:) = stat(k).lam(~stat(k).isoverlap)' * data(ipix,:);
       end
    end
    F(:,ix + (1:NT))    = Ftemp;

    Tneu                = StS\(S' * data);
    Ftemp2              = LtS * Tneu;
    Fneu(:,ix + (1:NT)) = Ftemp2;

%     Fneu(:,ix + (1:NT))     = m.LtS * Fdeconv(1+Nk:end, :); % estimated neuropil
%     F(:,ix + (1:NT))        = Fneu(:,ix + (1:NT)) + Fdeconv(1:Nk, :); % estimated ROI signal

    if ops.saveNeuropil
        Ntraces(:,ix + (1:NT)) = Tneu;
    end

    ix = ix + NT;
    if rem(ix, 3*NT)==0
        fprintf('Frame %d done in time %2.2f \n', ix, toc)
    end
end
fclose(fid);
%% add the means back in to both neuropil and total
data = my_conv2(mimg1, ops.sig, [1 2]);
data = bsxfun(@rdivide, data, m.sdmov);
data = single(reshape(data, [], 1));

scalefactors = nan(numel(stat),1);
Ftemp = zeros(Nk, 1, 'single');
for k = 1:Nk
    ipix = stat(k).ipix(~stat(k).isoverlap)';
    if ~isempty(ipix)
        Ftemp(k,:) = stat(k).lam(~stat(k).isoverlap)' * data(ipix,1);
        scalefactors(k) = mean(m.sdmov(ipix));
    end
end

Tneu                = StS\(S' * data);
Ftemp2              = LtS * Tneu;

Fneu     = bsxfun(@plus, Fneu, Ftemp2); % estimated neuropil
F        = bsxfun(@plus, F,    Ftemp);

Fneu     = bsxfun(@times, Fneu, scalefactors); % estimated neuropil
F        = bsxfun(@times, F,    scalefactors);

%%
% get activity stats
[stat, F, Fneu] = getActivityStats(ops, stat, F, Fneu);

%
csumNframes = [0 cumsum(ops.Nframes)];
Fcell       = cell(1, length(ops.Nframes));
FcellNeu    = cell(1, length(ops.Nframes));
for i = 1:length(ops.Nframes)
    Fcell{i}     = F(:, csumNframes(i) + (1:ops.Nframes(i)));
    FcellNeu{i}  = Fneu(:, csumNframes(i) + (1:ops.Nframes(i)));
end

if getOr(ops, 'saveNeuropil', 0)
    S = reshape(S, numel(ops.yrange), numel(ops.xrange), nBasis);
    save(sprintf('%s/NEU_%s_%s_plane%d.mat', ops.ResultsSavePath, ...
        ops.mouse_name, ops.date, ops.iplane),  'ops', 'S', 'Ntraces', '-v7.3')
end

author: cortex-lab


text:
function [tColor] = hyperTruecolor(rflFile, height, width, bands, rgbBands, strechtype)
% HYPERTRUECOLOR Returns a truecolor composite image, given a reflectance file
%   hyperTruecolor returns a truecolor composite image, given an AVIRIS reflectance file,
%  with an optional enhancement of the image using one of two contrast streches.
%
% Usage
%   [tColor] = hyperTruecolor(rflFile, height, width, bands)
%   [tColor] = hyperTruecolor(rflFile, height, width, bands, rgbBands)
%   [tColor] = hyperTruecolor(rflFile, height, width, bands, rgbBands, strechtype)
% Inputs
%   rflFile    - Path and filename of AVIRIS reflectance measurements
%   height     - Height in pixels
%   width      - Width in pixels
%   bands      - Number of bands
%   rgbBands   - [Red, green, blue] bands (1x3 vector)
%                [31,20,12] is the default if no argument is given
%   strechtype - truecolor enhancement contrast strech (string):
%                'none'         : default if no argument is given
%                'stretchlim'   : linear contrast stretch
%                'decorrstretch': decorrelation stretch
%                   (enhanced color separation across highly correlated channels)
% Outputs
%   tColor - Truecolor composite image

if nargin < 4 || nargin > 6
    help hyperTruecolor
    error('Incorrect usage, see function description above.')
end

if nargin == 4
    % RGB bands: [665.73, 557.07, 478.17] nm
    rgbBands   = [31,20,12];
    strechtype = 'none';
elseif nargin == 5
    strechtype = 'none';
end

h = height;
w = width;
p = bands;

% Read the 3 RGB bands
tColor = multibandread(rflFile, [h w p], 'int16', 0, 'bip', 'ieee-be', ...
            {'Row', 'Range', [1 h]}, {'Column', 'Range', [1 w]}, ...
            {'Band', 'Direct', rgbBands} );

% Normalize to proper reflectance units.
tColor = tColor ./ 1e4;

switch lower(strechtype)
    case 'none'
        return
    case 'stretchlim'
        % Apply a linear contrast stretch to the truecolor composite image
        tColor = imadjust(tColor, stretchlim(tColor));
    case 'decorrstretch'
        % Apply a decorrelation stretch to the truecolor composite image
        tColor = decorrstretch(tColor, 'Tol', 0.01);
    otherwise
        help hyperTruecolor
        error(sprintf('Unknown method %s. See function description above.', strechtype))
end

author: davidkun

repo: HyperSpectralToolbox


text:
function[data,h,hc]=provec(varargin)
%PROVEC  Generate progressive vector diagrams (simple and fancy).
%
%   Simple provecs:
%
%     PROVEC(DT,U,V) generates a simple progressive vector diagram plotting
%     CUMSUM(U*DT) vs CUMSUM(V*DT).  U and V are column vectors, or
%     matrices with time oriented in columns. DT is a scalar with units of
%     of hours, while U and V must have units of cm/s.
%
%     PROVEC(DT,CV), where CV=U+iV, also works.
%
%     INT=PROVEC(...) outputs the integrated dispacement INT in km.
%
%   Fancy provecs:
%
%     A fancy provec use SCATTER to plot the color and/or sizes of the
%     points according to another parameter, say density or temperature.
%     A colorbar is also plotted.
%
%     PROVEC(DT,U,V,C) uses C, of size SIZE(U), as the symbol color.
%
%     PROVEC(DT,U,V,C,S) also uses S, of size SIZE(U), as the symbol size.
%
%     PROVEC(DT,CV,C) and PROVEC(DT,CV,C,S) also work.
%
%     Note that since SCATTER is slow for large datasets, it is useful to
%     decimate the data after CUMSUMing but before plotting.  This is
%     accomplished using PROVEC(...,INDEX).  Then only the points
%     INT(INDEX,:) of the integrated trajectory are plotted.
%
%     [INT,H,HC]=PROVEC returns the intergrated displacement INT, the
%     handle H to the data, and the handle HC to the colorbar axis.
%
%   As an example,
%
%          load bravo94
%          th=100*detrend(bravo94.cat.th(:,3));
%          [int,h,hc]=provec(1,bravo94.rcm.cv(:,3),th,20+0*th,[1:10:4000]);
%          caxis([-8 8])
%
%   makes part of Figure 6b of Lilly and Rhines (2002) JPO.
%
%   Usage: int=provec(dt,cv);
%          [int,h,hc]=provec(dt,cv,c,index);
%          [int,h,hc]=provec(dt,cv,c,s,index);
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 1999--2017 J.M. Lilly --- type 'help jlab_license' for details

%          [int,h,hc]=provec(dcol,cv,c,s,index,deltat);
%

%   PROVEC(DCOL,...), where DCOL is a scalar or a row vector, specifies
%   that the columns of the data are to be offset by amount DCOL after
%   plotting, for both simple and fancy provecs.
%
%   If DCOL is a scalar, then successive columns after the first are offset
%   by amount -DCOL.  If DCOL is a row vector, then it specifies the offset
%   for each column of the data.

inunits=100;
outunits=1000;
factor=100*1000/3600;

bcolor=0;
defsize=5;

deltat=varargin{1};
varargin=varargin(2:end);
na=length(varargin);

%/********************************************************
% %Look for initial row vector
offs=0;
% if isrow(varargin{1}) || isscalar(varargin{1});
%   offs=varargin{1};
%   na=na-1;
%   varargin=varargin(2:end);
% end
% %\********************************************************

% if length(varargin{end})==1
%   deltat=varargin{end};
%   na=na-1;
% end

index=[];
if length(varargin{na})~=length(varargin{1})
  index=varargin{na};
  na=na-1;
end

c=[];
if isreal(varargin{1})
  data=varargin{1}+sqrt(-1)*varargin{2};
  if na>2
      c=varargin{3};
  end
  if na>3
      s=varargin{4};
  else
      s=defsize+zeros(size(data));
  end
else
  data=varargin{1};
  if na>1
      c=varargin{2};
  end
  if na>2
      s=varargin{3};
  else
      s=defsize+zeros(size(data));
  end
end

data=vswap(data,nan,0);
data=cumsum(data)*deltat/factor;
if ~(length(offs)==1&&allall(abs(offs)==0))
  if isscalar(offs)
    offs=-(0:1:size(data,2)-1)*offs;
  end
  if length(offs)~=size(data,2)
    error('Length of DCOL must equal number of columns of the data.')
  end
  for i=1:size(data,2)
      data(:,i)=data(:,i)+offs(i);
  end
end

if ~isempty(index)
  data=data(index,:);
  if ~isempty(c)
    c=c(index,:);
  end
  if ~isempty(s)
    s=s(index,:);
  end
end

if isempty(c)
	h=plot(data);
else
	bcolor=1;
        vcolon(data,s,c);
	if any(isnan(s))
	  error('S cannot contain NANs.')
	end
	if any(isnan(c))
	  error('C cannot contain NANs.')
	end
	h=scatter(real(data),imag(data),s,c,'filled');
end

xlabel('Displacement eastward (km)')
ylabel('Displacement northward (km)')

set(gca,'box','on')
set(gca,'dataaspectratio',[1 1 1])
pos=get(gca,'position');

%put a colorbar if needed
if bcolor
	ax=gca;
	hc=colorbar;
	%posc=get(hc,'position');
	%set(hc,'position',[posc(1) pos(2) posc(3) pos(4)])
	%the above doesn't take care of the relative size
	%problem--- instead, try to change aspect ratio
	axes(ax)
end

hold on

if nargout ==0
    clear data h hc
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [engine, loglik] = enter_evidence(engine, evidence, varargin)
% ENTER_EVIDENCE Add the specified evidence to the network (cond_gauss)
% [engine, loglik] = enter_evidence(engine, evidence, ...)
%
% evidence{i} = [] if if X(i) is hidden, and otherwise contains its observed value (scalar or column vector)

bnet = bnet_from_engine(engine);
ns = bnet.node_sizes(:);
observed = ~isemptycell(evidence);
onodes = find(observed);
hnodes = find(isemptycell(evidence));
engine.evidence = evidence;

% check there are no C->D links where C is hidden
pot_type = determine_pot_type(bnet, onodes);

dhid = myintersect(hnodes, bnet.dnodes);
S = prod(ns(dhid));
T = zeros(S,1);

N = length(bnet.dag);
mu = cell(1,N);
Sigma = cell(1,N);
cobs = myintersect(bnet.cnodes, onodes);
chid = myintersect(bnet.cnodes, hnodes);
ens = ns;
ens(cobs) = 0;
for j=chid(:)'
  mu{j} = zeros(ens(j), S);
  Sigma{j} = zeros(ens(j), ens(j), S);
end

for i=1:S
  dvals = ind2subv(ns(dhid), i);
  evidence(dhid) = num2cell(dvals);
  [sub_engine, loglik] = enter_evidence(engine.sub_engine, evidence);
  for j=chid(:)'
    m = marginal_nodes(sub_engine, j);
    mu{j}(:,i) = m.mu;
    Sigma{j}(:,:,i) = m.Sigma;
  end
  T(i) = exp(loglik);
end

[T, lik] = normalise(T);
loglik = log(lik);

engine.T = T;
engine.mu = mu;
engine.Sigma = Sigma;

dnodes = bnet.dnodes;
dobs = myintersect(dnodes, onodes);
ens(dobs) = 1;
engine.joint_dmarginal = dpot(dnodes, ens(dnodes), myreshape(engine.T, ens(dnodes)));

engine.onodes = onodes;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function onl_write_background(varargin)
% Periodically process data using a predictive model, and write results to some external device.
% onl_write_background(ResultWriter,StreamNames,Model,OutputFormat,UpdateFrequency,StartDelay,PredictorName,PredictAt,Verbose,EmptyResultValue)
%
% This is a convenience function which simplifies the definition of components which load and then
% periodically query a predictive model, in real time, and forward the results to some external
% device. The function is internally implemented using a timer that periodically triggers the
% computation of updated estimates, and their transfer to the data sink.
%
% In:
%   ResultWriter : Callback function that receives one or more BCI estimates and writes them to some
%                  external device. The format passed to the ResultWriter is according to OutputFormat.
%
%   StreamNames : optional names of stream data structures in the MATLAB base workspace to consider as
%                 possible data sources (previously created with onl_newstream); any stream that
%                 contains channels that are needed by the predictor will be linked to it (assuming
%                 that the choice of stream to use is not ambiguous).
%
%                 The identification of needed channels is primarily done on the basis of the channel
%                 labels -- if a stream has channels with labels that are required by a filter pipeline,
%                 it will be used as a source for this pipeline. The framework attempts to gracefully
%                 handle cases where a stream only provides a subset of the channels that were in the
%                 training set and the model only effectively operates on this subset via flt_selchans.
%
%   Model : A model data structure (as obtained from bci_train) according to which predictions shall be
%           made; typically this is a model struct, but for convenience it can be a file name,
%           variable name in the base workspace, or a cell array of {file name, variable name} to
%           refer to a variable inside a .mat file. The model is not modified by this function.
%           (default: 'lastmodel')
%
%   OutputFormat : Output data format, see onl_predict (default: 'distribution')
%
%   UpdateFrequency : Frequency at which the device should be queried, in Hz (default: 10)
%
%   StartDelay : Delay before real-time processing begins; grace period until user resources are
%                created (default: 1)
%
%   PredictorName : Name of a predictor data structure that shall be created in the MATLAB base
%                   workspace to hold the predictor's state. If a variable of this name already exists
%                   it will be overridden. (default: 'lastpredictor')
%
%   PredictAt : Predict at markers. If nonempty, this is a cell array of online target markers
%               relative to which predictions shall be made. If empty, predictions are always made
%               on the most recently added sample. (default: {})
%
%   Verbose : Verbose output. If false, the console output of the online pipeline will be suppressed.
%             (default: false)
%
%   EmptyResultValue : Empty-result value. This value is returned for predictions that yielded no
%                      result (e.g., due to an error or because not enough data was available).
%                      (default: NaN)
%
% Examples:
%   % after a predictive model has been learned using bci_train, and a data stream supplying raw
%   % data has been established, load the model into the online system and periodically send its
%   % outputs to a target destination
%   onl_write_background(@send_outputs_to_destination,'mystream')
%
%   % as before, but also specify a custom output format and a higher update frequency
%   onl_write_background(@send_outputs_to_destination,'mystream','lastmodel','expectation',25)
%
%   % as before, but pass all arguments by their short names
%   onl_write_background('ResultWriter',@send_outputs_to_destination,'MatlabStream','mystream','Model','lastmodel','OutputFormat','expectation','UpdateFrequency',25)
%
% See also:
%   onl_predict
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-01-18

% read options
arg_define(varargin, ...
    arg_norep({'result_writer','ResultWriter'},[],[],'Result writing callback. Callback function that receives one or more BCI estimates and writes them to some external device. The format passed to the ResultWriter is according to OutputFormat.'), ...
    arg({'in_stream','StreamName','MatlabStream'}, 'laststream',[],'Input stream name. This is the name of the stream data structure in the MATLAB base workspace that shall be read from. Can also be a cell array of stream names, if multiple, or empty if non-ambiguous.','typecheck',false,'shapecheck',false), ...
    arg({'pred_model','Model'}, 'lastmodel', [], 'Predictive model. A model data structure (as obtained from bci_train) according to which predictions shall be made; typically this is a model struct, but for convenience it can be a file name, variable name in the base workspace, or a cell array of {file name, variable name} to refer to a variable inside a .mat file. The model is not modified by this function.','type','expression'), ...
    arg({'out_form','OutputFormat'},'distribution',{'expectation','distribution','mode','raw'},'Format of the produced output values. Can be the expected value (posterior mean) of the target variable, or the distribution over possible target values (probabilities for each outcome, or parametric distribution), or the mode (most likely value) of the target variable.'), ...
    arg({'update_freq','UpdateFrequency'},10,[0 Inf],'Update frequency. This is the rate at which the outputs should be calculated.'), ...
    arg({'start_delay','StartDelay'}, 1, [0 Inf],'Start-up delay. Delay before real-time processing begins; grace period to initialize everything.'), ...
    arg({'pred_name','PredictorName'}, 'lastpredictor',[],'Name of new predictor. Name of a predictor data structure that shall be created in the MATLAB base workspace to hold the predictor''s state. If a variable of this name already exists it will be overridden.'), ...
    arg({'predict_at','PredictAt'}, {},[],'Predict at markers. If nonempty, this is a cell array of online target markers relative to which predictions shall be made. If empty, predictions are always made on the most recently added sample.','type','expression'), ...
    arg({'verbose','Verbose'}, false,[],'Verbose output. If false, the console output of the online pipeline will be suppressed.'), ...
    arg({'empty_result_value','EmptyResultValue'},NaN,[],'Empty-result value. This value is returned for predictions that yielded no result (e.g., due to an error or because not enough data was available).','type','expression'));

% input validation
if ischar(result_writer) && ~isempty(result_writer) %#ok<NODEF>
    if result_writer(1) ~= '@' && ~exist(result_writer,'file')
        error('The given ResultWriter argument (%s) is not a valid function name',result_writer); end
    result_writer = str2func(result_writer);
end
if ~isa(result_writer,'function_handle')
    error('The given ResultWriter argument must be a function handle (or function name), but was: %s',hlp_tostring(result_writer,10000)); end
if ~iscell(in_stream) %#ok<NODEF>
    in_stream = {in_stream}; end
for c=1:length(in_stream)
    stream_name = in_stream{c};
    if ~isvarname(stream_name)
        error('The given StreamName argument must be a valid variable name, but was: %s',hlp_tostring(stream_name,10000)); end
    try
        stream = evalin('base',stream_name);
    catch e
        error('Failed to look up stream named %s in MATLAB base workspace with error: %s',stream_name,e.message);
    end
    if ~isstruct(stream)
        error('The given data structure named %s in the MATLAB base workspace was expected to be a stream data structure, but was not a struct (wrong name?): %s',stream_name,hlp_tostring(stream,10000)); end
    if ~isfield(stream,'streamid')
        if isfield(stream,{'data','srate'})
            error('The given stream data structure named %s appears to be an EEGLAB data set struct but is not a stream (use onl_newstream to create a valid stream)',stream_name);
        else
            error('The given data structure named %s is not a valid stream (use onl_newstream to create a valid stream)',stream_name);
        end
    end
    streamids{c} = stream.streamid; %#ok<AGROW>
end

% create new predictor
predid = onl_newpredictor(pred_name,pred_model,in_stream,predict_at);

% create & start timer (which periodically writes to the stream)
start(timer('ExecutionMode','fixedRate', 'Name',[pred_name '_timer'], 'Period',1/update_freq, ...
    'StartDelay',start_delay, 'TimerFcn',@(timer_handle,varargin) write_data(pred_name,in_stream,out_form,result_writer,predid,streamids,timer_handle,verbose,empty_result_value)));

% background data writer
function write_data(predictor,stream_names,fmt,result_writer,pred_id,stream_ids,timer_handle,verbose,empty_result_value)
try
    % check if the stream and the predictor are still there
    for c=1:length(stream_names)
        s = evalin('base',stream_names{c});
        if s.streamid ~= stream_ids{c}
            error('Note: the stream named %s was recreated.',stream_names{c}); end
    end
    p = evalin('base',predictor);
    if p.predictorid ~= pred_id
        error('Note: the predictor named %s was recreated.',predictor); end
    % make a prediction
    y = onl_predict(predictor,fmt,~verbose,empty_result_value);
    % and write it out
    try
        result_writer(y);
    catch e
        disp('Error in result-writing function:');
        hlp_handleerror(e);
    end
catch e
    if ~strcmp(e.identifier,'MATLAB:UndefinedFunction')
        hlp_handleerror(e); end
    % stream or predictor have changed (e.g., replaced/deleted) --> stop timer
    stop(timer_handle);
    delete(timer_handle);
end

repo: REST


text:
function tanh_plot ( )

%*****************************************************************************80
%
%% TANH_PLOT makes a plot of certain tanh() functions.
%
%  Discussion:
%
%    The value is +1 if the number is positive or zero, and it is -1 otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2012
%
%  Author:
%
%    John Burkardt
%
  m = 101;
  x = ( linspace ( -1.0, +1.0, m ) )';

  u = zeros ( m, 6 );
  for j = 0 : 5
    s = 2^j;
    u(:,j+1) = tanh ( s * x(:) / 2.0 ) / ( tanh ( s * -1.0 / 2.0 ) );
  end

  plot ( x, u, 'Linewidth', 3 )
  grid on
  axis ( [-1, +1, -1.1, +1.1 ] )
  xlabel ( '<--- X --->' )
  ylabel ( '<--- U(X) ---> ' )
  title ( 'tanh ( 2^J * x / 2.0 ), J = 0, 1, ..., 5 scaled to be +1 at -1' )

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = year_length_common ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_COMMON returns the number of days in a Common year.
%
%  Discussion:
%
%    The "common" calendar is meant to be the calendar which is Julian up to
%    day JED = 2299160, and Gregorian from day JED = 2299161 and after.
%
%    If Y is 0, then the routine returns 0, reflecting the fact that
%    there was officially no year 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of
%    days in the year.
%
  if ( y == 0 )
    value = 0;
  elseif ( y == 1582 )
    value = 355;
  elseif ( year_is_leap_common ( y ) )
    value = 366;
  else
    value = 365;
  end

  return
end


text:
function ADEM_cost_SHC
% This demo illustrates the use of priors on the motion of hidden states as
% polices. It simulates exploration and exploitation using radial basis
% function attractors and auto-vitiative (self-destroying) attractors
% as the basis of the prior. These dynamics enforce exploration, under
% active inference. In turn, this exploration enables perceptual learning
% to associate attractors with changes in physiological states (cf,
% rewards). This can be exploited to by formal priors to ensure regions of
% physiological state-space are avoided.
% We look at this scheme using simulated pathology; first, we simulate a
% (neurodegenerative) reduction in log-precision (cf Parkinson's disease) on
% the motion of physical states.  This results in active inference with
% a loss of precise volitional movement and subsequent failure to optimise
% physiological states. Second, we look at the effects of precision on
% learning by increasing log-precision (cf, Addition) on the motion of
% physiological states. This results in a failure to learn and, again,
% subsequent failure to optimise physiological states.
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: ADEM_cost_SHC.m 7679 2019-10-24 15:54:07Z spm $

% switch for demo
%--------------------------------------------------------------------------
DEMO = 1;

% location and radius of attractors (A)
%--------------------------------------------------------------------------
global A
A.x  = [1 -1 -1  1;
        2 -2  2 -2];
A.d  = 1/4;              % s.d. of Gaussian radial basis function
A.u  = 1/8;              % threshold to induce rapid decay of x.a
A.q  = [1 2];            % indices of locations that increase x.q

% generative process
%==========================================================================

% level 1
%--------------------------------------------------------------------------
G(1).x.x = [0;0];
G(1).x.v = [1;0];
G(1).x.q = [1;1/2];

G(1).f   = inline('spm_cost_SHC_fxa(x,v,a,P)','x','v','a','P');
G(1).g   = inline('[x.x; x.v; x.q]','x','v','a','P');
G(1).V   = exp(16);                          % error precision
G(1).W   = exp(16);                          % error precision
G(1).U   = exp(4);                           % action precision

% level 2
%--------------------------------------------------------------------------
G(2).a   = sparse(2,1);                      % action
G(2).v   = 0;                                % inputs
G(2).V   = exp(16);
G        = spm_ADEM_M_set(G);

% generative model
%==========================================================================

% level 1
%--------------------------------------------------------------------------
M(1).x.x = G(1).x.x;
M(1).x.v = G(1).x.v;
M(1).x.q = G(1).x.q;
M(1).x.a = randn(4,1)/8;

M(1).f   = inline('spm_cost_SHC_fx(x,v,P)','x','v','P');
M(1).g   = inline('[x.x; x.v; x.q]','x','v','P');
M(1).pE  = speye(4,2);

M(1).V   = exp(8);
M(1).W   = diag(exp([[1 1 1 1]*4 [1 1]*6 [1 1 1 1]*4]));

% level 2 (no exogenous inputs in this simulation)
%--------------------------------------------------------------------------
M(2).v   = 0;
M(2).V   = exp(16);

% Integrate: active inference
%==========================================================================
M(1).E.nE = 1;
M(1).E.n  = 4;

N      = 128;
DEM.U  = sparse(N,1);
DEM.C  = sparse(N,1);
DEM.G  = G;
DEM.M  = M;

if DEMO
    load DEM_addiction
else
    DEM = spm_ADEM(DEM);
end

% show behavior
%==========================================================================

% overview
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1');
spm_DEM_qU(DEM.qU)

subplot(2,2,3)
spm_cost_SHC_path(DEM.pU,A)

% a closer look at physiology
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 2'); clf

subplot(2,1,1)
t    = 1:N;
plot(t,DEM.qU.x{1}(5:6,:)',t,t*0 + A.u,':'), hold on
plot(t,DEM.qU.x{1}(7:10,:)','-.'),           hold off
xlabel('time','Fontsize',14)
title('internal states','FontSize',16)
axis square, box off, set(gca,'XLim',[1 N])

subplot(2,2,3)
plot(t,DEM.qU.x{1}(5:6,:)',t,t*0 + A.u,':')
xlabel('time','Fontsize',14)
title('physiological states','FontSize',16)
axis square, box off, set(gca,'XLim',[1 N])

subplot(2,2,4)
plot(DEM.pU.x{1}(5,:),DEM.pU.x{1}(6,:)),               hold on
plot([A.u 1],[A.u A.u],'-.r',[A.u A.u],[A.u 1],'-.r'), hold off
xlabel('time','Fontsize',14)
title('physiological states','FontSize',16)
axis square, box off, axis([-.1 1.2 -.1 1.2])

% look at the effect of precision on inference (cf Parkinson's disease)
%==========================================================================

% initialize states
%--------------------------------------------------------------------------
DEM = spm_ADEM_update(DEM);

% simulate exposure with different log-precisions on physical motion
%--------------------------------------------------------------------------
WP    = [4 2 0];
if ~DEMO
    for i = 1:length(WP)
        DEM_P{i}        = DEM;
        W               = [[1 1 1 1]*WP(i) [1 1]*6 [1 1 1 1]*4];
        DEM_P{i}.M(1).W = diag(exp(W));
        DEM_P{i}        = spm_ADEM(DEM_P{i});
    end
end

% Graphics - path and physiology
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 3'); clf
for i = 1:3

    % path
    %----------------------------------------------------------------------
    spm_figure('GetWin','Figure 3');
    subplot(3,2,i*2 - 1)
    spm_cost_SHC_path(DEM_P{i}.pU,A)

    % and physiology
    %----------------------------------------------------------------------
    subplot(3,2,i*2)
    plot(t,DEM_P{i}.pU.x{1}(5:6,:)',t,t*0 + A.u,':'), hold on
    plot(t,sum(DEM_P{i}.pU.x{1}(1:2,:).^2)/8,'m:'),    hold off
    xlabel('time','Fontsize',12)
    title(sprintf('%s (%1.0f)','physiological states',WP(i)),'FontSize',16)
    axis square, box off, set(gca,'XLim',[1 N])
    drawnow

end

% Graphics - action and prediction error
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 4'); clf
for i = 1:3

    % path
    %----------------------------------------------------------------------
    subplot(3,2,2*i - 1)
    plot(t,DEM_P{i}.qU.a{2})
    xlabel('time','FontSize',12)
    title('action','FontSize',16)
    axis square, box off, set(gca,'XLim',[1 N])
    if i == 1; a = axis; end, axis(a)

    % and prediction error
    %----------------------------------------------------------------------
    subplot(3,2,2*i)
    plot(t,DEM_P{i}.qU.z{1}(1:4,:))
    xlabel('time','FontSize',12)
    title('sensory error','FontSize',16)
    axis square, box off, set(gca,'XLim',[1 N])
    if i == 1; aa = axis; end, axis(aa)
    drawnow
end

% Learning (at different levels of precision on physiological motion)
%==========================================================================

% Switch attractor for second physiological state and log-precision levels
%--------------------------------------------------------------------------
spm_figure('GetWin','DEM'); clf
A.q  = [1 3];
WA   = [4 8 12];

if ~DEMO

    for i = 1:length(WA)

        % Enable learning
        %------------------------------------------------------------------
        A.u              = 0;
        DEM_L{i}         = DEM;
        DEM_L{i}.M(1).pC = exp(8);
        W                = [[1 1 1 1]*4 [1 1]*WA(i) [1 1 1 1]*4];
        DEM_L{i}.M(1).W  = diag(exp(W));
        DEM_L{i}         = spm_ADEM(DEM_L{i});

        % replace prior with posterior and re-expose
        %------------------------------------------------------------------
        A.u              = 1/8;
        DEM_D{i}         = DEM_L{i};
        DEM_D{i}.M(1).pE = DEM_D{i}.qP.P{1};
        DEM_D{i}.M(1).pC = [];
        DEM_D{i}         = spm_ADEM(DEM_D{i});

    end

    % save DEM structures for future demonstrations
    %----------------------------------------------------------------------
    save DEM_addiction DEM DEM_P DEM_L DEM_D

end

% Graphics - optimal learning
%--------------------------------------------------------------------------
spm_figure('GetWin','DEM'); clf
spm_DEM_qU(DEM.qU,DEM.pU)

spm_figure('GetWin','Figure 5'); clf
spm_DEM_qP(DEM_L{1}.qP)

subplot(2,2,3)
spm_cost_SHC_path(DEM.pU,A)
title('Before','Fontsize',16)
subplot(2,2,4)
spm_cost_SHC_path(DEM_D{1}.pU,A)
title('After','Fontsize',16)
drawnow

% Comparison of learning over levels of log-precisions (cf. Addiction)
%==========================================================================
spm_figure('GetWin','Figure 6'); clf

for i = 1:3

    % path
    %----------------------------------------------------------------------
    subplot(3,2,i*2 - 1)
    spm_cost_SHC_path(DEM_D{i}.pU,A)
    title('After','FontSize',16)

    % and physiology
    %----------------------------------------------------------------------
    subplot(3,2,i*2)
    spm_plot_ci(DEM_L{i}.qP.P{1}(:),DEM_L{i}.qP.C)
    xlabel('parameter','FontSize',12)
    title(sprintf('%s (%1.0f)','learning',WA(i)),'FontSize',16)
    axis square, box off

end
drawnow

% and physiology
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 7'); clf
for i = 1:3

    % predicted motion of second physiological state
    %----------------------------------------------------------------------
    for j = 1:N
        x    = DEM_L{i}.qU.x{1}(:,j);
        f    = spm_cost_SHC_fx(spm_unvec(x,M(1).x),M(2).v,M(1).pE);
        p(j) = DEM_L{i}.qU.w{1}(6,j) + f.q(2);
        x    = DEM_L{i}.pU.x{1}(:,j);
        f    = spm_cost_SHC_fxa(spm_unvec(x,G(1).x),M(2).v,G(2).a,G(1).pE);
        q(j) = f.q(2);
    end
    subplot(3,2,2*i - 1)
    plot(t,p,t,q,'b:')
    xlabel('time','FontSize',12)
    title(sprintf('%s (%1.0f)','predicted motion',WA(i)),'FontSize',16)
    axis square, box off, set(gca,'XLim',[1 N])
    if i == 1; a = axis; end, axis(a)

    subplot(3,2,2*i)
    plot(t,DEM_L{i}.qU.w{1}(6,:))
    xlabel('time','FontSize',12)
    title('prediction error','FontSize',16)
    axis square, box off, set(gca,'XLim',[1 N])
    if i == 1; aa = axis; end, axis(aa)

end

author: spm


text:
function truth = isexpandable2row(obj, n)

truth = isscalar(obj) || (isvector(obj) && length(obj)==n);

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd