text:
function[varargout]=ridgemap(varargin)
%RIDGEMAP  Maps ridge quantities back onto the time series.
%
%   X=RIDGEMAP(M,XR,IR) where IR is a ridge index and XR is a quantity
%   along the ridge, maps the values of XR to their correct row locations
%   in a time series of length M, and returns the result in the array X.
%
%   If IR and XR contain L different ridges separated by NaNs, as output by
%   RIDGEWALK, then X is M x L with the values of XR from each ridge in a
%   separate column. Values not specified by the IR are left as NaNs.
%
%   [X1,X2,...,XM]=RIDGEMAP(M,X1R,X2R,...,XPR,IR) also works for any P
%   different ridge quantities X1R--XPR.
%
%   When using RIDGEWALK's joint ridges algorithm, in which some quantities
%   have more than one column, they should be passed to RIDGEMAP
%   individually, for example [X1,X2]=RIDGEMAP(M,X(:,1),X(:,2),IR).
%   __________________________________________________________________
%
%   Collapsing
%
%   X=RIDGEMAP(...'collapse') combines values from all the ridges using
%   a power-weighted mean.  Then X is a column vector of size M x 1.
%   __________________________________________________________________
%
%   Ridge multiplicity
%
%   [...,MULT]=RIDGEMAP returns the ridge multiplicity MULT after all the
%   expected output quantities.  MULT is a column vector with size M x 1.
%
%   The ridge multiplicity is the number of ridges present at each time.
%   __________________________________________________________________
%
%   See also RIDGEWALK.
%
%   'ridgemap --t' runs some a test.
%
%   Usage:   x=ridgemap(M,xr,ir);
%            [x,f]=ridgemap(M,xr,fr,ir);
%            [x,f]=ridgemap(M,xr,fr,ir,'collapse');
%            [x,mult]=ridgemap(M,xr,ir);
%            [x,f,mult]=ridgemap(M,xr,fr,ir);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2009--2019 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
    ridgemap_test,return
end

if ischar(varargin{end})
    str=varargin{end};
    varargin=varargin(1:end-1);
else
    str='all';
end

M=varargin{1}(1);
N=1;
bNinput=false;
if length(M)==2
    N=M(2);
    M=M(1);
    bNinput=true;
end

varargin=varargin(2:end);
if bNinput
    kr=varargin{end};
    ir=varargin{end-1};
    varargin=varargin(1:end-2);
else
    ir=varargin{end};
    kr=1+0*ir;
    varargin=varargin(1:end-1);
end

if ~isempty(ir)
    if iscell(ir)
        ir=cell2col(ir);
        for i=1:length(varargin)
            varargin{i}=cell2col(varargin{i});
        end
    end
    for i=1:length(varargin)
        [~,varargin{i}]=col2mat(ir,varargin{i});
    end
    ir=col2mat(ir);
    for i=1:length(varargin)
        varargout{i}=nan*zeros(M,size(ir,2));
%        size(varargout{i})
 %       size(ir)
        for k=1:size(ir,2)
           %ir(isfinite(ir(:,k)),k)
            varargout{i}(ir(isfinite(ir(:,k)),k),k)=varargin{i}(isfinite(ir(:,k)),k);
        end
    end
    mult=vsum(0+isfinite(varargout{1}),2);
else
    for i=1:length(varargin)
        varargout{i}=nan*zeros(M,1);
    end
    mult=zeros(M,1);
end

varargout{end+1}=mult;

% for i=1:length(varargin)
%     varargout{i}=vswap(varargout{i},inf+sqrt(-1)*inf,nan+sqrt(-1)*nan);
%     varargout{i}=vswap(varargout{i},inf,nan);
% end

%Calculate multiplicity
%mult=~isnan(vswap(varargout{1}(:,1,:),0,nan));
%The "+0" is to convert the logical into a numerical value

if strfind(str,'col')
    for i=length(varargin):-1:1
        if i==1
            varargout{i}=vsum(varargout{i},2);
        else
            varargout{i}=vmean(varargout{i},2,squared(varargout{1}));
        end
    end
end

for i=1:length(varargin)
    varargout{i}=squeeze(varargout{i});
    %anyany(isinf(varargout{i}))
end

function[]=ridgemap_test

load npg2006
use npg2006

%Decide on frequencies
fs=2*pi./(logspace(log10(10),log10(100),50)');

%Compute wavelet transforms using generalized Morse wavelets
[wx,wy]=wavetrans(real(cx),imag(cx),{1,2,4,fs,'bandpass'},'mirror');
[wp,wn]=vectmult(tmat,wx,wy);

%Form ridges of component time series
[wr,ir,jr,fr]=ridgewalk(dt,wn,fs);
[wa,fa,mult]=ridgemap(length(wn),wr,fr,ir);
reporttest('RIDGEMAP has one column per ridge, non-joint ridges',size(fa,2)==size(wa,2)&&size(fa,2)==length(find(~isfinite(ir))))

%[ir,jr,wpr,wnr,fpr,fnr]=ridgewalk(dt,wp,wn,fs,{0,0});
%[wpa,wna,fpa,fna,mult]=ridgemap(length(wn),wpr,wnr,fpr,fnr,ir);

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [cerrFeatS,pyFeatS] = computeRadiomicsFromCERRAndPyrad(planC,structNum,...
    cerrParamFilePath,pyParamFilePath)
%
% function [cerrFeatS,pyFeatS] = computeRadiomicsFromCERRAndPyrad(planC,structNum,...
%     cerrParamFilePath,pyParamFilePath)
%
% This function returns radiomics features from CERR and PyRadiomics.
% planC: CERR's planC data structure.
% structNum: Structure index in planC to compute radiomics.
% cerrParamFilePath: Radiomics settings for CERR.
% pyParamFilePath: Radiomics settings for PyRadiomics.
%
% This routine assumes that PyRadiomics and SciPy are already added to Python path.
% For example, on Windows OS
% P = py.sys.path;
% insert(P,int64(0),'C:\Program Files\Python37\Lib\site-packages\radiomics');
% P = py.sys.path;
% insert(P,int64(0),'C:\Program Files\Python37\Lib\site-packages\scipy');

%------------------------------------------------------------------------
% APA 01/08/2022

%% 1. Compute features using CERR
% cerrParamFilePath = fullfile(fileparts(fileparts(getCERRPath)),...
%             'Unit_Testing/settings_for_comparisons/cerrOrigNoInterp.json');
paramS = getRadiomicsParamTemplate(cerrParamFilePath);

% strC = {planC{indexS.structures}.structureName};
% structNum = getMatchingIndex(paramS.structuresC{1},strC,'exact');
scanNum = getStructureAssociatedScan(structNum,planC);
cerrFeatS = calcGlobalRadiomicsFeatures...
    (scanNum, structNum, paramS, planC);

filtName = 'Original';

%% 2. Compute features using Pyradiomics
% pyParamFilePath = fullfile(fileparts(fileparts(getCERRPath)),...
%             'Unit_Testing/settings_for_comparisons/pyOrigNoInterp.yaml');
% pyParamFilePath = 'L:\Aditya\forJung\TCIA_Jung\features_1_6_2022\pyOrigWithInterp.yaml';
pyCalcS = calcRadiomicsFeatUsingPyradiomics(planC,structNum,pyParamFilePath);

%Map to cerr fieldnames
pyFeatS = struct();

% First-order
pyFirstOrdFeatS = getPyradFeatDict(pyCalcS,{['original','_firstorder']});
pyFirstOrdFeatS = mapPyradFieldnames(pyFirstOrdFeatS,'original','firstorder');
pyFeatS.(filtName).firstOrderS = pyFirstOrdFeatS;

% GLCM
pyGlcmFeatS = getPyradFeatDict(pyCalcS,{['original','_glcm']});
pyGlcmFeatS = mapPyradFieldnames(pyGlcmFeatS,'original','glcm');
pyFeatS.(filtName).glcmFeatS = pyGlcmFeatS;

% GLRLM
pyGlrlmFeatS = getPyradFeatDict(pyCalcS,{['original','_glrlm']});
pyGlrlmFeatS = mapPyradFieldnames(pyGlrlmFeatS,'original','glrlm');
pyFeatS.(filtName).rlmFeatS = pyGlrlmFeatS;

% NGLDM
pyGldmFeatS = getPyradFeatDict(pyCalcS,{['original','_gldm']});
pyGldmFeatS = mapPyradFieldnames(pyGldmFeatS,'original','ngldm');
pyFeatS.(filtName).ngldmFeatS = pyGldmFeatS;

%GLSZM
pyGlszmFeatS = getPyradFeatDict(pyCalcS,{['original','_glszm']});
pyGlszmFeatS = mapPyradFieldnames(pyGlszmFeatS,'original','glszm');
pyFeatS.(filtName).szmFeatS = pyGlszmFeatS;

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function s = isstatv(pc)

%function s = isstatv(pc)
%
% Determines whether the set of partial correlations pc
% defines a stationary process
%
% It checks whether I-pc*pc'>0 for all pc's
%
% It does this by checking whether the singular values
% of all pc's are smaller than 1.
%

Lmax = size(pc,3)-1;
nons = 0;
for i = 1:Lmax,
   nons = nons+any(svd(pc(:,:,i+1))>1);
end
s = ~nons;

author: Sable

repo: mcbench-benchmarks


text:
function centroids = computeCentroids(X, idx, K)
%COMPUTECENTROIDS returs the new centroids by computing the means of the
%data points assigned to each centroid.
%   centroids = COMPUTECENTROIDS(X, idx, K) returns the new centroids by
%   computing the means of the data points assigned to each centroid. It is
%   given a dataset X where each row is a single data point, a vector
%   idx of centroid assignments (i.e. each entry in range [1..K]) for each
%   example, and K, the number of centroids. You should return a matrix
%   centroids, where each row of centroids is the mean of the data points
%   assigned to it.
%

% Useful variables
[m n] = size(X);

% You need to return the following variables correctly.
centroids = zeros(K, n);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every centroid and compute mean of all points that
%               belong to it. Concretely, the row vector centroids(i, :)
%               should contain the mean of the data points assigned to
%               centroid i.
%
% Note: You can use a for-loop over the centroids to compute this.
%

count = zeros(K, 1);
for i = 1 : m
    count(idx(i,1),1) = count(idx(i,1),1) + 1;
    centroids(idx(i,1),:) = centroids(idx(i,1),:) + X(i,:);
end

for i = 1 : K
    centroids(i,:) = centroids(i,:) ./ count(i,1);
end

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

end

author: UtkarshPathrabe


text:
%TestSuiteTest Unit tests for mtest command-line test runner.

classdef MtestTest < TestCaseInDir

   methods

       function self = MtestTest(name)
           self = self@TestCaseInDir(name, ...
               fullfile(fileparts(which(mfilename)), 'cwd_test'));
       end

      function test_noInputArgs(self)
          [T, did_pass] = evalc('mtest');
          % The cwd_test directory contains some test cases that fail,
          % so output of mtest should be false.
          assertFalse(did_pass);
      end

      function test_oneInputArg(self)
          [T, did_pass] = evalc('mtest(''testFoobar'')');
          % cwd_test/testFoobar.m is supposed to pass.
          assertTrue(did_pass);
      end

      function test_oneInputArgWithFilter_passing(self)
          [T, did_pass] = evalc('mtest(''TestCaseSubclass:testA'')');
          assertTrue(did_pass);
      end

      function test_oneInputArgWithFilter_failing(self)
          [T, did_pass] = evalc('mtest(''TestCaseSubclass:testB'')');
          assertFalse(did_pass);
      end

   end

end

author: ricedsp


text:
function node_num = grid_q9_node_num ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_Q9_NODE_NUM counts the nodes in a grid of Q9 elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    NELEMX * NELEMY.
%
%    Output, integer NODE_NUM, the number of nodes in the grid.
%
  node_num = ( 2 * nelemx + 1 ) * ( 2 * nelemy + 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function plotData(x, y)
%PLOTDATA Plots the data points x and y into a new figure
%   PLOTDATA(x,y) plots the data points and gives the figure axes labels of
%   population and profit.

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the training data into a figure using the
%               "figure" and "plot" commands. Set the axes labels using
%               the "xlabel" and "ylabel" commands. Assume the
%               population and revenue data have been passed in
%               as the x and y arguments of this function.
%
% Hint: You can use the 'rx' option with plot to have the markers
%       appear as red crosses. Furthermore, you can make the
%       markers larger by using plot(..., 'rx', 'MarkerSize', 10);

figure; % open a new figure window

plot(x, y, 'rx', 'MarkerSize', 10); % Plot the data
ylabel('Profit in $10,000s'); % Set the y􀀀axis label
xlabel('Population of City in 10,000s'); % Set the x􀀀axis label

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

end

author: vugsus

repo: coursera-machine-learning


text:
function prob_test089 ( )

%*****************************************************************************80
%
%% TEST089 tests HALF_NORMAL_MEAN, HALF_NORMAL_SAMPLE, HALF_NORMAL_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST089\n' );
  fprintf ( 1, '  For the Half Normal PDF:\n' );
  fprintf ( 1, '  HALF_NORMAL_MEAN computes the mean;\n' );
  fprintf ( 1, '  HALF_NORMAL_SAMPLE samples;\n' );
  fprintf ( 1, '  HALF_NORMAL_VARIANCE computes the variance.\n' );

  a = 0.0;
  b = 10.0;

  check = half_normal_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST089 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = half_normal_mean ( a, b );
  variance = half_normal_variance ( a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );
  fprintf ( 1, '  PDF mean =                    %14f\n', mean );
  fprintf ( 1, '  PDF variance =                %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = half_normal_sample ( a, b, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end


text:
function [ y2, m2, d2, f2 ] = ymdf_prev_islamic ( y1, m1, d1, f1 )

%*****************************************************************************80
%
%% YMDF_PREV_ISLAMIC returns the Islamic YMDF date of the previous day.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the YMDF date.
%
%    Output, integer Y2, M2, D2, real F2,
%    yesterday's YMDF date.
%
  y2 = y1;
  m2 = m1;
  d2 = d1 - 1;
  f2 = f1;

  [ y2, m2, d2 ] = day_borrow_islamic ( y2, m2, d2 );

  return
end


text:
function [mxc,mlat,OUT] = xc_stats(xc,xl,varargin)
% [mxc,mlat,OUT] = xc_stats(xc,xl,varargin)
% calculate t-statistic values for a k x k x n 3-D matrix of
% correlation (xc) and latency (xl) values across n subjects.
%
% tor wager

names = [];
if length(varargin) > 0, names = varargin{1};,end

 if isempty(xl), mlat = [];, end

% ------------------------------------------------------------
% stats on correlation values
% ------------------------------------------------------------

warning off, clear stelatency, clear tlat, clear stecorr, clear tcorr
% standard errors and t-values, element by element

% only do this on the upper triangle, to be faster

for j = 1:size(xc,2)-1
    for k = j + 1 : size(xc,2)

        if ~isempty(xl)
            tmp = squeeze(xl(j,k,:));
            z = tmp;
            [stelatency(j,k),tlat(j,k)] = ste(z);   % t-test on latencies; maybe should be poisson dist?
        end

        tmp = squeeze(xc(j,k,:));
        z = .5 * log( (1+tmp) ./ (1-tmp) );     % Fisher's r-to-z transform
        [stecorr(j,k),tcorr(j,k)] = ste(z);     % correl in rand fx analysis across ss
    end
end
warning on

% makes symmetric matrices from upper tri

if ~isempty(xl), stelatency(end+1,:) = 0; ,end
stecorr(end+1,:) = 0;
tcorr(end+1,:) = 0;
if ~isempty(xl),
    tlat(end+1,:) = 0;
    stelatency = stelatency + stelatency';
end

stecorr = stecorr + stecorr' + Inf * eye(size(stecorr,1));
tcorr = tcorr + tcorr';     % t-scores for significance/reliability of cross-correlations across subjects

if ~isempty(xl)
    tlat = tlat + tlat';        % same for time lags across subjects
end

% group means

if ~isempty(xl)
    mlat = mean(xl,3);          % mean of latency data (cross-lag, group average)
end
mxc = mean(xc,3);           % mean cross-correlations among k regions, group average

% ------------------------------------------------------------
% print and save output
% ------------------------------------------------------------

OUT.desc1 = 'xl is latency mtx for each subject, xc is individual correlation mtx';
if ~isempty(xl), OUT.xl = xl;, end
OUT.xc = xc;
OUT.desc2 = 'mlat and mxc: matrices of means across subjects for latency and correl'
OUT.desc3 = 'tlat/tcorr are t-values across z-transformed correlations';
OUT.mxc = mxc;

if ~isempty(xl)
    OUT.mlat = mlat;
    OUT.tlat = tlat;
    OUT.tcorr = tcorr;
end

% Output
if ~isempty(xl)
    fprintf(1,'Latency\n');
    str = correlation_to_text(mlat,Inf,names);
end

fprintf(1,'correlation T-values, p<.05 corrected, assuming normal distribution of latency diffs\n');
% Alpha correction - bonferroni.
numobs=(size(mxc,1)*(size(mxc,1)-1))/2;
corrp=1-(0.05/ (2 * (   numobs   )));       % 2-tailed corr p
crit_t = tinv_t(corrp,size(xc,1)-3);          % critical t-value for corrected significance
crit_tu = tinv_t(1-(.05/2),size(xc,1)-3);         % critical t-value for uncorrected significance

if ~isempty(xl)
    str = correlation_to_text(tlat,crit_t,names);
end

OUT.crit_t = crit_t;
OUT.crit_tu = crit_tu;
OUT.crit_tdesc = 'crit_t crit t value corrected, _tu uncorrected';
OUT.corrp = corrp;
OUT.corrpdesc = 'critical corrected p-value (bonferroni)';

fprintf(1,'Average correlation across subjects\n');
[str, sigu] = correlation_to_text(mxc,Inf,names);

fprintf(1,'Corr T-values, p<.05 uncorrected, 2-tailed, assuming independent samples over time\n');
[str, sigu] = correlation_to_text(tcorr,crit_tu,names);
OUT.sigmat_uncorrected = sigu;

fprintf(1,'Corr T-values, p<.05 corrected, 2-tailed, assuming independent samples over time\n');
[str, sigmat] = correlation_to_text(tcorr,crit_t,names);

OUT.sigmat = sigmat;

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
clear; close all; clc;

%% Path names

toolboxPath=fileparts(fileparts(fileparts(mfilename('fullpath'))));
docPath=fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'docs');
libPath=fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'lib');
docToolsPath=fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'docs','docTools');

