text:
classdef TREE2 < PROBLEM
% <multi> <real> <large/none> <constrained> <expensive/none>
% The time-varying ratio error estimation problem
% T --- 1000 --- Length of data (related to the number of variables)

%------------------------------- Reference --------------------------------
% C. He, R. Cheng, C. Zhang, Y. Tian, Q. Chen, and X. Yao, Evolutionary
% large-scale multiobjective optimization for ratio error estimation of
% voltage transformers, IEEE Transactions on Evolutionary Computation,
% 2020, 24(5): 868-881.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties(Access = private)
        Data;   % Dataset
        Mean;   % Mean values of the dataset
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            % Load data
            T = obj.ParameterSet(1000);
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'Dataset_TREE.mat'),'Dataset');
            obj.Data = Dataset.TREE2(1:min(max(T,10),end),:);
            % Set the numbers of objectives and decision variables
            K = 3;
            obj.M = 2;
            obj.D = T*K;
            % Set the upper and lower boundaries
            Lower    = zeros(T,K);
            Upper    = zeros(T,K);
            obj.Mean = zeros(T,K);
            for i = 1 : 3
                Lower(:,i)    = min(obj.Data(:,i:3:end),[],2)/1.01;
                Upper(:,i)    = max(obj.Data(:,i:3:end),[],2)/0.99;
                obj.Mean(:,i) = mean(obj.Data(:,i:3:end),2);
            end
            % The decision variables are the offset of the mean values of
            % the dataset
            obj.lower    = reshape(Lower-obj.Mean,1,[]);
            obj.upper    = reshape(Upper-obj.Mean,1,[]);
            obj.encoding = ones(1,obj.D);
        end
        %% Generate initial solutions
        function Population = Initialization(obj,N)
            if nargin < 2; N = obj.N; end
            PopDec = obj.Mean(:)'.*(rand(N,obj.D)*0.008-0.004);
            PopDec = min(max(PopDec,repmat(obj.lower,N,1)),repmat(obj.upper,N,1));
            Population = obj.Evaluation(PopDec);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            N      = size(PopDec,1);
            KP     = size(obj.Data,2);
            [T,K]  = size(obj.Mean);
            PopDec = reshape(PopDec,N,T,K) + repmat(reshape(obj.Mean,1,T,K),N,1,1);
            PopObj = zeros(N,obj.M);
            % First objective
            eA = abs(repmat(reshape(obj.Data,1,T,KP),N,1,1)./repmat(PopDec(:,:,1:3),1,1,KP/3)-1);
            PopObj(:,1) = sum(sum(eA,3),2);
            % Second objective
            Delta = std(eA(:,2:end,:)-eA(:,1:end-1,:),0,2);
            PopObj(:,2) = sum(reshape(Delta,N,[]),2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            N      = size(PopDec,1);
            [T,K]  = size(obj.Mean);
            PopDec = reshape(PopDec,N,T,K) + repmat(reshape(obj.Mean,1,T,K),N,1,1);
            PopCon = TREE_CalCon(PopDec,1);
        end
        %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            X = zeros(1,obj.D);
            X(1:2:end) = obj.lower(1:2:end);
            X(2:2:end) = obj.upper(2:2:end);
            R = obj.CalObj(X);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [w0] = mci_interp_init (Y,M)
% Linear interpolate to t=0
% FORMAT [w0] = mci_interp_init (Y,M)
%
% Y     Cell of data from multiple subjects
%       Y{n}.y, Y{n}.ind for n=1..N
% M     Model structure
%
% w0    [d x N] matrix of initial states
%       where d is number of states
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_interp_init.m 6548 2015-09-11 12:39:47Z will $

N=length(Y);
d=size(Y{1}.y,2);
doplot=0;

for n=1:N,
    for j=1:d,

        % Fit
        Nt=size(Y{n}.y,1);
        xd=[Y{n}.ind(:),ones(Nt,1)];
        yd=Y{n}.y(:,j);
        beta=pinv(xd)*yd;

        % Extrapolate
        xt=[[1:M.N]',ones(M.N,1)];
        yhat=xt*beta;

        if doplot
            figure;plot(Y{n}.ind,yd,'x');
            hold on; plot([1:M.N]',yhat,'r');
        end
        w0(j,n)=yhat(1);
    end
end

author: spm


text:
function c = getNamespace()
global namespace

% fprintf('In getNamespace:\n');

if length(namespace)==1
    c = namespace(1).name;
    return;
end

c = '';
s = dbstack();
if isempty(s)
    return;
end

for jj = length(s):-1:1
%     fprintf('  s(%d).file = %s:\n', jj, s(jj).file);
%     fprintf('  s(%d).name = %s:\n', jj, s(jj).name);

    try
        pname = filesepStandard(which(s(jj).file));
    catch
        continue
    end
    for ii = 1:length(namespace)
%         fprintf('  namespace = [%s, %s]\n', namespace(ii).name, namespace(ii).pname);
        if isRootpathsSame(pname, namespace(ii).pname)
            c = namespace(ii).name;
            return;
        end
    end
%     fprintf('\n');
end

% ---------------------------------------------------------
function b = isRootpathsSame(p1, p2)
b = false;
p1 = filesepStandard(p1);
p2 = filesepStandard(p2);
k = findstr(p1, p2);
if isempty(k)
    return;
end
b = true;

author: BUNPC


text:
function [ b, det ] = r8mat_inverse_3d ( a )

%*****************************************************************************80
%
%% R8MAT_INVERSE_3D inverts a 3 by 3 R8MAT using Cramer's rule.
%
%  Discussion:
%
%    If DET is zero, then A is singular, and does not have an
%    inverse.  In that case, B is simply set to zero, and a
%    message is printed.
%
%    If DET is nonzero, then its value is roughly an estimate
%    of how nonsingular the matrix A is.
%
%  Modified:
%
%    31 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(3,3), the matrix to be inverted.
%
%    Output, real B(3,3), the inverse of the matrix.
%
%    Output, real DET, the determinant of the matrix.
%

%
%  Compute the determinant of A
%
  det =   a(1,1) * ( a(2,2) * a(3,3) - a(2,3) * a(3,2) ) ...
        + a(1,2) * ( a(2,3) * a(3,1) - a(2,1) * a(3,3) ) ...
        + a(1,3) * ( a(2,1) * a(3,2) - a(2,2) * a(3,1) );
%
%  If the determinant is zero, bail out.
%
  if ( det == 0.0 )
    b(1:3,1:3) = 0.0;
    return
  end
%
%  Compute the entries of the inverse matrix using an explicit
%  formula.
%
  b(1,1) = + ( a(2,2) * a(3,3) - a(2,3) * a(3,2) ) / det;
  b(1,2) = - ( a(1,2) * a(3,3) - a(1,3) * a(3,2) ) / det;
  b(1,3) = + ( a(1,2) * a(2,3) - a(1,3) * a(2,2) ) / det;

  b(2,1) = - ( a(2,1) * a(3,3) - a(2,3) * a(3,1) ) / det;
  b(2,2) = + ( a(1,1) * a(3,3) - a(1,3) * a(3,1) ) / det;
  b(2,3) = - ( a(1,1) * a(2,3) - a(1,3) * a(2,1) ) / det;

  b(3,1) = + ( a(2,1) * a(3,2) - a(2,2) * a(3,1) ) / det;
  b(3,2) = - ( a(1,1) * a(3,2) - a(1,2) * a(3,1) ) / det;
  b(3,3) = + ( a(1,1) * a(2,2) - a(1,2) * a(2,1) ) / det;

author: johannesgerer

repo: jburkardt-m


text:
function out = map_feature(feat1, feat2)
% MAP_FEATURE    Feature mapping function for Exercise 5
%
%   map_feature(feat1, feat2) maps the two input features
%   to higher-order features as defined in Exercise 5.
%
%   Returns a new feature array with more features
%
%   Inputs feat1, feat2 must be the same size
%
% Note: this function is only valid for Ex 5, since the degree is
% hard-coded in.
    degree = 6;
    out = ones(size(feat1(:,1)));
    for i = 1:degree
        for j = 0:i
            out(:, end+1) = (feat1.^(i-j)).*(feat2.^j);
        end
    end

author: zellyn


text:
function err = getL2error3RT0(node,elem,sigma,sigmah,markedElem)
%% GETL2ERROR3RT0  L2 norm of RT0 element in 3D.
%
%  The input sigma can now just be function boundle. sigmah is the
%  flux through faces.
%  err = getL2error3RT0(node,elem,sigma,sigmah,markedElem)
%
% Example
%     node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];  % nodes
%     elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7]; % elements
%     elem = label3(node,elem);
%     [node,elem] = uniformbisect3(node,elem);
%     bdFace = setboundary3(node,elem,'Dirichlet','y~=-1','Neumann','y==-1');
%     maxIt = 3;
%     pde = mixBCdata3;
%     err = zeros(maxIt,1); N = zeros(maxIt,1);
%     for i =1:maxIt
%         [node,elem,~,bdFace] = uniformbisect3(node,elem,[],bdFace);
%         barycenter = 1/4.*(node(elem(:,1),:)+node(elem(:,2),:)+node(elem(:,3),:)+node(elem(:,4),:));
%         uexa = pde.exactu(barycenter);
%         [u,sigma] = Poisson3RT0(node,elem,pde,bdFace);
%         err(i) = getL2error3(node,elem,pde.exactu,u);
%         N(i) = length(u) + length(sigma) ;
%     end
%     r1 = showrate(N,err,2);
%     legend('||\sigma-\sigma_h||',['N^{' num2str(r1) '}'],...
%           'LOCATION','Best')
%
%% TODO: add jump coefficent and sigma is vector array.
%
% Created by Ming Wang at Jan 17, 2011, modified m-lint May 15, 2011.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Construct Data Structure
[elem2dof,dofSign] = dof3RT0(elem);
NT = size(elem,1);% Ndof = max(elem2dof(:)); %N = size(node,1);
[Dlambda,area] = gradbasis3(node,elem);
locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2];

%% Compute square of the L2 error element-wise
[lambda,w] = quadpts3(3); % quadrature order is 3
nQuad = size(lambda,1);
err = zeros(NT,1);
for p = 1:nQuad
    % quadrature points in the x-y-z coordinate
    pxy = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:) ...
        + lambda(p,4)*node(elem(:,4),:);
    sigmap = sigma(pxy);
    sigmahp = zeros(NT,3);
    for l = 1:4 % for each basis
        i = locFace(l,1); j = locFace(l,2); k = locFace(l,3);
        % phi_k = lambda_iRot_j - lambda_jRot_i;
        sigmahp = sigmahp + repmat(double(dofSign(:,l)).*sigmah(elem2dof(:,l)),1,3)*2.*...
                   (lambda(p,i)*cross(Dlambda(:,:,j),Dlambda(:,:,k)) + ...
                    lambda(p,j)*cross(Dlambda(:,:,k),Dlambda(:,:,i)) + ...
                    lambda(p,k)*cross(Dlambda(:,:,i),Dlambda(:,:,j)));
    end
    err = err + w(p)*sum((sigmap - sigmahp).^2,2);
end
err = err.*area;
% modify the error
err(isnan(err)) = 0;
if (nargin == 5) && ~isempty(markedElem)
    err = err(markedElem); % L2 err on some marked region
end
err = sqrt(sum(err));

repo: Aero-matlab


text:
function [y m]=convolution_(x,nx,h,nh)
m=nx(1)+nh(1):nh(end)+nx(end);
Lx=length(x);
Lh=length(h);
H=[h zeros(1,Lx)]'*ones(1,Lx);
H=H(1:end-Lx);
H=reshape(H,Lx+Lh-1,Lx);
y=(H*x')';
end

author: Sable

repo: mcbench-benchmarks


text:
function [x,options]=simps(fun,x,ix,options,vlb,vub,varargin)
%SIMPS (Strategy Simplex) finds the constrained minimum.
%
%   Syntax:
%   [X1,OPTIONS]=SIMPS(FUN,X0,IX,OPTIONS,VLB,VUB,P1,P2...)
%   where IX,OPTIONS,VLB,VUB,P1,P2... are optional.
%
%   Compare to Matlab's functions:
%   [X1,OPTIONS]=FMINS(FUN,X0,OPTIONS,[],P1,P2...)
%   [X1,OPTIONS]=CONSTR(FUN,X0,OPTIONS,VLB,VUB,[],P1,P2...) - OPTIM
%   toolbox
%
%   'FUN' - function to be minimized (usually an M-file: FUN.M).
%   It must be of the form: [F,G]=FUN(X,P1,P2...), where F returns the
%   value, and G returns the vector of constraints.
%   Constraints are satisfied at the point X if the respective
%   components of G are non-positive. For the non-constrained
%   minimization, use G=[].
%   The point X can be a [real] vector, matrix, or several dimensional
%   array. Both F and G must be well defined for all VLB<=X<=VUB (see
%   below).
%   P1,P2... are optional problem-dependent parameters, bypassed by
%   the minimizer SIMPS.
%
%   X0 - initial guess,
%   X1 - point of minimum as obtained;
%   both of the same type and dimension as X in the function 'FUN'.
%
%   IX - integer vector, subset of 1:PROD(SIZE(X0)) (which is also its
%   default value). The minimization is performed only against X0(IX),
%   while the other components of X0 are supposed to be fixed.
%   Let us denote LX=LENGTH(IX).
%
%   OPTIONS - see FOPTIONS. Only the following components are used:
%   Input:
%   OPTIONS(1)  - Whether to display intermediate results (1),
%                 or not (0, default).
%   OPTIONS(2)  - Relative X-termination tolerance (default: 1e-4).
%   OPTIONS(3)  - Relative F-termination tolerance (default: 1e-4).
%   OPTIONS(14) - Maximum number of function evaluations per each
%                 internal LX-dim. and 2-dim. simplex call
%                 (default: 100*LX).
%   Output:
%   OPTIONS(8)  - Function value at the point X1 (minimum).
%   OPTIONS(10) - Total number of function evaluations.
%
%   VLB,VUB - lower and upper bounds for X.
%   Both must be vectors of length PROD(SIZE(X0)) - total number of
%   components of X0.
%   Supply -INF, resp. +INF, (default) for unrestricted coordinates.
%
%   Comments:
%
%   - SIMPS is an iterative method, where each iteration consists of
%   one LX (full dim.) Nelder-Mead simplex call (direct search),
%   followed by LX+1 two dim. Nelder-Mead simplex calls (working on 2
%   dim. subspaces). These LX+1 two-dim. simplex calls provide initial
%   vertices V for the next full dim. simplex call (next iteration).
%
%   - The number of outer iterations is limited to LX or until the
%   longest dimension of simplex V (see above) falls below the
%   X-tolerance (OPTIONS(2)).
%   The internal LX and 2 dim. simplex calls are limited by proposed
%   F-tolerance or max. number of function evaluations (OPTIONS(2) and
%   OPTIONS(14)).
%
%   - The inner LX and 2 dim. Nelder-Mead minimizers are both realized
%   through the calls to the internal AMOEBA function, based
%   essentially on the code of the Matlab's FMINS non-linear simplex
%   implementation.
%
%   - The constraints are basically implemented by penalizing the
%   function (see the internal WRAPPER function): unsatisfied
%   constraints (-sum(G(find(G<0))), multiplied by factor 1e6, are
%   being added to the function value F.
%
%   - The index vector IX provides a useful and convenient method for
%   numerical experiments with complicated unknown minima, allowing to
%   user to perform searches limited to smaller dimensional subspaces.
%   This extra feature does not require run-time editing on the side
%   of target function 'FUN', as would be the case if certain
%   variables were fixed by use of optional parameters P1,P2...
%
%   - In contrast to the Matlab's FMINS, CONSTR and other minimizers,
%   both X and F termination tolerances are considered relative, not
%   absolute.
%   Nevertheless, X, F and G should be ideally of the same order of
%   magnitude (1), at least locally around the point of minimum.
%
%   - There is of course no guarantee that the output X1 represents
%   the point of global minimum as required.
%   However, the method (strategy) provides an advantage over the
%   plain nonlinear simplex, and it has been proved to be specially
%   useful for target functions with plenty of narrow local minima -
%   standard traps for analitically based minimizers.
%   The method is not limited to continuos functions and does not
%   require derivatives.
%
%   Authors:
%   Zeljko Bajzer (bajzer@mayo.edu) and Ivo Penzar (penzar@mayo.edu)
%   Mayo Clinic and Foundation, Rochester, Minnesota, USA
%   June, 1998.

if nargin<2, error('simps requires two input arguments'); end
if nargin<3, ix=1:prod(size(x)); end
ix=unique(ix); lx=length(ix);
if nargin<4, options=[]; end
if nargin<5, vlb=-Inf+zeros(size(1,lx)); end
if nargin<6, vub=Inf+zeros(size(1,lx)); end

options=foptions(options);
prnt=options(1);
if ~options(14)
 options(14)=50*lx;
end
tolx=options(2);
tolf=options(3);
maxfcalls=options(14);

varargs={vlb,vub,fun,varargin{:}};

xx=x(ix)';
if all(abs(xx)<eps)
 xx=xx+eps;
end
d=0.5;
% Guess goes to the centre:
%v=repmat((1-d/(lx+1))*xx,1,lx+1)+[diag(d*xx) zeros(lx,1)];
% Guess goes to a corner:
v=repmat(xx,1,lx+1)+[diag(d*xx) zeros(lx,1)];
fcalls=0; dv=Inf;
if prnt
 fprintf('%7s %10s %12s\n','f-COUNT','FUNCTION','x-TOL');
end
i=0;
while (i<lx)&(dv>=tolx)
 [y,fmin,tmp4]=amoeba(v,tolf,maxfcalls,...
                      x,ix,varargs{:});
 y=y'; x(ix)=y;
 fcalls=fcalls+tmp4;
 if prnt
  fprintf('%7d %10.5f\n',fcalls,fmin);
 end
 if lx==1
  break
 end
 if lx==2
  v(:,1)=y;
 else
  for j=1:lx-1
   j1=[j j+1];
   [y(j1),fmin,tmp4]=amoeba([xx(j1) y(j1) [xx(j1(1)) y(j1(2))]'],...
                            tolf,maxfcalls,...
                            x,ix(j1),varargs{:});
   x(ix)=y; v(:,j)=y;
   fcalls=fcalls+tmp4;
  end
 end
 j1=[1 lx];
 [y(j1),fmin,tmp4]=amoeba([xx(j1) y(j1) [xx(j1(1)) y(j1(2))]'],...
                          tolf,maxfcalls,...
                          x,ix(j1),varargs{:});
 x(ix)=y; v(:,lx)=y;
 fcalls=fcalls+tmp4;
 if lx>3
  j1=[2 lx-1];
 end
 [y(j1),fmin,tmp4]=amoeba([xx(j1) y(j1) [y(j1(1)) xx(j1(2))]'],...
                          tolf,maxfcalls,...
                          x,ix(j1),varargs{:});
 x(ix)=y; v(:,lx+1)=y;
 fcalls=fcalls+tmp4;
 dv=max(tol(min(v,[],2),max(v,[],2)));
 if prnt
  fprintf('%7d %10.5f %12.6f\n',fcalls,fmin,dv);
 end
 i=i+1;
end
if prnt
 fprintf('\n');
end

options(8)=fmin;
options(10)=fcalls;

return

function [x,fmin,fcalls]=amoeba(v,tolf,maxfcalls,varargin)
% Internal Nelder-Mead minimizer, based essentially on the code
% of the Matlab's FMINS non-linear simplex implementation.
%
% - Input and output parameters are changed
%   (e.g., initial simplex v is provided from outside);
% - There is no x-tolerance termination checking;
% - There is no display of intermediate results.

rho=1; chi=2; psi=0.5; sigma=0.5;

n=length(v(:,1));
onesn=ones(1,n);
two2np1=2:n+1;
one2n=1:n;
zerosn=zeros(1,n);
x=zerosn;

fv=zeros(1,n+1);
for j=1:n+1
 x(:)=v(:,j);
 fv(j)=wrapper(x,varargin{:});
end
[fv,jj]=sort(fv);
v=v(:,jj);
fmin=fv(1); fcalls=n+1;

while (fcalls<maxfcalls)&...
      (tol(fmin,fv(n+1))>tolf)
 shrink=0;
 xbar=sum(v(:,one2n),2)/n;
 xr=(1+rho)*xbar-rho*v(:,n+1);
 x(:)=xr;
 fxr=wrapper(x,varargin{:});
 fcalls=fcalls+1;

 if fxr<fv(1)
  xe=(1+rho*chi)*xbar-rho*chi*v(:,n+1);
  x(:)=xe;
  fxe=wrapper(x,varargin{:});
  fcalls=fcalls+1;
  if fxe<fxr
   v(:,n+1)=xe;
   fv(n+1)=fxe;
  else
   v(:,n+1)=xr;
   fv(n+1)=fxr;
  end

 elseif fxr<fv(n)
  v(:,n+1)=xr;
  fv(n+1)=fxr;

 else
  if fxr<fv(n+1)
   xc=(1+psi*rho)*xbar-psi*rho*v(:,n+1);
   x(:)=xc;
   fxc=wrapper(x,varargin{:});
   fcalls=fcalls+1;
   if fxc<=fxr
    v(:,n+1)=xc;
    fv(n+1)=fxc;
   else
    shrink=1;
   end
  else
   xcc=(1-psi)*xbar+psi*v(:,n+1);
   x(:)=xcc;
   fxcc=wrapper(x,varargin{:});
   fcalls=fcalls+1;
   if fxcc<fv(n+1)
    v(:,n+1)=xcc;
    fv(n+1)=fxcc;
   else
    shrink=1;
   end
  end
  if shrink
   for j=two2np1
    v(:,j)=v(:,1)+sigma*(v(:,j)-v(:,1));
    x(:)=v(:,j);
    fxcc=wrapper(x,varargin{:});
   end
   fcalls=fcalls+n;
  end
 end

 [fv,jj]=sort(fv);
 v=v(:,jj);
 fmin=fv(1);
end
x(:)=v(:,1);

return

function r=tol(a,b)
% Internal helper function, calculates the relative distances
% between points a and b, componentwise.

c=(abs(b)+abs(a))/2;
j=find(c<eps); c(j)=ones(size(j));
r=abs(b-a)./c;

return

function fmin=wrapper(x,x0,ix,vlb,vub,funfcn,varargin)
% Internal wrapper function, provides a proper interface for
% AMOEBA to the user's target function.

x1=zeros(1,prod(size(x0)));
x1(:)=x0;
x1(ix)=x;
x0(:)=max(min(x1,vub),vlb);

funfcn=fcnchk(funfcn,1+length(varargin));
[fmin,fconstr]=feval(funfcn,x0,varargin{:});

fconstr=[vlb-x1,x1-vub,fconstr];
fmin=fmin+1e6*sum(fconstr(find(fconstr>0)));

return


text:
function A = warmUpExercise()
%WARMUPEXERCISE Example function in octave
%   A = WARMUPEXERCISE() is an example function that returns the 5x5 identity matrix

A = [];
% ============= YOUR CODE HERE ==============
% Instructions: Return the 5x5 identity matrix
%               In octave, we return values by defining which variables
%               represent the return values (at the top of the file)
%               and then set them accordingly.

A = eye(5);

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

end

author: zlotus


text:
function [dtf, dtfvar, n] = ft_connectivity_dtf(inputdata, varargin)

% FT_CONNECTIVITY_DTF computes the directed transfer function.
%
% Use as
%   [d, v, n] = ft_connectivity_dtf(inputdata, ...)
%
% The input should be a spectral transfer matrix organized as
%   Nrpt x Nchan x Nchan x Nfreq (x Ntime)
% where Nrpt can be 1.
%
% The output represents
%   d = partial directed coherence matrix Nchan x Nchan x Nfreq (x Ntime).
%       If multiple observations in the input, the average is returned.
%   v = variance of d across observations.
%   n = number of observations.
%
% Typically, nrpt should be 1 where the spectral transfer matrix is computed across
% observations. When nrpt>1 and hasjack=true, the input is assumed to contain the
% leave-one-out estimates of the spectral transfer matrix, thus a more reliable
% estimate of the relevant quantities.
%
% Additional optional input arguments come as key-value pairs:
%   'hasjack'   = boolean, specifying whether the input contains leave-one-outs,
%                 required for correct variance estimate (default = false)
%   'crsspctrm' = matrix containing the cross-spectral density. If this
%                 matrix is defined, the function returns the ddtf, which
%                 requires an estimation of partial coherence from this matrix.
%   'invfun'    = 'inv' (default) or 'pinv', the function used to invert the
%                 crsspctrm matrix to obtain the partial coherence. Pinv is
%                 useful if the data are poorly-conditioned.
%   'feedback'  = 'none', 'text', 'textbar', 'dial', 'etf', 'gui' type of feedback showing progress of computation, see FT_PROGRESS
%
% See also CONNECTIVITY, FT_CONNECTIVITYANALYSIS

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

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

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

if ~isempty(powindx)
  % this error message is rather uninformative, but is kept for now for
  % backward compatibility reasons (i.e. it might exist when called from
  % ft_connectivityanalysis
  ft_error('linearly indexed data for DTF computation is at the moment not supported');
end

siz    = [size(inputdata) 1];
n      = siz(1);
outsum = zeros(siz(2:end));
outssq = zeros(siz(2:end));

% check the crsspctrm; if it is present, compute the partial coherence
if ~isempty(crsspctrm)
  assert(isequal(size(crsspctrm),size(inputdata)), 'the input data should be of the same size as the crsspctrm');
  fprintf('computing dDTF in the presence of a crsspctrm\n');

  % the crsspctrm allows for the partial coherence to be computed
  pdim   = prod(siz(4:end));
  tmpcrs = reshape(crsspctrm, [siz(1:3) pdim]);
  ft_progress('init', feedback, 'computing partial coherence...');
  for k = 1:n
    ft_progress(k/n, 'computing partial coherence for replicate %d from %d\n', k, n);
    tmp = reshape(tmpcrs(k,:,:,:), [siz(2:3) pdim]);
    for m = 1:pdim
      tmp(:,:,m) = invfun(tmp(:,:,m));
      tmp(:,:,m) = abs(tmp(:,:,m))./sqrt(abs(diag(tmp(:,:,m))*diag(tmp(:,:,m))'));
    end
    tmpcrs(k,:,:,:) = tmp;
  end
  ft_progress('close');
  crsspctrm = reshape(tmpcrs, siz);
end

% data should be represented as chan_chan_therest
for j = 1:n
  tmph   = reshape(inputdata(j,:,:,:,:), siz(2:end));
  if isempty(crsspctrm)
    % plain DTF
    den    = sum(abs(tmph).^2,2);
    tmpdtf = abs(tmph)./sqrt(repmat(den, [1 siz(2) 1 1 1]));
  else
    % dDTF
    tmpc   = reshape(crsspctrm(j,:,:,:,:), siz(2:end));
    den    = sum(sum(abs(tmph).^2,3),2);
    tmpdtf = abs(tmph)./sqrt(repmat(den, [1 siz(2) siz(4) 1 1 1]));
    tmpdtf = tmpdtf.*tmpc;
  end
  outsum = outsum + tmpdtf;
  outssq = outssq + tmpdtf.^2;
end
dtf = outsum./n;

if n>1 % FIXME this is strictly only true for jackknife, otherwise other bias is needed
  if hasjack
    bias = (n - 1).^2;
  else
    bias = 1;
  end
  dtfvar = bias.*(outssq - (outsum.^2)/n)./(n-1);
else
  dtfvar = [];
end

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

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function g = newDomain(g, newDom)
%NEWDOMAIN   Change of domain of a CHEBFUN.
%  NEWDOMAIN(G, DOM) returns the CHEBFUN G but moved to the domain DOM. This is
%  done with a linear map. DOM may be a vector of length G.ends, or a two-vector
%  (in which case all breakpoints are scaled by the same amount).

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

% [TODO]: Unbounded domains.

% Current breakpoints:
oldDom = g.domain;

if ( numel(newDom) == numel(oldDom) )
    % All new breakpoints are given!
elseif ( numel(newDom) == 2 )
    % Scale breakpoints:
    c = oldDom(1);
    d = oldDom(end);
    a = newDom(1);
    b = newDom(2);
    newDom = (b - a)*(oldDom - c)/(d - c) + a;
else
    error('CHEBFUN:CHEBFUN:newDomain:numints', 'Inconsistent domains.');
end

for k = 1:numel(g.funs)
    % Update the domains of each of the funs:
    g.funs{k} = changeMap(g.funs{k}, newDom(k:k+1));
end

% Update the CHEBFUN:
g.domain = newDom;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function d = bb_distance(bb1,bb2)
% Info

switch nargin
    case 1
        d = 1 - bb_overlap(bb1);
    case 2
        d = 1 - bb_overlap(bb1,bb2);
end

repo: OpenTLD


text:
function varargout = gui_eventinsertion(varargin)
% GUI_EVENTINSERTION MATLAB code for gui_eventinsertion.fig
%      GUI_EVENTINSERTION, by itself, creates a new GUI_EVENTINSERTION or raises the existing
%      singleton*.
%
%      H = GUI_EVENTINSERTION returns the handle to a new GUI_EVENTINSERTION or the handle to
%      the existing singleton*.
%
%      GUI_EVENTINSERTION('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GUI_EVENTINSERTION.M with the given input arguments.
%
%      GUI_EVENTINSERTION('Property','Value',...) creates a new GUI_EVENTINSERTION or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before gui_eventinsertion_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to gui_eventinsertion_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help gui_eventinsertion

% Last Modified by GUIDE v2.5 28-Oct-2010 04:42:36

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @gui_eventinsertion_OpeningFcn, ...
                   'gui_OutputFcn',  @gui_eventinsertion_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before gui_eventinsertion is made visible.
function gui_eventinsertion_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to gui_eventinsertion (see VARARGIN)

% Choose default command line output for gui_eventinsertion
handles.output = {};

% Update handles structure
guidata(hObject, handles);

% invoke the callback to init the dialog
edit1_Callback(handles.edit1,{},guidata(hObject));

% UIWAIT makes gui_eventinsertion wait for user response (see UIRESUME)
%uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = gui_eventinsertion_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% hack to get the guipanel to render on Macs...
if ismac
    p=get(handles.uipanel1,'Position');
    set(handles.uipanel1,'Position',p.*[ 1 1 0.99 1]);
    set(handles.uipanel1,'Position',p.*[ 1 1 1/0.99 1]);
end
uiwait(handles.figure1);
handles = guidata(hObject);

try
    % Get default command line output from handles structure
    varargout{1} = handles.output;
    close(handles.figure1);
catch
    varargout{1} = {};
end

% --- set number of insertion operations!

function edit1_Callback(hObject, eventdata, handles)
gui_check_scalar(hObject,'1');
num = max(1,round(str2num(get(hObject,'String'))));
set(hObject,'String',num2str(num));

if ~isfield(handles,'tabpanel1') || length(get(handles.tabpanel1,'Title')) ~= num
    % generate titles
    title = {};
    for k=1:num
        title{end+1} = ['Step ' num2str(k)]; end
    % save old parameters
    if ~isfield(handles,'opgrids')
        handles.opgrids = {}; end
    parameters = cell(1,num);
    for h=1:min(num,length(handles.opgrids))
        parameters{h} = arg_tovals(handles.opgrids{h}.GetPropertySpecification(),false); end
    % delete old propertygrids
    for h=1:length(handles.opgrids)
        delete(handles.opgrids{h}); end
    % (re-) create tabpanel
    if isfield(handles,'tabpanel1')
        item = num;
        delete(handles.tabpanel1);
    else
        item = 1;
    end
    handles.tabpanel1 = uitabpanel('active',item,'parent',handles.uipanel2,'title',title,'TitleBackgroundColor',[0.93,0.91,0.85]*0.66,'PanelBackgroundColor',[0.93,0.91,0.85],'frameBackgroundColor',[0.66,0.76,1],'TitleForegroundColor',[0 0 0]);
    %set(handles.tabpanel1,'SelectedItem',item);
    subhandles = getappdata(handles.tabpanel1,'panels');
    % create new propertygrids
    handles.opgrids = {};
    for h=1:length(subhandles)
        handles.opgrids{h} = arg_guipanel(subhandles(h), 'Function',@set_insert_markers, 'Parameters',parameters(h)); end
end
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% get the parameter specification for all steps
handles.output = {};
for h=1:length(handles.opgrids)
    handles.output{h} = arg_tovals(handles.opgrids{h}.GetPropertySpecification(),false); end
guidata(hObject,handles);
uiresume(handles.figure1);

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.output = {};
guidata(hObject,handles);
uiresume(handles.figure1);

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
env_doc set_insert_markers

% --- Executes during object creation, after setting all properties.
function uipanel1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to uipanel1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

repo: REST


text:
function normalizedU = hyperConvexHullRemoval(U,wavelengths)
%HYPERCONVEXHULLREMOVAL Performs spectral normalization via convex hull removal
%
% Usage
%   [ normalizedU ] = hyperConvexHullRemoval( U, wavelengths )
%
% Inputs
%   U - 2D HSI data (p x q)
%   wavelengths - Wavelength of each band (p x 1)
%
% Outputs
%   normalizedU - Data with convex hull removed (p x q)
%
% Author
%   Luca Innocenti
%
% References
%   Clark, R.N. and T.L. Roush (1984) Reflectance Spectroscopy: Quantitative
% Analysis Techniques for Remote Sensing Applications, J. Geophys. Res., 89,
% 6329-6340.

% Metadata and formatting
wavelengths = wavelengths(:);
p = length(wavelengths);
q = size(U,2);
U = U.';

U(:,1) = 0;
U(:,420) = 0;

normalizedU = zeros(q,420);

% The algorithm
for s = 1:q,
    rifl = U(s,:);
    k = convhull(wavelengths,rifl');
    c = [rifl(k); wavelengths(k)'];
    d = sortrows(c',2);

    xs = d(:,2);
    ys = d(:,1);
    [xsp, idx] = unique(xs);
    ysp = ys(idx);
    rifl_i = interp1(xsp,ysp,wavelengths');

    for t = 1:420,
        if rifl_i(t) ~= 0
            normalizedU(s,t) = rifl(t)/rifl_i(t);
        else
            normalizedU(s,t) = 1;
        end
    end
end

normalizedU = normalizedU.';

author: davidkun

repo: HyperSpectralToolbox


text:
function DEM = uplus(DEM)
DEM.Z = uplus(DEM.Z);

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
%Esto es el PURPOSE de la FUNCION que es sumar a + b
function [ c ] = funct1( a,b )
%FUNCT1 Summary of this function goes here
%   Detailed explanation goes here
%   Esta funcion suma a + b
c=a+b;
t=0:0.001:1;
plot(sin(t));
end


text:
function [mix, num_iter, ll] = gmmem_kpm(mix, x, varargin)
%GMMEM_KPM Like GMMEM, but with additional optional arguments
% function [mix, num_iter, ll] = gmmem_kpm(mix, x, varargin)
%
% Input:
% mix - structure created by gmminit or gmmem_multi_restart
% data - each row is an example
%
% Output:
% mix - modified structure
% num_iter - number of iterations needed to reach convergence
% ll - final log likelihood
%
% [ ... ] = gmmem_kpm(..., 'param1',val1, 'param2',val2, ...) allows you to
% specify optional parameter name/value pairs.
% Parameters are below [default value in brackets]
%
% 'max_iter' - maximum number of EM iterations [10]
% 'll_thresh' - change in log-likelihood threshold for convergence [1e-2]
% 'verbose' - 1 means display output while running [0]
% 'prior_cov' - this will be added to each estimated covariance
%               to prevent singularities  [1e-3*eye(d)]
% 'fn'        - this function, if non-empty,  will be called at every iteration
%               (e.g., to display the parameters as they evolve) [ [] ]
%               The fn is called as fn(mix, x, iter_num, fnargs).
%               It is also called before the iteration starts as
%               fn(mix, x, -1, fnargs), which can be used to initialize things.
% 'fnargs'    - additional arguments to be passed to fn [ {} ]
%
% Modified by Kevin P Murphy, 29 Dec 2002

% Check that inputs are consistent
errstring = consist(mix, 'gmm', x);
if ~isempty(errstring)
  error(errstring);
end

[ndata, xdim] = size(x);

[max_iter, ll_thresh, verbose, prior_cov, fn, fnargs] = ...
    process_options(varargin, ...
	'max_iter', 10, 'll_thresh', 1e-2, 'verbose', 1, ...
	'prior_cov', 1e-3*eye(xdim), 'fn', [], 'fnargs', {});

options = foptions;
if verbose, options(1)=1; else options(1)=-1; end
options(14) = max_iter;
options(3) = ll_thresh;

% Sort out the options
if (options(14))
  niters = options(14);
else
  niters = 100;
end

display = options(1);
test = 0;
if options(3) > 0.0
  test = 1;	% Test log likelihood for termination
end

check_covars = 0;
if options(5) >= 1
  if display >= 0
    disp('check_covars is on');
  end
  check_covars = 1;	% Ensure that covariances don't collapse
  MIN_COVAR = eps;	% Minimum singular value of covariance matrix
  init_covars = mix.covars;
end

mix0 = mix; % save init values for debugging

if ~isempty(fn)
  feval(fn, mix, x, -1, fnargs{:});
end

% Main loop of algorithm
for n = 1:niters

  % Calculate posteriors based on old parameters
  [post, act] = gmmpost(mix, x);

  % Calculate error value if needed
  if (display |  test)
    prob = act*(mix.priors)';
    % Error value is negative log likelihood of data
    e = - sum(log(prob + eps));
    if display > 0
      fprintf(1, 'Cycle %4d  Error %11.6f\n', n, e);
    end
    if test
      if (n > 1 & abs(e - eold) < options(3))
        options(8) = e;
	ll = -e;
	num_iter = n;
        return; %%%%%%%%%%%%%%%% Exit here if converged
      else
        eold = e;
      end
    end
  end

  if ~isempty(fn)
    feval(fn, mix, x, n, fnargs{:});
  end

  % Adjust the new estimates for the parameters
  new_pr = sum(post, 1);
  new_c = post' * x;

  % Now move new estimates to old parameter vectors
  mix.priors = new_pr ./ ndata;

  mix.centres = new_c ./ (new_pr' * ones(1, mix.nin));

  switch mix.covar_type
  case 'spherical'
    n2 = dist2(x, mix.centres);
    for j = 1:mix.ncentres
      v(j) = (post(:,j)'*n2(:,j));
    end
    mix.covars = ((v./new_pr) + sum(diag(prior_cov)))./mix.nin;
    if check_covars
      % Ensure that no covariance is too small
      for j = 1:mix.ncentres
        if mix.covars(j) < MIN_COVAR
          mix.covars(j) = init_covars(j);
        end
      end
    end
  case 'diag'
    for j = 1:mix.ncentres
      diffs = x - (ones(ndata, 1) * mix.centres(j,:));
      wts = (post(:,j)*ones(1, mix.nin));
      mix.covars(j,:) = sum((diffs.*diffs).*wts + prior_cov, 1)./new_pr(j);
    end
    if check_covars
      % Ensure that no covariance is too small
      for j = 1:mix.ncentres
        if min(mix.covars(j,:)) < MIN_COVAR
          mix.covars(j,:) = init_covars(j,:);
        end
      end
    end
  case 'full'
    for j = 1:mix.ncentres
      diffs = x - (ones(ndata, 1) * mix.centres(j,:));
      diffs = diffs.*(sqrt(post(:,j))*ones(1, mix.nin));
      mix.covars(:,:,j) = (diffs'*diffs + prior_cov)/new_pr(j);
    end
    if check_covars
      % Ensure that no covariance is too small
      for j = 1:mix.ncentres
        if min(svd(mix.covars(:,:,j))) < MIN_COVAR
          mix.covars(:,:,j) = init_covars(:,:,j);
        end
      end
    end
  case 'ppca'
    for j = 1:mix.ncentres
      diffs = x - (ones(ndata, 1) * mix.centres(j,:));
      diffs = diffs.*(sqrt(post(:,j))*ones(1, mix.nin));
      [mix.covars(j), mix.U(:,:,j), mix.lambda(j,:)] = ...
        ppca((diffs'*diffs)/new_pr(j), mix.ppca_dim);
    end
    if check_covars
      if mix.covars(j) < MIN_COVAR
        mix.covars(j) = init_covars(j);
      end
    end
    otherwise
      error(['Unknown covariance type ', mix.covar_type]);
  end
end

ll = sum(log(gmmprob(mix, x)));
num_iter = n;

%if (display >= 0)
%  disp('Warning: Maximum number of iterations has been exceeded');
%end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function varargout = clusterImg(imgin,Clusters,varargin)

% Copyright 2010 The MathWorks, Inc.
% [IDX,C,sumd,D] = clusterImg(imgin,nClusters,inMask,labels)
%
% This function implements kmeans clustering on an input RGB (m x n x 3)
% image. The user inputs at least two inputs: IMGIN and
% NCLUSTERS, and this
% function will step through an interactive color segmentation using kmeans
% clustering. It is interactive in that the user will be prompted to
% click-define the target colors in the segmentation. (There will be
% NCLUSTERS prompts.) Optionally, the user can also input a third argument
% containing a binary mask of the input image, in which case the program
% will operate on a subimage defined by the BOUNDING BOX of the (true)
% inMask.
%
% NOTE 1: Following segmentation (clustering), the user can EXPORT THE CLUSTERED
% IMAGE TO THE WORKSPACE, along with its corresponding colormap. Additionally,
% the user can opt to EXPORT A BINARY MASK OF A SINGLE COLOR CLUSTER TO THE
% WORKSPACE. Both of these options are initiated by press of a uicontrol
% pushbutton. The naming convention is automatic; the clustered image and
% corresponding colormap will be named CLUSTERED1 and CLUSTEREDMAP1,
% respectively. If those variables already exist in the workspace, they
% will NOT be overwritten. Rather, the names will be modified by
% incrementing the integer value on the end (e.g., CLUSTERED2,
% CLUSTERMAP2). Similarly, the binary mask on a cluster will be named
% CLUSTERMASK, CLUSTERMASK1, CLUSTERMASK2,....
%
% NOTE 2: To reconstruct the image from workspace variables, use the command:
%       imshow(clustered,clusterMap);
%
% ***************************************
% SYNTAX:
%
% [IDX,C,sumd,D] = clusterImg(imgin,nClusters,inMask,labels);
%
% ***************************************
% INPUTS:
%
% Imgin:
%     an m x n x 3 (RGB) image
%
% Clusters:
%     Either A) an integer number of target clusters. (Note
%     that you will be prompted to click-define--using
%     IMPIXEL syntax--at least one pixel in each target
%     color cluster; you probably want this to be a
%     "reasonable" number (more than 1, less than, say,
%     20); or B) an nx3 array of predefined colors. Using
%     option B will skip the interactive selection process
%     and use the colors in Clusters as the bin targets.
%
% inMask (OPTIONAL)
%     a binary mask of size m x n. If a mask is included, this function
%     will operate on a sub-image defined by IMCROP(imgin,BoundingBox),
%     where BoundingBox is the rectangular extent of the true mask.
%
% labels (OPTIONAL)
%     a cell array of labels for the colormap display of the
%     segmented image. Note that the number of strings must
%     be equal to the number of colors.
%
% ***************************************
% OUTPUTS:
%
% IDX, C, sumd, D
%     All as returned by KMEANS. Please refer to the help for that
%     function. Note, however, that the size of the outputs may be modified
%     by the inclusion of optional inMask input argument. (See NOTE 2
%     below.)
%
% clustered
%     an indexed image of size p x q x 1 (see NOTE 4 below), segmented using kmeans
%     clustering. (See description of optional inMask input argument). If no
%     mask is provided, the CLUSTERED image is m x n x 1.
%
% clusteredMap
%     an Clusters x 3 matrix of RGB colors. These colors comprise the mean
%     values returned by each IMPIXEL selection when the user defines the
%     target cluster colors.
%
% clusterMask
%     a binary mask of size p x q x 1 (see NOTE 2 below). Once the input has been
%     clustered, it will be displayed in a figure with a pushbutton prompt
%     to create a mask on a clustered color. If the button is pressed, the
%     user will be prompted to click-select a single cluster color (using
%     IMPIXEL syntax). CLUSTER_MASK will return a binary image with ones in
%     the location of the selected color, and zeros elsewhere.
%
% NOTE 3: OUTPUT ARGUMENTS [CLUSTERED, CLUSTEREDMAP, CLUSTERMASK] ARE
% PROMTED BY UICONTROLS, and do not require the presence of output
% arguments. Nomenclature is automatic, as described in NOTE 1 above.)
%
% NOTE 4: OUTPUT IMAGE (AND MAP) SIZES are p x q, if a binary mask is provided
% as an optional input, where p and q are defined by the extent of the
% bounding box of the (true) input mask. If no input mask is provided, p=m
% and q=n.
%
%
% ***************************************
% EXAMPLES:
%
% 1)
% clusterImg(imread('peppers.png'),6);
%     segments the image in 'peppers.png' into 6 color
%     clusters.
% 2)
% a=[0.2902    0.1647    0.2824;
%    0.9137    0.7922    0.6784;
%    0.3725    0.4392    0.1412;
%    0.5961    0.1529    0.2039;
%    1.0000    0.7843    0.0431;
%    0.9686    0.5294         0];
% clusterImg(imread('peppers.png'),a,[],{'One','Two','Three','Four','Five','Six'});
%
% 3)
% [IDX, C, sumd, D] = clusterImg(x, 5, x(:,:,1) < 100);
%     masks image x with "red plane less than 100", finds the bounding box
%     of the resulting (logically true) inMask, and segments the IMCROPped
%     image x(inMask) into 5 color clusters. Returns the variables IDX, C,
%     sumd, and D, as described by the documentation for KMEANS.
%
% ***************************************

% ***************************************
% DEPENDENCIES:
% Image Processing Toolbox, Statistics Toolbox
%
% ***************************************
% Written by Brett Shoelson, PhD
% brett.shoelson@mathworks.com
%
% Ver 2.0
% Revisions:
% Ver 1.1; 11/12/06. Bounding box detection on imcrop modified.
% Ver 1.2; 1/31.07. The second input argument has been
%    renamed from nClusters to Clusters to reflect the fact
%    that one can now either A) specify the number of clusters
%    and select them interactively (as before); or B)input an
%    nx3 array of colors, skipping the interactive
%    selection. ALSO, the function now takes an additional
%    optional input argument to specify the labels on the
%    colormap of the segmented image.
% Ver 2.0; 11/29/07. Now includes a pushbutton for the auto-generation of
%    all segmented colormap. Also, included description of LABELS as input
%    argument, and fixed incorrect capitalization of "clusterMaskn"
%    indicator.
% ***************************************

%Parse and initialize input and output arguments
if ndims(imgin)~=3
	error('Input image must be m x n x 3 (RGB)');
end
if nargin == 3
	inMask = varargin{1};
	if all(inMask(:))
		%NO OP: inMask is all 1's...no masking performed
		inMask = [];
	end
else
	inMask = [];
end

if ~isempty(inMask) %at least one excluded point; operate on masked sub-image
	[r,c]=find(inMask);
	bb = [min(c) min(r) max(c)-min(c) max(r)-min(r)];
	imgin = imcrop(imgin,bb);
	%bb = regionprops(bwlabel(inMask),'BoundingBox');
	%imgin = imcrop(imgin,bb.BoundingBox);
end

%Operate on image as double, regardless of input type
if ~isa(imgin,'double')
	imgin = im2double(imgin);
end

%Get individual RGB colors included in segmentation
rc = imgin(:,:,1);
gc = imgin(:,:,2);
bc = imgin(:,:,3);
tmpimg = imgin;
segcolors = double([rc(:) gc(:) bc(:)]);

tmpfig=figure('numbertitle','off','name','Cluster RGB Image','visible','off');
imshow(tmpimg);
set(tmpfig,'units','pixels','position',[0 0 1000 800]);
centerfig(tmpfig);
set(tmpfig,'visible','on');

if numel(Clusters)==1
	% Prompt for target cluster colors
	colors = zeros(Clusters,3);
	for ii=1:Clusters
		figure(tmpfig);
		title(sprintf('Select sample pixel(s) in cluster %d of %d:',ii,Clusters));
		tmp = impixel;
		colors(ii,:) = mean(tmp,1);
	end
elseif size(Clusters,2) == 3
	% User passed in the target colors; skip interactive
	% selection
	colors = Clusters;
	Clusters = size(colors,1);
end
title('Segmenting. Please wait.......');
drawnow;
%Operate on colors (and image) as doubles()
colors = im2double(colors);

% CALL KMEANS with appropritate output arguments.
if nargout < 2
	Idx = kmeans(segcolors,Clusters,'distance','sqEuclidean','start',colors,'emptyaction','singleton');
elseif nargout == 2
	[Idx,varargout{2}] = ...
	   kmeans(segcolors,Clusters,'distance','sqEuclidean','start',colors,'emptyaction','singleton');
elseif nargout == 3
	[Idx, varargout{2}, varargout{3}] = ...
		kmeans(segcolors,Clusters,'distance','sqEuclidean','start',colors,'emptyaction','singleton');
elseif nargout == 4
		[Idx, varargout{2}, varargout{3}, varargout{4}] = ...
		kmeans(segcolors,Clusters,'distance','sqEuclidean','start',colors,'emptyaction','singleton');
else
	error('Inappropriate number of outputs requested in ClusterImg.');
end
if nargout >= 2
	varargout{1} = Idx;
end

pixel_labels = reshape(Idx,size(rc));
imshow(pixel_labels,[]);
title('Image labeled by cluster index');
set(tmpfig,'units','pixels','position',[0 0 1000 800]);
centerfig(tmpfig);

colormap(colors);
tmp = colorbar;
set(tmp,'ytick',1:Clusters);
if nargin == 4
	set(tmp,'yticklabel',varargin{2});
end
uicontrol('style','pushbutton','string','Export Clustered Image and Map to Workspace','units','normalized',...
	'pos',[0.05 0.05 0.3 0.05],'callback',@exportMain);
uicontrol('style','pushbutton','string','Export Mask on Single Color','units','normalized',...
	'pos',[0.35 0.05 0.3 0.05],'callback',@exportClusterMask);
uicontrol('style','pushbutton','string','Export All Color Masks','units','normalized',...
	'pos',[0.65 0.05 0.3 0.05],'callback',{@exportAllClusterMasks,colors});

%Subfunctions local to clustering only. Isn't MATLAB wonderful?
	function exportMain(varargin)
		n = 0; tmp = 1;
		while tmp
			n = n + 1;
			tmp = evalin('base',['exist(''clustered', num2str(n), ''')']);
		end
		fprintf('Image written to clustered%d and associated colormap written to clusteredMap%d\n',n,n);
		fprintf('Display with:    imshow(clustered%d,clusteredMap%d);\n',n,n);
		assignin('base',['clustered' num2str(n)],pixel_labels);
		assignin('base',['clusteredMap' num2str(n)],colors);
	end

	function exportClusterMask(varargin)
		title('Click on any pixel in the desired region.');
		p=impixel;
		if ~all(p(:)==p(1))
			beep;
			disp('Select a SINGLE region!!!');
		else
			n = 0; tmp = 1;
			while tmp
				n = n + 1;
				tmp = evalin('base',['exist(''clusterMask', num2str(n), ''')']);
			end
			title('Image labeled by cluster index');
			fprintf('Mask of region %d written to clusterMask%d\n',p(1),n);
			assignin('base',['clusterMask' num2str(n)],pixel_labels==p(1));
		end
    end

    function exportAllClusterMasks(varargin)
        colors = varargin{3};
        n = 0; tmp = 1;
        while tmp
            n = n + 1;
            tmp = evalin('base',['exist(''clusterMask', num2str(n), ''')']);
        end
        for ii = 1:size(colors,1)
            fprintf('Mask of region %d written to clusterMask%d\n',ii,n+ii-1);
            assignin('base',['clusterMask' num2str(n+ii-1)],pixel_labels==ii);
        end
    end

end


text:
%% cell get_algorithms_list(int)
%
function [alg_list] = get_algorithm_list(id)
  L = load_algorithm_list();
  [m,~] = size(L{1,1});
  alg_list = {};
  k = 1;
  for i = 1 : m
    method_id_field = char(L{1,1}(i));
    method_id_list = strsplit(method_id_field,',');

    [~,s] = size(method_id_list);

    flag = 0;
    for j = 1 : s
      method_id = strtrim(char(method_id_list{1,j}));
      if(strcmp(method_id,'-'))
        alg_list{k,1} = '-';
        alg_list{k,2} = '-';
        alg_list{k,3} = 0;
        alg_list{k,4} = '-';
        k = k + 1;
        break;
      end
%       if(method_id(1) == '-')
%         list{i,1} = '-';
%         break;
%       end
      if(strcmp(id,method_id))
        flag = 1;
      end
    end

    if(flag == 0)
      continue;
    end

    algorithm_id = strtrim(char(L{1,2}(i)));
    algorithm_name = strtrim(char(L{1,3}(i)));

    algorithm = [algorithm_id ' ' algorithm_name];

    alg_list{k,1} = algorithm_id;
    alg_list{k,2} = algorithm_name;
    alg_list{k,3} = L{1,4}(i);
    alg_list{k,4} = algorithm;
    k = k + 1;
  end
end

author: andrewssobral

repo: lrslibrary


text:
%% Create dataset object
% Our bonds are located in a number of XLS-files.  Through a simple
% FOR-loop, we can grab all of the bonds and place them into a simgle
% workspace variable.
data = dataset;
for yearIdx = 1980:2000
    fileName = ['bonds' num2str(yearIdx) '.xls'];
    data = [data; dataset('xlsfile', fileName)];
end
% This single command is all we need to work with the data set.  There are,
% however, three improvements that we can make by applying our knowledge of
% the nature of this dataset.  These changes will make operations on the
% array have a simpler syntax and operate faster, and they will make the
% memory footprint smaller as an added bonus.

% Variable size is 1370656 bytes.
clearvars -except data
%% Tuning the dataset
%% 1. Nominal arrays
% Three fields (the ticker symbol, the description, and the currency) are
% filled with often-repeated, categorical data.  These should be recast
% into nominal arrays, which are more efficient at storing the repeated
% data and are also much faster for indexing and searching purposes.

data.ticker      = nominal(data.ticker);
data.description = nominal(data.description);
data.currency    = nominal(data.currency);
% Variable size is now 859844 bytes: 63% of the original size.

%% 2. Ordinal arrays
% The "index quality" field gives the bond ratings as assigned by S&P.
% These are similar to nominal array in that there are often-repeated
% labels, but furthermore there is now a definite order in that some
% ratings are better than others.  This is a perfect use for ordinal
% arrays.  All we must do is supply the ratings scale from lowest to
% highest (or else it will default to alphabetical order):

ratingsScale = ...
    {'BBB-', 'BBB', 'BBB+', 'A-', 'A', 'A+', 'AA-', 'AA', 'AA+', 'AAA'};
data.indexQuality = ordinal(data.indexQuality, [], ratingsScale);
% Variable size is now 701570 bytes: 51% of the original size.

%% 3. Change dates into MATLAB date numbers
% The three date columns are, by default, imported as character strings.
% Using the DATENUM command we can convert them to MATLAB's numeric data
% format, which allows for relational indexing.  It also uses less memory.

data.valDate      = datenum(data.valDate);
data.maturityDate = datenum(data.maturityDate);
% Some small savvy is required to properly handle the blank dates in the
% data.call_date set:
data.callDate     = cellfun(@datenan, data.callDate);
% Variable size is now 206982 bytes: only 15% of the original size!

%% 4. Load another dataset with sector information
% If we have data in different locations we can read in the data into
% another dataset and tune the dataset.
sectorInfo = dataset('xlsfile', 'sectorInfo.xls');

sectorInfo.ticker = nominal(sectorInfo.ticker);
sectorInfo.description = nominal(sectorInfo.description);
sectorInfo.sector = nominal(sectorInfo.sector);

%% 5. Join datasets together
% Once the dataset is tuned we can use the join function to join the two
% datasets based on a similar column or key.  Here we join the datasets
% based on the ticker to include sector information that was not in the
% original data we brought in.
dataFull = join(data, sectorInfo, 'key', 'ticker', 'RightVars', 'sector');

%% Some simple examples of data manipulation and indexing using datasets
% Now that we have a finely-tuned dataset, we can peform a number of
% intuitive indexing tasks as well as some operations.
%% 1. Searching for all instruments of a given ticker name
% First, let us list all of the ticker labels in our array.  This
% information is stored in the data.ticker nominal array.  We can list all
% categories in a nominal array with the GETLABELS function...

getlabels(dataFull.ticker)
%%
% ...and we can focus on a specific ticker (say, Good Foods) with array
% indexing:

dataFull( dataFull.ticker == 'GOOD', {'description', 'par', 'marketValue'} )

%% 2. Searching for all instruments of certain grades
% Thanks to ordinal arrays, we can search for all high- and prime-grade
% instruments with a single command:

dataFull( dataFull.indexQuality >= 'AA-', {'ticker', 'par', 'indexQuality'})

%% 3. Searching for all instruments of certain valuation dates
% With MATLAB serial dates, we can use relational operators to find all
% instruments between certain dates:
subset = dataFull( dataFull.valDate >= datenum('Dec-01-1980') ...
    & dataFull.valDate <= datenum('Jan-01-1981'), {'description', 'coupon', 'par'} )

%%
% Incidentally, it is straightforward to export a dataset to a variety of
% delimited files.  Here, we export the subset created above to a
% tab-delimited file:
export( subset, 'XLSFile', 'Dec1980.xls')

%% 4. A simple example of taking an average on a 2-month block of data
% Let us suppose we wanted to find the simple average coupon and par value
% of the instruments with a valuation date between 12/01/80 and 01/01/81,
% (the same date range as above):
avgCouponAndPar = datasetfun( @mean, dataFull( dataFull.valDate >= datenum('Dec-01-1980') ...
    & dataFull.valDate <= datenum('Jan-01-1981'), {'coupon', 'par'}) )

%% 5. A way to easily look at group statisics
% If we wanted to see the mean of the coupon for each of the instruments,
% we can use the GRPSTATS function:
avgCoupon = grpstats(dataFull, 'ticker', {@mean, @std}, 'DataVars', 'coupon')

%% 6. Other visualizations
% We can create a heatmap function that allows us to quickly look at how
% our portfolio is distributed across the different bonds.
generate_heatmap(dataFull);


text:
function cg_rc_test ( )

%*****************************************************************************80
%
%% CG_RC_TEST tests the CG_RC library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CG_RC_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the CG_RC library.\n' );

  cg_test01 ( );
  cg_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CG_RC_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [ n_data, a, b, n, fx ] = cos_power_int_values ( n_data )

%*****************************************************************************80
%
%% COS_POWER_INT_VALUES returns some values of the cosine power integral.
%
%  Discussion:
%
%    The function has the form
%
%      COS_POWER_INT(A,B,N) = Integral ( A <= T <= B ) ( cos(T) )^N dt
%
%    In Mathematica, the function can be evaluated by:
%
%      Integrate [ ( Cos[x] )^n, { x, a, b } ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0
%    before the first call.  On each call, the routine increments N_DATA by 1,
%    and returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real A, B, the limits of integration.
%
%    Output, integer N, the power.
%
%    Output, real FX, the function value.
%
  n_max = 11;

  a_vec = [ ...
     0.00, ...
     0.00, ...
     0.00, ...
     0.00, ...
     0.00, ...
     0.00, ...
     0.00, ...
     0.00, ...
     0.00, ...
     0.00, ...
     0.00 ];
  b_vec = [ ...
     3.141592653589793, ...
     3.141592653589793, ...
     3.141592653589793, ...
     3.141592653589793, ...
     3.141592653589793, ...
     3.141592653589793, ...
     3.141592653589793, ...
     3.141592653589793, ...
     3.141592653589793, ...
     3.141592653589793, ...
     3.141592653589793 ];
  fx_vec = [ ...
     3.141592653589793, ...
     0.0, ...
     1.570796326794897, ...
     0.0, ...
     1.178097245096172, ...
     0.0, ...
     0.9817477042468104, ...
     0.0, ...
     0.8590292412159591, ...
     0.0, ...
     0.7731263170943632 ];
  n_vec = [ ...
     0, ...
     1, ...
     2, ...
     3, ...
     4, ...
     5, ...
     6, ...
     7, ...
     8, ...
     9, ...
    10 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    b = 0.0;
    n = 0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    b = b_vec(n_data);
    n = n_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function sleep_stage = ECG_sleep_staging_prediction(ECGdata,fs,class)
% sleep stage prediction from ECG signal by model trained from SHHSv1 DB
%
% input:
%   ECGdata: ECG data
%   fs: sampling frequency of ECG
%   class: sleep staging type: class = 1, four types classification
%                              class = 4, two types classification
%
% output:
%   sleep_stage:  = 1: Wake,
%                   2: REM sleep,
%                   3: NREM light sleep,
%                   4: NREM deep sleep, when class = 1
%   sleep_stage:  = 1: NREM sleep,
%                   2: Wake+REM, when class = 4

classesn=[4 3 3 2 2];
numClasses=classesn(class);

% best parameters
filterSize=6;
balanced=1;
selectn = 3;

% 70 (10, sdnn), 73(14,avgsqi), 78(20, lfhf), 80(ac),81(dc),85(28, sampEn)
selects={[1:86],[1:69 71 72 74:77 79 82:84 86],[51:86],[51:69 71 72 74:77 79 82:84 86]};
used=selects{selectn};

load(['SHHS_ECG_model_' num2str(filterSize) '_' num2str(class) '_' num2str(balanced) '_' num2str(selectn) '_for_prediction'],'C','ka','parameters','hyperparameters','features_mean','features_std');

% calculate features

features=ECG_sleep_staging_features_extraction(ECGdata,fs);

if ~isempty(features)
    doTraining=false;

    dlXTest=features;

    for l=1:length(features_mean)
            dlXTest(l,:)=bsxfun(@minus, dlXTest(l,:),features_mean(l));
            dlXTest(l,:)=bsxfun(@rdivide, dlXTest(l,:),features_std(l));
    end

    dlXTest=dlarray(dlXTest(used,:));
%     dlYTest = dlarray(s.Ydata_all);

    dimension = size(dlXTest,1);

    dlXTest=reshape(dlXTest,dimension,1,[]);
    dlYPred = model(dlXTest,parameters,hyperparameters,doTraining);
    dlYPred = softmax(dlYPred,'DataFormat','CBT');

%     YPred = onehotdecode(dlYPred,classes,1);
%     YTest = onehotdecode(dlYTest,classes,1);
%
%     if isempty(predictions2)
%         predictions2=YPred(:);
%     else
%         predictions2 = [predictions2(:); YPred(:)];
%     end
%     yy=YPred == YTest;
%     predCorr2 = [predCorr2(:); yy(:)];

    [~,YPred]=max(dlYPred);
    YPred=extractdata(YPred);
    sleep_stage=YPred(:);

%     YTest=extractdata(dlYTest);
%     YTest=YTest(:);
%     YTest=YTest+1;
%
%     numFilters
%     C = confusionmat(YPred,YTest,'order',[numClasses:-1:1])
%     ka = kappa(C)
else
    sleep_stage = [];
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
%This folder contains the so-called CONTROL functions for bbci_apply.
%These functions transform the classifier output into the control signal
%(PACKET), that will be sent to the application via UDP. The PACKET is
%formatted as a variable/value list in a CELL.
%
%These functions have the following format:
%
% --- --- --- ---
%BBCI_CONTROL_XYZ - Generate control signal for application XYZ
%
%Synopsis:
%  [PACKET, STATE]= bbci_control_XYZ(CFY_OUT, STATE, EVENT_INFO, <PARAMS>)
%
%Arguments:
%  CFY_OUT - Output of the classifier
%  STATE - Internal state variable, which is empty in the first call of a
%      run of bbci_apply.
%  EVENT_INFO - Structure that specifies the event (fields 'time' and 'desc')
%      that triggered the evaluation of this control. Furthermore, EVENT_INFO
%      contains the whole marker queue in the field 'all_markers'.
%  PARAMS - Additional parameters that can be specified in bbci.control.param
%      are passed as further arguments.
%
%Output:
% PACKET: Variable/value list in a CELL defining the control signal that
%     is to be sent via UDP to the application.
% STATE: Updated internal state variable
% --- --- --- ---
%
%
%List of CONTROL functions (prefix bbci_control_ is left out):
% ERPSpeller: ERP-based Hex-o-Spell, one output for each complete sequence
% ERPSpellerBinary: ERP-based Hex-o-Spell, one outputs for each stimulus

author: bbci


text:
% Digital circuit sizing (vectorized)
% Boyd, Kim, Vandenberghe, and Hassibi, "A Tutorial on Geometric Programming"
% Written for CVX by Almir Mutapcic 02/08/06
% (a figure is generated)
%
% Solves the problem of choosing gate scale factors x_i to give
% minimum ckt delay, subject to limits on the total area and power.
%
%   minimize   D
%       s.t.   P <= Pmax, A <= Amax
%              x >= 1
%
% where variables are scale factors x.
%
% This code uses matrices in order to evaluate signal paths
% through the circuit (thus, it uses vectorize Matlab features).
% It is specific to the digital circuit shown in figure 4 (page 28)
% of GP tutorial paper.

% digital circuit shown in figure 4 (page 28) of GP tutorial paper
m = 7;  % number of cells
n = 8;  % number of edges
A = sparse(m,n);

% A is standard cell-edge incidence matrix of the circuit
% A_ij = 1 if edge j comes out of cell i, -1 if it comes in, 0 otherwise
  A(1,1) =     1;
  A(2,2) =     1;
  A(2,3) =     1;
  A(3,4) =     1;
  A(3,8) =     1;
  A(4,1) =    -1;
  A(4,2) =    -1;
  A(4,5) =     1;
  A(4,6) =     1;
  A(5,3) =    -1;
  A(5,4) =    -1;
  A(5,7) =     1;
  A(6,5) =    -1;
  A(7,6) =    -1;
  A(7,7) =    -1;
  A(7,8) =    -1;

% decompose A into edge outgoing and edge-incoming part
Aout = double(A > 0);
Ain = double(A < 0);

% problem constants
f = [1 0.8 1 0.7 0.7 0.5 0.5]';
e = [1 2 1 1.5 1.5 1 2]';
Cout6 = 10;
Cout7 = 10;

a     = ones(m,1);
alpha = ones(m,1);
beta  = ones(m,1);
gamma = ones(m,1);

% varying parameters for an optimal trade-off curve
N = 20;
Pmax = linspace(10,100,N);
Amax = [25 50 100];
min_delay = zeros(length(Amax),N);

disp('Generating the optimal tradeoff curve...')

for k = 1:length(Amax)
    fprintf( 'Amax = %d:\n', Amax(k) );
    for n = 1:N
        fprintf( '    Pmax = %6.2f: ', Pmax(n) );
        cvx_begin gp quiet
          % optimization variables
          variable x(m)                 % scale factors
          variable t(m)                 % arrival times

          % objective is the upper bound on the overall delay
          % and that is the max of arrival times for output gates 6 and 7
          minimize( max( t(6),t(7) ) )
          subject to
            % input capacitance is an affine function of sizes
            cin = alpha + beta.*x;

            % load capacitance is the input capacitance times the fan-out matrix
            % given by Fout = Aout*Ain'
            cload = (Aout*Ain')*cin;
            cload(6) = Cout6;          % load capacitance of the output gate 6
            cload(7) = Cout7;          % load capacitance of othe utput gate 7

            % delay is the product of its driving resistance R = gamma./x and cload
            d = cload.*gamma./x;

            % power and area definitions
            power = (f.*e)'*x;
            area = a'*x;

            % scale size, power, and area constraints
            x >= 1;
            power <= Pmax(n);
            area <= Amax(k);

            % create timing constraints
            % these constraints enforce t_j + d_j <= t_i over all gates j that drive gate i
            Aout'*t + Ain'*d <= Ain'*t;

            % for gates with inputs not connected to other gates we enforce d_i <= t_i
            d(1:3) <= t(1:3);
        cvx_end
        fprintf( 'delay = %3.2f\n', cvx_optval );
        min_delay(k,n) = cvx_optval;
    end
end

% plot the tradeoff curve
plot(Pmax,min_delay(1,:), Pmax,min_delay(2,:), Pmax,min_delay(3,:));
xlabel('Pmax'); ylabel('Dmin');
disp('Optimal tradeoff curve plotted.')

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [ x_est, x_cov ] = measurement_update(mem_x_pre, mem_cov_pre, meas, ...
    f_func_g, f_jacobian, f_hessian, meas_noise_cov, multi_noise_cov)

% MEASUREMENT_UPDATE: estimate state using second order extended kalman filter (SOEKF)
% Input:
%       mem_x_pre:      previous estimate, [m1, m2, alpha, l1, l2, (velo1, velo2)]
%                       5x1 for 'static', 7x1 for 'NCV'
%       mem_cov_pre:    previous covariance matrix
%                       5x5 for 'static', 7x7 for 'NCV'
%       meas:           measurement, 2x1
%       f_func_g:         quadratic function handle, output of GET_JACOBIAN_HESSIAN function
%       f_jacobian:       handle of Jacobian of func_g, output of GET_JACOBIAN_HESSIAN function
%       f_hessian:        handle of Hessians of func_g, output of GET_JACOBIAN_HESSIAN function
%       meas_noise_cov: covariance of measurement noise, 2x2
%       multi_noise_cov: covariance of multiplicative noise, 2x2, diag(h1_var, h2_var)
% Output:
%       x_est: estimated state
%       x_cov: covariance of state estimate

dim_state = numel(mem_x_pre);
nr_param = 5; % 5 paramtets: m1, m2, alpha, l1, l2

%% shift center to improve robustness
shifted = mem_x_pre(1:2);
meas = meas - shifted;
mem_x_pre_shifted = mem_x_pre;
mem_x_pre_shifted(1:2)=[0;0];

%% augment state and its covariance
X_est = [mem_x_pre_shifted;0;0;0;0];
X_cov = blkdiag(mem_cov_pre, multi_noise_cov, meas_noise_cov);

%% construct pseudo-measurement
z = [meas(1);meas(2);meas(1)^2;meas(2)^2; meas(1)*meas(2)];

%% Substitute quadratic function, Jacobian and Hessian matrices using current estimate
val_subs = [mem_x_pre_shifted(1:nr_param);0;0;0;0]; % velocity does not appear in measurement equation

subs_func_g = f_func_g(val_subs');
subs_jacobian = f_jacobian(val_subs');
subs_hessian = f_hessian(val_subs');

%% calculate variance of pseudo-measurement
cov_zz = zeros(nr_param, nr_param);
for i = 1:nr_param
    for j = 1:nr_param
        cov_zz(i, j) = subs_jacobian(i, :)*X_cov*subs_jacobian(j, :)'+ ...
            (1/2)*trace(subs_hessian(:, :, i)*X_cov*subs_hessian(:, :, j)*X_cov);

    end
end

%% calculate mean of pseudo-measurement
mean_z = zeros(1, nr_param);
for i = 1:nr_param
    mean_z(i) = subs_func_g(i) + (1/2)*trace(subs_hessian(:, :, i)*X_cov);
end

%% calculate cross-variance of pseudo-measurement with augmented estimate
cov_Xz = X_cov*subs_jacobian';

%% Kalman filter update
X_est = double(X_est + cov_Xz*cov_zz^-1 *( z - mean_z'));
X_cov = double(X_cov - cov_Xz*cov_zz^-1*cov_Xz');

X_est(1:2) = X_est(1:2) + shifted;

%% Truncate augmented estimate and covariance
        X_cov = (X_cov+X_cov')/2; % enforce covariance matrix symmetric
        x_est = X_est(1:dim_state);
        x_cov = X_cov(1:dim_state, 1:dim_state);
end

author: Fusion-Goettingen

repo: ExtendedObjectTracking


text:
function str= util_toString(var, varargin)
% UTIL_TOSTRING - Convert value of a variable into an (evaluatable) string
%
%Synopsis:
%  STR= util_toString(VAR, <OPT>
%
%Arguments:
%  VAR  -  Variable which is to be converted into a string
%  OPT  -  Struct or property/value list of optional properties:
%    'NumericFormat' [CHAR '%f']  Format string for printing numeric values
%    'MaxDim'        [INT  inf]  Print arrays up to this number of dimensions
%    'MAXNumel'      [INT  inf]  Print arrays up to this number of elements

props= {'NumericFormat'   '%f'    'CHAR'
        'MaxDim'          inf     'INT'
        'MaxNumel'        inf     'INT'
        'inrecursion'     0       'BOOL'
       };

if nargin==0,
  str= props; return
end

opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props);
if ~opt.inrecursion,
  opt_checkProplist(opt, props);
  opt.inrecursion= 1;
end

if numel(var)>opt.MaxNumel,
  str= sprintf('[%s %s]', str_vec2str(size(var), '%d', 'x'), ...
               upper(class(var)));
  return;
end

str= [];
if islogical(var)
  if prod(size(var))==1
    if var==0, str = 'false'; else str = 'true';end
  elseif isempty(var)
    str = '[]';
  elseif ndims(var)<=2 && size(var,2)==1
    str = sprintf('%s''',util_toString(var', opt));
  elseif (ndims(var)<=2 && size(var,1)==1) || ...
        (ndims(var)<=2 && opt.MaxDim>=2),
    str = '[';
    for i = 1:size(var,1)
      for j = 1:size(var,2)
%        str = [str , util_toString(var(i,j), opt),','];
        str = [str , util_toString(var(i,j), opt),' '];
      end
      str = [str(1:end-1),';'];
    end
    str = [str(1:end-1),']'];
  elseif ndims(var)<=opt.MaxDim,
    str = '[';
    nd = ndims(var);
    sv = size(var); sv(end)=[];
    va = permute(var,[nd,1:nd-1]);
    for i = 1:size(var,nd)
      v = va(i,:)';
      v = reshape(v,sv);
      str = [str util_toString(v, opt) ';'];
    end
    str = [str(1:end-1),']'];
  end
elseif isnumeric(var)
  if prod(size(var))==1
    if var==round(var),
      str = sprintf('%d',var);
    else
      str = sprintf(opt.NumericFormat, var);
    end
  elseif isempty(var)
    str = '[]';
  elseif ndims(var)<=2 && size(var,2)==1
    str = sprintf('%s''',util_toString(var', opt));
  elseif (ndims(var)<=2 && size(var,1)==1) || ...
        (ndims(var)<=2 && opt.MaxDim>=2),
    str = '[';
    for i = 1:size(var,1)
      for j = 1:size(var,2)
%        str = [str , util_toString(var(i,j), opt),','];
        str = [str , util_toString(var(i,j), opt),' '];
      end
      str = [str(1:end-1),';'];
    end
    str = [str(1:end-1),']'];
  elseif ndims(var)<=opt.MaxDim,
    str = '[';
    nd = ndims(var);
    sv = size(var); sv(end)=[];
    va = permute(var,[nd,1:nd-1]);
    for i = 1:size(var,nd)
      v = va(i,:)';
      v = reshape(v,sv);
      str = [str util_toString(v, opt) ';'];
    end
    str = [str(1:end-1),']'];
  end
elseif ischar(var)
  if ndims(var)<=2 & size(var,1)==1
    str = ['''',var,''''];
  elseif isempty(var)
    str = '''''';
  elseif ndims(var)<=2
    str = '[';
    for i = 1:size(var,1)
      str = [str util_toString(var(i,:), opt), ';'];
    end
    str = [str(1:end-1),']'];
  else
    str = '[';
    nd = ndims(var);
    sv = size(var); sv(end)=[];
    va = permute(var,[nd,1:nd-1]);
    for i = 1:size(var,nd)
      v = va(i,:)';
      v = reshape(v,sv);
      str = [str util_toString(v, opt) ';'];
    end
    str = [str(1:end-1),']'];
  end
elseif iscell(var)
   if prod(size(var))==1
    str = ['{',util_toString(var{1}, opt),'}'];
  elseif isempty(var)
    str = '{}';
  elseif ndims(var)<=2 & size(var,2)==1
    str = sprintf('%s''',util_toString(var', opt));
  elseif (ndims(var)<=2 && size(var,1)==1) || ...
        (ndims(var)<=2 && opt.MaxDim>=2),
    str = '{';
    for i = 1:size(var,1)
      for j = 1:size(var,2)
        str = [str , util_toString(var{i,j}, opt),','];
      end
      str = [str(1:end-1),';'];
    end
    str = [str(1:end-1),'}'];
  elseif ndims(var)<=opt.MaxDim,
    str = '{';
    nd = ndims(var);
    sv = size(var); sv(end)=[];
    va = permute(var,[nd,1:nd-1]);
    for i = 1:size(var,nd)
      v = va(i,:)';
      v = reshape(v,sv);
      str = [str util_toString(v, opt) ';'];
    end
    str = [str(1:end-1),'}'];
  end
elseif isstruct(var)
  if prod(size(var)) == 1
    str = 'struct(';
    a = fieldnames(var);
    for i = 1:length(a)
      str = [str, '''', a{i}, ''',', util_toString(getfield(var,a{i}), opt), ...
	     ','];
    end
    str = [str(1:end-~isempty(a)),')'];
  elseif isempty(var)
    str = 'struct([])';
  elseif ndims(var)<=2 && size(var,2)==1
    str = sprintf('%s''',util_toString(var', opt));
  elseif (ndims(var)<=2 && size(var,1)==1) || ...
        (ndims(var)<=2 && opt.MaxDim>=2),
    str = '[';
    for i = 1:size(var,1)
      for j = 1:size(var,2)
        str = [str , util_toString(var(i,j), opt),','];
      end
      str = [str(1:end-1),';'];
    end
    str = [str(1:end-1),']'];
  elseif ndims(var)<=opt.MaxDim,
    str = '[';
    nd = ndims(var);
    sv = size(var); sv(end)=[];
    va = permute(var,[nd,1:nd-1]);
    for i = 1:size(var,nd)
      v = va(i,:)';
      v = reshape(v,sv);
      str = [str util_toString(v, opt) ';'];
    end
    str = [str(1:end-1),']'];
  end
elseif isa(var, 'function_handle'),
  str= ['@' func2str(var)];
end

if isempty(str),
  str= sprintf('[%s %s]', str_vec2str(size(var), '%d', 'x'), ...
               upper(class(var)));
end


text:
function yourReturns = fName(arguments)
%Initialize a new vista function
%
%  fName(arguments)
%
%
% INPUTS
%
% RETURNS
%
% Web Resources
%
% Example:
%   vistaEditFunction('writeMe');
%
% Copyright Stanford team, mrVista, 2011

%% Check inputs

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%filename = 'CantileverSquareSYmmetricMesh';
%filename = 'CantileverSquareMedium';
%'CantileverSquareSYmmetricMesh';
%filename = 'CantileverSquare';
filename = 'CantileverSquareSmall';
%filename = 'CantileverSquareNew';
%'CantileverSquareNewFine';
%filename = 'Cantilever_quad_coarse';
%filename = 'LshapeTriFine';
%filename = 'LshapeTri';
%filename = 'LshapeTriSmall';
%filename = 'Lshape';
%filename = 'LshapeFine';
%filename = 'ArchTriFine';
%'Arch_quad_coarse';
%'Bridge_quad_coarse';
%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';
%filename = 'Bridge';

ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
cost = {'compliance'};
%cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%cost = {'stressNorm','compliance'};
%weights = [1,10];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'PDE';
%filterType = 'P1';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-5;
constr_initial = 1e-5;

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 16;
%
optimizer = 'DualNestedInPrimal';
%optimizer = 'AlternatingPrimalDual';
optimizerUnconstrained = 'PROJECTED GRADIENT';

%optimizer = 'MMA';
%optimizer = 'IPOPT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
%vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
vademecumFileName = 'SuperEllipseQOptAnalytic';

% designVariable = 'LevelSet';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';
% initial_case = 'full';% optimizerUnconstrained = 'SLERP';

% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';
% initial_case = 'full';
% rho0 = 0.3;

line_search_initiator = 'INCREASING LAST STEP';
incrementFactor = 1.05;
%

%kfrac = 2;
nsteps = 1;%17;

plotting = true;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 800;

%
% % %
% isDirichletPartX = @(x) x > -1e-12 & x < 0.1;
% isDirichletPart1 = @(y) y > 0.20 & y < 0.30;
% isDirichletPart2 = @(y) y > 0.70 & y < 0.80;
% isDirichletPartY = @(y) isDirichletPart1(y) | isDirichletPart2(y);
% isDirichletPart = @(x,y) isDirichletPartX(x) & isDirichletPartY(y);
%  isNeumannPartX = @(x) x > (0.96) & x < (1+1e-12);
%  isNeumannPartY = @(y) y > 0.35 & y < 0.5;
%
%  isCornerX = @(x) x > (0.4 - 0.05) & x < (0.4 +0.05);
%  isCornerY = @(y) y > (0.4 - 0.05) & y < (0.4 +0.05);
%  isCornerPart = @(x,y) isCornerX(x) & isCornerY(y);
% %
%  isNeumannPart = @(x,y) isNeumannPartX(x) & isNeumannPartY(y);
% % iNotOptimizable = @(coord) isDirichletPart(coord(:,1),coord(:,2)) & ~isNeumannPart(coord(:,1),coord(:,2));
%
% iNotOptimizable = @(coord) isCornerPart(coord(:,1),coord(:,2)) | isNeumannPart(coord(:,1),coord(:,2));
%
% costDomainNotOptimizable       = {iNotOptimizable};
% constraintDomainNotOptimizable = {[]};
%
% isDesignVariableFixed.nodes  = iNotOptimizable;
% isDesignVariableFixed.values = @(x) [m1*ones(size(x,1),1);m2*ones(size(x,1),1)];

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%> @brief Convert body struct to a spatial_v2 model
%> The purose of this is to facility an ease of conversion between
%> spatial_v2 and RBDL via files read in by yaml
function [model] = body_struct_to_spatial_model(robot, gravity)

if nargin < 2
    gravity = -9.81;
end

% dim = body_struct.dim; % @todo Add in option for planar
bodies = robot.Joints;
NB = length(bodies);
blankc = cell(NB, 1);
blankm = zeros(NB, 1);

model = struct();
% Fields for spatial_v2
model.type = 'spatial';
model.NB = NB;
model.parent = blankm;
model.jtype = blankc;
model.Xtree = blankc;
model.I = blankc;
model.mass = blankm;
model.gravity = [0, 0, gravity];

% Additional fields for kinematics
model.name = blankc;
model.axis_index = blankm;
model.is_revolute = blankm;
axes = 'xyz';

for i = 1:NB
    body = bodies(i);
    reference = body.Reference;
    lambda = getJointIndices(robot, reference.Name);

    model.parent(i) = lambda;

    index = find(body.Axis);
    assert(isscalar(index));
    axis_value = body.Axis(index);

    if strcmpi(body.Type,'revolute') || strcmpi(body.Type,'continuous')
        model.is_revolute(i) = true;
        prefix = 'R';
    else
        prefix = 'P';
    end
    if axis_value < 0
        jtype = ['-', prefix, axes(index)];
    else
        jtype = [prefix, axes(index)];
    end

    model.jtype{i} = jtype;
    model.Xtree{i} = plux(body.R, body.Offset);

    child_link_idx = getLinkIndices(robot, body.Child);

    if ~isnan(child_link_idx)
        link = robot.Links(child_link_idx);
        model.I{i} = mcI(link.Mass, link.Offset, link.Inertia);
        model.mass(i) = link.Mass;
    else
        model.I{i} = mcI(0, zeros(3,1), zeros(3));
        model.mass(i) = 0;
    end
    model.name{i} = body.Name;
    model.axis_index(i) = index;
end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [P,G] = fs_signrank(x_fname,y_fname,varargin)
%
%  Author: Martin Reuter
%
%  Extends Wilcoxon signed rank test for zero median to FS surfaces
%  If y_fname is passed, do a paired test, see fs_ranksum for independent
%  samples.
%
%  Can work in parallel mode (run 'matlabpool OPEN' first)
%
%  x_fname    file name of mgh file (stacked surfaces)
%  y_fname    optional: file name for second file to do a paired test
%
%  fs_signrank(...,'mask',MASK) selects a mask (e.g. cortex.label)
%
%  fs_signrank(...,'outsig',OUTSIG) specifies file to output p-values
%      output will be in FreeSurfer mgh format, signed - log10(p)
%
%  fs_signrank(...,'outgamma',OUTGAMMA specifies file to output gamma
%      output will be in FreeSurfer mgh format (see G below)
%
%  See signrank for more arguments ('method' and 'tail')
%
%  Return (only inside mask):
%
%  G          gamma: median(X)            for single sample
%                    median(X)-median(Y)  for paired samples
%
%  P          p-values of signed rank test
%

% parse and remove mask and outfile params:
oknames = {'mask' 'outsig' 'outgamma' 'tail' 'alpha' 'method'};
dflts   = {'' '' '' 'both' 0.05   '' };
[mask_fname , outsig, outgamma, tail] = internal.stats.parseArgs(oknames,dflts,varargin{:});
% remove mask and outfile:
del = zeros(size(varargin));
for i = 1:length(varargin)
  if (strcmp(varargin{i},'mask'))
    del(i) = 1;
    del(i+1) = 1;
  end
  if strcmp(varargin{i},'outsig')
    del(i) = 1;
    del(i+1) = 1;
  end
  if strcmp(varargin{i},'outgamma')
    del(i) = 1;
    del(i+1) = 1;
  end
end
varargin = varargin(find(del == 0));
%tail = getParamVal(tail,{'both' 'right' 'left'},'''tail''');
if (strcmp(tail,'both'))
  tail = 0;
elseif (strcmp(tail,'left'))
  tail = -1;
elseif (strcmp(tail,'right'))
  tail = 1;
else
  error('tail can only be "both","left" or "right"');
end

% read X from file:
if ( ischar(x_fname) && exist(x_fname,'file') )
  [X,Mx] = fs_read_Y(x_fname);
  nsubj = size(X,1);
  nvert = size(X,2);
else
  error(['Can not load ' x_fname ' as an mgh or mgz file']);
end

% if passed read Y from file:
Y = [];
if (  ~ isempty(y_fname) )
  if ischar(y_fname) && exist(y_fname,'file')
    [Y,My] = fs_read_Y(y_fname);
    if (size(X) ~= size(Y))
      error('Dimensions of X and Y do not agree');
    end
  else
      error(['Can not load ' y_fname ' as an mgh or mgz file']);
  end
end

% read mask from file:
mask = [1:nvert];
outside = [];
if ( ~ isempty(mask_fname) )
  mask = fs_read_label(mask_fname);
  zz = zeros(1,nvert);
  zz(mask) = ones(1,length(mask));
  outside = find(zz==0);
end

% process
n = length(mask);
Pt = ones(1,n);
if (isempty(Y))
  parfor i = 1:n
    Pt(i) = signrank(X(:,mask(i)),'',varargin{:});
  end
  G = median(X);
else
  parfor i = 1:n
    idx = mask(i);
    Pt(i) = signrank(X(:,idx),Y(:,idx),varargin{:});
  end
  G = median(X)-median(Y);
end
P = ones(1,nvert);
P (mask) = Pt;
G(outside) = zeros(1,length(outside));

% write (on full surface)
if ( ~ isempty(outsig) )
  Pl10 = -log10(P) .* sign(G);
  Mx.volsz(4) = 1;
  fs_write_Y(Pl10,Mx,outsig);
end

if ( ~ isempty(outgamma) )
  Mx.volsz(4) = 1;
  fs_write_Y(G,Mx,outgamma);
end

% return values (only inside mask):
P = Pt;
G = G(mask);
meanG = mean(G);


text:
function [um] = ft2um(ft)
% Convert length from feet to micrometers.
% Chad Greene 2012
um = ft*304800;


text:
classdef SuperpixelLSC < handle
    %SUPERPIXELLSC  Class implementing the LSC (Linear Spectral Clustering) superpixels algorithm
    %
    % As described in [LiCVPR2015LSC].
    %
    % LSC (Linear Spectral Clustering) produces compact and uniform
    % superpixels with low computational costs. Basically, a normalized cuts
    % formulation of the superpixel segmentation is adopted based on a
    % similarity metric that measures the color similarity and space proximity
    % between image pixels. LSC is of linear computational complexity and high
    % memory efficiency and is able to preserve global properties of images.
    %
    % ## References
    % [LiCVPR2015LSC]:
    % > Zhengqin Li and Jiansheng Chen. "Superpixel Segmentation using Linear
    % > Spectral Clustering". IEEE Conference on Computer Vision and Pattern
    % > Recognition (CVPR), June 2015.
    %
    % See also: cv.SuperpixelLSC.SuperpixelLSC, cv.SuperpixelSLIC,
    %  cv.SuperpixelSEEDS, superpixels
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = SuperpixelLSC(img, varargin)
            %SUPERPIXELLSC  Class implementing the LSC (Linear Spectral Clustering) superpixels
            %
            %     obj = cv.SuperpixelLSC(img)
            %     obj = cv.SuperpixelLSC(img, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ Image to segment.
            %
            % ## Options
            % * __RegionSize__ Chooses an average superpixel size measured in
            %   pixels. default 10
            % * __Ratio__ Chooses the enforcement of superpixel compactness
            %   factor of superpixel. default 0.075
            %
            % The function initializes a SuperpixelLSC object for the input
            % image. It sets the parameters of superpixel algorithm, which
            % are: `RegionSize` and `Ratio`. It preallocate some buffers for
            % future computing iterations over the given image. An example of
            % LSC is ilustrated in the following picture.
            % For enanched results it is recommended for color images to
            % preprocess image with little gaussian blur with a small 3x3
            % kernel and additional conversion into CieLAB color space.
            %
            % ![image](https://docs.opencv.org/3.3.1/superpixels_lsc.png)
            %
            % See also: cv.SuperpixelLSC.iterate
            %
            this.id = SuperpixelLSC_(0, 'new', img, varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.SuperpixelLSC
            %
            if isempty(this.id), return; end
            SuperpixelLSC_(this.id, 'delete');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.SuperpixelLSC.empty, cv.SuperpixelLSC.load
            %
            SuperpixelLSC_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.SuperpixelLSC.clear, cv.SuperpixelLSC.load
            %
            b = SuperpixelLSC_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.SuperpixelLSC.load
            %
            SuperpixelLSC_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.SuperpixelLSC.save
            %
            SuperpixelLSC_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.SuperpixelLSC.save, cv.SuperpixelLSC.load
            %
            name = SuperpixelLSC_(this.id, 'getDefaultName');
        end
    end

    %% SuperpixelLSC
    methods
        function num = getNumberOfSuperpixels(this)
            %GETNUMBEROFSUPERPIXELS  Calculates the actual amount of superpixels on a given segmentation computed and stored in object
            %
            %     num = obj.getNumberOfSuperpixels()
            %
            % ## Output
            % * __num__ actual amount of superpixels.
            %
            % See also: cv.SuperpixelLSC.iterate
            %
            num = SuperpixelLSC_(this.id, 'getNumberOfSuperpixels');
        end

        function iterate(this, varargin)
            %ITERATE  Calculates the superpixel segmentation on a given image with the initialized parameters in the object
            %
            %     obj.iterate()
            %     obj.iterate('OptionName',optionValue, ...)
            %
            % ## Options
            % * __NumIterations__ Number of iterations. Higher number improves
            %   the result. default 10
            %
            % The function computes the superpixels segmentation of an image
            % with the parameters initialized with the constructor. The
            % algorithms starts from a grid of superpixels and then refines
            % the boundaries by proposing updates of edges boundaries.
            %
            % This function can be called again without the need of
            % initializing the algorithm with the constructor This save the
            % computational cost of allocating memory for all the structures
            % of the algorithm.
            %
            % See also: cv.SuperpixelLSC.getLabels
            %
            SuperpixelLSC_(this.id, 'iterate', varargin{:});
        end

        function labels = getLabels(this)
            %GETLABELS  Returns the segmentation labeling of the image
            %
            %     labels = obj.getLabels()
            %
            % ## Output
            % * __labels__ Return a `int32` integer array containing the
            %   labels of the superpixel segmentation. The labels are in the
            %   range `[0, obj.getNumberOfSuperpixels()]`.
            %
            % The function returns an image with the labels of the superpixel
            % segmentation. The labels are in the range
            % `[0, obj.getNumberOfSuperpixels()]`.
            %
            % Each label represents a superpixel, and each pixel is assigned
            % to one superpixel label.
            %
            % See also: cv.SuperpixelLSC.iterate
            %
            labels = SuperpixelLSC_(this.id, 'getLabels');
        end

        function img = getLabelContourMask(this, varargin)
            %GETLABELCONTOURMASK  Returns the mask of the superpixel segmentation stored in object
            %
            %     img = obj.getLabelContourMask()
            %     img = obj.getLabelContourMask('OptionName',optionValue, ...)
            %
            % ## Output
            % * __img__ Return `logical` image mask where 1 indicates that the
            %   pixel is a superpixel border, and 0 otherwise.
            %
            % ## Options
            % * __ThickLine__ If false, the border is only one pixel wide,
            %   otherwise all pixels at the border are masked. default true
            %
            % The function return the boundaries of the superpixel
            % segmentation.
            %
            % See also: cv.SuperpixelLSC.iterate, boundarymask
            %
            img = SuperpixelLSC_(this.id, 'getLabelContourMask', varargin{:});
            img = (img == 255);  % fg:uint8(255), bg:uint8(0)
        end

        function enforceLabelConnectivity(this, varargin)
            %ENFORCELABELCONNECTIVITY  Enforce label connectivity
            %
            %     obj.enforceLabelConnectivity()
            %     obj.enforceLabelConnectivity('OptionName',optionValue, ...)
            %
            % ## Options
            % * __MinElementSize__ The minimum element size in percents that
            %   should be absorbed into a bigger superpixel. Given resulted
            %   average superpixel size valid value should be in 0-100 range,
            %   25 means that less then a quarter sized superpixel should be
            %   absorbed, this is default. default 20
            %
            % The function merge component that is too small, assigning the
            % previously found adjacent label to this component. Calling this
            % function may change the final number of superpixels.
            %
            % See also: cv.SuperpixelLSC.iterate
            %
            SuperpixelLSC_(this.id, 'enforceLabelConnectivity', varargin{:});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv