text:
classdef EvolveProcessCapture
    % Shows the evolution process of spectral unmixing algorithms.
    properties
      ScatterOpts
      AbundOpts
      FramesScatter
      FramesAbund
    end
    methods
        function obj = EvolveProcessCapture()
        end

        function obj = begin(obj,max_iter)
            scatter_opts = [];
            scatter_opts.hs = [];
            scatter_opts.update = 0;

            abund_opts = [];
            abund_opts.imhs = [];
            abund_opts.update = 0;

            frames_scatter(max_iter) = struct('cdata',[],'colormap',[]);
            frames_abund(max_iter) = struct('cdata',[],'colormap',[]);

            obj.ScatterOpts = scatter_opts;
            obj.AbundOpts = abund_opts;
            obj.FramesScatter = frames_scatter;
            obj.FramesAbund = frames_abund;
        end

        function obj = end(obj,end_iter)
            obj.FramesScatter = obj.FramesScatter(1:end_iter);
            obj.FramesAbund = obj.FramesAbund(1:end_iter);
        end

        function obj = updateEndmVar(obj,Y,A,rows,cols,w_jk,mu_jk,sigma_jk,names,iter)
            obj = updateScatterPlotEndmVar(obj,Y,w_jk,mu_jk,sigma_jk,names,iter);
            obj = updateAbundanceMaps(obj,A,rows,cols,iter);

            drawnow
        end

        function obj = updateEndmFixed(obj,Y,A,rows,cols,R,names,iter)
            obj = updateScatterPlotEndmFixed(obj,Y,R,names,iter);
            obj = updateAbundanceMaps(obj,A,rows,cols,iter);

            drawnow
        end

        function obj = updateScatterPlotEndmVar(obj,Y,w_jk,mu_jk,sigma_jk,names,iter)
            scatter_opts = obj.ScatterOpts;
            frames_scatter = obj.FramesScatter;

            [~,hs,axes_h,scatter_fh] = endmember_scatter_plot_end_var( ...
                Y,w_jk,mu_jk,sigma_jk,names,scatter_opts);
            title(axes_h, ['Scatter plot of iteration ',num2str(iter)]);
            frames_scatter(iter) = getframe(scatter_fh);

            if ~scatter_opts.update
                pos = get(scatter_fh,'Position');
                pos(1) = pos(1) - pos(3) / 2 - 1;
                set(scatter_fh,'Position',pos);
            end
            scatter_opts.update = 1;
            scatter_opts.hs = hs;
            scatter_opts.axes_h = axes_h;
            scatter_opts.fig_h = scatter_fh;

            obj.ScatterOpts = scatter_opts;
            obj.FramesScatter = frames_scatter;
        end

        function obj = updateScatterPlotEndmFixed(obj,Y,R,names,iter)
            scatter_opts = obj.ScatterOpts;
            frames_scatter = obj.FramesScatter;

            [~,hs,axes_h,scatter_fh] = endmember_scatter_plot( ...
                Y,R,names,scatter_opts);
            title(axes_h, ['Scatter plot of iteration ',num2str(iter)]);
            frames_scatter(iter) = getframe(scatter_fh);

            if ~scatter_opts.update
                pos = get(scatter_fh,'Position');
                pos(1) = pos(1) - pos(3) / 2 - 1;
                set(scatter_fh,'Position',pos);
            end
            scatter_opts.update = 1;
            scatter_opts.hs = hs;
            scatter_opts.axes_h = axes_h;
            scatter_opts.fig_h = scatter_fh;

            obj.ScatterOpts = scatter_opts;
            obj.FramesScatter = frames_scatter;
        end

        function obj = updateAbundanceMaps(obj,A,rows,cols,iter)
            abund_opts = obj.AbundOpts;
            frames_abund = obj.FramesAbund;

            [imhs,abund_fh] = show_abundances(A,rows,cols,'',0,0,abund_opts);
            set(abund_fh, 'name', ['Abundances of iteration ',num2str(iter)]);
            frames_abund(iter) = getframe(abund_fh);

            if ~abund_opts.update
                pos = get(abund_fh,'Position');
                pos(1) = pos(1) + pos(3) / 2 + 1;
                set(abund_fh,'Position',pos);
            end
            abund_opts.update = 1;
            abund_opts.imhs = imhs;
            abund_opts.fig_h = abund_fh;

            obj.AbundOpts = abund_opts;
            obj.FramesAbund = frames_abund;
        end
    end
end

author: zhouyuanzxcv

repo: Hyperspectral


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================

% function  [Jc,para,dJ,H] = PIRBFGSobjFctn(T,Rc,omega,m,beta,M,wRef,xc,wc)
%
% Objective Function for Parametric Image Registration with lBFGS
%
% computes J(wC) = D(T(Y(wc)),R) + S(wc), where
%
% yc       = y(wc,xc) = trafo(wc,xc)
% Tc       = T(yc) = inter(T,omega,yc),
% D(Tc,Rc) = distance(Tc,Rc,omega,m)
% S(wc)    = 0.5*(wc-wRef)'*M*(wc-wRef);
%
% Input:
%   T       coefficients for template image
%   Rc      sampled reference image
%   omega   spatial domain
%   m       number of discretization points
%   beta    adding beta*I to approximation of Hessian
%   xc      discretization of Omega
%   wc      current parameters
% Output:
%  Jc       current function value J(wc)
%  para     struct {Tc=T(y(wc)), Rc, omega, m, yc=y(wc,xc), Jc}, for plots
%  dJ       gradient of J
%  H        approximation to Hessian of J
%==============================================================================

function [Jc,para,dJ,H] = PIRBFGSobjFctn(T,Rc,omega,m,beta,M,wRef,xc,wc)

if nargin == 0,

  help(mfilename)
  runMinimalExample;
  Jc = 'endMinimalExample';
  return

elseif ~exist('wc','var') || isempty(wc),

  % if wc is not an input argument, reports status
  if nargout == 1, Jc = 'PIR';  return; end;

  % report current settings
  dimstr  = @(m) sprintf('[%s]',sprintf(' %d',m));
  wc      = trafo('w0');

  FAIRmessage('Parametric Image Registration with BFGS','-');
  fprintf('J(wc)=D(T(y(wc)),R) + (wc-wRef)''*M*(wc-wRef) != min\n');
  fprintf('%-20s : %s\n','m',dimstr(m));
  fprintf('%-20s : %s\n','omega',dimstr(omega));
  fprintf('%-20s : %s\n','IMAGE MODEL',imgModel);
  fprintf('%-20s : %s\n','DISTANCE',distance);
  fprintf('%-20s : %s\n','TRAFO',trafo);
  fprintf('%-20s : %s\n','length(wc)',num2str(length(wc)));
  fprintf('%-20s : %s\n','REGULARIZATION',...
    sprintf('M is %d-by-%d, beta=%s',size(M),num2str(beta)));
  FAIRmessage('-');

  return;
end;

% do the work ------------------------------------------------------------
doDerivative = (nargout>2);            % flag for necessity of derivatives

% compute transformation, distance, and regularization and combine these
[yc,dy]    = trafo(wc,center(xc,m),'doDerivative',doDerivative);
[Tc,dT]    = imgModel(T,omega,yc,'doDerivative',doDerivative);
[Dc,rc,dD] = distance(Tc,Rc,omega,m,'doDerivative',doDerivative);

% add regularization
if isempty(M) || (all(size(M)==1) && (M==0)),
  Sc = 0;
  dS = 0;
  M  = 0;
else
  dS = (wc-wRef)'*M;
  Sc = 0.5*dS*(wc-wRef);
end;

Jc = Dc + Sc;

% collect variables for plots
para = struct('Tc',Tc,'Rc',Rc,'omega',omega,'m',m,'yc',yc,'Jc',Jc);

if ~doDerivative, return; end;
dD = dD * dT;                    % multiply outer and inner derivatives
if size(dD,2) == size(dy,1),     % generic case, dy comes complete
  dJ = dD*dy + dS;
else                             % tricky case,  dy comes sparse
  n  = size(dy{1},1);
  dJ = reshape(dy{1}'*reshape(dD*dT,n,[]),1,[]) + dS;
%   dr = [dr(:,1:n)*dy{1},dr(:,n+1:2*n)*dy{1},dr(:,2*n+1:3*n)*dy{1}];
end;
if nargout<4, return; end;

% approximation to Hessian
if M == 0,
  H = 1;
else
  H = M;
end;

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

function runMinimalExample

setup2DhandData;

% extract data of level 4
level = 4;
omega = ML{level}.omega;
m     = ML{level}.m;

imgModel('reset','imgModel','linearInter');
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);
xc    = getCellCenteredGrid(omega,m);
Rc    = imgModel(R,omega,xc);

% initialize distance measure
distance('reset','distance','SSD');

% initialize the transformation and a starting guess
% trafo('reset','trafo','affine3Dsparse');
trafo('reset','trafo','affine2D');
w0 = trafo('w0');

%------------------------------------------------------------------------------
% build objective function
% note: T  is data for template image
%       Rc is sampled reference image
%       optional Tikhonov-regularization is disabled by setting m = [], wRef = []
%       beta = 0, M = [], wRef = []:
%       disables additional regularization of Hessian approximation
beta = 0; M = []; wRef = [];
fctn = @(wc) PIRobjFctn(T,Rc,omega,m,beta,M,wRef,xc,wc);
fctn([]);   % report status
checkDerivative(fctn,w0+rand(size(w0)));

% setup plots and initialize
FAIRplots('reset','mode','PIR-objective','omega',omega,'m',m,'fig',1,'plots',1);
FAIRplots('init',struct('Tc',T,'Rc',R,'omega',omega,'m',m));

%% -- solve the optimization problem on one level
OPTpara = FAIRcell2struct(optPara('lBFGS','solver',''));
[wc,his] = lBFGS(fctn,w0,OPTpara{:});
return;

% %------------------------------------------------------------------------------
% %% finally: run the MultiLevel Non-Parametric Image Registration
% [wc,his] = MLPIR(ML);
%==============================================================================


text:
classdef MeasurementListX < BaseX
% MeasurementListX class
%
% Summary of MeasurementListX:
% Class implementation of the primitive Measurement List type.

    properties (Dependent)
        Measurements
        NumMeasurements
    end

    properties
        Vectors
        Models
        Tags
        Timestamp
        Metadata
    end

    methods
        function this = MeasurementListX(varargin)
        % Measurement Base constructor
        %
        % Parameters
        % ----------
        % measurements: vector of MeasurementX objects

            if nargin>1
                this.Vectors = varargin{1};
                this.Timestamp = varargin{2};
                if (nargin>2)
                    this.Tags = varargin{3};
                end
            elseif nargin && ~isempty(varargin{1})
                if isa(varargin{1},'MeasurementListX')
                    this.Vectors = varargin{1}.Vectors;
                    this.Tags = varargin{1}.Tags;
                    this.Models = varargin{1}.Models;
                    this.Timestamp = varargin{1}.Timestamp;
                    this.Metadata = varargin{1}.Metadata;
%                     if(isfield(varargin{1},'Tags'))
%                         this.Tags
                elseif isa(varargin{1},'MeasurementX')
                    measurements = varargin{1};
                    this.Vectors = [measurements.Vector];
                    this.Tags = [measurements.Tag];
                    this.Models = [measurements.Model];
                    this.Metadata = [measurements.Metadata];
                    this.Timestamp = measurements.Timestamp;
                else
                    this.Vectors = varargin{1};
                    this.Timestamp = [];
                end
            else
                this.Vectors = [];
                this.Timestamp = [];
            end
        end

        function measurements = getMeasurements(this, indices)
            numMeasurements = size(this.Vectors,2);
            measurements = MeasurementX.empty();
            if any(indices<1) || any(indices>numMeasurements)
                error('Index out of range');
            end
            for k=indices
                if ~isempty(this.Models)
                    model = this.Models(k);
                else
                    model = [];
                end
                if ~isempty(this.Tags)
                    tag = this.Tags(k);
                else
                    tag = [];
                end
                if ~isempty(this.Metadata)
                    metadata = this.Metadata(k);
                else
                    metadata = struct();
                end
                measurements(end+1) = MeasurementX(this.Vectors(:,k),...
                                               this.Timestamp,...
                                               model,...
                                               tag);
                measurements(end).Metadata = metadata;
            end
        end

        function numMeasurements = get.NumMeasurements(this)
            numMeasurements = size(this.Vectors,2);
        end

        function measurements = get.Measurements(this)
            if(isempty(this.Vectors))
                measurements = MeasurementX.empty();
            else
                for k = 1 : this.NumMeasurements
                    if ~isempty(this.Models)
                        model = this.Models(k);
                    else
                        model = [];
                    end
                    if ~isempty(this.Tags)
                        tag = this.Tags(k);
                    else
                        tag = [];
                    end
                    if ~isempty(this.Metadata)
                        metadata = this.Metadata(k);
                    else
                        metadata = struct();
                    end
                    measurements(k) = MeasurementX(this.Vectors(:,k),...
                                                   this.Timestamp,...
                                                   model,...
                                                   tag);
                    measurements(k).Metadata = metadata;
                end
            end
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function H = m_tricontour(tri,p,z,N,C)
%   H = m_tricontour(tri,p,z,N,C)
%   N-> number of contour levels
%   C-> linespec color
global MAP_PROJECTION

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
    disp('No Map Projection initialized - call M_PROJ first!');
    return;
end

[X,Y] = m_ll2xy(p(:,1),p(:,2),'clip','on');
hold on;
if nargin < 5
    H = tricontour(tri,X,Y,z,N);
else
    H = tricontour(tri,X,Y,z,N,C);
end

end

author: CHLNDDEV


text:
function [dWdB, dWdD, dWdS, dWdalpha, dWdgParam, dWdn] = gpsimMapWGradient(model, ...
                                                  k)
% GPSIMMAPWGRADIENT Compute the gradients of W with respect to the parameters of the k-th gene.
% FORMAT
% DESC computes the gradients of W with respect to the parameters of the k-th gene given
% Gaussian process for use in a single input motif protein network.
% ARG model : the model for which the log likelihood is computed.
% ARG k : the k-th gene.
% RETURN dWdB : the gradients of W w.r.t the model paramter Bk.
% RETURN dWdD : the gradients of W w.r.t the model paramter Dk.
% RETURN dWdS : the gradients of W w.r.t the model paramter Sk.
% RETURN dWdalpha : the gradients of W w.r.t the model paramter alpha_k.
% RETURN dWdgParam : the gradients of W w.r.t the model paramter gamma_k
% RETURN dWdn : the gradients of W w.r.t the general noise variance of the k-gene.
%
% SEEALSO : gpsimMapCreate, gpsimMapLogLikelihood,
% gpsimMapGradient, gpsimMapFunctionalLogLikeGradients
%
% COPYRIGHT : Pei Gao, Magnus Rattray and Neil D. Lawrence, 2008

% GPSIM

intPoints = model.times_index(1)+1:(model.numMapPts);
step2 = model.step*model.step;
S2 = model.S.*model.S;
numData = length(model.t);

[w1, w2] = size(model.W);

dWdB = zeros(w1, w2);
dWdD = zeros(w1, w2);
dWdS = zeros(w1, w2);
dWdalpha = [];
dWdgParam = [];
dWdn = [];
dWdalpha = [];
if model.ngParam > 0
  ngParamk = model.ngParam/model.numGenes;
  dWdgParam= zeros(w1, w2, ngParamk);
  gInd = k;
else
  gInd = 1;
end

% check if it's multiple g(f).
if isfield(model, 'isGroupNonlinearity') && strcmp(model.nonLinearity{k}, ...
                                                   'repression')
  dWdalpha = zeros(w1, w2);
end

if isfield(model, 'includeNoise') && model.includeNoise
  noiseMat = ones(numData, 1)*model.noiseVar;
  yvar = model.yvar + noiseMat;
  dWdn = zeros(w1, w2);
else
  yvar = model.yvar;
end

for p = intPoints
    for i=1:numData
      arg = model.t(i)-model.mapt(p);
      if arg >= 0
        ind = i + (k-1)*numData;
        beta_ik=1/yvar(ind);

        [dxdB dxdD dxdS dxdalpha dxdgParam] = gpsimXGradient(model, i, k);

        dWdB(p, p)=dWdB(p, p)+beta_ik*model.g_grad2(p,gInd)*dxdB* ...
                exp(-model.D(k)*arg+log(model.S(k)) +log(model.step));

        if isfield(model, 'isGroupNonlinearity')
          if strcmp(model.nonLinearity{k}, 'repression')
            dWdalpha(p,p) = dWdalpha(p, p)+beta_ik*model.g_grad2(p,gInd)* ...
                dxdalpha*exp(-model.D(k)*arg+log(model.S(k)) +log(model.step));
          end
        end

        factor = model.ypred(model.times_index(i), k)-model.y(ind);

        dWdD(p, p) = dWdD(p, p)+model.step*beta_ik*model.g_grad2(p,gInd)* ...
                (dxdD-factor*arg)*exp(-model.D(k)*arg+log(model.S(k))) ;

        dWdS(p, p) = dWdS(p, p)+model.step*beta_ik*exp(-model.D(k)* ...
                           arg)*(dxdS*model.S(k)*model.g_grad2(p,gInd)+ ...
                                 factor*model.g_grad2(p,gInd));

        if model.ngParam > 0
          for gParamInd = 1:ngParamk
            dWdgParam(p, p, gParamInd)= dWdgParam(p, p, gParamInd)+ ...
                model.step*beta_ik*exp(-model.D(k)*arg)* ...
                (dxdgParam(gParamInd)*model.S(k)*model.g_grad2(p, ...
                gInd)+ factor*model.S(k)*model.dggrad2(p,gInd));
          end
        end

        if isfield(model,'includeNoise') && model.includeNoise
          dWdn(p, p) = dWdn(p, p)-2*sqrt(model.noiseVar(k))* ...
              beta_ik^2*factor*model.g_grad2(p,gInd)*exp(-model.D(k)*arg+ ...
                                                         log(model.step)+ ...
                                                         log(model.S(k)));
        end
      end
    end
end

for p = intPoints
  for q = intPoints
    for i = 1:numData
      arg1 = model.t(i)-model.mapt(p);
      arg2 = model.t(i)-model.mapt(q);
      if arg1 >= 0 && arg2 >= 0
        ind = i + (k-1)*numData;
        beta_ik = 1/yvar(ind);

        dWdD(p,q) = dWdD(p,q)-beta_ik*model.g_grad(p,gInd)* ...
            model.g_grad(q,gInd)*(arg1+arg2)*exp(-model.D(k)*(arg1+arg2)+ ...
                                                 log(S2(k))+log(step2));

        dWdS(p,q) = dWdS(p,q)+2*beta_ik*model.S(k)* model.g_grad(q,gInd)* ...
            model.g_grad(p,gInd) * exp(-model.D(k)*(arg1+arg2)+log(step2));

        if model.ngParam > 0
          for gParamInd = 1:ngParamk
            dWdgParam(p, q, gParamInd)= dWdgParam(p, q, gParamInd)+beta_ik*(model.dggrad(q,gInd)*model.g_grad(p,gInd)+model.g_grad(q,gInd)*model.dggrad(p,gInd))*exp(-model.D(k)*(arg1+arg2)+log(S2(k))+log(step2));
          end
        end

        if isfield(model, 'includeNoise') && model.includeNoise
          dWdn(p, q) = dWdn(p, q)-2*sqrt(model.noiseVar(k))* ...
              beta_ik^2*model.g_grad(p,gInd)*model.g_grad(q,gInd)*exp(- ...
                    model.D(k)*(arg1+arg2)+log(step2)+log(S2(k)));
        end
      end
    end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcFlipV(Temp,Event,handles)
handles=guidata(handles.MU_matrix_display);

if ~isempty(handles.V.Segs)
    choice = questdlg('Segmentation mask is detected, transform operation will reset them, preceed?','Mask Reset', ...
                      'No, go save mask','Yes','No, go save mask');
    if isempty(choice)
        warndlg('Transform is cancelled.');
        return;
    end
    % Handle response
    switch choice
        case 'No, go save mask'
            warndlg('Save your mask before transform.');
            return;
    end

    handles.Mask=handles.Mask*0;
    handles.V.Segs=[];
end

% close 3D slicer
global Figure_handles
if isfield(Figure_handles,'MU_display2')
    slicer_display_handles=guidata(Figure_handles.MU_display2);
    if Figure_handles.MU_display == slicer_display_handles.Parent
        close(Figure_handles.MU_display2);
    end
end

handles.TMatrix=flipdim(handles.TMatrix,1);
handles.Mask=flipdim(handles.Mask,1);

MergeM=get(handles.Matrix_name_edit,'String');
set(handles.Matrix_name_edit,'String',[MergeM '_fpv']);

% update current display matrix
handles=MU_update_image(handles.Matrix_display_axes,{handles.TMatrix,handles.Mask},handles,0);
guidata(handles.MU_matrix_display, handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function d = dirac(f, varargin)
%DIRAC    Dirac delta function.
% D = DIRAC(F) returns a CHEBFUN D which is zero on the domain of the CHEBFUN F
% except at the simple roots of F, where it is infinite.
%
% DIRAC(F, N) is the nth derivative of DIRAC(F).
%
% DIRAC(F) is not defined if F has a zero of order greater than one within the
% domain of F.
%
% If F has break-points, they should not coincide with the roots of F. However,
% F can have simple roots at either end points of its domain.
%
% See also HEAVISIDE.

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

% Empty argument:
if ( isempty(f) )
    d = f;
    return
end

% Deal with quasimatrices:
if ( numColumns(f) > 1 )
    f = cheb2cell(f);
    for k = 1:numel(f)
        f{k} = dirac(f{k}, varargin{:});
    end
    d = horzcat(f{:});
    return
end

% Handle the case for derivatives of delta function:
if ( nargin > 1 )
    if ( nargin > 2 )
        error('CHEBFUN:CHEBFUN:dirac:dirac', 'Too many input arguments.');
    end

    % Order of the derivative of dirac delta function:
    n = varargin{1};
    if ( ~isnumeric(n) || n < 0 || round(n) ~= n || ~isscalar(n) )
        error('CHEBFUN:CHEBFUN:dirac:dirac', ...
            'Order of the derivative must be be a non-negative integer.');
    end

    if ( n == 0 ) % Trivial case
        d = dirac(f);
        return
    else
        d = diff(dirac(f), n);
        return
    end
end

% Set a tolerance and get the domain of f:
tol = eps;
dom = f.domain;
a = dom(1);
b = dom(end);
vscl = vscale(f);

% Extract the 'normal' roots of f:
r = roots(f, 'nojump', 'nozerofun');
r = sort(r(:));

% Check roots at the end points of f:
if ( isempty(r) )
    % If there are no roots, still check roots at the end points:
    if ( abs(feval(f, a, 'right')) < 100*tol*vscl )
        rootA = 1;
        r = [r; a];
    else
        rootA = 0;
    end

    if ( abs(feval(f, b, 'left')) < 100*tol*vscl )
        rootB = 1;
        r = [r; b];
    else
        rootB = 0;
    end
else
    % If there are roots, check if they are at the end points:
    if ( r(1) > a )
        rootA = 0;
    elseif ( abs(feval(f, a, 'right')) < 100*tol*vscl )
        rootA = 1;
        if ( r(1) ~= a )
            r = [a ; r];
        end
    end
    if ( r(end) < b )
        rootB = 0;
    elseif ( abs(feval(f, b, 'left')) < 100*tol*vscl )
        rootB = 1;
        if ( r(end) ~= b )
            r = [r ; b];
        end
    end
end

% Initialize a zero CHEBFUN:
d = chebfun(0, [a, b]);

% If there is no root of F within the domain or at the end points, return with a
% zero CHEBFUN:
if ( isempty( r ) )
    return
end

% Check if any of the roots is not simple by looking at the derivative of F:
fp = diff(f);
fpVals = feval(fp, r);

% Check root order for interior break-points:
if ( any(abs(fpVals) < 100*tol*vscale(fp)) )
    error('CHEBFUN:CHEBFUN:dirac:dirac', ...
        'Function has a root which is not simple');
else
    % Place deltas with appropriate scaling at interior roots.
    deltaMag = 1./abs(fpVals);
end

% Use half of the strength if there is a root at the end point of the
% domain of the input CHEBFUN and update the pointValues:
pointValues = [0; 0];
if ( rootA )
    deltaMag(1) = deltaMag(1)/2;
    pointValues(1) = sign(deltaMag(1))*inf;
end
if ( rootB )
    deltaMag(end) = deltaMag(end)/2;
    pointValues(2) = sign(deltaMag(end))*inf;
end

% Call the DELTAFUN constructor directly:
data.deltaMag = deltaMag.';
data.deltaLoc = r.';
d.funs{1} = deltafun(d.funs{1}, data);
d.pointValues = pointValues;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function value = ch_indexi ( s, c )

%*****************************************************************************80
%
%% CH_INDEXI is the (case insensitive) first occurrence of a character in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be searched.
%
%    Input, character C, the character to be searched for.
%
%    Output, integer VALUE, the location of the first occurrence of C
%    (upper or lowercase), or 0 if C does not occur.
%
  value = 0;

  for i = 1 : length ( s )

    if ( ch_eqi ( s(i), c ) )
      value = i;
      return
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%--- help for ts.step_dummy ---
%
%  Implements a step dummy in the time series
%
%  Args:
%
%     start_date: start date of the time series
%
%     end_date: end date of the time series
%
%     dummy_start_date: start date of the step dummy
%
%  Returns:
%     :
%
%     - **db** (ts object): a time series with zeros from the first observation to the date
%       before the start of the dummy and ones from then on.
%
%

author: jmaih


text:
function z = mutInfo(x, y)
% Compute mutual information I(x,y) of two discrete variables x and y.
% Input:
%   x, y: two integer vector of the same length
% Output:
%   z: mutual information z=I(x,y)
% Written by Mo Chen (sth4nth@gmail.com).
assert(numel(x) == numel(y));
n = numel(x);
x = reshape(x,1,n);
y = reshape(y,1,n);

l = min(min(x),min(y));
x = x-l+1;
y = y-l+1;
k = max(max(x),max(y));

idx = 1:n;
Mx = sparse(idx,x,1,n,k,n);
My = sparse(idx,y,1,n,k,n);
Pxy = nonzeros(Mx'*My/n); %joint distribution of x and y
Hxy = -dot(Pxy,log2(Pxy));

Px = nonzeros(mean(Mx,1));
Py = nonzeros(mean(My,1));

% entropy of Py and Px
Hx = -dot(Px,log2(Px));
Hy = -dot(Py,log2(Py));
% mutual information
z = Hx+Hy-Hxy;
z = max(0,z);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
classdef AltitudeTermCondition < AbstractEventTerminationCondition
    %AltitudeTermCondition Summary of this class goes here
    %   Detailed explanation goes here

    properties
        altitude(1,1) double = 0; %km
        bodyInfo KSPTOT_BodyInfo
    end

    methods
        function obj = AltitudeTermCondition(altitude)
            obj.altitude = altitude;
        end

        function evtTermCondFcnHndl = getEventTermCondFuncHandle(obj)
            evtTermCondFcnHndl = @(t,y) obj.eventTermCond(t,y);
        end

        function initTermCondition(obj, initialStateLogEntry)
            obj.bodyInfo = initialStateLogEntry.centralBody;
        end

        function name = getName(obj)
            name = sprintf('Altitude (%.3f km)', obj.altitude);
        end

        function tf = shouldBeReinitOnRestart(obj)
            tf = true;
        end

        function params = getTermCondUiStruct(obj)
            params = struct();

            params.paramName = 'Altitude';
            params.paramUnit = 'km';
            params.useParam = 'on';
            params.useStages = 'off';
            params.useTanks = 'off';
            params.useEngines = 'off';
            params.useStopwatches = 'off';

            params.value = obj.altitude;
            params.refStage = LaunchVehicleStage.empty(1,0);
            params.refTank = LaunchVehicleEngine.empty(1,0);
            params.refEngine = LaunchVehicleEngine.empty(1,0);
            params.refStopwatch = LaunchVehicleStopwatch.empty(1,0);
        end

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

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

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end
    end

    methods(Static)
        function termCond = getTermCondForParams(paramValue, stage, tank, engine, stopwatch)
            termCond = AltitudeTermCondition(paramValue);
        end
    end

    methods(Access=private)
        function [value,isterminal,direction] = eventTermCond(obj, t,y)
            rVect = y(1:3);
            vVect = y(4:6);
            cartElem = CartesianElementSet(t, rVect(:), vVect(:), obj.bodyInfo.getBodyCenteredInertialFrame());
            geoElem = cartElem.convertToFrame(obj.frame).convertToGeographicElementSet();

            actualAltitude = geoElem.alt;

            value = actualAltitude - obj.altitude;
            isterminal = 1;
            direction = 0;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function output=SSMultibandKamath02(signal,fs,IS)

% OUTPUT=SSMULTIBANDKAMATH02(S,FS,IS)
% Multi-band Spectral subtraction [Kamath2002]
% subtraction with adjusting subtraction factor. the adjustment is
% according to local a postriori SNR and the frequency band.
% S is the noisy signal, FS is the sampling frequency and IS is the initial
% silence (noise only) length in seconds (default value is .25 sec)
%
% April-04
% Esfandiar Zavarehei

if (nargin<3 | isstruct(IS))
    IS=.25; %seconds
end
W=fix(.025*fs); %Window length is 25 ms
nfft=W;
SP=.4; %Shift percentage is 40% (10ms) %Overlap-Add method works good with this value(.4)
wnd=hamming(W);

% IGNORE THIS SECTION FOR CAMPATIBALITY WITH ANOTHER PROGRAM FROM HERE.....
if (nargin>=3 & isstruct(IS))%This option is for compatibility with another programme
    W=IS.windowsize
    SP=IS.shiftsize/W;
    nfft=IS.nfft;
    wnd=IS.window;
    if isfield(IS,'IS')
        IS=IS.IS;
    else
        IS=.25;
    end
end
% .......IGNORE THIS SECTION FOR CAMPATIBALITY WITH ANOTHER PROGRAM T0 HERE

NIS=fix((IS*fs-W)/(SP*W) +1);%number of initial silence segments
Gamma=2;%Magnitude Power (1 for magnitude spectral subtraction 2 for power spectrum subtraction)

y=segment(signal,W,SP,wnd);
Y=fft(y,nfft);
YPhase=angle(Y(1:fix(end/2)+1,:)); %Noisy Speech Phase
Y=abs(Y(1:fix(end/2)+1,:)).^Gamma;%Specrogram
numberOfFrames=size(Y,2);
FreqResol=size(Y,1);

N=mean(Y(:,1:NIS)')'; %initial Noise Power Spectrum mean

NoiseCounter=0;
NoiseLength=9;%This is a smoothing factor for the noise updating

Beta=.03;
minalpha=1;
maxalpha=5;
minSNR=-5;
maxSNR=20;
alphaSlope=(minalpha-maxalpha)/(maxSNR-minSNR);
alphaShift=maxalpha-alphaSlope*minSNR;

BN=Beta*N;

%Delta is a frequency dependent coefficient
Delta=1.5*ones(size(BN));
Delta(1:fix((-2000+fs/2)*FreqResol*2/fs))=2.5; %if the frequency is lower than FS/2 - 2KHz
Delta(1:fix(1000*FreqResol*2/fs))=1; %if the frequency is lower than 1KHz

for i=1:numberOfFrames
    [NoiseFlag, SpeechFlag, NoiseCounter, Dist]=vad(Y(:,i).^(1/Gamma),N.^(1/Gamma),NoiseCounter); %Magnitude Spectrum Distance VAD
    if SpeechFlag==0
        N=(NoiseLength*N+Y(:,i))/(NoiseLength+1); %Update and smooth noise
        BN=Beta*N;
    end

    SNR=10*log(Y(:,i)./N);
    alpha=alphaSlope*SNR+alphaShift;
    alpha=max(min(alpha,maxalpha),minalpha);

    D=Y(:,i)-(Delta.*alpha.*N); %Nonlinear (Non-uniform) Power Specrum Subtraction

    X(:,i)=max(D,BN); %if BY>D X=BY else X=D which sets very small values of subtraction result to an attenuated
                      %version of the input power spectrum.
end

output=OverlapAdd2(X.^(1/Gamma),YPhase,W,SP*W);

function ReconstructedSignal=OverlapAdd2(XNEW,yphase,windowLen,ShiftLen);

%Y=OverlapAdd(X,A,W,S);
%Y is the signal reconstructed signal from its spectrogram. X is a matrix
%with each column being the fft of a segment of signal. A is the phase
%angle of the spectrum which should have the same dimension as X. if it is
%not given the phase angle of X is used which in the case of real values is
%zero (assuming that its the magnitude). W is the window length of time
%domain segments if not given the length is assumed to be twice as long as
%fft window length. S is the shift length of the segmentation process ( for
%example in the case of non overlapping signals it is equal to W and in the
%case of %50 overlap is equal to W/2. if not givven W/2 is used. Y is the
%reconstructed time domain signal.
%Sep-04
%Esfandiar Zavarehei

if nargin<2
    yphase=angle(XNEW);
end
if nargin<3
    windowLen=size(XNEW,1)*2;
end
if nargin<4
    ShiftLen=windowLen/2;
end
if fix(ShiftLen)~=ShiftLen
    ShiftLen=fix(ShiftLen);
    disp('The shift length have to be an integer as it is the number of samples.')
    disp(['shift length is fixed to ' num2str(ShiftLen)])
end

[FreqRes FrameNum]=size(XNEW);

Spec=XNEW.*exp(j*yphase);

if mod(windowLen,2) %if FreqResol is odd
    Spec=[Spec;flipud(conj(Spec(2:end,:)))];
else
    Spec=[Spec;flipud(conj(Spec(2:end-1,:)))];
end
sig=zeros((FrameNum-1)*ShiftLen+windowLen,1);
weight=sig;
for i=1:FrameNum
    start=(i-1)*ShiftLen+1;
    spec=Spec(:,i);
    sig(start:start+windowLen-1)=sig(start:start+windowLen-1)+real(ifft(spec,windowLen));
end
ReconstructedSignal=sig;

function [NoiseFlag, SpeechFlag, NoiseCounter, Dist]=vad(signal,noise,NoiseCounter,NoiseMargin,Hangover)

%[NOISEFLAG, SPEECHFLAG, NOISECOUNTER, DIST]=vad(SIGNAL,NOISE,NOISECOUNTER,NOISEMARGIN,HANGOVER)
%Spectral Distance Voice Activity Detector
%SIGNAL is the the current frames magnitude spectrum which is to labeld as
%noise or speech, NOISE is noise magnitude spectrum template (estimation),
%NOISECOUNTER is the number of imediate previous noise frames, NOISEMARGIN
%(default 3)is the spectral distance threshold. HANGOVER ( default 8 )is
%the number of noise segments after which the SPEECHFLAG is reset (goes to
%zero). NOISEFLAG is set to one if the the segment is labeld as noise
%NOISECOUNTER returns the number of previous noise segments, this value is
%reset (to zero) whenever a speech segment is detected. DIST is the
%spectral distance.
%Saeed Vaseghi
%edited by Esfandiar Zavarehei
%Sep-04

if nargin<4
    NoiseMargin=3;
end
if nargin<5
    Hangover=8;
end
if nargin<3
    NoiseCounter=0;
end

FreqResol=length(signal);

SpectralDist= 20*(log10(signal)-log10(noise));
SpectralDist(find(SpectralDist<0))=0;

Dist=mean(SpectralDist);
if (Dist < NoiseMargin)
    NoiseFlag=1;
    NoiseCounter=NoiseCounter+1;
else
    NoiseFlag=0;
    NoiseCounter=0;
end

% Detect noise only periods and attenuate the signal
if (NoiseCounter > Hangover)
    SpeechFlag=0;
else
    SpeechFlag=1;
end

function Seg=segment(signal,W,SP,Window)

% SEGMENT chops a signal to overlapping windowed segments
% A= SEGMENT(X,W,SP,WIN) returns a matrix which its columns are segmented
% and windowed frames of the input one dimentional signal, X. W is the
% number of samples per window, default value W=256. SP is the shift
% percentage, default value SP=0.4. WIN is the window that is multiplied by
% each segment and its length should be W. the default window is hamming
% window.
% 06-Sep-04
% Esfandiar Zavarehei

if nargin<3
    SP=.4;
end
if nargin<2
    W=256;
end
if nargin<4
    Window=hamming(W);
end
Window=Window(:); %make it a column vector

L=length(signal);
SP=fix(W.*SP);
N=fix((L-W)/SP +1); %number of segments

Index=(repmat(1:W,N,1)+repmat((0:(N-1))'*SP,1,W))';
hw=repmat(Window,1,N);
Seg=signal(Index).*hw;

author: Sable

repo: mcbench-benchmarks


text:
function I2 = undistortImage(I,K,kc)
% function to undo lens distortion with perspective camera intrinsics K and kc
% This is the model from Oulu university and described clearly in Bouguet's Matlab camera calibration toolbox
% (http://www.vision.caltech.edu/bouguetj/calib_doc/). This is done using
% Bouguet's toolbox (rect). You need to download it.

if ~exist('rect.m','file')
    error('To undistort an image you need to download Bouguets Matlab camera calibration toolbox (http://www.vision.caltech.edu/bouguetj/calib_doc/) and add it to the path.')
end
if isa(I,'uint8')
   I = im2double(I);
end
for i=1:size(I,3)
    I2(:,:,i) =rect((I(:,:,i)),eye(3),[K(1,1),K(2,2)],[K(1,3),K(2,3)],kc,K(1,2)/K(1,1),K);
end
I2 = uint8(I2*255);

author: tobycollins

repo: IPPE

save_path: github-repos/MATLAB/tobycollins-IPPE


text:
classdef MaF1 < PROBLEM
% <multi/many> <real> <large/none>
% Inverted DTLZ1

%------------------------------- Reference --------------------------------
% R. Cheng, M. Li, Y. Tian, X. Zhang, S. Yang, Y. Jin, and X. Yao, A
% benchmark test suite for evolutionary many-objective optimization,
% Complex & Intelligent Systems, 2017, 3(1): 67-81.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = obj.M + 9; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            g      = sum((PopDec(:,obj.M:end)-0.5).^2,2);
            PopObj = repmat(1+g,1,obj.M) - repmat(1+g,1,obj.M).*fliplr(cumprod([ones(size(g,1),1),PopDec(:,1:obj.M-1)],2)).*[ones(size(g,1),1),1-PopDec(:,obj.M-1:-1:1)];
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = 1 - UniformPoint(N,obj.M);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a = linspace(0,1,10)';
                R = {1-a*a',1-a*(1-a'),1-(1-a)*ones(size(a'))};
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function g = dexpKernGradient(kern, x1, varargin)

% DEXPKERNGRADIENT Gradient of the double exponential kernel's parameters.
%
% FORMAT
% DESC computes the gradient of functions with respect to the double
% exponential kernel's parameters. As well as the kernel structure
% and the input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the relevant
% elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations for which the gradients are being
% computed in the form of a design matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in x1.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix in the form of a design matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix in the form of a design matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as x1 and the same number of columns
% as x2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO dexpKernParamInit, kernGradient, dexpKernDiagGradient, kernGradX
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin < 4
    x2 = x1;
else
    x2 = varargin{1};
end

[K, sK, n1] = dexpKernCompute(kern, x1, x2);

g(1) = kern.variance * sum(sum(varargin{end} .* sK ...
    .* (1 - n1*kern.decay)));
g(2) = kern.decay * sum(sum(varargin{end} .* sK));

%/~
if any(isnan(g))
  warning('g is NaN')
end
%~/


text:
function [net, state] = DnCNN_train(sigma_min,sigma_max, net, varargin)

%    The function automatically restarts after each training epoch by
%    checkpointing.
%
%    The function supports training on CPU or on one or more GPUs
%    (specify the list of GPU IDs in the `gpus` option).

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

%%%-------------------------------------------------------------------------
%%% solvers: SGD(default) and Adam with(default)/without gradientClipping
%%%-------------------------------------------------------------------------

%%% solver: Adam
%%% opts.solver = 'Adam';
opts.beta1   = 0.9;
opts.beta2   = 0.999;
opts.alpha   = 0.01;
opts.epsilon = 1e-8;

%%% solver: SGD
opts.solver = 'SGD';
opts.learningRate = 0.01;
opts.weightDecay  = 0.001;
opts.momentum     = 0.9 ;

%%% GradientClipping
opts.gradientClipping = false;
opts.theta            = 0.005;

%%% specific parameter for Bnorm
opts.bnormLearningRate = 0;

%%%-------------------------------------------------------------------------
%%%  setting for simplenn
%%%-------------------------------------------------------------------------

opts.conserveMemory = true;
opts.mode = 'normal';
opts.cudnn = true ;
opts.backPropDepth = +inf ;
opts.skipForward = false;
opts.numSubBatches = 1;
%%%-------------------------------------------------------------------------
%%%  setting for model
%%%-------------------------------------------------------------------------

opts.batchSize = 128 ;
opts.gpus = [];
opts.numEpochs = 300 ;
opts.modelName   = 'model';
opts.expDir = fullfile('data',opts.modelName) ;
opts.numberData   = 1;
opts.DataDir      = '../../../LDAMP_TensorFlow/TrainingData/traindata_50_128.mat';
opts.ValDataDir          = '../../../LDAMP_TensorFlow/TrainingData/valdata_50_128.mat';

%%%-------------------------------------------------------------------------
%%%  update settings
%%%-------------------------------------------------------------------------

opts = vl_argparse(opts, varargin);
opts.numEpochs = numel(opts.learningRate);

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end

%%%-------------------------------------------------------------------------
%%%  Initialization
%%%-------------------------------------------------------------------------

net = vl_simplenn_tidy(net);    %%% fill in some eventually missing values
net.layers{end-1}.precious = 1;
vl_simplenn_display(net, 'batchSize', opts.batchSize) ;

state.getBatch = getBatch ;

%%%-------------------------------------------------------------------------
%%%  Train and Test
%%%-------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf([opts.modelName,'-epoch-%d.mat'], ep));

start = findLastCheckpoint(opts.expDir,opts.modelName) ;
if start >= 1
    fprintf('%s: resuming by loading epoch %d', mfilename, start) ;
    load(modelPath(start), 'net') ;
    net = vl_simplenn_tidy(net) ;
end

%%% load training and validation data

opts.DataPath    = fullfile(opts.DataDir);
TrainData = load(opts.DataPath) ;
opts.train = find(TrainData.set==1);

opts.ValDataPath    = fullfile(opts.ValDataDir);
ValData = load(opts.ValDataPath) ;
opts.test = find(ValData.set==0);

bad_epochs=0;
max_bad_epochs=2;
for epoch = start+1 : opts.numEpochs
    %%% Train for one epoch.

    %Stop training if three epochs produced no improvement in validation
    %error
    if bad_epochs>=max_bad_epochs
        if opts.learningRate(epoch-max_bad_epochs)==.001
            load(modelPath(epoch-max_bad_epochs-1), 'net')%load the last best epoch
            save(modelPath(epoch-1), 'net')%save the just-loaded model as the best, in case it needs to be loaded because errors don't improve
            save(fullfile(opts.expDir, sprintf([opts.modelName,'-best_epoch',num2str(epoch-1)])),'net');
            opts.learningRate(epoch:end)=.0001;
            bad_epochs=0;
        elseif opts.learningRate(epoch-max_bad_epochs)==.0001
            load(fullfile(opts.expDir, sprintf([opts.modelName,'-best_epoch',num2str(epoch-max_bad_epochs-1)])),'net');
            opts.learningRate(epoch:end)=.00001;
            bad_epochs=0;
        else
            break;
        end
    end

    state.epoch = epoch ;
    state.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate)));
    opts.thetaCurrent = opts.theta(min(epoch, numel(opts.theta)));
    if numel(opts.gpus) == 1
        net = vl_simplenn_move(net, 'gpu') ;
    end

    if epoch==start+1
        %Determine Initial Validation Error
        tic
        state.test = opts.test(randperm(numel(opts.test))) ; %%% shuffle
        [net, state, min_val_loss] = process_epoch(sigma_min,sigma_max,net, state, ValData, opts, 'test');
        net.layers{end}.class =[];
        toc
    end

    %Perform Training
    t0=cputime;
    state.train = opts.train(randperm(numel(opts.train))) ; %%% shuffle
    [net, state, ~] = process_epoch(sigma_min,sigma_max,net, state, TrainData, opts, 'train');
    net.layers{end}.class =[];
    time_taken=cputime-t0;

    %Determine Validation Error
    state.test = opts.test(randperm(numel(opts.test))) ; %%% shuffle
    [net, state, loss] = process_epoch(sigma_min,sigma_max,net, state, ValData, opts, 'test');
    net.layers{end}.class =[];

    net = vl_simplenn_move(net, 'cpu');

    %Save model
    save(modelPath(epoch), 'net')

    %Save as best model if validation error is low
    if loss<min_val_loss
        save(fullfile(opts.expDir, sprintf([opts.modelName,'-best_epoch',num2str(epoch)])),'net');
        min_val_loss=loss;
        bad_epochs=0;
    else
        bad_epochs=bad_epochs+1;
    end
end

%%%-------------------------------------------------------------------------
function  [net, state, total_l2_loss] = process_epoch(sigma_min,sigma_max,net, state, Data, opts, mode)
%%%-------------------------------------------------------------------------

if strcmp(mode,'train')

    switch opts.solver

        case 'SGD' %%% solver: SGD
            for i = 1:numel(net.layers)
                if isfield(net.layers{i}, 'weights')
                    for j = 1:numel(net.layers{i}.weights)
                        state.layers{i}.momentum{j} = 0;
                    end
                end
            end

        case 'Adam' %%% solver: Adam
            for i = 1:numel(net.layers)
                if isfield(net.layers{i}, 'weights')
                    for j = 1:numel(net.layers{i}.weights)
                        state.layers{i}.t{j} = 0;
                        state.layers{i}.m{j} = 0;
                        state.layers{i}.v{j} = 0;
                    end
                end
            end

    end

end

subset = state.(mode) ;
num = 0 ;
res = [];
total_l2_loss=0;
for t=1:opts.batchSize:numel(subset)

    for s=1:opts.numSubBatches
        % get this image batch
        batchStart = t + (s-1);
        batchEnd = min(t+opts.batchSize-1, numel(subset)) ;
        batch = subset(batchStart : opts.numSubBatches : batchEnd) ;
        num = num + numel(batch) ;
        if numel(batch) == 0, continue ; end

        [inputs,labels] = state.getBatch(sigma_min,sigma_max,Data, batch) ;
%         truth=inputs-labels;
%         imshow(truth(:,:,1,1),[]);

        if numel(opts.gpus) == 1
            inputs = gpuArray(inputs);
            labels = gpuArray(labels);
        end

        if strcmp(mode, 'train')
            dzdy = single(1);
            evalMode = 'normal';%%% forward and backward (Gradients)
        else
            dzdy = [] ;
            evalMode = 'test';  %%% forward only
%             evalMode = 'normal';  %%% forward and backward (Gradients)%Forces the network not to remove batch normalization.
        end

        net.layers{end}.class = labels ;
        res = vl_simplenn(net, inputs, dzdy, res, ...
            'accumulate', s ~= 1, ...
            'mode', evalMode, ...
            'conserveMemory', opts.conserveMemory, ...
            'backPropDepth', opts.backPropDepth, ...
            'cudnn', opts.cudnn) ;

    end

    if strcmp(mode, 'train')
        [state, net] = params_updates(state, net, res, opts, opts.batchSize) ;
    end

    lossL2 = gather(res(end).x) ;
    %%%--------add your code here------------------------
    if mod(floor(t/opts.batchSize),60)==1%show figure every 20 batches
        k=2;
        label_hat=res(end-1).x;
        label_hat_k=squeeze(label_hat(:,:,1,k));
        label_k=squeeze(labels(:,:,1,k));
        sig_hat_k=squeeze(inputs(:,:,1,k))-label_hat_k;
        sig_k=squeeze(inputs(:,:,1,k))-label_k;
        figure(1);
        subplot(3,2,1);imshow(inputs(:,:,1,k),[]);title('Noisy Input');
        subplot(3,2,3);imshow(label_hat_k,[]);title('Noise Estimate');
        subplot(3,2,4);imshow(label_k,[]);title('Noise Truth');
        subplot(3,2,5);imshow(sig_hat_k,[]);title('Signal Estimate');
        subplot(3,2,6);imshow(sig_k,[]);title('Signal Truth');
        Signals=inputs-labels;
        Signal_Estimates=inputs-res(end-1).x;
        MyLossL2=vl_nnloss(Signals,Signal_Estimates);%confirm this value is the same as lossL2
    end

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

    fprintf('%s: epoch %02d dataset %02d: %3d/%3d:', mode, state.epoch, mod(state.epoch,opts.numberData), ...
        fix((t-1)/opts.batchSize)+1, ceil(numel(subset)/opts.batchSize)) ;
    fprintf('error: %f \n', lossL2) ;
    total_l2_loss=total_l2_loss+lossL2;
end

%%%-------------------------------------------------------------------------
function [state, net] = params_updates(state, net, res, opts, batchSize)
%%%-------------------------------------------------------------------------

switch opts.solver

    case 'SGD' %%% solver: SGD

        for l=numel(net.layers):-1:1
            for j=1:numel(res(l).dzdw)
                if j == 3 && strcmp(net.layers{l}.type, 'bnorm')
                    %%% special case for learning bnorm moments
                    thisLR = net.layers{l}.learningRate(j) - opts.bnormLearningRate;
                    net.layers{l}.weights{j} = vl_taccum(...
                        1 - thisLR, ...
                        net.layers{l}.weights{j}, ...
                        thisLR / batchSize, ...
                        res(l).dzdw{j}) ;

                else
                    thisDecay = opts.weightDecay * net.layers{l}.weightDecay(j);
                    thisLR = state.learningRate * net.layers{l}.learningRate(j);

                    if opts.gradientClipping
                        theta = opts.thetaCurrent/thisLR;
                        state.layers{l}.momentum{j} = opts.momentum * state.layers{l}.momentum{j} ...
                            - thisDecay * net.layers{l}.weights{j} ...
                            - (1 / batchSize) * gradientClipping(res(l).dzdw{j},theta) ;
                        net.layers{l}.weights{j} = net.layers{l}.weights{j} + ...
                            thisLR * state.layers{l}.momentum{j} ;
                    else
                        state.layers{l}.momentum{j} = opts.momentum * state.layers{l}.momentum{j} ...
                            - thisDecay * net.layers{l}.weights{j} ...
                            - (1 / batchSize) * res(l).dzdw{j} ;
                        net.layers{l}.weights{j} = net.layers{l}.weights{j} + ...
                            thisLR * state.layers{l}.momentum{j} ;
                    end
                end
            end
        end

    case 'Adam'  %%% solver: Adam

        for l=numel(net.layers):-1:1
            for j=1:numel(res(l).dzdw)

                if j == 3 && strcmp(net.layers{l}.type, 'bnorm')

                    %%% special case for learning bnorm moments
                    thisLR = net.layers{l}.learningRate(j) - opts.bnormLearningRate;
                    net.layers{l}.weights{j} = vl_taccum(...
                        1 - thisLR, ...
                        net.layers{l}.weights{j}, ...
                        thisLR / batchSize, ...
                        res(l).dzdw{j}) ;
                else
                    thisLR = state.learningRate * net.layers{l}.learningRate(j);
                    state.layers{l}.t{j} = state.layers{l}.t{j} + 1;
                    t = state.layers{l}.t{j};
                    alpha = thisLR;
                    lr = alpha * sqrt(1 - opts.beta2^t) / (1 - opts.beta1^t);

                    state.layers{l}.m{j} = state.layers{l}.m{j} + (1 - opts.beta1) .* (res(l).dzdw{j} - state.layers{l}.m{j});
                    state.layers{l}.v{j} = state.layers{l}.v{j} + (1 - opts.beta2) .* (res(l).dzdw{j} .* res(l).dzdw{j} - state.layers{l}.v{j});

                    if opts.gradientClipping
                        theta = opts.thetaCurrent/lr;
                        net.layers{l}.weights{j} = net.layers{l}.weights{j} - lr * gradientClipping(state.layers{l}.m{j} ./ (sqrt(state.layers{l}.v{j}) + opts.epsilon),theta);
                    else
                        net.layers{l}.weights{j} = net.layers{l}.weights{j} - lr * state.layers{l}.m{j} ./ (sqrt(state.layers{l}.v{j}) + opts.epsilon);
                    end

                    % net.layers{l}.weights{j} = weightClipping(net.layers{l}.weights{j},2); % gradually clip the weights

                end
            end
        end
end

%%%-------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir,modelName)
%%%-------------------------------------------------------------------------
list = dir(fullfile(modelDir, [modelName,'-epoch-*.mat'])) ;
tokens = regexp({list.name}, [modelName,'-epoch-([\d]+).mat'], 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

%%%-------------------------------------------------------------------------
function A = gradientClipping(A, theta)
%%%-------------------------------------------------------------------------
A(A>theta)  = theta;
A(A<-theta) = -theta;

%%%-------------------------------------------------------------------------
function A = weightClipping(A, theta)
%%%-------------------------------------------------------------------------
A(A>theta)  = A(A>theta) -0.0005;
A(A<-theta) = A(A<-theta)+0.0005;

%%%-------------------------------------------------------------------------
function fn = getBatch
%%%-------------------------------------------------------------------------
fn = @(sigma_min,sigma_max,x,y) getSimpleNNBatch(sigma_min,sigma_max,x,y);

%%%-------------------------------------------------------------------------
function [inputs,labels] = getSimpleNNBatch(sigma_min,sigma_max,Data, batch)
%%%-------------------------------------------------------------------------
inputs = Data.inputs(:,:,:,batch);
rng('shuffle');
mode = randperm(8);
inputs = data_augmentation(inputs, mode(1));
sigma_vec=sigma_min+rand(numel(batch),1)*(sigma_max-sigma_min);%amount of random noise to apply to each image
sigma_array=ones(size(inputs));
for i=1:numel(batch)
    sigma_array(:,:,1,i)=sigma_vec(i);
end
labels  = sigma_array/255.*randn(size(inputs),'single');
inputs = inputs + labels;

function image = data_augmentation(image, mode)

if mode == 1
    return;
end

if mode == 2 % flipped
    image = flipud(image);
    return;
end

if mode == 3 % rotation 90
    image = rot90(image,1);
    return;
end

if mode == 4 % rotation 90 & flipped
    image = rot90(image,1);
    image = flipud(image);
    return;
end

if mode == 5 % rotation 180
    image = rot90(image,2);
    return;
end

if mode == 6 % rotation 180 & flipped
    image = rot90(image,2);
    image = flipud(image);
    return;
end

if mode == 7 % rotation 270
    image = rot90(image,3);
    return;
end

if mode == 8 % rotation 270 & flipped
    image = rot90(image,3);
    image = flipud(image);
    return;
end

author: ricedsp


text:
classdef ShapesPrinter < CompositeResultsPrinter

    properties (Access = private)
        allShapes
        printableIndex
    end

    methods (Access = public)

        function obj = ShapesPrinter(d)
            obj.simulationStr = 'ShapePrinters';
            obj.obtainAllShapes(d);
            obj.obtainPrintableIndex();
            obj.init(d);
        end

    end

    methods (Access = protected)

        function createPrinters(obj,d)
            pI = obj.printableIndex;
            aS = obj.allShapes;
            nShapes = numel(pI);
            obj.printers = cell(nShapes,1);
            for ishape = 1:nShapes
                index = obj.printableIndex(ishape);
                shape = aS{index};
                d.fieldToPrint = shape.createPrintVariables();
                p = ShapePrinter(d);
                obj.printers{ishape} = p;
            end
        end

        function storeFieldsToPrint(obj,d)
            obj.obtainAllShapes(d);
            obj.obtainPrintableIndex();
            for iprinter = 1:numel(obj.printers)
                index = obj.printableIndex(iprinter);
                shape = obj.allShapes{index};
                fP = shape.addPrintableVariables();
                d.fieldToPrint = fP;
                p = obj.printers{iprinter};
                p.storeFieldsToPrint(d);
            end
        end

        function createHeadPrinter(obj,d,dh)
            for iprinter = 1:numel(obj.printers)
                index = obj.printableIndex(iprinter);
                shape = obj.allShapes{index};
                if obj.getHasGaussData()
                  if ~isfield(d,'quad')
                      d.quad = shape.getQuad();
                  else
                      if isempty(d.quad)
                        d.quad = shape.getQuad();
                      end
                  end
                else
                    d.quad = [];
                end
                p = obj.printers{iprinter};
                p.createHeadPrinter(d,dh);
                if p.getHasGaussData()
                    h = p.getHeadPrinter();
                    obj.headPrinter = h;
                    return
                else
                    h = p.getHeadPrinter();
                    obj.headPrinter = h;
                end
            end
        end

    end

    methods (Access = private)

        function obtainAllShapes(obj,d)
            obj.obtainCostShapes(d);
            obj.obtainConstraintShapes(d);
        end

        function obtainCostShapes(obj,d)
            nShape = numel(d.cost.shapeFunctions);
            for ishape = 1:nShape
                obj.allShapes{ishape} = d.cost.shapeFunctions{ishape};
            end
        end

        function obtainConstraintShapes(obj,d)
            nShape = numel(d.constraint.shapeFunctions);
            for ishape = 1:nShape
                obj.allShapes{end + ishape} = d.constraint.shapeFunctions{ishape};
            end
        end

        function obtainPrintableIndex(obj)
            aShapes = obj.allShapes;
            nAllShapes = numel(aShapes);
            iprint = 0;
            for ishape = 1:nAllShapes
                shape = aShapes{ishape};
                if obj.isShapePrintable(shape)
                    iprint = iprint + 1;
                    obj.printableIndex(iprint) = ishape;
                end
            end
        end

        function itIs = isShapePrintable(obj,shapeFunc)
            shapeName = class(shapeFunc);
            printingShapes = {'ShFunc_NonSelfAdjoint_Compliance',...
                'ShFunc_Compliance', 'ShFunc_Compliance',...
                'ShFunc_Chomog_alphabeta','ShFunc_Chomog_EnforceCh_CCstar_L2',...
                'ShFunc_Chomog_fraction','ShFunc_Perimeter',...
                'ShFunc_StressNorm','ShFunc_StressNorm2',...
                'ShFunc_StressNorm3'};
            itIs = any(strcmp(printingShapes,shapeName));
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% rundesignsim
clear;
figure;
index = 1;
name = input('Enter the name of the output variable and files ');
eval([name ' = zeros(10,10);'])

% set initial frequency of rests
	H = findobj('Tag', 'E2');
	freqConditions = str2num(get(H(1), 'String'));
	freqConditions = num2str([freqConditions(1:4) .01]);
	set(H,'String',freqConditions);

	% set initial ISI
	H = findobj('Tag', 'E4'); % ISI
	ISI = str2num(get(H(1), 'String'));
	ISI = num2str(1);
	set(H,'String',ISI);

myISI = .1:.2:4;
myFC = .01:.1:.81;

for RR = 1:size(myISI,2)
	% set ISI
	H = findobj('Tag', 'E4'); % ISI
	ISI = str2num(get(H(1), 'String'));
	ISI = num2str(myISI(RR));
	set(H(1),'String',ISI);

	for CC = 1:size(myFC,2)
		disp(['Starting model ' num2str(index)])

		% set rest frequency
		H = findobj('Tag', 'E2');
		freqConditions = str2num(get(H(1), 'String'));
		freqConditions(end) = myFC(CC);
		freqConditions(1:end-1) = (1 - myFC(CC)) / (size(freqConditions,2)-1)
		freqConditions = num2str(freqConditions);
		set(H(1),'String',freqConditions);

		disp(['FREQC and ISI are:	' get(H(1),'String') ' and ' get(findobj('Tag','E4'),'String')])
		designsim_gui_script

		H = findobj('Tag','SIMplot');
		if ~isempty(H)
			figure(H)
			close
		end

		disp(['mean is ' num2str(mean(SIM.rnd.t))])
		eval([name '(RR,CC) = mean(SIM.rnd.t);']);
		eval([name 'HP(RR,CC) = mean(SIM.rnd.HPt);']);
		eval([name 'HL(RR,CC) = mean(SIM.rnd.HLt);']);
		eval([name '(1:10,1:10)'])
		eval(['save ' name num2str(index) ' ' name ' SIM'])
		index = index + 1;
	end
end

eval(['save ' name])
eval(['surf(' name ')']); colormap(copper)
eval(['set(gca,''XTick'',1:size(' name ',2))'])
set(gca,'XTickLabel',myFC)
eval(['set(gca,''YTick'',1:size(' name ',1))'])
set(gca,'YTickLabel',myISI)
xlabel('Proportion of blank intervals','FontSize',18)
ylabel('ITI in s','FontSize',18)
title('Surface map of ITI vs. proportion of blank intervals','FontSize',18)

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function jac = p01_jac ( option, nvar, x )

%*****************************************************************************80
%
%% P01_JAC evaluates the jacobian for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the argument of the jacobian.
%
%    Output, real JAC(NVAR-1,NVAR), the jacobian matrix evaluated
%    at X.  The NVAR-th row is not set by this routine.
%
  jac = zeros ( nvar, nvar );

  jac(1,1) = 1.0;
  jac(2,1) = 1.0;

  jac(1,2) = ( - 3.0 * x(2) + 10.0 ) * x(2) -  2.0;
  jac(2,2) = (   3.0 * x(2) +  2.0 ) * x(2) - 14.0;
%
%  Get the starting point
%
  y = p01_start ( option, nvar );
%
%  Get the function value at the starting point
%
  gy = p01_gx ( y );

  jac(1,3) = gy(1);
  jac(2,3) = gy(2);

  return
end


text:
function [ vX ] = ProxHuberLossBoyd( vY, paramDelta, paramLambda )
% ----------------------------------------------------------------------------------------------- %
% [ vX ] = ProxHuberLossBoyd( vY, paramDelta, paramLambda )
%   Solves the Proximal Operator of the Huber Loss Function:
%   $$ \arg \min_{x} \frac{1}{2} {\left| x - y \right\|}_{2}^{2} + \lambda {H}_{\delta} \left( x \right) $$
%   Where {H}_{\delta} \left( x \right) is the Huber Loss Function.
% Input:
%   - vY            -   Input Vector.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - paramDelta    -   Parameter Delta.
%                       The Delta Parameter of the Huber Loss Function.
%                       This is the value the Huber Loss changes from
%                       L2 Norm to the L1 Norm.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: (0, inf).
%   - paramLambda   -   Parameter Lambda.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: (0, inf).
% Output:
%   - vX            -   Output Vector.
%                       The solution of the Proximal Operator.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  Huber Loss (Wikipedia) - https://en.wikipedia.org/wiki/Huber_loss.
%   2.  Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers (See Huber Fitting).
%   3.  Proximal Operator of the Huber Loss Function - https://math.stackexchange.com/questions/3589025.
% Remarks:
%   1.  In the reference by Boyd they use the $ \rho = 1 / \lambda $ form of
%       the Proximal Operator. Hence the adoption of the scaling. Also the
%       book use Huber Loss Function with $ \delta = 1 $.
% TODO:
%   1.  U.
% Release Notes:
%   -   1.0.000     21/03/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

hProxL1 = @(vX, paramLambda) sign(vX) .* max(abs(vX) - paramLambda, 0);
vX = ((1 / (1 + paramLambda)) * vY) + (paramDelta * (paramLambda / (1 + paramLambda)) * hProxL1((vY / paramDelta), 1 + paramLambda));

end

author: RoyiAvital

repo: StackExchangeCodes


text:
%% OPT13_RUN
%
%  Modified:
%
%    12 February 2008
%
   %---------------------------------------------------------------------
   %  Use NEWTON method.
   %---------------------------------------------------------------------
   fprintf('---------------------------------------------------------\n')
   fprintf('Test case 13 as a standard minimization problem.\n' );
   fprintf('Use Newton''s method.\n' );
   fprintf('The function is badly scaled.\n' );
   fprintf('There is a local minimizer at X=(0.285,0.279), F(X)=5.92\n' );
   fprintf('There is a global minimizere at X=-21.02,36.76), F(X)=0.0\n' );
   fprintf('---------------------------------------------------------\n')

   fname = 'opt13_fgh';
   options = [];
   options.max_iterations     = 30;
   options.method             = 'newton';

   fprintf('Newton:\n')

   x0 = [ 1; 1 ];
   x = entrust(fname, x0, options);

   fprintf('Newton''s method produced  (%10.7e,%10.7e)\n\n',x(1),x(2))
   f = opt13_fgh ( x, 'f' );
   fprintf('Value of F(X) = %f\n', f );

   %---------------------------------------------------------------------
   %  Use GAUSS-NEWTON method on least squares problem.
   %---------------------------------------------------------------------
   fprintf('---------------------------------------------------------\n')
   fprintf('Running testcase_13 as least squares problem: \n')
   fprintf('---------------------------------------------------------\n')
   fname = 'opt13_rj';
   options = [];
   options.verbose            = 0;
   options.method             = 'gauss_newton';
   options.step_tolerance     = 1.e-15;
   options.globalization      = 'none';
   options.gradient_tolerance = 1.e-10;
   options.max_iterations     = 10000;
   options.max_fevals         = 10000;

   x0 = [ 1; 1 ];
   x = entrust(fname, x0, options);
   fprintf('Gauss-Newton produced  (%10.7e, %10.7e)\n\n',x(1),x(2))
   [ res, jac ] = opt13_rj ( x, 'f' );
   fprintf('Norm of RES(X) = %f\n', norm ( res ) );


text:
function c = nancov(x,varargin)
%NANCOV Covariance matrix, ignoring NaNs.
%   C = NANCOV(X), if X is a vector, returns the sample variance of the
%   values in X, treating NaNs as missing values.  For matrices, where
%   each row is an observation and each column a variable, NANCOV(X) is
%   the covariance matrix computing using rows of X that do not contain
%   any NaN values.  NANCOV(X,Y), where X and Y are matrices with
%   the same number of elements, is equivalent to NANCOV([X(:) Y(:)]).
%
%   NANCOV(X) or NANCOV(X,Y) normalizes by (N-1) if N>1, where N is the
%   number of observations after removing missing values.  This makes
%   NANCOV(X) the best unbiased estimate of the covariance matrix if the
%   observations are from a normal distribution. For N=1, COV normalizes
%   by N.
%
%   NANCOV(X,1) or NANCOV(X,Y,1) normalizes by N and produces the second
%   moment matrix of the observations about their mean.  NANCOV(X,Y,0) is
%   the same as NANCOV(X,Y), and NANCOV(X,0) is the same as NANCOV(X).
%
%   C = NANCOV(...,'pairwise') computes C(I,J) using rows with no NaN
%   values in columns I or J.  The result may not be a positive definite
%   matrix. C = NANCOV(...,'complete') is the default, and it omits rows
%   with any NaN values, even if they are not in column I or J.
%
%   The mean is removed from each column before calculating the
%   result.
%
%   Example:  Generate random data having non-zero covariance between
%             column 4 and the other columns.
%       x = randn(30,4);       % uncorrelated data
%       x(:,4) = sum(x,2);     % introduce correlation
%       x(2,3) = NaN;          % introduce one missing value
%       c = nancov(x)          % compute sample covariance
%
%   Class support for inputs X,Y:
%      float: double, single
%
%   See also COV, VAR, NANVAR.

%   Copyright 1984-2011 The MathWorks, Inc.
%   $Revision: 1.1.8.3 $  $Date: 2011/02/09 19:35:29 $

if nargin<1
   error(message('stats:nancov:NotEnoughInputs'));
end

% Should we ignore NaNs by complete rows or pairwise?
dopairwise = false;
if numel(varargin)>0
   temp = varargin{end};
   if ischar(temp)
      j = find(strncmpi(temp, {'pairwise' 'complete'},length(temp)));
      if isempty(j)
         error(message('stats:nancov:InvalidArg', temp));
      end
      dopairwise = (j==1);
      varargin(end) = [];
   end
end

% Should we use the mle (divide by N) or unbiased estimate (N-1)?
domle = false;
if numel(varargin)>0
   temp = varargin{end};
   if isequal(temp,0) || isequal(temp,1)
      domle = (temp==1);
      varargin(end) = [];
   end
end

if numel(varargin)>1
   error(message('stats:nancov:TooManyArgs'));
end

scalarxy = false; % nancov(scalar,scalar) is an ambiguous case
if numel(varargin)>0
   y = varargin{1};

   % Two inputs, convert to equivalent single input
   x = x(:);
   y = y(:);
   if length(x)~=length(y)
      error(message('stats:nancov:XYmismatch'));
   end
   scalarxy = isscalar(x) && isscalar(y);
   x = [x y];
elseif ndims(x)>2
   error(message('stats:nancov:InputDim'));
end

if isvector(x) && ~scalarxy
  x = x(:);
end

xnan = isnan(x);
[m,n] = size(x);

if isempty(x);
  if (m==0 && n==0)
      c = NaN(class(x));
  else
      c = NaN(n,class(x));
  end
  return;
end

if ~dopairwise || ~any(any(xnan))    % no need to do pairwise
   nanrows = any(xnan,2);
   if any(nanrows)
       x = x(~nanrows,:);
   end
   c = localcov(x,domle);
else                                         % pairwise with some NaNs
   % Compute variance using complete data separately by column
   c = zeros(n,class(x));
   x(xnan) = 0;
   colsize = sum(~xnan,1);
   xmean = sum(x,1) ./ max(1,colsize);
   xmean(colsize==0) = NaN;
   xc = x - repmat(xmean,m,1);
   xc(xnan) = 0;
   xvar = sum(xc.^2,1);
   if domle
      denom = colsize;
   else
      denom = max(0,colsize-1);
   end
   xvar(denom>1) = xvar(denom>1) ./ denom(denom>1);
   xvar(denom==0) = NaN;
   c(1:n+1:end) = xvar;

   % Now compute off-diagonal entries
   jk = 1:2;
   for j = 2:n
      jk(1) = j;
      for k=1:j-1
         jk(2) = k;
         rowsjk = ~any(xnan(:,jk),2);
         njk = sum(rowsjk);
         if njk<=1
            cjk = NaN;
         else
            cjk = localcov(x(rowsjk,jk),domle);
            cjk = cjk(1,2);
         end
         c(j,k) = cjk;
      end
   end
   c = c + tril(c,-1)';
end

% ------------------------------------------------
function [c,n] = localcov(x,domle)
%LOCALCOV Compute cov with no error checking and assuming NaNs are removed

[m,n] = size(x);
if domle
   denom = m;
else
   denom = max(0,m-1);
end

if m==1   % and doing mle, be sure to get exact 0
   c = zeros(n,class(x));
elseif denom==0
   c = NaN(n,class(x));
else
   xc = x - repmat(mean(x),m,1);
   c = xc' * xc / denom;
end

author: marianux

repo: ecg-kit

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


text:
function a = wishart_sample ( m, df, sigma )

%*****************************************************************************80
%
%% WISHART_SAMPLE samples the Wishart distribution.
%
%  Discussion:
%
%    This function requires functions from the PDFLIB and RNGLIB libraries.
%
%    The "initialize()" function from RNGLIB must be called before using
%    this function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 July 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Patrick Odell, Alan Feiveson,
%    A numerical procedure to generate a sample covariance matrix,
%    Journal of the American Statistical Association,
%    Volume 61, Number 313, March 1966, pages 199-203.
%
%    Stanley Sawyer,
%    Wishart Distributions and Inverse-Wishart Sampling,
%    Washington University,
%    30 April 2007, 12 pages.
%
%  Parameters:
%
%    Input, integer M, the order of the matrix.
%
%    Input, integer DF, the number of degrees of freedom.
%    M <= DF.
%
%    Input, real SIGMA(M,M), the covariance matrix, which should be
%    a symmetric positive definite matrix.
%
%    Output, real A(M,M), the sample matrix from the Wishart distribution.
%
  if ( df < m )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'WISHART_SAMPLE - Fatal error!\n' );
    fprintf ( 1, '  DF = %d < M = %d.\n', df, m );
    error ( 'WISHART_SAMPLE - Fatal error!\n' );
  end
%
%  Get R, the upper triangular Cholesky factor of SIGMA.
%
  r = chol ( sigma );
%
%  Get AU, a sample from the unit Wishart distribution.
%
  au = wishart_unit_sample ( m, df );
%
%  Construct the matrix A = R' * AU * R.
%
  a = r' * au * r;

  return
end


text:
function K = sdlfmvXsdlfmvKernCompute(sdlfmvKern1, sdlfmvKern2, t1, t2, covIC)

% SDLFMVXSDLFMVKERNCOMPUTE Compute a cross kernel between two SDLFMV kernels.
% FORMAT
% DESC computes cross kernel terms between the velocity of two switching
% dynamical LFM kernels for the multiple output kernel.
% ARG sdlfmvKern1 : the kernel structure associated with the velocity of the
% first SDLFM kernel.
% ARG sdlfmvKern2 : the kernel structure associated with the velocity of the
% second SDLFM kernel.
% ARG t : inputs for which kernel is to be computed.
% ARG covIC : covariance for the initial conditions
% RETURN K : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between the velocity of two SDLFM kernels
% for the multiple output kernel.
% ARG sdlfmvKern1 : the kernel structure associated with the velocity of
% the first SDLFM kernel.
% ARG sdlfmvKern2 : the kernel structure associated with the velocity of
% the second SDLFM kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covIC : covariance for the initial conditions
% RETURN K : block of values from kernel matrix.
%
% SEEALSO : sdlfmvKernParamInit, sdlfmvKernCompute, sdlfmKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 4
    t2 = t1;
end

K = sdlfmXsdlfmKernCompute(sdlfmvKern1, sdlfmvKern2, t1, t2, covIC, 'VelVel');


text:
function [z,output] = minf_lbfgsdl(f,g,z0,options)
%MINF_LBFGSDL Minimize a function by L-BFGS with dogleg trust region.
%   [z,output] = minf_lbfgsdl(f,g,z0) starts at z0 and attempts to find a
%   local minimizer of the real-valued function f(z). The input variables z
%   may be a scalar, vector, matrix, tensor or even a (nested) cell array
%   of tensors and its contents may be real or complex.
%
%   If f(x) is a function of real variables x, the function g(x) should
%   compute the partial derivatives of f with respect to the real variables
%   x, i.e. g(xk) := df(xk)/dx. If f(z) is a function of complex variables
%   z, the function g(z) should compute two times the partial derivative
%   of f with respect to conj(z) (treating z as constant), i.e. g(zk) :=
%   2*df(zk)/d(conj(z)) = 2*conj(df(zk)/dz). If g is the empty matrix [],
%   the real gradient or scaled conjugate cogradient is approximated with
%   finite differences. The output of the function g(z) may have the same
%   structure as z (although this is not necessary). The structure output
%   returns additional information:
%
%      output.alpha      - The plane search step lengths in every
%                          iteration, if a plane search is selected.
%      output.delta      - The trust region radius at every step attempt.
%      output.fevals     - The total number of function calls.
%      output.fval       - The value of the objective function f in every
%                          iteration.
%      output.gevals     - The total number of gradient calls.
%      output.info       - The circumstances under which the procedure
%                          terminated:
%                             1: Objective function tolerance reached.
%                             2: Step size tolerance reached.
%                             3: Maximum number of iterations reached.
%      output.infops     - The circumstances under which the plane search
%                          terminated in every iteration.
%      output.iterations - The number of iterations.
%      output.relfval    - The difference in objective function value
%                          between every two successive iterates, relative
%                          to its initial value.
%      output.relstep    - The step size relative to the norm of the
%                          current iterate in every iteration.
%      output.rho        - The trustworthiness at every step attempt.
%
%   minf_lbfgsdl(f,g,z0,options) may be used to set the following options:
%
%      options.Delta =            - The initial trust region radius.
%      0.3*max(1,norm(z0))
%      options.Display = 10       - Displays the objective function value,
%                                   its difference with the previous
%                                   iterate relative to the first iterate
%                                   and the relative step size each
%                                   options.Display iterations. Set to 0 to
%                                   disable.
%      options.M =                - The number of updates to store.
%      min(30,length(z0))
%      options.MaxIter = 500      - The maximum number of iterations.
%      options.PlaneSearch        - The plane search used to minimize the
%      = false                      objective function in the plane spanned
%                                   by the steepest descent direction and
%                                   the Gauss-Newton step. Disables dogleg
%                                   trust region strategy. The method
%                                   should have the function signature
%                                   options.PlaneSearch(F,dF,z,p1,p2, ...
%                                   state,options.PlaneSearchOptions).
%      options.PlaneSearchOptions - The options structure passed to the
%                                   plane search search routine.
%      options.TolFun = 1e-6      - The tolerance for output.relfval.
%      options.TolX = 1e-8        - The tolerance for output.relstep.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Unconstrained
%       optimization of real functions in complex variables", SIAM J. Opt.,
%       Vol. 22, No. 3, 2012, pp. 879-898.

% Store the structure of the input space and evaluate the gradient.
dim = structure(z0);
fval = f(z0);
if ~isa(g,'function_handle') && isempty(g)
    grad = serialize(deriv(f,z0,fval));
else
    grad = serialize(g(z0));
end
z0 = serialize(z0);

% Check the options structure.
if nargin < 4, options = struct; end
if ~isfield(options,'Delta'), options.Delta = 0.3*max(1,norm(z0)); end
if ~isfield(options,'Display'), options.Display = 10; end
if ~isfield(options,'M'), options.M = min(30,length(z0)); end
if ~isfield(options,'MaxIter'), options.MaxIter = 500; end
if ~isfield(options,'PlaneSearch'), options.PlaneSearch = false; end
if ~isfield(options,'PlaneSearchOptions')
    options.PlaneSearchOptions = struct;
end
if ~isfield(options,'TolFun'), options.TolFun = 1e-6; end
if ~isfield(options,'TolX'), options.TolX = 1e-8; end

% Initialize the algorithm.
S = zeros(numel(z0),options.M);
Y = zeros(numel(z0),options.M);
a = zeros(1,options.M);
r = zeros(1,options.M);
m = 0;
midx = 0;

% L-BFGS with dogleg trust region.
output.alpha = [];
output.delta = options.Delta;
output.fevals = 1;
output.fval = fval;
output.gevals = 1;
output.info = false;
output.infops = [];
output.iterations = 0;
output.relfval = [];
output.relstep = [];
output.rho = [];
while ~output.info

    % Compute the quasi-Newton step pqn = -H*grad.
    pqn = -grad;
    for i = 1:m
        a(i) = r(midx(i))*real(S(:,midx(i))'*pqn);
        pqn = pqn-a(i)*Y(:,midx(i));
    end
    if m > 0
        y1 = Y(:,midx(1));
        y1y1 = y1'*y1;
        gamma = y1y1*r(midx(1));
        pqn = 1/gamma*pqn;
    end
    for i = m:-1:1
        b = r(midx(i))*real(Y(:,midx(i))'*pqn);
        pqn = pqn+(a(i)-b)*S(:,midx(i));
    end

    % Approximate the Cauchy point pcp = -alpha*grad, where alpha is equal
    % to arg min m(-alpha*grad) and m(p) is a second order model of f at z.
    gg = grad'*grad;
    if m == 0
        alpha = 1;
    else
        s1 = S(:,midx(1));
        gBg = gg-real(grad'*s1)^2/(s1'*s1)+real(grad'*y1)^2/y1y1;
        gBg = gamma*gBg;
        alpha = gg/gBg;
    end

    % Plane search in the plane spanned by {pqn,pcp}.
    if isa(options.PlaneSearch,'function_handle')
        if output.iterations == 0, pcp = zeros(size(grad));
        else pcp = -alpha*grad; end
        state = output; state.grad = grad;
        [alpha,outputps] = options.PlaneSearch( ...
            f,g,deserialize(z0,dim),deserialize(pqn,dim), ...
            deserialize(pcp,dim),state,options.PlaneSearchOptions);
        output.alpha(:,end+1) = alpha;
        if length(alpha) < 3, alpha(3) = 1; end
        p = alpha(1)*pqn+alpha(2)*pcp;
        z = deserialize(alpha(3)*(z0+p),dim);
        relstep = norm(p)/norm(z0); if isnan(relstep), relstep = 0; end
        if isfield(outputps,'fval'), fval = outputps.fval;
        else fval = f(z); end
        if isfield(outputps,'info')
            output.infops(end+1) = outputps.info;
        end
        rho = 1;
    else
        rho = -inf;
    end

    % Dogleg trust region.
    normpqn = norm(pqn);
    while rho <= 0

        % Compute the dogleg step p.
        delta = output.delta(end);
        if normpqn <= delta
            p = pqn;
            dfval = -0.5*real(grad'*pqn);
        elseif abs(alpha)*sqrt(gg) >= delta
            p = (-delta/sqrt(gg))*grad;
            dfval = delta*(sqrt(gg)-0.5*delta/alpha);
        else
            bma = pqn+alpha*grad; bmabma = bma'*bma;
            a = -alpha*grad; aa = alpha^2*gg;
            c = real(a'*bma);
            if c <= 0
                beta = (-c+sqrt(c^2+bmabma*(delta^2-aa)))/bmabma;
            else
                beta = (delta^2-aa)/(c+sqrt(c^2+bmabma*(delta^2-aa)));
            end
            p = a+beta*bma;
            dfval = 0.5*alpha*(1-beta)^2*gg- ...
                    0.5*beta *(2-beta)*real(grad'*pqn);
        end

        % Compute the trustworthiness rho.
        if dfval > 0
            z = deserialize(z0+p,dim);
            fval = f(z);
            rho = (output.fval(end)-fval)/dfval;
            if isnan(rho), rho = -inf; end
            output.rho(end+1) = rho;
            output.fevals = output.fevals+1;
        end

        % Update trust region radius delta.
        if rho > 0.5
            output.delta(end+1) = max(delta,2*norm(p));
        else
            sigma = (1-0.25)/(1+exp(-14*(rho-0.25)))+0.25;
            if normpqn < sigma*delta && rho < 0
                e = ceil(log2(normpqn/delta)/log2(sigma));
                output.delta(end+1) = sigma^e*delta;
            else
                output.delta(end+1) = sigma*delta;
            end
        end

        % Check for convergence.
        relstep = norm(p)/norm(z0); if isnan(relstep), relstep = 0; end
        if rho <= 0 && relstep <= options.TolX
            output.rho(end+1) = rho;
            fval = output.fval(end);
            z = deserialize(z0,dim);
            break;
        end

    end

    % Save current state.
    if rho > 0
        z0 = serialize(z);
        grad1 = grad;
    end

    % Evaluate the gradient and update step information.
    if rho > 0
        if isa(options.PlaneSearch,'function_handle') && ...
           output.iterations >= 1 && isfield(outputps,'grad')
            grad = outputps.grad;
        elseif ~isa(g,'function_handle') && isempty(g)
            grad = serialize(deriv(f,z,fval));
        else
            grad = serialize(g(z));
        end
        s = p;
        y = grad-grad1;
        sy = real(y'*s);
        if sy > 0
            m = min(m+1,options.M);
            midx = [midx(1)+1:-1:1,m:-1:midx(1)];
            S(:,midx(1)) = s;
            Y(:,midx(1)) = y;
            r(:,midx(1)) = 1/sy;
        end
    end

    % Update the output structure.
    output.fval(end+1) = fval;
    output.gevals = output.gevals+1;
    output.iterations = output.iterations+1;
    output.relfval(end+1) = ...
        abs(diff(output.fval(end:-1:end-1)))/abs(output.fval(1));
    output.relstep(end+1) = relstep;
    if output.relfval(end) <= options.TolFun, output.info = 1; end
    if output.relstep(end) <= options.TolX, output.info = 2; end
    if output.iterations >= options.MaxIter, output.info = 3; end

    % Display progress.
    if options.Display > 0 && (output.iterations == 1 || output.info || ...
       mod(output.iterations,options.Display) == 0)
        if output.iterations == 1
            bold = '%s';
            [~,~,~,~,v] = regexp(version('-release'),'([0-9]+)([ab])');
            if usejava('Desktop') && str2double(v{1}{1}) > 2011 || ...
               (str2double(v{1}{1}) == 2011 && strcmpi(v{1}{2},'b'))
                bold = '%s';
            end
        end
        if output.iterations == 1 || ...
           mod(output.iterations,15*options.Display) == 0
            fprintf('\n%7s%s','',sprintf(bold,'fval'));
            fprintf('%13s%s','',sprintf(bold,'relfval'));
            fprintf('%10s%s','',sprintf(bold,'relstep'));
            if isa(options.PlaneSearch,'function_handle')
                fprintf('%10s%s','',sprintf(bold,'alpha'));
            else
                fprintf('%10s%s','',sprintf(bold,'delta'));
                fprintf('%8s%s','',sprintf(bold,'rho'));
            end
            fprintf('\n%21s%9s = %4.e %6s = %4.e\n\n','=1/2*norm(F)^2', ...
                    'TolFun',options.TolFun,'TolX',options.TolX);
        end
        if output.iterations == 1
            fprintf('%4i: % 14.8e |\n',0,output.fval(1));
        end
        if isa(options.PlaneSearch,'function_handle')
            stralpha = [repmat('%10.4e ',1,size(output.alpha,1)) '\n'];
            fprintf(['%4i: % 14.8e | %14.8e | %14.8e | ' stralpha], ...
                    output.iterations,output.fval(end), ...
                    output.relfval(end),output.relstep(end), ...
                    abs(output.alpha(:,end)));
        else
            fprintf(['%4i: % 14.8e | %14.8e | %14.8e | ' ...
                     '%10.4e | %10.4e\n'],...
                    output.iterations,output.fval(end), ...
                    output.relfval(end),output.relstep(end), ...
                    output.delta(end),output.rho(end));
        end
    end

end

% Display termination message.
if options.Display > 0
    ahref = '\n%s\n\n';
    x = round(linspace(0,output.iterations,min(500,output.iterations)));
    if length(bold) > 2
        ahref = sprintf(['\n%%s\n\n'],mat2str(x'), ...
            mat2str([output.fval(x+1)' [nan output.relfval(x(2:end))]' ...
                    [nan output.relstep(x(2:end))]'],3));
    end
    switch output.info
        case 1, fprintf(ahref,'Objective function tolerance reached.');
        case 2, fprintf(ahref,'Step size tolerance reached.');
        case 3, fprintf(ahref,'Maximum number of iterations reached.');
    end
end

end

function [z,offset] = deserialize(z,dim,offset)
    if iscell(dim)
        v = z;
        z = cell(size(dim));
        if nargin < 3, offset = 0; end
        for i = 1:numel(z)
            if iscell(dim{i})
                [z{i},offset] = deserialize(v,dim{i},offset);
            else
                n = prod(dim{i}(:));
                z{i} = reshape(v(offset+(1:n)),dim{i});
                offset = offset+n;
            end
        end
    elseif ~isempty(dim)
        z = reshape(z,dim);
    end
end

function z = serialize(z)
    if iscell(z)
        for i = find(cellfun(@iscell,z(:).'))
            z{i} = serialize(z{i});
        end
        s = cellfun(@numel,z(:)); o = [0; cumsum(s)];
        c = z; z = zeros(o(end),1);
        for i = 1:length(s), z(o(i)+(1:s(i))) = c{i}(:); end
    else
        z = z(:);
    end
end

function dim = structure(z)
    if iscell(z)
        dim = cellfun(@size,z,'UniformOutput',false);
        for i = find(cellfun(@iscell,z(:).'))
            dim{i} = structure(z{i});
        end
    else
        dim = size(z);
        if numel(z) == dim(1), dim = []; end
    end
end

author: andrewssobral

repo: lrslibrary