text:
function [Di_new,delet] = FDDL_UpdateDi (A,X,index,classn,Fish_ipts,Fish_par)
% ========================================================================
% Dictionary updating of FDDL, Version 1.0
% Copyright(c) 2011  Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang
% All Rights Reserved.
%
% -----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is here
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
%
% This is an implementation of the algorithm for updating the
% dictionary of FDDL (fix the coefficient)
%
% Please refer to the following paper
%
% Meng Yang, Lei Zhang, Xiangchu Feng, and David Zhang,"Fisher Discrimination
% Dictionary Learning for Sparse Representation", In IEEE Int. Conf. on
% Computer Vision, 2011.
% Meng Yang, Lei Zhang, Jian yang and David Zhang, "Metaface learning for
% sparse representation based face recognition", In ICIP, 2010.
%----------------------------------------------------------------------
%
%  Inputs :   (1) A :    the training data
%             (2) X :    the coefficient matrix of the training data
%             (3) index:   the label of the class being processed
%             (4) classn:  the total number of classes
%             (5) Fish_ipts
%                        . D  the dictionary in the last interation
%                        . trls  the labels of training data
%             (6) Fish_par
%                        .dls   the labels of the dicitonary atoms
%
% Outputs:    (1) Di_new :  the updated dictionary of the index-th class
%             (2) delet:    the indication of deleted atoms
%
%--------------------------------------------------------------------

D    =  Fish_ipts.D;
tau3 =  1;
tau2 =  1;
trls =  Fish_ipts.trls;
dls =  Fish_par.dls;
delet = [];

Do = D(:,dls~=index);
Di = D(:,dls==index);
Ai = A(:,trls==index);
Ao = A(:,trls~=index);
Xi = X(:,trls==index);
Xo = X(:,trls~=index);

Xi_i = Xi(dls==index,:);
Xi_o = Xi(dls~=index,:);
Xo_i = Xo(dls==index,:);
Xo_o = Xo(dls~=index,:);
% X_i  = X(dls==index,:);

Zi = Ai-Do*Xi_o;
Zo = Ao-Do*Xo_o;

for i = 1: size(Di,2)
    Yi = Zi - Di*Xi_i + Di(:,i)*Xi_i(i,:);
    Ui = Ai - Di*Xi_i + Di(:,i)*Xi_i(i,:);
%     Ua = A  - Di*X_i+Di(:,i)*X_i(i,:);
    Vo = Zo - Di*Xo_i + Di(:,i)*Xo_i(i,:);

    UaXti = zeros(size(Yi*(Xi_i(i,:))'));
    for t_i = 1:classn
        if t_i~=index
        Xt_i = X(dls==index,trls==t_i);
        UaXti = UaXti + (0  - Di*Xt_i + Di(:,i)*Xt_i(i,:))*(Xt_i(i,:))';
        end
    end

    tem1 = -Yi*(Xi_i(i,:))' - (tau2)*Ui*(Xi_i(i,:))' - Vo*(Xo_i(i,:))' - tau3*UaXti;
%     tem1 = - (tau2)*Ui*(Xi_i(i,:))';
%     tem1 = -Yi*(Xi_i(i,:))' - (tau2)*Ui*(Xi_i(i,:))'  - tau3*Ua*(Xo_i(i,:))';
    tem  = -tem1;
    if norm(tem,2)<1e-6
        Di(:,i) = zeros(size(tem));
        delet = [delet i];
    else
        Di(:,i) = tem./norm(tem,2);
    end
end

Di_new = Di;

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
%APPEND_PDFS Appends/concatenates multiple PDF files
%
% Example:
%   append_pdfs(output, input1, input2, ...)
%   append_pdfs(output, input_list{:})
%   append_pdfs test.pdf temp1.pdf temp2.pdf
%
% This function appends multiple PDF files to an existing PDF file, or
% concatenates them into a PDF file if the output file doesn't yet exist.
%
% This function requires that you have ghostscript installed on your
% system. Ghostscript can be downloaded from: http://www.ghostscript.com
%
% IN:
%    output - string of output file name (including the extension, .pdf).
%             If it exists it is appended to; if not, it is created.
%    input1 - string of an input file name (including the extension, .pdf).
%             All input files are appended in order.
%    input_list - cell array list of input file name strings. All input
%                 files are appended in order.

% Copyright: Oliver Woodford, 2011

% Thanks to Reinhard Knoll for pointing out that appending multiple pdfs in
% one go is much faster than appending them one at a time.

% Thanks to Michael Teo for reporting the issue of a too long command line.
% Issue resolved on 5/5/2011, by passing gs a command file.

% Thanks to Martin Wittmann for pointing out the quality issue when
% appending multiple bitmaps.
% Issue resolved (to best of my ability) 1/6/2011, using the prepress
% setting

function append_pdfs(varargin)
% Are we appending or creating a new file
append = exist(varargin{1}, 'file') == 2;
if append
    output = [tempname '.pdf'];
else
    output = varargin{1};
    varargin = varargin(2:end);
end
% Create the command file
cmdfile = [tempname '.txt'];
fh = fopen(cmdfile, 'w');
fprintf(fh, '-q -dNOPAUSE -dBATCH -sDEVICE=pdfwrite -dPDFSETTINGS=/prepress -sOutputFile="%s" -f', output);
fprintf(fh, ' "%s"', varargin{:});
fclose(fh);
% Call ghostscript
ghostscript(['@"' cmdfile '"']);
% Delete the command file
delete(cmdfile);
% Rename the file if needed
if append
    movefile(output, varargin{1});
end

author: opencobra

repo: cobratoolbox


text:
function M = fit_metric(fDir)

% fit_metric - fit a 2D Riemannian metric
%
%   M = fit_metric(fDir);
%
%   fDir(h) gives the value of the metric in direction h (a 1x2 vector).
%   M is a 2x2 symmetric matrix, that fit the measurements, i.e. such that
%       fDir(h)^2 ~ h*M*h'
%
%   Copyright (c) 2014 Gabriel Peyre

a = fDir([1 0])^2;
b = fDir([0 1])^2;
c = (fDir([1 1])^2-a-b)/2;
M = [a c; c b];

[U,S] = eig(M); S = diag(S);
if min(S)<0
    warning('Non positive matrix, clamping singular values');
    M = U*diag(max(S,0))*U';
end

author: gpeyre


text:
function [block, mfg, mbg] = amrf_seg(im, blocklen, overlap)
% AMRF_SEG  Adaptive Markov Random Field segmentation.
%
%   Note: This function is in development, and can produce "blocky"
%   segmentations.
%
%   This function splits an image into blocks, estimates typical background
%   and foreground intensities, and uses them to seed a Markov Random Field
%   segmentation algorithm from the Insight Toolbox (itk::MRFImageFilter).
%
% BW = amrf_seg(IM, BLOCKLEN)
%
%   IM is a 3D array with grayscale intesity values, forming a 3D image
%   (e.g. Magnetic Resonance). The background is assumed to be light
%   (higher intensity values), and the foreground dark (lower intensity
%   values).
%
%   IM is split into blocks to estimate the typical foreground and
%   background intensity values.
%
%   BLOCKLEN is a 3-vector with the number of rows, columns and slices in
%   each block. If an integer number of blocks does not fit in the image,
%   the last blocks in each dimension are cropped.
%
%   BW is a binary segmentation of IM. Foreground voxels = 1, and
%   background voxels = 0.
%
% [..., MFG, MBG] = amrf_seg(..., OVERLAP)
%
%   To improve the smoothness of the segmentation result, we don't want
%   sudden changes in the MUFG, MUGB between neighbouring blocks.
%
%   OVERLAP is a 3-vector with the number of rows, columns and slices added
%   on each side of each block so that it overlaps with its neighbours. By
%   default, OVERLAP = [0 0 0] (no overlap is used). Note that large
%   overlaps will increase the memory used by this function substantially.
%
%   MFG, MBG are 3D arrays with the typical (modes) foreground and
%   background intensities for each of the estimation blocks.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2014 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2, 3);
nargoutchk(0, 3);

% defaults
if (nargin < 3 || isempty(overlap))
    overlap = [0 0 0];
end

% image dimensions
NR = size(im, 1);
NC = size(im, 2);
NS = size(im, 3);

%% Split image into overlapping blocks and estimate foreground and
%% background

% "botton left" index of every block in the image, without overlap. If the
% last block starts at the right edge of the image, it has size zero. Thus,
% if we have any block starting there, we remove it with setdiff()
br0 = setdiff(1:blocklen(1):NR, NR);
bc0 = setdiff(1:blocklen(2):NC, NC);
bs0 = setdiff(1:blocklen(3):NS, NS);

% "top right" index of every block in the image, without overlap (taking
% care that the last boxes don't go outside the image)
tr0 = min(NR, br0 + blocklen(1) - 1);
tc0 = min(NC, bc0 + blocklen(2) - 1);
ts0 = min(NS, bs0 + blocklen(3) - 1);

% the last blocks have to cover to the end of the image
tr0(end) = NR;
tc0(end) = NC;
ts0(end) = NS;

% extend the blocks with the overlap margins
br = max(1, br0 - overlap(1));
bc = max(1, bc0 - overlap(2));
bs = max(1, bs0 - overlap(3));
tr = min(NR, tr0 + overlap(1));
tc = min(NC, tc0 + overlap(2));
ts = min(NS, ts0 + overlap(3));

% split the image into blocks. This is necessary if we want to process the
% blocks in parallel. We create a vector array of blocks instead of a 3D
% array of blocks because parfor cannot be nested
block = cell(1, length(br) * length(bc) * length(bs));
for ri = 1:length(br)
    for ci = 1:length(bc)
        for si = 1:length(bs)

            block{sub2ind([length(br) length(bc) length(bs)], ri, ci, si)} ...
                = im(br(ri):tr(ri), bc(ci):tc(ci), bs(si):ts(si));

        end
    end
end

% parallel processing of the blocks
mfg = zeros(length(br), length(bc), length(bs));
mbg = zeros(length(br), length(bc), length(bs));
parfor I = 1:numel(block)

    [mfg(I), mbg(I)] = im_modes(block{I});

end

% clear memory
clear block

% reshape outputs as 3D arrays
mfg = reshape(mfg, length(br), length(bc), length(bs));
mbg = reshape(mbg, length(br), length(bc), length(bs));

%% Interpolate typical values of foreground and background intensity that
%% couldn't be estimated in the previous phase

% values that need to be interpolated
fgidx = isnan(mfg);
bgidx = isnan(mbg);

if (nnz(fgidx) || nnz(bgidx))

    % centers of each block without overlap
    [gr0, gc0, gs0] ...
        = ndgrid((br0 + tr0)*0.5, (bc0 + tc0)*0.5, (bs0 + ts0)*0.5);

    % scale the center coordinates to be within the [0, 1] cube, to avoid
    % ill-conditioned L matrices in the thin-plate spline. Note that
    % because of the scale-invariance of the TPS, this doesn't change the
    % result
    K = max([(br0(end) + tr0(end))*0.5, ...
        (bc0(end) + tc0(end))*0.5, ...
        (bs0(end) + ts0(end))*0.5]);
    gr0 = gr0 / K;
    gc0 = gc0 / K;
    gs0 = gs0 / K;

    % interpolate foreground missing values
    mfg(fgidx) = pts_tps_map([gr0(~fgidx), gc0(~fgidx), gs0(~fgidx)], ...
        mfg(~fgidx), [gr0(fgidx), gc0(fgidx), gs0(fgidx)]);

    % interpolate background missing values
    mbg(bgidx) = pts_tps_map([gr0(~bgidx), gc0(~bgidx), gs0(~bgidx)], ...
        mbg(~bgidx), [gr0(bgidx), gc0(bgidx), gs0(bgidx)]);

end

%% Split the image into non-overlapping blocks, and segment them using a
%% Markov Random Field method

% split the image into blocks, this time without overlap
block = cell(1, length(br) * length(bc) * length(bs));
for ri = 1:length(br0)
    for ci = 1:length(bc0)
        for si = 1:length(bs0)

            block{sub2ind([length(br0) length(bc0) length(bs0)], ri, ci, si)} ...
                = im(br0(ri):tr0(ri), bc0(ci):tc0(ci), bs0(si):ts0(si));

        end
    end
end

% Markov Random Field segmentation of the blocks
parfor I = 1:numel(block)

    block{I} = itk_imfilter('mrf', block{I}, [mfg(I) mbg(I)]);

end

% reshape the vector of blocks as 3D array
block = reshape(block, length(br), length(bc), length(bs));

% format output as segmentation, and invert so that tissue=1, background=0
block = ~cell2mat(block);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%%***********************************************************************
%% validate: validate data
%%
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%***********************************************************************

   function  [blk,At,C,b,dim,nnblk,parbarrier] = ...
              validate(blk,At,C,b,par,parbarrier);

   if (nargin >= 5)
      spdensity = par.spdensity;
   else
      spdensity = 0.4;
   end
%%
   if ~iscell(blk);
      error('validate: blk must be a cell array'); end;
   if (size(blk,2) < 2)
      error('validate: blk must be a cell array with at least 2 columns');
   end
   if ~iscell(At) | ~iscell(C);
      error('validate: At, C must be cell arrays');
   end
   if (size(At,1) ~= size(blk,1))
      if (size(At,2) == size(blk,1));
         At = At';
      else
         error('validate: size of At is not compatible with blk');
      end
   end
   if (size(C,1) ~= size(blk,1))
      if (size(C,2) == size(blk,1))
         C = C';
      else
         error('validate: size of C is not compatible with blk');
      end
   end
   if (min(size(b)) > 1); error('validate: b must be a vector'); end
   if (size(b,1) < size(b,2)); b = b'; end
   m = length(b);
%%
%%-----------------------------------------
%% validate blk, At, C
%%-----------------------------------------
%%
   for p = 1:size(blk,1)
      if (size(blk{p,2},1) > size(blk{p,2},2))
         blk{p,2} = blk{p,2}';
      end
      pblk = blk(p,:);
      n = sum(pblk{2});
      numblk = length(pblk{2});
      if strcmp(pblk{1},'s');
         m1 = size(At{p,1},2);
         n2 = sum(pblk{2}.*pblk{2});  n22 = sum(pblk{2}.*(pblk{2}+1))/2;
         ntotal(p) = n22;
         if ~all(size(C{p}) == n)
            error('validate: blk and C are not compatible');
         end
         if (norm(C{p}-C{p}',inf) > 1e-13*norm(C{p},inf));
            error('validate: C is not symmetric');
         end
         if (size(At{p,1}) == [m1, n22] & m1~=n22);
            At{p,1} = At{p,1}';
         end
         if (~isempty(At{p,1})) & (size(At{p,1},1) ~= n22)
            error('validate: blk and At not compatible');
         end
         if (nnz(At{p,1}) < spdensity*n22*m1)
            if ~issparse(At{p,1}); At{p,1} = sparse(At{p,1}); end
         end
         if (length(pblk) > 2) %% for low rank constraints
            len = sum(pblk{3});
            if (size(pblk{1,3},2) < size(pblk{1,3},1))
                blk{p,3} = blk{p,3}';
            end
            if any(size(At{p,2})- [n,len])
               error(' low rank structure specified in blk and At not compatible')
            end
            if (length(At(p,:)) > 2)  & ~isempty(At{p,3})
               if all(size(At{p,3},2)-[1,4])
                  error(' low rank structure in At{p,3} not specified correctly')
               end
               if (size(At{p,3},2) == 1)
                  if (size(At{p,3},1) < size(At{p,3},2)); At{p,3} = At{p,3}'; end
                  lenn = length(At{p,3});
                  constrnum = mexexpand(pblk{3},[1:length(pblk{3})]');
                  At{p,3} = [constrnum, [1:lenn]', [1:lenn]', At{p,3}];
               elseif (size(At{p,3},2) == 4)
                  dd = At{p,3};
                  [dummy,idxsort] = sort(dd(:,1));
                  dd = dd(idxsort,:);
                  lenn = size(dd,1);
                  idxD = [0; find(diff(dd(:,1))); lenn];
                  ii = zeros(lenn,1); jj = zeros(lenn,1);
                  ss = [0,cumsum(pblk{3})];
                  for k = 1:length(pblk{3})
                     idx = [idxD(k)+1 : idxD(k+1)];
                     ii(idx) = dd(idx,2)+ss(k); %% convert to cumulative indexing
                     jj(idx) = dd(idx,3)+ss(k);
                  end
                  At{p,3} = [dd(:,1),ii,jj,dd(:,4)];
               end
            else
               constrnum = mexexpand(pblk{3},[1:length(pblk{3})]');
               At{p,3} = [constrnum,  [1:len]', [1:len]', ones(len,1)];
            end
         end
         if (nnz(C{p}) < spdensity*n2) | (numblk > 1);
            if ~issparse(C{p}); C{p} = sparse(C{p}); end;
         else
            if issparse(C{p}); C{p} = full(C{p}); end;
         end
      elseif strcmp(pblk{1},'q') | strcmp(pblk{1},'l') | strcmp(pblk{1},'u');
         ntotal(p) = n;
         if (min(size(C{p})) ~= 1 | max(size(C{p})) ~= n);
            error(['validate: blk and C are not compatible']);
         end;
         if (size(C{p},1) < size(C{p},2)); C{p} = C{p}'; end
         if (size(At{p,1}) == [m, n] & m~=n);
            At{p,1} = At{p,1}';
         end
         if ~all(size(At{p,1}) == [n,m]);
            error('validate: blk and At not compatible');
         end
         if ~issparse(At{p,1});
            At{p,1} = sparse(At{p,1});
         end
         if (nnz(C{p}) < spdensity*n);
            if ~issparse(C{p}); C{p} = sparse(C{p}); end;
         else
            if issparse(C{p}); C{p} = full(C{p}); end;
         end;
      else
         error(' blk: some fields are not specified correctly');
      end
   end
   if (sum(ntotal) < m)
      error(' total dimension of C should be > length(b)');
   end
%%
%%-----------------------------------------
%% problem dimension
%%-----------------------------------------
%%
   dim = zeros(1,4);  nnblk = zeros(1,2);
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      if strcmp(pblk{1},'s')
         dim(1) = dim(1) + sum(pblk{2});
         nnblk(1) = nnblk(1) + length(pblk{2});
         nn(p) = sum(pblk{2});
      elseif strcmp(pblk{1},'q')
         dim(2) = dim(2) + sum(pblk{2});
         nnblk(2) = nnblk(2) + length(pblk{2});
         nn(p) = length(pblk{2});
      elseif strcmp(pblk{1},'l')
         dim(3) = dim(3) + sum(pblk{2});
         nn(p) = sum(pblk{2});
      elseif strcmp(pblk{1},'u')
         dim(4) = dim(4) + sum(pblk{2});
         nn(p) = sum(pblk{2});
      end
   end
%%
%%-----------------------------------------
%% validate parbarrier
%%-----------------------------------------
%%
    if (nargin == 6)
       if ~iscell(parbarrier);
	  if (length(parbarrier) == size(blk,1))
	     tmp = parbarrier;
             clear parbarrier;
	     parbarrier = cell(size(blk,1),1);
	     for p = 1:size(blk,1)
		parbarrier{p} = tmp(p);
             end
          end
       end
       if (size(parbarrier,2) > size(parbarrier,1))
          parbarrier = parbarrier';
       end
       for p = 1:size(blk,1)
          pblk = blk(p,:);
          if (size(parbarrier{p},1) > size(parbarrier{p},2))
             parbarrier{p} = parbarrier{p}';
          end
          len = length(parbarrier{p});
          if strcmp(pblk{1},'s') | strcmp(pblk{1},'q')
             if (len == 1)
                parbarrier{p} = parbarrier{p}*ones(1,length(pblk{2}));
	     elseif (len == 0)
                parbarrier{p} = zeros(1,length(pblk{2}));
	     elseif (len ~= length(pblk{2}))
                error('blk and parbarrier not compatible');
             end
          elseif strcmp(pblk{1},'l')
             if (len == 1)
                parbarrier{p} = parbarrier{p}*ones(1,sum(pblk{2}));
	     elseif (len == 0)
		parbarrier{p} = zeros(1,sum(pblk{2}));
             elseif (len ~= sum(pblk{2}))
                error('blk and parbarrier not compatible');
             end
          elseif strcmp(pblk{1},'u')
             parbarrier{p}= zeros(1,sum(pblk{2}));
          end
       end
    end
%%***********************************************************************

repo: REST


text:
function ch = digit_oct_to_ch ( i )

%*****************************************************************************80
%
%% DIGIT_OCT_TO_CH returns the character representation of an octal digit.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the integer, between 0 and 7.
%
%    Output, character CH, the octal representation of the integer.
%
  if ( 0 <= i && i <= 7 )
    ch = char ( '0' + i );
  else
    ch = '?';
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = inputindexGUI(varargin)

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @inputindexGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @inputindexGUI_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 inputindexGUI is made visible.
function inputindexGUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = [];

try
      prompt    = char(varargin{1});
      num_lines = varargin{2};
      def       = varargin{3};
catch
      prompt    = 'Enter value';
      num_lines = 1;
      def       = 1;
end

set(handles.uipanel_prompt,'Title', prompt)
set(handles.popupmenu_value,'String', cellstr(num2str(([1:num_lines]'))))
set(handles.popupmenu_value,'Value', def)

%
% Color GUI
%
handles = painterplab(handles);

%
% Set font size
%
handles = setfonterplab(handles);

% Update handles structure
guidata(hObject, handles);

% help
% helpbutton

% UIWAIT makes inputindexGUI wait for user response (see UIRESUME)
uiwait(handles.gui_chassis);

% -----------------------------------------------------------------------
function varargout = inputindexGUI_OutputFcn(hObject, eventdata, handles)

varargout{1} = handles.output;
% The figure can be deleted now
delete(handles.gui_chassis);
pause(0.1)

% -----------------------------------------------------------------------
function pushbutton_CANCEL_Callback(hObject, eventdata, handles)
handles.output = [];
% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

% -----------------------------------------------------------------------
function pushbutton_OK_Callback(hObject, eventdata, handles)

value = get(handles.popupmenu_value,'Value'); % JLC, 30 Aug 2012
handles.output = {value};

% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);
% -----------------------------------------------------------------------
function popupmenu_value_Callback(hObject, eventdata, handles)

% -----------------------------------------------------------------------
function popupmenu_value_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% -----------------------------------------------------------------------
function uipanel_prompt_CreateFcn(hObject, eventdata, handles)

% -----------------------------------------------------------------------
function gui_chassis_CloseRequestFcn(hObject, eventdata, handles)

if isequal(get(handles.gui_chassis, 'waitstatus'), 'waiting')
        %The GUI is still in UIWAIT, us UIRESUME
        handles.output = '';
        %Update handles structure
        guidata(hObject, handles);
        uiresume(handles.gui_chassis);
else
        % The GUI is no longer waiting, just close it
        delete(handles.gui_chassis);
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function [F,FF] = icml2011_mohsst5_parse_results(filename, N, burnin)

F_mean = 0;
FF = 0;

for n=(burnin+1):N
  fprintf('%d/%d\n', n, N);
  name = sprintf('%s_F%d',filename,n);
  load(name);
  F_mean = (F + (n-burnin-1)*F_mean) / (n-burnin);
  FF = (F.*F + (n-burnin-1)*FF) / (n-burnin);
end

F = F_mean;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function wfiltdtinfo(dw,varargin)
%WFILTDTINFO Plots dual-tree filters info
%   Usage: wfiltdtinfo(dw);
%
%   Input parameters:
%         dw     : Wavelet dual-tree filterbank
%
%   `wfiltdtinfo(w)` plots impulse responses, frequency responses and
%   approximation of the scaling and of the wavelet function(s) associated
%   with the dual-tree wavelet filters defined by *w* in a single figure.
%   Format of *dw* is the same as in |dtwfb|.
%
%   The figure is organized as follows:
%
%   First row shows impulse responses of the first (real) tree.
%
%   Second row shows impulse responses of the second (imag) tree.
%
%   Third row contains plots of real (green), imaginary (red) and absolute
%   (blue) parts of approximation of scaling and wavelet function(s).
%
%   Fourth and fifth row show magnitude and phase frequency responses
%   respectivelly of filters from rows 1 and 2 with matching colors.
%
%   Optionally it is possible to define scaling of the y axis of the
%   frequency seponses. Supported are:
%
%   'db','lin'
%       dB or linear scale respectivelly. By deault a dB scale is used.
%
%   Examples:
%   ---------
%   :::
%      wfiltdtinfo('qshift4');
%

% AUTHOR: Zdenek Prusa

complainif_notenoughargs(nargin,1,'WFILTDTINFO');

definput.flags.freqzscale = {'db','lin'};
[flags]=ltfatarghelper({},definput,varargin);

[dwstruct,info] = dtwfbinit({'strict',{dw,6,'dwt'}});
dw = info.dw;

filtNo = size(dw.g,1);
grayLevel = [0.6,0.6,0.6];
clf;

colorAr ={repmat('rbk',1,filtNo),repmat('cmg',1,filtNo)};

for ii=1:2
    subplot(5,filtNo,filtNo*(ii-1)+1);
    title(sprintf('Scaling imp. response, tree %i',ii));
    loAna = dw.g{1,ii}.h;
    loShift = -dw.g{1,ii}.offset;
    xvals = -loShift + (0:length(loAna)-1);
    hold on;
    if ~isempty(loAna(loAna==0))
       stem(xvals(loAna==0),loAna(loAna==0),'Color',grayLevel);
    end

    loAnaNZ = find(loAna~=0);
    stem(xvals(loAnaNZ),loAna(loAnaNZ),colorAr{ii}(1));
    axis tight;
    hold off;
end

for ii=1:2
    for ff=2:filtNo
        subplot(5,filtNo,ff+filtNo*(ii-1));
        title(sprintf('Wavelet imp. response no: %i, tree %i',ff-1,ii));
        filtAna = dw.g{ff,ii}.h;
        filtShift = -dw.g{ff,ii}.offset;
        xvals = -filtShift + (0:length(filtAna)-1);
        filtNZ = find(filtAna~=0);
        hold on;

        if ~isempty(filtAna(filtAna==0))
           stem(xvals(filtAna==0),filtAna(filtAna==0),'Color',grayLevel);
        end

        stem(xvals(filtNZ),filtAna(filtNZ),colorAr{ii}(ff));
        axis tight;
        hold off;
    end
end

L =  wfbtlength(1024,dwstruct,'per');
Lc = wfbtclength(L,dwstruct,'per');
c = wavpack2cell(zeros(sum([Lc;Lc(end:-1:1)]),1),...
                 [Lc;Lc(end:-1:1)]);
c{1}(1) = 1;
sfn = idtwfb(c,dwstruct,L);

subplot(5,filtNo,[2*filtNo+1]);
xvals = ((-floor(L/2)+1:floor(L/2)).');
plot(xvals,fftshift([abs(sfn),real(sfn),imag(sfn)],1));
axis tight;
title('Scaling function');

%legend({'abs','real','imag'},'Location','south','Orientation','horizontal')

for ff=2:filtNo
   subplot(5,filtNo,[2*filtNo+ff]);

   c{ff-1}(1) = 0;
   c{ff}(1) = 1;
   wfn = idtwfb(c,dwstruct,L);

   plot(xvals,fftshift([abs(wfn),real(wfn),imag(wfn)],1));
   axis tight;
   %legend({'abs','real','imag'},'Location','south','Orientation','horizontal')
   title(sprintf('Wavelet function: %i',ff-1));
end

subplot(5,filtNo,3*filtNo + (1:filtNo) );
title('Magnitude frequency response');
maxLen=max(cellfun(@(gEl) numel(gEl.h),dw.g));
Ls = nextfastfft(max([maxLen,1024]));
H = filterbankfreqz(dw.g(:),[dw.a(:);dw.a(:)],Ls);

%[H] = wtfftfreqz(w.g);
if flags.do_db
   plotH = 20*log10(abs(H));
elseif flags.do_lin
   plotH = abs(H);
else
    error('%s: Unknown parameter',upper(mfilaname));
end
xVals = linspace(0,1,numel(H(:,1)));
hold on;
for ii=1:2
for ff=1:filtNo
   plot(xVals,plotH(:,ff+(ii-1)*filtNo),colorAr{ii}(ff));
   axis tight;
end
end
if flags.do_db
    ylim([-30,max(plotH(:))])
end
ylabel('|\itH|[dB]');
xlabel('\omega [-]')
hold off;

subplot(5,filtNo,4*filtNo + (1:filtNo) );
title('Phase frequency response');
hold on;
for ii=1:2
for ff=1:filtNo
   plot(xVals,unwrap(angle((H(:,ff+(ii-1)*filtNo))))/pi,colorAr{ii}(ff));
   axis tight;
end
end
ylabel('arg H(\omega)[\pi rad]');
xlabel('\omega [-]')

axis tight;
% plot(unwrap(angle([H])));
% axis tight;
hold off;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%DOUBLE Dataset to double conversion
%
%	D = double(A)
%
% Converts a dataset A to a double D, which is just the set of datavectors.

author: marianux

repo: ecg-kit

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


text:
function pass = test_sum( pref )
% Test with function cos(cos(lam)*sin(th))

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e4*pref.techPrefs.chebfuneps;

%% Integrate over r

% Example 1
f = ballfun(@(r,lam,th)r.*cos(lam).*sin(th), 'spherical');
g = sum(f, 1);
exact = spherefun(@(lam,th)cos(lam).*sin(th)/4);
pass(1) = norm(g-exact) < tol;

% Example 2
f = ballfun(@(r,lam,th)1, 'spherical');
g = sum(f, 1);
exact = spherefun(@(lam,th)1/3);
pass(2) = norm(g-exact) < tol;

%% Integrate over lambda

% Example 3
f = ballfun(@(r,lam,th)(r.*sin(lam).*sin(th)).^2, 'spherical');
g = sum(f, 2);
exact = diskfun(@(th,r)pi*r.^2.*sin(th).^2,'polar');
pass(3) = norm(g-exact) < tol;

%% Integrate over theta

% Example 4
f = ballfun(@(r,lam,th)r.*cos(lam).*sin(th), 'spherical');
g = sum(f, 3);
exact = diskfun(@(lam,r)r.*cos(lam)*pi/2,'polar');
pass(4) = norm(g-exact) < tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function seed = van_der_corput_seed_get ( )

%*****************************************************************************80
%
%% VAN_DER_CORPUT_SEED_GET gets the seed of the van der Corput sequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Johannes van der Corput,
%    Verteilungsfunktionen I & II,
%    Nederl. Akad. Wetensch. Proc.,
%    Volume 38, 1935, pages 813-820, pages 1058-1066.
%
%  Parameters:
%
%    Output, integer SEED, the current value for the seed.
%
%    Global, integer VAN_DER_CORPUT_SEED, the current seed.
%
  global van_der_corput_BASE
  global van_der_corput_SEED
%
 seed = van_der_corput_SEED;

  return
end


text:
function d = p09_fd ( p )

%*****************************************************************************80
%
%% P09_FD is a signed distance function for problem 9.
%
%  Discussion:
%
%    Notice that the V1 and V2 arrays, which define hexagons, contain SEVEN
%    points, not six.  This is necessary in order for DPOLY to carry out its
%    calculations correctly.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    11 March 2006
%
%  Parameters:
%
%    Input, real P, one or more points.
%
%    Output, real D, the signed distance of each point to the boundary of the region.
%
  v1 = [ 0.15, 0.75;
         0.20, 0.6634;
         0.30, 0.6634;
         0.35, 0.75;
         0.30, 0.8366;
         0.20, 0.8366;
         0.15, 0.75 ];

  d1 = dpoly ( p, v1 );

  v2 = [ 0.50, 0.40;
         0.55, 0.3134;
         0.65, 0.3134;
         0.70, 0.40;
         0.65, 0.4866;
         0.55, 0.4866;
         0.50, 0.40 ];

  d2 = dpoly ( p, v2 );

  d3 = dunion ( d1, d2 );

  d4 = drectangle ( p, 0.0, 1.0, 0.0, 1.0 );

  d = ddiff ( d4, d3 );

  return
end


text:
function fx = p34_fun ( n, x )

%*****************************************************************************80
%
%% P34_FUN evaluates the integrand for problem 34.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    ( 10 * x - 1 ) * ( 10 * x - 1.1 ) * ( 10 * x - 1.2 ) * ( 10 * x - 1.3 )
%
%  Exact Integral:
%
%    1627879 / 1500
%
%  Approximate Integral (20 digits):
%
%    1085.2526666666666666...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Hermann Engels,
%    Numerical Quadrature and Cubature,
%    Academic Press, 1980.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  fx = ( 10.0 * x - 1.0 ) .* ( 10.0 * x - 1.1 ) .* ( 10.0 * x - 1.2 ) ...
    .* ( 10.0 * x - 1.3 );

  return
end


text:
function [y,mm]=v_momfilt(x,r,w,m)
%V_MOMFILT calculates moments of a signal using a sliding window Y=(X,R,W,M)
%
% Inputs: x    is the input signal
%         r    is a list of moments to calculate
%              (+ve = relative to mean, -ve = relative to zero)
%         w    is the window (or just the length of a Hamming window)
%              Note: If the window is asymmetric, you should be aware that it gets
%              flipped in the convolution process
%         m    is the sample of w to use as the centre [default=ceil(length(w)/2+0.5)]
%
%         mm   the actual value of m used. Output point y(i) is based on x(i+m-w:i+m-1).
%
% Example:
%  To calculate a running kurtosis using a Hamming window of length 30:
%             y=v_momfilt(x,[2 4],30); k=y(:,2)./y(:,1).^2

%	   Copyright (C) Mike Brookes 2007
%      Version: $Id: v_momfilt.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin < 3
   w=hamming(length(x));
elseif  length(w)==1
   w=hamming(w);
end
lw=length(w);
w=w(:);                 % force to a column vector
if nargin < 4
   m=(1+lw)/2;
end
m=max(round(m),1);
mm=m;
r=round(r(:))';             % force integer row vector of moments

lx=prod(size(x));
lxx=lx+m-1;
xx=zeros(lxx,1);
xx(1:lx)=x(:);              % extend with zeros so filter() works correctly

cw=cumsum(w);
sw=cw(end);
y0=repmat(sw,lxx,1);
lxw=min(lxx,lw);
y0(1:lxw)=cw(1:lxw);
y0(lx+1:lx+m-1)=y0(lx+1:lx+m-1)-cw(1:m-1);      % equivalent to y0=filter(w,1,xx^0);
yd=y0(m:end);
yd(abs(yd)<eps)=1;

nr=length(r);
wlx=ones(lx,1);
wlxx=ones(lxx,1);
y=zeros(lx,nr);
mr=max(abs(r));                         % max moment to calculate
mk=zeros(1,mr);                       % list of moments
mk(-r(r<0))=1;                         % choose the moments we need to calculate
maxr=max(r);
if maxr>1
    mk(1:maxr)=1;
end
ml=find(mk>0);
lml=length(ml);
if lml
    mlx=mk;
    mlx(ml)=1:lml;                        % mapping from moment into ml
    wlml=ones(1,lml);
    xm=filter(w,1,xx(:,wlml).^ml(wlxx,:));    % calculate all the moments
    xm=xm(m:end,:)./yd(:,wlml);                             % remove the useless start values and normalize
end
fr=find(r<0);
if length(fr)
    y(:,fr)=xm(:,mlx(-r(fr)));    % zero-centred moments
end
fr=find(r==0);                              % 0'th moment
if length(fr)
    y(:,fr)=1;
end
fr=find(r==1);                              % 1'st moment about mean
if length(fr)
    y(:,fr)=0;
end
fr=find(r==2);                              % 2'nd moment about mean
if length(fr)
    yfr=xm(:,2)-xm(:,1).^2;
    y(:,fr)=yfr(:,ones(1,length(fr)));
end
if maxr>2
    mon=[1 -1];
    bc=[1 -2 1];
    am=zeros(lx,maxr);
    am(:,1)=xm(:,1);                            % copy the mean across
    ml=2:maxr;
    wlml=ones(1,length(ml));
    am(:,2:end)=xm(:,wlml).^ml(ones(lx,1),:);              % calculate powers of the mean
    for i=3:maxr
        bc=conv(bc,mon);                        % calculate binomial coefficients
        fr=find(r==i);
        if length(fr)
            yfr=xm(:,i)+sum(xm(:,i-1:-1:2).*am(:,1:i-2).*bc(wlx,2:i-1),2)+am(:,i)*(bc(i)+bc(i+1));
            y(:,fr)=yfr(:,ones(1,length(fr)));
        end
    end
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
classdef (ConstructOnLoad, Abstract)ZmapCatalogAddon < matlab.mixin.Copyable
    % ZMAPCATALOGADDON is a base class for add-ons to a zmap catalog. For example, Moment Tensors
    properties(Abstract, Constant)
        Type
    end

    methods
        function obj = ZmapCatalogAddon()
        end

        function subsetInPlace(obj, idx)
            pef = obj.possibly_empty_fields();
            for ii =1: numel(pef)
                fn = pef{ii};
                if ~isempty(obj.(fn))
                    obj.(fn) = obj.(fn)(idx,:);
                end
            end
        end

        function newobj = subset(obj, idx)
            newobj=copy(obj);
            newobj.subsetInPlace(idx)
        end

        function disp(obj)
            display_helper(obj)
        end

    end
    methods (Abstract, Static, Hidden)
    	pef = possibly_empty_fields()
        s = display_order()
    end
end

author: CelsoReyes


text:
%%*****************************************************************************
%% HSDsqlp: solve an semidefinite-quadratic-linear program
%%   by infeasible path-following method on the homogeneous self-dual model.
%%
%%  [obj,X,y,Z,info,runhist] =
%%      HSDsqlp(blk,At,C,b,OPTIONS,X0,y0,Z0,kap0,tau0,theta0);
%%
%%  Input: blk: a cell array describing the block diagonal structure of SQL data.
%%          At: a cell array with At{p} = [svec(Ap1) ... svec(Apm)]
%%         b,C: data for the SQL instance.
%%     OPTIONS: a structure that specifies parameters required in HSDsqlp.m,
%%              (if it is not given, the default in sqlparameters.m is used).
%%
%%  (X0,y0,Z0): an initial iterate (if it is not given, the default is used).
%%  (kap0,tau0,theta0): initial parameters (if not given, the default is used).
%%
%%  Output: obj  = [<C,X> <b,y>].
%%          (X,y,Z): an approximately optimal solution or a primal or dual
%%                   infeasibility certificate.
%%          info.termcode = termination-code
%%          info.iter     = number of iterations
%%          info.obj      = [primal-obj, dual-obj]
%%          info.cputime  = total-time
%%          info.gap      = gap
%%          info.pinfeas  = primal_infeas
%%          info.dinfeas  = dual_infeas
%%          runhist.pobj    = history of primal objective value.
%%          runhist.dobj    = history of dual   objective value.
%%          runhist.gap     = history of <X,Z>.
%%          runhist.pinfeas = history of primal infeasibility.
%%          runhist.dinfeas = history of dual   infeasibility.
%%          runhist.cputime = history of cputime spent.
%%----------------------------------------------------------------------------
%%  The OPTIONS structure specifies the required parameters:
%%      vers  gam  predcorr  expon  gaptol  inftol  steptol
%%      maxit  printlevel  ...
%%      (all have default values set in sqlparameters.m).
%%
%%*************************************************************************
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*************************************************************************

  function [obj,X,y,Z,info,runhist] = ...
            HSDsqlpmain(blk,At,C,b,par,X0,y0,Z0,kap0,tau0,theta0);

%%
%%-----------------------------------------
%% get parameters from the OPTIONS structure.
%%-----------------------------------------
%%
   global spdensity  solve_ok  printlevel
   global schurfun   schurfun_par
%%
   randstate = rand('state');  randnstate = randn('state');
   rand('state',0);   randn('state',0);
%%
   matlabversion = par.matlabversion;
   vers          = par.vers;
   predcorr      = par.predcorr;
   gam           = par.gam;
   expon         = par.expon;
   gaptol        = par.gaptol;
   inftol        = par.inftol;
   steptol       = par.steptol;
   maxit         = par.maxit;
   printlevel    = par.printlevel;
   stoplevel     = par.stoplevel;
   scale_data    = par.scale_data;
   spdensity     = par.spdensity;
   rmdepconstr   = par.rmdepconstr;
   cachesize     = par.cachesize;
   smallblkdim   = par.smallblkdim;
   schurfun      = par.schurfun;
   schurfun_par  = par.schurfun_par;
   ublksize      = par.ublksize;
%%
   tstart = cputime;
   X = X0; y = y0; Z = Z0;
   for p = 1:size(blk,1)
      if strcmp(blk{p,1},'u'); Z{p} = zeros(blk{p,2},1); end
   end
%%
%%-----------------------------------------
%% convert unrestricted blk to linear blk.
%%-----------------------------------------
%%
   ublkidx = zeros(size(blk,1),1);
   Cpert = zeros(size(blk,1),1); Cnew = C;
   perturb_C = 1;
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      n = sum(pblk{2});
      tmp = max(1,norm(C{p},'fro'))/sqrt(n);
      if strcmp(pblk{1},'u')
         if (printlevel); fprintf(' *** convert ublk to linear blk'); end
         ublkidx(p) = 1;
         n = 2*pblk{2};
         blk{p,1} = 'l'; blk{p,2} = n;
         if (perturb_C); Cpert(p) = 1e-2*tmp; end
         C{p}  = [C{p}; -C{p}];
         At{p} = [At{p}; -At{p}];
         Cnew{p} = C{p} + Cpert(p)*ones(n,1);
         X{p}  = 1+rand(n,1); %% do not add a factor of n
         Z{p}  = 1+rand(n,1); %%
      elseif strcmp(pblk{1},'s')
         if (perturb_C); Cpert(p) = 1e-3*tmp; end
         Cnew{p} = C{p} + Cpert(p)*speye(n);
      else
         if (perturb_C); Cpert(p) = 1e-3*tmp; end
         Cnew{p} = C{p} + Cpert(p)*ones(n,1);
      end
   end
%%
%%-----------------------------------------
%% check if the matrices Ak are
%% linearly independent.
%%-----------------------------------------
%%
   m0 = length(b);
   [At,b,y,indeprows,depconstr,feasible,AAt] = ...
   checkdepconstr(blk,At,b,y,rmdepconstr);
   if (~feasible)
      msg = 'SQLP is not feasible';
      if (printlevel); fprintf('\n %s',msg); end
      return;
   end
   par.depconstr = depconstr;
%%
   normC = zeros(length(C),1);
   for p = 1:length(C); normC(p) = max(max(abs(C{p}))); end
   normC = 1+max(normC);
   normb = 1+max(abs(b));
   nn = ops(C,'getM');
   m = length(b);
   if (nargin <= 8) | (isempty(kap0) | isempty(tau0) | isempty(theta0))
      if (max([ops(At,'norm'),ops(C,'norm'),norm(b)]) > 1e6)
         kap0 = 10*blktrace(blk,X,Z);
      else
         kap0 = blktrace(blk,X,Z);
      end
      tau0 = 1; theta0 = 1;
   end
   kap = kap0; tau = tau0; theta = theta0;
%%
   normX0 = ops(X0,'norm')/tau; normZ0 = ops(Z0,'norm')/tau;
   bbar  = (tau*b-AXfun(blk,At,[],X))/theta;
   ZpATy = ops(Z,'+',Atyfun(blk,At,[],[],y));
   Cbar = ops(ops(ops(tau,'*',C),'-',ZpATy),'/',theta);
   gbar = (blktrace(blk,C,X)-b'*y+kap)/theta;
   abar = (blktrace(blk,X,Z)+tau*kap)/theta;
   for p = 1:size(blk,1);
      pblk = blk(p,:);
      if strcmp(pblk{1},'s')
         At{p} = [At{p}, -svec(pblk,Cnew{p},1), svec(pblk,Cbar{p},1)];
      else
         At{p} = [At{p}, -Cnew{p}, Cbar{p}];
      end
   end
   Bmat = [sparse(m,m), -b, bbar; b', 0, gbar; -bbar', -gbar, 0];
   em1 = zeros(m+2,1); em1(m+1) = 1;
   em2 = zeros(m+2,1); em2(m+2) = 1;
   par.Umat = [[b;0;0], [bbar;gbar;0], em1, em2];
   par.m = m;
   par.diagAAt = [full(diag(AAt)); 1; 1];
%%
%%-----------------------------------------
%% find the combined list of non-zero
%% elements of Aj, j = 1:k, for each k.
%%-----------------------------------------
%%
   par.numcolAt = length(b)+2;
   [At,C,Cnew,X,Z,par.permA,par.invpermA,par.permZ] = ...
    HSDsortA(blk,At,C,Cnew,[b;0;0],X,Z);
   [par.isspA,par.nzlistA,par.nzlistAsum,par.isspAy,par.nzlistAy] = ...
    nzlist(blk,At,par);
%%
%%-----------------------------------------
%% initialization
%%-----------------------------------------
%%
   y2 = [y; tau; theta];
   AX = AXfun(blk,At,par.permA,X);
   rp = [zeros(m,1); kap; -abar] - AX - Bmat*y2;
   Rd = ops(Atyfun(blk,At,par.permA,par.isspAy,-y2),'-',Z);
   trXZ = blktrace(blk,X,Z);
   mu  = (trXZ+kap*tau)/(nn+1);
   obj = [blktrace(blk,C,X), b'*y]/tau;
   gap = trXZ/tau^2;
   relgap = gap/(1+mean(abs(obj)));
   ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,[y;0;0]));
   ZpATynorm = ops(ZpATy,'norm');
   prim_infeas = norm(b - AX(1:m)/tau)/normb;
   dual_infeas = ops(ops(C,'-',ops(ZpATy,'/',tau)),'norm')/normC;
   infeas = max(prim_infeas,dual_infeas);
   pstep = 0; dstep = 0; pred_convg_rate = 1; corr_convg_rate = 1;
   prim_infeas_bad = 0;
   termcode  = -6;
   msg = [];
   runhist.pobj    = obj(1);
   runhist.dobj    = obj(2);
   runhist.gap     = gap;
   runhist.relgap  = relgap;
   runhist.pinfeas = prim_infeas;
   runhist.dinfeas = dual_infeas;
   runhist.infeas  = infeas;
   runhist.cputime = cputime-tstart;
   runhist.step    = 0;
   ttime.preproc   = runhist.cputime;
   ttime.pred = 0; ttime.pred_pstep = 0; ttime.pred_dstep = 0;
   ttime.corr = 0; ttime.corr_pstep = 0; ttime.corr_dstep = 0;
   ttime.pchol = 0; ttime.dchol = 0; ttime.misc = 0;
%%
%%-----------------------------------------
%% display parameters, and initial info
%%-----------------------------------------
%%
   if (printlevel >= 2)
      fprintf('\n********************************************');
      fprintf('***********************\n');
      fprintf('   SDPT3: homogeneous self-dual path-following algorithms');
      fprintf('\n********************************************');
      fprintf('***********************\n');
      [hh,mm,ss] = mytime(ttime.preproc);
      if (printlevel>=3)
         fprintf(' version  predcorr  gam  expon\n');
         if (vers == 1); fprintf('   HKM '); elseif (vers == 2); fprintf('    NT '); end
         fprintf('     %1.0f      %4.3f   %1.0f\n',predcorr,gam,expon);
         fprintf('it  pstep dstep p_infeas d_infeas  gap')
         fprintf('     mean(obj)    cputime\n');
         fprintf('------------------------------------------------');
         fprintf('-------------------\n');
         fprintf('%2.0f  %4.3f %4.3f %2.1e %2.1e',0,0,0,prim_infeas,dual_infeas);
         fprintf('  %2.1e %- 7.6e  %s:%s:%s',gap,mean(obj),hh,mm,ss);
         fprintf(' %2.1e %2.1e %2.1e',kap,tau,theta);
      end
   end
%%
%%---------------------------------------------------------------
%% start main loop
%%---------------------------------------------------------------
%%
   EE = ops(blk,'identity');
   normE = ops(EE,'norm'); Zpertold = 1;
   [Xchol,indef(1)] = blkcholfun(blk,X);
   [Zchol,indef(2)] = blkcholfun(blk,Z);
   if any(indef)
      msg = 'Stop: X, Z are not both positive definite';
      if (printlevel); fprintf('\n  %s\n',msg); end
      info.termcode = -3;
      info.msg1 = msg;
      return;
   end
%%
   breakyes = 0; dy = zeros(length(b),1); dtau = 0; dtheta = 0;
   for iter = 1:maxit;

       update_iter = 0; pred_slow = 0; corr_slow = 0; step_short = 0;
       tstart = cputime;
       time = zeros(1,11);
       time(1) = cputime;
       par.tau = tau;
       par.kap = kap;
       par.theta = theta;
       par.mu   = mu;
       par.iter = iter;
       par.y    = y;
       par.dy2  = [dy; dtau; dtheta];
       par.rp   = rp;
       par.ZpATynorm = ZpATynorm;
%%
%%--------------------------------------------------
%% perturb C associated with ublk
%%--------------------------------------------------
%%
       if (perturb_C)
          Cpertold = Cpert;
          for p = 1:size(blk,1)
             pblk = blk(p,:);
             n = sum(pblk{2});
             tmp = max(1,norm(C{p},'fro'))/sqrt(n);
             if (max(relgap,infeas) < 1e-6)
                if (norm(X{p},'fro') < 1e2); const=0.2; else; const=0.3; end
                Cpert(p) = max(const*Cpert(p),1e-10*tmp);
             elseif (max(relgap,infeas) < 1e-2)
                if (norm(X{p},'fro') < 1e2); const=0.4; else; const=0.5; end
                Cpert(p) = max(const*Cpert(p),1e-8*tmp);
   	     else
	        Cpert(p) = max(0.9*Cpert(p),1e-6*tmp);
             end
             Cpert = min(Cpert,Cpertold);
             if (prim_infeas < min([0.1*dual_infeas, 1e-7*runhist.pinfeas(1)])) ...
      	        & (iter > 1 & dual_infeas > 0.8*runhist.dinfeas(iter-1) & relgap < 1e-4)
                Cpert(p) = 0.5*Cpert(p);
             end
             if (dual_infeas < min([0.1*prim_infeas, 1e-7*runhist.dinfeas(1)])) ...
   	        & (iter > 1 & prim_infeas > 0.8*runhist.pinfeas(iter-1) & relgap < 1e-4)
                Cpert(p) = 0.5*Cpert(p);
             end
	     if (max(relgap,1e-2*infeas) < 1e-6 & relgap < 0.1*infeas)
                Cpert(p) = 0.5*Cpert(p);
             end
             if (prim_infeas < min([1e-4*dual_infeas,1e-7]) & theta < 1e-6) ...
	        | (prim_infeas < 1e-4 & theta < 1e-10)
                Cpert(p) = 0.1*Cpert(p); 0;
             end
             if (dual_infeas < min([1e-4*prim_infeas,1e-7]) & theta < 1e-6) ...
	        | (dual_infeas < 1e-4 & theta < 1e-10)
                Cpert(p) = 0.1*Cpert(p); 0;
             end
             if strcmp(pblk{1},'s')
                Cnew{p} = C{p} + Cpert(p)*speye(n);
                At{p}(:,par.invpermA(p,end-1)) = -svec(pblk,Cnew{p},1);
             else
                Cnew{p} = C{p} + Cpert(p)*ones(n,1);
                At{p}(:,par.invpermA(p,end-1)) = -Cnew{p};
             end
          end
          maxCpert(iter) = max(Cpert);
          fprintf(' %2.1e',max(Cpert));
          if (iter > 10 & norm(diff(maxCpert([iter-3,iter]))) < 1e-13)
             Cpert = 0.5*Cpert;
             maxCpert(iter) = max(Cpert);
          end
          AX = AXfun(blk,At,par.permA,X);
          rp = [zeros(m,1); kap; -abar] - AX - Bmat*y2;
          Rd = ops(Atyfun(blk,At,par.permA,par.isspAy,-y2),'-',Z);
       end
%%---------------------------------------------------------------
%% predictor step.
%%---------------------------------------------------------------
%%
       if (predcorr)
          sigma = 0;
       else
          sigma = 1-0.9*min(pstep,dstep);
          if (iter == 1); sigma = 0.5; end;
       end
       sigmu = sigma*mu;

       invXchol = cell(size(blk,1),1);
       invZchol = ops(Zchol,'inv');
       if (vers == 1);
          [par,dX,dy,dZ,coeff,L,hRd] = ...
           HSDHKMpred(blk,At,par,rp,Rd,sigmu,X,Z,invZchol);
       elseif (vers == 2);
          [par,dX,dy,dZ,coeff,L,hRd] = ...
           HSDNTpred(blk,At,par,rp,Rd,sigmu,X,Z,Zchol,invZchol);
       end
       if (solve_ok <= 0)
          msg = 'Stop: difficulty in computing predictor directions';
          if (printlevel); fprintf('\n  %s',msg); end
          runhist.cputime(iter+1) = cputime-tstart;
          termcode = -4;
          break;
       end
       time(2) = cputime;
       ttime.pred = ttime.pred + time(2)-time(1);
%%
%%-----------------------------------------
%% step-lengths for predictor step
%%-----------------------------------------
%%
      if (gam == 0)
         gamused = 0.9 + 0.09*min(pstep,dstep);
      else
         gamused = gam;
      end
      kapstep = max( (par.dkap<0)*(-kap/(par.dkap-eps)), (par.dkap>=0)*1e6 );
      taustep = max( (par.dtau<0)*(-tau/(par.dtau-eps)), (par.dtau>=0)*1e6 );
      [Xstep,invXchol] = steplength(blk,X,dX,Xchol,invXchol);
      time(3) = cputime;
      Zstep = steplength(blk,Z,dZ,Zchol,invZchol);
      time(4) = cputime;
      pstep = min(1,gamused*min([Xstep,Zstep,kapstep,taustep]));
      dstep = pstep;
      kappred = kap + pstep*par.dkap;
      taupred = tau + pstep*par.dtau;
      trXZpred = trXZ + pstep*blktrace(blk,dX,Z) + dstep*blktrace(blk,X,dZ) ...
                 + pstep*dstep*blktrace(blk,dX,dZ);
      mupred  = (trXZpred + kappred*taupred)/(nn+1);
      mupredhist(iter) = mupred;
      ttime.pred_pstep = ttime.pred_pstep + time(3)-time(2);
      ttime.pred_dstep = ttime.pred_dstep + time(4)-time(3);
%%
%%-----------------------------------------
%%  stopping criteria for predictor step.
%%-----------------------------------------
%%
      if (min(pstep,dstep) < steptol) & (stoplevel)
         msg = 'Stop: steps in predictor too short';
         if (printlevel)
            fprintf('\n  %s',msg);
            fprintf(': pstep = %3.2e,  dstep = %3.2e',pstep,dstep);
         end
         runhist.cputime(iter+1) = cputime-tstart;
         termcode = -2;
         breakyes = 1;
      end
      if (iter >= 2)
         idx = [max(2,iter-2) : iter];
         pred_slow = all(mupredhist(idx)./mupredhist(idx-1) > 0.4);
         idx = [max(2,iter-5) : iter];
         pred_convg_rate = mean(mupredhist(idx)./mupredhist(idx-1));
         pred_slow = pred_slow + (mupred/mu > 5*pred_convg_rate);
      end
      if (~predcorr)
         if (max(mu,infeas) < 1e-6) & (pred_slow) & (stoplevel)
            msg = 'Stop: lack of progress in predictor';
            if (printlevel)
               fprintf('\n  %s',msg);
               fprintf(': mupred/mu = %3.2f, pred_convg_rate = %3.2f.',...
               mupred/mu,pred_convg_rate);
            end
            runhist.cputime(iter+1) = cputime-tstart;
            termcode = -1;
            breakyes = 1;
         else
            update_iter = 1;
         end
      end
%%
%%---------------------------------------------------------------
%% corrector step.
%%---------------------------------------------------------------
%%
      if (predcorr) & (~breakyes)
         step_pred = min(pstep,dstep);
         if (mu > 1e-6)
            if (step_pred < 1/sqrt(3));
               expon_used = 1;
            else
               expon_used = max(expon,3*step_pred^2);
            end
         else
            expon_used = max(1,min(expon,3*step_pred^2));
         end
         sigma = min( 1, (mupred/mu)^expon_used );
         sigmu = sigma*mu;
%%
         if (vers == 1)
            [par,dX,dy,dZ] = HSDHKMcorr(blk,At,par,rp,Rd,sigmu,hRd,...
             dX,dZ,coeff,L,X,Z);
         elseif (vers == 2)
            [par,dX,dy,dZ] = HSDNTcorr(blk,At,par,rp,Rd,sigmu,hRd,...
             dX,dZ,coeff,L,X,Z);
         end
         if (solve_ok <= 0)
            msg = 'Stop: difficulty in computing corrector directions';
            if (printlevel); fprintf('\n  %s',msg); end
            runhist.cputime(iter+1) = cputime-tstart;
            termcode = -4;
            break;
         end
         time(5) = cputime;
         ttime.corr = ttime.corr + time(5)-time(4);
%%
%%-----------------------------------
%% step-lengths for corrector step
%%-----------------------------------
%%
         if (gam == 0)
            gamused = 0.9 + 0.09*min(pstep,dstep);
         else
            gamused = gam;
         end
         kapstep = max( (par.dkap<0)*(-kap/(par.dkap-eps)), (par.dkap>=0)*1e6 );
         taustep = max( (par.dtau<0)*(-tau/(par.dtau-eps)), (par.dtau>=0)*1e6 );
         Xstep = steplength(blk,X,dX,Xchol,invXchol);
         time(6) = cputime;
         Zstep = steplength(blk,Z,dZ,Zchol,invZchol);
         time(7) = cputime;
         pstep = min(1,gamused*min([Xstep,Zstep,kapstep,taustep]));
         dstep = pstep;
         kapcorr = kap + pstep*par.dkap;
         taucorr = tau + pstep*par.dtau;
         trXZcorr = trXZ + pstep*blktrace(blk,dX,Z) + dstep*blktrace(blk,X,dZ)...
                    + pstep*dstep*blktrace(blk,dX,dZ);
         mucorr = (trXZcorr+kapcorr*taucorr)/(nn+1);
         ttime.corr_pstep = ttime.corr_pstep + time(6)-time(5);
         ttime.corr_dstep = ttime.corr_dstep + time(7)-time(6);
%%
%%-----------------------------------------
%%  stopping criteria for corrector step
%%-----------------------------------------
%%
         if (iter >= 2)
            idx = [max(2,iter-2) : iter];
            corr_slow = all(runhist.gap(idx)./runhist.gap(idx-1) > 0.8);
            idx = [max(2,iter-5) : iter];
            corr_convg_rate = mean(runhist.gap(idx)./runhist.gap(idx-1));
            corr_slow = corr_slow + (mucorr/mu > max(min(1,5*corr_convg_rate),0.8));
         end
	 if (max(mu,infeas) < 1e-6) & (iter > 10) & (stoplevel) ...
            & (corr_slow & mucorr/mu > 1.0)
            msg = 'Stop: lack of progress in corrector';
   	    if (printlevel)
               fprintf('\n  %s',msg);
               fprintf(': mucorr/mu = %3.2f, corr_convg_rate = %3.2f',...
               mucorr/mu,corr_convg_rate);
            end
            runhist.cputime(iter+1) = cputime-tstart;
            termcode = -1;
            breakyes = 1;
         else
            update_iter = 1;
         end
      end
%%
%%---------------------------------------------------------------
%% udpate iterate
%%---------------------------------------------------------------
%%
      indef = [1 1];
      if (update_iter)
         for t = 1:5
            [Xchol,indef(1)] = blkcholfun(blk,ops(X,'+',dX,pstep)); time(8) = cputime;
            if (predcorr); ttime.pchol = ttime.pchol + time(8)-time(7);
            else;          ttime.pchol = ttime.pchol + time(8)-time(4);
            end
            if (indef(1)); pstep = 0.8*pstep; else; break; end
         end
	 if (t > 1); pstep = gamused*pstep; end
	 for t = 1:5
            [Zchol,indef(2)] = blkcholfun(blk,ops(Z,'+',dZ,dstep)); time(9) = cputime;
            if (predcorr); ttime.dchol = ttime.dchol + time(9)-time(8);
            else;          ttime.dchol = ttime.dchol + time(9)-time(4);
            end
            if (indef(2)); dstep = 0.8*dstep; else; break; end
         end
	 if (t > 1); dstep = gamused*dstep; end
         AdX = AXfun(blk,At,par.permA,dX);
         AXtmp = AX(1:m) + pstep*AdX(1:m); tautmp = par.tau+pstep*par.dtau;
         prim_infeasnew = norm(b-AXtmp/tautmp)/normb;
         if any(indef)
            msg = 'Stop: X, Z not both positive definite';
            if (printlevel); fprintf('\n  %s',msg); end
   	    termcode = -3;
            breakyes = 1;
         elseif (prim_infeasnew > max([1e-8,relgap,10*infeas])) ...
	    | (prim_infeasnew > max([1e-4,20*prim_infeas]) & (infeas < 1e-2)) ...
            | (prim_infeasnew > max([1e-6,3*prim_infeas,10*dual_infeas]) ...
            & max([relgap,dual_infeas]) < 1e-5)
            if (stoplevel) & (max(pstep,dstep)<=1) & (kap < 1e-3)
               msg = 'Stop: primal infeas has deteriorated too much';
               if (printlevel); fprintf('\n  %s, %2.1e',msg,prim_infeasnew); end
               termcode = -7;
               breakyes = 1;
            end
	 end
         if (~breakyes)
            X = ops(X,'+',dX,pstep);
            y = y + dstep*dy;  Z = ops(Z,'+',dZ,dstep);
            theta = max(0, theta + pstep*par.dtheta);
            kap = kap + pstep*par.dkap;
            if (tau + pstep*par.dtau > theta)
               tau = tau + pstep*par.dtau;
            end
         end
      end
%%
%%--------------------------------------------------
%% perturb Z: do this step before checking for break
%%--------------------------------------------------
      perturb_Z = 1;
      if (~breakyes) & (perturb_Z)
         trXZtmp = blktrace(blk,X,Z);
         trXE  = blktrace(blk,X,EE);
         Zpert = max(1e-12,0.2*min(relgap,prim_infeas)).*normC./normE;
         Zpert = min(Zpert,0.1*trXZtmp./trXE);
         Zpert = min([1,Zpert,1.5*Zpertold]);
         if (infeas < 1e-2)
            Z = ops(Z,'+',EE,Zpert);
            [Zchol,indef(2)] = blkcholfun(blk,Z);
            if any(indef(2))
               msg = 'HSDsqlp stop: Z not positive definite';
               if (printlevel); fprintf('\n  %s',msg); end
               termcode = -3;
               breakyes = 1;
            end
         end
         Zpertold = Zpert;
      end
%%
%%---------------------------------------------------------------
%% compute rp, Rd, infeasibities, etc.
%%---------------------------------------------------------------
%%
      y2 = [y; tau; theta];
      AX = AXfun(blk,At,par.permA,X);
      rp = [zeros(m,1); kap; -abar] - AX - Bmat*y2;
      Rd = ops(Atyfun(blk,At,par.permA,par.isspAy,-y2),'-',Z);
      trXZ = blktrace(blk,X,Z);
      mu  = (trXZ+kap*tau)/(nn+1);
      obj = [blktrace(blk,C,X), b'*y]/tau;
      gap = trXZ/tau^2;
      relgap = gap/(1+mean(abs(obj)));
      ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,[y;0;0]));
      prim_infeas = norm(b-AX(1:m)/tau)/normb;
      dual_infeas = ops(ops(C,'-',ops(ZpATy,'/',tau)),'norm')/normC;
      infeas = max(prim_infeas,dual_infeas);
      runhist.pobj(iter+1)    = obj(1);
      runhist.dobj(iter+1)    = obj(2);
      runhist.gap(iter+1)     = gap;
      runhist.relgap(iter+1)  = relgap;
      runhist.pinfeas(iter+1) = prim_infeas;
      runhist.dinfeas(iter+1) = dual_infeas;
      runhist.infeas(iter+1)  = infeas;
      runhist.cputime(iter+1) = cputime-tstart;
      runhist.step(iter+1)    = min(pstep,dstep);
      time(10) = cputime;
      ttime.misc = ttime.misc + time(10)-time(9);
      [hh,mm,ss] = mytime(sum(runhist.cputime));
      if (printlevel>=3)
         fprintf('\n%2.0f  %4.3f %4.3f',iter,pstep,dstep);
         fprintf(' %2.1e %2.1e  %2.1e',prim_infeas,dual_infeas,gap);
         fprintf(' %- 7.6e  %s:%s:%s',mean(obj),hh,mm,ss);
         fprintf(' %2.1e %2.1e %2.1e',kap,tau,theta);
      end
%%
%%--------------------------------------------------
%% check convergence.
%%--------------------------------------------------
%%
      ZpATynorm = ops(ZpATy,'norm');
      if (obj(2) > 0); homRd = (ZpATynorm/tau)/obj(2);   else; homRd = inf; end
      if (obj(1) < 0); homrp = (norm(AX(1:m))/tau)/(-obj(1)); else; homrp = inf; end
      if (ops(X,'norm')/tau > 1e15*normX0 | ops(Z,'norm')/tau > 1e15*normZ0)
         termcode = 3;
         breakyes = 1;
      end
      if (homRd < min(1e-6,1e-2*sqrt(max([infeas,relgap]*inftol)))) ...
         | (homRd < 10*tau & tau < 1e-7)
         termcode = 1;
         breakyes = 1;
      end
      if (homrp < min(1e-6,1e-2*sqrt(max([infeas,relgap]*inftol)))) ...
         | (homrp < 10*tau & tau < 1e-7)
         termcode = 2;
         breakyes = 1;
      end
      if (max(relgap,infeas) < gaptol)
         msg = sprintf('Stop: max(relative gap, infeasibilities) < %3.2e',gaptol);
         if (printlevel); fprintf('\n  %s',msg); end
         termcode = 0;
         breakyes = 1;
      end
      if (stoplevel)
         min_prim_infeas = min(runhist.pinfeas(1:iter));
         prim_infeas_bad = prim_infeas_bad + (prim_infeas > ...
            max(1e-10,5*min_prim_infeas) & (min_prim_infeas < 1e-2));
         if (mu < 1e-6)
            idx = [max(1,iter-1): iter];
         elseif (mu < 1e-3);
            idx = [max(1,iter-2): iter];
         else
            idx = [max(1,iter-3): iter];
         end
         idx2 = [max(1,iter-4): iter];
         gap_ratio2 = runhist.gap(idx2+1)./runhist.gap(idx2);
         gap_slowrate = min(0.8,max(0.6,2*mean(gap_ratio2)));
         gap_ratio = runhist.gap(idx+1)./runhist.gap(idx);
         if (infeas < 1e-4 | prim_infeas_bad) & (relgap < 1e-3) ...
            & (iter > 5) & (prim_infeas > (1-pstep/2)*runhist.pinfeas(iter))
            gap_slow = all(gap_ratio > gap_slowrate) & (relgap < 1e-3);
            min_pinfeas = min(runhist.pinfeas);
            const = 0.1;
            if (relgap < const*max(prim_infeas,dual_infeas)) ...
               & ((runhist.step(iter+1) < 0.5) | ...
                  (prim_infeas > 10*min_pinfeas & min_pinfeas < 1e-6)) ...
               & (dual_infeas > 0.8*runhist.dinfeas(iter) | (dual_infeas < 1e-2*gaptol))
               msg = 'Stop: relative gap < infeasibility';
               if (printlevel); fprintf('\n  %s',msg); end
               termcode = -1;
               breakyes = 1;
            elseif (gap_slow) & (infeas > 0.8*runhist.infeas(iter)) ...
               & (theta < 1e-8)
               msg = 'Stop: progress is too slow';
               if (printlevel); fprintf('\n  %s',msg); end
               termcode = -5;
               breakyes = 1;
            end
         elseif (prim_infeas_bad) & (iter >50) & all(gap_ratio > gap_slowrate)
            msg = 'Stop: progress is bad';
            if (printlevel); fprintf('\n  %s',msg); end
            termcode = -5;
            breakyes = 1;
         elseif (infeas < 1e-8) & (gap > 1.2*mean(runhist.gap(idx)))
            msg = 'Stop: progress is bad*';
            if (printlevel); fprintf('\n  %s',msg); end
            termcode = -5;
            breakyes = 1;
         end
         if (max([relgap,infeas]) < 1e-3) & (iter > 10) ...
            & (runhist.pinfeas(iter+1) > 0.9*runhist.pinfeas(max(1,iter-5))) ...
            & (runhist.dinfeas(iter+1) > 0.9*runhist.dinfeas(max(1,iter-5)))
            msg = 'Stop: progress is bad**';
            if (printlevel); fprintf('\n  %s',msg); end
            termcode = -5;
            breakyes = 1;
         end
         if (min(runhist.infeas) < 1e-4 | prim_infeas_bad) ...
            & (max(runhist.infeas) > 1e-4) & (iter > 5)
            relgap2 = abs(diff(obj))/(1+mean(abs(obj)));
            if (relgap2 < 1e-3);
               step_short = all(runhist.step([iter:iter+1]) < 0.1) ;
            elseif (relgap2 < 1)
               idx = [max(1,iter-3): iter+1];
               step_short = all(runhist.step(idx) < 0.05);
            else
               step_short = 0;
            end
            if (step_short)
               msg = 'Stop: steps too short consecutively';
               if (printlevel); fprintf('\n  %s',msg); end
               termcode = -5;
               breakyes = 1;
            end
         end
	 if (infeas > 100*max(1e-12,min(runhist.infeas)) & relgap < 1e-4)
            msg = 'Stop: infeas has deteriorated too much';
            if (printlevel); fprintf('\n  %s, %3.1e',msg,infeas); end
            X = ops(X,'-',dX,pstep);
            y = y - dstep*dy;  Z = ops(Z,'-',dZ,dstep);
            kap = kap - pstep*par.dkap; tau = tau - pstep*par.dtau;
            theta = theta - pstep*par.dtheta;
            prim_infeas = runhist.pinfeas(iter); dual_infeas = runhist.dinfeas(iter);
            gap = runhist.gap(iter); relgap = runhist.relgap(iter);
            obj = [runhist.pobj(iter), runhist.dobj(iter)];
            termcode = -7;
            breakyes = 1;
         end
         if (iter > 3 & iter < 20) & (max(runhist.step(max(1,iter-3):iter+1)) < 1e-3) ...
            & (infeas > 1) & (min(homrp,homRd) > 1000*inftol)
            if (stoplevel >= 2)
               msg = 'Stop: steps too short consecutively';
               if (printlevel); fprintf('\n  %s',msg); end
               termcode = -5;
               breakyes = 1;
            end
         end
      end
      if (breakyes > 0.5); break; end
   end
%%---------------------------------------------------------------
%% end of main loop
%%---------------------------------------------------------------
%%
   if (termcode == -6)
      msg = 'Stop: maximum number of iterations reached';
      if (printlevel); fprintf('\n  %s',msg); end
   end
%%
%%---------------------------------------------------------------
%% produce infeasibility certificates if appropriate
%%---------------------------------------------------------------
%%
   X = ops(X,'/',tau); y = y/tau; Z = ops(Z,'/',tau);
   if (iter >= 1)
      param.obj         = obj;
      param.relgap      = relgap;
      param.prim_infeas = prim_infeas;
      param.dual_infeas = dual_infeas;
      param.ZpATynorm   = ZpATynorm/tau;
      param.inftol      = inftol;
      param.m0          = m0;
      param.indeprows   = indeprows;
      param.termcode    = termcode;
      param.AX          = AX(1:m)/tau;
      param.normX0      = normX0;
      param.normZ0      = normZ0;
      param.printlevel  = printlevel;
      [X,y,Z,resid,reldist,param,msg2] = ...
      HSDsqlpmisc(blk,At,C,b,X,y,Z,par.permZ,param);
      termcode = param.termcode;
   end
%%
%%---------------------------------------------------------------
%% recover unrestricted blk from linear blk
%%---------------------------------------------------------------
%%
   for p = 1:size(blk,1)
      if (ublkidx(p) == 1)
         n = blk{p,2}/2;
         X{p} = X{p}(1:n)-X{p}(n+[1:n]);
         Z{p} = Z{p}(1:n);
      end
   end
%%
%%---------------------------------------------------------------
%%  print summary
%%---------------------------------------------------------------
%%
   dimacs = [prim_infeas; 0; dual_infeas; 0];
   dimacs = [dimacs; [-diff(obj); gap]/(1+sum(abs(obj)))];
   info.dimacs   = dimacs;
   info.termcode = termcode;
   info.iter     = iter;
   info.obj      = obj;
   info.gap      = gap;
   info.relgap   = relgap;
   info.pinfeas  = prim_infeas;
   info.dinfeas  = dual_infeas;
   info.cputime  = sum(runhist.cputime);
   info.resid    = resid;
   info.reldist  = reldist;
   info.normX    = ops(X,'norm');
   info.normy    = norm(y);
   info.normZ    = ops(Z,'norm');
   info.normA    = ops(At,'norm');
   info.normb    = norm(b);
   info.normC    = ops(C,'norm');
   info.msg1     = msg;
   info.msg2     = msg2;
%%
   sqlpsummary(info,ttime,[],printlevel);
   rand('state',randstate);
%%*****************************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function cursor_demo()
% CURSOR_DEMO - Demo of the MYCURSORS function

% Pictures for the toolbar buttons
pics = CreatePics();

% Create the figure and controls
handles.Fig    = figure('Tag','Fig','Name','Cursors Demo','NumberTitle','off','MenuBar','none','Toolbar','none','Position',[200 200 550 550],'Color',[1 1 1]);
handles.Axe    = axes('Parent',handles.Fig,'Tag','Axe','Units','pixels','Position',[50 200 450 325]);
handles.TbCurs = uitoolbar('Parent',handles.Fig,'Tag','TbCurs');

handles.PtAdd1  = uipushtool('Parent',handles.TbCurs,'Tag','PtAdd1','CData',pics.Plus,'ClickedCallback',@PtAdd1_ClickedCallback,'TooltipString','Add cursor');
handles.PtDel1  = uipushtool('Parent',handles.TbCurs,'Tag','PtDel1','CData',pics.Minus,'ClickedCallback',@PtDel1_ClickedCallback,'TooltipString','Remove last cursor');
handles.PtVal1  = uipushtool('Parent',handles.TbCurs,'Tag','PtVal1','CData',pics.Val,'ClickedCallback',@PtVal1_ClickedCallback,'TooltipString','Display cursors positions');

handles.PtAdd2  = uipushtool('Parent',handles.TbCurs,'Tag','PtAdd2','CData',pics.PlusR,'ClickedCallback',@PtAdd2_ClickedCallback,'TooltipString','Add cursor');
handles.PtDel2  = uipushtool('Parent',handles.TbCurs,'Tag','PtDel2','CData',pics.MinusR,'ClickedCallback',@PtDel2_ClickedCallback,'TooltipString','Remove last cursor');
handles.PtVal2  = uipushtool('Parent',handles.TbCurs,'Tag','PtVal2','CData',pics.ValR,'ClickedCallback',@PtVal2_ClickedCallback,'TooltipString','Display cursors positions');

handles.StVal1  = uicontrol('Parent',handles.Fig,'Tag','StVal1','style','text','units','pixels','position',[ 10 10 210 150],'string','','HorizontalAlignment','left','FontName','Courier New','BackgroundColor',[1 1 1],'ForegroundColor','k');
handles.StVal2  = uicontrol('Parent',handles.Fig,'Tag','StVal2','style','text','units','pixels','position',[230 10 210 150],'string','','HorizontalAlignment','left','FontName','Courier New','BackgroundColor',[1 1 1],'ForegroundColor','r');

% Plot data
plot(sin(-20:.01:20),'Parent',handles.Axe);

% Init Cursor structure
mycurs1 = [];
mycurs2 = [];

%*******************************************************************************

% --- FIRST CURSOR SET ----------------------------------------------

  function PtAdd1_ClickedCallback(varargin)
    % Adds a cursor on the graph

    % If no cursor was created, simply initializes one
    if isempty(mycurs1)
      mycurs1 = mycursors(handles.Axe,'k'); % Create the first cursor and store the structure of function handles
    else
      mycurs1.add(); % Add a cursor to the graph
    end

  end

  function PtDel1_ClickedCallback(varargin)
    % Removes the last cursor added

    % Execute code only if at least one cursor has already been added
    if ~isempty(mycurs1)
      mycurs1.off('last');
    end

  end

  function PtVal1_ClickedCallback(varargin)
    % Displays the positions of all cursors

    % Execute code only if at least one cursor has already been added
    if ~isempty(mycurs1)

      % Get cursors positions
      data = mycurs1.val();

      % Create a matrix with cursor indices on the 1st line and cursor positions
      % on the 2nd line
      data = [1:length(data);data(:)'];

      % Display info in the static text below the axe
      set(handles.StVal1,'string',sprintf('  Cursor %u: %10.2f\n',data(:)))

    end

  end

% --- SECOND CURSOR SET ---------------------------------------------

  function PtAdd2_ClickedCallback(varargin)
    % Adds a cursor on the graph

    % If no cursor was created, simply initializes one
    if isempty(mycurs2)
      mycurs2 = mycursors(handles.Axe,'r'); % Create the first cursor and store the structure of function handles
    else
      mycurs2.add(); % Add a cursor to the graph
    end

  end

  function PtDel2_ClickedCallback(varargin)
    % Removes the last cursor added

    % Execute code only if at least one cursor has already been added
    if ~isempty(mycurs2)
      mycurs2.off('last');
    end

  end

  function PtVal2_ClickedCallback(varargin)
    % Displays the positions of all cursors

    % Execute code only if at least one cursor has already been added
    if ~isempty(mycurs2)

      % Get cursors positions
      data = mycurs2.val();

      % Create a matrix with cursor indices on the 1st line and cursor positions
      % on the 2nd line
      data = [1:length(data);data(:)'];

      % Display info in the static text below the axe
      set(handles.StVal2,'string',sprintf('  Cursor %u: %10.2f\n',data(:)))

    end

  end

% --- OTHER NESTED FUNCTIONS ----------------------------------------

  function pics = CreatePics()
    % Generate pictures matrices to use in the toolbar

    % + sign
    pics.Plus  = repmat([ ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN   0   0   0   0   0   0   0   0   0   0 NaN NaN NaN ; ...
      NaN NaN NaN   0   0   0   0   0   0   0   0   0   0 NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN   ...
      ],[1,1,3]);

    % - sign
    pics.Minus = repmat([ ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN   0   0   0   0   0   0   0   0   0   0 NaN NaN NaN ; ...
      NaN NaN NaN   0   0   0   0   0   0   0   0   0   0 NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN   ...
      ],[1,1,3]);

    % Matrix
    pics.Val = repmat([ ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN   0   0 NaN NaN   0   0   0 NaN NaN   0   0   0 NaN NaN ; ...
      NaN   0 NaN   0 NaN NaN NaN NaN   0 NaN NaN NaN NaN   0 NaN NaN ; ...
      NaN NaN NaN   0 NaN NaN   0   0   0 NaN NaN NaN   0   0 NaN NaN ; ...
      NaN NaN NaN   0 NaN NaN   0 NaN NaN NaN NaN NaN NaN   0 NaN NaN ; ...
      NaN NaN NaN   0 NaN NaN   0   0   0 NaN NaN   0   0   0 NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN   0 NaN NaN NaN NaN   0   0   0 NaN NaN   0   0   0 NaN NaN ; ...
      NaN   0 NaN   0 NaN NaN   0 NaN NaN NaN NaN   0 NaN NaN NaN NaN ; ...
      NaN   0   0   0 NaN NaN   0   0   0 NaN NaN   0   0   0 NaN NaN ; ...
      NaN NaN NaN   0 NaN NaN NaN NaN   0 NaN NaN   0 NaN   0 NaN NaN ; ...
      NaN NaN NaN   0 NaN NaN   0   0   0 NaN NaN   0   0   0 NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ; ...
      NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN   ...
      ],[1,1,3]);

    pS = [16 16 3];
    col = [1 0 0]; % Red

    % + sign in red
    pics.PlusR = pics.Plus;
    pics.PlusR(pics.PlusR==0) = 1;
    pics.PlusR = reshape( repmat(col,pS(1)*pS(2),1) .* reshape(pics.PlusR,pS(1)*pS(2),3) , pS );

    % - sign in red
    pics.MinusR = pics.Minus;
    pics.MinusR(pics.MinusR==0) = 1;
    pics.MinusR = reshape( repmat(col,pS(1)*pS(2),1) .* reshape(pics.MinusR,pS(1)*pS(2),3) , pS );

    % Matrix in red
    pics.ValR = pics.Val;
    pics.ValR(pics.ValR==0) = 1;
    pics.ValR = reshape( repmat(col,pS(1)*pS(2),1) .* reshape(pics.ValR,pS(1)*pS(2),3) , pS );

  end

end

author: Sable

repo: mcbench-benchmarks


text:
function F = set(varargin)
%SET OBSOLETE

disp('SET has been considered obsolete for many years, and the time has come...');
disp('Update your code. https://yalmip.github.io/tutorial/basics/');
error('Obsolete command')

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef prtUiClassLabeler < hgsetget

    properties

        useClassifier = true;

        dataSet
        classifier = prtClassMap;
        handles

        onClickDisplayFunction = @(ds, obsInd, featureInds)plotAsTimeSeries(ds.retainObservations(obsInd));
    end
    methods
        function self = prtUiClassLabeler(varargin)
            self = prtUtilAssignStringValuePairs(self, varargin{:});

            if isempty(self.dataSet)
                error('dataSet must be defined on input');
            end

            if isempty(self.dataSet.Y)
                self.dataSet.Y = nan(self.dataSet.nObservations,1);
            end

            self.init();
        end
        function init(self)
            self.handles.mainFigure = figure;
            gca;
            if self.useClassifier
                self.handles.classifierFigure = figure;
            end

            self.udpate();
        end
        function udpate(self)
            axesObj = findobj(self.handles.mainFigure,'type','axes');
            %oldView = view(axesObj);
            oldXlim = xlim(axesObj);
            oldYlim = ylim(axesObj);
            oldZlim = zlim(axesObj);
            figure(self.handles.mainFigure)
            explore(self.dataSet, @(ds,obsInd,featureInd)self.onClick(ds,obsInd,featureInd));

            %view(axesObj,oldView);
            xlim(axesObj,oldXlim);
            ylim(axesObj,oldYlim);
            zlim(axesObj,oldZlim);

            if self.useClassifier && self.dataSet.nClasses > 1 && all(self.dataSet.nObservationsByClass > 2*self.dataSet.nFeatures)
                self.classifier = train(self.classifier, self.dataSet);
                figure(self.handles.classifierFigure)
                plot(self.classifier);
            end
        end

        function onClick(self,ds,obsInd, featureInd)
            self.handles.secondaryFigure = figure;
            self.handles.secondaryFigure.KeyPressFcn = @(h,e)self.onClickKeypress(h,e,obsInd);
            self.handles.secondaryFigure.WindowStyle = 'modal';

            self.onClickDisplayFunction(ds, obsInd, featureInd);

        end
        function onClickKeypress(self, h, e, obsInd)
            key = e.Key;
            switch key
                case {'0','1','2','3','4','5','6','7','8','9'}
                    newY = str2double(key);
                    self.dataSet.Y(obsInd) = newY;

                    close(self.handles.secondaryFigure);
                    self.udpate();
                otherwise
                    title('Invalid key')
            end

        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
%--- help for setup_restrictions ---
%
%  INTERNAL FUNCTION: sets up linear, nonlinear and general restrictions
%
%  ::
%
% 	obj=setup_general_restrictions(obj)
%
%  Args:
%
%     obj (rise | dsge | rfvar | svar): RISE object
%
%  Returns:
%     :
%
%     - **obj** [rise|dsge|rfvar|svar]: RISE object
%
%  See also:
%     - rise_generic.setup_linear_restrictions,
%     - rise_generic.setup_nonlinear_restrictions,
%     - rise_generic.setup_general_restrictions
%
%

author: jmaih


text:
function [v,psnrall] = gapwnnm_int( y, opt )
%GAPWNNM_INT Integrated GAP and WNNM for snapshot compressive imaging.
%   Note integrated version of GAP-WNNM instead of embedded for
%   acceleration of WNNM with less block-match process for multiple
%   iterations.
%   Reference(s)
%     [1] Y. Liu, X. Yuan, J. Suo, D.J. Brady, and Q. Dai, Rank Minimization
%         for Snapshot Compressive Imaging, preprint, 2018.
%   Code credit
%     Xin Yuan, Bell Labs, xyuan@bell-labs.com, initial version Jul 2,
%       2015.
%     Yang Liu, Tsinghua University, y-liu16@mails.tsinghua.edu.cn, last
%       update Jul 13, 2018.
%   See also GAPDENOISE.
if nargin<2
    opt = [];
end
% [0] default parameter configuration, to be specified
A  = @(x) M_func(x);
At = @(z) Mt_func(z);
if isfield(opt,'Mfunc'),  A  = @(x) opt.Mfunc(x);  end
if isfield(opt,'Mtfunc'), At = @(z) opt.Mtfunc(z); end
% Phisum   = ;
% v0       = At(y); % start point (initialization of iteration)
lambda   = 1;     % correction coefficiency
maxiter  = 300;     % maximum number of iteration
acc      = 1;       % enable acceleration
sigma    = 10/255;  % noise deviation
flag_iqa = true;    % flag of showing image quality assessments
save_iter_image = false; % flag of saving all the images of the iteration
                         %  process (true only necessary)
                         %  require assigned directory opt.iter_image_dir

if isfield(opt,'Phisum'),     Phisum = opt.Phisum;   end
if isfield(opt,'v0'),             v0 = opt.v0;       end
if isfield(opt,'lambda'),   lambda   = opt.lambda;   end
if isfield(opt,'maxiter'),   maxiter = opt.maxiter;  end
if isfield(opt,'acc'),           acc = opt.acc;      end
if isfield(opt,'sigma'),       sigma = opt.sigma;    end
if isfield(opt,'flag_iqa'), flag_iqa = opt.flag_iqa; end
if isfield(opt,'save_iter_image'), save_iter_image = opt.save_iter_image; end

if ~exist('v0','var') || isempty(v0)
    v0 = At(y); % start point (initialization of iteration)
end
y1 = zeros(size(y),'like',y);
psnrall = []; % return empty with no ground truth
% ssimall = []; % return empty with no ground truth
% [1] start iteration
v = v0; % initialization
% opt.sigma = sigma;
k = 1; % current number of iteration
nonlocalarr = uint32.empty;
vindarr     = uint8.empty;
for isig = 1:length(maxiter)
    opt.sigma = sigma(isig);
    for iter = 1:maxiter(isig)
        for ii = 1:1
            % [1.1] Euclidean projection
            yb = A(v);
            if acc % enable acceleration
                y1 = y1+(y-yb);
                v = v+lambda*(At((y1-yb)./Phisum)); % v=v+lambda*(At*A)^-1*At*dy
            else
                v = v+lambda*(At((y-yb)./Phisum));
            end
        end
        % [1.2] Denoising to match the video prior
        % WNNM video denoising (MATLAB-style matrix version)
        % v = wnnm_vdenoise(v,[],opt); % opt.sigma
        noisyv = v;
        para = vdefparaconf(opt); % default parameter configuration

        [nrow,ncol,nframe] = size(noisyv); % grayscale video (color todo)
        % [1] pre-calculation of the indexes of the neighbors within the search
        %     window
        [neighborindarr,neighbornumarr,selfindarr] = neighborind([nrow ncol],para);

        % [2] WNNM denoisng for several iterations
        estv = noisyv;
        for iit = 1:para.iternum
            % correction between adjacent iterations
            if ~para.adaptboost %
                estv = estv + para.delta*(noisyv-estv);
            else % adaptive boosting for WNNM-based denoising
                Eestv = mean(abs(estv(:)).^2);
                Enoise  = para.abeta*abs(para.nsigma^2-var(noisyv(:)-estv(:)));
                rho = sqrt(Eestv)/(sqrt(Eestv)+sqrt(max(Eestv-Enoise,0)));
                fprintf('    Iteration % 2d, rho = %.3f.\n',iit,rho);
                estv = estv + (1-rho)*(noisyv-estv);
            end
            % use all frames for denoising
            curvall = zeros(nrow,ncol,nframe,nframe);
            curfall = zeros(nrow,ncol,nframe,nframe);

            % [2.1] video to patches
            [rawpatchmat,nsigmamat] = v2patch(estv,noisyv,para);
            blockmatch_period = para.blockmatch_period;
            % inner loop to reuse the block-matching results
                if mod(k,blockmatch_period) == 1 % block matching periodically
                    nonlocalarr = uint32.empty;
                    vindarr     = uint8.empty;
                    parfor (iframe = 1:nframe,nframe) % maximum nframe parpool for optimal
                    % [2.2] calculate the patches with non-local similarity for each
                    %       key patch
                    [curnonlocalarr,curvindarr] = vblockmatch(iframe,rawpatchmat,...
                        neighborindarr,neighbornumarr,selfindarr,para);
                    nonlocalarr(:,:,iframe) = curnonlocalarr;
                    vindarr(:,:,iframe) = curvindarr;
                    end
                end
                if iit==1 % initial noise level of each patch
                    nsigmamat = para.nsigma*ones(size(nsigmamat));
                end
            % frame-wise denosing using parfor instead of for to get accelerated
            % performance, requiring Parallel Computing Toolbox.
            parfor (iframe = 1:nframe,nframe) % maximum nframe parpool for optimal
            % for iframe = 1:nframe % maximum nframe parpool for optimal
                % use all frames for denoising
                % inner loop to reuse the block-matching results
                % if mod(k,blockmatch_period) == 1 % block matching periodically
                %     % [2.2] calculate the patches with non-local similarity for each
                %     %       key patch
                %     [nonlocalarr(:,:,iframe),vindarr(:,:,iframe)] = vblockmatch(iframe,rawpatchmat,...
                %         neighborindarr,neighbornumarr,selfindarr,para);
                % end
                curnonlocalarr = nonlocalarr(:,:,iframe);
                curvindarr = vindarr(:,:,iframe);
                % [2.3] patch estimation by means of WNNM
                [estpatchmat,frqpatchmat] = vpatchestimate(iframe,curnonlocalarr,...
                    curvindarr,rawpatchmat,nsigmamat,selfindarr,para);
                % [2.4] aggregate overlapped patches to the whole image
                [curv,curf] = patch2v(estpatchmat,frqpatchmat,size(noisyv),...
                    para.patchsize);
                % use all frames for denoising
                curvall(:,:,:,iframe) = curv;
                curfall(:,:,:,iframe) = curf;
            end
            % use all frames for denoising
            v_ = sum(curvall,ndims(curvall));
            f_ = sum(curfall,ndims(curfall));

            estv = v_./(f_+eps);

            if mod(iit-1,para.innerloop)==0
                % % [2.2] calculate the patches with non-local similarity for
                % %       each key patch
                % [nonlocalarr,vindarr] = vblockmatch(cframe,rawpatchmat,...
                %     neighborindarr,neighbornumarr,selfindarr,para);
                % less non-local patches with lower noise level
                para.patchnum = para.patchnum-10;
            end
            % [1.4] save and show intermediate results of psnr and ssim
            if flag_iqa && isfield(opt,'orig') && (~isempty(opt.orig))
                psnrall(k) = psnr(double(estv),double(opt.orig)); % record all psnr
                % ssimall(k) = ssim(double(estv),opt.orig); % record all ssim
                if (mod(k,5)==0)
                    fprintf('  GAP-%s iteration % 4d, sigma %.1f, PSNR %2.2f dB.\n',...
                        upper(opt.denoiser),k,opt.sigma*255,psnrall(k));
                end
                if save_iter_image % save all the images of the iteration process
                    MAXB = opt.MAXB;
                    % save each frame in the assigned directory
                    nim = size(estv,ndims(estv)); % number of frames in the video
                    for iim = 1:nim
                        im = estv(:,:,iim);
                        impsnr = psnr(im,opt.orig(:,:,iim));
                        imssim = ssim(im,opt.orig(:,:,iim));
                        if k == 1 % mkdir for the first iteration
                            subdir = sprintf('%s/frame%02d',opt.iter_image_dir,iim); % frame-wise
                            if ~exist(subdir,'dir')
                                mkdir(subdir);
                            end
                        end
                        imwrite(uint8(im*MAXB),sprintf('%s/frame%02d/frame%02d_iter%03d_sigma%.1f_psnr%2.2f_ssim%.4f.png',...
                            opt.iter_image_dir,iim,iim,k,opt.sigma*MAXB,impsnr,imssim));
                    end
                end
            end
            k = k+1;
        end % WNNM loop [iternum]
        v = estv;
        % % [1.3] update noise standard deviation
        % nsigma = eta*sqrt(abs(sigma^2-var(v(:)-v0(:))));
        % opt.sigma = nsigma;

    end % GAP loop [maxiter]
end % sigma loop [length(sigma)]

end

repo: DeSCI


text:
function plugin_deactivate(foldername);

    % get plugin path
    % ---------------
    fullpluginfolder = fullfile(fileparts(which('eeglab.m')), 'plugins', foldername);
    if ~exist(fullpluginfolder)
        error([ 'Could not find folder ' foldername 'in extension folder' ]);
    end;

    disp('Removing extension from path');
    allPaths = path;
    indSep   = [0 find(allPaths == ':') length(allPaths)+1];
    for index = 1:length(indSep)-1
        tmpPath = allPaths(indSep(index)+1:indSep(index+1)-1);
        if ~isempty(strfind(tmpPath, fullpluginfolder))
            rmpath(tmpPath);
            fprintf('Removing path %s\n', tmpPath);
        end;
    end;

    disp([ 'Moving extension ' foldername ' to deactivatedplugins folder' ]);
    fulldeactivatedpluginfolder = fullfile(fileparts(which('eeglab.m')), 'deactivatedplugins');
    if ~exist(fulldeactivatedpluginfolder), mkdir(fulldeactivatedpluginfolder); end;
    try
        movefile(fullpluginfolder, fulldeactivatedpluginfolder);
    catch
        eeglab_error;
    end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function EEG = est_aamp(varargin)

% Obtain the instantaneous analytic amplitude within a given band of a collection of
% processes in EEG.CAT.srcdata. This uses the hilbert transform and EEGLAB's eegfilt().
%
% See Also: eegfilt(), hilbert(), hlp_aamp()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual.
%   Available at: http://www.sccn.ucsd.edu/wiki/SIFT/
%
%
% Author: Tim Mullen 2012, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

g = arg_define([0 1], varargin, ...
        arg_norep({'EEG'},mandatory),...
        arg({'ampband','AmplitudePassBand'},[80 150],[],'The [lo hi] pass-band to use for amplitude (Hz)','shape','row'), ...
        arg_subtoggle({'normalize','NormalizeData'},{'method', {'ensemble'}},@pre_normData,'Data normalization. Normalize analytic amplitude trials across time, ensemble, or both','cat','Normalization'), ...
        arg({'plot','PlotData'},false,[],'Plot amplitude envelope'), ...
        arg({'verb','Verbosity'},true,[],'Verbose output') ...
    );

% yeild EEG structure to workspace
data = hlp_splitstruct(g,{'EEG'});
arg_toworkspace(data);
clear data;

%% hilbert transform the data
AA = hlp_aamp(EEG, 'AmplitudePassBand', g.ampband, 'verb', g.verb);

% normalize analytic amplitude, if desired
if g.normalize.arg_selection
    AA = pre_normData('data',AA,g.normalize,'verb',g.verb);
end

% concatenate analytic amplitudes to end of dataset
EEG.CAT.srcdata = cat(1, EEG.CAT.srcdata, AA);

% update EEG datastructure to reflect new AA 'channels'
EEG.CAT.curComponentNames = EEG.CAT.curComponentNames(EEG.CAT.curComps);
EEG.data       = cat(1, EEG.data, AA);
EEG.CAT.curComps = [EEG.CAT.curComps EEG.CAT.curComps+EEG.CAT.nbchan];
EEG.nbchan     = 2*EEG.CAT.nbchan;
EEG.CAT.nbchan = 2*EEG.CAT.nbchan;

N = length(EEG.CAT.curComponentNames);
for k=1:N
    EEG.CAT.curComponentNames{k+N} = [EEG.CAT.curComponentNames{k} '_amp'];
end
% EEG.CAT.curComponentNames = cellstr(num2str(EEG.CAT.curComps'))';
N = length(EEG.chanlocs);
for k=1:N
    EEG.chanlocs(k+N) = EEG.chanlocs(k);
    EEG.chanlocs(k+N).labels = [EEG.chanlocs(k).labels ' amp'];
end

EEG.CAT.datamode = 'AnalyticAmplitude';
EEG.CAT.ampband = g.ampband;

if g.plot
    eegplot(EEG.CAT.srcdata,'srate',EEG.srate,'title',sprintf('AmplitudeModulation [%0.5g %0.5g] Hz',g.ampband(1),g.ampband(2)));
    ax = findobj(gcf,'tag','eegaxis');
    set(ax,'YTickLabel',flipud([EEG.CAT.curComponentNames,' ']'));
end


text:
%  Create a new page
%
%  ::
%
%     document.newpage();
%
%  Args:
%     none
%
%


text:
%This Matlab script can be used to reproduce Figure 7.20 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%% Define parameters
M = 64; %Number of antennas
Mprime = 30; %Number of
beta = [0,10^(-5/10),1];
IT = 10000;

%Derived parameters
betad = (M-beta*Mprime)/(M-Mprime);

%% Simulate channels
correlationAA = zeros(numel(beta),IT);
correlationAB = zeros(numel(beta),IT);

%Go through all beta values
for b = 1:numel(beta)
    R_A = diag([beta(b)*ones(1,Mprime), betad(b)*ones(1,M-Mprime)]);
    R_B = diag([betad(b)*ones(1,M-Mprime), beta(b)*ones(1,Mprime)]);
    R_A2 = sqrt(R_A);
    R_B2 = sqrt(R_B);

    %Compute UE correlation metric in (7.22)
    for it = 1:IT

        hA1 = R_A2*1/sqrt(2)*(randn(M,1) + 1i*randn(M,1));
        hA2 = R_A2*1/sqrt(2)*(randn(M,1) + 1i*randn(M,1));
        hB1 = R_B2*1/sqrt(2)*(randn(M,1) + 1i*randn(M,1));
        correlationAA(b,it) = 10*log10((abs(hA1'*hA2)/(norm(hA1)*norm(hA2)))^2);
        correlationAB(b,it) = 10*log10((abs(hA1'*hB1)/(norm(hA1)*norm(hB1)))^2);

    end
end

%% Plot the simulation results
figure;
hold on; box on;
nth = 1000;
Colors = {'r','b'};
Markers = {'square', 'o'};

plot(-100,-100,'Color', 'k', 'LineStyle', '-');
plot(-100,-100,'Color', 'k', 'LineStyle', '--');
plot(-100,-100,'Color','k', 'LineStyle', 'none', 'Marker', '*');
plot(-100,-100,'Color', Colors{2}, 'LineStyle', 'none', 'Marker', Markers{2});
plot(-100,-100,'Color', Colors{1}, 'LineStyle', 'none', 'Marker', Markers{1});

[y,x] = ecdf(correlationAA(end,:));
plot(x,y,'Color','k', 'LineWidth', 1);
plot(x(1:nth:end),y(1:nth:end),'Color','k', 'LineStyle', 'none', 'Marker', '*');

for b = 1:numel(beta)-1
    [y,x] = ecdf(correlationAA(b,:));
    plot(x,y,'Color', Colors{b}, 'LineWidth', 1);
end

for b = 1:numel(beta)-1
    [y,x] = ecdf(correlationAB(b,:));
    plot(x,y,'--','Color', Colors{b}, 'LineWidth', 1)
end

for b = 1:numel(beta)-1
    [y,x] = ecdf(correlationAA(b,:));
    plot(x(1:nth:end),y(1:nth:end),'LineStyle', 'none','Color', Colors{b}, 'LineWidth', 1, 'Marker', Markers{b});
end

for b = 1:numel(beta)-1
    [y,x] = ecdf(correlationAB(b,:));
    plot(x(1:nth:end),y(1:nth:end),'LineStyle', 'none','Color', Colors{b}, 'LineWidth', 1, 'Marker', Markers{b})
end

xlim([-40,-10])
legend('Same region', 'Different regions','\beta=1', '\beta=-5 dB', '\beta=0','Location','NorthWest')
xlabel('Average UE correlation [dB]')
ylabel('CDF')

author: emilbjornson

repo: massivemimobook