pathNames={docToolsPath,libPath,docPath};

%%

fileExtension='.m'; %Extension to remove boilerplate from

footerTargetStart='% _*GIBBON footer text*_ '; %Target for removal

%N.B. Removal is up to end from the start so make sure the target is
%appropriately set!!!!!

%% Removing footer text

numPaths=numel(pathNames);
for q_path=1:1:numPaths
    pathName=pathNames{q_path};
    files = dir(fullfile(pathName,'*.m'));
    files={files(1:end).name};
    files=sort(files(:));
    numFiles=numel(files);
    for q_file=1:1:numFiles
        fileName=fullfile(pathName,files{q_file});
        [T_now]=txtfile2cell(fileName);
        targetStartIndex = find(strcmp(footerTargetStart,T_now));
        if ~isempty(targetStartIndex)
            targetStartIndex=targetStartIndex(end); %Keep last occurance
            T_now=T_now(1:targetStartIndex-1-1); %NB -1 is used to remove %% above target
            cell2txtfile(fileName,T_now,0,0);
        end
    end
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%TestSuiteTest Unit tests for runtests command-line test runner.

classdef RuntestsTest < TestCaseInDir

   methods

       function self = RuntestsTest(name)
           self = self@TestCaseInDir(name, ...
               fullfile(fileparts(which(mfilename)), 'cwd_test'));
       end

      function test_noInputArgs(self)
          [T, did_pass] = evalc('runtests');
          % The cwd_test directory contains some test cases that fail,
          % so output of runtests should be false.
          assertFalse(did_pass);
      end

      function test_Verbose(self)
          [T, did_pass] = evalc('runtests(''-verbose'')');
          assertFalse(did_pass);
      end

      function test_oneInputArg(self)
          [T, did_pass] = evalc('runtests(''testFoobar'')');
          % cwd_test/testFoobar.m is supposed to pass.
          assertTrue(did_pass);
      end

      function test_verboseThenTestName(self)
          [T, did_pass] = evalc('runtests(''-verbose'', ''.'')');
          assertFalse(did_pass);
      end

      function test_testNameThenVerbose(self)
          [T, did_pass] = evalc('runtests(''.'', ''-verbose'')');
          assertFalse(did_pass);
      end

      function test_oneInputArgWithFilter_passing(self)
          [T, did_pass] = evalc('runtests(''TestCaseSubclass:testA'')');
          assertTrue(did_pass);
      end

      function test_oneInputArgWithFilter_failing(self)
          [T, did_pass] = evalc('runtests(''TestCaseSubclass:testB'')');
          assertFalse(did_pass);
      end

      function test_oneDirname(self)
          [T, did_pass] = evalc('runtests(''../dir1'')');
          assertTrue(did_pass);

          [T, did_pass] = evalc('runtests(''../dir2'')');
          assertFalse(did_pass);
      end

      function test_twoDirnames(self)
          [T, did_pass] = evalc('runtests(''../dir1'', ''../dir2'')');
          assertFalse(did_pass);
      end

      function test_packageName(self)
          [T, did_pass] = evalc('runtests(''xunit.mocktests'')');
          assertTrue(did_pass);
      end

      function test_noTestCasesFound(self)
          assertExceptionThrown(@() runtests('no_such_test'), ...
              'xunit:runtests:noTestCasesFound');
      end

      function test_optionStringsIgnored(self)
          % Option string at beginning.
          [T, did_pass] = evalc('runtests(''-bogus'', ''../dir1'')');
          assertTrue(did_pass);

          % Option string at end.
          [T, did_pass] = evalc('runtests(''../dir2'', ''-bogus'')');
          assertFalse(did_pass);
      end

      function test_logfile(self)
          name = tempname;
          command = sprintf('runtests(''../dir1'', ''-logfile'', ''%s'')', name);
          [T, did_pass] = evalc(command);
          assertTrue(did_pass);
          assertTrue(exist(name, 'file') ~= 0);
          delete(name);
      end

      function test_logfileWithNoFile(self)
          assertExceptionThrown(@() runtests('../dir1', '-logfile'), ...
              'xunit:runtests:MissingLogfile');
      end

      function test_logfileWithNoWritePermission(self)
          assertExceptionThrown(@() runtests('../dir1', '-logfile', ...
              'C:\dir__does__not__exist\foobar.txt'), ...
              'xunit:runtests:FileOpenFailed');
      end

      function test_namesInCellArray(self)
          [T, did_pass] = evalc('runtests({''TestCaseSubclass:testA''})');
          assertTrue(did_pass);

          [T, did_pass] = evalc('runtests({''TestCaseSubclass:testA'', ''TestCaseSubclass:testB''})');
          assertFalse(did_pass);
      end

   end

end


text:
%==========================================================================
% This is the testing code of a special case of SRMD (scale factor = 1) for real image <denoising & deblurring>.
% There are two models, "SRMDx1_gray.mat" for grayscale image, "SRMDx1_color.mat"
% for color image. The models can do:
%   1. Deblurring. (The kernel is assumed to be Gaussian-like!!! For other kernels, you should re-train the model!)
%      there are two types of kernels,
%      including isotropic Gaussian (width range: [0.1, 3]),
%      anisotropic Gaussian ([0.5, 8]).
%   2. Denoising. the noise level range is [0, 75].
%      For denoising only, set "kerneltype = 1; kernelwidth = 0.1." (i.e., delta kernel)
%
%==========================================================================
% The basic idea of SRMD is to learn a CNN to infer the MAP of general SISR (with special case of sf=1), i.e.,
% solve x^ = arg min_x 1/(2 sigma^2) ||kx - y||^2 + lamda \Phi(x)
% via x^ = CNN(y,k,sigma;\Theta) or x^ = CNN(y,kernel,noiselevel;\Theta).
%
% There involves two important factors, i.e., blur kernel (k; kernel) and noise
% level (sigma; nlevel).
%
% For more information, please refer to the following paper.
%    @article{zhang2017learningsrmd,
%    title={Learning a Single Convolutional Super-Resolution Network for Multiple Degradations},
%    author={Kai, Zhang and Wangmeng, Zuo and Lei, Zhang},
%    year={2017},
%    }
%
% If you have any question, please feel free to contact with <Kai Zhang (cskaizhang@gmail.com)>.
%
% This code is for research purpose only.
%
% by Kai Zhang (Nov, 2017)
%==========================================================================

% clear; clc;
format compact;

addpath('utilities');
imageSets    = {'Hilbert','Set5','Set14','BSD68','BSD100','Urban100'}; % testing dataset

%%======= ======= ======= degradation parameter settings ======= ======= =======

% For real image 'Hilbert', some examples of degradation setting are given as follows.
% sf = 1; nlevel =  8; kerneltype = 1; kernelwidth = 0.2; % denoising only
% sf = 1; nlevel = 10; kerneltype = 1; kernelwidth = 0.5; % denoising and sharpening(deblurring)

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

%% select testing dataset, use GPU or not, ...
setTest      = imageSets([1]); %
showResult   = 1; % 1, show results; 2, save restored images
pauseTime    = 1;
useGPU       = 1; % 1 or 0, true or false
method       = 'SRMD';
folderTest   = 'testsets';
folderResult = 'results';
if ~exist(folderResult,'file')
    mkdir(folderResult);
end

%% scale factor (it is fixed to 1)

sf          = 1; %{1}

%% load model with scale factor sf
folderModel = 'models';
load(fullfile(folderModel,['SRMDx',int2str(sf),'_gray.mat']));
%net.layers = net.layers(1:end-1);
net = vl_simplenn_tidy(net);
if useGPU
    net = vl_simplenn_move(net, 'gpu') ;
end

%% degradation parameter (noise level and kernel) setting
%############################# noise level ################################
% noise level, from a range of [0, 75]

nlevel     = 10;  % [0, 75]

kerneltype = 1;  % {1}

%############################### kernel ###################################
% there are tree types of kernels, including isotropic Gaussian,
% anisotropic Gaussian, and estimated kernel k_b for isotropic Gaussian k_d
% under direct downsampler (x2 and x3 only).

if kerneltype == 1
    % type 1, isotropic Gaussian---although it is a special case of anisotropic Gaussian.
    kernelwidth = 0.5; % from a range of [0.1, 3]. set kernelwidth from (0.001, 0.2) to generate delta kernel (no blur)
    kernel = fspecial('gaussian',15, kernelwidth); % Note: the kernel size is fixed to 15X15.
    tag    = ['_',method,'_x',num2str(sf),'_itrG_',int2str(kernelwidth*10),'_nlevel_',int2str(nlevel)];

elseif kerneltype == 2
    % type 2, anisotropic Gaussian
    nk     = 1;%randi(size(net.meta.AtrpGaussianKernel,4)); % randomly select one
    kernel = net.meta.AtrpGaussianKernel(:,:,:,nk);
    tag    = ['_',method,'_x',num2str(sf),'_atrG_',int2str(nk),'_nlevel_',int2str(nlevel)];

end

%##########################################################################

surf(kernel) % show kernel
view(45,55);
title('Assumed kernel');
xlim([1 15]);
ylim([1 15]);
pause(2)
close;

%% for degradation maps
global degpar;
degpar = single([net.meta.P*kernel(:); nlevel(:)/255]);

for n_set = 1 : numel(setTest)

    %% search images
    setTestCur = cell2mat(setTest(n_set));
    disp('--------------------------------------------');
    disp(['    ----',setTestCur,'-----Super-Resolution-----']);
    disp('--------------------------------------------');
    folderTestCur = fullfile(folderTest,setTestCur);
    ext                 =  {'*.jpg','*.png','*.bmp'};
    filepaths           =  [];
    for i = 1 : length(ext)
        filepaths = cat(1,filepaths,dir(fullfile(folderTestCur, ext{i})));
    end

    %% prepare results
    folderResultCur = fullfile(folderResult, [setTestCur,tag]);
    if ~exist(folderResultCur,'file')
        mkdir(folderResultCur);
    end

    %% perform denoising or/and deblurring (only support Gaussian-like kernel)
    for i = 1 : length(filepaths)

        input  = imread(fullfile(folderTestCur,filepaths(i).name));
        %label  = modcrop(label, 2);
        [h,w,C]   = size(input);
        input     = im_pad(input);
        if C == 3
            input = rgb2gray(input); % input = rgb2ycbcr(input); input = input(:,:,1); % another option
        end

        [~,imageName,ext] = fileparts(filepaths(i).name);
        input1 = im2single(input);
        %tic
        if useGPU
            input1 = gpuArray(input1);
        end
        res = vl_srmd(net, input1,[],[],'conserveMemory',true,'mode','test','cudnn',true);
        %res = vl_srmd_concise(net, input1); % a concise version of "vl_srmd".
        %res = vl_srmd_matlab(net, input1); % You should also set "useGPU = 0;" and comment "net = vl_simplenn_tidy(net);"

        output = im2uint8(gather(res(end).x));

        %toc;
        output = im_crop(output,h,w);
        input  = im_crop(input,h,w);
        %output2 = 0.8*output + 0.2*input; % add noise and structure back to make the output more visually plausible. or GAN?

        disp([setTestCur,'    ',int2str(i),'    ',filepaths(i).name]);
        if showResult
            imshow(cat(2,input,output));
            drawnow;
            title(['Denoising and deblurring     ',filepaths(i).name],'FontSize',12)
            pause(pauseTime)
            imwrite(output,fullfile(folderResultCur,[imageName,'_x',int2str(sf),'.png']));% save results

        end
    end
end

author: cszn

repo: SRMD

save_path: github-repos/MATLAB/cszn-SRMD


text:
function a = dorr_inverse ( alpha, n )

%*****************************************************************************80
%
%% DORR_INVERSE returns the inverse of the DORR matrix.
%
%  Discussion:
%
%    The DORR matrix is a special case of the TRIV matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    CM daFonseca, J Petronilho,
%    Explicit Inverses of Some Tridiagonal Matrices,
%    Linear Algebra and Its Applications,
%    Volume 325, 2001, pages 7-21.
%
%  Parameters:
%
%    Input, real ALPHA, the parameter.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(N,N), the inverse of the matrix.
%

%
%  Form the three diagonals.
%
  x = zeros ( n - 1, 1 );
  y = zeros ( n,     1 );
  z = zeros ( n - 1, 1 );

  for i = 1 : n
    for j = 1 : n

      if ( i <= floor ( ( n + 1 ) / 2 ) )

        if ( j == i - 1 )
          x(i-1) = - alpha * ( n + 1 )^2;
        elseif ( j == i )
          y(i) = 2.0 * alpha * ( n + 1 )^2 + 0.5 * ( n + 1 ) - i;
        elseif ( j == i + 1 )
          z(i) = - alpha * ( n + 1 )^2 - 0.5 * ( n + 1 ) + i;
        end

      else

        if ( j == i - 1 )
          x(i-1) = - alpha * ( n + 1 )^2 + 0.5 * ( n + 1 ) - i;
        elseif ( j == i )
          y(i) = 2.0 * alpha * ( n + 1 )^2 - 0.5 * ( n + 1 ) + i;
        elseif ( j == i + 1 )
          z(i) = - alpha * ( n + 1 )^2;
        end

      end

    end
  end
%
%  Now evaluate the inverse.
%
  a = zeros ( n, n );

  d = zeros ( n, 1 );
  e = zeros ( n, 1 );

  d(n) = y(n);
  for i = n - 1 : -1 : 1
    d(i) = y(i) - x(i) * z(i) / d(i+1);
  end

  e(1) = y(1);
  for i = 2 : n
    e(i) = y(i) - x(i-1) * z(i-1) / e(i-1);
  end

  for i = 1 : n
    for j = 1 : i
      a(i,j) = r8_mop ( i + j ) * prod ( x(j:i-1) ) ...
        * prod ( d(i+1:n) ) / prod ( e(j:n) );
    end
    for j = i + 1 : n
      a(i,j) = r8_mop ( i + j ) * prod ( z(i:j-1) ) ...
        * prod ( d(j+1:n) ) / prod ( e(i:n) );
    end
  end

  return
end


text:
function [masks,keepLocs] = maskGaussians( siz, M, width, offset, show )
% Divides a volume into softly overlapping gaussian windows.
%
% Return M^nd masks each of size siz.  Each mask represents a symmetric
% gaussian window, the locations are evenly spaced throughout the array of
% size siz.  For example, if M=2, then along each dimension d the location
% of each mask is either 1/4 or 3/4 of siz(d) and likewise if M=3 that mask
% is at 1/6,3/6, or 5/6 of siz(d). For higher M the locations are:
% 1/2M,3/2M,...,M-1/2M. See examples below to visualize the masks.
%
% The std of each gaussian is set to be equal to the spacing between two
% adjacent masks multiplied by width.  Reducing the widths of the gaussians
% causes there to be less overlap between masks, but if the width is
% reduced too far certain pixels in between masks receive very little
% weight.  A desired property of the masks is that their coverage (the
% total weight placed on the pixel by all the masks) is approximately
% constant.  Typically, we settle for having the coverage be monotonically
% decreasing as we move away from the center.  (In reality the coverage
% oscilates as we move past peaks, it's just that the oscillations tend to
% be small).  The default value of the width is .6, which minimizes overlap
% while still providing good overall coverage. Values lower tend to produce
% noticeable oscillations in coverage.  offset in (-.5,1) controls the
% spacing of the locations.  Essentially, a positive offset moves the
% locations away from the center of the array and a negative offset moves
% the windows away from the center.  Using a positive offset gives better
% coverage to areas near the borders.
%
% USAGE
%  [masks,keepLocs] = maskGaussians( siz, M, [width], [offset], [show] )
%
% INPUTS
%  siz         - dimensions of each mask
%  M           - # mask locations along each dim [either scalar or vector]
%  width       - [.6] widths of the gaussians
%  offset      - [.1] spacing of mask centers; in (-.5,1)
%  [show]      - [0] figure to use for display (no display if==0) (nd<=3)
%
% OUTPUTS
%  masks       - [see above] array of size [siz x M^nd]
%  keepLocs    - logical array of all locs where masks is nonzero
%
% EXAMPLE
%  masks = maskGaussians( 100, 10, .6, -.1, 1 );  %1D
%  masks = maskGaussians( [35 35], 3, .6, .1, 1 );  %2D
%  masks = maskGaussians( [35 35 35], [2 2 4], .6, .1, 1 ); %3D
%
% See also HISTCIMLOC, MASKCIRCLE
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

nd = length(siz);
if( nargin<3 || isempty(width)); width = .6; end;
if( nargin<4 || isempty(offset)); offset = .1; end;
if( nargin<5 || isempty(show) || nd>3 ); show = 0; end;

%%% uses a cache since this is slow but often called with same inputs.
persistent cache; if( isempty(cache) ); cache=simpleCache('init'); end;
key = [nd siz M width offset];
[found,val] = simpleCache( 'get', cache, key );
if( found ) %%% get masks and keepLocs from cache
  [masks,keepLocs] = deal(val{:});
else %%% create masks and keepLocs
  [M,er] = checkNumArgs( M, [1 nd], 0, 2 ); error(er);
  inds = {':'}; inds = inds(:,ones(1,nd));
  if( offset<=-.5 || offset>=1 ); error('offset must be in (-.5,1)'); end;

  %%% the covariance of each window
  spacing = (siz*(1+2*offset))./M;
  sigmas = spacing * width;
  C = diag(sigmas.^2);

  %%% create each mask
  masks = zeros( [siz,prod(M)] );
  for c=1:prod(M)
    sub = ind2sub2( M, c );
    mus = (sub-.5).* spacing + .5-offset*siz;
    masks(inds{:},c) = filterGauss( siz, mus, C );
  end
  keepLocs = masks>1e-7;

  %%% place into cache
  cache = simpleCache( 'put', cache, key, {masks,keepLocs} );
end;

%%% optionally display
if( show )
  if( nd==1 )
    figure(show); clf; plot( masks );
    figure(show+1); clf; plot(sum( masks,nd+1 ));
    a=axis; a(3)=0; axis(a);
    title('coverage');
 elseif( nd==2)
    figure(show); clf; montage2( masks );
    figure(show+1); clf; im(sum( masks,nd+1));
    title('coverage');
  elseif( nd==3)
    figure(show); clf; montage2( masks );
    figure(show+1); clf; montage2(sum( masks,nd+1) );
    title('coverage');
  end
end;

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function tcell_flow_movie ( )

%*****************************************************************************80
%
%% TCELL_FLOW_MOVIE animates the TCELL data.
%
%  Discussion:
%
%    This MATLAB script file reads the T-Cell flow data:
%
%      geometry (XY values at nodes, assumed to be in 'xy.txt')
%      flow (UV values at nodes, in a sequence of files starting with 'up001.txt')
%
%    and plots the velocity vectors (U,V)(X,Y), saving each plot and
%    making an animation.
%
%    The file plots either the velocity vector field, or the velocity
%    direction field, depending on the value of the internal logical
%    parameter "normalized".
%
%    The MATLAB routine quiver internally scales the vectors, but this can
%    be adjusted by using a value of SCALE that is not 1.
%
%    For the unnormalized case, the velocity field is computed to a fixed
%    scale determined by finding the largest velocity vector over time.
%
%    This routine requires some auxiliary routines in order to "increment"
%    the name of the current velocity data file to get the name of the next
%    one, with the main such routine being called FILE_NAME_INC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 July 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TCELL_FLOW_MOVIE:\n' );
  fprintf ( 1, '  Create a movie of flow in the TCELL.\n' );

  FALSE = 0;
  TRUE = 1;
  scale = 1.0;
  normalized = TRUE;

  nframes = 500;
%
%  Get the XY coordinates of the nodes.
%
  load xy.txt;

  x = xy(:,1);
  y = xy(:,2);
%
%  Compute the thinning index.
%
  thin_factor = 4;
  thin_dex = thin_index ( x, y, thin_factor );
  thin_num = length ( thin_dex );
  x = x(thin_dex);
  y = y(thin_dex);
%
%  For unnormalized plots, you need to make sure that a fixed scale is preserved from
%  step to step.  The only way I can see to do this requires that I determine the
%  maximum velocity over all time, and then append one extra node and velocity
%  to the data structure.
%
  if ( normalized == FALSE )

    upfile = 'up000.txt';
    vnorm_max = 0.0;

    for i = 1 : nframes
      i
      upfile = file_name_inc ( upfile );
      uv = load ( upfile );
      u = uv(thin_dex,1);
      v = uv(thin_dex,2);
      norm = sqrt ( u.^2 + v.^2 );
      vnorm_max = max ( vnorm_max, max ( norm ) );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Maximum visible velocity magnitude is %f\n', vnorm_max );

  end
%
%  Set the coordinates of the boundary lines, and an extra set of
%  lines that will be invisible, but which make some space for the
%  file name to be displayed within the plot area.
%
  bx1 = [ 0.00, 1.00, 1.00, 0.00, 0.00, 0.00 ];
  by1 = [ 1.00, 1.00, 0.99, 0.99, 0.99, 1.00 ];

  bx2 = [ 0.00, 0.25, 0.25, 0.75, 0.75, 1.00, 1.00, 0.74, 0.74, ...
          0.26, 0.26, 0.00, 0.00];
  by2 = [ 0.50, 0.50, 0.00, 0.00, 0.50, 0.50, 0.51, 0.51, 0.01, ...
          0.01, 0.51, 0.51, 0.50];

  bx3 = [ -0.10,  1.10, 1.10, -0.10, -0.10 ];
  by3 = [ -0.10, -0.10, 1.10,  1.10, -0.10 ];
%
%  Set the name of the 0th (nonexistent) velocity file.
%  All the velocity files will have this format, with the
%  numeric part of the name incremented to get the next one.
%  In particular, the first velocity file is called UP001.TXT.
%
  upfile = 'up000.txt';

  for i = 1 : nframes

    i
    upfile = file_name_inc ( upfile );
    uv = load ( upfile );
    u = uv(thin_dex,1);
    v = uv(thin_dex,2);

    if ( normalized == TRUE )
      norm = sqrt ( u.^2 + v.^2 );
      nonzero = find ( norm ~= 0.0 );
      u(nonzero) = u(nonzero) ./ norm(nonzero);
      v(nonzero) = v(nonzero) ./ norm(nonzero);
      vnorm_max = 1.0;
    end

    x(thin_num+1) = 0.00;
    y(thin_num+1) = 1.05;
    u(thin_num+1) = vnorm_max;
    v(thin_num+1) = 0.0;

    quiver ( x, y, u, v, scale )
%
%  Draw the boundary lines, and the invisible bounding line.
%
    line ( bx1, by1, 'color', 'r' )
    line ( bx2, by2, 'color', 'r' )
    line ( bx3, by3, 'color', 'w' )
    axis equal
    if ( normalized == TRUE )
      title ( 'T-Cell Direction Field' )
    else
      title ( 'T-Cell Flow Field' )
    end
    text ( 0.420, 1.03, upfile )

    my_frames(:,i) = getframe;

  end
%
%  Display the movie a given number of times.
%
  movie ( my_frames, 1 )

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TCELL_FLOW_MOVIE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function file_name = file_name_inc ( file_name )

%*****************************************************************************80
%
%% FILE_NAME_INC generates the next filename in a series.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP!         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the string to be incremented.
%
%    Output, string FILE_NAME, the incremented string.
%
  lens = length ( file_name );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_NAME_INC - Fatal error!\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILE_NAME_INC - Fatal error!' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = file_name(i);

    if ( '0' <= c & c <= '8' )

      change = change + 1;

      c = c + 1;

      file_name(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';

      file_name(i) = c;

    end

  end

  if ( change == 0 )
    file_name = ' ';
  end

  return
end
function thin_dex = thin_index ( x, y, thin_factor )

%*****************************************************************************80
%
%%  THIN_INDEX determines thinning indices for a X, Y data.
%
%  Discussion:
%
%    A set of X, Y data is given, that is presumably, not too far off
%    from being on a rectangular grid.
%
%    The input value of THIN_FACTOR indicates by how much the data should
%    be thinned.
%
%    The X and Y ranges are computed, and only those data points are
%    retained for which both X and Y lie in an appropriate subrange.
%
%    For instance, a THIN_FACTOR of 2 would essentially save data
%    that lay in the black squares of a checkerboard.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(NODE_NUM), Y(NODE_NUM), the X and Y coordinates
%    of the nodes.
%
%    Input, integer THIN_FACTOR, the thinning factor.
%
%    Output, integer THIN_DEX(NODE_NUM), contains in (1:THIN_NUM) the
%    indices into X and Y of the vectors to be retained after thinning.
%
  TRUE = 1;
  FALSE = 0;

  node_num = length ( x );

  x_unique_num = 0;

  for i = 1 : node_num

    unique = TRUE;

    for j = 1 : x_unique_num
      if ( x(i) == x_unique(j) )
        unique = FALSE;
        break;
      end
    end

    if ( unique )
      x_unique_num = x_unique_num + 1;
      x_unique(x_unique_num) = x(i);
    end

  end

  sort ( x_unique );

  y_unique_num = 0;

  for i = 1 : node_num

    unique = TRUE;

    for j = 1 : y_unique_num
      if ( y(i) == y_unique(j) )
        unique = FALSE;
        break;
      end
    end

    if ( unique )
      y_unique_num = y_unique_num + 1;
      y_unique(y_unique_num) = y(i);
    end

  end

  sort ( y_unique );

  thin_num = 0;

  for i = 1 : node_num

    for j = 1 : x_unique_num-1
      if ( x_unique(j) <= x(i) & x(i) <= x_unique(j+1) )
        x_bin = j;
        break;
      end
    end

    for j = 1 : y_unique_num-1
      if ( y_unique(j) <= y(i) & y(i) <= y_unique(j+1) )
        y_bin = j;
        break;
      end
    end

    if ( mod ( y_bin, thin_factor ) == thin_factor / 2 & ...
         mod ( x_bin, thin_factor ) == thin_factor / 2 )

      thin_num = thin_num + 1;
      thin_dex(thin_num) = i;

    end

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TCELL_FLOW_MOVIE\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Problem 8
% Solve the difference equation
% y[n]+1.5y[n-1]+0.5y[n-2]=x[n]+x[n-1] , x[n]=0.8^n u[n]

syms n z Y
x=0.8^n;
X=ztrans(x,z);
X1=z^(-1)*X;
Y1=z^(-1)*Y;
Y2=z^(-2)*Y;
G=Y+1.5*Y1+0.5*Y2-X-X1;
SOL=solve(G,Y);
y=iztrans(SOL,n)

% a)
n_s=0:20;
y_s=subs(y,n,n_s);
stem(n_s,y_s);
legend('Solution y[n]')
xlim([-.5 20.5])
ylim([0 1.1])

% b)
xn=x;
xn_1=0.8^(n-1);
yn=y;
yn_1=subs(y,n,n-1);
yn_2=subs(y,n,n-2);
test=yn+1.5*yn_1+0.5*yn_2-xn-xn_1
simplify(test)


text:
function x = power( x, y )

% POWER     Matrix power, z = x.^y

if isa(y,'tfocs_tuple')
	x.value_ = cellfun( @power, x.value_, y.value_, 'UniformOutput', false );
elseif isscalar(y)
	x.value_ = cellfun( @power, x.value_, {y}, 'UniformOutput', false );
else
	x.value_ = cellfun( @power, x.value_, {y}, 'UniformOutput', false );
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function installUnfold(flatSubdir,query)
%
% function installUnfold([flatSubdir],[query])
%
% Rebuild coords after a new unfold.
%
%    flatSubdir: used to specify the flat subdirectory in case there is
%    more than one unfold available. Default: prompt user to choose from a menu.
%
% Note: could automatically rebuild Gray/corAnals by xforming from inplane, but we
% would have to do this for each dataType.
%
% 12/8/00 Wandell and Huk
% djh, 2/14/2001
%    updated to use hiddenFlat instead of opening a window

if ~exist('flatSubdir','var')
	flatSubdir = getFlatSubdir;
end

if ~exist('query','var')
    query = 1;
end

if query
    resp = questdlg(['Delete ',flatSubdir,'/anat, coords, corAnal, ' ...
                     'and parameter map files and close flat window(s)?']);
else
    resp = 'Yes';
end

%delete coords and rebuild
if strcmp(resp,'Yes')
    closeAllFlatWindows(flatSubdir);
    cleanFlat(flatSubdir);

    % Open hidden gray structure because we need the gray coords
    hiddenGray = initHiddenGray;

    % Load flat file, compute & save coords
    hiddenFlat = initHiddenFlat(flatSubdir);

    % Compute and save flat anat file
    hiddenFlat = loadAnat(hiddenFlat);
    disp('Rebuilt Flat coords and anat files. Deleted old corAnals.');

else
    disp('Coords and corAnals not deleted.');

end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function d = bayesgauss(X, CA, MA, P)
%BAYESGAUSS Bayes classifier for Gaussian patterns.
%   D = BAYESGAUSS(X, CA, MA, P) computes the Bayes decision
%   functions of the n-dimensional patterns in the rows of X.
%   CA is an array of size n-by-n-by-W containing W covariance
%   matrices of size n-by-n, where W is the number of classes.
%   MA is an array of size n-by-W, whose columns are the corres-
%   ponding mean vectors. A cov. matrix and a mean vector must be
%   specified for each class, even is some are equal.  X is of size
%   K-by-n, where K is the number of patterns to be classified. P is
%   a 1-by-W array, containing the probabilities of occurrence of
%   each class.  If P is not included in the argument, the classes
%   are assumed to be equally likely.
%
%   D, is a column vector of length K. Its ith element is the class
%   number assigned to the ith vector in X during classification.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.10 $  $Date: 2004/12/15 20:15:38 $

d = [ ]; % Initialize d.
error(nargchk(3, 4, nargin)) % Verify correct no. of inputs.
n = size(CA, 1);             % Dimension of patterns.

% Protect against the possibility that the class number is
% included as an (n+1)th element of the vectors.
X = double(X(:, 1:n));
W = size(CA, 3); % Number of pattern classes.
K = size(X, 1);  % Number of patterns to classify.
if nargin == 3
   P(1:W) = 1/W; % Classes assumed equally likely.
else
   if sum(P) ~= 1
      error('Elements of P must sum to 1.');
   end
end
% Compute the determinants.
for J = 1:W
   DM(J) = det(CA(:, :, J));
end

% Evaluate the decision functions. Note the use of
% function mahalanobis discussed in Section 12.2.
MA = MA'; % Organize the mean vectors as rows.
for J = 1:W
   C = CA(:,:,J);
   M = MA(J,:);
   k1 = log(P(J));
   k2 = 0.5*log(DM(J));
   L(1:K,1) = k1;
   DET(1:K,1) = k2;
   if P(J) == 0;
      D(1:K, J) = -inf;
   else
      D(:, J) = L - DET - 0.5*mahalanobis(X, C, M);
   end
end

% Find the coordinates of the maximum value in each row. These maxima
% maxima give the class of each pattern:
[i, j] = find(D==repmat(max(D')',1,size(D,2)));
% Re-use X. It now contains the max value along each column.
X = [i j];
% Eliminate multiple classifications of the same patterns. Since
% the class assignment when two or more decision functions give
% the same value is arbitrary, we need to keep only one.
X = sortrows(X);
[b, m] = unique(X(:,1));
X = X(m, :);
% X is now sorted, with the 2nd column giving the class of the
% pattern no. in the 1st col.;  i.e., X(j, 1) refers to the jth
% input pattern, and X(j, 2) is its the class number.

% Output the result of classification. d is a col. vector with
% length equal to the total no. of inout patterns. The elements of
% d are the classes into which the patterns were classified.
d = X(:,2);

repo: weiyanmin


text:
function x0 = bsv_crossing ( a, b, n, x, u )

%*****************************************************************************80
%
%% BSV_CROSSING estimates the location X0 where U(X0) = 0.
%
%  Discussion:
%
%    This function is intended for a special case, in which we are solving
%    the Burgers equation over an interval [A,B], with positive boundary
%    condition at one end and negative at the other.  In that case, the
%    solution must change be zero at some point X0 in the interval.
%
%    We assume the solution is discretized by a piecewise linear function.
%    We use binary search to locate consecutive indices I and I+1 so that
%    U(I) and U(I+1) are of opposite signs.  We use linear interpolation to
%    estimate the location of X0 between X(I) and X(I+1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the left and right endpoints.
%
%    Input, integer N, the number of nodes to use between A and B.
%
%    Input, real X(N), the node coordinates.
%
%    Input, real U(N), the computed discretized solution.
%    It must be the case that U(1) and U(N) are of opposite sign.
%
%    Output, real X0, a point where the piecewise linear approximation to
%    the solution is zero.
%
  if ( u(1) == 0.0 )
    x0 = x(1);
    return
  end

  if ( u(n) == 0.0 )
    x0 = x(n);
    return
  end

  if ( r8_sign ( u(1) ) == r8_sign ( u(n) ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BSV_CROSSING - Fatal error!\n' );
    fprintf ( 1, '  We require opposite signs for U(1) and U(N).\n' );
    error ( 'BSV_CROSSING - Fatal error!\n' );
  end
%
%  Set the initial change of sign indices.
%
  i = 1;
  k = n;
%
%  Do a binary search for the smallest change of sign interval.
%
  while ( 1 < k - i )
%
%  Try halfway.
%
    j = floor ( ( i + k ) / 2 );

    if ( u(j) == 0.0 )
      x0 = x(j)
      return
    end

    if ( r8_sign ( u(j) ) == r8_sign ( u(i) ) )
      i = j;
    else
      k = j;
    end

  end
%
%  The change of sign interval indices are I and K=I+1.
%  Now interpolate to get X0.
%
  x0 = x(i) + u(i) * ( x(i) - x(k) ) / ( u(k) - u(i) );

  return
end


text:
%% DNN Face Detection and Recognition
% This tutorial will show us how to run deep learning models, with face
% detection and face recognition models pipeline.
%
% Sources:
%
% * <https://docs.opencv.org/3.4.1/d5/d86/tutorial_dnn_javascript.html>
% * <https://github.com/opencv/opencv/blob/3.4.1/samples/dnn/js_face_recognition.html>
%
%% Face detection
% Face detection network gets BGR image as input and produces set of bounding
% boxes that might contain faces. All that we need is just select the boxes
% with a strong confidence.
%
% Face detector is based on SSD framework (Single Shot MultiBox Detector),
% using a reduced ResNet-10 model.
%
%% Face recognition
% Network is called <https://github.com/cmusatyalab/openface OpenFace>. Face
% recognition model receives RGB face image of size |96x96|. Then it returns
% 128-dimensional unit vector that represents input face as a point on the
% unit multidimensional sphere. So difference between two faces is an angle
% between two output vectors.
%
%% Code
% Start the demo, then press "Add a person" to name a person that is
% recognized as an unknown one.
%

%%
% import deep learning models
[netDet, blobDetOpts] = ResNetSSD_FaceDetector();
[netRec, blobRecOpts] = OpenFace_Embedding();
assert(~netDet.empty() && ~netRec.empty());

%%
% options
confThreshold = 0.5;  % minimum confidence for face detection
scoreThreshold = 0.5; % minimum score for face recognition

%%
% prepare video input
cap = cv.VideoCapture();
pause(1);
assert(cap.isOpened(), 'Could not initialize capturing');
frame = cap.read();
assert(~isempty(frame), 'Could not read frame');

%%
% prepare figure
hThumb = [];
hImg = imshow(frame);
hFig = ancestor(hImg, 'figure');
setappdata(hFig, 'flag',false);
hBtn = uicontrol('Parent',hFig, 'Style','pushbutton', ...
    'Position',[20 20 100 20], 'String','Add a person', ...
    'Callback',@(~,~) setappdata(hFig, 'flag',true));

%%
% dataset of extracted faces and corresponding names
vecs = zeros(128,0,'single');        % feature vectors
names = cell(1,0);
faces = zeros([96 96 3 0], 'uint8'); % only needed for visualization

%%
% main loop that receives a frames from a camera and makes a recognition of
% every detected face on the frame
while ishghandle(hImg)
    % read frame
    frame = cap.read();
    if isempty(frame), break; end
    out = frame;

    % detect faces
    rects = detectFaces(frame, netDet, blobDetOpts, confThreshold);
    for i=1:size(rects,1)
        % preprocess face
        rect = rects(i,:);
        face = alignFace(cv.Rect.crop(frame, rect));

        % recognize face
        vec = face2vec(face, netRec, blobRecOpts);
        [name, score] = recognizeFace(vec, vecs, names, scoreThreshold);

        % show detection and prediction
        out = insertAnnotation(out, rect, sprintf('%s (%.2f)', name, score), ...
            'Color',name2clr(name, names), 'TextColor',[255 255 255], ...
            'Thickness',2, 'FontScale',0.9);
    end

    % update plot
    set(hImg, 'CData',out);
    drawnow;

    % check if add-a-person button is pressed
    flag = ishghandle(hFig) && getappdata(hFig, 'flag');
    if flag && ~isempty(rects)
        setappdata(hFig, 'flag',false);

        % prompt for name
        name = inputdlg('Enter person name:');
        name = strtrim(name{1});

        % face representation as feature vector
        rect = rects(1,:);
        face = alignFace(cv.Rect.crop(frame, rect));
        vec = face2vec(face, netRec, blobRecOpts);

        % store
        vecs(:,end+1) = vec;
        names{end+1} = name;

        % visualize face + name
        face = cv.resize(face, [96 96]);
        face = cv.putText(face, name, [5 25], ...
            'Color',name2clr(name, names), 'FontScale',0.6, 'Thickness',2);
        faces(:,:,:,end+1) = face;

        % show montage of tracked people
        if ishghandle(hThumb)
            clf(hThumb, 'reset');
            figure(hThumb)
        else
            hThumb = figure;
        end
        montage(faces, 'Size',[NaN 2]);
        movegui(hThumb, 'east')
    end
end
cap.release();

%%
% Helper functions

function [rects, confs] = detectFaces(img, net, blobOpts, thresh)
    %DETECTFACES  Run face detection network to detect faces on input image
    %
    % You may play with input blob sizes to balance detection quality and
    % efficiency. The bigger input blob the smaller faces may be detected.
    %

    % detect faces
    net.setInput(cv.Net.blobFromImages(flip(img,3), blobOpts{:}));
    dets = net.forward();

    % SSD output is 1-by-1-by-ndetections-by-7
    % d = [img_id, class_id, confidence, left, top, right, bottom]
    dets = permute(dets, [3 4 2 1]);

    % filter out weak detections
    if nargin < 4, thresh = 0.5; end
    idx = (dets(:,2) == 1 & dets(:,3) > thresh);  % 0: background, 1: face
    dets = dets(idx,:);

    % adjust relative coordinates to image size
    sz = [size(img,2) size(img,1)];
    dets(:,4:7) = bsxfun(@times, dets(:,4:7), [sz sz]);

    % output detections (clamp coords and remove small and out-of-bound rects)
    rects = cv.Rect.from2points(dets(:,4:5), dets(:,6:7));
    rects = cv.Rect.intersect(rects, [0 0 sz]);
    idx = (cv.Rect.area(rects) >= 10);
    rects = rects(idx,:);
    confs = dets(idx,3);
end

function img = alignFace(img)
    %ALIGNFACE  Align face to make the eyes and bottom lip appear in the same location
    %
    % OpenFace expects faces to be aligned, it uses Dlib.
    %

    %TODO: not implemented, maybe we could port this:
    % https://www.pyimagesearch.com/2017/05/22/face-alignment-with-opencv-and-python/

    %TODO: we could also use facial landmarks from opencv_contrib face module
    % (cv.Facemark and cv.FacemarkKazemi)
end

function vec = face2vec(img, net, blobOpts)
    %FACE2VEC  Get 128 floating points feature vector
    %
    % Run face recognition network to receive 128-dimensional unit feature
    % vector from input face image.
    %

    net.setInput(cv.Net.blobFromImages(img, blobOpts{:}));
    vec = net.forward();
    vec = vec(:);
end

function [name, score] = recognizeFace(vec, vecs, names, thresh)
    %RECOGNIZEFACE  Perform face recognition
    %
    % Match a new feature vector with registered ones. Return a name of the
    % best matched person.
    %
    % (NOTE: For more advanced usage, we could train an SVM classifier)
    %
    % See also: pdist2
    %

    if nargin < 4, thresh = 0.5; end
    name = 'unknown';
    score = -1;

    if ~isempty(vecs)
        scores = vec.' * vecs;  % dot-product of vec against each vecs(:,i)
        [s, idx] = max(scores);
        if s > thresh
            name = names{idx};
            score = s;
        end
    end
end

function clr = name2clr(name, names)
    clrs = round(255 * lines(7));
    idx = find(strcmp(name, names));
    if isempty(idx)
        clr = [128 128 128];
    else
        idx = rem(idx - 1, 7) + 1;
        clr = clrs(idx,:);
    end
end

function img = insertAnnotation(img, rect, str, varargin)
    % See also: insertObjectAnnotation, insertShape, insertText
    p = inputParser();
    p.addParameter('Alpha', 0.6);
    p.addParameter('Thickness', 1);
    p.addParameter('Color', [255 255 0]);
    p.addParameter('TextColor', [0 0 0]);
    p.addParameter('FontFace', 'HersheySimplex');
    p.addParameter('FontScale', 0.4);
    p.addParameter('AntiAlias', true);
    p.addParameter('Shape', 'rectangle');
    p.parse(varargin{:});
    opts = p.Results;
    opts.Shape = validatestring(opts.Shape, {'rectangle','circle'});
    thick = 1;

    [sz,b] = cv.getTextSize(str, 'Thickness',thick, ...
        'FontFace',opts.FontFace, 'FontScale',opts.FontScale);
    txt_rect = [rect(1), rect(2)-sz(2)-b, sz(1), sz(2)+b];
    txt_orig = [rect(1), rect(2)-b];

    if opts.AntiAlias
        alias = {'LineType','AA'};
    else
        alias = {'LineType',8};
    end

    overlay = img;
    if strcmp(opts.Shape, 'rectangle')
        overlay = cv.rectangle(overlay, rect, ...
            'Color',opts.Color, 'Thickness',opts.Thickness, alias{:});
    else
        c = rect(1:2) + rect(3:4)/2;
        r = max(rect(3:4)/2);
        overlay = cv.circle(overlay, c, r, ...
            'Color',opts.Color, 'Thickness',opts.Thickness, alias{:});
    end
    overlay = cv.rectangle(overlay, txt_rect, ...
        'Color',opts.Color, 'Thickness','Filled', alias{:});
    if opts.Thickness > 1
        overlay = cv.rectangle(overlay, txt_rect, ...
            'Color',opts.Color, 'Thickness',opts.Thickness, alias{:});
    end
    overlay = cv.putText(overlay, str, txt_orig, ...
        'FontFace',opts.FontFace, 'FontScale',opts.FontScale, ...
        'Color',opts.TextColor, 'Thickness',thick, alias{:});

    img = cv.addWeighted(img,1-opts.Alpha, overlay,opts.Alpha, 0);
end

%%
% Pretrained models

function dname = get_dnn_dir(dname)
    %GET_DNN_DIR  Path to model files, and show where to get them if missing

    dname = fullfile(mexopencv.root(), 'test', 'dnn', dname);
    b = isdir(dname);
    if ~b
        % display help of calling function
        % (assumed to be a local function in current file)
        st = dbstack(1);
        help([mfilename() filemarker() st(1).name])
    end
    assert(b, 'Missing model: %s', dname);
end

function [net, blobOpts] = ResNetSSD_FaceDetector()
    %RESNETSSD_FACEDETECTOR  face detector based on SSD framework with reduced ResNet-10 backbone
    %
    % homepage = https://github.com/opencv/opencv/blob/3.4.0/samples/dnn/face_detector/how_to_train_face_detector.txt
    %
    % ## Model
    %
    % file = test/dnn/ResNetSSD_FaceDetector/deploy.prototxt
    % url  = https://github.com/opencv/opencv/raw/3.4.0/samples/dnn/face_detector/deploy.prototxt
    % hash = 006BAF926232DF6F6332DEFB9C24F94BB9F3764E
    %
    % ## Weights (FP32)
    %
    % file = test/dnn/ResNetSSD_FaceDetector/res10_300x300_ssd_iter_140000.caffemodel
    % url  = https://github.com/opencv/opencv_3rdparty/raw/dnn_samples_face_detector_20170830/res10_300x300_ssd_iter_140000.caffemodel
    % hash = 15aa726b4d46d9f023526d85537db81cbc8dd566
    % size = 10.1 MB
    %
    % ## Weights (FP16)
    %
    % file = test/dnn/ResNetSSD_FaceDetector/res10_300x300_ssd_iter_140000_fp16.caffemodel
    % url  = https://github.com/opencv/opencv_3rdparty/raw/dnn_samples_face_detector_20180205_fp16/res10_300x300_ssd_iter_140000_fp16.caffemodel
    % hash = 31fc22bfdd907567a04bb45b7cfad29966caddc1
    % size = 5.1 MB
    %

    dname = get_dnn_dir('ResNetSSD_FaceDetector');
    if true
        % fp32
        net = cv.Net('Caffe', ...
            fullfile(dname, 'deploy.prototxt'), ...
            fullfile(dname, 'res10_300x300_ssd_iter_140000.caffemodel'));
    else
        % fp16: weights converted to half-precision floats (2x less model size)
        net = cv.Net('Caffe', ...
            fullfile(dname, 'deploy.prototxt'), ...
            fullfile(dname, 'res10_300x300_ssd_iter_140000_fp16.caffemodel'));
    end
    blobOpts = {'SwapRB',false, 'Crop',false, 'Size',[300 300], 'Mean',[104 117 123]};
end

function [net, blobOpts] = OpenFace_Embedding()
    %OPENFACE  OpenFace embedding for face recognition
    %
    % homepage = https://cmusatyalab.github.io/openface/
    %
    % ## Model + Weights
    %
    % file = test/dnn/OpenFace/nn4.small2.v1.t7
    % url  = https://storage.cmusatyalab.org/openface-models/nn4.small2.v1.t7
    % hash = ac8161a4376fb5a79ceec55d85bbb57ef81da9fe
    % size = 30 MB
    %

    dname = get_dnn_dir('OpenFace');
    net = cv.Net('Torch', fullfile(dname, 'nn4.small2.v1.t7'));
    blobOpts = {'SwapRB',false, 'Crop',false, 'Size',[96 96], 'ScaleFactor',1/255};
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function dice = dice_labels(lname1, lname2)
% d = dice_labels(lname1, lname2)
%
% computes Dice coefficient of the two given label files.
%

%
% dice_labels.m
%
% Original Author: Nick Schmansky
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.2 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

label1 = read_label('',lname1);
label2 = read_label('',lname2);

label1Size = size(label1,1);
label2Size = size(label2,1);

hits = 0;
for i=1:label1Size
   x1 = label1(i,2);
   y1 = label1(i,3);
   z1 = label1(i,4);
   for j=1:label2Size
     x2 = label2(j,2);
     y2 = label2(j,3);
     z2 = label2(j,4);
     if ((x1==x2) && (y1==y2) && (z1==z2))
       hits = hits+1;
       break;
     end
   end
end

dice = (2 * hits) / (label1Size + label2Size);


text:
function q_current = coppelia_get_joint_positions(coppelia)
    joint_handles = coppelia.robot.j_handles;

    q_current = [];
    for i=1:length(joint_handles)
        [error, value]=coppelia.sim.simxGetJointPosition(coppelia.clientID, joint_handles(i), coppelia.sim.simx_opmode_buffer);
        q_current = [q_current value];
    end
    q_current = q_current(:);
end

repo: ARTE


text:
classdef dm_converter_mpc2_legacy < mp.dm_converter_mpc2
%MP.DM_CONVERTER_MPC2_LEGACY  MATPOWER data model converter for MATPOWER case v2.

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

    methods
        function dm = legacy_user_mod_inputs(obj, dm, mpopt, dc)
            %% pre-process input related to user vars, constraints, costs

            %% create (read-only) copies of individual fields for convenience
            mpc = dm.source;
            [~, bus, gen, ~, ~, A, l, u, mpopt, ...
                N, fparm, H, Cw, z0, zl, zu, ~] = opf_args(mpc, mpopt);

            %% data dimensions
            nb   = size(bus, 1);    %% number of buses
            ng   = size(gen, 1);    %% number of dispatchable injections
            nlin = size(A, 1);      %% number of linear user constraints
            nw = size(N, 1);        %% number of general cost vars, w
            if dc
                %% reduce A and/or N from AC dimensions to DC dimensions, if needed
                if nlin || nw
                    acc = [nb+(1:nb) 2*nb+ng+(1:ng)];   %% Vm and Qg columns
                    if nlin && size(A, 2) >= 2*nb + 2*ng
                        %% make sure there aren't any constraints on Vm or Qg
                        if any(any(A(:, acc)))
                            error('mp.dm_converter_mpc2_legacy/legacy_user_mod_inputs: attempting to solve DC OPF with user constraints on Vm or Qg');
                        end
                        A(:, acc) = [];         %% delete Vm and Qg columns
                    end
                    if nw && size(N, 2) >= 2*nb + 2*ng
                        %% make sure there aren't any costs on Vm or Qg
                        if any(any(N(:, acc)))
                            [ii, jj] = find(N(:, acc));
                            ii = unique(ii);    %% indices of w with potential non-zero cost terms from Vm or Qg
                            if any(Cw(ii)) || (~isempty(H) && any(any(H(:, ii))))
                                error('mp.dm_converter_mpc2_legacy/legacy_user_mod_inputs: attempting to solve DC OPF with user costs on Vm or Qg');
                            end
                        end
                        N(:, acc) = [];         %% delete Vm and Qg columns
                    end
                end
                nv = 0;
                nq = 0;
            else
                nv = nb;
                nq = ng;
            end

            %% get number of user vars, check consistency
            nx = nb+nv + ng+nq;  %% number of standard OPF control variables
            if nlin
                nz = size(A, 2) - nx;   %% number of user z variables
                if nz < 0
                    error('mp.dm_converter_mpc2_legacy/legacy_user_mod_inputs: user supplied A matrix must have at least %d columns.', nx);
                end
            else
                nz = 0;               %% number of user z variables
                if nw                 %% still need to check number of columns of N
                    if size(mpc.N, 2) ~= nx;
                        error('mp.dm_converter_mpc2_legacy/legacy_user_mod_inputs: user supplied N matrix must have %d columns.', nx);
                    end
                end
            end

            %% package up parameters
            z = struct( 'nz', nz, ...           %% num user variables
                        'z0', z0, ...
                        'zl', zl, ...
                        'zu', zu    );
            lin = struct(   'nlin', nlin, ...   %% num user linear constraints
                            'A', A, ...
                            'l', l, ...
                            'u', u  );
            cost =  struct( 'nw', nw, ...       %% num user cost rows
                            'N', N, ...
                            'Cw', Cw    );
            nlc = obj.legacy_user_nln_constraints(dm, mpopt);

            if ~isempty(H)
                cost.H = H;
            end
            if ~isempty(fparm)
                cost.dd = fparm(:, 1);
                cost.rh = fparm(:, 2);
                cost.kk = fparm(:, 3);
                cost.mm = fparm(:, 4);
            end
            dm.userdata.legacy_opf_user_mods = struct( ...
                'z', z, ...
                'lin', lin, ...
                'nlc', {nlc}, ...
                'cost', cost );
        end

        function uc = legacy_user_nln_constraints(obj, dm, mpopt)
            mpc = dm.source;

            %% check for user-defined nonlinear constraints
            nnle = 0;   %% number of nonlinear user-defined equality cons
            nnli = 0;   %% number of nonlinear user-defined inequality cons
            if isfield(mpc, 'user_constraints')
                if isfield(mpc.user_constraints, 'nle')
                    for k = 1:length(mpc.user_constraints.nle)
                        nnle = nnle + mpc.user_constraints.nle{k}{2};
                    end
                end
                if isfield(mpc.user_constraints, 'nli')
                    for k = 1:length(mpc.user_constraints.nli)
                        nnli = nnli + mpc.user_constraints.nli{k}{2};
                    end
                end
            end

            %% initialize cell array for add_nln_constraint() args
            uc = cell(nnle+nnli, 1);

            %% user-defined nonlinear equalities
            if nnle
                for k = 1:length(mpc.user_constraints.nle)
                    nlc = mpc.user_constraints.nle{k};
                    fcn  = eval(['@(x)' nlc{3} '(x, nlc{6}{:})']);
                    hess = eval(['@(x, lam)' nlc{4} '(x, lam, nlc{6}{:})']);
                    uc{k} = {nlc{1:2}, 1, fcn, hess, nlc{5}};
                end
            end

            %% user-defined nonlinear inequalities
            if nnli
                for k = 1:length(mpc.user_constraints.nli)
                    nlc = mpc.user_constraints.nli{k};
                    fcn  = eval(['@(x)' nlc{3} '(x, nlc{6}{:})'])
                    hess = eval(['@(x, lam)' nlc{4} '(x, lam, nlc{6}{:})'])
                    uc{nnle+k} = {nlc{1:2}, 0, fcn, hess, nlc{5}};
                end
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% Calculate one-shot classification error rate
function run_classification

    classdir = 'model_refits';
    load('items_classification','nrun','ntrain','ntest','cell_Y');

    trainset = num2cell(vec(1:ntrain));
    testset = num2cell(vec(1:ntest));

    fprintf(1,'one-shot classification results\n');
    perror = zeros(nrun,1);
    for r=1:nrun
        Y = cell_Y{r};
        fscore = @(itrain,itest) fclassify(itrain,itest,r,classdir);
        perror(r) = myclassify(trainset,testset,fscore,Y,'score');
        fprintf(1,' run %d (error %s%%)\n',r,num2str(perror(r),3));
    end
    fprintf(1,'average error: %s%%\n',num2str(mean(perror),3));

end

%
% Bayesian classification score
%
%  log P(I_T|I_C) + log(I_C|I_T) - log(I_C)
%   for training image I_C
%   and test image I_T
%
% Input
%  itrain: train index
%  itest: test index
%  irun: run index
%  classdir: file directory for "crossFit" files
%
% Output
%  log_score: [scalar] score
%
function log_score = fclassify(itrain,itest,irun,classdir)

    srun = num2str(irun);
    strain = num2str(itrain);
    stest = num2str(itest);

    fn_test_to_train = ['run',srun,'_fit_test' ,stest,'_to_image_train',strain];
    fn_train_to_test = ['run',srun,'_fit_train',strain,'_to_image_test',stest];

    % load file optimizing P(I_T|I_C)
    if ~exist(fullfile(classdir,[fn_train_to_test,'.mat']),'file')
       fprintf(1,'Please download pre-computed model results to use this feature. Program quiting...\n');
       assert false;
    end
    load(fullfile(classdir,fn_train_to_test),'fit_score','prior_score');
    pair.fit_score = fit_score;
    pair.prior_score = prior_score;
    pair_fit_train_to_test = pair;
    clear pair

    % load file optimizing P(I_C|I_T)
    load(fullfile(classdir,fn_test_to_train),'fit_score','prior_score');
    pair.fit_score = fit_score;
    pair.prior_score = prior_score;
    pair_fit_test_to_train = pair;
    clear pair

    % compute score
    [log_P_IT_given_IC,prior_scores] = log_post_pred(pair_fit_train_to_test);
    log_P_IC_given_IT = log_post_pred(pair_fit_test_to_train);
    log_P_IC = logsumexp(prior_scores(:));
    log_score = log_P_IC_given_IT + log_P_IT_given_IC - log_P_IC;

end

%
% Log-posterior predictive score, log P(I_2 | I_1), using discrete approximation
%
% Input
%  pair: structure
%
% Output
%  logscore: log conditional probability
%  prior_scores:  [k x 1] log joint probability (unnormalized) weights for
%       each parse of I_1
%
function [logscore,prior_scores] = log_post_pred(pair)

    % extra info
    logfit = pair.fit_score;
    prior_scores = pair.prior_score;

    % normalize weights
    logwt = prior_scores - logsumexp(prior_scores(:));

    % combine weights with fit term
    logv = logfit(:) + logwt;
    logscore = logsumexp(logv);

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


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

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

% Last Modified by GUIDE v2.5 01-May-2013 21:32:35

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @MotList_OpeningFcn, ...
                   'gui_OutputFcn',  @MotList_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 MotList_figure is made visible.
function MotList_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to MotList_figure (see VARARGIN)

handles.Simuh=varargin{1};
Category=dir([handles.Simuh.MRiLabPath filesep 'Config' filesep 'Mot']);
ind=1;
for i=1:length(Category)
    if ~strcmp(Category(i).name,'.') & ~strcmp(Category(i).name,'..') & Category(i).isdir
        handles.CategoryDir{ind,1}=Category(i).name;
        ind=ind+1;
    end
end
set(handles.Category_listbox,'String',handles.CategoryDir);

handles.accept = 0;
% Choose default command line output for MotList_figure
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = MotList_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on selection change in Mot_listbox.
function Mot_listbox_Callback(hObject, eventdata, handles)
% hObject    handle to Mot_listbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns Mot_listbox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from Mot_listbox

if ~isempty(get(handles.MotPath_text,'String'))
    return;
end

if strcmp(handles.MotList,'Empty')
    return;
end

handles.MotXMLFile=[handles.Simuh.MRiLabPath filesep ...
                    'Config' filesep ...
                    'Mot' filesep ...
                    handles.CategoryDir{get(handles.Category_listbox,'Value')} filesep ...
                    handles.MotList{get(handles.Mot_listbox,'Value')} filesep...
                    handles.MotList{get(handles.Mot_listbox,'Value')} '.xml'];

handles.MotXMLDir=[handles.Simuh.MRiLabPath filesep ...
                    'Config' filesep ...
                    'Mot' filesep ...
                    handles.CategoryDir{get(handles.Category_listbox,'Value')} filesep ...
                    handles.MotList{get(handles.Mot_listbox,'Value')}];

guidata(hObject, handles);

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

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

% --- Executes on selection change in Category_listbox.
function Category_listbox_Callback(hObject, eventdata, handles)
% hObject    handle to Category_listbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns Category_listbox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from Category_listbox

if ~isempty(get(handles.MotPath_text,'String'))
    return;
end

if strcmp(handles.CategoryDir,'Empty')
    return;
end

set(handles.Mot_listbox,'Enable','on');

Mot=dir([handles.Simuh.MRiLabPath filesep ...
          'Config' filesep ...
          'Mot' filesep ...
          handles.CategoryDir{get(handles.Category_listbox,'Value')}]);
ind=1;
handles.MotList=[];
for i=1:length(Mot)
    if ~strcmp(Mot(i).name,'.') & ~strcmp(Mot(i).name,'..')
        handles.MotList{ind,1}=Mot(i).name;
        ind=ind+1;
    end
end
if ind==1
    handles.MotList='Empty';
end
set(handles.Mot_listbox,'String',handles.MotList);
set(handles.Mot_listbox,'Value',1);

guidata(hObject, handles);

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

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

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

if isempty(get(handles.MotPath_text,'String'))
    try
        if strcmp(handles.MotList{get(handles.Mot_listbox,'Value')},'Empty')
            DoUpdateInfo(handles.Simuh,'No Motion is selected!');
            return;
        end
    catch me
        DoUpdateInfo(handles.Simuh,'No Motion is selected!');
        return;
    end
end

%----update Motion selection
handles.Simuh.MotXMLFile=handles.MotXMLFile;
handles.Simuh.MotXMLDir=handles.MotXMLDir;
Mot=DoParseXML(handles.MotXMLFile);
handles.Simuh.MotStruct=Mot;
Mot_text=uicontrol(handles.Simuh.MotSel_uipanel,'Style','text');
set(Mot_text,'FontWeight','bold','FontSize',10,'Position',[0,0,140,18],...
    'String', Mot.Attributes(1).Value,'TooltipString',handles.MotXMLFile);
%----end
guidata(handles.Simuh.SimuPanel_figure,handles.Simuh);
DoUpdateInfo(handles.Simuh,'Mot file was successfully loaded!');
handles.accept = 1;
guidata(hObject, handles);
close(handles.MotList_figure);

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

% --- Executes when user attempts to close MotList_figure.
function MotList_figure_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to MotList_figure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

DoDisableButton([],[],handles.Simuh);
if handles.accept == 0
    DoUpdateInfo(handles.Simuh,'Mot file loading was cancelled!');
end
% Hint: delete(hObject) closes the figure
delete(handles.MotList_figure);

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

[filename,pathname,filterindex]=uigetfile({'Mot*.xml','XML-files (*.xml)'},'MultiSelect','off');
if filename~=0
    handles.MotXMLFile=[pathname filename];
    handles.MotXMLDir=pathname(1:end-1);
    set(handles.MotPath_text,'String',[pathname filename]);
    % Add searchig path
    path(path,handles.MotXMLDir);

else
    errordlg('No Mot was loaded!');
    return;
end
guidata(hObject, handles);

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

if isfield(handles.Simuh,'MotXMLFile')
    handles.Simuh=rmfield(handles.Simuh,'MotXMLFile');
    delete(get(handles.Simuh.MotSel_uipanel,'Children'));
    guidata(handles.Simuh.SimuPanel_figure,handles.Simuh);
end
DoUpdateInfo(handles.Simuh,'Use default Motion!');
handles.accept = 1;
guidata(hObject, handles);
close(handles.MotList_figure);

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [ FRF, FBB ] = OMP( Fopt, NRF, At )

FRF = [];
Fres = Fopt;
for k = 1:NRF
    PU = At' * Fres;
%     [aa,bb] = max(diag(PU * PU'));
    [aa,bb] = max(sum( abs(PU).^2, 2 ));
    FRF = [FRF , At(:,bb)];
    FBB = pinv(FRF) * Fopt; %use pseudoinverse to avoid the inverse of a possible singular matrix
    Fres = (Fopt - FRF * FBB) / norm(Fopt - FRF * FBB,'fro');
end

end

author: yuxianghao


text:
function [u,s,v] = svdsim(a,tol)
%SVDSIM  simple SVD program
%
% A simple program that demonstrates how to use the
% QR decomposition to perform the SVD of a matrix.
% A may be rectangular and complex.
%
% usage: [U,S,V]= SVDSIM(A)
%     or      S = SVDSIM(A)
%
% with A = U*S*V' , S>=0 , U'*U = Iu  , and V'*V = Iv
%
% The idea is to use the QR decomposition on A to gradually "pull" U out from
% the left and then use QR on A transposed to "pull" V out from the right.
% This process makes A lower triangular and then upper triangular alternately.
% Eventually, A becomes both upper and lower triangular at the same time,
% (i.e. Diagonal) with the singular values on the diagonal.
%
% Matlab's own SVD routine should always be the first choice to use,
% but this routine provides a simple "algorithmic alternative"
% depending on the users' needs.
%
%see also: SVD, EIG, QR, BIDIAG, HESS
%

% Paul Godfrey
% October 23, 2006

if ~exist('tol','var')
   tol=eps*1024;
end

%reserve space in advance
sizea=size(a);
loopmax=100*max(sizea);
loopcount=0;

% or use Bidiag(A) to initialize U, S, and V
u=eye(sizea(1));
s=a';
v=eye(sizea(2));

Err=realmax;
while Err>tol & loopcount<loopmax ;
%   log10([Err tol loopcount loopmax]); pause
    [q,s]=qr(s'); u=u*q;
    [q,s]=qr(s'); v=v*q;

% exit when we get "close"
    e=triu(s,1);
    E=norm(e(:));
    F=norm(diag(s));
    if F==0, F=1;end
    Err=E/F;
    loopcount=loopcount+1;
end
% [Err/tol loopcount/loopmax]

%fix the signs in S
ss=diag(s);
s=zeros(sizea);
for n=1:length(ss)
    ssn=ss(n);
    s(n,n)=abs(ssn);
    if ssn<0
       u(:,n)=-u(:,n);
    end
end

if nargout<=1
   u=diag(s);
end

return


text:
% Copyright (c) Technische Universität München, 2012
% All rights reserved.
%
% Author: Clemens Hage
% Contact: hage@tum.de
% Date: 2012-11-12
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% This script runs a quick numerical demo on Robust PCA (batch mode with
% incomplete observation)

% Dimension of the data set
m = 400;

% Number of samples
n = 400;

% Dimension of the underlying subspace
k = 40; % i.e. relative rank of 0.1

% Relative density of outliers
rho = 0.1; % i.e. 10% of the matrix entries are corrupted

% Max. magnitude of the outliers
amplitude = 5;

% Create test data
[X,L_0]=testdata(m,n,k,rho,amplitude,0);

support=0.5; % 50% of the entries are observed

% Create support set

% number of observed entries
nentries = ceil(support * m * n);

% random positions
indices = randperm(m*n,nentries);

% convert to indices
[i,j]=ind2sub([m n],indices);

% create sparse matrix
Omega = sparse(i,j,true(nentries,1),m,n,nentries);

% run Robust PCA
L=robustpca_batch(X,k,'atansquare','L_0',L_0,'Omega',Omega);

% Compute relative subspace reconstruction error
err = norm(Omega.*L-Omega.*L_0,'fro') / norm(Omega.*L_0,'fro');
display(['Relative subspace reconstruction error (robust method): ' num2str(err)])

% % Compare to common PCA approach
X_sampled = Omega.*X;
[U,~,~]=svds(X_sampled,k);

L_PCA = U * (U' * X_sampled);

err_PCA = norm(Omega.*L_PCA-Omega.*L_0,'fro') / norm(Omega.*L_0,'fro');
display(['Relative subspace reconstruction error (common PCA): ' num2str(err_PCA)])

author: andrewssobral

repo: lrslibrary


text:
function demoMertonModel()
%%DEMOMERTONMODEL Demonstrates scenario simulation of two independent
%             Merton processes using the strongStochTaylorStepJump
%             function. A plot comparing the approximated paths to the true
%             path is generated. In order to demonstrate ideal results, the
%             same collection of driving processes are used for both the
%             true path and the simulations. The first line may be
%             uncommented and modified if reproducible plots are desired.
%
%The Merton model is discussed in detail in Chapter 1.7, pages 50-52 and
%Chapter 9.6, page 414 of [1]. See equations 9.6.3 and 9.6.4 and the
%discussion on page 414. Note we simplified the jump process term's
%coefficient to c(y)=psi*y as suggested in Chapter 7.1, page 313.
%
%REFERENCES:
%[1] Platen, Eckhard, and Nicola Bruti-Liberati. Numerical solution of
%    stochastic differential equations with jumps in finance. Vol. 64.
%    Springer Science & Business Media, 2010.
%
%July 2019 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    %rng(16,'twister'); % Make the example reproducible
    close all;

    %Define Parameters
    X0=[1;1]; % Initial path value
    mue=[0.05;0.05]; % Drift mean
    sig=[0.2;0.2]; % Wiener process standard deviations
    psi=[-0.25;-0.25]; % Jump intensity coefficients
    T=10; % End time starting from 0
    lambda=0.3; % Poisson jump frequency for unit time interval

    %Define Discretization Coefficients and Derivatives
    % Assume a time invariant SDE of the form:
    % dY_{n+1} = a(Y_n) dt + b(Y_n) dW + c(Y_n) dP
    a = @(x) diag(mue)*x;
    B = @(x) diag(sig)*diag(x);
    c = @(x) diag(psi)*x;
    pBpy = cat(3,sig(1)*[1 0; 0 0],sig(2)*[0 0; 0 1]);
    pcpy = diag(psi);

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

    %%%%%%%%%%%%%%%%%%%% Exact Solution %%%%%%%%%%%%%%%%%%%%%
    %The Merton model is discussed in detail on p.50-52 and p.414 of [1].
    % We simplified the compound jump process to have constant jump size.
    % See equations 9.6.3 and 9.6.4.

    dt = 0.001; % True path time step
    t = 0:dt:T;
    X = zeros([2,length(t)]); % A vector to contain the true path values
    X(:,1) = X0; % Initial condition

    %Generate Poisson samples as 2 independent compound Poisson processes
    CP = [0,genMarkedPointProc(lambda,dt,0+dt,T,0)];

    %Simulate 2 independent Wiener Processes
    W = [[0;0],sqrt(dt).*randn([2,length(t)-1])];
    W = cumsum(W,2);

    %Compute Merton Process Values using Exact Solution from [1]
    for i = 2:length(X)
        X(:,i) = X0.*exp((mue-0.5*sig.^2)*t(i)+sig.*W(:,i)+psi*CP(i));
    end

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

    %%%%%%%%%%%%%%%%%%%% Approximations %%%%%%%%%%%%%%%%%%%%%

    %Define Approximation Parameters
    dtapprox = 0.5; % Approximation time step.
    tapprox = 0:dtapprox:T;
    numruns = 50; % Number of times to perform the random simulation steps.

    XEuler = zeros([2,length(tapprox)]);
    XEuler(:,1) = X0;
    XTaylor = zeros([2,length(tapprox)]);
    XTaylor(:,1) = X0;
    XDerFree = zeros([2,length(tapprox)]);
    XDerFree(:,1) = X0;

    for i = 2:length(tapprox)
        [deltaP,Wj] = getDrivers(W,CP,dtapprox,dt,i);

        %Get Euler Approximation
        for run = 1:numruns
            XEuler(:,i) = XEuler(:,i) + strongStochTaylorStepJump(...
                XEuler(:,i-1),a,B,@(x,k,idx)c(x),lambda,dtapprox,1,[],[],deltaP,Wj);
        end
        XEuler(:,i) = XEuler(:,i)/numruns;

        %Get Taylor Approximation
        for run = 1:numruns
            XTaylor(:,i) = XTaylor(:,i) + strongStochTaylorStepJump(...
                XTaylor(:,i-1),a,B,c,lambda,dtapprox,3,pBpy,pcpy,deltaP,Wj,tapprox(i));
        end
        XTaylor(:,i) = XTaylor(:,i)/numruns;

        %Get Derivative Free Taylor Approximation
        for run = 1:numruns
            XDerFree(:,i) = XDerFree(:,i) + strongRungeKStepJump(...
                XDerFree(:,i-1),a,B,c,lambda,dtapprox,3,deltaP,Wj,tapprox(i));
        end
        XDerFree(:,i) = XDerFree(:,i)/numruns;
    end

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

    %%%%%%%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%%

    %Plot the Paths
    figure('Position',[10 10 900 600])
    plot3(t,X(1,:),X(2,:)); label11='Exact';
    hold on
    plot3(tapprox,XEuler(1,:),XEuler(2,:)); label12='Euler';
    plot3(tapprox,XTaylor(1,:),XTaylor(2,:)); label13='Taylor';
    plot3(tapprox,XDerFree(1,:),XDerFree(2,:)); label14='DF';
    legend(label11,label12,label13,label14)
    xlabel('Time (t)','FontSize',16)
    ylabel('Position ($X_t^{(1)}$ or $\tilde{X}_t^{(1)})$','Interpreter','latex','FontSize',16)
    zlabel('Position ($X_t^{(2)}$ or $\tilde{X}_t^{(2)})$','Interpreter','latex','FontSize',16)
    title('A Simulated Merton Model','FontSize',16)
    hold off

    %Plot the Errors
    figure('Position',[10 10 900 600])
    errorEuler = sum((X(:,1:floor(dtapprox/dt):end)-XEuler).^2)/2;
    errorTaylor = sum((X(:,1:floor(dtapprox/dt):end)-XTaylor).^2)/2;
    errorDerFree = sum((X(:,1:floor(dtapprox/dt):end)-XDerFree).^2)/2;
    plot(tapprox,errorEuler); label21='Euler';
    hold on
    plot(tapprox,errorTaylor); label22='Taylor';
    plot(tapprox,errorDerFree); label23='DF';
    legend(label21,label22,label23)
    xlabel('Time (t)','FontSize',16)
    ylabel('$\frac{1}{2}\left[(X_t^{(1)}-\tilde{X}_t^{(1)})^2+(X_t^{(2)}-\tilde{X}_t^{(2)})^2\right]$','Interpreter','latex','FontSize',16)
    title('Mean Squared Error for Each Time Step','FontSize',16)
    hold off
end

function [deltaP,Wj]=getDrivers(W,CP,dtapprox,dt,iter)
%%GETDRIVERS A helper function to automate retrieval of the correct
%            interval of values from the driving processes used for the
%            generation of the true path.
%
%INPUTS: W An mXlength(t) vector containing realized Wiener processes in
%          each row.
%       CP A 1Xlength(t) row vector containing a realized Poisson process.
% dtapprox The coarse approximation time step (assumed to be constant).
%       dt The fine true path time step (assumed to be constant).
%     iter The current iteration index (assumed indices are positive and
%          begin at 1).
%
%OUTPUTS: deltaP The number of jumps which occurred in CP during the time
%                interval.
%             Wj An mXdeltaP+2 vector which contains the values from W
%                which most closely occured at the same time as each jump
%                counted in deltaP. This assumes the distribution of jumps
%                over the time interval are uniformly random.
%
%July 2019 Codie T. Lewis, Naval Research Laboratory, Washington D.C.

    if(iter>1)
        currentIdx=floor(dtapprox/dt)*(iter-1)+1;
        prevIdx=floor(dtapprox/dt)*(iter-2)+1;
        deltaP=CP((iter-1)*floor(dtapprox/dt)+1)-CP((iter-2)*floor(dtapprox/dt)+1);
        if(deltaP~=0)
            jumpIdx=floor((prevIdx+rand([deltaP,1])*dtapprox/dt));
            Wj=W(:,[prevIdx;jumpIdx;currentIdx]);
        else
            Wj=W(:,[prevIdx;currentIdx]);
        end
    else
        error('Iter must be >1.')
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary