text:
function [Fnew,Qnew,Hnew,Rnew,Bnew,Params] = KalmanFilterX_LearnEM_Mstep_new2()

    if(nargin<5)
        Bnew = 0;
    end

    y       = cell2mat(cellfun(@(x)x.y,est_f,'un',0));    % Observations {N} x (yDim x 1)
    u       = cell2mat(cellfun(@(x)x.u,est_f,'un',0));    % Control inputs {N{ x (uDim x 1)
    V_f     = cellfun(@(x)x.P,est_f,'un',0);    % Filtered co-variences {N} x (xDim x xDim)
    x_s     = cell2mat(cellfun(@(x)x.x,est_s,'un',0));    % Smoothed estimates {N} x (xDim x 1)
    V_s     = cellfun(@(x)x.P,est_s,'un',0);    % Smoothed co-variences {N} x (xDim x xDim)
    K       = cellfun(@(x)x.K,est_f,'un',0);    % Filtered Kalman gains {N} x (xDim x yDim)
    [xDim,N] = size(x_s);
    est_s{N}.C = est_s{N-1}.C;
    C       = cellfun(@(x)x.C,est_s,'un',0);    % Smoothed Kalman gains {N} x (xDim x yDim)

    Vt_tm1  = zeros(xDim,xDim,N);     % V_{t,t-1}
    Pt      = zeros(xDim,xDim,N);     % P_{t}
    Pt_tm1  = zeros(xDim,xDim,N);     % P_{t,t-1}

    Vt_tm1(:,:,N) = F*V_f{N-1} - K{N}*H*F*V_f{N-1} ;
    Pt_tm1(:,:,N) = Vt_tm1(:,:,N) + x_s(:,N)*x_s(:,N-1)';
    Pt(:,:,N)     = V_s{N} + x_s(:,N)*x_s(:,N)';
    Pt(:,:,N-1)     = V_s{N-1} + x_s(:,N-1)*x_s(:,N-1)';

    Params.loglik = computeEMLoglik(est_f{1}.x, est_f{1}.P, x_s,y,F,H,Q,R,B,u);

    for k = N-1:-1:2
        Vt_tm1(:,:,k) = V_f{k}*C{k-1}' + C{k}*(Vt_tm1(:,:,k+1) - F*V_f{k})*C{k-1}';
        Pt_tm1(:,:,k) = Vt_tm1(:,:,k) + x_s(:,k)*x_s(:,k-1)';
        Pt(:,:,k-1)   = V_s{k-1} + x_s(:,k-1)*x_s(:,k-1)';
    end

%     %% Backward Rauch Recursion %Equation 2
%     [xhat, VtT, J, Pt2, Vtt1T, Ptt1] = Shuang_KF_Smooth(cell2mat(cellfun(@(x)x.x,est_f,'un',0)),cell2mat(cellfun(@(x)x.P,est_f,'un',0)),...
%             cell2mat(cellfun(@(x)x.PPred,est_f,'un',0)), F, H, B, [u(2:end),0], K{N});
%     sum(xhat == x_s)
%     test = permute(Pt_tm1,[2 3 1]);
%     sum(Ptt1 == test)
%     test = permute(Vt_tm1,[2 3 1]);
%     %sum(Vtt1T == test)
%     sum(Vtt1T - test<eps)
%     sum(cell2mat(C(1,1:end-1)) == J)
%     %cell2mat(C(1,1:end-1)) == J
%     test = permute(Pt,[2 3 1]);
%     sum(Pt2 == test)

%     Hnew = sum(y(:,1:N).*xhat(:,1:N))/sum(Pt2(1:N));
%
%     %Update noise covariance
%     Rnew = 1/N*sum(y.*y - Hnew*xhat.*y);
%
%     %Update dynamics matrix
%     Fnew = sum(Ptt1(2:N)-B*u(2:N).*xhat(1:N-1))/sum(Pt2(1:N-1));
%
%     %new version of B and Q
%     Bnew = sum((xhat(2:N)-F*xhat(1:N-1)).*u(2:N))/sum((u(2:N).^2));
%
%     %Update state noise covariance
%     Qnew = 1/(N-1)*sum(Pt2(2:N) + xhat(2:N).*u(2:N)*B - F*Ptt1(2:N) + F*xhat(1:N-1).*u(2:N)*B - B*u(2:N).*xhat(2:N) + B*(u(2:N).^2)*B);

    %Update output matrix
    Hnew = y*x_s'/sum(Pt,3);

    %Update noise covariance
    Rnew = 1/N*( y*y' - y*(x_s'*H') - (H*x_s)*y' + H*sum(Pt,3)*H'); %sum( - H*xhat.*y);

    %Update dynamics matrix
    Fnew = (sum(Pt_tm1(:,:,2:N),3)-(B*u(:,2:N))*x_s(:,1:N-1)') / sum(Pt(:,:,1:N-1),3);

    %new version of B and Q
    Bnew = (x_s(:,2:N)*u(:,2:end)' - F*x_s(:,1:N-1)*u(:,2:N)') / (u(:,2:N)*u(:,2:N)');

    %Update state noise covariance
    Qnew = 1/(N-1)*(sum(Pt(:,:,2:end),3) - sum(Pt_tm1(:,:,2:end),3)*F' - F*sum(Pt_tm1(:,:,2:end),3)' + F*sum(Pt(:,:,1:N-1),3)*F' ...
         + F*(x_s(:,1:N-1)*u(:,2:N)')*B' + B*(u(:,2:end)*x_s(:,1:N-1)')*F'- x_s(:,2:N)*u(:,2:N)'*B' - B*(u(:,2:N)*x_s(:,2:N)')+ B*(u(:,2:N)*u(:,2:N)')*B');

    Params.Vt_tm1 = Vt_tm1;
    Params.Pt_tm1 = Pt_tm1;
    Params.Pt = Pt;

end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [W, H, w, h, x, y] = floorplan(adj_H, adj_V, rho, Amin, l, u )

% Computes a minimum-perimeter bounding box subject to positioning constraints
%
% Inputs:
%      adj_H,adj_V: adjacency matrices
%      Amin:        minimum spacing: w_i * h_i >= Amin
%      rho:         boundaries: rho <= x_i <= W-rho, rho <= y_i <= H-rho
%      l, u:        aspect ratio constraints: l_i <= h_i/w_i <= u_i
% Only adj_H and adj_V are required; the rest are optional. If n is the
% number of cells, then adj_H and adj_V must be nxn matrices, and Amin,
% l, and u must be vectors of length n. rho must be a scalar. The default
% values of rho and Amin are 0.
% Joelle Skaf - 12/04/05

if nargin < 2
    error('Insufficient number of input arguments');
end

[n1, n2] = size(adj_H);
[m1, m2] = size(adj_V);

if n1~=n2
    error('Input adjacency matrix for horizontal graph must be square');
end

if m1~=m2
    error('Input adjacency matrix for horizontal graph must be square');
end

if n1~=m1
    error('Input adjacency matrices must be of the same size');
end

n = n1;                     % number of cells

if nargin <3
    rho = 0;
end

if nargin <4
    Amin = zeros(1,n);
else
    if min(size(Amin)) ~=1
        error('Amin should be a vector');
    end
    if max(size(Amin)) ~= n
        error('Amin should have the same length as the input graphs');
    end
    if size(Amin,1)~=1
        Amin = Amin';
    end
end

if nargin == 5
    if min(size(1)) ~= 1
        error('l must be a vector');
    end
    if max(size(l)) ~= n
        error('the vector l must have same length as the input graphs');
    end
    if size(l,1) == 1
        l = l';
    end
end

if nargin == 6
    if min(size(1)) ~= 1
        error('u must be a vector');
    end
    if max(size(u)) ~= n
        error('the vector u must have same length as the input graphs');
    end
    if size(u,1) == 1
        u = u';
    end
end

if nargin < 6
    u = [];
end
if nargin < 5
    l = [];
end

% verifying that there is a directed path between any pair of cells in at
% least one of the 2 graphs

paths_H = adj_H;
paths_V = adj_V;
temp_H = adj_H^2;
temp_V = adj_V^2;
while (sum(temp_H(:))>0)
    paths_H = paths_H + temp_H;
    temp_H = temp_H*adj_H;
end
while (sum(temp_V(:))>0)
    paths_V = paths_V + temp_V;
    temp_V = temp_V*adj_V;
end

hh = paths_H + paths_H';
vv = paths_V + paths_V';
p = hh+vv+eye(n);
all_paths = p>0;
if sum(all_paths(:)) ~= n^2
    error('There must be a directed graph between every pair of cells in one or the other input graphs');
end

par_H = sum(adj_H,2);               % number of parents of each node in H
par_V = sum(adj_V,2);               % number of parents of each node in V
chi_H = sum(adj_H);                 % number of children of each node in H
chi_V = sum(adj_V);                 % number of children of each node in V

% find the root(s) for each tree
roots_H = find(par_H==0);
roots_V = find(par_V==0);

% find all non-root nodes for each tree
nodes_H = find(par_H>0);
nodes_V = find(par_V>0);

% find leaf(s) for each tree
leafs_H = find(chi_H==0);
leafs_V = find(chi_V==0);

cvx_begin quiet
        variables x(n) y(n) w(n) h(n) W H
        minimize ( W + H )
        w >= 0;
        h >= 0;
        x(leafs_H) >= rho;
        y(leafs_V) >= rho;
        x(roots_H) + w(roots_H) + rho <= W;
        y(roots_V) + h(roots_V) + rho <= H;
        for i=1:length(nodes_H)
            node = nodes_H(i);
            c = adj_H(node,:);
            prnt = find(c>0)';
            m = length(prnt);
            x(node) + w(node) + rho <= x(prnt);
        end

        for i=1:length(nodes_V)
            node = nodes_V(i);
            c = adj_V(node,:);
            prnt = find(c>0)';
            m = length(prnt);
            y(node) + h(node) + rho <= y(prnt);
        end

        if sum(size(u))~= 0
            h <= u.*w;
        end
        if sum(size(l))~= 0
            h >= l.*w;
        end
        w' >= quad_over_lin([Amin.^.5;zeros(1,n)],h');
cvx_end

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% STD_GETTRIALSIND - return the index of the trials that comply with the
%                      defined values from trialinfo
%
% Usage:
%   >>  trialindsx = std_gettrialsind('testnewstudyfile.mat', 'type', { 'C' 'X' 'Y' }, ...
%                                     'load', [0 1 2 3 4 5 6], 'uncertainty1',[1 2])
%   >>  trialindsx = std_gettrialsind('testnewstudyfile.mat', 'load', [0 1 2 3 4 5 6])
%   >>  [trialindsx cont] = std_gettrialsind('testnewstudyfile.mat', 'load', [0 1 2 3 4 5 6], ...
%                                      'duration', '0<100', 'rt', '')
%
% Inputs:
% filename - Filename of the study file, the structure contained in that
% file or, the trialinfo structure for that specific subject
% Beware of the value type: cell or string is categorical string variables,
% numerical array is categorical numerical variables, empty string or string
% containing "<" character is for continuous variables. For continuous variables,
% a range may be entered. Enter empty for the whole value range.
%
% varargin - Pairs of inputs ,varname(string),varvalues . (see Usage)
%
% Optional inputs:
%   ['var']    - [cell or array] ['var'] is the name of one independent
%                variable followed by a value or a list of value for this
%                variable. To obtain values for a continuous variable,
%                use '' for all the values or a range '0<100' to get values
%                within a given range (see example below).
%
% Outputs:
%   trialind   - trial indices
%   cont       - optional continuous variable value for these indices
%
% Example:
%   >>  trialindsx = std_gettrialsind('testnewstudyfile.mat', 'type', { 'C' 'X' 'Y' }, ...
%                                     'load', [0 1 2 3 4 5 6], 'uncertainty1',[1 2])
%   >>  trialindsx = std_gettrialsind('testnewstudyfile.mat', 'load', [0 1 2 3 4 5 6])
%   >>  [trialindsx cont] = std_gettrialsind('testnewstudyfile.mat', 'load', [0 1 2 3 4 5 6], ...
%                                      'duration', '0<100', 'rt', '')
% Authors: Arnaud Delorme
%          Ramon Martinez-Cancino
%
% Copyright (C) 2015  Ramon Martinez-Cancino, UCSD, INC, SCCN
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [trialindsx, eventvals] = std_gettrialsind(filename,varargin)

trialindsx = [];
eventvals = [];
if nargin < 1
    help std_gettrialsind;
    return;
end

% Input stuff
try
    if length(varargin) == 1
        varargin = varargin{:}; % Call from std_readfile or other func
    end
    varsin = varargin;
    if ~isempty( varargin )
        for i = 1:2:numel(varsin)
            queryvars.(varsin{i}) = varsin{i+1};
        end
    else queryvars = []; end
catch
    disp('std_gettrialsind() error: calling convention {''key'', value, ... } error'); return;
end

% Checking if first entry is string or struct
if isstruct(filename)
    if isfield(filename,'trialinfo')
        trialinfo = filename.trialinfo;
    else
        trialinfo = filename;
        %warning('std_gettrialsind: First input assumed as trialinfo structure');
    end
else
    % Loading file
    trialinfo = load(filename,'-mat', 'trialinfo');
    trialinfo = trialinfo.trialinfo;
end

% --- Query ---
varnames = fieldnames(queryvars);
hits = zeros(length(trialinfo),length(varnames));

for iVar = 1 :  length(varnames)
    dattrials = {trialinfo.(varnames{iVar})};
    indvarvals = queryvars.(varnames{iVar});

    % case of continuous variable
    if ischar(indvarvals) && (isempty(indvarvals) || any(indvarvals == '<'))
        if ischar(dattrials{1})
            error(sprintf('Type error - excepting string numerical for field %s in data file', varnames{iVar}));
        end
        if isempty(indvarvals)
            hits(:,iVar) = ~cellfun(@isempty, dattrials);
        else % test for range
            indGreater = find(indvarvals == '<');
            if length(indGreater) ~= 1, error('Undefined range for continuous var'); end
            if any(indvarvals == '='), error('<= or =< not allowed for continuous var'); end
            lowerBound = str2num(indvarvals(1:indGreater-1));
            upperBound = str2num(indvarvals(indGreater+1:end));
            hits(:,iVar) = cellfun(@(x)(~isempty(x)&&x>lowerBound&&x<upperBound), dattrials);
        end
        dattrials(~hits(:,iVar)) = { NaN };
        eventvals(end+1,:) = [ dattrials{:} ];
    else
        % case of categorical variable
        if ischar(dattrials{1})
            if ischar(indvarvals) indvarvals = { indvarvals }; end
            if ~iscell(indvarvals)
                error(sprintf('Type error - excepting numerical values for field %s', varnames{iVar}));
            end
            for iVal = 1:length(indvarvals) % programmed for speed - AD
                hits(strmatch(indvarvals{iVal}, dattrials, 'exact'),iVar) = 1;
                % older version not compatible with old Matlab
                % hits(:,iVar) = hits(:,iVar) || strncmp(indvarvals{iVal}, dattrials, 100)';
            end
        else
            if iscell(indvarvals)
                error(sprintf('Type error - excepting string values for field %s', varnames{iVar}));
            end
            excl = find(cellfun(@isempty, dattrials));
            if ~(isempty(excl))
                [dattrials{excl}] = deal(nan);
            end
            dattrials = [ dattrials{:} ];
            for iVal = 1:length(indvarvals) % programmed for speed - AD
                hits(:,iVar) = hits(:,iVar) | [ dattrials == indvarvals(iVal) ]';
            end
        end
    end
end

% Retrieving overlapped values
hits = sum(hits,2);
trialindsx = find(hits == length(varnames));
if ~isempty(eventvals)
    eventvals  = eventvals(:,trialindsx);
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [intEst,absErrEst,exitCode]=integral1DAdaptiveCC(func,bounds,RelTol,AbsTol,NPowMin,NPowMax)
%%INTEGRAL1DADAPTIVECC Perform adaptive numerical integration over a
%           function that takes a scalar input and returns a scalar or
%           vector quantity. The integration is performed using
%           Clenshaw-Curtis points. The function ClenshawCurtisPoints1D is
%           used to get points of two orders and the difference in the
%           integral values from each is the error estimate. If the error
%           is too high, then the next higher pair of orders is used (the
%           next one being double the previous order). Note
%           that Clenshaw-Curtis points of length 2*N have all the
%           quadrature points of length N embedded, so by using powers of
%           2 for the lengths, we can reuse previously computed function
%           values.
%
%INPUTS: f The handle to the function. f(x) must be able to take a 1XN
%          vector of points and return the numDimXN set of function values
%          at those points.
%   bounds A 2X1 or 1X2 vector such that bounds(1) is the lower bound of
%          integration and bounds(2) is the upper bound of integration. it
%          is not required that bounds(2)>bounds(1). Integration bounds
%          must be finite. The default if omitted or an empty matrix is
%          passed is [-1;1].
%   RelTol The maximum relative error tolerance allowed. If omitted or an
%          empty matrix is passed, the default value of 1e-8 is used. When
%          f returns a vector, this is applied to each element separately.
%   AbsTol The absolute error tolerance allowed. If omitted or an empty
%          matrix is passed, the default value of 1e-11 is used. When
%          f returns a vector, this is applied to each element separately.
%  NPowMin The minimum order of the Chebyshev polynomials to use is
%          2^NPowMin. NPowMin is an integer>=1. The default if omitted or
%          an empty matrix is passed is 2.
%  NPowMax The maximum order of the Chebyshev polynomials that will be
%          considered is 2^NPowMax. The default if omitted or
%          an empty matrix is passed is 11, referring to order 2048.
%          NPowMax must be greater than NPowMin.
%
%OUTPUTS: intEst The estimate of the integral.
%      absErrEst An estimate of the error of the integral. This is the
%                difference of the values at two different orders.
%       exitCode A value indicating the status of the algorithm on
%                termination. Possible values are
%                0 Termination occurred due to the absolute or relative
%                  error tolerances being fulfilled.
%                1 Termination occurred due to reaching NPowMax.
%                2 A singularity of a NaN was encountered during a step;
%                  the returned values of intEst and totalError are from
%                  the previous step and do not meet the error criterion.
%                3 A singularity or NaN was encountered. The values of
%                  intEst and totalError are most likely invalid.
%
%The initial Chebyshev interpolation attempted will be order 2^NPowMin and
%an approximation of order 2^(NPowMin+1) will be used to get the error
%estimate. After that, the order of the pairs attempted is increased until
%convergence or the last error estimate was obtained having order NPowMax
%without convergence.
%
%EXAMPLE:
%Here, we compare the integral of a nonlinear function that is rather peaky
%in the beginning with the Chebyshev approximated integral.
% f=@(x)(x.*exp(-x.^2-x/2));
% bounds=[0;16];
% [intEst,absErrEst,exitCode]=integral1DAdaptiveCC(f,bounds);
% explicitSol=(1/8)*(4-4/exp(264)+exp(1/16)*sqrt(pi)*(erf(1/4)-erf(65/4)));
% relErr=abs(intEst-explicitSol)/abs(explicitSol)
%The relative error will be on the order of 1e-16, indicating a good
%agreement between the explicit an numerically determined integrals.
%
%May 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(NPowMax))
    NPowMax=11;
end

if(nargin<5||isempty(NPowMin))
    NPowMin=2;
end

if(NPowMax<=NPowMin)
    error('It is required that NPowMax>NPowMin.');
end

if(nargin<4||isempty(AbsTol))
    AbsTol=1e-11;
end

if(nargin<3||isempty(RelTol))
    RelTol=1e-8;
end

if(any(~isfinite(bounds(:))))
    error('The bounds must be finite. Perform a change of variables to handle infinite bounds.');
end

if(nargin<2||isempty(bounds)||(bounds(1)==-1&&bounds(2)==1))
    %Integration bounds are -1 to 1.
    f=func;
    df=1;
else
    %Integration bounds are given in bounds. We have to map them to (-1,1).
    LB=bounds(1);
    UB=bounds(2);
    f=@(x)(func(((UB-LB)/2).*x+(UB+LB)/2));
    df=(UB-LB)/2;
end
NCur=NPowMin;
[xi,w]=ClenshawCurtisPoints1D(2^NCur);
fLo=f(xi);
numDim=size(fLo,1);
%The integral value with the lower polynomial order approximation.
ILow=sum(fLo.*w',2);

%If a non-finite term was encountered.
if(~isfinite(ILow))
    absErrEst=Inf;
    exitCode=3;
    intEst=ILow*df;
    return;
end

absErrEstPrev=Inf;
NCur=NCur+1;
while(NCur<=NPowMax)
    highPow=2^NCur;
    [xi,w]=ClenshawCurtisPoints1D(2^NCur);
    %The Clenshaw-Curtis points are nested, so we do not need to reevaluate
    %the function at the points it was already evaluated at (every other
    %point).
    numPts=highPow+1;
    fHigh=zeros(numDim,numPts);
    fHigh(:,1+(0:2:highPow))=fLo;
    fHigh(:,1+(1:2:(highPow-1)))=f(xi(1+(1:2:(highPow-1))));

    IHigh=sum(fHigh.*w',2);

    %If a non-finite term was encountered.
    if(~isfinite(ILow))
        absErrEst=absErrEstPrev;
        exitCode=2;
        intEst=ILow*df;
        return;
    end

    %Check convergence
    absErrEst=abs(IHigh-ILow);
    if(all(absErrEst<=AbsTol)||all(absErrEst<=RelTol*abs(IHigh)))
        exitCode=0;
        intEst=IHigh*df;
        return
    end

    absErrEstPrev=absErrEst;
    ILow=IHigh;
    fLo=fHigh;
    NCur=NCur+1;
end

%If the maximum power was reached.
exitCode=1;
intEst=IHigh*df;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function nii=jnii2nii(jnii, varargin)
%
%    nii=jnii2nii(jnii)
%       or
%    nii=jnii2nii(jniifile)
%    jnii2nii(jniifile, niifile)
%
%    Covert a JNIfTI file or data structure to a NIfTI-1/2 structure or file
%
%    This function is compatible with both MATLAB and GNU Octave.
%    It accepts .jnii and .bnii input files
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        jnii: a JNIfTI data structure (a struct with NIFTIHeader and NIFTIData fields);
%              if jnii is a string, it represents a JNIfTI file (.jnii/.bnii)
%        niifile: if the 2nd parameter is given as a file name, the converted nifti data
%              will be save as a nii file with filename specified by niifile.
%              if the filename in niifile contains .gz, the file will be compressed using
%              the zmat toolbox.
%
%    output:
%        nii: is the converted nifti-1/2 data structure, it contains the below subfields
%          nii.img: the data volume read from the nii file
%          nii.hdr: extended raw file header, a structure that is byte-wise compatible with a
%                   nifti-1 - in this case, typecast(nii.hdr,'uint8') must be 348+4=352 bytes,
%                       including the raw nifti-1 hdr header (348 bytes) plus the 4-byte
%                       extension flags), or
%                   nifti-2 - in this case, typecast(nii.hdr,'uint8') must be 540+4=544 bytes,
%                       including the raw nifti-2 hdr header (540 bytes) plus the 4-byte
%                       extension flags)
%              if one run nii.hdr.extension=[]; the resulting struct is 348/540-byte in length
%              nii.hdr key subfileds include
%
%              sizeof_hdr: must be 348 (for NIFTI-1) or 540 (for NIFTI-2)
%              dim: short array, dim(2: dim(1)+1) defines the array size
%              datatype: the type of data stored in each voxel
%              bitpix: total bits per voxel
%              magic: must be 'ni1\0' or 'n+1\0' for NIFTI-1 data, and 'ni2\0' or 'n+2\0' for NIFTI-2 data
%
%              For the detailed nii header, please see
%              https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h
%
%    dependency:
%
%          To load a JNIfTI file with compression or niifile ends with (.nii.gz/.hdr.gz/.img.gz),
%          one must install the ZMat Toolbox (http://github.com/fangq/zmat) and
%          JSONLab Toolbox (http://github.com/fangq/jsonlab);
%
%    this file is part of JNIfTI specification: https://github.com/fangq/jnifti
%
%    License: Apache 2.0, see https://github.com/fangq/jnifti for details
%

if(nargin<=0)
    help jnii2nii
    return;
end

if(~isstruct(jnii))
    jnii=loadjnifti(jnii);
end

if(~(isfield(jnii,'NIFTIHeader') && isfield(jnii,'NIFTIData')))
    error('input must be a valid JNIfTI structure (needs both NIFTIHeader and NIFTIData subfields)');
end

niiformat='nifti1';

if((isfield(jnii.NIFTIHeader,'NIIFormat') && ismember(jnii.NIFTIHeader.NIIFormat(1:3),{'ni2','n+2'})) || max(jnii.NIFTIHeader.Dim)>=2^32)
    niiformat='nifti2';
end

nii.hdr=nifticreate(jnii.NIFTIData, niiformat);
nii.img=jnii.NIFTIData;

if(isfield(jnii.NIFTIHeader,'NIIHeaderSize'))
    nii.hdr.sizeof_hdr	=bytematch(jnii.NIFTIHeader, 'NIIHeaderSize', nii.hdr.sizeof_hdr);
end
if(isfield(nii.hdr,'data_type'))
    nii.hdr.data_type	=bytematch(jnii.NIFTIHeader, 'A75DataTypeName', nii.hdr.data_type);
    nii.hdr.db_name	=bytematch(jnii.NIFTIHeader, 'A75DBName', nii.hdr.db_name);
    nii.hdr.extents	=bytematch(jnii.NIFTIHeader, 'A75Extends', nii.hdr.extents);
    nii.hdr.session_error	=bytematch(jnii.NIFTIHeader, 'A75SessionError', nii.hdr.session_error);
    nii.hdr.regular	=bytematch(jnii.NIFTIHeader, 'A75Regular', nii.hdr.regular);
end

dim_info=bitor(uint8(jnii.NIFTIHeader.DimInfo.Freq),bitshift(uint8(jnii.NIFTIHeader.DimInfo.Phase),3));
dim_info=bitor(dim_info,bitshift(uint8(jnii.NIFTIHeader.DimInfo.Slice),6));
nii.hdr.dim_info=cast(dim_info, class(nii.hdr.dim_info));

nii.hdr.dim(1)=cast(length(jnii.NIFTIHeader.Dim),class(nii.hdr.dim));
nii.hdr.dim(2:1+length(jnii.NIFTIHeader.Dim))	=bytematch(jnii.NIFTIHeader, 'Dim', nii.hdr.dim(2:1+length(jnii.NIFTIHeader.Dim)));
nii.hdr.intent_p1	=bytematch(jnii.NIFTIHeader, 'Param1', nii.hdr.intent_p1);
nii.hdr.intent_p2	=bytematch(jnii.NIFTIHeader, 'Param2', nii.hdr.intent_p2);
nii.hdr.intent_p3	=bytematch(jnii.NIFTIHeader, 'Param3', nii.hdr.intent_p3);

if(ischar(jnii.NIFTIHeader.Intent))
    jnii.NIFTIHeader.Intent=niicodemap('intent',jnii.NIFTIHeader.Intent);
end
nii.hdr.intent_code	=bytematch(jnii.NIFTIHeader, 'Intent', nii.hdr.intent_code);
if(ischar(jnii.NIFTIHeader.DataType))
    jnii.NIFTIHeader.DataType=niicodemap('datatype',jnii.NIFTIHeader.DataType);
end

nii.hdr.datatype	=bytematch(jnii.NIFTIHeader, 'DataType', nii.hdr.datatype);
nii.hdr.bitpix	=bytematch(jnii.NIFTIHeader, 'BitDepth', nii.hdr.bitpix);
nii.hdr.slice_start	=bytematch(jnii.NIFTIHeader, 'FirstSliceID', nii.hdr.slice_start);
nii.hdr.pixdim(1)=cast(length(jnii.NIFTIHeader.VoxelSize),class(nii.hdr.pixdim));
nii.hdr.pixdim(2:2+nii.hdr.dim(1)-1)	=bytematch(jnii.NIFTIHeader, 'VoxelSize', nii.hdr.pixdim(2:2+nii.hdr.dim(1)-1));
nii.hdr.vox_offset	=bytematch(jnii.NIFTIHeader, 'NIIByteOffset', nii.hdr.vox_offset);
nii.hdr.scl_slope	=bytematch(jnii.NIFTIHeader, 'ScaleSlope', nii.hdr.scl_slope);
nii.hdr.scl_inter	=bytematch(jnii.NIFTIHeader, 'ScaleOffset', nii.hdr.scl_inter);
nii.hdr.slice_end	=bytematch(jnii.NIFTIHeader, 'LastSliceID', nii.hdr.slice_end);

if(ischar(jnii.NIFTIHeader.SliceType))
    jnii.NIFTIHeader.SliceType=niicodemap('slicetype',jnii.NIFTIHeader.SliceType);
end
nii.hdr.slice_code	=bytematch(jnii.NIFTIHeader, 'SliceType', nii.hdr.slice_code);
if(ischar(jnii.NIFTIHeader.Unit.L))
    jnii.NIFTIHeader.Unit.L=niicodemap('unit',jnii.NIFTIHeader.Unit.L);
end
if(ischar(jnii.NIFTIHeader.Unit.T))
    jnii.NIFTIHeader.Unit.T=niicodemap('unit',jnii.NIFTIHeader.Unit.T);
end

xyzt_units=bitor(uint8(jnii.NIFTIHeader.Unit.L),uint8(jnii.NIFTIHeader.Unit.T));
nii.hdr.xyzt_units=cast(xyzt_units, class(nii.hdr.xyzt_units));

nii.hdr.cal_max	=bytematch(jnii.NIFTIHeader, 'MaxIntensity', nii.hdr.cal_max);
nii.hdr.cal_min	=bytematch(jnii.NIFTIHeader, 'MinIntensity', nii.hdr.cal_min);
nii.hdr.slice_duration	=bytematch(jnii.NIFTIHeader, 'SliceTime', nii.hdr.slice_duration);
nii.hdr.toffset	=bytematch(jnii.NIFTIHeader, 'TimeOffset', nii.hdr.toffset);
if(isfield(nii.hdr,'glmax'))
    nii.hdr.glmax	=bytematch(jnii.NIFTIHeader, 'A75GlobalMax', nii.hdr.glmax);
    nii.hdr.glmin	=bytematch(jnii.NIFTIHeader, 'A75GlobalMin', nii.hdr.glmin);
end

nii.hdr.descrip	=bytematch(jnii.NIFTIHeader, 'Description', nii.hdr.descrip);
nii.hdr.aux_file	=bytematch(jnii.NIFTIHeader, 'AuxFile', nii.hdr.aux_file);

nii.hdr.qform_code	=bytematch(jnii.NIFTIHeader, 'QForm', nii.hdr.qform_code);
nii.hdr.sform_code	=bytematch(jnii.NIFTIHeader, 'SForm', nii.hdr.sform_code);
nii.hdr.quatern_b	=bytematch(jnii.NIFTIHeader, 'Quatern.b', nii.hdr.quatern_b);
nii.hdr.quatern_c	=bytematch(jnii.NIFTIHeader, 'Quatern.c', nii.hdr.quatern_c);
nii.hdr.quatern_d	=bytematch(jnii.NIFTIHeader, 'Quatern.d', nii.hdr.quatern_d);
nii.hdr.qoffset_x	=bytematch(jnii.NIFTIHeader, 'QuaternOffset.x', nii.hdr.qoffset_x);
nii.hdr.qoffset_y	=bytematch(jnii.NIFTIHeader, 'QuaternOffset.y', nii.hdr.qoffset_y);
nii.hdr.qoffset_z	=bytematch(jnii.NIFTIHeader, 'QuaternOffset.z', nii.hdr.qoffset_z);
nii.hdr.srow_x	=cast(jnii.NIFTIHeader.Affine(1,:), class(nii.hdr.srow_x));
nii.hdr.srow_y	=cast(jnii.NIFTIHeader.Affine(2,:), class(nii.hdr.srow_y));
nii.hdr.srow_z	=cast(jnii.NIFTIHeader.Affine(3,:), class(nii.hdr.srow_z));

nii.hdr.intent_name	=bytematch(jnii.NIFTIHeader, 'Name', nii.hdr.intent_name);
%nii.hdr.magic	=bytematch(jnii.NIFTIHeader, 'NIIFormat', nii.hdr.magic);

if(isfield(jnii.NIFTIHeader,'NIIExtender'))
    nii.hdr.extension	=bytematch(jnii.NIFTIHeader, 'NIIExtender', nii.hdr.extension);
end
if(isfield(jnii.NIFTIHeader,'NIIQfac_'))
    nii.hdr.pixdim(1)	=bytematch(jnii.NIFTIHeader, 'NIIQfac_', nii.hdr.pixdim(1));
end
if(isfield(jnii.NIFTIHeader,'NIIUnused_'))
    nii.hdr.reserved	=bytematch(jnii.NIFTIHeader, 'NIIUnused_', nii.hdr.reserved);
end

if(isfield(jnii,'NIFTIExtension') && iscell(jnii.NIFTIExtension))
    nii.extension=jnii.NIFTIExtension;
    if(nii.hdr.extension(1)~=length(jnii.NIFTIExtension))
        nii.hdr.extension(1)=length(jnii.NIFTIExtension);
	warning('header extension count does not match the extension data, force update');
    end
end

if(nargin>=2 && ischar(varargin{1}))
    savenifti(nii.img, varargin{1}, nii.hdr);
end

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

function dat=bytematch(jobj, key, orig)
dtype=class(orig);
dat=orig;
if(isfield(jobj,key))
    dat=cast(jobj.(key),dtype);
else
    dat=cast(0,dtype);
end
if(length(dat)<length(orig))
    dat(length(orig))=cast(0,dtype);
end
dat=dat(1:length(orig));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% ACOTH  Inverse hyperbolic cotangent.
%    ACOTH(X) is the inverse hyperbolic cotangent of the elements of X.
%
%    Class support for input X:
%       float: double, single
%
%    See also COTH.
%
%    Reference page in Doc Center
%       doc acoth
%
%    Other functions named acoth
%
%       codistributed/acoth    gpuArray/acoth    sym/acoth    ts/acoth
%

author: jmaih


text:
function generate_rb_regressor(path_to_urdf)
% ---------------------------------------------------------------------
% The function generates regressor for UR10E. It assumes that links of
% the robot are rigid bodies. Thus the inverse dynamics can be written
% in linear-in-parameters from: tau = Y(q, dq, ddq) pi
% ---------------------------------------------------------------------

% Parse urdf to get robot description
ur10 = parse_urdf(path_to_urdf);

% Create symbolic generilized coordiates, their first and second deriatives
q_sym = sym('q%d',[6,1],'real');
qd_sym = sym('qd%d',[6,1],'real');
q2d_sym = sym('q2d%d',[6,1],'real');

% ------------------------------------------------------------------------
% Getting gradient of energy functions, to derive dynamics
% ------------------------------------------------------------------------
T_pk = sym(zeros(4,4,6)); % transformation between links
w_kk(:,1) = sym(zeros(3,1)); % angular velocity k in frame k
v_kk(:,1) = sym(zeros(3,1)); % linear velocity of the origin of frame k in frame k
g_kk(:,1) = sym([0,0,9.81])'; % vector of graviatational accelerations in frame k
p_kk(:,1) = sym(zeros(3,1)); % origin of frame k in frame k

for i = 1:6
    jnt_axs_k = str2num(ur10.robot.joint{i}.axis.Attributes.xyz)';
    % Transformation from parent link frame p to current joint frame
    rpy_k = sym(str2num(ur10.robot.joint{i}.origin.Attributes.rpy));
    R_pj = RPY(rpy_k);
    R_pj(abs(R_pj)<sqrt(eps)) = sym(0); % to avoid numerical errors
    p_pj = str2num(ur10.robot.joint{i}.origin.Attributes.xyz)';
    T_pj = sym([R_pj, p_pj; zeros(1,3), 1]); % to avoid numerical errors
    % Tranformation from joint frame of the joint that rotaties body k to
    % link frame. The transformation is pure rotation
    R_jk = Rot(q_sym(i),sym(jnt_axs_k));
    p_jk = sym(zeros(3,1));
    T_jk = [R_jk, p_jk; sym(zeros(1,3)),sym(1)];
    % Transformation from parent link frame p to current link frame k
    T_pk(:,:,i) = T_pj*T_jk;
    z_kk(:,i) = sym(jnt_axs_k);

    w_kk(:,i+1) = T_pk(1:3,1:3,i)'*w_kk(:,i) + sym(jnt_axs_k)*qd_sym(i);
    v_kk(:,i+1) = T_pk(1:3,1:3,i)'*(v_kk(:,i) + cross(w_kk(:,i),sym(p_pj)));
    g_kk(:,i+1) = T_pk(1:3,1:3,i)'*g_kk(:,i);
    p_kk(:,i+1) = T_pk(1:3,1:3,i)'*(p_kk(:,i) + sym(p_pj));

    beta_K(i,:) = [sym(0.5)*w2wtlda(w_kk(:,i+1)),...
                   v_kk(:,i+1)'*vec2skewSymMat(w_kk(:,i+1)),...
                   sym(0.5)*v_kk(:,i+1)'*v_kk(:,i+1)];
    beta_P(i,:) = [sym(zeros(1,6)), g_kk(:,i+1)',...
                   g_kk(:,i+1)'*p_kk(:,i+1)];
end

% ---------------------------------------------------------------------
% Dynamic regressor of the full paramters
% ---------------------------------------------------------------------
beta_Lf = [beta_K(1,:) - beta_P(1,:), beta_K(2,:) - beta_P(2,:),...
         beta_K(3,:) - beta_P(3,:), beta_K(4,:) - beta_P(4,:),...
         beta_K(5,:) - beta_P(5,:), beta_K(6,:) - beta_P(6,:)];
dbetaLf_dq = jacobian(beta_Lf,q_sym)';
dbetaLf_dqd = jacobian(beta_Lf,qd_sym)';
tf = sym(zeros(6,60));
for i = 1:6
   tf = tf + diff(dbetaLf_dqd,q_sym(i))*qd_sym(i)+...
                diff(dbetaLf_dqd,qd_sym(i))*q2d_sym(i);
end
Y_f = tf - dbetaLf_dq;

% Generate function from a symbolic expression for the regressor
matlabFunction(Y_f,'File','autogen/standard_regressor_UR10E',...
               'Vars',{q_sym,qd_sym,q2d_sym});

author: shamilmamedov


text:
function x = p29_start ( n )

%*****************************************************************************80
%
%% P29_START returns a starting point for optimization for problem 29.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = [ -0.5, +0.25 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef UniformD
%%UNIFORMD Functions to handle the scalar and multivariate (hyper-
%     rectangularly bounded) uniform distributions.
%Implemented methods are: mean, cov, PDF, CDF (scalar only), rand, entropy
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)
function val=mean(bounds)
%%MEAN Obtain the mean of the uniform distribution
%
%INPUTS: bounds A 2XnumDim matrix of the minimum and maximum bounds of each
%               of the numDim dimensions of the uniform distribution.
%               bounds(1,:) is the set of minimum bounds and bounds(2,:) is
%               the set of maximum bounds.
%
%OUTPUTS: val  The numDimX1 mean od the uniform distribution.
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=(bounds(2,:)+bounds(1,:))'/2;
end

function val=cov(bounds)
%%COV Obtain the covariance matrix of the uniform distribution (the
%     variance if scalar).
%
%INPUTS: bounds A 2XnumDim matrix of the minimum and maximum bounds of each
%               of the numDim dimensions of the uniform distribution.
%               bounds(1,:) is the set of minimum bounds and bounds(2,:) is
%               the set of maximum bounds.
%
%OUTPUTS: val  The covariance matrix of the uniform distribution. This is a
%              diagonal matrix.
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=diag((1/12)*(bounds(2,:)-bounds(1,:)).^2);
end

function vals=PDF(z,bounds)
%%PDF Evaluate a scalar or multivariate uniform  PDF at a certain point
%     given the bounds of the PDF.
%
%INPUTS: z The points at which the PDF should be evaluated. If the PDF is
%          multivariate, then this is a column vector. If evaluation at
%          multiple points are desired, then this is a matrix with each
%          column being the a point (a vector).
%   bounds A 2XnumDim matrix of the minimum and maximum bounds of each of
%          the numDim dimensions of the uniform distribution. bounds(1,:)
%          is the set of minimum bounds and bounds(2,:) is the set of
%          maximum bounds.
%
%OUTPUTS: vals The scalar value of the uniform PDF with the given bounds.
%              If multiple points are passed (z is a matrix), then val is a
%              column vector.
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    numDim=size(z,1);
    numPoints=size(z,2);
    val=1/prod(bounds(2,:)-bounds(1,:));
    vals=repmat(val,[numPoints,1]);

    sel=false(1,numPoints);
    for curDim=1:numDim
        sel=sel|z(curDim,:)<bounds(1,curDim)|z(curDim,:)>bounds(2,curDim);
    end
    vals(sel)=0;
end

function val=CDF(z,bounds)
%%CDF Evaluate a scalar uniform cumulative distribution function (CDF) at a
%     certain point given the bounds of the distribution.
%
%INPUTS: z The scalar point(s) at which the CDF should be evaluated.
%   bounds A 2X1 vector of the minimum and maximum bounds of the scalar
%          uniform distribution. bounds(1) is the minimum bound and
%          bounds(2) is the maximum bound.
%
%OUTPUTS: val The scalar value(s) of the uniform CDF with the given
%             bounds at point(s) in z.
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=(z-bounds(1))/(bounds(2)-bounds(1));
    val(z<=bounds(1))=0;
    val(z>=bounds(2))=1;
end

function x=rand(N,bounds)
%%RAND Generate multivariate uniform random variables with a given set of
%      bounds
%
%INPUTS: N The scalar number of random variables to generate.
%  bounds A 2XnumDim matrix of the minimum and maximum integer bounds of
%         each of the numDim dimensions of the uniform distribution.
%         bounds(1,:) is the set of minimum bounds and bounds(2,:) is
%         the set of maximum bounds.
%
%OUTPUT: x A sample of an numDimXN set of random uniform variables with the
%          given bounds.
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    xDim=size(bounds,2);

    x=(bounds(2,:)-bounds(1,:))'.*rand(xDim,N)+bounds(1,:)';
end

function entropyVal=entropy(bounds)
%%ENTROPY Obtain the differential entropy of the multivariate uniform
%         distribution given in nats. The differential entropy of a
%         continuous distribution is entropy=-int_x p(x)*log(p(x)) dx where
%         the integral is over all values of x. Units of nats mean that the
%         natural logarithm is used in the definition. Unlike the Shannon
%         entropy for discrete variables, the differential entropy of
%         continuous variables can be both positive and negative.
%
%INPUTS: bounds A 2XnumDim matrix of the minimum and maximum integer bounds
%               of each of the numDim dimensions of the uniform
%               distribution. bounds(1,:) is the set of minimum bounds and
%               bounds(2,:) is the set of maximum bounds.
%
%OUTPUTS: entropyVal The value of the differential entropy in nats.
%
%Differential entropy is defined in Chapter 8 of [1].
%
%REFERENCES:
%[1] T. M. Cover and J. A. Thomas, Elements of Information Theory, 2nd ed.
%    Hoboken, NJ: Wiley-Interscience, 2006.
%
%April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    entropyVal=log(prod(bounds(2,:)-bounds(1,:)));
end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
function output = callsdpnal(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

% Bounded variables converted to constraints
if ~isempty(ub)
    [F_struc,K] = addStructureBounds(F_struc,K,ub,lb);
end

[blk,A,C,b,oldKs]=sedumi2sdpt3(F_struc(:,1),F_struc(:,2:end),c,K,1);

if options.savedebug
    ops = options.sdpnal;
    save sdpnaldebug blk A C b ops -v6
end

if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
solvertime = tic;
if options.verbose==0
    evalc('[obj,X,s,y,Z,Z2,y2,v,info,runhist] = sdpnalplus(blk,A,C,b,[],[],[],[],[],options.sdpnal);');
else
    [obj,X,s,y,Z,Z2,y2,v,info,runhist] = sdpnalplus(blk,A,C,b,[],[],[],[],[],options.sdpnal);
end
solvertime = toc(solvertime);

% Create YALMIP dual variable and slack
Dual = [];
Slack = [];
top = 1;
if K.f>0
    Dual = [Dual;X{top}(:)];
    Slack = [Slack;Z{top}(:)];
    top = top+1;
end
if K.l>0
    Dual = [Dual;X{top}(:)];
    Slack = [Slack;Z{top}(:)];
    top = top + 1;
end
if any(K.q)
    Dual = [Dual;X{top}(:)];
    Slack = [Slack;Z{top}(:)];
    top = top + 1;
end
if any(K.s)
    % Messy format in SDPT3 to block and sort small SDPs
    u = blk(:,1);
    u = find([u{:}]=='s');
    s = 1;
    for top = u
        ns = blk(top,2);ns = ns{1};
        k = 1;
        for i = 1:length(ns)
            Xi{oldKs(s)} = X{top}(k:k+ns(i)-1,k:k+ns(i)-1);
            Zi{oldKs(s)} = Z{top}(k:k+ns(i)-1,k:k+ns(i)-1);
            s = s + 1;
            k = k+ns(i);
        end
    end
    for i = 1:length(Xi)
        Dual = [Dual;Xi{i}(:)];
        Slack = [Slack;Zi{i}(:)];
    end
end
Primal = -y;  % Primal variable in YALMIP

% No error code available
if isfield(info,'termcode')
    switch info.termcode
        case -1
            problem = 4;
        case -2
            problem = 3;
        case 0
            problem = 0;
        case 1
            problem = 5;
        case 2
            problem = 3;
        otherwise
            problem = 11;
    end
else
    if isfield(info,'msg')
        if isequal(info.msg,'maximum iteration reached')
            problem = 3;
        else
            problem = 9;
        end
    end
end

if options.savesolveroutput
    solveroutput.obj = obj;
    solveroutput.X = X;
    solveroutput.y = y;
    solveroutput.Z = Z;
    solveroutput.info = info;
    solveroutput.runhist = runhist;
 else
    solveroutput = [];
end

if options.savesolverinput
    solverinput.blk = blk;
    solverinput.A   = A;
    solverinput.C   = C;
    solverinput.b   = b;
    solverinput.options   = options.sdpnal;
else
    solverinput = [];
end

% Standard interface
output = createOutputStructure(Primal,Dual,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

function [F_struc,K] = deblock(F_struc,K);
X = any(F_struc(end-K.s(end)^2+1:end,:),2);
X = reshape(X,K.s(end),K.s(end));
[v,dummy,r,dummy2]=dmperm(X);
blks = diff(r);

lint = F_struc(1:end-K.s(end)^2,:);
logt = F_struc(end-K.s(end)^2+1:end,:);

newlogt = [];
for i = 1:size(logt,2)
    temp = reshape(logt(:,i),K.s(end),K.s(end));
    temp = temp(v,v);
    newlogt = [newlogt temp(:)];
end
logt = newlogt;

pattern = [];
for i = 1:length(blks)
    pattern = blkdiag(pattern,ones(blks(i)));
end

F_struc = [lint;logt(find(pattern),:)];
K.s(end) = [];
K.s = [K.s blks];
K.m = blks;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function result = computeMutualCoherence(A)
% Computes mutual coherence
%
% USAGE:
%
%    result = computeMutualCoherence(A)
%
% INPUT:
%    A:         input matrix
%
% OUTPUT:
%    result:    mutual coherence
%
% ..
%     [M N] = size(A);
%     if (N<2)
%         disp('error - input contains only one column');
%         u=NaN;   beep;    return
%     end
%
%     % normalize the columns
%     nn = sqrt(sum(A.*conj(A),1));
%     if ~all(nn)
%         disp('error - input contains a zero column');
%         u=NaN;   beep;    return
%     end
%     nA = bsxfun(@rdivide,A,nn);  % nA is a matrix with normalized columns
%     result = max(max(triu(abs((nA')*nA),1)));

    X = A/(diag(sqrt(diag(A'*A))));
    mcoh = abs(X'*X);
    V = diag(mcoh);
    V = -diag(V,0);
    mcoh = mcoh+V;
    result = max(max(mcoh));
end

author: opencobra

repo: cobratoolbox


text:
function c = cat(varargin)

% c = cat(c1,c2,c3,...)
% Concatenate the traces in two or more correlation objects. c1,c2,c3,...
% will be combined in the order listed. Only the four data descriptor
% fields (TRIG,START,WAVES,FS) will be filled in the output object. There
% is not enough information to fill in the CORRELATION, LAG, STAT, LINK, or
% CLUSTER fields. These fields will be empty. Data in c1,c2,c3,... must
% have the same number of samples and sampling frequency.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

if nargin <= 1
    error('Not enough inputs');
end

for i = 1:nargin
    if ~isa(varargin{i},'correlation')
        error('All inputs must be correlation objects');
    end
    if get(varargin{i},'Fs') ~= get(varargin{1},'Fs')
        error('All inputs must have the same sampling frequency');
    end
    if size(get(varargin{i},'DATA_LENGTH'),1) ~= size(get(varargin{1},'DATA_LENGTH'),1)
        error('All inputs must have the same trace length');
    end
end

c = varargin{1};
c.C = [];
c.L = [];
c.stat = [];
c.link = [];
c.clust = [];

for i = 2:nargin
    cin = varargin{i};
    c.W = cat(1,c.W,cin.W);
    c.trig = cat(1,c.trig,cin.trig);
end

author: geoscience-community-codes

repo: GISMO


text:
function varargout = coeffs2( f, m, n )
%COEFFS2   Fourier--Chebyshev coefficients of a DISKFUN.
%
%   X = COEFFS2( F ) returns the coefficients of the DISKFUN F in the
%   Fourier--Chebyshev bases, where the columns are the Chebyshev
%   coefficients and the rows are Fourier coefficients.
%
%   [C, D, R] = COEFFS2( F ) returns a low rank approximation to the
%   coefficients.
%
%   X = COEFFS2(F, M, N) returns bivariate coefficients with N Chebyshev
%   coefficients in the radial direction and M Fourier coefficients in
%   angular direction.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

%empty check
if isempty(f)
    varargout = {};
    return
end

% Calculate the CDR decomposition:
[C, D, R] = cdr( f );

if nargin == 1
    % Find the coefficients of each slice:
    U = C.coeffs;
    R = R.coeffs;
else
    if nargin == 2
        n = m;
    end
    % Find the coefficients of each slice:
    U = chebtech2.alias(C.coeffs,n);
    R = trigtech.alias(R.coeffs,m);
end

% Prepare the output. Keep in low rank form if nargin > 1.
if ( nargout <= 1 )
    varargout = {U*D*R.'};
elseif ( nargout <= 3 )
    varargout = {U, D, R};
else
    error('DISKFUN:COEFFS:NARGOUT',...
            'Too many output arguments')
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function b = r8mat_power ( n, a, npow )

%*****************************************************************************80
%
%% R8MAT_POWER computes a nonnegative power of an R8MAT.
%
%  Discussion:
%
%    The algorithm is:
%
%      B = I
%      do NPOW times:
%        B = A * B
%      end
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real A(N,N), the matrix to be raised to a power.
%
%    Input, integer NPOW, the power to which A is to be raised.
%    N must be nonnegative.
%
%    Output, real B(N,N), the value of A**NPOW.
%
  if ( npow < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_POWER - Fatal error!\n' );
    fprintf ( 1, '  Input value of NPOW < 0.\n' );
    error ( 'R8MAT_POWER - Fatal error!' );
  end

  b = r8mat_identity ( n );

  for ipow = 1 : npow
    b(1:n,1:n) = a(1:n,1:n) * b(1:n,1:n);
  end

  return
end


text:
%BLOCKLOGDISTRIBUTION
%
%   LogBS = BlockLogDistribution(V, G, F, A) returns the log of a
%   block-sampling array (which contains the log-unnormalized-probabilities of
%   selecting each label for the block), given variables V to block-sample in
%   network G with factors F and current assignment A.  Note that the variables
%   in V must all have the same dimensionality.
%
%   Input arguments:
%   V -- an array of variable names
%   G -- the graph with the following fields:
%     .names - a cell array where names{i} = name of variable i in the graph
%     .card - an array where card(i) is the cardinality of variable i
%     .edges - a matrix such that edges(i,j) shows if variables i and j
%              have an edge between them (1 if so, 0 otherwise)
%     .var2factors - a cell array where var2factors{i} gives an array where the
%              entries are the indices of the factors including variable i
%   F -- a struct array of factors.  A factor has the following fields:
%       F(i).var - names of the variables in factor i
%       F(i).card - cardinalities of the variables in factor i
%       F(i).val - a vectorized version of the CPD for factor i (raw probability)
%   A -- an array with 1 entry for each variable in G s.t. A(i) is the current
%       assignment to variable i in G.
%
%   Each entry in LogBS is the log-probability that that value is selected.
%   LogBS is the P(V | X_{-v} = A_{-v}, all X_i in V have the same value), where
%   X_{-v} is the set of variables not in V and A_{-v} is the corresponding
%   assignment to these variables consistent with A.  In the case that |V| = 1,
%   this reduces to Gibbs Sampling.  NOTE that exp(LogBS) is not normalized to
%   sum to one at the end of this function (nor do you need to worry about that
%   in this function).
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function LogBS = BlockLogDistribution(V, G, F, A)
if length(unique(G.card(V))) ~= 1
    disp('WARNING: trying to block sample invalid variable set');
    return;
end

% d is the dimensionality of all the variables we are extracting
d = G.card(V(1));

LogBS = zeros(1, d);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE
% Compute LogBS by multiplying (adding in log-space) in the correct values from
% each factor that includes some variable in V.
%
% NOTE: As this is called in the innermost loop of both Gibbs and Metropolis-
% Hastings, you should make this fast.  You may want to make use of
% G.var2factors, repmat,unique, and GetValueOfAssignment.
%
% Also you should have only ONE for-loop, as for-loops are VERY slow in matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

factors = G.var2factors{V(1)};
for i = 2:length(V)
	factors =unique([factors,G.var2factors{V(i)}]);
end

fp = F(factors(1));
for i = 2:length(factors)
	fp = FactorProduct(fp,F(factors(i)));
end

x = repmat(A,d,1);
x(:,V) = repmat([1:d]',1,length(V));
for i = 1:d
	LogBS(i) = log(GetValueOfAssignment(fp,x(i,:),1:length(A)));
end

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

% Re-normalize to prevent underflow when you move back to probability space
LogBS = LogBS - min(LogBS);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
% test_polygeom.m - test polygeom
%   area, centroid, perimeter and area moments of polygonal outline
% H.J. Sommer III - 02.05.14 - tested under MATLAB v5.2

clear

% constants
d2r = pi / 180;

% 3x5 test rectangle with long axis at 30 degrees
% area=15, x_cen=3.415, y_cen=6.549, perimeter=16
% I1=11.249, I2=31.247, J=42.496
x = [ 2.000  0.500  4.830  6.330 ]';
y = [ 4.000  6.598  9.098  6.500 ]';

% get geometry
[ geom, iner, cpmo ] = polygeom( x, y );

% show results
area = geom(1);
x_cen = geom(2);
y_cen = geom(3);
perimeter = geom(4);
disp( [ ' ' ] )
disp( [ '3x5 test rectangle with long axis at 30 degrees' ] )
disp( [ ' ' ] )
disp( [ '      area     x_cen     y_cen     perim' ] )
disp( [ area  x_cen  y_cen  perimeter ] )

I1 = cpmo(1);
angle1 = cpmo(2);
I2 = cpmo(3);
angle2 = cpmo(4);
disp( [ ' ' ] )
disp( [ '        I1        I2' ] )
disp( [ I1 I2 ] )
disp( [ '    angle1    angle2' ] )
disp( [ angle1/d2r angle2/d2r ] )

% plot outline
xplot = [ x ; x(1) ];
yplot = [ y ; y(1) ];
rad = 10;
x1 = [ x_cen-rad*cos(angle1)  x_cen+rad*cos(angle1) ];
y1 = [ y_cen-rad*sin(angle1)  y_cen+rad*sin(angle1) ];
x2 = [ x_cen-rad*cos(angle2)  x_cen+rad*cos(angle2) ];
y2 = [ y_cen-rad*sin(angle2)  y_cen+rad*sin(angle2) ];
plot( xplot,yplot,'b', x_cen,y_cen,'ro', ...
      x1,y1,'g:', x2,y2,'g:'  )
axis( [ 0  rad  0  rad ] )
axis square

author: Sable

repo: mcbench-benchmarks


text:
function [Xf]=nspd(A)

% function [Xf]=nspd(A)
% checks whether a matrix A is symmetric positive definite (SPD), and if not, computes the nearest SPD matrix to A.
% typically, A is either SPD, or close to SPD but not recognised as such because of a very small numerical perturbation
% this functions helps to fix the problem by producing a matrix Xf that is for sure SPD
% it uses the methodology proposed by N.J Highams (1988): 'Computing a nearest symmetric positive semidefinite matrix',Linear Algebra and its Applications, No. 103
% inputs:  - matrix A: must be a square matrix
% outputs: - matrix 'Xf': the closest symmetric positive definite matrix to A

% first, check if the matrix is SPD by attempting to take the choleski factor
[~,test1]=chol(A,'lower');

% if it is SPD in the first place (test value=0), no need to apply any change
if test1<1
Xf=A;

% if it is not SPD, fix it
elseif test1>=1

% obtain the matrix Xf by applying theorem 2.1

% create B, as defined in theorem 2.1
B=(A+A')/2;
% produce the polar decomposition of B
% first, implement a singular value decomposition of B
[U,S,V]=svd(B);
% then recover the polar matrix H from the singular value decomposition formula: H=V*S*V'
H=V*S*V';
% finally, recover Xf from the formula in theorem 2.1: Xf=(B+H)/2
Xf=(B+H)/2;
% ensure symmetry

Xf=(Xf+Xf')/2;

% now the matrix should be symmetric positive definite: test again
[~,test2]=chol(Xf,'lower');

   % if the test is passed and Choleski factorisation was possible (test2=0), then Xf needs no further work;
   if test2<1

   % if the test is not passed, it must be because of some very small numerical disturbance
   % hence keep modifying Xf very slightly until it becomes SPD
   elseif test2>=1
      n=1;
      while test2>=1
      mineig=min(eig(Xf));
      Xf=Xf+(-mineig*n.^2+eps(mineig)+0.0000000000000000001)*eye(size(Xf));
      % test again
      [~,test2]=chol(Xf,'lower');
      n=n+1;
      end

   end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function varargout = process_source_flat( varargin )
% PROCESS_SOURCE_FLAT: Convert an unconstrained source file into a flat map.
%
% USAGE:  OutputFiles = process_source_flat('Run', sProcess, sInput)
%          ResultsMat = process_source_flat('Compute', ResultsMat, Method, Field)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2013-2015

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % ===== PROCESS =====
    % Description the process
    sProcess.Comment     = 'Unconstrained to flat map';
    sProcess.Category    = 'File';
    sProcess.SubGroup    = 'Sources';
    sProcess.Index       = 337;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation?highlight=(Unconstrained+to+flat+map)#Z-score';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'results'};
    sProcess.OutputTypes = {'results'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % === SELECT METHOD
    sProcess.options.label1.Comment = ['Converts unconstrained source files (3 values per vertex: x,y,z)
' 10 ...
                                       'to simpler files with only one value per vertex.

' 10 ...
                                       'Method used to perform this conversion:'];
    sProcess.options.label1.Type    = 'label';
    sProcess.options.method.Comment = {'Norm: sqrt(x^2+y^2+z^2)', ...
                                       'PCA: First mode of svd(x,y,z)'};
    sProcess.options.method.Type    = 'radio';
    sProcess.options.method.Value   = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInput) %#ok<DEFNU>
    OutputFiles = {};
    % Get options
    switch(sProcess.options.method.Value)
        case 1, Method = 'rms';  fileTag = 'norm';
        case 2, Method = 'pca';  fileTag = 'pca';
    end

    % ===== PROCESS INPUT =====
    % Load the source file
    ResultsMat = in_bst_results(sInput.FileName, 1);
    % Error: Not an unconstrained model
    if (ResultsMat.nComponents == 1) || (ResultsMat.nComponents == 2)
        bst_report('Error', sProcess, sInput, 'The input file is not an unconstrained source model.');
        return;
    end
    % Compute flat map
    ResultsMat = Compute(ResultsMat, Method);

    % ===== SAVE FILE =====
    % File tag
    if ~isempty(strfind(sInput.FileName, '_abs_zscore'))
        FileType = 'results_abs_zscore';
    elseif ~isempty(strfind(sInput.FileName, '_zscore'))
        FileType = 'results_zscore';
    else
        FileType = 'results';
    end
    % Get study description
    sStudy = bst_get('Study', sInput.iStudy);
    % Output filename
    OutputFiles{1} = bst_process('GetNewFilename', bst_fileparts(sStudy.FileName), [FileType '_' fileTag]);
    % Save on disk
    bst_save(OutputFiles{1}, ResultsMat, 'v6');
    % Register in database
    db_add_data(sInput.iStudy, OutputFiles{1}, ResultsMat);
end

%% ===== COMPUTE =====
function ResultsMat = Compute(ResultsMat, Method, Field)
    % Field to process
    if (nargin < 3) || isempty(Field)
        Field = 'ImageGridAmp';
    end
    if (nargin < 2) || isempty(Method)
        Method = 'rms';
    end
    % Process differently the different types of source spaces
    switch (ResultsMat.nComponents)
        % Mixed source models
        case 0
            % Apply orientations for each region independently
            [ResultsMat.(Field), ResultsMat.GridAtlas] = bst_source_orient([], ResultsMat.nComponents, ResultsMat.GridAtlas, ResultsMat.(Field), Method);        % Constrained source models
        case 1
            % The input file is not an unconstrained source model, nothing to do
        case 2
            error('Not supported.');
        % Unconstrained source models
        case 3
            % Apply source orientations
            ResultsMat.(Field) = bst_source_orient([], ResultsMat.nComponents, [], ResultsMat.(Field), Method);
            % Set the number components
            ResultsMat.nComponents = 1;
    end
    % Save new file function
    ResultsMat.Function = Method;

    % File tag
    switch(Method)
        case 'rms',  fileTag = 'norm';
        case 'pca',  fileTag = 'pca';
    end
    % Reset the data file initial path
    if ~isempty(ResultsMat.DataFile)
        ResultsMat.DataFile = file_win2unix(file_short(ResultsMat.DataFile));
    end
    % Add comment
    ResultsMat.Comment = [ResultsMat.Comment ' | ' fileTag];
    % Add history entry
    ResultsMat = bst_history('add', ResultsMat, 'flat', ['Convert unconstrained sources to a flat map with option: ' fileTag]);
end

author: brainstorm-tools


text:
function F = diag(f,varargin)
%DIAG(F)   Diagonal of a DISKFUN.
%   G = DIAG(F) returns the CHEBFUN representing g(r) = f(pi/4, r), where r
%   is the radial variable and is taken across the whole disk: -1 <= r <= 1
%
%   G = diag(F,T) returns the CHEBFUN representing g(r) = f(T, r), where
%  -pi < T < pi and r is the radial variable.
%
% See also DISKFUN/TRACE.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( isempty( f ) )
    F = chebfun;
    return
end

% When c = 0, choose the diagonal radial slice for t = pi/4
if (nargin<2)
    F = feval(f, pi/4,':', 'polar');
else
    c = varargin{1};
    if abs(c) >= pi
        error('CHEBFUN:DISKFUN:diag:angleOutOfRange','Diagonal parameter must be between -pi and pi.');
    else
        F = feval(f, varargin{1},':', 'polar');
    end
end

end


text:
function [output,timing] = global_solve_upper(p,p_original,x,options,uppersolver,timing)

if ~isempty(p.binary_variables)
    local_gave_good = find(abs(x(p_original.binary_variables)-fix(x(p_original.binary_variables)))< options.bnb.inttol);
    p.lb(p.binary_variables(local_gave_good)) = fix(x(p.binary_variables(local_gave_good)));
    p.ub(p.binary_variables(local_gave_good)) = fix(x(p.binary_variables(local_gave_good)));
    for i = 1:3
        p = update_eval_bounds(p);
        p = propagate_bounds_from_equalities(p);
        p = update_monomial_bounds(p);
    end
end

if ~isempty(p_original.integer_variables)
    local_gave_good = find(abs(x(p.integer_variables)-fix(x(p.integer_variables)))< options.bnb.inttol);
    p.lb((p.integer_variables(local_gave_good))) = fix(x(p.integer_variables(local_gave_good)));
    p.ub((p.integer_variables(local_gave_good))) = fix(x(p.integer_variables(local_gave_good)));
    for i = 1:3
        p = update_eval_bounds(p);
        p = propagate_bounds_from_equalities(p);
        p = update_monomial_bounds(p);
    end
end

% The bounds and relaxed solutions have been computed w.r.t to the relaxed
% bilinear model. We only need the original bounds and variables.
p.lb = p.lb(1:length(p_original.c));
p.ub = p.ub(1:length(p_original.c));

x = x(1:length(p_original.c));
%
p_upper = p_original;

p_upper.lb = p.lb;
p_upper.ub = p.ub;

% KNITRO does not like fixed binary/integer variables, so we remove them
fixed_binary = find(p_upper.lb(p_upper.binary_variables) == p_upper.ub(p_upper.binary_variables));
fixed_integer = find(p_upper.lb(p_upper.integer_variables) == p_upper.ub(p_upper.integer_variables));
p_upper.binary_variables(fixed_binary) = [];if length(p_upper.binary_variables)==0;p_upper.binary_variables = [];end
p_upper.integer_variables(fixed_integer) = [];if length(p_upper.integer_variables)==0;p_upper.integer_variables = [];end

% Pick an initial point (this can be a bit tricky...)
lbinfbounds = find(isinf(p.lb));
ubinfbounds = find(isinf(p.ub));
switch p.options.bmibnb.localstart
    case 'relaxed'
        p_upper.x0 = x;
    case 'boxcenter'
        p_upper.x0 = x;
        p_upper.x0(lbinfbounds)=0;
        p_upper.x0(ubinfbounds)=0;
    case 'zero'
        p_upper.x0 = zeros(length(p.lb),1);
    otherwise
        error('Unknown local starting point option');
end

% Find x(a)*x(b) where x(a) or x(b) is fixed, and remove from Q
[a,b,kk] = find(triu(2*p_upper.Q));
for i = 1:length(a)
    if p.lb(a(i)) == p.ub(a(i))
        p_upper.c(b(i)) = p_upper.c(b(i)) + p.lb(a(i))*kk(i);
        p_upper.Q(a(i),b(i))=0;
        p_upper.Q(b(i),a(i))=0;
    elseif p.lb(b(i)) == p.ub(b(i))
        p_upper.c(a(i)) = p_upper.c(a(i)) + p.lb(b(i))*kk(i);
        p_upper.Q(a(i),b(i))=0;
        p_upper.Q(b(i),a(i))=0;
    end
end

% Save time
p_upper.options.saveduals = 0;

% Solve upper bounding problem
p_upper.options.usex0 = 1;
tstart = tic;
try
    output = feval(uppersolver,p_upper);
catch
    output.Primal = zeros(length(p_upper.lb),1);
    output.Problem = -1;
end
if isempty(output.Primal)
     output.Primal = zeros(length(p_upper.lb),1);
end

timing.uppersolve = timing.uppersolve + toc(tstart);

% Project into the box (numerical issue)
output.Primal(output.Primal<p_upper.lb) = p_upper.lb(output.Primal<p_upper.lb);
output.Primal(output.Primal>p_upper.ub) = p_upper.ub(output.Primal>p_upper.ub);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function condition_test02 ( )

%*****************************************************************************80
%
%% CONDITION_TEST02 tests CONDITION_SAMPLE1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 April 2012
%
%  Author:
%
%    John Burkardt
%
  m_test = [ 10, 1000, 100000 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CONDITION_TEST02\n' );
  fprintf ( 1, '  For a matrix in general storage,\n' );
  fprintf ( 1, '  CONDITION_SAMPLE1 estimates the L1 condition number using sampling.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix                 Samples Order   Condition        Estimate\n' );
%
%  Combinatorial matrix.
%
  name = 'Combinatorial';
  n = 4;
  alpha = 2.0;
  beta = 3.0;
  a = combin ( alpha, beta, n );
  a_inverse = combin_inverse ( alpha, beta, n );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1  = a_norm_l1 * a_inverse_norm_l1;
  fprintf ( 1, '\n' );
  for i = 1 : 3
    m = m_test(i);
    cond = condition_sample1 ( n, a, m );
    fprintf ( 1, '  %20s  %8d  %4d  %14.6g  %14.6g\n', name, m, n, cond_l1, cond );
  end
%
%  CONEX1
%
  name = 'CONEX1';
  n = 4;
  alpha = 100.0;
  a = conex1 ( alpha );
  a_inverse = conex1_inverse ( alpha );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1  = a_norm_l1 * a_inverse_norm_l1;
  fprintf ( 1, '\n' );
  for i = 1 : 3
    m = m_test(i);
    cond = condition_sample1 ( n, a, m );
    fprintf ( 1, '  %20s  %8d  %4d  %14.6g  %14.6g\n', name, m, n, cond_l1, cond );
  end
%
%  CONEX2
%
  name = 'CONEX2';
  n = 3;
  alpha = 100.0;
  a = conex2 ( alpha );
  a_inverse = conex2_inverse ( alpha );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1  = a_norm_l1 * a_inverse_norm_l1;
  fprintf ( 1, '\n' );
  for i = 1 : 3
    m = m_test(i);
    cond = condition_sample1 ( n, a, m );
    fprintf ( 1, '  %20s  %8d  %4d  %14.6g  %14.6g\n', name, m, n, cond_l1, cond );
  end
%
%  CONEX3
%
  name = 'CONEX3';
  n = 5;
  a = conex3 ( n );
  a_inverse = conex3_inverse ( n );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1  = a_norm_l1 * a_inverse_norm_l1;
  fprintf ( 1, '\n' );
  for i = 1 : 3
    m = m_test(i);
    cond = condition_sample1 ( n, a, m );
    fprintf ( 1, '  %20s  %8d  %4d  %14.6g  %14.6g\n', name, m, n, cond_l1, cond );
  end
%
%  CONEX4
%
  name = 'CONEX4';
  n = 4;
  a = conex4 ( );
  a_inverse = conex4_inverse ( );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1  = a_norm_l1 * a_inverse_norm_l1;
  fprintf ( 1, '\n' );
  for i = 1 : 3
    m = m_test(i);
    cond = condition_sample1 ( n, a, m );
    fprintf ( 1, '  %20s  %8d  %4d  %14.6g  %14.6g\n', name, m, n, cond_l1, cond );
  end
%
%  KAHAN
%
  name = 'KAHAN';
  n = 4;
  alpha = 0.25;
  a = kahan ( alpha, n, n );
  a_inverse = kahan_inverse ( alpha, n );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1  = a_norm_l1 * a_inverse_norm_l1;
  fprintf ( 1, '\n' );
  for i = 1 : 3
    m = m_test(i);
    cond = condition_sample1 ( n, a, m );
    fprintf ( 1, '  %20s  %8d  %4d  %14.6g  %14.6g\n', name, m, n, cond_l1, cond );
  end
%
%  Random
%
  seed = 123456789;

  for j = 1 : 5
    name = 'RANDOM';
    n = 4;
    [ a, seed ] = r8mat_uniform_01 ( n, n, seed );
    a_inverse = inv ( a );
    a_norm_l1 = norm ( a, 1 );
    a_inverse_norm_l1 = norm ( a_inverse, 1 );
    cond_l1  = a_norm_l1 * a_inverse_norm_l1;
    fprintf ( 1, '\n' );
    for i = 1 : 3
      m = m_test(i);
      cond = condition_sample1 ( n, a, m );
      fprintf ( 1, '  %20s  %8d  %4d  %14.6g  %14.6g\n', name, m, n, cond_l1, cond );
    end
  end

  return
end


text:
function s = char(sR)
% convert spherical region to char

switch length(sR.N)
  case 0
    s = 'full sphere';
  case 1
    if sR.isUpper
      s = 'upper hemisphere';
    else
      s = 'lower hemisphere';
    end
  otherwise
    s = 'sector';
end

end

author: mtex-toolbox

repo: mtex

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


text:
function [pvec, pstruct] = tapas_logrt_linear_whatworld_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2014 Christoph Mathys, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pvec    = NaN(1,length(ptrans));
pstruct = struct;

pvec(1)     = ptrans(1);         % be0
pstruct.be0 = pvec(1);
pvec(2)     = ptrans(2);         % be1
pstruct.be1 = pvec(2);
pvec(3)     = ptrans(3);         % be2
pstruct.be2 = pvec(3);
pvec(4)     = ptrans(4);         % be3
pstruct.be3 = pvec(4);
pvec(5)     = exp(ptrans(5));    % ze
pstruct.ze  = pvec(5);

return;

author: translationalneuromodeling

repo: tapas


text:
classdef PTKViewer < GemFigure
    % PTKViewer. A standalone image viewer for showing 3D images slice-by-slice
    %
    %     PTKViewer uses MimViewerPanel to create a visualisation window for the
    %     supplied image, which can be a PTKImage or raw data. If a raw data
    %     matrix is supplied, the type argument can be supplied to ensure the
    %     image is displayed as expected.
    %
    %     Syntax:
    %         obj = PTKViewer;
    %         obj = PTKViewer(image);
    %         obj = PTKViewer(image, image_type);
    %         obj = PTKViewer(image, image_type, reporting);
    %         obj = PTKViewer(image, overlay_image);
    %         obj = PTKViewer(image, overlay_image, reporting);
    %
    %     Examples
    %     --------
    %
    %         % Displays my_image, where my_image is a 3D image matrix or a PTKImage
    %         PTKViewer(my_image);
    %
    %         % Displays my_image as a greyscale image, where my_image is a 3D image.
    %         % If my_Image is a PTKImage, the type will be ignored
    %         PTKViewer(my_image, PTKImageType.Grayscale);
    %
    %         % Displays my_image in the background with overlay_image as an overlay
    %         viewer = PTKViewer(my_image, overlay_image);
    %
    %         % Displays my_image as a greyscale image, and overlay_image as an indexed colour image
    %         viewer = PTKViewer(my_image, PTKImageType.Grayscale);
    %         viewer.ViewerPanelHandle.OverlayImage = overlay_image; % Must be a PTKImage
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ViewerPanelHandle
    end

    properties (SetAccess = private)
        Image
        Overlay
    end

    methods
        function obj = PTKViewer(image, overlay_or_type, reporting)

            if nargin < 1
                image = [];
            end

            % The second argument can be the type of the image, or an overlay image
            if nargin < 2
                overlay = [];
                type = [];
            else
                if isa(overlay_or_type, 'PTKImageType')
                    overlay = [];
                    type = overlay_or_type;
                elseif isa(overlay_or_type, 'PTKImage')
                    overlay = overlay_or_type;
                    type = [];
                end
            end

            if nargin < 3
                reporting = CoreReportingDefault;
            end

            title = 'PTK Viewer';

            if isa(image, 'PTKImage')
                image_handle = image;
                if ~isempty(image.Title)
                    title = image.Title;
                end
            else
                if isempty(type)
                    image_handle = PTKImage(image);
                else
                    image_handle = PTKImage(image, type);
                end
            end

            if isa(overlay, 'PTKImage')
                overlay_handle = overlay;
                if isempty(title) && ~isempty(image.Title)
                    title = image.Title;
                end
            else
                overlay_handle = PTKImage(overlay);
            end

            if isempty(title)
                title = 'PTK Viewer';
            end

            % Call the base class to initialise the hidden window
            obj = obj@GemFigure(title, [], reporting);

            obj.ViewerPanelHandle = MimViewerPanelWithControlPanel(obj);
            obj.AddChild(obj.ViewerPanelHandle);

            obj.Image = image_handle;
            obj.Overlay = overlay_handle;

            obj.Resize([100 50 700 600]);

            % Create the figure
            obj.Show;

        end

        function CreateGuiComponent(obj, position)
            CreateGuiComponent@GemFigure(obj, position);

            obj.ViewerPanelHandle.BackgroundImage = obj.Image;
            obj.ViewerPanelHandle.OverlayImage = obj.Overlay;
            obj.ViewerPanelHandle.SliceNumber = max(1, round(obj.Image.ImageSize/2));
        end

        function Resize(obj, position)
            Resize@GemFigure(obj, position);

            parent_width_pixels = max(1, position(3));
            parent_height_pixels = max(1, position(4));
            panel_position = [1 1 parent_width_pixels parent_height_pixels];
            obj.ViewerPanelHandle.Resize(panel_position);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function RS_PF(names, endo, ind_deval, actualdata, forecasts, ind_feval, rs_pf_windowSize, pdate)
%% Rossi-Sekhposyan (JAE,2016) Fluctuation Rationality Test

% Window size for Rossi-Sekhposyan (JAE,2016) fluctuation rationality test,
% used in computing series of Wald statistics
% rs_pf_windowSize = fT-round(0.2*fT);

ind_var=find(strcmp((names(1,2:end)),endo(ind_deval(1))));

y = actualdata(ind_var,:)'-forecasts(ind_feval(1),:)'; % construct forecast errors
fT=length(y);
x = [ones(fT,1),forecasts(ind_feval(1),:)'];                 % construct matrix of regressors

result = bear.RS_PF_Rationality(rs_pf_windowSize,y,x,pdate);         % perform rationality test

ptruncdate = result.ptruncdate; % truncated date vector
resultt    = result.resultt;    % test statistic
cvones     = result.cvones;     % critical values

%% Plots

% Actual Data vs Median Forecasts
figure('Tag','BEARresults');
subplot(2,2,1);
plot(pdate,actualdata(ind_var,:), 'LineWidth',2)
hold on
plot(pdate,forecasts(ind_feval(1),:),':k','LineWidth',2)
hold off
hleg=legend('Actual','Median Forecast','Location','best');
set(hleg,'Fontsize',8, 'Fontname','Palatino Linotype');
title(strcat([endo(ind_feval) 'Actual and Forecasted Values']),'Fontsize',8,'Fontname','Times New Roman');
axis tight
set(gca,'Fontname','Times New Roman');

% Plot Test Statistic
subplot(2,2,2);
plot(ptruncdate,resultt,ptruncdate,cvones,'r:','LineWidth',3);
xlim([ptruncdate(1) ptruncdate(end)])
ylabel('Test Statistic','fontsize',8,'Fontname','Palatino Linotype');
xlabel('Test rejects forecast rationality if test line is outside critical value line.','fontsize',7,'Fontname','Palatino Linotype');
hleg=legend('Fluctuation Test','5% Critical Value','Location','best');
set(hleg,'Fontsize',8, 'Fontname','Palatino Linotype');
title('Rossi-Sekhposyan(JAE,2016) Fluctuation Rationality Test','Fontsize',8,'Fontname','Palatino Linotype');
set(gca,'Fontname','Palatino Linotype');


text:
%% Image Segmentation with Distance Transform and Watershed Algorithm
% Sample code showing how to segment overlapping objects using Laplacian
% filtering, in addition to Watershed and Distance Transformation.
%
% In this demo you will learn how to:
%
% * Use the OpenCV function |cv.filter2D| in order to perform some laplacian
%   filtering for image sharpening
% * Use the OpenCV function |cv.distanceTransform| in order to obtain the
%   derived representation of a binary image, where the value of each pixel is
%   replaced by its distance to the nearest background pixel
% * Use the OpenCV function |cv.watershed| in order to isolate objects in the
%   image from the background
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/d2/dbd/tutorial_distance_transform.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/ImgTrans/imageSegmentation.cpp>
% * <https://docs.opencv.org/3.2.0/d3/db4/tutorial_py_watershed.html>
%

%%
% load image
fname = fullfile(mexopencv.root(), 'test', 'cards.png');
if exist(fname, 'file') ~= 2
    disp('Downloading image...')
    url = 'https://cdn.rawgit.com/opencv/opencv/3.2.0/samples/data/cards.png';
    urlwrite(url, fname);
end
src = cv.imread(fname, 'Color',true);
imshow(src), title('source image')

%%
% change the background from white to black, since that will help later to
% extract better results during the use of Distance Transform
mask = all(src == 255, 3);
src(repmat(mask,[1 1 3])) = 0;
imshow(src), title('black background image')

%%
% create a kernel that we will use for accuting/sharpening our image
% (an approximation of second derivative, a quite strong kernel)
if true
    kernel = [1 1 1; 1 -8 1; 1 1 1];
else
    kernel = [0 1 0; 1 -4 1; 0 1 0];
end

%%
% do the laplacian filtering, with intermediate calculations performed in
% floating-point precision (given the kernel used, Laplacian image will likely
% have negative values which would otherwise get truncated too soon)
if true
    imgLaplacian = cv.filter2D(src, kernel, 'DDepth','single');
else
    imgLaplacian = cv.Laplacian(src, 'DDepth','single');
end
src = uint8(single(src) - imgLaplacian);
imshow(src), title('sharpened image')

%%
% create binary image from source image
bw = cv.cvtColor(src, 'RGB2GRAY');
bw = cv.threshold(bw, 'Otsu');
imshow(bw), title('binary')

%%
% perform the distance transform algorithm, and normalize it to [0,1] range
% so we can visualize and threshold it
D = cv.distanceTransform(bw);
D = cv.normalize(D, 'Alpha',0, 'Beta',1, 'NormType','MinMax');
imshow(D), title('distance transform image')

%%
% threshold distance image to obtain the peaks and dilate it a bit
% (these will be the markers for the foreground objects)
D = cv.threshold(D, 0.4, 'MaxValue',1);
D = cv.dilate(D);
imshow(D), title('peaks')
if ~mexopencv.isOctave()
    %HACK: transparency not yet implemented in Octave
    image('CData',src, 'AlphaData',0.3)
end

%%
% find total markers
D = uint8(D * 255);
contours = cv.findContours(D, 'Mode','External', 'Method','Simple');
contours = cellfun(@(C) cat(1,C{:}), contours, 'UniformOutput',false);
N = numel(contours);

%%
% create the marker image for the watershed algorithm and draw the markers
% (seed regions are marked with positive labels 1,2,3,etc. everything else 0)
markers = zeros(size(D), 'int32');

% foreground markers
for i=1:N
    markers = cv.drawContours(markers, contours, 'ContourIdx',i-1, ...
        'Color',i, 'Thickness','Filled');
end

% background marker
markers = cv.circle(markers, [5 5], 3, 'Color',N+1, 'Thickness','Filled');

%%
% perform the watershed algorithm
markers = cv.watershed(src, markers);

%%
% in output, pixels are assigned to regions (1,2,etc.), boundaries set to -1
L = uint8(markers);  % -1,1,2,..,N,N+1 -> 0,1,2,..,N,N+1
if ~mexopencv.isOctave() && mexopencv.require('stats')
    %HACK: tabulate in Octave behaves differently
    tabulate(L(:));  % 0=boundaries, 1:N=regions, N+1=background
end

%%
% show segmented objects color-coded according to labels
clr = [1 1 1; hsv(N); 0 0 0];
rgb = uint8(255 * ind2rgb(L, clr));
imshow(rgb), title('segmented image')

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function pdf = uniform_pdf ( x, a, b )

%*****************************************************************************80
%
%% UNIFORM_PDF evaluates the Uniform PDF.
%
%  Discussion:
%
%    The Uniform PDF is also known as the "Rectangular" or "de Moivre" PDF.
%
%  Formula:
%
%    PDF(X)(A,B) = 1 / ( B - A ) for A <= X <= B
%               = 0 otherwise
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameters of the PDF.
%    A < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < a | b < x )
    pdf = 0.0;
  else
    pdf = 1.0 / ( b - a );
  end

  return
end


text:
function ll = modelPointLogLikelihood(model, varargin)

% MODELPOINTLOGLIKELIHOOD Compute the log likelihood of a given point.
% FORMAT
% DESC computes the log likelihood of the given model.
% ARG model : the model for which the log likelihood is to be
% computed.
% ARG P1, P2 ... : additional arguments as required.
% RETURN ll : the log likelihood of the given data point.
%
% SEEALSO : modelLogLikelihood, modelCreate
%
% COPYRIGHT : Neil D. Lawrence, 2006

% MLTOOLS

fhandle = str2func([model.type 'PointLogLikelihood']);
ll = fhandle(model, varargin{:});

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function Population = WOF_optimiseByNSGAII(GlobalDummy, Population, evaluations, isDummy)
% -----------------------------------------------------------------------
%  Copyright (C) 2020 Heiner Zille
%
%  This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%  International License. (CC BY-NC-SA 4.0). To view a copy of this license,
%  visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or see the
%  pdf-file "License-CC-BY-NC-SA-4.0.pdf" that came with this code.
%
%  You are free to:
%  * Share ? copy and redistribute the material in any medium or format
%  * Adapt ? remix, transform, and build upon the material
%  Under the following terms:
%  * Attribution ? You must give appropriate credit, provide a link to the
%     license, and indicate if changes were made. You may do so in any reasonable
%     manner, but not in any way that suggests the licensor endorses you or your use.
%  * NonCommercial ? You may not use the material for commercial purposes.
%  * ShareAlike ? If you remix, transform, or build upon the material, you must
%    distribute your contributions under the same license as the original.
%  * No additional restrictions ? You may not apply legal terms or technological
%    measures that legally restrict others from doing anything the license permits.
%
%  Author of this Code:
%   Heiner Zille <heiner.zille@ovgu.de> or <heiner.zille@gmail.com>
%
%  This code is based on the following publications:
%
%  1) Heiner Zille
%     "Large-scale Multi-objective Optimisation: New Approaches and a Classification of the State-of-the-Art"
%     PhD Thesis, Otto von Guericke University Magdeburg, 2019
%     http://dx.doi.org/10.25673/32063
%
%  2) Heiner Zille and Sanaz Mostaghim
%     "Comparison Study of Large-scale Optimisation Techniques on the LSMOP Benchmark Functions"
%     IEEE Symposium Series on Computational Intelligence (SSCI), IEEE, Honolulu, Hawaii, November 2017
%     https://ieeexplore.ieee.org/document/8280974
%
%  3) Heiner Zille, Hisao Ishibuchi, Sanaz Mostaghim and Yusuke Nojima
%     "A Framework for Large-scale Multi-objective Optimization based on Problem Transformation"
%     IEEE Transactions on Evolutionary Computation, Vol. 22, Issue 2, pp. 260-275, April 2018.
%     http://ieeexplore.ieee.org/document/7929324
%
%  4) Heiner Zille, Hisao Ishibuchi, Sanaz Mostaghim and Yusuke Nojima
%     "Weighted Optimization Framework for Large-scale Mullti-objective Optimization"
%     Genetic and Evolutionary Computation Conference (GECCO), ACM, Denver, USA, July 2016
%     http://dl.acm.org/citation.cfm?id=2908979
%
%  This file is intended to work with the PlatEMO framework version 2.5.
%  Date of publication of this code: 06.04.2020
%  Last Update of this code: 06.04.2020
%  A newer version of this algorithm may be available. Please contact the author
%  or see http://www.ci.ovgu.de/Research/Codes.html.
%
% The files may have been modified in Feb 2021 by the authors of the Platemo framework to work with the Platemo 3.0 release.
% -----------------------------------------------------------------------

    maximum = currentEvaluations(GlobalDummy, isDummy) + evaluations;

    [~,FrontNo,CrowdDis] = WOF_NSGAIIEnvironmentalSelection(Population,GlobalDummy.N);

    %% Optimization
    while currentEvaluations(GlobalDummy, isDummy) < maximum
        MatingPool = TournamentSelection(2,GlobalDummy.N,FrontNo,-CrowdDis);
        if isDummy == true
            NewDec      = WOF_GA(GlobalDummy, Population(MatingPool));
            numberOfNewSolutions = size(NewDec,1);
            Offspring = [];
            for j = 1:numberOfNewSolutions
                Offspring   = [Offspring, WOF_WeightIndividual(NewDec,GlobalDummy)];
            end
        else
            Offspring   = OperatorGA(GlobalDummy, Population(MatingPool));
        end
        [Population,FrontNo,CrowdDis] = WOF_NSGAIIEnvironmentalSelection([Population,Offspring],GlobalDummy.N);

    end
end

function e = currentEvaluations(GlobalDummy, isDummy)
    if isDummy == true
        e = GlobalDummy.Global.FE;
    else
        e = GlobalDummy.FE;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function opt = cplex_options(overrides, mpopt)
%CPLEX_OPTIONS  Sets options for CPLEX.
%
%   OPT = CPLEX_OPTIONS
%   OPT = CPLEX_OPTIONS(OVERRIDES)
%   OPT = CPLEX_OPTIONS(OVERRIDES, FNAME)
%   OPT = CPLEX_OPTIONS(OVERRIDES, MPOPT)
%
%   Sets the values for the options struct normally passed to
%   CPLEXOPTIMSET.
%
%   Inputs are all optional, second argument must be either a string
%   (FNAME) or a struct (MPOPT):
%
%       OVERRIDES - struct containing values to override the defaults
%       FNAME - name of user-supplied function called after default
%           options are set to modify them. Calling syntax is:
%                   MODIFIED_OPT = FNAME(DEFAULT_OPT);
%       MPOPT - MATPOWER options struct, uses the following fields:
%           opf.violation    - used to set opt.simplex.tolerances.feasibility
%           verbose          - used to set opt.barrier.display,
%               opt.conflict.display, opt.mip.display, opt.sifting.display,
%               opt.simplex.display, opt.tune.display
%           cplex.lpmethod   - used to set opt.lpmethod
%           cplex.qpmethod   - used to set opt.qpmethod
%           cplex.opts       - struct containing values to use as OVERRIDES
%           cplex.opt_fname  - name of user-supplied function used as FNAME,
%               except with calling syntax:
%                   MODIFIED_OPT = FNAME(DEFAULT_OPT, MPOPT);
%           cplex.opt        - numbered user option function, if and only if
%               cplex.opt_fname is empty and cplex.opt is non-zero, the value
%               of cplex.opt_fname is generated by appending cplex.opt to
%               'cplex_user_options_' (for backward compatibility with old
%               MATPOWER option CPLEX_OPT).
%
%   Output is an options struct to pass to CPLEXOPTIMSET.
%
%   There are multiple ways of providing values to override the default
%   options. Their precedence and order of application are as follows:
%
%   With inputs OVERRIDES and FNAME
%       1. FNAME is called
%       2. OVERRIDES are applied
%   With inputs OVERRIDES and MPOPT
%       1. FNAME (from cplex.opt_fname or cplex.opt) is called
%       2. cplex.opts (if not empty) are applied
%       3. OVERRIDES are applied
%
%   Example:
%
%   If cplex.opt = 3, then after setting the default CPLEX options,
%   CPLEX_OPTIONS will execute the following user-defined function
%   to allow option overrides:
%
%       opt = cplex_user_options_3(opt, mpopt);
%
%   The contents of cplex_user_options_3.m, could be something like:
%
%       function opt = cplex_user_options_3(opt, mpopt)
%       opt.threads          = 2;
%       opt.simplex.refactor = 1;
%       opt.timelimit        = 10000;
%
%   For details on the available options, see the "Parameters of CPLEX"
%   section of the CPLEX documentation at:
%
%       http://www.ibm.com/support/knowledgecenter/SSSA5P
%
%   See also CPLEXLP, CPLEXQP, MPOPTION.

%   MP-Opt-Model
%   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%%-----  initialization and arg handling  -----
%% defaults
verbose = 1;
feastol = 1e-6;
fname   = '';

%% second argument
if nargin > 1 && ~isempty(mpopt)
    if ischar(mpopt)        %% 2nd arg is FNAME (string)
        fname = mpopt;
        have_mpopt = 0;
    else                    %% 2nd arg is MPOPT (MATPOWER options struct)
        have_mpopt = 1;
        %% (make default opf.violation correspond to default CPLEX feastol)
        feastol  = mpopt.opf.violation/5;
        verbose  = mpopt.verbose;
        if isfield(mpopt.cplex, 'opt_fname') && ~isempty(mpopt.cplex.opt_fname)
            fname = mpopt.cplex.opt_fname;
        elseif mpopt.cplex.opt
            fname = sprintf('cplex_user_options_%d', mpopt.cplex.opt);
        end
    end
else
    have_mpopt = 0;
end

%%-----  set default options for CPLEX  -----
opt = struct( ...
    'simplex', struct('tolerances', struct('feasibility', feastol)), ...
    'output', struct('clonelog', -1) ...
);

%% printing
if verbose > 2
    opt.display = 'iter';
elseif verbose > 1
    opt.display = 'on';
elseif verbose > 0
    opt.display = 'off';
end

%% solution algorithm
if have_mpopt
    opt.lpmethod = mpopt.cplex.lpmethod;
    opt.qpmethod = mpopt.cplex.qpmethod;
end

%%-----  call user function to modify defaults  -----
if ~isempty(fname)
    if have_mpopt
        opt = feval(fname, opt, mpopt);
    else
        opt = feval(fname, opt);
    end
end

%%-----  apply overrides  -----
if have_mpopt && isfield(mpopt.cplex, 'opts') && ~isempty(mpopt.cplex.opts)
    opt = nested_struct_copy(opt, mpopt.cplex.opts);
end
if nargin > 0 && ~isempty(overrides)
    opt = nested_struct_copy(opt, overrides);
end

%--------------------------  Default Options Struct  --------------------------
% as returned by ...
%   >> opt = cplexoptimset('cplex')
%
%   opt =
%       advance:        1
%       barrier:        [1x1 struct]
%           algorithm:      0
%           colnonzeros:    0
%           convergetol:    1.0000e-08
%           crossover:      0
%           display:        1
%           limits:         [1x1 struct]
%               corrections:    -1
%               growth:         1.0000e+12
%               iteration:      9.2234e+18
%               objrange:       1.0000e+20
%           ordering:       0
%           qcpconvergetol: 1.0000e-07
%           startalg:       1
%       clocktype:      2
%       conflict:       [1x1 struct]
%           display:        1
%       diagnostics:    'off'
%       emphasis:       [1x1 struct]
%           memory:         0
%           mip:            0
%           numerical:      0
%       exportmodel:    ''
%       feasopt:        [1x1 struct]
%           mode:           0
%           tolerance:      1.0000e-06
%       lpmethod:       0
%       mip:            [1x1 struct]
%           cuts:           [1x1 struct]
%               cliques:        0
%               covers:         0
%               disjunctive:    0
%               flowcovers:     0
%               gomory:         0
%               gubcovers:      0
%               implied:        0
%               mcfcut:         0
%               mircut:         0
%               pathcut:        0
%               zerohalfcut:    0
%           display:        2
%           interval:       0
%           limits:         [1x1 struct]
%               aggforcut:      3
%               auxrootthreads: 0
%               cutpasses:      0
%               cutsfactor:     4
%               eachcutlimit:   2.1000e+09
%               gomorycand:     200
%               gomorypass:     0
%               nodes:          9.2234e+18
%               polishtime:     0
%               populate:       20
%               probetime:      1.0000e+75
%               repairtries:    0
%               solutions:      9.2234e+18
%               strongcand:     10
%               strongit:       0
%               submipnodelim:  500
%               treememory:     1.0000e+75
%           ordertype:      0
%           polishafter:    [1x1 struct]
%               absmipgap:      0
%               mipgap:         0
%               nodes:          9.2234e+18
%               solutions:      9.2234e+18
%               time:           1.0000e+75
%           pool:           [1x1 struct]
%               absgap:         1.0000e+75
%               capacity:       2.1000e+09
%               intensity:      0
%               relgap:         1.0000e+75
%               replace:        0
%           strategy:       [1x1 struct]
%               backtrack:      0.9999
%               bbinterval:     7
%               branch:         0
%               dive:           0
%               file:           1
%               fpheur:         0
%               heuristicfreq:  0
%               kappastats:     0
%               lbheur:         0
%               miqcpstrat:     0
%               nodeselect:     1
%               order:          1
%               presolvenode:   0
%               probe:          0
%               rinsheur:       0
%               search:         0
%               startalgorithm: 0
%               subalgorithm:   0
%               variableselect: 0
%           tolerances:     [1x1 struct]
%               absmipgap:      1.0000e-06
%               integrality:    1.0000e-05
%               lowercutoff:    -1.0000e+75
%               mipgap:         1.0000e-04
%               objdifference:  0
%               relobjdifference: 0
%               uppercutoff:    1.0000e+75
%       output:         [1x1 struct]
%           clonelog:       1
%           intsolfileprefix: ''
%           mpslong:        1
%           writelevel:     0
%       parallel:       0
%       preprocessing:  [1x1 struct]
%           aggregator:     -1
%           boundstrength:  -1
%           coeffreduce:    -1
%           dependency:     -1
%           dual:           0
%           fill:           10
%           linear:         1
%           numpass:        -1
%           presolve:       1
%           qpmakepsd:      1
%           reduce:         3
%           relax:          -1
%           repeatpresolve: -1
%           symmetry:       -1
%       qpmethod:       0
%       read:           [1x1 struct]
%           apiencoding:    ''
%           constraints:    30000
%           datacheck:      0
%           fileencoding:   'ISO-8859-1'
%           nonzeros:       250000
%           qpnonzeros:     5000
%           scale:          0
%           variables:      60000
%       sifting:        [1x1 struct]
%           algorithm:      0
%           display:        1
%           iterations:     9.2234e+18
%       simplex:        [1x1 struct]
%           crash:          1
%           dgradient:      0
%           display:        1
%           limits:         [1x1 struct]
%               iterations:     9.2234e+18
%               lowerobj:       -1.0000e+75
%               perturbation:   0
%               singularity:    10
%               upperobj:       1.0000e+75
%           perturbation:   [1x1 struct]
%               indicator:      0
%               constant:       1.0000e-06
%           pgradient:      0
%           pricing:        0
%           refactor:       0
%           tolerances:     [1x1 struct]
%               feasibility:    1.0000e-06
%               markowitz:      0.0100
%               optimality:     1.0000e-06
%       solutiontarget: 0
%       threads:        0
%       timelimit:      1.0000e+75
%       tune:           [1x1 struct]
%           display:        1
%           measure:        1
%           repeat:         1
%           timelimit:      10000
%       workdir:        '.'
%       workmem:        128

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower