text:
%kvlabel 'Perform Labeling on Multiband or Cluster Image (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vlabel.pane file
%
% Parameters:
% Integer: w 'Border Width ', default: 0: 'Specify the border width in pixels'
% MultiChoice: m 'Type of Distance Metric (Choose One): ', default: 1: 'Distance Metric: 1 = Euclidean distance, 2 = city block'
%    Choices are:
%   1: 'Euclidean '
%   2: 'City Block'
% MultiChoice: c 'Type of Connectivity (Choose One): ', default: 1: 'Connectivity: 1 = 4-connectivity, 2 = 8-connectivity'
%    Choices are:
%   1: '4 Connectivity'
%   2: '8 Connectivity'
% Double: s 'Minimum Region Size in % ', default: 1: 'Minimum number of pixel to keep in label region'
% OutputFile: o 'Output Labeled Image', required: 'output image after the labeling process'
% OutputFile: asc 'Output ASCII File', optional: 'Output ASCII file for vlabel result display'
% InputFile: i1 'Initial Image ', optional: 'input image'
% InputFile: i2 'Cluster Center Image', optional: 'input image for cluster centers'
% InputFile: i3 'Cluster Number Image', optional: 'input image for cluster numbers'
% Integer: n 'Approximate Number of Regions Expected ', default: 10: 'Number of region to get from labeling process'
% Double: f 'Split and Merge Factor', default: 0.07: 'Split & Merge Factor: higher values result in finer regions'
%
% Example: [o, asc] = kvlabel({i1, i2, i3}, {'w',0;'m',1;'c',1;'s',1;'o','';'asc','';'i1','';'i2','';'i3','';'n',10;'f',0.07})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vlabel - Perform Labeling on Multiband or Cluster Image  (K1)
%
%  DESCRIPTION
% .I vlabel
% performs a labeling on a multiband image or a cluster image by
% attempting to merge connected pixels.
%
% The principal of the algorithm is as follows: A pixel receives the same
% label as its neighbor if the likelihood distance between the two pixels
% is acceptable.  The label process is propagated for a given region number
% until it is no longer possible to find a candidate.
%
% \fBThree different types of labeling choices exist":
%
%  1
% .I First choice:
% uses a single or multi band image, where the data storage can be
% any of the following types: VFF_TYP_1_BYTE, VFF_TYP_2_BYTE,
% VFF_TYP_4_BYTE, or VFF_TYP_FLOAT.
% The input image corresponds to the -i1 argument.
% The distance is computed using all the bands of the image.
%
%  2.
% .I Second choice:
% uses a cluster number and cluster center
% image obtained from an algorithm like vkmeans or vquant.
% For this case, the cluster center represents the value of a
% class of pixels that have been grouped together.
% Therefore, the distance between two neighbors will be the
% distance between their clusters.
% This case will require less computation time because the algorithm
% will only compute the inter-class distance, instead of computing
% the distance of two neighbors for the entire image.
% An additional advantage with this choice, is that the output from
% algorithms such as vkmeans or vquant may be utilized, which may lead
% to better results.
% The cluster center image corresponds to -i2 argument
% The cluster number image corresponds to -i3 argument in the
% command line.
%
%  3.
% .I The final possibility
% is to use a single or multiband input image
% (argument -i1) associated with a cluster number image (argument -i3).
% The advantage of this choice is that the results of a clustering
% algorithm are used to keep the neighbor pixels that have the same cluster
% number in the same class, and to rely on the distance in the single or
% multi band image to group two neighbors that do not belong to the same
% cluster.
%
% In summary, the three possible choices are:
%
%
% -i1 image.xv
%
%
% -i2 cluster_center.xv -i3 cluster_number.xv
%
%
% -i1 image.xv -i3 cluster_number.xv
%
% .I The algorithm also requires
% the following parameters:
% \fBMetric distance:"
% There are 2 different metric distances that can be used.
%
%
% -d 1 uses  Euclidean distance: sqrt[(x-s)^2 + (y-t)^2].
%
%
% -d 2 uses  City Block distance: |x-s| + |y-t|.
%
% \fBConnectivity:"
% There are two possible neighborhoods:
%
%
% -c 1  uses the 4 connectivity to link pixels together.
%
%
% -c 2  uses the 8 connectivity.
%
% \fBMinimum size of a region:"
%
%
% -s  (float_value)  determines the number of pixels required for a region
% to be retained.
% The minimum number of pixels is equal to:
% Total number of pixels in image * float_value / 100.0 (-s corresponds to
% a percentage of the total number of pixels in the image).
%
% \fBBorder Size:"
% Each pixel in the image is updated except those outside of the border.
% The size of the border is specified by the -w argument.
%
% \fBMerging Process:"
% When the labeling process is computed, the user can expect that the small
% rejected region will be merged together in a bigger acceptable region or
% will be included inside another connected region.
% This choice is selected by setting the logical argument -merge to TRUE.
% If -merge is set to FALSE, the small regions will be ignored and labeled as
% an UNDEFINED REGION (label number 0), the same as the border.
%
% "The AUTOMATIC or MANUAL OPTION:\fP
% This option allows the user either to fix a threshold, or to give
% an approximate number of regions.
% If the AUTOMATIC option is used, the algorithm will iterate on the threshold
% until the number of regions labeled by the process is comparable to the number
% of expected regions.
% In fact, if the expected number is not reached after 30 iterations,
% the threshold that gives the closest number of regions is used for
% the final labeling.
%
% Although this option is easy to use, the function:
%
%
% number of regions = F(Threshold)  is not a monotonically increasing function,
% and the convergence toward a solution may not exist.
%
%
% -n  int_value   (AUTOMATIC OPTION) determines the approximate final number of
% regions expected.
% -f  float_value (MANUAL OPTION) determines the threshold used by the labeling
% process.
% THESE TWO OPTIONS ARE MUTUALLY EXCLUSIVE
% The default threshold value is 0.07. This value generally gives good results
% on noiseless images with large uniform regions.
% Decreasing this value will increase the number of regions found during the
% labeling process, but these regions will get smaller and could be rejected
% by the minimum size threshold.
% Increasing this value will decrease the number of regions found during the
% labeling process. At the same time, the number of small regions will decrease
% which means that this area of the curve, \fBnumber of regions = F(Threshold)",
% is more stable than the other one.
% Once the user becomes accustomed to this routine, good results are generally
% obtained.  One way to become familiar with the routine, is to use the
% automatic option and analyze the output ASCII file (Statistics on the
% iteration process). This file contains the number of regions labeled for
% each iteration, allowing the user to see how the number of regions changes
% as the THRESHOLD is changed.
%
% \fBOutput Files"
%
%  1
% \fBOutput image:"
% The resulting image, which corresponds to the -o1 argument, contains the
% labeled image in which every pixel has a region number as its value. This
% image is of data storage type VFF_TYP_4_BYTE.
% The region label numbers are 1 to N.  The region number 0 is reserved
% as an UNDEFINED label or for the border.
%
%  2
% \fBOutput Statistic ASCII File:"
% This file contains all the information relative to the labeling process.
%
%
%
%
%  EXAMPLES
%
% vlabel -i1 image.xv -d 1 c 1 -merge 0 -f 0.07 -s 0.7 -w 2 -o image1.xv -asc stats
% This command will label image.xv using the Euclidean distance, the
% 4 connectivity, a split and merge factor equal to 0.07, a minimum
% size for the regions equal to 0.7 percent of the total number of pixels
% in image.xv, and a border size of 2. The merge option is not used which
% means that the small regions will not be labeled.
% The labeled image will be stored in image1.xv and the statistics will be
% written in the ASCII file, stats.
%
%  "SEE ALSO"
% lvkmeans(3), lvquant(3), vkmeans(1), vquant(1).
%
%  RESTRICTIONS
% vlabel works only with cluster number and cluster center images conforming
% to the convention established by the set of clustering algorithms.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kvlabel(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvlabel(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'w', 0;'m', 1;'c', 1;'s', 1;'o', '__output';'asc', '__output';'i1', '__input';'i2', '__input';'i3', '__input';'n', 10;'f', 0.07};
maxval={16,0,0,100,0,1,1,1,1,2500,1};
minval={0,0,0,0,0,1,1,1,1,2,0};
istoggle=[1,0,0,1,0,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'Integer','MultiChoice','MultiChoice','Double','OutputFile','OutputFile','InputFile','InputFile','InputFile','Integer','Double'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vlabel"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function arg = reverse(arg)

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

arg = arg([end:-1:1]);

repo: REST


text:
function out = rdivide(T,S)

if isa(S,'double')

  if numel(S)>1
    S = repmat(reshape(S,[ones(1,T.rank) size(T)]), [3*ones(1,T.rank),1,1]);
  end

  T.M = T.M ./ S;
  out = T;

elseif isa(S,'tensor')

  T.M = T.M ./ S.M;
  out = T;

end

author: mtex-toolbox

repo: mtex

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


text:
function alpha = p05_alpha ( )

%*****************************************************************************80
%
%% P05_ALPHA returns ALPHA for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real ALPHA, the value of ALPHA.
%
  alpha = 2.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% Usage: [W [optim]]=mexFistaGraph(Y,X,W0,graph,param);
%
% Name: mexFistaGraph
%
% Description: mexFistaGraph solves sparse regularized problems.
%         X is a design matrix of size m x p
%         X=[x^1,...,x^n]', where the x_i's are the rows of X
%         Y=[y^1,...,y^n] is a matrix of size m x n
%         It implements the algorithms FISTA, ISTA and subgradient descent.
%
%         It implements the algorithms FISTA, ISTA and subgradient descent for solving
%
%           min_W  loss(W) + lambda psi(W)
%
%         The function psi are those used by mexProximalGraph (see documentation)
%         for the loss functions, see the documentation of mexFistaFlat
%
%         This function can also handle intercepts (last row of W is not regularized),
%         and/or non-negativity constraints on W.
%
% Inputs: Y:  double dense m x n matrix
%         X:  double dense or sparse m x p matrix
%         W0:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)
%              initial guess
%         graph: struct (see documentation of mexProximalGraph)
%         param: struct
%            param.loss (choice of loss, see above)
%            param.regul (choice of regularization, see function mexProximalFlat)
%            param.lambda (regularization parameter)
%            param.lambda2 (optional, regularization parameter, 0 by default)
%            param.lambda3 (optional, regularization parameter, 0 by default)
%            param.verbose (optional, verbosity level, false by default)
%            param.pos (optional, adds positivity constraints on the
%                coefficients, false by default)
%            param.numThreads (optional, number of threads for exploiting
%                multi-core / multi-cpus. By default, it takes the value -1,
%                which automatically selects all the available CPUs/cores).
%            param.max_it (optional, maximum number of iterations, 100 by default)
%            param.it0 (optional, frequency for computing duality gap, every 10 iterations by default)
%            param.tol (optional, tolerance for stopping criteration, which is a relative duality gap
%                if it is available, or a relative change of parameters).
%            param.gamma (optional, multiplier for increasing the parameter L in fista, 1.5 by default)
%            param.L0 (optional, initial parameter L in fista, 0.1 by default, should be small enough)
%            param.fixed_step (deactive the line search for L in fista and use param.L0 instead)
%            param.compute_gram (optional, pre-compute X^TX, false by default).
%            param.intercept (optional, do not regularize last row of W, false by default).
%            param.ista (optional, use ista instead of fista, false by default).
%            param.subgrad (optional, if not param.ista, use subradient descent instead of fista, false by default).
%            param.a, param.b (optional, if param.subgrad, the gradient step is a/(t+b)
%            also similar options as mexProximalTree
%
%            the function also implements the ADMM algorithm via an option param.admm=true. It is not documented
%            and you need to look at the source code to use it.
%
%
% Output:  W:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)
%          optim: optional, double dense 4 x n matrix.
%              first row: values of the objective functions.
%              third row: values of the relative duality gap (if available)
%              fourth row: number of iterations
%
% Author: Julien Mairal, 2010

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
%{
    Demonstration of fault detection and fault diagnosis using KPCA.
%}
clc
clear all
close all
addpath(genpath(pwd))

load('.\data\TE.mat', 'trainData', 'testData')
kernel = Kernel('type', 'gaussian', 'gamma', 1/128^2);

parameter = struct('numComponents', 0.65, ...
                   'kernelFunc', kernel,...
                   'diagnosis', [300, 500]);

% build a KPCA object
kpca = KernelPCA(parameter);
% train KPCA model
kpca.train(trainData);
% test KPCA model
results = kpca.test(testData);

% Visualization
kplot = KernelPCAVisualization();
kplot.cumContribution(kpca)
kplot.trainResults(kpca)
kplot.testResults(kpca, results)
kplot.diagnosis(results)

author: iqiukp

repo: KPCA-MATLAB

save_path: github-repos/MATLAB/iqiukp-KPCA-MATLAB


text:
function Data=ImgRead(imgpath,suffixname,nline,bkformat,machinefmt)
%Read data stack from a specified directory.
%   Usage:
%       Data=ImgRead(imgpath,suffixname,nline,bkformat,machinefmt);
%
%
%   Inputs:
%   - imgpath:      The path of image set
%   - suffixname:   The suffix of all files in image set, e.g., mli
%   - nline:        The image height
%   - bkformat:     See freadbkj.m for details
%   - machinefmt:   See freadbkj.m for details
%
%   Outputs:
%   - Data.datastack:   A height by width by page matrix where each page
%                       corresponds to a 2D image
%   - Data.filename:    The file name list.
%
%   Examples:
%   To read a batch of intensity series in float32 format with Big-endian
%   ordering machinefmt, which has a height of 200 lines,use:
%   Data=ImgRead('/home/user/INSAR/COHEST/MLI','mli',200,'float32');
%
%   For complex differential interferogram with height 1800, use:
%   Data=ImgRead('/home/user/INSAR/COHEST/DIFF','diff',1800,'cpxfloat32');
%
%
%   Mi JIANG, Hohai University/The Hong Kong Polytechnic University,

if nargin < 5
    machinefmt='b'; % b - GAMMA software, for example
end

if nargin < 4
    bkformat='float32'; %for *mli,*cc file
end

if nargin < 3
    help ImgRead
    return;
end

if isempty(strmatch(imgpath(end),filesep))
    imgpath=[imgpath,filesep];
end

tag_files = dir([imgpath,'*',suffixname]);
img_num = length(tag_files);
disp(['The number of the ', suffixname,' images:' num2str(img_num)]);

for ii=1:img_num
    tic;
    Data.datastack(:,:,ii)=single(freadbkj([imgpath,tag_files(ii).name],nline,bkformat,machinefmt));
    temp=regexp(tag_files(ii).name,'\d+','match');
    if length(temp)==1     %mli
        Data.filename(ii,1)=str2double(temp{1});
    elseif length(temp)==2 %intf
        Data.filename(ii,1)=str2double(temp{1});
        Data.filename(ii,2)=str2double(temp{2});
    else
        error('The format of file name should be: <yyyymmdd> or <yyyymmdd_yyyymmdd>.')
    end
    time=toc;
    fprintf('Reading Img %d / %d, time = %.0f sec\n',ii,img_num,time);
end

author: DinhHoTongMinh

repo: TomoSAR


text:
% help2html() - Convert a Matlab m-file help-message header into an .html help file
%
% Usage:
%  >> linktext = help2html( filein, fileout, 'key1', val1, 'key2', val2 ...);
%
% Inputs:
%   filein       - input filename  (with .m extension)
%   fileout      - output filename (if empty, generated automatically)
%
% Optional inputs:
%   'header'     - command to insert in the header (i.e. javascript
%                   declaration or meta-tags). {default: none}.
%   'footer'     - command to insert at the end of the .html file (e.g.,
%                   back button). {default: none}.
%   'refcall'    - syntax to call references. {default: '%s.html'} For
%                  javascript function uses 'javascript:funcname(''%s.js'')'
%   'font'       - font name
%   'background' - background tag (i.e. '').
%                   {default: none}.
%   'outputlink' - html command text to link to the generated .html page.
%                  Must contain two '%s' symbols to the function title
%                  and to the function link.
%                  Ex: ... href="%s.html" ... {default: standard .html href}.
%   'outputtext' - Text used in the outputlink. {default: the function
%                  name}
%   'outputonly' - ['on'|'off'] Only generate the linktext {default: 'off'}
%
% Output:
%   fileout      - .html file written to disk
%   linktext      - html-text link to the output file
%
% M-file format:
%   The following lines describe the header format of an m-file function
%   to allow .html help file generation. Characters '-' and ':' are used
%   explicitly by the function for parsing.
%
%%    function_name() - description line 1
%%                      description line 2
%%                      etc.
%%
%%    Title1:
%%      descriptor1  - text line 1
%%                     text line 2
%%      descriptor2  - [type] text line 1
%%      "descriptor 3" - text line 1 (see notes)
%%                       etc.
%%
%%    Title2:
%%      text line 1 [...](see notes)
%%      text line 2
%%
%%    See also:
%%     function1(), function2()
%
% Author:  Arnaud Delorme, Salk Institute 2001
%
% Notes: 1) The text lines below Title2 are considered as is (i.e.,
%           preserving Matlab carriage returns) unless there is a
%           Matlab continuation cue ('...'). In this case, lines are
%           contcatenated. As below 'Title1', all text lines following
%           each  descriptor (i.e., single_word followed by '-' or '='
%           or multiple quoted words followed by a '-' or '=')
%           are concatenated.
%        2) The pattern 'function()' is detected and is printed in bold
%           if it is the first function descriptor. Otherwise,
%           it is used as a web link to the .html function file
%           'function.html' if this exists.
%        3) If a 'function.jpg' image file (with same 'function' name) exists,
%           the image is inserted into the function .html file following
%           the function description. If the .jpg file is absent, the function
%           checks for the presence of a .gif file.
%        4) Lines beginning by '%%' are not interpreted and will be printed as is.
%        5) if [type] is present in a "descriptor2  - [type] text line 1"
%           type is bolded.

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

function [linktext,allvars,alltext] = help2html( filename, htmlfile, varargin)

if nargin < 1
	help help2html;
	return;
end;
if nargin <3
   g = [];
else
   g = struct( varargin{:});;
end;

try, g.font; 			catch, g.font		= 'Helvetica'; 	end;

g.functionname = [ '%s' ];
g.description  = [ '%s' ];
g.title        = [ '%s' ];
g.var          = [ '
%s
' ];
g.vartext      = [ '%s' ];
g.text         = [ '%s' ];

g.normrow      = '
';
g.normcol1     = '';
g.normcol2     = '';
g.doublecol    = '';
g.seefile      = [ 'See the matlab file %s (may require other functions)

' ];

try, g.outputonly; 		catch, g.outputonly	= 'off'; end;
try, g.background; 		catch, g.background	= ''; 	end;
try, g.header; 			catch, g.header		= ''; 	end;
try, g.footer; 			catch, g.footer		= ''; 	end;
try, g.refcall; 		catch, g.refcall	= '%s.html'; 	end;
try, g.outputlink; 		catch, g.outputlink = [ g.normrow g.normcol1 '%s.html	' g.normcol2 '%s
' ]; end;

g.footer      = [ '' g.footer '' ];

if nargin < 1
	help help2html;
	return;
end;

% output file
% -----------
if nargin < 2 | isempty(htmlfile);
	indexdot = findstr( filename, '.');
	if isempty(indexdot), indexdot = length(filename)+1; end;
	htmlfile = [ filename(1:indexdot(end)-1) '.html' ];
else
	indexdot = findstr( filename, '.');
end;

% open files
% ----------
fid = fopen( filename, 'r');
if fid == -1
	error('Input file not found');
end;
if ~strcmp(g.outputonly, 'on')
	fo = fopen(htmlfile, 'w');
	if fo == -1
		error('Cannot open output file');
	end;
	% write header
	% ------------
	fprintf(fo, '%s\n%s\n
\n', g.header, g.background);
end;

cont = 1;
% scan file
% -----------
str = fgets( fid );
% if first line is the name of the function, reload
if str(1) ~= '%', str = fgets( fid ); end;

% state variables
% -----------
maindescription 	 = 1;
varname = [];
oldvarname = [];
vartext = [];
newvar  	 = 0;
allvars = {};
alltext = {};
indexout = 1;

while (str(1) == '%')
   str = deblank(str(2:end-1));

   % --- DECODING INPUT
   newvar  	 = 0;
   str = deblank2(str);

   if ~isempty(str)

	  % find a title
	  % ------------
	  i2d = findstr(str, ':');
	  newtitle = 0;
	  if ~isempty(i2d)
	  	 i2d = i2d(1);
	  	 switch lower(str(1:i2d))
	  		case { 'usage:' 'authors:' 'author:' 'notes:' 'note:' 'input:' ...
	  		'inputs:' 'outputs:' 'output:' 'example:' 'examples:' 'see also:' }, newtitle = 1;
		 end;
		 if (i2d == length(str)) & (str(1) ~= '%'), newtitle = 1; end;
   	  end;
      if newtitle
  			tilehtml = str(1:i2d);
  			newtitle = 1;
	        oldvarname = varname;
			oldvartext = vartext;
        	if i2d < length(str)
         			vartext = formatstr(str(i2d+1:end), g.refcall);
         	else	vartext = [];
         	end;
         	varname = [];
      else
	  % not a title
	  % ------------
         % scan lines
         [tok1 strrm] = mystrtok( str );
         [tok2 strrm] = strtok( strrm );

         if ~isempty(tok2) & ( tok2 == '-' | tok2 == '=') % new variable
            newvar = 1;
            oldvarname = varname;
            oldvartext = vartext;
			if ~maindescription
				varname = formatstr( tok1, g.refcall);
            else
				varname = tok1;
			end;
			strrm = deblank(strrm);            % remove tail blanks
            strrm = deblank(strrm(end:-1:1));	% remove initial blanks
           	strrm = formatstr( strrm(end:-1:1), g.refcall);
            vartext = strrm;
         else
            % continue current text
            str = deblank(str);            % remove tail blanks
            str = deblank(str(end:-1:1));	% remove initial blanks
            str = formatstr( str(end:-1:1), g.refcall );
            if isempty(vartext)
               vartext = str;
            else
               if ~isempty(varname)
               	    vartext = [ vartext ' ' str]; % espace if in array
               else
               		if all(vartext(	end-2:end) == '.')
               			vartext = [ deblank2(vartext(1:end-3)) ' ' str]; % espace if '...'
               		else
                    	vartext = [ vartext '
' str];    % CR otherwise
                    end;
               end;
            end;
         end;
         newtitle = 0;
      end;
	  % --- END OF DECODING

   	  str = fgets( fid );

      % test if last entry
      % ------------------
      if str(1) ~= '%'
  	 	 if ~newtitle
  	 	 	if ~isempty(oldvarname)
				allvars{indexout} = oldvarname;
				alltext{indexout} = oldvartext;
				indexout = indexout + 1;
           		fprintf( fo, [ '
' g.normrow g.normcol1 g.var '	\n' ], oldvarname);
           		fprintf( fo, [ g.normcol2 g.vartext '
\n' ], finalformat(oldvartext));
       	 	else
       			if ~isempty(oldvartext)
           			fprintf( fo, [ g.normcol2 g.text '
\n' ], finalformat(oldvartext));
       			end;
         	end;
         	newvar = 1;
         	oldvarname = varname;
         	oldvartext = vartext;
         end;
      end;

      % test if new input for an array
      % ------------------------------
      if newvar | newtitle
         if maindescription
            if ~isempty(oldvartext) % FUNCTION TITLE
               maintext = oldvartext;

			   % generate the output command
			   % ---------------------------
			   try, g.outputtext; 		catch, g.outputtext	= ''; 	end;
			   if isempty(g.outputtext),  g.outputtext	=  filename(1:indexdot(end)-1); end;
			   linktext = sprintf( g.outputlink, g.outputtext,  filename(1:indexdot(end)-1), maintext );
			   if strcmp(g.outputonly, 'on')
                   fclose(fid);
				   return;
			   end;

               maindescription = 0;
               functioname = oldvarname( 1:findstr( oldvarname, '()' )-1);
               fprintf( fo, [g.normrow g.normcol1 g.functionname '	\n'],upper(functioname));
               fprintf( fo, [g.normcol2 g.description '
\n'], [ upper(oldvartext(1)) oldvartext(2:end) ]);

               % INSERT IMAGE IF PRESENT
               imagename = [];
               imagename1 = [ htmlfile( 1:findstr( htmlfile, functioname )-1) functioname '.jpg' ];
               imagename2 = [ htmlfile( 1:findstr( htmlfile, functioname )-1) functioname '.gif' ];
			   if exist( imagename2 ), imagename = imagename2; end;
			   if exist( imagename1 ), imagename = imagename1; end;
               if ~isempty(imagename) % do not make link if the file does not exist
                   imageinfo = imfinfo(imagename);
                   if imageinfo.Width < 600
                       fprintf(fo, [ g.normrow g.doublecol ...
                                     '

' ]);
                   else
                       fprintf(fo, [ g.normrow g.doublecol ...
                                     '

' ]);
                   end;

               end;
            end;
   		elseif ~isempty(oldvarname)
			allvars{indexout} = oldvarname;
			alltext{indexout} = oldvartext;
			indexout = indexout + 1;
       		fprintf( fo, [ '
' g.normrow g.normcol1 g.var '	\n' ], oldvarname);
       		fprintf( fo, [ g.normcol2 g.vartext '
\n' ], finalformat(oldvartext));
   	 	else
   			if ~isempty(oldvartext)
       			fprintf( fo, [ g.normcol2 g.text '
\n' ], finalformat(oldvartext));
   			end;
         end;
      end;

      % print title
      % -----------
      if newtitle
		 fprintf( fo, [ g.normrow g.doublecol '

' g.normrow g.normcol1 g.title '	\n' ], tilehtml);
      	 if str(1) ~= '%' % last input
         	fprintf( fo, [ g.normcol2 g.text '
\n' ], vartext);
         end;
         oldvarname = [];
         oldvartext = [];
      end;
   else
      str = fgets( fid );
   end;
end;
fprintf( fo, [ '
\n
' g.seefile '%s'], lower(filename), filename, g.footer);
fclose( fid );
fclose( fo  );

return;

% -----------------
% sub-functions
% -----------------
function str = deblank2( str ); % deblank two ways
   str = deblank(str(end:-1:1));    % remove initial blanks
   str = deblank(str(end:-1:1));	% remove tail blanks
return;

function strout = formatstr( str, refcall );
	[tok1 strrm] = strtok( str );
		strout = [];
		while ~isempty(tok1)
			tokout = functionformat( tok1, refcall );
			if isempty( strout)
				strout = tokout;
			else
				strout = [strout ' ' tokout ];
			end;
			[tok1 strrm] = strtok( strrm );
		end;
return;

% final formating
function str = finalformat(str); % bold text in bracket if just in the beginning
    tmploc = sort(union(find(str == '['), find(str == ']')));
    if ~isempty(tmploc) & str(1) == '['
        if mod(length(tmploc),2) ~= 0, str, error('Opening but no closing bracket'); end;
        tmploc = tmploc(1:2);
        str = [ str(1:tmploc(1)) '' str(tmploc(1)+1:tmploc(2)-1) '' str(tmploc(2):end) ];
    end;

    %tmploc = find(str == '"');
    %if ~isempty(tmploc)
    %    if mod(length(tmploc),2) ~= 0, str, error('Opening but no closing parenthesis'); end;
    %    for index = length(tmploc):-2:1
    %        str = [ str(1:tmploc(index-1)-1) '' str(tmploc(index-1)+1:tmploc(index)-1) '' str(tmploc(index)+1:end) ];
    %    end;
    %end;

function tokout = functionformat( tokin, refcall );
	tokout = tokin;	% default
	[test, realtokin, tail, beg] = testfunc1( tokin );
	if ~test,  [test, realtokin, tail] = testfunc2( tokin ); end;
	if test
		i1 = findstr( refcall, '%s');
		i2 = findstr( refcall(i1(1):end), '''');
		if isempty(i2) i2 = length( refcall(i1(1):end) )+1; end;
		filename  = [ realtokin refcall(i1+2:i1+i2-2)]; % concatenate filename and extension
		%disp(filename)
		if exist( filename ) % do not make link if the file does not exist
			tokout =  sprintf( [ beg '%s' tail ' ' ], realtokin, realtokin );
		end;
	end;
return;

function [test, realtokin, tail, beg] = testfunc1( tokin ) % test if is string is 'function()[,]'
	test = 0; realtokin = ''; tail = ''; beg = '';
	if ~isempty( findstr( tokin, '()' ) )
		if length(tokin)<3, return; end;
		realtokin = tokin( 1:findstr( tokin, '()' )-1);
		if length(realtokin) < (length(tokin)-2) tail = tokin(end); else tail = []; end;
		test = 1;
		if realtokin(1) == '(', realtokin = realtokin(2:end); beg = '('; end;
		if realtokin(1) == ',', realtokin = realtokin(2:end); beg = ','; end;
	end;
return;

function [test, realtokin, tail] = testfunc2( tokin ) % test if is string is 'FUNCTION[,]'
	test = 0; realtokin = ''; tail = '';
	if all( upper(tokin) == tokin)
		if tokin(end) == ','
			realtokin = tokin(1:end-1);
			tail = ',';
		else
			realtokin = tokin;
		end;
		testokin = realtokin;
		testokin(findstr(testokin, '_')) = 'A';
		testokin(findstr(testokin, '2')) = 'A';
		if all(double(testokin) > 64) & all(double(testokin) < 91)
			test = 1;
		end;
		realtokin = lower(realtokin);
	end;
return;

function [tok, str] = mystrtok(str)

    [tok str] = strtok(str);
    if tok(1) == '"'
        while tok(end) ~= '"'
            [tok2 str] = strtok(str);
            if isempty(tok2), tok, error('can not find closing quote ''"'' in previous text'); end;
            tok = [tok ' ' tok2];
        end;
    end;

author: PatternRecognition

repo: OpenBMI


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [value] = BinarySABR_3(f, k, t, a, b, r, n)
% the value of a digital option within a SABR model
    epsilon = 1e-004;
    kp = k + epsilon;
    km = k - epsilon;
    sigmap = svol_2(a,b,r,n,f,kp,t);
    sigmam = svol_2(a,b,r,n,f,km,t);
    d1p = 1./(sigmap .* sqrt(t)).*(log(f./kp)+0.5*sigmap.^2*t);
    d2p = 1./(sigmap .* sqrt(t)).*(log(f./kp)-0.5*sigmap.^2*t);
    d1m = 1./(sigmam .* sqrt(t)).*(log(f./km)+0.5*sigmam.^2*t);
    d2m = 1./(sigmam .* sqrt(t)).*(log(f./km)-0.5*sigmam.^2*t);
    pp = f * normcdf(d1p) - kp .* normcdf(d2p);
    pm = f * normcdf(d1m) - km.* normcdf(d2m);
    value = (pm-pp) / (2*epsilon);
end

author: Sable

repo: mcbench-benchmarks


text:
function biharmF = biharmonic(f)
%BIHARMONIC   Biharmonic operator of a DISKFUN.
%   B = BIHARMONIC(F) returns a DISKFUN representing the biharmonic
%   operator applied to F.
%
% See also DISKFUN/BIHARM.

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

% Use the implementation in biharm.
biharmF = biharm(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function coords = roiGrowCoords(ui);
% Grow Coords from a seed in mrViewer.
%
% coords = roiGrowCoords(ui);
%
% Get a point from the user from which to grow a 'blob' of
% interconnected points in 3D, and return a 3xN list of coordinates
% for those points (in the base MR data).
% ras 09/2005.
figure(ui.fig);
tmp = get(ui.panels.display,'BackgroundColor');
set(ui.panels.display,'BackgroundColor','y');
hmsg = mrMessage('Select point from which to grow.');
pt = round(ginput(1));
close(hmsg);
set(ui.panels.display,'BackgroundColor',tmp);

% figure out what image was clicked on
imgNum = find(ui.display.axes==gca);
if isempty(imgNum)
    warning('Didn''t click on mrViewer display. No points added.')
    coords = [];
    return
end

% find seed location in the data coords
imgCoords = mrViewGet(ui,'DataCoords',imgNum);
imgSz = size(ui.display.images{imgNum});
ind = sub2ind(imgSz,pt(2),pt(1));
seed = round(imgCoords(:,ind));

% First-pass method of growing a blob efficiently:
% (1) Create an initial box, centered on seed
% (2) Restrict the box to coords which pass overlay thresholds
% (3) Find contiguous blob which contains seed
% (4) Test if box contains whole blob: if not, repeat 1-4 iteratively
hmsg = msgbox('Growing ROI...');
boxSize = 10; stepSize = 10;
coords = [];
while 1
    % (1) create box (get coords specifying box in data)
    ymin = max(1,seed(1)-boxSize);
    ymax = min(ui.mr.dims(1),seed(1)+boxSize);
    xmin = max(1,seed(2)-boxSize);
    xmax = min(ui.mr.dims(2),seed(2)+boxSize);
    zmin = max(1,seed(3)-boxSize);
    zmax = min(ui.mr.dims(3),seed(3)+boxSize);
    [X Y Z] = meshgrid(xmin:xmax,ymin:ymax,zmin:zmax);
    boxCoords = [Y(:) X(:) Z(:)]'; clear X Y Z
    mask = logical(zeros(ui.mr.dims(1:3)));

    % (2) restrict box to overlays
    boxCoords = mrViewRestrict(ui,boxCoords);
    ok = sub2ind(size(mask),boxCoords(1,:),boxCoords(2,:),boxCoords(3,:));
    mask(round(ok)) = 1;

    % (3) Find contiguous blob which contains seed
    L = bwlabeln(mask, 6); % integer label matrix of 6-connected blobs
    seedLabel = L(seed(1),seed(2),seed(3));
    if seedLabel==0, return;   end;  % no data so quit
    blob = (L==seedLabel); % binary matrix where blob is 1
    [i1, i2, i3] = ind2sub(size(mask),find(blob>0.5));
    coords = [i1 i2 i3]';

    % (4) Test if the box contains the entire blob
    % (Or is otherwise flush with the bounds of the data)
    inBoundsY=(all(i1>ymin)|ymin==1) & (all(i1<ymax)|ymax==ui.mr.dims(1));
    inBoundsX=(all(i2>xmin)|xmin==1) & (all(i2<xmax)|xmax==ui.mr.dims(2));
    inBoundsZ=(all(i3>zmin)|zmin==1) & (all(i3<zmax)|zmax==ui.mr.dims(3));
    if inBoundsX & inBoundsY & inBoundsZ
        break;
    else
        boxSize = boxSize + stepSize;
    end
end
close(hmsg);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [J grad] = nnCostFunction(nn_params, ...
                                   input_layer_size, ...
                                   hidden_layer_size, ...
                                   num_labels, ...
                                   X, y, lambda)
%NNCOSTFUNCTION Implements the neural network cost function for a two layer
%neural network which performs classification
%   [J grad] = NNCOSTFUNCTON(nn_params, hidden_layer_size, numn,/_labels, ...
%   X, y, lambda) computes the cost and gradient of the neural network. The
%   parameters for the neural network are "unrolled" into the vector
%   nn_params and need to be converted back into the weight matrices.
%
%   The returned parameter grad should be a "unrolled" vector of the
%   partial derivatives of the neural network.
%

% Reshape nn_params back into the parameters Theta1 and Theta2, the weight matrices
% for our 2 layer neural network
Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), ...
                 hidden_layer_size, (input_layer_size + 1));

Theta2 = reshape(nn_params((1 + (hidden_layer_size * (input_layer_size + 1))):end), ...
                 num_labels, (hidden_layer_size + 1));

% Setup some useful variables
m = size(X, 1);

% You need to return the following variables correctly
J = 0;
Theta1_grad = zeros(size(Theta1));
Theta2_grad = zeros(size(Theta2));

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the code by working through the
%               following parts.
%
% Part 1: Feedforward the neural network and return the cost in the
%         variable J. After implementing Part 1, you can verify that your
%         cost function computation is correct by verifying the cost
%         computed in ex4.m
%
% Part 2: Implement the backpropagation algorithm to compute the gradients
%         Theta1_grad and Theta2_grad. You should return the partial derivatives of
%         the cost function with respect to Theta1 and Theta2 in Theta1_grad and
%         Theta2_grad, respectively. After implementing Part 2, you can check
%         that your implementation is correct by running checkNNGradients
%
%         Note: The vector y passed into the function is a vector of labels
%               containing values from 1..K. You need to map this vector into a
%               binary vector of 1's and 0's to be used with the neural network
%               cost function.
%
%         Hint: We recommend implementing backpropagation using a for-loop
%               over the training examples if you are implementing it for the
%               first time.
%
% Part 3: Implement regularization with the cost function and gradients.
%
%         Hint: You can implement this around the code for
%               backpropagation. That is, you can compute the gradients for
%               the regularization separately and then add them to Theta1_grad
%               and Theta2_grad from Part 2.
%
% FORWARD PROPOGATION
a1 = [ones(m,1) X];
z2 = (a1*Theta1');
a2 = [ones(size(z2,1),1) sigmoid(z2)];
h_theta = sigmoid(a2*Theta2');
a3 = h_theta;
y_matrix = eye(num_labels)(y,:);
J = (-sum(sum(y_matrix.*log(h_theta))) - sum(sum((1-y_matrix).*(log(1-h_theta)))))/m;

% REGULARIZATION
regularization_term = (lambda/(2*m))*((sum(sum(Theta1(:,2:end).^2))) + sum(sum(Theta2(:,2:end).^2)));
J = J + regularization_term;

% BACK PROPOGATION
d3 = a3 - y_matrix;                                             % has same dimensions as a3
d2 = (d3*Theta2).*[ones(size(z2,1),1) sigmoidGradient(z2)];     % has same dimensions as a2

D1 = d2(:,2:end)' * a1;    % has same dimensions as Theta1
D2 = d3' * a2;    % has same dimensions as Theta2

Theta1_grad = Theta1_grad + (1/m) * D1;
Theta2_grad = Theta2_grad + (1/m) * D2;

% REGULARIZATION OF THE GRADIENT

Theta1_grad(:,2:end) = Theta1_grad(:,2:end) + (lambda/m)*(Theta1(:,2:end));
Theta2_grad(:,2:end) = Theta2_grad(:,2:end) + (lambda/m)*(Theta2(:,2:end));
% -------------------------------------------------------------

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

% Unroll gradients
grad = [Theta1_grad(:) ; Theta2_grad(:)];

end

author: anirudhjayaraman

repo: Machine-Learning


text:
function [  Z, H, dnorm, AC, MIhat  ] = seminmf( X, k, varargin )
% Matrix sizes
% X: m x n
% Z: m x num_of_components
% H: num_of_components x num_of_components

% Process optional arguments
pnames = {'z0' 'h0' 'bUpdateH' 'maxiter' 'TolFun' 'bUpdateZ' 'verbose'};

% Do SVD initialisation of the init components

[z0, h0] = NNDSVD(abs(X), k);
%rng(0);
%h0 = abs(rand(k , size(X, 2)));
%z0 = X * pinv(h0) + eps;

% z0 = rand(size(X, 1), k);

dflts  = {z0, h0, 1, 300,  1e-5, 1, 1};

[Z, H, bUpdateH, max_iter, tolfun, bUpdateZ, verbose] = ...
        internal.stats.parseArgs(pnames,dflts,varargin{:});

for i = 1:max_iter

    if bUpdateZ
        Z = X * pinv(H);
    end

    A = Z' * X;
    Ap = (abs(A)+A)./2;
    An = (abs(A)-A)./2;

    B = Z' * Z;
    Bp = (abs(B)+B)./2;
    Bn = (abs(B)-B)./2;

    if bUpdateH
        H = H .* sqrt((Ap + Bn * H) ./ (An + Bp * H + eps));
    end

    if mod(i, 10) == 0 || mod(i+1, 10) == 0

        s = X - Z * H;
        dnorm = sqrt(sum(s(:).^2));
        % dnorm = norm(gX - Z * H, 'fro');

        if mod(i+1, 10) == 0
            dnorm0 = dnorm;
            continue
        end

        if mod(i, 100) == 0 && verbose
            display(sprintf('...Semi-NMF iteration #%d out of %d, error: %f\n', i, max_iter, dnorm));
        end

        if 0 && exist('dnorm0')
            assert(dnorm <= dnorm0, sprintf('Rec. error increasing! From %f to %f. (%d)', dnorm0, dnorm, k));
        end

        % Check for convergence
        if exist('dnorm0') && dnorm0-dnorm <= tolfun*max(1,dnorm0)
            if verbose
                display(sprintf('Stopped at %d: dnorm: %f, dnorm0: %f', i, dnorm, dnorm0));
            end
            break;
        end

    end
end

author: andrewssobral

repo: lrslibrary


text:
function args = cvx_parse
try
    global cvx___ %#ok
    oargs = cvx___.args{1};
    nargs = length( oargs );
    if isequal( cvx___.args{2}, 'dual' ),
        nvars = nargs;
        rpat = '^([a-zA-Z]\w*)({.*})?$';
        darg = {};
    else
        nvars = cvx___.args{2};
        darg = cvx___.args{3};
        rpat = '^([a-zA-Z]\w*)(\(.*\))?$';
        if nvars < 0, nvars = nargs; end
    end
    args = struct( 'name', cell( 1, nargs ), 'args', cell( 1, nargs ) );
    for k = 1 : nargs,
        tok = regexp( oargs{k}, rpat, 'tokens' );
        if isempty( tok ),
            if k <= nvars, type = 'variable'; else type = 'structure'; end
            cvx_throw( 'Invalid %s specification: %s', type, oargs{k} );
        end
        tok = tok{1};
        nam = tok{1};
        if length(cvx___.problems) > 1 && ( nvars < 0 || k <= nvars ),
            if ~isvarname( nam ),
                cvx_throw( 'Invalid variable name: %s', nam );
            elseif nam(end) == '_',
                cvx_throw( 'Invalid variable name: %s\n   Variables ending in underscores are reserved for internal use.', nam );
            elseif isfield( cvx___.reswords, nam ),
                if cvx___.reswords.(nam) == 'S',
                    cvx_throw( 'Invalid variable name: %s\n   This is a reserved word in CVX.\n   Trying to declare a structured matrix? Use the VARIABLE keyword instead.', nam );
                else
                    cvx_throw( 'Invalid variable name: %s\n   This is a reserved word in CVX.', nam );
                end
            end
            switch exist( nam ), %#ok
                case {0,1},
                case 5,
                    cvx_throw( 'Variable name "%s" is the name of a built-in MATLAB function.\nPlease choose a different name.', nam );
                case 8,
                    cvx_throw( 'Variable name "%s" is the name of an existing MATLAB class.\nPlease choose a different name.', nam );
                otherwise,
                    mpath = which( nam );
                    if ~isempty( nam ),
                        if strncmp( nam, matlabroot, length(matlabroot) ),
                            cvx_throw( 'Variable name "%s" is the name of an existing MATLAB function or directory:\n    %s\nPlease choose a different name.', nam, mpath );
                        elseif strncmp( mpath, cvx___.where, length(cvx___.where) ),
                            cvx_throw( 'Variable name "%s" matches the name of a CVX function or directory:\n    %s\nPlease choose a different name.', nam, mpath );
                        else
                            warning( 'Variable name "%s" matches the name of an function or directory:\n    %s\nThis may cause unintended behavior with CVX models.\nPlease consider moving this file or choosing a different variable name.', nam, mpath );
                        end
                    end
            end
        end
        if length(tok) < 2 || isempty( tok{2} ),
            if k <= nvars, arg = darg;
            else arg = {}; end
        else
            try
                if k <= nvars,
                    arg = evalin( 'caller', [ '[', tok{2}(2:end-1), ']' ] );
                    arg = cvx_get_dimlist( arg );
                else
                    arg = evalin( 'caller', [ '{', tok{2}(2:end-1), '}' ] );
                end
            catch exc
                cvx_throw( 'Error attempting to evaluate arguments of: %s\n   %s', oargs{k}, exc.message );
            end
        end
        args(k).name = nam;
        args(k).args = arg;
        args(k).orig = oargs{k};
    end
    cvx___.args = {};
catch exc
    cvx___.args = {};
    rethrow( exc );
end

author: yu-jiang

repo: radpbook

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


text:
function []=mavardisp(beta_median,beta_std,beta_lbound,beta_ubound,psi_median,psi_std,psi_lbound,psi_ubound,sigma_median,X,Z,Y,n,m,p,k1,k3,q1,q2,T,ar,lambda1,lambda2,lambda3,lambda4,lambda5,IRFt,beta_gibbs,endo,exo,psi_gibbs,startdate,enddate,stringdates1,decimaldates1,datapath)

% function []=mavardisp(beta_median,beta_std,beta_lbound,beta_ubound,psi_median,psi_std,psi_lbound,psi_ubound,sigma_median,X,Z,Y,n,m,p,k1,k3,q1,q2,T,ar,lambda1,lambda2,lambda3,lambda4,lambda5,IRFt,beta_gibbs,endo,exo,psi_gibbs,startdate,enddate,stringdates1,decimaldates1,datapath)
% displays estimation results for the MABVAR model on Matlab prompt, creates a copy of these results on the text file results.txt
% inputs:  - vector 'beta_median': median value of the posterior distribution of beta
%          - vector 'beta_std': standard deviation of the posterior distribution of beta
%          - vector 'beta_lbound': lower bound of the credibility interval of beta
%          - vector 'beta_ubound': upper bound of the credibility interval of beta
%          - vector 'psi_median': median value of the posterior distribution of psi
%          - vector 'psi_std': standard deviation of the posterior distribution of psi
%          - vector 'psi_lbound': lower bound of the credibility interval of psi
%          - vector 'psi_ubound': upper bound of the credibility interval of psi
%          - vector 'sigma_median': median value of the posterior distribution of sigma (vectorised)
%          - matrix 'X': matrix of regressors for the VAR model (defined in 3.5.10)
%          - matrix 'Z': matrix of exogenous regressors for the MABVAR model (defined in 3.5.10)
%          - matrix 'Y': matrix of regressands for the VAR model (defined in 3.5.10)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k1': number of endogenous coefficients to estimate for each equation in the MABVAR model (defined p 77 of technical guide)
%          - integer 'k3': number of exogenous coefficients to estimate for each equation in the reformulated MABVAR model (defined p 77 of technical guide)
%          - integer 'q1': total number of endogenous coefficients to estimate in the MABVAR model (defined p 77 of technical guide)
%          - integer 'q2': total number of exogenous coefficients to estimate in the MABVAR model (defined p 77 of technical guide)
%          - integer 'T': number of sample time periods (defined p 7 of technical guide)
%          - scalar 'lambda1': overall tightness hyperparameter (defined p 16 of technical guide)
%          - scalar 'lambda2': cross-variable weighting hyperparameter(defined p 16 of technical guide)
%          - scalar 'lambda3': lag decay hyperparameter (defined p 16 of technical guide)
%          - scalar 'lambda4': exogenous variable tightness hyperparameter (defined p 17 of technical guide)
%          - scalar 'lambda5': block exogeneity shrinkage hyperparameter (defined p 32 of technical guide)
%          - integer 'IRFt': determines which type of structural decomposition to apply (none, Choleski, triangular factorization)
%          - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - cell 'endo': list of endogenous variables of the model
%          - matrix 'psi_gibbs': record of the gibbs sampler draws for the psi vector
%          - cell 'exo': list of exogenous variables of the model
%          - string 'startdate': start date of the sample
%          - string 'enddate': end date of the sample
%          - cell 'stringdates1': date strings for the sample period
%          - vector 'decimaldates1': dates converted into decimal values, for the sample period
%          - string 'datapath': user-supplied path to excel data spreadsheet
% outputs: none

% before displaying and saving the results, start estimating the evaluation measures for the model

% obtain first a point estimate betatilde of the VAR coefficients related to endogenous variables
% this is simply beta_median, which is both the mean and the median of the normal distribution
betatilde=beta_median;
Btilde=reshape(betatilde,k1,n);
% generate U, using Btilde, from (3.5.15)
% U=eye(n*m);
%    for jj=1:p
%    U=[U;kron(eye(m),Btilde((jj-1)*n+1:jj*n,:)')];
%    end

% obtain then a point estimate psitilde of the VAR coefficients related to exogenous variables
% this is simply psi_median, which is both the mean and the median of the normal distribution
psitilde=psi_median;

% % combine the two to obtain vec(DELTA'), using (3.5.14)
% vecdeltap=U*psitilde;
% % recover delta
% deltap=reshape(vecdeltap,n,k3);
% DELTA=deltap';

% use this estimate to produce predicted values for the model, following (3.5.9)
Ytilde=X*Btilde+Z*DELTA;
% then produce the corresponding residuals, using (3.5.9)
EPStilde=Y-Ytilde;

% check first whether the model is stationary, using (1.9.1)
[stationary eigmodulus]=bear.macheckstable(Btilde,n,p);

% Compute then the sum of squared residuals
% compute first the RSS matrix, defined in (1.9.5)
RSS=EPStilde'*EPStilde;
% retain only the diagonal elements to get the vector of RSSi values
rss=diag(RSS);

% Go on calculating R2
% generate Mbar
Mbar=eye(T)-ones(T,T)/T;
% then compute the TSS matrix, defined in (1.9.8)
TSS=Y'*Mbar*Y;
% generate the R2 matrix in (1.9.9)
R2=eye(n)-RSS./TSS;
% retain only the diagonal elements to get the vector of R2 values
r2=diag(R2);

% then calculate the adjusted R2, using (1.9.11)
R2bar=eye(n)-((T-1)/(T-(k1+m)))*(eye(n)-R2);
% retain only the diagonal elements to get the vector of R2bar values
r2bar=diag(R2bar);

% now start displaying and saving the results

% preliminary task: create and open the txt file used to save the results

filelocation=datapath;
filelocation=[filelocation '\results.txt'];
fid=fopen(filelocation,'wt');

% print toolbox header

fprintf('%s\n','');
fprintf(fid,'%s\n','');

fprintf('%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
fprintf(fid,'%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
fprintf('%s\n','%                                                                                %%');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    BAYESIAN ESTIMATION, ANALYSIS AND REGRESSION (BEAR) TOOLBOX                                         %');
fprintf(fid,'%s\n','%    BAYESIAN ESTIMATION, ANALYSIS AND REGRESSION (BEAR) TOOLBOX                                     %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    This statistical package has been developed by the external developments division of the ECB.       %');
fprintf(fid,'%s\n','%    This statistical package has been developed by the external developments division of the ECB.   %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    Authors:                                                                                %');
fprintf(fid,'%s\n','%    Authors:                                                                                %');
fprintf('%s\n','%    Romain Legrand  (Romain Legrand <b00148883@essec.edu>)                                              %');
fprintf(fid,'%s\n','%    Romain Legrand  (Romain Legrand <b00148883@essec.edu>)                                          %');
fprintf('%s\n','%    Alistair Dieppe (adieppe@worldbank.org)                                                             %');
fprintf(fid,'%s\n','%    Alistair Dieppe (adieppe@worldbank.org)                                                         %');
fprintf('%s\n','%    Björn van Roye  (Bjorn.van_Roye@ecb.europa.eu)                                                      %');
fprintf(fid,'%s\n','%    Björn van Roye  (Bjorn.van_Roye@ecb.europa.eu)                                                  %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    Version 4.4                                                                                %');
fprintf(fid,'%s\n','%    Version 4.4                                                                                %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    The authors are grateful to Paolo Bonomolo, Marta Banbura, Martin Bruns, Fabio Canova,              %');
fprintf(fid,'%s\n','%    The authors are grateful to Paolo Bonomolo, Marta Banbura, Martin Bruns, Fabio Canova,          %');
fprintf('%s\n','%    Matteo Ciccarelli, Marek Jarocinski, Niccolo Battistini, Gabriel Bobeica                            %');
fprintf(fid,'%s\n','%    Matteo Ciccarelli, Marek Jarocinski, Niccolo Battistini, Gabriel Bobeica                        %');
fprintf('%s\n','%    Michele Lenza, Chiara Osbat, Mirela Miescu, Gary Koop, Giorgio Primiceri                            %');
fprintf(fid,'%s\n','%    Michele Lenza, Chiara Osbat, Mirela Miescu, Gary Koop, Giorgio Primiceri,                       %');
fprintf('%s\n','%    Michal Rubaszek, Barbara Rossi, Ben Schumann, Peter Welz, Hugo Vega de la Cruz and Francesca Loria. %');
fprintf(fid,'%s\n','%  Michal Rubaszek, Barbara Rossi, Ben Schumann, Peter Welz, Hugo Vega de la Cruz and Francesca Loria%');
fprintf('%s\n','%    valuable input and advice which contributed to improve the quality of this work.                    %');
fprintf(fid,'%s\n','%  valuable input and advice which contributed to improve the quality of this work.                  %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%   These programmes are the responsibilities of the authors and not of the ECB and the Worldbank.       %');
fprintf(fid,'%s\n','%   These programmes are the responsibilities of the authors and not of the ECB and the Worldbank.   %');
fprintf('%s\n','%   Errors and ommissions remain those of the authors.                                                   %');
fprintf(fid,'%s\n','%   Errors and ommissions remain those of the authors.                                               %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    Please do not use or quote this work without permission.                                            %');
fprintf(fid,'%s\n','%    Please do not use or quote this work without permission.                                        %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
fprintf(fid,'%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');

% print then estimation results

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

toolboxinfo='BEAR toolbox estimates';
fprintf('%s\n',toolboxinfo);
fprintf(fid,'%s\n',toolboxinfo);

time=clock;
datestring=datestr(time);
dateinfo=['Date: ' datestring(1,1:11) '   Time: ' datestring(1,13:17)];
fprintf('%s\n',dateinfo);
fprintf(fid,'%s\n',dateinfo);

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

VARtypeinfo='Mean-adjusted BVAR';
fprintf('%s\n',VARtypeinfo);
fprintf(fid,'%s\n',VARtypeinfo);

if IRFt==1
SVARinfo='structural decomposition: none';
elseif IRFt==2
SVARinfo='structural decomposition: choleski factorisation';
elseif IRFt==3
SVARinfo='structural decomposition: triangular factorisation';
elseif IRFt==4
SVARinfo='structural decomposition: sign restrictions';
end
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);

temp='endogenous variables: ';
for ii=1:n
temp=[temp ' ' endo{ii,1} ' '];
end
endoinfo=temp;
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

temp='exogenous variables: constant';
% if m>1
%    for ii=1:m-1
%    temp=[temp ' ' exo{ii,1} ' '];
%    end
% end
exoinfo=temp;
fprintf('%s\n',exoinfo);
fprintf(fid,'%s\n',exoinfo);

sampledateinfo=['estimation sample: ' startdate '-' enddate];
fprintf('%s\n',sampledateinfo);
fprintf(fid,'%s\n',sampledateinfo);

samplelengthinfo=['sample size (omitting initial conditions): ' num2str(T)];
fprintf('%s\n',samplelengthinfo);
fprintf(fid,'%s\n',samplelengthinfo);

laginfo=['number of lags included in regression: ' num2str(p)];
fprintf('%s\n',laginfo);
fprintf(fid,'%s\n',laginfo);

hyperparam1='hyperparameters:';
fprintf('%s\n',hyperparam1);
fprintf(fid,'%s\n',hyperparam1);

hyperparam2=['autoregressive coefficient (ar):          ' num2str(ar)];
fprintf('%s\n',hyperparam2);
fprintf(fid,'%s\n',hyperparam2);

hyperparam3=['overall tightness (lambda1):              ' num2str(lambda1)];
fprintf('%s\n',hyperparam3);
fprintf(fid,'%s\n',hyperparam3);

hyperparam4=['cross-variable weighting (lambda2):       ' num2str(lambda2)];
fprintf('%s\n',hyperparam4);
fprintf(fid,'%s\n',hyperparam4);

hyperparam5=['lag decay (lambda3):                      ' num2str(lambda3)];
fprintf('%s\n',hyperparam5);
fprintf(fid,'%s\n',hyperparam5);

hyperparam6=['exogenous variable tightness (lambda4):   ' num2str(lambda4)];
fprintf('%s\n',hyperparam6);
fprintf(fid,'%s\n',hyperparam6);

hyperparam7=['block exogeneity shrinkage (lambda5):     ' num2str(lambda5)];
fprintf('%s\n',hyperparam7);
fprintf(fid,'%s\n',hyperparam7);

% display coefficient estimates

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

coeffinfo=['VAR coefficients (beta): posterior estimates'];
fprintf('%s\n',coeffinfo);
fprintf(fid,'%s\n',coeffinfo);

for ii=1:n

fprintf('%s\n','');
fprintf(fid,'%s\n','');
if ii~=1
fprintf('%s\n','');
fprintf(fid,'%s\n','');
end

endoinfo=['Endogenous: ' endo{ii,1}];
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

coeffheader=fprintf('%25s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');
coeffheader=fprintf(fid,'%25s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');

% handle the endogenous
   for jj=1:n
      for kk=1:p
      values=[beta_median((ii-1)*k1+n*(kk-1)+jj,1) beta_std((ii-1)*k1+n*(kk-1)+jj,1) beta_lbound((ii-1)*k1+n*(kk-1)+jj,1) beta_ubound((ii-1)*k1+n*(kk-1)+jj,1)];
      fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1},'(-',int2str(kk),')'),values);
      fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1},'(-',int2str(kk),')'),values);
      end
   end

% handle the exogenous
% display the results related to the constant
values=[psi_median(ii,1) psi_std(ii,1) psi_lbound(ii,1) psi_ubound(ii,1)];
fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n','Constant',values);
fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n','Constant',values);
   % if there is no other exogenous, stop here
   if m==1
   % if there are other exogenous, display their results
   elseif m>1
      for jj=1:m-1
      values=[psi_median(n*jj+ii,1) psi_std(n*jj+ii,1) psi_lbound(n*jj+ii,1) psi_ubound(n*jj+ii,1)];
      %fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
     %fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
     end
   end

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display evaluation measures
rssinfo=['Sum of squared residuals: ' num2str(rss(ii,1),'%.2f')];
fprintf('%s\n',rssinfo);
fprintf(fid,'%s\n',rssinfo);

r2info=['R-squared: ' num2str(r2(ii,1),'%.3f')];
fprintf('%s\n',r2info);
fprintf(fid,'%s\n',r2info);

adjr2info=['adj. R-squared: ' num2str(r2bar(ii,1),'%.3f')];
fprintf('%s\n',adjr2info);
fprintf(fid,'%s\n',adjr2info);

end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display VAR stability results
eigmodulus=reshape(eigmodulus,p,n);
stabilityinfo1=['Roots of the characteristic polynomial (modulus):'];
fprintf('%s\n',stabilityinfo1);
fprintf(fid,'%s\n',stabilityinfo1);
for ii=1:p
temp=num2str(eigmodulus(ii,1),'%.3f');
   for jj=2:n
   temp=[temp,'  ',num2str(eigmodulus(ii,jj),'%.3f')];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end
if stationary==1;
stabilityinfo2=['No root lies outside the unit circle.'];
stabilityinfo3=['The estimated VAR model satisfies the stability condition'];
fprintf('%s\n',stabilityinfo2);
fprintf(fid,'%s\n',stabilityinfo2);
fprintf('%s\n',stabilityinfo3);
fprintf(fid,'%s\n',stabilityinfo3);
else
stabilityinfo2=['Warning: at leat one root lies on or outside the unit circle.'];
stabilityinfo3=['The estimated VAR model will not be stable'];
fprintf('%s\n',stabilityinfo2);
fprintf(fid,'%s\n',stabilityinfo2);
fprintf('%s\n',stabilityinfo3);
fprintf(fid,'%s\n',stabilityinfo3);
end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display posterior for sigma
sigmainfo=['sigma (residual covariance matrix): posterior estimates'];
fprintf('%s\n',sigmainfo);
fprintf(fid,'%s\n',sigmainfo);
% calculate the (integer) length of the largest number in sigma, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(sigma_median))))));
% add a separator, a potential minus sign, and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:n
temp=[];
   for jj=1:n
   % convert matrix entry into string
   number=num2str(sigma_median(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

fclose(fid);

% Finally, display the results in terms of graph

% graph the empirical posterior distribution of the VAR parameters
% first reshape the gibbs sampler output to be compatible with the way Matlab uses subplot
% first concatenate psi_gibbs on beta_gibbs to obtain a vector similar to beta in the classical BVAR model
temp=[];
for ii=1:n
temp=[temp;beta_gibbs((ii-1)*k1+1:ii*k1,:)];
   for jj=1:m
   temp=[temp;psi_gibbs((jj-1)*m+ii,:)];
   end
end
beta_reshape=temp;
[beta_swap]=bear.betaswap(beta_reshape,n,m,p,k1+m);
plotvar=1;
plotlag=1;
plotexo=1;
plotconst=1;
% then plot the figure
%{
postdis=figure;
set(postdis,'Color',[0.9 0.9 0.9]);
set(postdis,'name','posterior ditribution of VAR coefficients');
for ii=1:(q1+q2)
subplot(n,k1+m,ii)
hist(beta_swap(ii,:),20);
histproperties=findobj(gca,'Type','patch');
set(histproperties,'FaceColor',[0.7 0.78 1],'EdgeColor',[0 0 0],'LineWidth',0.1);
set(gca,'FontName','Times New Roman');
   % top labels for endogenous
   if ii<=n*p
   temp=[endo{plotvar,1} '(-' num2str(plotlag) ')'];
   title(temp,'FontWeight','normal');
      if plotlag<p
      plotlag=plotlag+1;
      elseif plotlag==p
      plotlag=1;
      plotvar=plotvar+1;
      end
   end
   % top labels for exogenous
   %if ii==k1+1
   %title('Constant','FontWeight','normal');
   %end
   %if ii>k1+1 && ii<=k1+m
   %title(exo{plotexo,1},'FontWeight','normal');
   %plotexo=plotexo+1;
   %end
   % side labels
   %if rem((ii-1)/(k1+m),1)==0
   ylabel(endo{(ii-1)/(k1+m)+1,1},'FontWeight','normal');
   %end
end
%}

% plot actual vs. fitted
actualfitted=figure;
set(actualfitted,'Color',[0.9 0.9 0.9]);
set(actualfitted,'name','model estimation: actual vs fitted')
ncolumns=ceil(n^0.5);
nrows=ceil(n/ncolumns);
for ii=1:n
subplot(nrows,ncolumns,ii)
hold on
plot(decimaldates1,Y(:,ii),'Color',[0 0 0],'LineWidth',2);
plot(decimaldates1,Ytilde(:,ii),'Color',[1 0 0],'LineWidth',2);
hold off
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
title(endo{ii,1},'FontName','Times New Roman','FontSize',10);
   if ii==1
   plotlegend=legend('actual','fitted');
   set(plotlegend,'FontName','Times New Roman');
   end
end

% plot the residuals
residuals=figure;
set(residuals,'Color',[0.9 0.9 0.9]);
set(residuals,'name','model estimation: residuals')
for ii=1:n
subplot(nrows,ncolumns,ii)
plot(decimaldates1,EPStilde(:,ii),'Color',[0 0 0],'LineWidth',2)
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
title(endo{ii,1},'FontName','Times New Roman','FontSize',10);
end

% finally, save the results on excel
bear.data.excelrecord2fcn(T, n, endo, stringdates1, Y, Ytilde, pref, EPStilde)

author: european-central-bank

repo: BEAR-toolbox


text:
%% Example 2.7: Hyperbolic Conservation Laws in 2D
%
% In this example we will solve the 2D scalar conservation law
%
% $$ u_t + f(u)_x + g(u)_y = 0, \qquad u(x,y,0) = u_0(x,y) $$
%
% by dimensional splitting. To this end, we use two substeps,
%
% $$ S_x(t): \quad v_t + f(v)_x = 0, \qquad v(x,0)=v_0(x)$$
%
% $$ S_y(t): \quad w_t + g(w)_y = 0, \qquad w(y,0)=w_0(y)$$
%
% and construct the approximate solution from the formula
%
% $$u(x,t)\approx [S_y(\Delta t)\circ S_x(\Delta t) ]^n u_0(x)$$
%
% The 1-D hyperbolic steps will be solved using the Lax-Friedrichs scheme

%% Initial setup
N = 256;
h = 2*pi/N;
x =-pi+(0:N)*h;  x=0.5*(x(1:end-1)+x(2:end));
y = x;
[X,Y] = ndgrid(x,y);
u0 = exp( -4*sin(X/2).^2 - 4*sin(Y/2).^2 );
h=path; path(h,'../Example2_5');

%% Flux functions
% As a concrete example, we choose the Buckley-Leverett flux in the
% x-direction and a convex flux in the y-direction
%
% $$f(u) = u^2 / (u^2 + (1-u)^2), \qquad g(u) = u(1-u)$$
%
% This system can be seen as a simplified model of two-phase flow in a
% periodic porous medium under the influence of convective and gravity
% forces (convective in the x-direction and gravity in the y-direction).
% The two fluxes are shown in the plot below
u = linspace(0,1,101);
plot(u, fflux(u), u, gflux(u)), axis tight, legend('f(u)','g(u)',2)

%% Evolution of the solution
T = 1.5;
nsplit = 50;
u = dimsplit('fflux','gflux', u0, x, y, T, nsplit);
I = 1:10:nsplit+1;
t = linspace(0,T,nsplit+1);
c = linspace(0,1,11);
for i=1:6
   subplot(2,3,i),
   contourf(x,y,u(:,:,I(i))',c), caxis([0 1]) %, colorbar('horiz')
   %surf(x,y,u(:,:,I(i))'), shading interp, view(3)
   title(['Time: t=' num2str(t(I(i)))]);
end
%%
% To understand the evolution of the initial blob, we discuss the movement
% in the two spatial directions separately. In the x-direction, the
% nonconvex Buckley-Leverett flux will turn the right-hand side of the
% symmetric blob into a leading shock, followed by a rarefaction. Likewise,
% the left-hand side of the blob turns into a shock followed by a
% rarefaction wave. In the y-direction, the convex flux gives a shock along
% the upper part of the blob and a rarefaction wave along the lower part.

%% Compare different time steps
% Having established the dynamics of the problem, we can investigate the
% performance of our numerical method for different choices of the
% splitting step. To this end, we increase time interval to [0,8].
T = 8;
c = linspace(0,.25,11);
for n=1:4,
	nsplit=4*2.^(n-1);
	u=dimsplit('fflux','gflux',u0,x,y,T,nsplit);
	subplot(2,2,n)
	contourf(x,y,u(:,:,nsplit+1)',c), caxis([0 .25]);
	xlabel('x'), ylabel('y','Rotation',0),  colorbar
	title([num2str(nsplit) ' steps']);
end;
path(h);
%%
% It is quite amazing how accurate one can capture the dynamics of this
% flow with only four splitting steps. Increasing the number of steps,
% improves small-scale features of the solution, but does not change the
% qualitative behavior of the approximate solution.

repo: Aero-matlab


text:
classdef MimShowMarkers < MimGuiPlugin
    % MimShowMarkers. Gui Plugin for entering marker edit mode
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimEnterMarkerMode is a Gui Plugin for showing or hiding markers
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Show Markers'
        SelectedText = 'Hide Markers'
        ToolTip = 'Enters marker edit mode'
        Category = 'Marker display'
        Visibility = 'Dataset'
        Mode = 'Markers'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Location = 1
        Icon = 'markers.png'
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.ShowMarkers = ~gui_app.ImagePanel.ShowMarkers;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded;
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.ShowMarkers;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function value = torus_volume_3d ( r1, r2 )

%*****************************************************************************80
%
%% TORUS_VOLUME_3D returns the volume of a torus in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that:
%
%      ( SQRT ( X**2 + Y**2 ) - R1 )**2 + Z**2 = R2**2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R1, R2, the two radii that define the torus.
%
%    Output, real TORUS_VOLUME_3D, the volume of the torus.
%
  value = 2.0 * pi * pi * r1 * r2 * r2;

  return
end


text:
function p = addMonomialCuts(p)

if any(p.originalModel.variabletype==3)
    singleMonomial =  sum(p.originalModel.monomtable |  p.originalModel.monomtable,2)==1;
    p_cut = emptyNumericalModel;
    monomials = find(p.originalModel.variabletype==3);
    for i = 1:length(monomials)
        monom_index = monomials(i);
        monom = p.originalModel.monomtable(monomials(i),:);
        monom_variable = find(monom);
        if length(monom_variable)==1
            n = monom(monom_variable);
            L = p.lb(monom_variable);
            U = p.ub(monom_variable);
            if ~isinf(L) && ~isinf(U)
                if even(n)
                    M = (L+U)/2;
                    if L <= 0 && U >= 0 && M~=0
                        if M < 0
                            [Ax,Ay,b,K] = convexhullConvex(L,M,0,U,L^n,M^n,0,U^n,n*L^(n-1),n*M^(n-1),0,n*U^(n-1));
                        else
                            [Ax,Ay,b,K] = convexhullConvex(L,0,M,U,L^n,0,M^n,U^n,n*L^(n-1),0,n*M^(n-1),n*U^(n-1));
                        end
                    else
                        [Ax,Ay,b,K] = convexhullConvex(L,M,U,L^n,M^n,U^n,n*L^(n-1),n*M^(n-1),n*U^(n-1));
                    end
                    p_cut.F_struc(end+1:end+length(b),1) = b;
                    p_cut.F_struc(end-length(b)+1:end,1+monom_variable) = -Ax;
                    p_cut.F_struc(end-length(b)+1:end,1+monom_index) = -Ay;
                    p_cut.K.l = p_cut.K.l+length(b);
                else
                    if p.lb(monom_variable)<0 &  p.ub(monom_variable)>0 & ~isinf(p.lb(monom_variable)) & ~isinf(p.ub(monom_variable))

                        % Line between lower bound and tangent intersection
                        r = zeros(1,n+1);r(1)=n-1;r(2)=-L*n;r(end)=L^n;
                        r = roots(r);
                        %r = r(min(find(r==real(r))));
                        r = max(r(imag(r)==0));
                        if r >= U
                            r = U;
                            fprim = (U^n-L^n)/(U-L);
                        else
                            fprim = n*r^(n-1);
                        end
                        p_cut.F_struc(end+1,1) = -L^n+L*fprim;
                        p_cut.F_struc(end,1+monom_index)=1;
                        p_cut.F_struc(end,1+monom_variable)=-fprim;
                        p_cut.K.l = p_cut.K.l+1;

                        % Line between upper bound and tangent intersection
                        r = zeros(1,n+1);r(1)=n-1;r(2)=-U*n;r(end)=U^n;
                        r = roots(r);
                        r = min(r(imag(r)==0));
                        if r <= L
                            r = L;
                            fprim = (U^n-L^n)/(U-L);
                        else
                            fprim = n*r^(n-1);
                        end
                        p_cut.F_struc(end+1,1) = U^n-U*fprim;
                        p_cut.F_struc(end,1+monom_index)=-1;
                        p_cut.F_struc(end,1+monom_variable)=fprim;
                        p_cut.K.l = p_cut.K.l+1;
                    end
                end
            end
        end
    end
    p = mergeNumericalModels(p,p_cut);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function wavelet_test07 ( )

%*****************************************************************************80
%
%% WAVELET_TEST07 tests DAUB14_TRANSFORM and DAUB14_TRANSFORM_INVERSE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WAVELET_TEST07\n' );
  fprintf ( 1, '  DAUB14_TRANSFORM computes the DAUB14 transform of a vector.\n' );
  fprintf ( 1, '  DAUB14_TRANSFORM_INVERSE inverts it.\n' );
%
%  Random data.
%
  n = 16;
  seed = 123456789;
  [ u, seed ] = r8vec_uniform_01 ( n, seed );

  v = daub14_transform ( n, u );

  w = daub14_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)       D14(U)(i)  D14inv(D14(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Constant signal.
%
  n = 8;
  u(1:n) = 1.0;

  v = daub14_transform ( n, u );

  w = daub14_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D14(U)(i)  D14inv(D14(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Linear signal.
%
  n = 16;
  a_first = 1.0;
  a_last = n;
  u = linspace ( a_first, a_last, n );

  v = daub14_transform ( n, u );

  w = daub14_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D14(U)(i)  D14inv(D14(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Quadratic data.
%
  n = 8;
  u(1) = 25.0;
  u(2) = 16.0;
  u(3) = 9.0;
  u(4) = 4.0;
  u(5) = 1.0;
  u(6) = 0.0;
  u(7) = 1.0;
  u(8) = 4.0;

  v = daub14_transform ( n, u );

  w = daub14_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D14(U)(i)  D14inv(D14(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end

  return
end


text:
function gX = ardKernDiagGradX(kern, X)

% ARDKERNDIAGGRADX Gradient of ARD kernel's diagonal with respect to X.
% FORMAT
% DESC computes the gradient of the diagonal of the pre-built RBF and linear ARD kernel matrix with
% respect to the elements of the design matrix given in X.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG X : the input data in the form of a design matrix.
% RETURN gX : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : ardKernParamInit, kernDiagGradX, ardkernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004

% KERN

gX = 2*kern.linearVariance*X.*repmat(kern.inputScales, [size(X, 1) 1]);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%ISVEC Test if vector
%
% ISVEC(V) is true (1) if the argument V is a 3-vector, else false (0).
%
% ISVEC(V, L) is true (1) if the argument V is a vector of length L,
% either a row- or column-vector.  Otherwise false (0).
%
% Notes::
% - Differs from MATLAB builtin function ISVECTOR, the latter returns true
%   for the case of a scalar, ISVEC does not.
% - Gives same result for row- or column-vector, ie. 3x1 or 1x3 gives true.
%
% See also ISHOMOG, ISROT.

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function h = isvec(v, l)
    if nargin == 1
            l = 3;
    end
    d = size(v);
    h = logical( length(d) == 2 && min(d) == 1 && numel(v) == l );

author: Allopart

repo: rbpf-gmapping


text:
function engine = jtree_mnet_inf_engine(model, varargin)
% JTREE_MNET_INF_ENGINE Junction tree inference engine for Markov nets
% engine = jtree_inf_engine(mnet, ...)
%

% set default params
N = length(mnet.graph);
root = N;

engine = init_fields;
engine = class(engine, 'jtree_mnet_inf_engine', inf_engine(bnet));

onodes = bnet.observed;
if is_mnet(bnet)
  MG = bnet.graph;
else
  error('should be a mnet')
end

%[engine.jtree, dummy, engine.cliques, B, w, elim_order, moral_edges, fill_in_edges, strong] = ...
%    dag_to_jtree(bnet, onodes, stages, clusters);

porder = determine_elim_constraints(bnet, onodes);
strong = ~isempty(porder);
ns = bnet.node_sizes(:);
ns(onodes) = 1; % observed nodes have only 1 possible value
[engine.jtree, root2, engine.cliques, B, w] = ...
    graph_to_jtree(MG, ns, porder, stages, clusters);

engine.cliques_bitv = B;
engine.clique_weight = w;
C = length(engine.cliques);
engine.clpot = cell(1,C);

% Compute the separators between connected cliques.
[is,js] = find(engine.jtree > 0);
engine.separator = cell(C,C);
for k=1:length(is)
  i = is(k); j = js(k);
  engine.separator{i,j} = find(B(i,:) & B(j,:)); % intersect(cliques{i}, cliques{j});
end

% A node can be a member of many cliques, but is assigned to exactly one, to avoid
% double-counting its CPD. We assign node i to clique c if c is the "lightest" clique that
% contains i's family, so it can accomodate its CPD.

engine.clq_ass_to_node = zeros(1, N);
for i=1:N
  %c = clq_containing_nodes(engine, family(bnet.dag, i));
  clqs_containing_family = find(all(B(:,family(bnet.dag, i)), 2)); % all selected columns must be 1
  c = clqs_containing_family(argmin(w(clqs_containing_family)));
  engine.clq_ass_to_node(i) = c;
end

% Make the jtree rooted, so there is a fixed message passing order.
if strong
  % the last clique is guaranteed to be a strong root
  engine.root_clq = length(engine.cliques);
else
  % jtree_dbn_inf_engine requires the root to contain the interface.
  % This may conflict with the strong root requirement! *********** BUG *************
  engine.root_clq = clq_containing_nodes(engine, root);
  if engine.root_clq <= 0
    error(['no clique contains ' num2str(root)]);
  end
end

[engine.jtree, engine.preorder, engine.postorder] = mk_rooted_tree(engine.jtree, engine.root_clq);

% collect
engine.postorder_parents = cell(1,length(engine.postorder));
for n=engine.postorder(:)'
  engine.postorder_parents{n} = parents(engine.jtree, n);
end
% distribute
engine.preorder_children = cell(1,length(engine.preorder));
for n=engine.preorder(:)'
  engine.preorder_children{n} = children(engine.jtree, n);
end

%%%%%%%%

function engine = init_fields()

engine.jtree = [];
engine.cliques = [];
engine.separator = [];
engine.cliques_bitv = [];
engine.clique_weight = [];
engine.clpot = [];
engine.clq_ass_to_node = [];
engine.root_clq = [];
engine.preorder = [];
engine.postorder = [];
engine.preorder_children = [];
engine.postorder_parents = [];
engine.maximize = [];
engine.evidence = [];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% This matlab codes implement the adaptive residual subsampling
% method for radial basis function 1-D initial-boundary value problem.
%
% Implementation : Method of lines
%                  Time-stepping: 4th order Runge-Kutta method
%                  Space adaptivity: Adaptive residual subsampling for RBFs
%
% For reference, see:
% Adaptive residual subsampling methods for radial basis function
% interpolation and collocation problems. submitted to Computers Math. Appl
%
% Tobin A. Driscoll and Alfa R.H. Heryudono    05/14/2006
% MATLAB 7 is recommended.
%
% Test problem:
% Burger's Equation
% epsilon.u_xx - u.u_x = u_t, 0 < x < 1
% u(0,t)=u(1,t)=0.

Tfin = 1;
dt = Tfin/100;

epsilon = 1e-3; theta = 1e-4;
theta = [theta theta/1000];
initcond = @(x) sin(2*pi*x) + 0.5*sin(pi*x);
[x,c,u0] = coarserefine(initcond,[0 1],theta,13,0.75);

Tdone = 0;
han = plot(x,u0,'-ko',x,-1*x.^0,'ko','MarkerFaceColor','k','MarkerSize',2);
tp = title('','erasemode','xor'); hold on;
xlabel('x');ylabel('u(x,t)');
axis([0 1 -1 1.5])

while Tdone < Tfin
N = length(x);

[A,D1,D2] = deal(zeros(N));
for j=1:N
    [A(:,j),D1(:,j),D2(:,j)] = mq(x,x(j),c(j));
end

lambda = A\eye(N);
D1 = D1*lambda; D1([1 N],:)=[]; D1(:,[1 N])=[];
D2 = D2*lambda; D2([1 N],:)=[]; D2(:,[1 N])=[];

% Solve system of ODEs
options = odeset('RelTol',1e-5,'AbsTol',1e-8,'Jacobian',@Jburgers,'OutputFcn', ...
                  @(t,y,flag,varargin) burgersplot(t,y,flag,varargin,x,han));

set(tp,'string',sprintf('T = %.3f,    N = %3i.',Tdone,N));
[T,U] = ode15s(@burgers,[Tdone Tdone+dt],u0(2:N-1),options,epsilon,D1,D2);

disp('-----------------------');
[x,c,u0] = coarserefine(@(xp)predictor(xp,x,[0;U(end,:)';0]),[0 1],theta,13,0.75);
Tdone = Tdone + dt;
end


text:
function V = par_read_volume(info)
% Function for reading the volume of a Philips Par / Rec  MR V4.* file
%
% volume = par_read_volume(file-header)
%
% examples:
% 1: info = par_read_header()
%    V = par_read_volume(info);
%    imshow(squeeze(V(:,:,round(end/2),1)),[]);
%
% 2: V = par_read_volume('test.par');

if(~isstruct(info)), info=par_read_header(info); end

% Open file
fid=fopen(info.FilenameREC','rb','ieee-le');
% Skip header
fseek(fid,0,'bof');

datasize=prod(info.Dimensions)*info.BitDepth/8;

% Read the Data
switch(info.BitDepth)
    case 8
        info.DataType='char';
    case 16
        info.DataType='short';
    case 32
        info.DataType='float';
    case 64
        info.DataType='double';
end

switch(info.DataType)
    case 'char'
        V = int8(fread(fid,datasize,'char=>int8'));
    case 'uchar'
        V = uint8(fread(fid,datasize,'uchar=>uint8'));
    case 'short'
        V = int16(fread(fid,datasize,'short=>int16'));
    case 'ushort'
        V = uint16(fread(fid,datasize,'ushort=>uint16'));
    case 'int'
        V = int32(fread(fid,datasize,'int=>int32'));
    case 'uint'
        V = uint32(fread(fid,datasize,'uint=>uint32'));
    case 'float'
        V = single(fread(fid,datasize,'float=>single'));
    case 'double'
        V = double(fread(fid,datasize,'double=>double'));
end

fclose(fid);
V = reshape(V,info.Dimensions);


text:
clf reset;
echo on

% This script demonstrates the use of the RSOM.
clc;

% load the example dissimilarity data
load exampleDissimilarity.mat;

% display the eigenvalue spectrum
[V, D] = eig(Dissim);
eVals = diag(D);
figure; bar(eVals);

pause % Strike any key to continues...
clc;

% init RSOM
sMap = rsom_lininit(Dissim,  [10 10]);

% train the RSOM
sMap = rsom_batchtrain(sMap, Dissim);

% Display the U-Matrix
figure;
rsom_show(sMap, Dissim);

% The U-Matrix shows clearly the structure of the data, i.e. that two
% clusters are available

pause % Strike any key to continues...
clc;

% do linear embedding of the distance matrix into a 3 dimensional space
x = cmdscale(Dissim.^(1/2));
x = x(:,1:3);

% plot the resulting data
h = figure; hold on;
plot3(x(1:100, 1), x(1:100, 2), x(1:100, 3), '.b');
plot3(x(101:200, 1), x(101:200, 2), x(101:200, 3), '.r');

% Since approximated vectorial data are available now, we can compute
% the (approximated) neuron positions and plot them
Neurons = sMap.cCodebook * x;
figure(h);

% create a som struct
sMapSOM = som_map_struct(3, sMap.topol);
som_grid(sMapSOM,'Coord',Neurons);

echo off;

author: ilarinieminen

repo: SOM-Toolbox


text:
function [f, grad] = DL1c_regularization(Knots, kspc, D, csqrt)
    gvol = prod(kspc);
    k = 1;
    ksz = size(Knots(:,:,:, 1));
    Nd = size(Knots, 4);
    f = 0;
    grad = zeros(size(Knots));
    p = 1;
%     p = 0.85;

    for nn = 1 : size(Knots, 5)
        for i = 1 : Nd
            tmp = D * fl(Knots(:,:,:, i, nn));

            t2 = sqrt(tmp(:).^2 + csqrt);
            f = f + sum(t2.^p);

%             grad = D' * ((abs(tmp)+deps).^(p-1) .* sign(tmp)) * p;

%             gr = D' * fl(tmp ./ t2);
            gr = D' * (tmp .* (t2.^(p-2)) * p);

            grad(:,:,:, i, nn) = grad(:,:,:, i, nn) + reshape(gr, ksz);
        end
    end
    f = f * gvol;
    grad = grad * gvol;
end

repo: pTVreg


text:
function f=fitmultiplegauss_test(x, varargin)
% f=loglikGaP(x, varargin)
% complete log likellihood of the GaP model
% Vxt = varargin{1};      %data
% sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
% peval = varargin{5}; %parameters
% x(1:end-2*peval.ncomp) is Hkt

Vxyt = varargin{1};      %data
sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
peval = varargin{3}; %parameters

Vxt=reshape(Vxyt,peval.nx*peval.ny,1);
sigpsf_vec=repmat(sigpsf,peval.ncomp,1); %all psfs same sigma
a_vec=1./(sigpsf_vec.^2*2*pi); % all normalised to 1
cxy_vec=[x(1:peval.ncomp)'+1, x(peval.ncomp+1:end)'+1]; %different notation of the dip_image/new version

Wxkpix=gauss2dmultislice([peval.nx, peval.ny, peval.ncomp], cxy_vec, sigpsf_vec, a_vec);
Wxkpix=normalize(Wxkpix); %normalize PSFs to 1
Wxk=reshape(Wxkpix,peval.nx*peval.ny, peval.ncomp);

factor=(sum(Vxt(:))/sum(Wxk(:)))/peval.ncomp;
Hkt=factor*ones(peval.ncomp,1); %just sum the components

P=Wxk*Hkt+peval.bg; %current approximation

%Poisson contribution
t1=Vxt.*log(P) - P;
%Gamma contribution
% t2=(alpha-1)*log(Hkt)-1/beta*Hkt-alpha*log(beta)-log(gamma(alpha));

f=sum(t1(:));%+sum(t2(:)));
f=-f; %conjugate gradient is mimimizing!
end

author: aludnam

repo: MATLAB


text:
% Test file for trigtech/qr.m

function pass = test_qr(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
x = 2 * rand(100, 1) - 1;

testclass = trigtech();

%%
% Do a few spot-checks.

f = testclass.make(@(x) exp(sin(pi*x)), [], pref);
pass(1:2) = test_one_qr(f, x);
pass(3:4) = test_one_qr_with_perm(f, x);

f = testclass.make(@(x) [exp(sin(pi*x)) 3./(4-cos(pi*x))], [], pref);
pass(5:6) = test_one_qr(f, x);
pass(7:8) = test_one_qr_with_perm(f, x);

f = testclass.make(@(x) [ones(size(x)) sin(pi*x) cos(pi*x) sin(2*pi*x) ...
                         cos(2*pi*x) sin(3*pi*x)], [], pref);
pass(9:10) = test_one_qr(f, x);
pass(11:12) = test_one_qr_with_perm(f, x);

f = testclass.make(@(x) [3./(4-exp(1i*pi*x)) exp(sin(pi*x)) cos(3*pi*x)], ...
[], pref);
pass(13:14) = test_one_qr(f, x);
pass(15:16) = test_one_qr_with_perm(f, x);

%%
% Check that the 'vector' flag works properly.
N = size(f, 2);
[Q1, R1, E1] = qr(f, []);
[Q2, R2, E2] = qr(f, 'vector');
err = E1(E2, :) - eye(N);
pass(17) = all(err(:) == 0);

%%
% Check a rank-deficient problem:
% [TODO]: Is this correct?
%
% Rank deficient QR factorizations fail due to the bug
% in issue #1441. These tests are disabled until the bug
% is addressed.
f = testclass.make(@(x) [cos(pi*x) cos(pi*x) cos(pi*x)], [], pref);
[Q, R] = qr(f, []);
Q = simplify(Q,100*eps);
pass(18) = all(size(Q) == 3) && all(size(R) == 3);
% pass(18) = 1;
I = eye(3);
pass(19) = norm(innerProduct(Q, Q) - I, inf) < ...
10*max(vscale(f)*eps);
pass(19) = 1;
% These tests should be reverted once issue #1441 is
% fixed.

%%
% Check that the vscale comes out with the correct size for QR of an
% array-valued chebtech.
f = testclass.make(@(x) [sin(pi*x) cos(pi*x) cos(2*pi*x)], [], pref);
[Q, R] = qr(f, []);
pass(20) = isequal(size(vscale(Q)), [1 3]);

end

% Tests the QR decomposition for a CHEBTECH object F using a grid of points X
% in [-1  1] for testing samples.
function result = test_one_qr(f, x)
    N = size(f, 2);
    [Q, R] = qr(f);

    % Check orthogonality.
    ip = innerProduct(Q, Q);
    result(1) = max(max(abs(ip - eye(N)))) < 10*max(vscale(f)*eps);

    % Check that the factorization is accurate.
    err = Q*R - f;
    result(2) = norm(feval(err, x), inf) < 100*max(vscale(f)*eps);
end

% Same as the previous function but this time uses the QR factorization with
% permutations.
function result = test_one_qr_with_perm(f, x)
    N = size(f, 2);
    [Q, R, E] = qr(f);

    % Check orthogonality.
    ip = innerProduct(Q, Q);
    result(1) = max(max(abs(ip - eye(N)))) < 10*max(vscale(f)*eps);

    % Check that the factorization is accurate.
    err = Q*R - f*E;
    result(2) = norm(feval(err, x), inf) < 100*max(vscale(f)*eps);
end


text:
function [GAmp,GTime]=GxTrapezoid(p)

tStart=p.tStart; % start time ms
tEnd=p.tEnd;     % end time ms
tRamp=p.tRamp;   % ramp duration time ms
sRamp=p.sRamp;   % ramp steps
GxAmp=p.GxAmp;   % Gx amplitude
Duplicates=max(1,p.Duplicates);
DupSpacing=max(0,p.DupSpacing);

[GAmp,GTime]=StdTrap(tStart-tRamp, ...
                     tEnd+tRamp,   ...
                     tStart,               ...
                     tEnd,                 ...
                     GxAmp,max(2,sRamp),2,max(2,sRamp));

[GTime,m,n]=unique(GTime);
GAmp=GAmp(m);

% Create Duplicates
if Duplicates~=1 & DupSpacing ~=0
    GAmp=repmat(GAmp,[1 Duplicates]);
    TimeOffset = repmat(0:DupSpacing:(Duplicates-1)*DupSpacing,[length(GTime) 1]);
    GTime=repmat(GTime,[1 Duplicates]) + (TimeOffset(:))';
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab