text:
function [hpol,radlab_h] = polarLabels(theta,rho,line_style,labelrotate,rlabshift,clear)
% [hpol,radlab_h] = polarLabels(theta,rho,line_style,labelrotate,rlabshift,clear)
%
% polarLabels   Polar coordinate plot. (with added options)
%   polarLabels(THETA, RHO) makes a plot using polar coordinates of
%   the angle THETA, in radians, versus the radius RHO.
%   polarLabels(THETA,RHO,S) uses the linestyle specified in string S.
%   See PLOT for a description of legal linestyles.
%
%   See also PLOT, LOGLOG, SEMILOGX, SEMILOGY.
%
% Additions
%   Plot is labelled +/-180 instead of 0-360
%   Additional input options
%       labelrotate     value (in degrees) to add to angle labels to facilitate rotation of LABEL
%                       def = 0
%       rlabshift       value to add to radial labels to facilitate impression of negative values
%                       def = 0
%       clear           make polar background clear
%                       def = false (0)
%
% Ex. polarlabels([0:360]*pi/180,abs(sin([0:360]*pi/180))*10,'b',-90,-10)
%
%   Additional output option
%       [hpol,radlab_h] = polarLabels(theta,rho,line_style,labelrotate,rlabshift,clear)
%       radlab_h           vector of handles to radial text labels
%
%       Backward compatible with standard POLAR
%       For multiple plots (hold on) should only need to be used with the first one

%   Revised by BFGK 26-mars-2002
%    For better access to labels
%   Revised by BFGK 27-mai-2002
%    Fix to work with full 360 plot
%   Revised by BFGK 22-mai-2003
%    Added "clear" option for special background application
%    Tidied up nargin section

if nargin < 1
    error('Requires at least 2 input arguments.')
elseif nargin == 2
    if isstr(rho)
        line_style = rho;
        rho = theta;
        [mr,nr] = size(rho);
        if mr == 1
            theta = 1:nr;
        else
            th = (1:mr)';
            theta = th(:,ones(1,nr));
        end
    else
        line_style = 'auto';
    end
elseif nargin == 1
    line_style = 'auto';
    rho = theta;
    [mr,nr] = size(rho);
    if mr == 1
        theta = 1:nr;
    else
        th = (1:mr)';
        theta = th(:,ones(1,nr));
    end
end

if ~isequal(size(theta),size(rho))
    error('THETA and RHO must be the same size.');
end

% Dealing with additional arguments
if nargin < 4,  labelrotate = 0 ;    end
if nargin < 5,  rlabshift = 0   ;    end
if nargin < 6,  clear = 0       ;    end

if isstr(theta) | isstr(rho) | isstr(labelrotate) | isstr(rlabshift)
    error('Input arguments must be numeric.');
end

% get hold state
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;

% get x-axis text color so grid is in same color
tc = get(cax,'xcolor');
ls = get(cax,'gridlinestyle');

% Hold on to current Text defaults, reset them to the
% Axes' font attributes so tick marks use them.
fAngle  = get(cax, 'DefaultTextFontAngle');
fName   = get(cax, 'DefaultTextFontName');
fSize   = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
fUnits  = get(cax, 'DefaultTextUnits');
set(cax, 'DefaultTextFontAngle',  get(cax, 'FontAngle'), ...
    'DefaultTextFontName',   get(cax, 'FontName'), ...
    'DefaultTextFontSize',   get(cax, 'FontSize'), ...
    'DefaultTextFontWeight', get(cax, 'FontWeight'), ...
    'DefaultTextUnits','data')

% only do grids if hold is off
if ~hold_state

% make a radial grid
    hold on;
    maxrho = max(abs(rho(:)));
    hhh=plot([-maxrho -maxrho maxrho maxrho],[-maxrho maxrho maxrho -maxrho]);
    set(gca,'dataaspectratio',[1 1 1],'plotboxaspectratiomode','auto')
    v = [get(cax,'xlim') get(cax,'ylim')];
    ticks = sum(get(cax,'ytick')>=0);
    delete(hhh);
% check radial limits and ticks
    rmin = 0; rmax = v(4); rticks = max(ticks-1,2);
    if rticks > 5   % see if we can reduce the number
        if rem(rticks,2) == 0
            rticks = rticks/2;
        elseif rem(rticks,3) == 0
            rticks = rticks/3;
        end
    end

% define a circle
    th = 0:pi/50:2*pi;
    xunit = cos(th);
    yunit = sin(th);
% now really force points on x/y axes to lie on them exactly
    inds = 1:(length(th)-1)/4:length(th);
    xunit(inds(2:2:4)) = zeros(2,1);
    yunit(inds(1:2:5)) = zeros(3,1);
% plot background if necessary
    if (~isstr(get(cax,'color')) & ~clear),
       patch('xdata',xunit*rmax,'ydata',yunit*rmax, ...
             'edgecolor',tc,'facecolor',get(gca,'color'),...
             'handlevisibility','off');
    end

% draw radial circles
    c82 = cos(82*pi/180);
    s82 = sin(82*pi/180);
    rinc = (rmax-rmin)/rticks;
    cnt = 1;
    for i=(rmin+rinc):rinc:rmax
        hhh = plot(xunit*i,yunit*i,ls,'color',tc,'linewidth',1,...
                   'handlevisibility','off');
        radlab_h(cnt) = text((i+rinc/20)*c82,(i+rinc/20)*s82, ...
            ['  ' num2str(i+rlabshift)],'verticalalignment','bottom',...
            'handlevisibility','off');
        cnt = cnt+1;
    end
    set(hhh,'linestyle','-') % Make outer circle solid

% plot spokes
    th = (1:6)*2*pi/12;
    cst = cos(th); snt = sin(th);
    cs = [-cst; cst];
    sn = [-snt; snt];
    plot(rmax*cs,rmax*sn,ls,'color',tc,'linewidth',1,...
         'handlevisibility','off')

% annotate spokes in degrees
    rt = 1.1*rmax;
    for i = 1:length(th)
        angTextNum = (i*30)+labelrotate    ;
        text(rt*cst(i),rt*snt(i),int2str(angTextNum),...
             'horizontalalignment','center',...
             'handlevisibility','off');
        if i == length(th)
            loc = int2str(0+labelrotate);
        else
            angTextNum = 180+(i*30)+labelrotate    ;
            if angTextNum > 180, angTextNum = angTextNum - 360; end
            loc = int2str(angTextNum);
        end
        text(-rt*cst(i),-rt*snt(i),loc,'horizontalalignment','center',...
             'handlevisibility','off')
    end

% set view to 2-D
    view(2);
% set axis limits
    axis(rmax*[-1 1 -1.15 1.15]);
end

% Reset defaults.
set(cax, 'DefaultTextFontAngle', fAngle , ...
    'DefaultTextFontName',   fName , ...
    'DefaultTextFontSize',   fSize, ...
    'DefaultTextFontWeight', fWeight, ...
    'DefaultTextUnits',fUnits );

% transform data to Cartesian coordinates.
xx = rho.*cos(theta);
yy = rho.*sin(theta);

% plot data on top of grid
if strcmp(line_style,'auto')
    q = plot(xx,yy);
else
    q = plot(xx,yy,line_style);
end
if nargout > 0
    hpol = q;
end
if ~hold_state
    set(gca,'dataaspectratio',[1 1 1]), axis off; set(cax,'NextPlot',next);
end
set(get(gca,'xlabel'),'visible','on')
set(get(gca,'ylabel'),'visible','on')


text:
function opts = getdata(dataformat)

for i = 1:200
  imshow(reshape(images(:,i),144,176),[]);
  pause(0.1);
end

if strcmp(dataformat,'surveillance-video-Hall')
    load Hall_airport_1000_1496_497_144_176_gray.mat;
    D = images(:,1:200); [n1,n2] = size(images); % imn1 = 144; imn2 = 176;
    opts.D = D; opts.mu = norm(D)/1.25;
    Xs = D; Ys = D;

elseif strcmp(dataformat,'surveillance-video-Campus-color')
%     load Campus_trees_1000_1993_994_128_160_gray.mat;
    load Campus_trees_1000_1996_997_128_160_R.mat; DR = images(:,181:500);
    load Campus_trees_1000_1995_996_128_160_G.mat; DG = images(:,181:500);
    load Campus_trees_1000_1994_995_128_160_B.mat; DB = images(:,181:500);
    D = [DR, DG, DB];
    [n1,n2] = size(D); % imn1 = 128; imn2 = 160;
    opts.D = D; opts.mu = norm(D)/1.25;
    Xs = D; Ys = D;
end

opts.Xs = Xs;  opts.Ys = Ys;
opts.n1 = n1; opts.n2 = n2;
opts.sigma = 1e-6; opts.maxitr = 500; opts.rho = 1/sqrt(n1);
opts.eta_mu = 2/3; opts.eta_sigma = 2/3;
opts.muf = 1e-6;
opts.sigmaf = 1e-6;
opts.epsilon = 1e-7;
opts.sv = 100;

author: andrewssobral

repo: lrslibrary


text:
function [m,v] = spm_csd_chf(P,M,U)
% Characteristic (expected) frequency of a NMM
% FORMAT [G,w] = spm_csd_chf(P,M,U)
%
% P - parameters
% M - neural mass model structure
% U - trial-specific effects
%
% m - expected frequency
% v - dispersion
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_csd_chf.m 4718 2012-04-19 15:34:45Z karl $

% compute log-spectral density (without noise)
%==========================================================================
P.a(1) = -0;
P.a(2) = -32;
P.b    = P.b - 32;
P.c    = P.c - 32;

[G,w]  = spm_csd_mtf(P,M,U);

% compute moments (treating abs(G{i}) as a density)
%--------------------------------------------------------------------------
for i = 1:length(G)
    p     = abs(G{i});
    p     = p/sum(p);
    m(i)  = w'*p;
    v(i)  = ((w - m(i)).^2)'*p;
end

% average over trials
%--------------------------------------------------------------------------
m    = mean(m);
v    = mean(v);

author: spm


text:
classdef NSGAIIDTI < ALGORITHM
% <multi> <real/integer/label/binary/permutation> <constrained/none> <robust>
% NSGA-II of Deb's type I robust version

%------------------------------- Reference --------------------------------
% K. Deb and H. Gupta, Introducing robustness in multi-objective
% optimization, Evolutionary Computation, 2006, 14(4): 463-494.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

	methods
       function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();
            [PopObjV,PopConV]    = MeanEffective(Problem,Population);
            [~,FrontNo,CrowdDis] = EnvironmentalSelection(Population,Problem.N,PopObjV,PopConV);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,FrontNo,-CrowdDis);
                Offspring  = OperatorGA(Problem,Population(MatingPool),{1,10,1,50});
                [OffObjV,OffConV] = MeanEffective(Problem,Offspring);
                [Population,FrontNo,CrowdDis,PopObjV,PopConV] = EnvironmentalSelection([Population,Offspring],Problem.N,[PopObjV;OffObjV],[PopConV;OffConV]);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% Alignment of the Crystal Axes
%
% Default is $\vec c$ axis of highest symmetry.
%
% TODO: Explain the default setting in more detail.
%
%% Switching between different Alignment Options
%
% Since, especialy for lower symmetry groups, different conventions for
% aligning the crystal axes are used it might be necessary to transform
% data, e.g, orientations or tensors, from one convention into another.
% This can be done using the command <tensor.transformReferenceFrame.html
% transformReferenceFrame> as it illustrated below.
%
% First we import the stiffness tensor Forsterite with respect to the axes
% alignment

cs = crystalSymmetry('mmm',[4.7646 10.2296 5.9942],'mineral','Olivin');

% import some stiffness tensor
fname = fullfile(mtexDataPath,'tensor','Olivine1997PC.GPa');
C = stiffnessTensor.load(fname,cs)

plot(C)

%%
% Let us now consider a different setup of the Forsterite symmetry, where
% the $\vec a$ axis is the longest and the $\vec c$-axis is the shortest.

cs_new = crystalSymmetry('mmm',[10.2296 5.9942 4.7646],'mineral','Olivin')

%%
% In order to represent the stiffness tensor |C| with respect to this
% setupt we use the command <tensor.transformReferenceFrame.html
% transformReferenceFrame>.

C_new = C.transformReferenceFrame(cs_new)

nextAxis
plot(C_new)

author: mtex-toolbox

repo: mtex

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


text:
function[newPop] = normGeomSelect(oldPop,options)
% NormGeomSelect is a ranking selection function based on the normalized
% geometric distribution.
%
% function[newPop] = normGeomSelect(oldPop,options)
% newPop  - the new population selected from the oldPop
% oldPop  - the current population
% options - options to normGeomSelect [gen probability_of_selecting_best]

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 1, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

q=options(2); 				% Probability of selecting the best
e = size(oldPop,2); 			% Length of xZome, i.e. numvars+fit
n = size(oldPop,1); 			% Number of individuals in pop
newPop = zeros(n,e); 			% Allocate space for return pop
fit = zeros(n,1); 			% Allocates space for prob of select
x=zeros(n,2); 			        % Sorted list of rank and id
x(:,1) =[n:-1:1]'; 			% To know what element it was
[y x(:,2)] = sort(oldPop(:,e)); 	% Get the index after a sort
r = q/(1-(1-q)^n); 			% Normalize the distribution, q prime
fit(x(:,2))=r*(1-q).^(x(:,1)-1); 	% Generates Prob of selection
fit = cumsum(fit); 			% Calculate the cumulative prob. func
rNums=sort(rand(n,1)); 			% Generate n sorted random numbers
fitIn=1; newIn=1; 			% Initialize loop control
while newIn<=n 				% Get n new individuals
  if(rNums(newIn)<fit(fitIn))
    newPop(newIn,:) = oldPop(fitIn,:); 	% Select the fitIn individual
    newIn = newIn+1; 			% Looking for next new individual
  else
    fitIn = fitIn + 1; 			% Looking at next potential selection
  end
end

author: HuangCongQing


text:
function h=plotgauss2d(mu, Sigma)
% PLOTGAUSS2D Plot a 2D Gaussian as an ellipse with optional cross hairs
% h=plotgauss2(mu, Sigma)
%

h = plotcov2(mu, Sigma);
return;

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

% PLOTCOV2 - Plots a covariance ellipse with major and minor axes
%            for a bivariate Gaussian distribution.
%
% Usage:
%   h = plotcov2(mu, Sigma[, OPTIONS]);
%
% Inputs:
%   mu    - a 2 x 1 vector giving the mean of the distribution.
%   Sigma - a 2 x 2 symmetric positive semi-definite matrix giving
%           the covariance of the distribution (or the zero matrix).
%
% Options:
%   'conf'    - a scalar between 0 and 1 giving the confidence
%               interval (i.e., the fraction of probability mass to
%               be enclosed by the ellipse); default is 0.9.
%   'num-pts' - the number of points to be used to plot the
%               ellipse; default is 100.
%
% This function also accepts options for PLOT.
%
% Outputs:
%   h     - a vector of figure handles to the ellipse boundary and
%           its major and minor axes
%
% See also: PLOTCOV3

% Copyright (C) 2002 Mark A. Paskin

function h = plotcov2(mu, Sigma, varargin)

if size(Sigma) ~= [2 2], error('Sigma must be a 2 by 2 matrix'); end
if length(mu) ~= 2, error('mu must be a 2 by 1 vector'); end

[p, ...
 n, ...
 plot_opts] = process_options(varargin, 'conf', 0.9, ...
					'num-pts', 100);
h = [];
holding = ishold;
if (Sigma == zeros(2, 2))
  z = mu;
else
  % Compute the Mahalanobis radius of the ellipsoid that encloses
  % the desired probability mass.
  k = conf2mahal(p, 2);
  % The major and minor axes of the covariance ellipse are given by
  % the eigenvectors of the covariance matrix.  Their lengths (for
  % the ellipse with unit Mahalanobis radius) are given by the
  % square roots of the corresponding eigenvalues.
  if (issparse(Sigma))
    [V, D] = eigs(Sigma);
  else
    [V, D] = eig(Sigma);
  end
  % Compute the points on the surface of the ellipse.
  t = linspace(0, 2*pi, n);
  u = [cos(t); sin(t)];
  w = (k * V * sqrt(D)) * u;
  z = repmat(mu, [1 n]) + w;
  % Plot the major and minor axes.
  L = k * sqrt(diag(D));
  h = plot([mu(1); mu(1) + L(1) * V(1, 1)], ...
	   [mu(2); mu(2) + L(1) * V(2, 1)], plot_opts{:});
  hold on;
  h = [h; plot([mu(1); mu(1) + L(2) * V(1, 2)], ...
	       [mu(2); mu(2) + L(2) * V(2, 2)], plot_opts{:})];
end

h = [h; plot(z(1, :), z(2, :), plot_opts{:})];
if (~holding) hold off; end

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

% CONF2MAHAL - Translates a confidence interval to a Mahalanobis
%              distance.  Consider a multivariate Gaussian
%              distribution of the form
%
%   p(x) = 1/sqrt((2 * pi)^d * det(C)) * exp((-1/2) * MD(x, m, inv(C)))
%
%              where MD(x, m, P) is the Mahalanobis distance from x
%              to m under P:
%
%                 MD(x, m, P) = (x - m) * P * (x - m)'
%
%              A particular Mahalanobis distance k identifies an
%              ellipsoid centered at the mean of the distribution.
%              The confidence interval associated with this ellipsoid
%              is the probability mass enclosed by it.  Similarly,
%              a particular confidence interval uniquely determines
%              an ellipsoid with a fixed Mahalanobis distance.
%
%              If X is an d dimensional Gaussian-distributed vector,
%              then the Mahalanobis distance of X is distributed
%              according to the Chi-squared distribution with d
%              degrees of freedom.  Thus, the Mahalanobis distance is
%              determined by evaluating the inverse cumulative
%              distribution function of the chi squared distribution
%              up to the confidence value.
%
% Usage:
%
%   m = conf2mahal(c, d);
%
% Inputs:
%
%   c    - the confidence interval
%   d    - the number of dimensions of the Gaussian distribution
%
% Outputs:
%
%   m    - the Mahalanobis radius of the ellipsoid enclosing the
%          fraction c of the distribution's probability mass
%
% See also: MAHAL2CONF

% Copyright (C) 2002 Mark A. Paskin

function m = conf2mahal(c, d)

m = chi2inv(c, d); % matlab stats toolbox

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ftps2] = mGal2ftps2(mGal)
% Convert acceleration from milligals to feet per second-squared
% Chad A. Greene 2012
ftps2 = mGal*3.28083333e-5;
end

author: Sable

repo: mcbench-benchmarks


text:
function line = fitLine3d(points)
%FITLINE3D Fit a 3D line to a set of points.
%
%   LINE = fitLine3d(PTS)
%
%   Example
%     pts = randn(300, 3);
%     pts = transformPoint3d(pts, createScaling3d([6 4 2]));
%     pts = transformPoint3d(pts, createRotationOx(pi/6));
%     pts = transformPoint3d(pts, createRotationOy(pi/4));
%     pts = transformPoint3d(pts, createRotationOz(pi/3));
%     pts = transformPoint3d(pts, createTranslation3d([5 4 3]));
%     elli = inertiaEllipsoid(pts);
%     figure; drawPoint3d(pts); axis equal;
%     hold on; drawEllipsoid(elli, ...
%         'drawEllipses', true, 'EllipseColor', 'b', 'EllipseWidth', 3);
%     line = fitLine3d(pts);
%     drawLine3d(line, 'color', 'm', 'LineWidth', 4);
%
%   See also
%     lines3d, inertiaEllipsoid, fitPlane
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2012-11-11,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

% number of points
n = size(points, 1);

% compute centroid
center = mean(points);

% compute the covariance matrix
covPts = cov(points)/n;

% perform a principal component analysis with 2 variables,
% to extract inertia axes
[U, S] = svd(covPts);

% sort axes from greater to lower
[dummy, ind] = sort(diag(S), 'descend'); %#ok<ASGLU>

% format U to ensure first axis points to positive x direction
U = U(ind, :);
if U(1,1) < 0
    U = -U;
    % keep matrix determinant positive
    U(:,3) = -U(:,3);
end

line = [center U(:,1)'];

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ix] = spm_fieldindices(X,varargin)
% Returns the indices of fields in a structure (and vice versa)
% FORMAT [i]     = spm_fieldindices(X,field1,field2,...)
% FORMAT [field] = spm_fieldindices(X,i)
%
% X         - structure
% field1,.. - fields
%
% i         - vector of indices or fieldname{s}
%
% Note: Fields are returned in column order of X, regardless of the order
% of fields specified in the input.
%
%__________________________________________________________________________
% Copyright (C) 2010-2013 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fieldindices.m 7031 2017-02-27 19:46:48Z karl $

% if varargin is a vector simply return fieldnames
%--------------------------------------------------------------------------
if isnumeric(varargin{1})
    if numel(varargin{1}) > 1
        for j = 1:length(varargin{1})
            ix{j} = spm_fieldindices(X,varargin{1}(j));
        end
        return
    end
end

% create structure of zeros
%--------------------------------------------------------------------------
X0    = spm_zeros(X);
ix    = spm_vec(X0);

% and add one to specified fields
%--------------------------------------------------------------------------
for i = 1:length(varargin)

    if ischar(varargin{i})

        field = varargin{i};
        x  = X0;
        try
            % field is a field name
            %--------------------------------------------------------------
            f  = x.(field);
            f  = spm_unvec(spm_vec(f) + 1,f);
            x.(field) = f;
            ix = ix + spm_vec(x);

        catch
            try
                % field is an expression
                %----------------------------------------------------------
                eval(['x.' field ' = x.' field ' + 1;']);
                ix = ix + spm_vec(x);
            end
        end

        % or return the name of the field
        %----------------------------------------------------------------------
    elseif isnumeric(varargin{i})

        ind    = varargin{i};
        x      = ix;
        x(ind) = 1;
        x      = spm_unvec(x,X);
        name   = fieldnames(x);

        for j = 1:length(name)
            if any(spm_vec(x.(name{j})))
                if iscell(x.(name{j}))
                    for k = 1:numel(x.(name{j}))
                        xv = spm_vec(x.(name{j}){k});
                        if any(xv)
                            if isstruct(x.(name{j}){k})
                                s     = spm_fieldindices(x.(name{j}){k},find(xv));
                                ix    = sprintf('%s{%i}.%s',name{j},k,s);
                            else
                                [p,q] = find(x.(name{j}){k});
                                ix    = sprintf('%s{%i}(%i,%i)',name{j},k,p,q);
                            end
                            return
                        end
                    end
                elseif isnumeric(x.(name{j}))
                    s   = size(x.(name{j}));
                    if numel(s) < 3
                        if min(s) == 1
                            p    = find(x.(name{j}));
                            ix   = sprintf('%s(%i)',name{j},p);
                            return
                        else
                            [p,q] = find(x.(name{j}));
                            ix    = sprintf('%s(%i,%i)',name{j},p,q);
                            return
                        end
                    else
                        [p,q,r] = ind2sub(s,find(x.(name{j})));
                        ix      = sprintf('%s(%i,%i,%i)',name{j},p,q,r);
                        return
                    end
                end
            end
        end
    end
end

% find indices
%--------------------------------------------------------------------------
ix = find(ix);


text:
%DEMO_LGCP  Demonstration for a log Gaussian Cox process
%           with inference via EP or Laplace approximation
%
%  Description
%    Log Gaussian Cox process (LGCP) is a model for non-homogeneous
%    point-process in which the log intensity is modelled using
%    Gaussian process. LGCP can be modelled using log GP and
%    Poisson observation model in a discretized space.
%
%    The model constructed is as follows:
%
%    The number of occurrences of the realised point pattern within cell w_i
%
%         y_i ~ Poisson(y_i| |w_i|exp(f_i))
%
%    where |w_i| is area of cell w_i and f_i is the log intensity.
%
%    We place a zero mean Gaussian process prior for f =
%    [f_1, f_2,...,f_n] ~ N(0, K),
%
%    where K is the covariance matrix, whose elements are given as
%    K_ij = k(x_i, x_j | th). The function k(x_i, x_j | th) is
%    covariance function and th its parameters. We place a
%    prior for parameters, p(th).
%
%    The inference is conducted via EP or Laplace, where we find
%    Gaussian approximation for p(f| th, data), where th is the
%    maximum a posterior (MAP) estimate for the parameters.
%
%  See also  LGCP, DEMO_SPATIAL2
%
% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% =====================================
% 1D-example
% =====================================
fprintf(['Coal disaster data with EP integration over the latent values\n'])

% Coal disaster data
S = which('demo_lgcp');
L = strrep(S,'demo_lgcp.m','demodata/coal.txt');
x = load(L);
xt = (1850:1963)';
[p1,pq1] = lgcp(x,xt,'gpcf',@gpcf_exp);

figure()
hp=patch([xt; xt(end:-1:1)],[pq1(:,1); pq1(end:-1:1,2)],[.9 .9 .9]);
set(hp,'edgecolor',[.9 .9 .9])
xlim([min(x) max(x)])
line(xt,p1,'linewidth',2);
line([x x],[5 5.3],'color','k')
line(xlim,[5.15 5.15],'color','k')
xlim([1850 1963])
ylim([0 5.29])
title('The coal mine disaster data, estimated intensity, and 90% interval')
xlabel('Year')
ylabel('Intensity')

% =====================================
% 2D-example
% =====================================
fprintf(['Redwood data with Laplace integration over the latent\n' ...
         'values and MAP estimate for the parameters\n'])

% Redwood data
S = which('demo_lgcp');
L = strrep(S,'demo_lgcp.m','demodata/redwoodfull.txt');
x=load(L);
x1min=min(x(:,1));x1max=max(x(:,1));
x2min=min(x(:,2));x2max=max(x(:,2));
[xt1,xt2]=meshgrid(linspace(x1min,x1max,100),...
                   linspace(x2min,x2max,100));
xt=[xt1(:) xt2(:)];
p2 = lgcp(x,xt,'range',[0 1 0 1],'latent_method','Laplace','gridn',20);

figure()
G=zeros(size(xt1));
G(:)=p2;
pcolor(xt1,xt2,G);
shading flat
colormap('jet')
cx=caxis;
cx(1)=0;
caxis(cx);
colorbar
h=line(x(:,1),x(:,2),'marker','.','linestyle','none','color','k','markersize',10);
colorbar
axis square
set(gca,'xtick',[0:.2:1],'ytick',[0:.2:1])
title('Redwood data and intensity estimate')

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function test_glm
%Validate calculation of general linear model
%
%  test_glm()
%
% Tests: applyGlm, computeContrastMap2
%
% INPUTS
%  No inputs
%
% RETURNS
%  No returns
%
% Example: test_glm()
%
% See also MRVTEST
%
% Copyright Stanford team, mrVista, 2013

%% Initialize the key variables and data path
% Data directory (where the mrSession file is located)
dataDir = mrtInstallSampleData('functional','vwfaLoc');

% This is the validation file
storedGLM = mrtGetValididationData('glm');

% These are the items we stored in the validation file
%
% val.dim             = size(contrastMap);
% val.contrastMean    = nanmean(contrastMap(:));
% val.contrastMed     = median(contrastMap(:));
% val.contrastMax     = max(contrastMap(:));
% val.contrastVoc     = contrastMap(40,40,10);
%
% save(vFile, '-struct', 'val')

%% Retain original directory, change to data directory
curDir = pwd;
cd(dataDir);

% There can be several data types - name the one you want to use for
% computing GLM.
dataType = 'MotionComp';

% We will compute a GLM for 2 scans ...
scans = 1:2;

%% Get data structure and calculate glm

% open a session
vw = initHiddenInplane;

% set to motion comp dataTYPES
vw = viewSet(vw, 'current data type', dataType);

% get the GLM params
params = er_getParams(vw, scans(1), dataType);

% run the GLM
vw = applyGlm(vw, dataType, scans, params, 'GLM');

% Compute a contrast map from command line
stim     = er_concatParfiles(vw);
active   = find(strcmpi(stim.condNames, 'Word'));
control  = find(strcmpi(stim.condNames, 'Fix'));
saveName = 'WordVFix';
vw       = computeContrastMap2(vw, active, control, saveName);

contrastMap = viewGet(vw, 'map', 1); contrastMap = contrastMap{1};

%% Return to original directory
cd(curDir)

%% Validate the results

assertEqual(storedGLM.dim, size(contrastMap));
assertElementsAlmostEqual(storedGLM.contrastMean, nanmean(contrastMap(:)), 'relative', .00001);
assertElementsAlmostEqual(storedGLM.contrastMed,  median(contrastMap(:)), 'relative', .00001);
assertElementsAlmostEqual(storedGLM.contrastMax,  max(contrastMap(:)), 'relative', .00001);
assertElementsAlmostEqual(storedGLM.contrastVoc,  contrastMap(40,40,10), 'relative', .00001);

mrvCleanWorkspace;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function SimRndResults = SimRnd(Model, RndParam, Opt)
%VaryRndParam Multi Voxel simulation of normally distributed parameters

fields = fieldnames(RndParam);
n = length(RndParam.(fields{1})); % number of voxels

% Create waitbar
h = waitbar(0, sprintf('Data 0/%0.0f',n), 'Name', 'Simulating data',...
    'CreateCancelBtn', 'if ~strcmp(get(gcbf,''Name''),''canceling...''), setappdata(gcbf,''canceling'',1); set(gcbf,''Name'',''canceling...''); else delete(gcbf); end');
setappdata(h,'canceling',0)
setappdata(0,'Cancel',0);

tic;
for ii = 1:n
    for ix=1:length(Model.xnames)
        x(ix) = RndParam.(Model.xnames{ix})(ii);
    end
    Fit = Model.Sim_Single_Voxel_Curve(x,Opt,0);

    if ii==1 % initialize structure at first Sim
        fields = fieldnames(Fit);
        for vox = 1:length(fields)
            SimRndResults.(fields{vox}) = nan(n,1);
        end
    end

    for jj = 1:length(fields)
        SimRndResults.(fields{jj})(ii) = Fit.(fields{jj})(1);
    end

    % Update waitbar
    if getappdata(h,'canceling');  break;  end
    waitbar(ii/n,h,sprintf('Data %0.0f/%0.0f',ii,n));
end

delete(h);
SimRndResults.time = toc
SimRndResults.fields = fields;
SimRndResults = AnalyzeResults(RndParam, SimRndResults);
end

function Results = AnalyzeResults(Input, Results)
Fields = intersect(fieldnames(Input), fieldnames(Results));
for ii = 1:length(Fields)
    n = length(Input.(Fields{ii}));
    Results.Error.(Fields{ii})    = Results.(Fields{ii}) - Input.(Fields{ii}) ;
    Results.PctError.(Fields{ii}) = 100*(Results.(Fields{ii}) - Input.(Fields{ii})) ./ Input.(Fields{ii});
    Results.MPE.(Fields{ii})      = 100/n*sum((Results.(Fields{ii}) - Input.(Fields{ii})) ./ Input.(Fields{ii}));
    Results.RMSE.(Fields{ii})     = sqrt(sum((Results.(Fields{ii}) - Input.(Fields{ii})).^2 )/n);
    Results.NRMSE.(Fields{ii})    = Results.RMSE.(Fields{ii}) / (max(Input.(Fields{ii})) - min(Input.(Fields{ii})));
end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function plotTree(nodes, cats, subTree, figLabelHierarchy)
% plotTree(nodes, cats, subTree, figLabelHierarchy)
%
% Plot a tree using the results of the CocoStuffClasses.getClassHierarchyX() functions.
%
% subTree: (optional) +-1 for left/right sub tree, 0 for the entire tree
% figLabelHierarchy: (optional) handle to a figure
%
% Copyright by Holger Caesar, 2017

% By default we plot the entire tree
if ~exist('subTree', 'var')
    subTree = 0;
end

% Create figure if necessary
if ~exist('figLabelHierarchy', 'var')
    figLabelHierarchy = figure();
end

% Check that tree is binary at the top node
firstChildren = find(nodes == 1);
assert(numel(firstChildren) == 2);

% Get only relevant nodes and cats
if subTree ~= 0
    % Find descendents of the specified startTreeInd node
    sel = false(size(nodes));
    if subTree == -1
        sel(firstChildren(1)) = true;
    elseif subTree == 1
        sel(firstChildren(2)) = true;
    end
    while true
        oldSel = sel;
        sel = sel | ismember(nodes, find(sel));
        if isequal(sel, oldSel)
            break;
        end
    end
    nodes = nodes(sel);
    cats = cats(sel);

    % Remap nodes in 0:x range
    map = false(max(nodes), 1);
    map(unique(nodes)) = true;
    map = cumsum(map)-1;
    nodes = map(nodes);
end

% Plot them
ax = axes('Parent', figLabelHierarchy, 'Units', 'Norm');
axis(ax, 'off');
treeplot(nodes');
moveLeft = 0.08;
if subTree == -1
    set(ax, 'Position', [0-moveLeft,   0, 0.5+moveLeft, 1]);
elseif subTree == 1
    set(ax, 'Position', [0.5-moveLeft, 0, 0.5+moveLeft, 1]);
end
[xs, ys] = treelayout(nodes);

% Set appearance settings and show labels
isLeaf = ys == min(ys);
textInner = text(xs(~isLeaf) + 0.01, ys(~isLeaf) - 0.025, cats(~isLeaf), 'VerticalAlignment', 'Bottom', 'HorizontalAlignment', 'right'); %#ok<NASGU>
textLeaf  = text(xs( isLeaf) - 0.01, ys( isLeaf) - 0.02,  cats( isLeaf), 'VerticalAlignment', 'Bottom', 'HorizontalAlignment', 'left'); %#ok<NASGU>
set(ax, 'XTick', [], 'YTick', [], 'Units', 'Normalized');
ax.XLabel.String = '';
axis off;

% Rotate view
camroll(90);

author: nightrome


text:
% check whether a variable is in GPU or not
function output = IsInGPU(x)
if strcmpi(class(x), 'gpuArray')
    output = 1;
else
    output = 0;
end
end

author: singaxiong

repo: SignalGraph


text:
function [output] = F_exp(input_layer)
input = input_layer.a;
output = exp(input);
end


text:
%%******************************************************************
%% randlowranksdp.m : creates random feasible SDP problems where the
%%                constraint matrices are low-rank matrices of the
%%                form V*diag(d)*V'.
%%
%% [blk2,At2,C,b,blk,At] = randlowranksdp(n,m1,m2,r);
%%
%% blk2,At2: data with low-rank structure coded.
%% blk, At:  data without taking low-rank structure into account.
%%
%% n = size of the sdp variable
%% m1 = number of general constraints
%% m2 = number of low-rank constraints
%% r = rank of each constraint matrix.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function [blk2,At2,C,b,blk,At] = randlowranksdp(n,m1,m2,r)

  blk = cell(1,2);
  if (m1 > 0)
     [blk,At0,C,b0] = randsdp(n,[],[],m1);
  else
     b0 = []; At0 = cell(1);
  end
  if (m2 > 0)
     if (nargout > 4)
        [blk2,At2,C,b2,blk,At1] = randlowranksdpfun(n,m2,r);
     else
        [blk2,At2,C,b2] = randlowranksdpfun(n,m2,r);
     end
  else
     b2 = []; At1 = cell(1); blk2 = blk;
  end
  b = [b0; b2];
  if (nargout > 4)
     At{1} = [At0{1}, At1{1}];
  end
  At2{1,1} = At0{1};
%%******************************************************************
%%******************************************************************
  function [blk2,At2,C,b,blk,At] = randlowranksdpfun(n,m,r)

  randn('state',0);

  blk{1,1} = 's'; blk{1,2} = n;
  blk2{1,1} = 's'; blk2{1,2} = n; blk2{1,3} = r*ones(1,m);
  %%
  %% construct data with low rank structure
  %%
  At2 = cell(1,3);
  b = zeros(m,1);
  X0 = randn(n); X0 = X0*X0'; X0 = 0.5*(X0+X0');
  ss = [0,cumsum(blk2{1,3})];
  V = randn(n,m*r);
  dd = [];
  for k = 1:length(blk2{1,3})
     idx = [ss(k)+1 : ss(k+1)];
     len = blk2{1,3}(k);
     Dk = randn(len,len); Dk = 0.5*(Dk+Dk');
     [ii,jj,vv] = find(Dk);
     numnz = length(ii);
     dd = [dd; k*ones(numnz,1),ii,jj,vv]; %% each row has the form [constr,i,j,val]
     tmp1 = X0*V(:,idx);
     tmp2 = V(:,idx)*Dk;
     b(k) = sum(sum(tmp1.*tmp2));
  end
  At2{1,2} = V; At2{1,3} = dd;
  Z0 = X0;
  y0 = randn(m,1);
  Aty = At2{1,2}*spdiags(mexexpand(blk2{1,3},y0),0,r*m,r*m)*At2{1,2}';
  C{1} = Z0 + norm(Aty,'fro')*speye(n,n) + Aty;
  C{1} = 0.5*(C{1}+C{1}');
  %%
  %% construct data without exploiting low rank structure
  %%
  if (nargout > 4)
     idxD = [0; find(diff(dd(:,1))); size(dd,1)];
     for k = 1:m
        idx = [ss(k)+1 : ss(k+1)];
        Vk = At2{1,2}(:,idx);
        len = blk2{1,3}(k);
        idx2 = [idxD(k)+1:idxD(k+1)];
        Dk = spconvert([dd(idx2,2:4); len,len,0]);
        A{k} = Vk*Dk*Vk';
     end
     At = svec(blk,A);
  end
%%******************************************************************

author: shamilmamedov


text:
function model = robOneDynamicsCreate(q, d, latentVals)

% ROBONEDYNAMICSCREATE Create the dynamics model.

% FGPLVM

if(q~=2)
  error('Robot One Dynamics designed for 2-D latent spaces.');
end
if(d~=q)
  error('Input dimension must equal output dimension');
end
model.a = 20;
model.sigma2 = (pi/16).^2;
model.mixTheta = 0.8;
model.mixR = 0.5;

model = robOneDynamicsSetLatentValues(model, latentVals);

% set b (the scale of the gamma distributions) so that a/b is the average jump.
aveR = mean(model.r);
model.b = model.a/aveR;
model.type = 'robOneDynamics';

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function loc = match_sorted_rows(a,s)
%MATCH_SORTED_ROWS   Location of matches in sorted rows.
% MATCH_SORTED_ROWS(A,S) returns LOC such that S(LOC(i),:) = A(i,:).
% If A(i,:) is not in S, LOC(i) = 0.
% S must be sorted and unique.
%
% This function is a special case of ISMEMBER_SORTED_ROWS.

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

[tf,loc] = ismember_sorted_rows(a,s);


text:
function test_failed = test_libltfat_ifft(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

[~,~,enuminfo]=libltfatprotofile;
phaseconv = enuminfo.ltfat_phaseconvention;

fftwflags = struct('FFTW_MEASURE',0,'FFTW_ESTIMATE',64,'FFTW_PATIENT',32,'FFTW_DESTROY_INPUT',1,...
    'FFTW_UNALIGNED',2,'FFTW_EXHAUSTIVE',8,'FFTW_PRESERVE_INPUT',16);

Larr  = [350 350   9   1];
Warr  = [  1   3   3   1];

    for idx = 1:numel(Larr)
        L = Larr(idx);
        W = Warr(idx);

        f = randn(L,W,flags.complexity) + 1i*randn(L,W,flags.complexity);
        fin = complex2interleaved(f);
        fPtr = libpointer(dataPtr,fin);

        c = cast(randn(L,W)+1i*randn(L,W),flags.complexity);
        cout = complex2interleaved(c);
        coutPtr = libpointer(dataPtr,cout);

        truec = ifft(f)*L;

        funname = makelibraryname('ifft',flags.complexity,0);
        status = calllib('libltfat',funname,fPtr,L,W,coutPtr);

        res = norm(truec - interleaved2complex(coutPtr.Value),'fro');
        [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
        fprintf(['ifft   L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

        % With plan
        c = cast(randn(L,W)+1i*randn(L,W),flags.complexity);
        cout = complex2interleaved(c);
        coutPtr = libpointer(dataPtr,cout);

        plan = libpointer();
        funname = makelibraryname('ifft_init',flags.complexity,0);
        statusInit = calllib('libltfat',funname,L,W,fPtr,coutPtr, fftwflags.FFTW_MEASURE, plan);

        funname = makelibraryname('ifft_execute',flags.complexity,0);
        statusExecute = calllib('libltfat',funname,plan);

        res = norm(truec - interleaved2complex(coutPtr.Value),'fro');
        [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
        fprintf(['ifft L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

        c = cast(randn(L,W)+1i*randn(L,W),flags.complexity);
        cout = complex2interleaved(c);
        coutPtr = libpointer(dataPtr,cout);

        funname = makelibraryname('ifft_execute_newarray',flags.complexity,0);
        statusExecute = calllib('libltfat',funname,plan,fPtr,coutPtr);

        res = norm(truec - interleaved2complex(coutPtr.Value),'fro');
        [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
        fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

        funname = makelibraryname('ifft_done',flags.complexity,0);
        statusDone = calllib('libltfat',funname,plan);

        %%%%%% Inplace
        c = f;
        cout = complex2interleaved(c);
        coutPtr = libpointer(dataPtr,cout);

        plan = libpointer();
        funname = makelibraryname('ifft_init',flags.complexity,0);
        statusInit = calllib('libltfat',funname,L,W, coutPtr, coutPtr, fftwflags.FFTW_MEASURE, plan);

        funname = makelibraryname('ifft_execute',flags.complexity,0);
        statusExecute = calllib('libltfat',funname,plan);

        res = norm(truec - interleaved2complex(coutPtr.Value),'fro');
        [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
        fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

        c = f;
        cout = complex2interleaved(c);
        coutPtr = libpointer(dataPtr,cout);

        funname = makelibraryname('ifft_execute_newarray',flags.complexity,0);
        statusExecute = calllib('libltfat',funname,plan,coutPtr,coutPtr);

        res = norm(truec - interleaved2complex(coutPtr.Value),'fro');
        [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
        fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

        funname = makelibraryname('ifft_done',flags.complexity,0);
        statusDone = calllib('libltfat',funname,plan);

    end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [result_for_kitti_eval] = generateDataForKITTIEvaluationFromFile(result_for_kitti_eval, fno, dnos, detections)

    sz = size(result_for_kitti_eval,1);

    for i = 1:length(dnos)

        result_for_kitti_eval(sz+i, 1)   = fno;
        result_for_kitti_eval(sz+i, 2)   = dnos(i);
        result_for_kitti_eval(sz+i, 3:6) = detections{fno}(i, 1:4);
        result_for_kitti_eval(sz+i, 7)   = 0.9;

    end

end

repo: MOTBeyondPixels


text:
classdef OptimizerNullSpace < Optimizer

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

    properties (Access = private)
        tau
        lineSearchTrials
        lineSearch
        costOld
        upperBound
        lowerBound
        tol = 1e-8
        nX
        hasConverged
        acceptableStep
        oldDesignVariable
        oldCost
        incrementalScheme
        hasFinished
        mOld
        meritNew
        nConstr
        meritGradient

        globalCost
        globalConstraint
        globalCostGradient
        globalMerit
        globalLineSearch
        globalDual
        globalDesignVar
    end

    methods (Access = public)

        function obj = OptimizerNullSpace(cParams)
            obj.initOptimizer(cParams);
            obj.init(cParams);
            obj.outputFunction.monitoring.create(cParams);
            obj.createPrimalUpdater(cParams);
            obj.createDualUpdater(cParams);
            obj.prepareFirstIter();
        end

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

    end

    methods(Access = private)

        function init(obj,cParams)
            obj.upperBound             = cParams.uncOptimizerSettings.ub;
            obj.lowerBound             = cParams.uncOptimizerSettings.lb;
            obj.cost                   = cParams.cost;
            obj.constraint             = cParams.constraint;
            obj.designVariable         = cParams.designVar;
            obj.dualVariable           = cParams.dualVariable;
            obj.incrementalScheme      = cParams.incrementalScheme;
            obj.nConstr                = cParams.constraint.nSF;
            obj.nX                     = length(obj.designVariable.value);
            obj.maxIter                = cParams.maxIter;
            obj.hasConverged           = false;
            obj.nIter                  = 0;
        end

        function prepareFirstIter(obj)
            obj.cost.computeFunctionAndGradient();
            obj.costOld = obj.cost.value;
            obj.designVariable.updateOld();
            obj.dualVariable.value = zeros(obj.nConstr,1);
        end

        function obj = update(obj)
            x0 = obj.designVariable.value;
            obj.saveOldValues(x0);
            if obj.nIter ~= 0
                obj.dualUpdater.update();
            end
            obj.mOld = obj.computeMeritFunction(x0);
            obj.calculateInitialStep();
            obj.acceptableStep   = false;
            obj.lineSearchTrials = 0;
            DJ = obj.cost.gradient;
            Dg = obj.constraint.gradient;
            while ~obj.acceptableStep
                obj.computeMeritGradient(DJ,Dg);
                x = obj.updatePrimal();
                obj.designVariable.update(x);
                obj.dualUpdater.update();
                obj.checkStep(x,x0);
            end
            obj.updateOldValues(x);
        end

        function obj = calculateInitialStep(obj)
            if obj.nIter == 0
                obj.cost.computeFunctionAndGradient();
                obj.constraint.computeFunctionAndGradient();
                x       = obj.designVariable.value;
                l       = obj.dualVariable.value;
                DJ      = obj.cost.gradient;
                Dg      = obj.constraint.gradient;
                aJ      = 1;
                DmF     = aJ*(DJ + Dg*l);
                factor  = 1;
                obj.primalUpdater.computeFirstStepLength(DmF,x,factor);
            else
                factor = 1.2;
                obj.primalUpdater.increaseStepLength(factor);
            end
        end

        function displayIter(obj,x)
            m = obj.designVariable.mesh;
            bm = m.createBoundaryMesh();
            s.backgroundMesh = m;
            s.boundaryMesh   = bm;
            um = UnfittedMesh(s);
            um.compute(x);
            figure()
            um.plot();
        end

        function x = updatePrimal(obj)
            x       = obj.designVariable.value;
            g       = obj.meritGradient;
            x       = obj.primalUpdater.update(g,x);
        end

        function computeMeritGradient(obj,DJ,Dg)
            l       = obj.dualVariable.value;
            aJ      = 1;
            DmF     = aJ*(DJ + Dg*l);
            obj.meritGradient = DmF;
        end

        function checkStep(obj,x,x0)
            mNew = obj.computeMeritFunction(x);
            if obj.nIter == 0 %&& mNew <= obj.mOld
                obj.acceptableStep = true;
                obj.meritNew       = mNew;
                obj.dualUpdater.updateOld();
            end
            if mNew < obj.mOld
                obj.acceptableStep = true;
                obj.meritNew = mNew;
                obj.dualUpdater.updateOld();
            elseif obj.primalUpdater.isTooSmall()
                error('Convergence could not be achieved (step length too small)')
            else
                obj.primalUpdater.decreaseStepLength();
                obj.designVariable.update(x0);
                obj.lineSearchTrials = obj.lineSearchTrials + 1;
            end
        end

        function mF = computeMeritFunction(obj,x)
            obj.designVariable.update(x)
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();
            J  = obj.cost.value;
            h  = obj.constraint.value;
            l  = obj.dualVariable.value;
            aJ = 1;
            AJ = aJ*(J + l'*h);
            mF = AJ;
        end

        function obj = saveOldValues(obj,x)
            obj.designVariable.update(x);
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();
            obj.oldCost            = obj.cost.value;
            obj.oldDesignVariable  = x;
        end

        function obj = updateOldValues(obj,x)
            obj.designVariable.update(x);
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();
        end

        function obj = checkConvergence(obj)
           if abs(obj.meritNew - obj.mOld) < obj.tol && obj.checkConstraint()
               obj.hasConverged = true;
           else

           end

        end

        function obj = updateMonitoring(obj)
            s.nIter            = obj.nIter;
            s.tau              = obj.primalUpdater.tau;
            s.lineSearch       = obj.lineSearch;
            s.lineSearchTrials = obj.lineSearchTrials;
            s.oldCost          = obj.oldCost;
            s.hasFinished      = obj.hasFinished;
            s.meritNew         = obj.meritNew;
            obj.outputFunction.monitoring.compute(s);
        end

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

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

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

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

        function saveVariablesForAnalysis(obj)
            i                           = obj.nIter + 1;
            obj.globalCost(i)           = obj.cost.value;
            obj.globalConstraint(:,i)   = obj.constraint.value;
            obj.globalCostGradient(i)   = norm(obj.cost.gradient);
%             obj.globalMerit(i)          = obj.meritNew;
%             obj.globalLineSearch(i)     = obj.primalUpdater.tau;
            obj.globalDual(:,i)         = obj.dualVariable.value;
            obj.globalDesignVar(:,i)    = obj.designVariable.value;
            if obj.hasConverged
                c = obj.globalCost;
                h = obj.globalConstraint;
                g = obj.globalCostGradient;
%                 m = obj.globalMerit;
%                 t = obj.globalLineSearch;
                d = obj.globalDual;
                v = obj.globalDesignVar;
                save('NullSpaceCant04.mat',"c","g","h","d","v");
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [M_est,U_est,V_est,L1_error] = RobustApproximation_M_UV_TraceNormReg(M,W,r,lambda,rho,maxIterIN,signM)
%% Robust low-rank matrix approximation with missing data and outliers
% min |W.*(M-E)|_1 + lambda*|V|_*
% s.t., E = UV, U'*U = I
%
%Input:
%   M: m*n data matrix
%   W: m*n indicator matrix, with '1' means 'observed', and '0' 'missing'.
%   r: the rank of r
%   lambda: the weighting factor of the trace-norm regularization, 1e-3 in default.
%   rho: increasing ratio of the penalty parameter mu, usually 1.05.
%   maxInterIN: maximum iteration number of inner loop, usually 100.
%   signM: if M>= 0, then signM = 1, otherwise, signM = 0;
%Output:
%   M_est: m*n full matrix, such that M_est = U_est*V_est
%   U_est: m*r matrix
%   V_est: r*n matrix
%   L1_error: the L1-norm error of observed data only.
%% Normalization
scale = max(max(abs(M)));
M = M/scale;

%% In-default parameters
[m n] = size(M); %matrix dimension
if nargin < 6
    maxIterIN = 100;
end
if nargin < 5
    rho = 1.05;
end
if nargin < 4
    lambda = 1e-3;
end
if nargin < 3
    disp('Please input the data matrix M, the indicator W and the rank r, and try again.');
end

maxIterOUT = 5000;
max_mu = 1e20;
mu = 1e-6;
M_norm = norm(M,'fro');
tol = 1e-8*M_norm;

cW = ones(size(W)) - W; %the complement of W.
display = 1; %display progress
%% Initializing optimization variables as zeros
E = zeros(m,n);
U = zeros(m,r);
V = zeros(r,n);
Y = zeros(m,n); %lagrange multiplier
%% Start main outer loop
iter_OUT = 0;
objs=[];
while iter_OUT < maxIterOUT
    iter_OUT = iter_OUT + 1;

    itr_IN = 0;
    obj_pre = 1e20;
    %start inner loop
    while itr_IN < maxIterIN
        %update U
        temp = (E + Y/mu)*V';

        [Us,sigma,Ud] = svd(temp,'econ'); % stable
        %[Us,sigma,Ud] = svdecon(temp); % fastest

        U = Us*Ud';

        %update V
        temp = U'*(E + Y/mu);

        [Vs,sigma,Vd] = svd(temp,'econ'); % stable
        %[Vs,sigma,Vd] = svdecon(temp); % fastest

        sigma = diag(sigma);
        svp = length(find(sigma > lambda/mu));
        if svp >= 1
            sigma = sigma(1:svp) - lambda/mu;
        else
            svp = 1;
            sigma = 0;
        end
        V = Vs(:,1:svp)*diag(sigma)*Vd(:,1:svp)';
        sigma0 = sigma;

        UV = U*V;

        %update E
        temp1 = UV - Y/mu;
        temp = M-temp1;
        E = max(0,temp - 1/mu) + min(0,temp + 1/mu);
        E = (M-E).*W + temp1.*cW;

        if signM > 0
            E(E<0) = 0;
        end

        %evaluate current objective
        obj_cur = sum(sum(abs(W.*(M-E)))) + lambda*sum(sigma0) + sum(sum(Y.*(E-UV))) + mu/2*norm(E-UV,'fro')^2;

        %check convergence of inner loop
        if abs(obj_cur - obj_pre) <= 1e-8*abs(obj_pre)
            break;
        else
            obj_pre = obj_cur;
            itr_IN = itr_IN + 1;
        end
    end

    leq = E - UV;
    stopC = norm(leq,'fro');
    if display
        obj = sum(sum(abs(W.*(M-UV)))) + lambda*sum(sigma0);
        objs = [objs,obj];
    end
    if display && (iter_OUT==1 || mod(iter_OUT,50)==0 || stopC<tol)
        disp(['iter ' num2str(iter_OUT) ',mu=' num2str(mu,'%2.1e') ...
            ',obj=' num2str(obj) ',stopALM=' num2str(stopC,'%2.3e')]);
    end
    if stopC<tol
        break;
    else
        %update lagrage multiplier
        Y = Y + mu*leq;
        %update penalty parameter
        mu = min(max_mu,mu*rho);
    end
end

%% Denormalization
U_est = sqrt(scale)*U; V_est = sqrt(scale)*V;
M_est = U_est*V_est;
L1_error = sum(sum(abs(W.*(scale*M-M_est))));

end

author: andrewssobral


text:
function f = constructor(f, op, varargin)
%CONSTRUCTOR   Main CHEBFUN3T constructor.
%   Classical full tensor approach for 3D functions. No low-rank technique
%   is involved here. This is experimental code, not for ordinary use.
%
% See also CHEBFUN3.

%   The 'trig' flag is NOT implemented.

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

% Parse the inputs:
[op, dom, pref, vectorize] = parseInputs(op, varargin{:});

% Set preferences:
tech        = pref.tech();
tpref       = tech.techPref;
grid = tpref.minSamples;
maxSample = 363;
% Since 363*363*363=(6916)^2, this is equivalent to working with a matrix
% of size 6916 * 6916 which needs 380MB of memory.
pseudoLevel = pref.cheb3Prefs.chebfun3eps;
prefStruct = pref.cheb2Prefs;
passSampleTest = prefStruct.sampleTest;

[xx, yy, zz] = ndgrid(dom(1:2), dom(3:4), dom(5:6));
A = op(xx, yy, zz);
if ( isscalar(A) )
    op = @(x,y,z) op(x,y,z) + 0*x + 0*y + 0*z;
end

%%
isHappy = 0;
m = grid; n = grid; p = grid;
for i=1:10
    %% Compute values at 3D chebpts
    out = tech.tensorGrid([m, n, p], dom);
    xx = out{1};
    yy = out{2};
    zz = out{3};
    F = op(xx,yy,zz);

    % Does the function blow up or evaluate to NaN?:
    vscale = max(abs(F(:)));
    if ( isinf(vscale) )
        error('CHEBFUN:CHEBFUN3T:constructor:inf', ...
            'Function returned INF when evaluated');
    elseif ( any(isnan(F(:)) ) )
        error('CHEBFUN:CHEBFUN3T:constructor:nan', ...
            'Function returned NaN when evaluated');
    end

    %% Call 3D vals2coeffs
    coeffs3D = chebfun3t.vals2coeffs(F);
    % happinessCheck
    [isHappyX, isHappyY, isHappyZ, cutoffX2, cutoffY2, cutoffZ2] = ...
        happinessCheck3D(coeffs3D, pref);
    isHappy = isHappyX & isHappyY & isHappyZ;

    %% PHASE II
    % This phase refines the grid ONLY in "unhappy" directions.
    % We don't refine here if ALL the directions are nonhappy (as mentioned
    % in the following). If, however, less than 3 directions are nonhappy,
    % we refine in those directions only.

    m2 = m; n2 = n; p2 = p;
    while (~isHappy)
        % We need to refine the grid. We refine in 2 different ways: If all
        % 3 directions are unhappy, the new grid is bigger in each
        % direction by a factor of sqrt(2). If 1 or 2 of the directions are
        % unhappy, we refine (only those unhappy directions) by a factor of
        % 2.

        if (~isHappyX && ~isHappyY && ~isHappyZ)
            % All directions are unhappy. To avoid getting an unncecessarily
            % big tensor, refine differently.
            break;
        end
        if (~isHappyX)
            m2 = gridRefinePhase2(m2, pref);
        end
        if (~isHappyY)
            n2 = gridRefinePhase2(n2, pref);
        end
        if (~isHappyZ)
            p2 = gridRefinePhase2(p2, pref);
        end
        out = tech.tensorGrid([m2, n2, p2], dom);
        xx = out{1};
        yy = out{2};
        zz = out{3};
        F = op(xx,yy,zz);
        coeffs3D = chebfun3t.vals2coeffs(F);
        [isHappyX, isHappyY, isHappyZ, cutoffX2, cutoffY2, cutoffZ2] = ...
            happinessCheck3D(coeffs3D,pref);
        isHappy = isHappyX & isHappyY & isHappyZ;
    end

     if ( isHappy )
         % Chop the tail:
         coeffs3D = coeffs3D(1:cutoffX2, 1:cutoffY2, 1:cutoffZ2);

         % Construct a CHEBFUN3T object:
         f.coeffs = coeffs3D;
         f.vscale = max(abs(F(:)));
         f.domain = dom;

         % Step 2: Apply a SAMPLETEST:
         [m2, n2, p2] = size(coeffs3D);
         if ( passSampleTest )
             tol2 = getTol3D(xx, yy, zz, F, length(xx), dom, pseudoLevel);
             pass = sampleTest(f, op, tol2, vectorize);
             if ( pass ) % :)
                 break
             end
         end
     end
    m = gridRefinePhase1(m);
    n = gridRefinePhase1(n);
    p = gridRefinePhase1(p);
end

% Reached maximum allowed tensor size and still unhappy?
if ( i == 10 && ~isHappy )
    % Throw a warning and construct from the latest unhappy approximation
    warning('CHEBFUN:CHEBFUN3T:constructor:notResolved', ...
                'Unresolved with maximum CHEBFUN3T length: %u.', maxSample);
    f.coeffs = coeffs3D;
    f.vscale = max(abs(F(:)));
    f.domain = dom;
end

end

%%
function [isHappyX, isHappyY, isHappyZ, cutoffX2, cutoffY2, cutoffZ2] = ...
    happinessCheck3D(simple_3D_coeffs, pref)
% An alternative technique: check whether all the rows, columns and tubes have been resolved.
% Somewhat similar to Chebfun2 constructor (but with all the rows, cols and tubes instead of the pivot ones)

% TODO: THIS DOES NOT PROPERLY WORK WITH TRIGs.

colChebtech = sum(chebfun3t.unfold(abs(simple_3D_coeffs), 1), 2);
fCol = chebtech2({[], colChebtech});
[isHappyX, cutoffX2] = happinessCheck(fCol, [], [], [], pref);

rowChebtech = sum(chebfun3t.unfold(abs(simple_3D_coeffs), 2), 2);
fRow = chebtech2({[], rowChebtech});
[isHappyY, cutoffY2] = happinessCheck(fRow, [], [], [], pref);

tubeChebtech = sum(chebfun3t.unfold(abs(simple_3D_coeffs), 3), 2);
fTube = chebtech2({[], tubeChebtech});
[isHappyZ, cutoffZ2] = happinessCheck(fTube, [], [], [], pref);
end
%%
function grid = gridRefinePhase1(grid)
% Grid refinement strategy in Phase 1. It does the same for all TECHS.
grid = floor(sqrt(2)^(floor(2*log2(grid)) + 1)) + 1;
end

%%
function grid = gridRefinePhase2(grid, pref)
% Grid refinement strategy for tech in Phase 2

% What tech am I based on?:
tech = pref.tech();

% What is the next grid size?
if ( isa(tech, 'chebtech2') )
    % Double sampling on tensor grid:
    grid = 2*grid-1;
elseif ( isa(tech, 'trigtech') )
    % Double sampling on tensor grid:
    grid = 2^(floor(log2(grid)+1));
elseif ( isa(tech, 'chebtech1') )
    grid = 3 * grid;
else
    error('CHEBFUN:CHEBFUN3:constructor:gridRefinePhase2:techType', ...
        'Technology is unrecognized.');
end

end
%%

function tol = getTol3D(xx, yy, zz, vals, grid, dom, pseudoLevel)
%See https://github.com/chebfun/chebfun/issues/1491
[m,n,p] = size(vals);
% Remove some edge values so that df_dx, df_dy and df_dz have the same size.
df_dx = diff(vals(:,1:n-1,1:p-1),1,1) ./ diff(xx(:,1:n-1,1:p-1),1,1); % xx changes in the first mode.
df_dy = diff(vals(1:m-1,:,1:p-1),1,2) ./ diff(yy(1:m-1,:,1:p-1),1,2); % yy changes row-wise (2nd mode).
df_dz = diff(vals(1:m-1,1:n-1,:),1,3) ./ diff(zz(1:m-1,1:n-1,:),1,3); % zz changes tube-wise (3rd mode).

J = max(max(abs(df_dx),abs(df_dy)), abs(df_dz));
Jac_norm = max(J(:)); % An approximation for the norm of the Jacobian over the whole domain.
vscale = max(abs(vals(:)));
tol = grid^(4/5) * max( abs(dom(:) ) ) * max(Jac_norm, vscale) * pseudoLevel;
end
%%
function [op, dom, pref, vectorize] = parseInputs(op, varargin)
vectorize = 0;
pref = chebfunpref();

if ( isa(op, 'char') )     % CHEBFUN3T( CHAR )
    op = str2op(op);
end

if nargin == 1
    dom = [-1 1 -1 1 -1 1];
elseif ( ( (nargin == 2) && ~any(strcmpi(varargin{1}, {'trig', 'periodic'})) ))
    % DOMAIN is specified
    dom = varargin{1};
elseif ( ( (nargin == 2) && any(strcmpi(varargin{1}, {'trig', 'periodic'})) ))
    % TECH is specified
    pref.tech = @trigtech;
    dom = [-1 1 -1 1 -1 1];
elseif ( ( (nargin == 3) && any(strcmpi(varargin{2}, {'trig', 'periodic'})) ))
    % DOMAIN and TECH are specified
    dom = varargin{1};
    pref.tech = @trigtech;
elseif ( ( (nargin == 3) && any(strcmpi(varargin{1}, {'trig', 'periodic'})) ))
    % TECH and DOMAIN are specified
    pref.tech = @trigtech;
    dom = varargin{3};
elseif ( (nargin == 3) && strcmpi(varargin{1}, 'eps') )
    % EPS is specified
    pref.chebfuneps = varargin{2};
    dom = [-1 1 -1 1 -1 1];
elseif ( (nargin == 4) && strcmpi(varargin{1}, 'eps') ...
        && any(strcmpi(varargin{3}, {'trig', 'periodic'})) )
    % EPS and TECH are specified
    pref.tech = @trigtech;
    pref.chebfuneps = varargin{3};
    dom = [-1 1 -1 1 -1 1];
elseif ( (nargin == 4) && strcmpi(varargin{2}, 'eps') && ...
        any(strcmpi(varargin{1}, {'trig', 'periodic'})) )
    % TECH and EPS are specified
    pref.tech = @trigtech;
    pref.chebfuneps = varargin{3};
    dom = [-1 1 -1 1 -1 1];
elseif ( (nargin == 4) && strcmpi(varargin{2}, 'eps') && ...
        ~any(strcmpi(varargin{1}, {'trig', 'periodic'})) )
    % DOMAIN and EPS are specified
    pref.chebfuneps = varargin{3};
    dom = varargin{1};
elseif ( (nargin == 4) && strcmpi(varargin{1}, 'eps') && ...
        ~any(strcmpi(varargin{3}, {'trig', 'periodic'})) )
    % EPS and DOMAIN are specified
    pref.chebfuneps = varargin{2};
    dom = varargin{3};
elseif ( (nargin == 5) && any(strcmpi(varargin{2}, {'trig', 'periodic'}))...
        &&  strcmpi(varargin{3}, 'eps') )
    % DOMAIN, TECH and EPS are specified
    dom = varargin{1};
    pref.tech = @trigtech;
    pref.chebfuneps = varargin{4};
elseif ( (nargin == 5) && strcmpi(varargin{2}, 'eps') ...
        && any(strcmpi(varargin{4}, {'trig', 'periodic'})) )
    % DOMAIN, EPS and TECH are specified
    dom = varargin{1};
    pref.chebfuneps = varargin{3};
    pref.tech = @trigtech;
elseif ( (nargin == 5) && any(strcmpi(varargin{1}, {'trig', 'periodic'}))...
        &&  strcmpi(varargin{2}, 'eps') )
    % TECH, EPS and DOMAIN are specified
    pref.tech = @trigtech;
    pref.chebfuneps = varargin{3};
    dom = varargin{4};
elseif ( (nargin == 5) && any(strcmpi(varargin{1}, {'trig', 'periodic'}))...
        &&  strcmpi(varargin{3}, 'eps') )
    % TECH, DOMAIN and EPS are specified
    pref.tech = @trigtech;
    dom = varargin{2};
    pref.chebfuneps = varargin{4};
elseif ( (nargin == 5) && any(strcmpi(varargin{4}, {'trig', 'periodic'}))...
        &&  strcmpi(varargin{1}, 'eps') )
    % EPS, DOMAIN and TECH are specified
    pref.chebfuneps = varargin{2};
    dom = varargin{3};
    pref.tech = @trigtech;
elseif ( (nargin == 5) && any(strcmpi(varargin{3}, {'trig', 'periodic'}))...
        &&  strcmpi(varargin{1}, 'eps') )
    % EPS, TECH and DOMAIN are specified
    pref.chebfuneps = varargin{2};
    pref.tech = @trigtech;
    dom = varargin{4};
else
    error('CHEBFUN:CHEBFUN3T:constructor', ...
                    'Domain not valid or fully determined.');
end
end
%%
function op = str2op( op )
% OP = STR2OP(OP), finds independent variables in a string and returns an op
% handle than can be evaluated.

vars = symvar(op); % Independent variables
if ( numel(vars) > 3)
    error('CHEBFUN:CHEBFUN3T:constructor:str2op:depvars', ...
        'Too many independent variables in string input.');
else
    op = eval(['@(' vars{1} ',' vars{2} ',' vars{3} ')' op]);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Jacobian of the state transition function for the UNGM-model.
%
% Copyright (C) 2007 Jouni Hartikainen
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function df = ungm_df_dx(x,param)
df = 0.5+25*(1-x.^2)./((1+x.^2).^2);

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
fprintf('Spectral Clustering: Adding directories to current path...\n');

% make sure we add the correct folders even if this file is
% not called from the current folder
fileName = mfilename();
filePath = mfilename('fullpath');
filePath = filePath(1:end-size(fileName, 2));

% Add folders to current path
path(genpath([filePath 'files']), path);

fprintf('Spectral Clustering: Starting GUI...\n');

% Open main GUI
guiMain;

fprintf('Spectral Clustering: Ready.\n\n');

% clear variables
clearvars fileName filePath

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function bw = audfiltbw(fc,varargin)
%AUDFILTBW  Bandwidth of auditory filter
%   Usage: bw = audfiltbw(fc)
%
%   `audfiltbw(fc)` returns the critical bandwidth of the auditory filter
%   at center frequency *fc* defined in equivalent rectangular bandwidth.
%   The function uses the relation:
%
%   .. bw = 24.7 + fc/9.265
%
%   .. math::  bw = 24.7 + \frac{fc}{9.265}
%
%   as estimated in Glasberg and Moore (1990). This function is also used
%   when the original ERB scale ('erb83') is chosen.
%
%   `audfiltbw(fc,'bark')` returns the critical bandwidth at *fc* according
%   to the Bark scale using the relation:
%
%   .. bw = 25 + 75 ( 1+1.4*10^{-6} fc^2 )^0.69
%
%   .. math::  bw = 25 + 75 ( 1+1.4\times10^{-6} fc^2 )^{0.69}
%
%   as estimated by Zwicker and Terhardt (1980).
%
%   For the scales 'mel', 'mel1000', 'log10' and 'semitone', no critical
%   bandwidth function is usually given. Following the example of the
%   equivalent rectangular bandwidth (associated with the ERB scale), we
%   use the derivative of the inverse of the scale function F_{Scale},
%   evaluated at F_{Scale}(fc), i.e.
%
%   .. bw = (F_{scale}^{-1})'(F_{scale}(fc))
%
%   .. math::  bw = (F_{scale}^{-1})'(F_{scale}(fc))
%
%   See also: freqtoerb, erbspace, freqtoaud, audtofreq
%
%   References: glasberg1990daf zwickerterhardt80

%   AUTHOR : Peter L. Søndergaard

% ------ Checking of input parameters ---------

% narginchk(1,1);
if nargin<1
  error('%s: Too few input parameters.',upper(mfilename));
end;

if ~isnumeric(fc) || any(fc(:)<0)
  error('AUDFILTBW: fc must be non-negative.');
end;

definput.flags.audscale={'erb','erb83','bark','mel','mel1000','log10','semitone'};
[flags,kv]=ltfatarghelper({},definput,varargin);

% ------ Computation --------------------------

% FIXME: What is the upper frequency for which the estimation is valid?
if flags.do_erb || flags.do_erb83
    bw = 24.7 + fc/9.265;
end

if flags.do_bark
    bw = 25 + 75*(1 + 1.4E-6*fc.^2).^0.69;
end

if flags.do_mel
    bw = log(17/7)*(700+fc)/1000;
end

if flags.do_mel1000
    bw = log(2)*(1+fc/1000);
end;

if flags.do_log10 || flags.do_semitone
    bw = fc;
end


text:
function xval = r8vec_even_select ( n, xlo, xhi, ival )

%*****************************************************************************80
%
%% R8VEC_EVEN_SELECT returns the I-th of N evenly spaced values in [ XLO, XHI ].
%
%  Discussion:
%
%    XVAL = ( (N-IVAL) * XLO + (IVAL-1) * XHI ) / dble ( N - 1 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of values.
%
%    Input, real XLO, XHI, the low and high values.
%
%    Input, integer IVAL, the index of the desired point.
%    IVAL is normally between 1 and N, but may be any
%    integer value.
%
%    Output, real XVAL, the IVAL-th of N evenly spaced values
%    between XLO and XHI.
%
%    Unless N = 1, X(1) = XLO and X(N) = XHI.
%
%    If N = 1, then X(1) = 0.5*(XLO+XHI).
%
  if ( n == 1 )

    xval = 0.5 * ( xlo + xhi );

  else

    xval = ( ( n - ival     ) * xlo   ...
           + (     ival - 1 ) * xhi ) ...
           / ( n        - 1 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [g1, g2] = disimXdisimKernGradient(disimKern1, disimKern2, t1, t2, covGrad)

% DISIMXDISIMKERNGRADIENT Compute a cross gradient between two DISIM kernels.
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel
% between two disim kernels for the multiple output kernel.
% ARG disimKern1 : the kernel structure associated with the first DISIM
% kernel.
% ARG disimKern2 : the kernel structure associated with the second DISIM
% kernel.
% ARG t : inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see disimKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see disimKernExtractParam.
%
% FORMAT
% DESC computes cross kernel terms between two DISIM kernels for
% the multiple output kernel.
% ARG disimKern1 : the kernel structure associated with the first DISIM
% kernel.
% ARG disimKern2 : the kernel structure associated with the second DISIM
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see disimKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see disimKernExtractParam.
%
% SEEALSO : multiKernParamInit, multiKernCompute, disimKernParamInit, disimKernExtractParam
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Antti Honkela, 2007-2009

% KERN

arg{1}=t1;
if nargin < 5
  covGrad = t2;
  t2 = t1;
else
  arg{2}=t2;
end
if size(t1, 2) > 1 || size(t2, 2) > 1
  error('Input can only have one column');
end
if disimKern1.inverseWidth ~= disimKern2.inverseWidth
  error('Kernels cannot be cross combined if they have different inverse widths.')
end
if disimKern1.di_decay ~= disimKern2.di_decay
  error('Kernels cannot be cross combined if they have different driving input decays.');
end
if disimKern1.di_variance ~= disimKern2.di_variance
  error('Kernels cannot be cross combined if they have different driving input variances.');
end
if disimKern1.rbf_variance ~= disimKern2.rbf_variance
  error('Kernels cannot be cross combined if they have different RBF variances.');
end

delta = disimKern1.di_decay;
D1 = disimKern1.decay;
D2 = disimKern2.decay;

l = sqrt(2/disimKern1.inverseWidth);
[h1, dh1_ddelta, dh1_dD1, dh1_dD2, dh1_dl] = disimComputeH(t1, t2, delta, D1, D2, l);
[hp1, dhp1_ddelta, dhp1_dD1, dhp1_dD2, dhp1_dl] = disimComputeHPrime(t1, t2, delta, D1, D2, l);

% Avoid making the expensive call twice with the same arguments
if ((length(t1) == length(t2)) && all(t1 == t2) && ...
    (D1 == D2)),
  h2 = h1;
  dh2_ddelta = dh1_ddelta;
  dh2_dD2 = dh1_dD1;
  dh2_dD1 = dh1_dD2;
  dh2_dl = dh1_dl;

  hp2 = hp1;
  dhp2_ddelta = dhp1_ddelta;
  dhp2_dD2 = dhp1_dD1;
  dhp2_dD1 = dhp1_dD2;
  dhp2_dl = dhp1_dl;
else,
  [h2, dh2_ddelta, dh2_dD2, dh2_dD1, dh2_dl] = disimComputeH(t2, t1, delta, D2, D1, l);
  [hp2, dhp2_ddelta, dhp2_dD2, dhp2_dD1, dhp2_dl] = disimComputeHPrime(t2, t1, delta, D2, D1, l);
end

dK_ddelta = dh1_ddelta + dh2_ddelta' + dhp1_ddelta + dhp2_ddelta';
dK_dD1 = dh1_dD1 + dh2_dD1' + dhp1_dD1 + dhp2_dD1';
dK_dD2 = dh1_dD2 + dh2_dD2' + dhp1_dD2 + dhp2_dD2';
dK_dl = dh1_dl + dh2_dl' + dhp1_dl + dhp2_dl';

C0 = disimKern1.di_variance;
C1 = sqrt(disimKern1.variance);
C2 = sqrt(disimKern2.variance);
C3 = disimKern1.rbf_variance;
K = h1 + h2' + hp1 + hp2';
K = 0.5*K*sqrt(pi);
var2 = C0*C1*C2*C3;
dk_ddelta = (sum(sum(covGrad.*dK_ddelta)))*0.5*sqrt(pi)*l*var2;
dk_dD1 = (sum(sum(covGrad.*dK_dD1)))*0.5*sqrt(pi)*l*var2;
dk_dD2 = (sum(sum(covGrad.*dK_dD2)))*0.5*sqrt(pi)*l*var2;
dk_dl = sum(sum(covGrad.*(dK_dl*0.5*sqrt(pi)*l + K)))*var2;
K = l*K;
dk_dC0 = C1*C2*C3*sum(sum(covGrad.*K));
dk_dC1 = C0*C2*C3*sum(sum(covGrad.*K));
dk_dC2 = C0*C1*C3*sum(sum(covGrad.*K));
dk_dC3 = C0*C1*C2*sum(sum(covGrad.*K));

dk_dDIVariance = dk_dC0;
dk_dDisim1Variance = dk_dC1*0.5/C1;
dk_dDisim2Variance = dk_dC2*0.5/C2;
dk_dRBFVariance = dk_dC3;

dk_dinvWidth = -0.5*sqrt(2)/(disimKern1.inverseWidth* ...
                             sqrt(disimKern1.inverseWidth))*dk_dl;

K = var2*K;

if isfield(disimKern1, 'gaussianInitial') && disimKern1.gaussianInitial && ...
  isfield(disimKern2, 'gaussianInitial') && disimKern2.gaussianInitial,
  if disimKern1.initialVariance ~= disimKern2.initialVariance
    error('Kernels cannot be cross combined if they have different initial variances.');
  end

  dim1 = size(t1, 1);
  dim2 = size(t2, 1);
  t1Mat = t1(:, ones(1, dim2));
  t2Mat = t2(:, ones(1, dim1))';

  the_rest = (exp(- delta * t1Mat) - exp(- D1 * t1Mat)) ./ (D1 - delta) .* ...
      (exp(- delta * t2Mat) - exp(- D2 * t2Mat)) ./ (D2 - delta);

  dk_dinitVariance = sum(sum((sqrt(disimKern1.variance) * ...
			      sqrt(disimKern2.variance) * ...
			      the_rest) .* covGrad));
  dk_dDisim1Variance = dk_dDisim1Variance + ...
      sum(sum((.5 ./ sqrt(disimKern1.variance) * ...
	       disimKern1.initialVariance * sqrt(disimKern2.variance) * the_rest) .* covGrad));
  dk_dDisim2Variance = dk_dDisim2Variance + ...
      sum(sum((.5 ./ sqrt(disimKern2.variance) * ...
	       disimKern1.initialVariance * sqrt(disimKern1.variance) * the_rest) .* covGrad));

  dk_dD1 = dk_dD1 + ...
	   sum(sum((disimKern1.initialVariance * ...
		    sqrt(disimKern1.variance) * sqrt(disimKern2.variance) * ...
		    (t1Mat * (D1 - delta).*exp(-D1*t1Mat) - exp(-delta*t1Mat) + exp(-D1*t1Mat)) ./ (D1-delta).^2 .* ...
		    (exp(- delta * t2Mat) - exp(- D2 * t2Mat)) ./ (D2 - delta)).*covGrad));

  dk_dD2 = dk_dD2 + ...
	   sum(sum((disimKern1.initialVariance * ...
		    sqrt(disimKern1.variance) * sqrt(disimKern2.variance) * ...
		    (t2Mat * (D2 - delta).*exp(-D2*t2Mat) - exp(-delta*t2Mat) + exp(-D2*t2Mat)) ./ (D2-delta).^2 .* ...
		    (exp(- delta * t1Mat) - exp(-D1 * t1Mat)) ./ (D1 - delta)).*covGrad));

  dk_ddelta = dk_ddelta + ...
      sum(sum((disimKern1.initialVariance * ...
	       sqrt(disimKern1.variance) * sqrt(disimKern2.variance) * ...
	       ((-t2Mat * (D2 - delta).*exp(-delta*t2Mat) + exp(-delta*t2Mat) - exp(-D2*t2Mat)) ./ (D2-delta).^2 .* ...
	       (exp(- delta * t1Mat) - exp(-D1 * t1Mat)) ./ (D1 - delta) + ...
		(-t1Mat * (D1 - delta).*exp(-delta*t1Mat) + exp(-delta*t1Mat) - exp(-D1*t1Mat)) ./ (D1-delta).^2 .* ...
		(exp(- delta * t2Mat) - exp(-D2 * t2Mat)) ./ (D2 - delta))).*covGrad));

  g1 = [dk_ddelta dk_dinvWidth dk_dDIVariance dk_dD1 dk_dDisim1Variance dk_dRBFVariance dk_dinitVariance];
  g2 = [0 0 0 dk_dD2 dk_dDisim2Variance 0 0];
else
  % only pass the gradient with respect to the inverse width to one
  % of the gradient vectors ... otherwise it is counted twice.
  g1 = [dk_ddelta dk_dinvWidth dk_dDIVariance dk_dD1 dk_dDisim1Variance dk_dRBFVariance];
  g2 = [0 0 0 dk_dD2 dk_dDisim2Variance 0];
end


text:
%   AUTHOR:
%       Boguslaw Obara, http://boguslawobara.net/
%% Read image
im = imread('triple.tif');
im = im<10;
%% Skeleton
im = bwmorph(im, 'thin', inf);
%% Hit or Miss
out1 = BOHitOrMiss(im, 'end');
out2 = BOHitOrMiss(im, 'triple');
%% Plot
ims = im;
ims = ims + 2*out1;
ims = ims + 8*out2;
imagesc(ims);
%%


text:
function pass = test_zerothOrder(~)
% TEST_ZEROTHORDER   Solve problems without any differentiation/integration

%% Do a manual Newton iteration to check everything is working at lower levels
x = chebfun('x');
% RHS
f = sin(x)+2;
% Initial guess
u = 0*x;
op = @(u) u.^2 + sin(u) + exp(u) - f;
p = cheboppref();
p.discretization = @chebcolloc2;
normdel = inf;
while normdel > 1e-10
    uad = adchebfun(u);
    res = op(uad);
    del = linsolve(linop(res.jacobian), res.func, p);
    normdel = norm(del);
    u = chebfun(u);
    u = u - del;
end

pass(1) = normdel < 1e-10;

%% Iteration with CHEBOPs
x = chebfun('x');
f = sin(4*x) + 2;
N = chebop(@(u) 4*u.^2 - x+.1+sin(u));
u = N\f;
pass(2) = norm(N(u) - f) < 1e-10;

%% Coupled systems
x = chebfun('x');
f = sin(4*x) + 2;
N = chebop(@(x,u,v) [u.^2 - x+.1+sin(v); v + exp(v) + u]);
[u, v] = N\[f; f];
pass(3) = norm(N(u,v) - [f;f]) < 1e-10;

end


text:
function [varargout]=gpatch(varargin)

% function [h]=gpatch(F,V,C,CE,A,L)
% ------------------------------------------------------------------------
% This function is a short-hand version of the |patch| command. The inputs
% for |gpatch| are the faces (F), the vertices (V), the color description
% (C), the edge color description CE, the transparancy (A), and the edge
% width (L).
% The color data descriptions C (or equivalently CE for edges) can be:
% 1) A string such as 'g' for green
% 2) A triplet of RGD values e.g. [1 0 0] is blue
% 3) A nx1 or a mx1 array of colormapped colors (where n=size(F,1) or
% m=size(V,1))
% 4) (simiarl to 3) A nx3 or a mx3 RGB color value array for the faces or
% vertices respectively.
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2017
% 2018/02/07 Added support for colormapped edges
%------------------------------------------------------------------------

switch nargin
    case 1
        error('Not enough input arguments, provide at least faces and vertices');
    case 2
        F=varargin{1};
        V=varargin{2};
        C='g';
        CE='k';
        A=1;
        L=[];
    case 3
        F=varargin{1};
        V=varargin{2};
        C=varargin{3};
        CE='k';
        A=1;
        L=[];
    case 4
        F=varargin{1};
        V=varargin{2};
        C=varargin{3};
        CE=varargin{4};
        A=1;
        L=[];
    case 5
        F=varargin{1};
        V=varargin{2};
        C=varargin{3};
        CE=varargin{4};
        A=varargin{5};
        L=[];
    case 6
        F=varargin{1};
        V=varargin{2};
        C=varargin{3};
        CE=varargin{4};
        A=varargin{5};
        L=varargin{6};
    otherwise
        error('Wrong number of input arguments');
end

if isempty(C)
    C='g';
end

if isempty(CE)
    C='k';
end

if isa(F,'cell') %Assume all entries are cells defining multiple patch data sets
    for q=1:1:numel(F)
        f=F{q};

        if isa(V,'cell')
            v=V{q};
        else
            v=V;
        end

        if isa(C,'cell')
            c=C{q};
        else
            c=C;
        end

        if isa(CE,'cell')
            ce=CE{q};
        else
            ce=CE;
        end

        if isa(A,'cell')
            a=A{q};
        else
            a=A;
        end

        hp(q)=plotPatch(f,v,c,ce,a,L);
    end
else
    hp=plotPatch(F,V,C,CE,A,L);
end

if nargout==1
    varargout{1}=hp;
end

end

%%
function hp=plotPatch(F,V,C,CE,A,L)

% hf=gcf;
% if isempty(hf.Children)
%     gca;
%     view(3);
% end

argInPatch.Faces=F;
argInPatch.Vertices=V;
argInPatch.EdgeColor=CE;

if ischar(C) %Plain single color
    argInPatch.FaceColor=C;
    if strcmp(C,'kw')
        argInPatch.FaceColor=grayColor(0.5);
    end
    if strcmp(C,'rw')
        argInPatch.FaceColor=[1 0.5 0.5];
    end
    if strcmp(C,'gw')
        argInPatch.FaceColor=[0.5 1 0.5];
    end
    if strcmp(C,'bw')
        argInPatch.FaceColor=[0.5 0.5 1];
    end
    if strcmp(C,'yw')
        argInPatch.FaceColor=[1 1 0.5];
    end
    if strcmp(C,'cw')
        argInPatch.FaceColor=[0.5 1 1];
    end
    if strcmp(C,'mw')
        argInPatch.FaceColor=[1 0.5 1];
    end
elseif size(C,2)==1
    argInPatch.FaceColor='flat';
    argInPatch.CData=double(C);
elseif size(C,2)==3 && size(C,1)==1 %Assume single RGB level
    argInPatch.FaceColor=double(C);
elseif size(C,2)==3 && size(C,1)>1 %Assume RGB array
    argInPatch.FaceColor='flat';
    argInPatch.FaceVertexCData=double(C);
else
    error('Invalid face-vertex color data input');
end

if ischar(CE) %Plain single color
    argInPatch.EdgeColor=CE;
    if strcmp(CE,'kw')
        argInPatch.EdgeColor=grayColor(0.5);
    end
    if strcmp(CE,'rw')
        argInPatch.EdgeColor=[1 0.5 0.5];
    end
    if strcmp(CE,'gw')
        argInPatch.EdgeColor=[0.5 1 0.5];
    end
    if strcmp(CE,'bw')
        argInPatch.EdgeColor=[0.5 0.5 1];
    end
    if strcmp(CE,'yw')
        argInPatch.EdgeColor=[1 1 0.5];
    end
    if strcmp(CE,'cw')
        argInPatch.EdgeColor=[0.5 1 1];
    end
    if strcmp(CE,'mw')
        argInPatch.EdgeColor=[1 0.5 1];
    end
elseif size(CE,2)==1
    if size(CE,1)>1
        if size(CE,1)==size(F,1)
            [CE]=faceToVertexMeasure(F,V,CE);
            argInPatch.EdgeColor='flat';
            argInPatch.CData=double(CE);
        end
        if size(CE,1)==size(V,1)
            argInPatch.EdgeColor='flat';
            argInPatch.CData=double(CE);
        end
    else
        argInPatch.EdgeColor='flat';
        argInPatch.CData=double(CE)*ones(size(V,1),1);
    end
elseif size(CE,2)==3 && size(CE,1)==1 %Assume single RGB level
    argInPatch.EdgeColor=double(CE);
elseif size(CE,2)==3 && size(CE,1)>1 %Assume RGB array
    argInPatch.EdgeColor='flat';
    argInPatch.FaceVertexCData=double(CE);
else
    error('Invalid edge color data input');
end

if numel(A)==1 %Plain single alpha
    argInPatch.FaceAlpha=A;
elseif size(A,2)==1 %Alpha mapping
    argInPatch.FaceAlpha='flat';
    argInPatch.FaceVertexAlphaData=A;
else
    error('Invalid alpha data input');
end

if ~isempty(L)
    argInPatch.LineWidth=L;
end

hp=patch(argInPatch);

end

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

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function string=implode(pieces,delimiter)
% IMPLODE    Joins strings with delimiter in between.
%
% IMPLODE(PIECES,DELIMITER) returns a string containing all the
% strings in PIECES joined with the DELIMITER string in between.
%
% :Inputs:
%
%   **PIECES:**
%        the pieces of string to join (cell array), each cell is a piece
%
%   **DELIMITER:**
%        the delimiter string to put between the pieces (string)
%
% :Output:
%
%   **STRING:**
%        all the pieces joined with the delimiter in between (string)
%
% :Examples:
% ::
%
%    PIECES = {'ab','c','d','e fgh'}
%    DELIMITER = '->'
%    STRING = IMPLODE(PIECES,DELIMITER)
%    STRING = ab->c->d->e fgh
%
% :See also:
% EXPLODE, STRCAT
%
% ..
%   Created: Sara Silva (sara@itqb.unl.pt) - 2002.08.25
% ..

if isempty(pieces) % no pieces to join, return empty string
   string='';

else % no need for delimiters yet, so far there's only one piece
   string=pieces{1};
end

l=length(pieces);
p=1;
while p<l % more than one piece to join with the delimiter, the interesting case
   p=p+1;
   string=[string delimiter pieces{p}];
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%%**********************************************************
%% SDPT3data_SEDUMIdata: convert SQLP data in SDPT3 format to
%%                       SeDuMi format
%%
%% [At,b,c,K] = SDPT3data_SEDUMIdata(blk,AAt,CC,bb);
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%**********************************************************

function [At,b,c,K] = SDPT3data_SEDUMIdata(blk,AAt,CC,bb)

c = []; At = [];
b = bb;
mm = length(bb);
%%
if (~iscell(CC))
    Ctmp = CC; clear CC; CC{1} = Ctmp;
end
%%
%% extract unrestricted blk
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if (p==1); K.f = []; end
    if strcmp(pblk{1},'u')
        K.f = [K.f, pblk{2}];
        At  = [At; AAt{p}]; %#ok
        c   = [c; CC{p}]; %#ok
    end
end
K.f = sum(K.f);
%%
%% extract linear blk
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if (p==1); K.l = []; end
    if strcmp(pblk{1},'l')
        K.l = [K.l, pblk{2}];
        At  = [At; AAt{p,1}]; %#ok
        c   = [c; CC{p,1}]; %#ok
    end
end
K.l = sum(K.l);
%%
%% extract second order cone blk
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if (p==1); K.q = []; end
    if strcmp(pblk{1},'q')
        K.q = [K.q, pblk{2}];
        At  = [At; AAt{p,1}]; %#ok
        c   = [c; CC{p,1}]; %#ok
    end
end
%%
%% extract rotated cone blk
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if (p==1); K.r = []; end
    if strcmp(pblk{1},'r')
        K.r = [K.r, pblk{2}];
        At  = [At; AAt{p,1}]; %#ok
        c   = [c; CC{p,1}]; %#ok
    end
end
%%
%% extract semidefinite cone blk
%%
for p = 1:size(blk,1)
    if (p==1); K.s = []; end
    pblk = blk(p,:);
    if strcmp(pblk{1},'s')
        K.s = [K.s, pblk{2}];
        ss = [0,cumsum(pblk{2})];
        idxstart = [0,cumsum(pblk{2}.*pblk{2})];
        numblk = length(pblk{2});
        nnzA = nnz(AAt{p,1});
        II = zeros(2*nnzA,1);
        JJ = zeros(2*nnzA,1);
        VV = zeros(2*nnzA,1);
        m2 = size(AAt{p,1},2);
        if (length(pblk) > 2)
            rr = [0, cumsum(pblk{3})];
            dd = AAt{p,3};
            idxD = [0; find(diff(dd(:,1))); size(dd,1)];
        end
        count = 0;
        for k = 1:mm
            if (k<= m2);
                Ak = smat(pblk,AAt{p,1}(:,k),1);
            else
                idx = rr(k)+1 : rr(k+1);
                Vk  = AAt{p,2}(:,idx);
                len = pblk{3}(k);
                if (size(dd,2) == 4)
                    idx2 = idxD(k)+1:idxD(k+1);
                    Dk = spconvert([dd(idx2,2:4); len,len,0]);
                elseif (size(dd,2) == 1);
                    Dk = spdiags(dd(idx),0,len,len);
                end
                Ak = Vk*Dk*Vk';
            end
            for tt = 1:numblk
                if (numblk > 1)
                    idx = ss(tt)+1: ss(tt+1);
                    Aksub = full(Ak(idx,idx));
                else
                    Aksub = Ak;
                end
                tmp = Aksub(:);
                nzidx = find(tmp);
                len = length(nzidx);
                II(count+1:count+len,1) = idxstart(tt)+nzidx;
                JJ(count+1:count+len,1) = k*ones(length(nzidx),1);
                VV(count+1:count+len,1) = tmp(nzidx);
                count = count + len;
            end
        end
        II = II(1:count);
        JJ = JJ(1:count);
        VV = VV(1:count);
        At = [At; spconvert([II,JJ,VV; sum(pblk{2}.*pblk{2}), mm, 0])]; %#ok
        Cp = CC{p};
        ctmp = [];
        for tt = 1:numblk
            if (numblk > 1)
                idx = ss(tt)+1: ss(tt+1);
                Csub = full(Cp(idx,idx));
            else
                Csub = Cp;
            end
            ctmp = [ctmp; Csub(:)]; %#ok
        end
        c = [c; ctmp]; %#ok
    end
end
%%**********************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function r8_epsilon_test ( )

%*****************************************************************************80
%
%% R8_EPSILON_TEST tests R8_EPSILON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_EPSILON_TEST\n' );
  fprintf ( 1, '  R8_EPSILON produces the R8 roundoff unit.\n' );
  fprintf ( 1, '\n' );

  r = r8_epsilon ( );
  fprintf ( 1, '  R = R8_EPSILON()         = %e\n', r );

  s = ( 1.0 + r ) - 1.0;
  fprintf ( 1, '  ( 1 + R ) - 1            = %e\n', s );

  s = ( 1.0 + ( r / 2.0 ) ) - 1.0;
  fprintf ( 1, '  ( 1 + (R/2) ) - 1        = %e\n', s );

  return
end


text:
function result = cartesianProduct(sets)
    c = cell(1, numel(sets));
    [c{:}] = ndgrid( sets{:} );
    result = cell2mat( cellfun(@(v)v(:), c, 'UniformOutput',false) );
end

author: Arrowstar


text:
classdef Filter < dagnn.Layer
  properties
    pad = [0 0 0 0]
    stride = [1 1]
    dilate = [1 1]
  end
  methods
    function set.pad(obj, pad)
      if numel(pad) == 1
        obj.pad = [pad pad pad pad] ;
      elseif numel(pad) == 2
        obj.pad = pad([1 1 2 2]) ;
      else
        obj.pad = pad ;
      end
    end

    function set.stride(obj, stride)
      if numel(stride) == 1
        obj.stride = [stride stride] ;
      else
        obj.stride = stride ;
      end
    end

    function set.dilate(obj, dilate)
      if numel(dilate) == 1
        obj.dilate = [dilate dilate] ;
      else
        obj.dilate = dilate ;
      end
    end

    function kernelSize = getKernelSize(obj)
      kernelSize = [1 1] ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      ks = obj.getKernelSize() ;
      ke = (ks - 1) .* obj.dilate + 1 ;
      outputSizes{1} = [...
        fix((inputSizes{1}(1) + obj.pad(1) + obj.pad(2) - ke(1)) / obj.stride(1)) + 1, ...
        fix((inputSizes{1}(2) + obj.pad(3) + obj.pad(4) - ke(2)) / obj.stride(2)) + 1, ...
        1, ...
        inputSizes{1}(4)] ;
    end

    function rfs = getReceptiveFields(obj)
      ks = obj.getKernelSize() ;
      ke = (ks - 1) .* obj.dilate + 1 ;
      y1 = 1 - obj.pad(1) ;
      y2 = 1 - obj.pad(1) + ke(1) - 1 ;
      x1 = 1 - obj.pad(3) ;
      x2 = 1 - obj.pad(3) + ke(2) - 1 ;
      h = y2 - y1 + 1 ;
      w = x2 - x1 + 1 ;
      rfs.size = [h, w] ;
      rfs.stride = obj.stride ;
      rfs.offset = [y1+y2, x1+x2]/2 ;
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function auto_cell = automobile_read ( input_filename )

%*****************************************************************************80
%
%% AUTOMOBILE_READ reads data from the automobile dataset.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    A Frank, A Asuncion,
%    UCI Machine Learning Repository,
%    http://archive.ics.uci.edu/ml,
%    School of Information and Computer Science,
%    University of California, Irvine, California.
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, cell array AUTO_CELL(205,26), the information from the file.
%
  auto_cell = cell ( 205, 26 );
%
%  Open the file.
%
  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'AUTOMOBILE_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'AUTOMOBILE_READ - Error!' );
  end
%
%  Read the records.
%
  i = 0;

  while ( 1 )

    line = fgets ( input_unit );

    if ( line == -1 )
      break
    end
%
%  Read the I-th record.
%
    i = i + 1;
%
%  Replace commas by spaces, so the line consists of 26 "words".
%
    line = s_replace_ch ( line, ',', ' ' );
%
%  Replace ? by -1, so missing values are marked.
%
    line = s_replace_ch_by_s ( line, '?', '-1' );
%
%  Read the 26 words, using an integer format for the first word,
%  a real format for others, and a string format for others.
%
    done = 1;

    for j = 1 : 26
      [ w, done ] = word_next_read ( line, done );
      if ( j == 1 )
        auto_cell{i,j} = sscanf ( w, '%d' );
      elseif ( j == 2 || j == 10 || j == 12 || j == 13 || j == 14 || j == 17 || ...
               j == 19 || j == 20 || j == 21 || j == 22 || j == 23 || j == 24 || ...
               j == 25 || j == 26 )
        auto_cell{i,j} = sscanf ( w, '%f' );
      else
        auto_cell{i,j} = w;
      end
    end

  end

  fclose ( input_unit );

  if ( i ~= 205 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'AUTOMOBILE_READ - Warning!\n' );
    fprintf ( 1, '  The data file should contain 205 records.\n' );
    fprintf ( 1, '  This program found %d records.\n', i );
  end

  return
end
function s = s_replace_ch ( s, c1, c2 )

%*****************************************************************************80
%
%% S_REPLACE_CH replaces all occurrences of one character by another.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, character S(*), the string.
%
%    Input, character C1, C2, the character to be replaced, and the
%    replacement character.
%
  for i = 1 : length ( s )
    if ( s(i) == c1 )
      s(i) = c2;
    end
  end

  return
end
function s_new = s_replace_ch_by_s ( s_old, c, s_rep )

%*****************************************************************************80
%
%% S_REPLACE_CH_BY_S replaces all occurrences of a character by a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S_OLD, the string to be modified.
%
%    Input, character C, the character to be replaced.
%
%    Input, string S_REP, the replacement string.
%
%    Output, string S_NEW, the modified string.
%
  n_old = length ( s_old );
  n_rep = length ( s_rep );

  i_new = 1;

  for i_old = 1 : n_old

    if ( s_old(i_old) ~= c )
      s_new(i_new) = s_old(i_old);
      i_new = i_new + 1;
    else
      s_new(i_new:i_new+n_rep-1) = s_rep;
      i_new = i_new + n_rep;
    end

  end

  return
end
function [ word, done ] = word_next_read ( s, done )

%*****************************************************************************80
%
%% WORD_NEXT_READ "reads" words from a string, one at a time.
%
%  Special cases:
%
%    The following characters are considered to be a single word,
%    whether surrounded by spaces or not:
%
%      " ( ) { } [ ]
%
%    Also, if there is a trailing comma on the word, it is stripped off.
%    This is to facilitate the reading of lists.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 December 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, a string, presumably containing words
%    separated by spaces.
%
%    Input, logical DONE.
%    TRUE, if we are beginning a new string;
%    FALSE, if we are continuing to process the current string.
%
%    Output, string WORD.
%    If DONE is FALSE, then WORD contains the "next" word read.
%    If DONE is TRUE, then WORD is blank, because there was no more to read.
%
%    Output, logical DONE.
%      FALSE if another word was read,
%      TRUE if no more words could be read.
%
  persistent lenc;
  persistent next;

  tab = char ( 9 );
%
%  We "remember" LENC and NEXT from the previous call.
%
%  An input value of DONE = TRUE signals a new line of text to examine.
%
  if ( done )

    next = 1;
    done = 0;
    lenc = length ( s );

    if ( lenc <= 0 )
      done = 1;
      word = ' ';
      return
    end

  end
%
%  Beginning at index NEXT, search the string for the next nonblank,
%  which signals the beginning of a word.
%
  ilo = next;
%
%  ...S(NEXT:) is blank.  Return with WORD = ' ' and DONE = TRUE.
%
  while ( 1 )

    if ( lenc < ilo )
      word = ' ';
      done = 1;
      next = lenc + 1;
      return
    end
%
%  If the current character is blank, skip to the next one.
%
    if ( s(ilo) ~= ' ' && s(ilo) ~= tab )
      break
    end

    ilo = ilo + 1;

  end
%
%  ILO is the index of the next nonblank character in the string.
%
%  If this initial nonblank is a special character,
%  then that's the whole word as far as we're concerned,
%  so return immediately.
%
  if ( s(ilo) == '"' || ...
       s(ilo) == '(' || ...
       s(ilo) == ')' || ...
       s(ilo) == '{' || ...
       s(ilo) == '}' || ...
       s(ilo) == '[' || ...
       s(ilo) == ']' )

    word = s(ilo);
    next = ilo + 1;
    return

  end
%
%  Now search for the last contiguous character that is not a
%  blank, TAB, or special character.
%
  next = ilo + 1;

  while ( next <= lenc )

    if ( s(next) == ' ' )
      break;
    elseif ( s(next) == tab )
      break;
    elseif ( s(next) == '"' )
      break;
    elseif ( s(next) == '(' )
      break;
    elseif ( s(next) == ')' )
      break;
    elseif ( s(next) == '{' )
      break;
    elseif ( s(next) == '}' )
      break;
    elseif ( s(next) == '[' )
      break;
    elseif ( s(next) == ']' )
      break;
    end

    next = next + 1;

  end

  if ( s(next-1) == ',' )
    word = s(ilo:next-2);
  else
    word = s(ilo:next-1);
  end

  return
end