text:
function []=panel3disp(n,N,m,p,k,T,Yi,Xi,Units,endo,exo,const,beta_gibbs,beta_median,beta_std,beta_lbound,beta_ubound,sigma_gibbs,sigma_median,D_estimates,gamma_estimates,lambda1,startdate,enddate,forecast_record,forecast_estimates,Fcperiods,stringdates3,Fstartdate,Fcenddate,Feval,Fcomp,data_endo_c,data_endo_c_lags,data_exo_c,It,Bu,IRF,IRFt,pref,names)

% start displaying and saving the general results

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

filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'wt');

% print toolbox header

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

% print the list of contributors
bear.printcontributors(fid);

% 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='Panel VAR: random effect model (Zellner and Hong)';
fprintf('%s\n',VARtypeinfo);
fprintf(fid,'%s\n',VARtypeinfo);

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

temp='units: ';
for ii=1:N
temp=[temp ' ' Units{ii,1} ' '];
end
unitinfo=temp;
fprintf('%s\n',unitinfo);
fprintf(fid,'%s\n',unitinfo);

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: ';
if const==0 && m==0
temp=[temp ' none'];
elseif const==1 && m==1
temp=[temp ' constant '];
elseif const==0 && m>0
   for ii=1:m-1
   temp=[temp ' ' exo{ii,1} ' '];
   end
elseif const==1 && m>1
temp=[temp ' constant '];
   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=['overall tightness (lambda1):              ' num2str(lambda1)];
fprintf('%s\n',hyperparam2);
fprintf(fid,'%s\n',hyperparam2);

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

% everything else is then unit-specific, hence loop over units
for hh=1:N

% first estimate the in-sample evaluation criteria
% obtain predicted values
B_median(:,:,hh)=reshape(beta_median(:,:,hh),k,n);
Yp(:,:,hh)=Xi(:,:,hh)*B_median(:,:,hh);
% then produce the corresponding residuals, using (1.9.4)
EPS(:,:,hh)=Yi(:,:,hh)-Yp(:,:,hh);

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

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

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

% now start the display phase: start with the name of the unit

unitinfo=['Unit: ' Units{hh,1}];
fprintf('%s\n',unitinfo);
fprintf(fid,'%s\n',unitinfo);

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

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

% loop over the equations(that is, endogenous variables) of the model
for ii=1:n
fprintf('%s\n','');
fprintf(fid,'%s\n','');

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)*k+n*(kk-1)+jj,1,hh) beta_std((ii-1)*k+n*(kk-1)+jj,1,hh) beta_lbound((ii-1)*k+n*(kk-1)+jj,1,hh) beta_ubound((ii-1)*k+n*(kk-1)+jj,1,hh)];
      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
   % if there is no constant:
   if const==0
      % if there is no exogenous at all, obvioulsy, don't display anything
      if m==0
      % if there is no constant but some other exogenous, display them
      else
         for jj=1:m
         values=[beta_median(ii*k-m+jj,1,hh) beta_std(ii*k-m+jj,1,hh) beta_lbound(ii*k-m+jj,1,hh) beta_ubound(ii*k-m+jj,1,hh)];
         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
   % if there is a constant
   else
   % display the results related to the constant
         values=[beta_median(ii*k-m+1,1,hh) beta_std(ii*k-m+1,1,hh) beta_lbound(ii*k-m+1,1,hh) beta_ubound(ii*k-m+1,1,hh)];
         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
      else
         for jj=1:m-1
         values=[beta_median(ii*k-m+jj+1,1,hh) beta_std(ii*k-m+jj+1,1,hh) beta_lbound(ii*k-m+jj+1,1,hh) beta_ubound(ii*k-m+jj+1,1,hh)];
         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
   end

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

% display evaluation measures

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

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

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

end

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

% display VAR stability results

% check whether the model is stationary
[stationary eigmodulus]=bear.checkstable(beta_median(:,:,hh),n,p,k);
eigmodulus=reshape(eigmodulus,p,n);
stabilityinfo1=['Roots of the characteristic polynomial (modulus):'];
fprintf('%s\n',stabilityinfo1);
fprintf(fid,'%s\n',stabilityinfo1);
for jj=1:p
temp=num2str(eigmodulus(jj,1),'%.3f');
   for kk=2:n
   temp=[temp,'  ',num2str(eigmodulus(jj,kk),'%.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.'];
fprintf('%s\n',stabilityinfo2);
fprintf(fid,'%s\n',stabilityinfo2);
stabilityinfo3=['The estimated VAR model will not be stable'];
fprintf('%s\n',stabilityinfo3);
fprintf(fid,'%s\n',stabilityinfo3);
end

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

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

% then display the results for D and gamma, if a structural decomposition was selected

if IRF==1 && IRFt~=1
svarinfo1=['D (structural decomposition matrix): posterior estimates'];
fprintf('%s\n',svarinfo1);
fprintf(fid,'%s\n',svarinfo1);

% recover D
D=reshape(D_estimates(:,:,hh),n,n);
% calculate the (integer) length of the largest number in D, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(D))))));
% 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(D(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

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

svarinfo2=['gamma (structural disturbances covariance matrix): posterior estimates'];
fprintf('%s\n',svarinfo2);
fprintf(fid,'%s\n',svarinfo2);

% recover gamma
gamma=reshape(gamma_estimates(:,:,hh),n,n);
% calculate the (integer) length of the largest number in D, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(gamma))))));
% 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(gamma(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

end

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

% finally, estimate and print the forecast evaluation criteria
% first note that forecast evaluation can only be conducted if it is activated, and if there is some observable data after the beginning of the forecast
if Feval==1 && Fcomp==1

% compute forecast evaluation
[RMSE MAE MAPE Ustat CRPS_estimates S1_estimates S2_estimates]=bear.panelfeval(n,p,k,beta_gibbs(:,:,hh),sigma_gibbs(:,:,hh),forecast_record(:,:,hh),forecast_estimates(:,:,hh),Fcperiods,data_endo_c(:,:,hh),data_endo_c_lags(:,:,hh),data_exo_c,const,It,Bu);
% then display and save the results
bear.panelfprint(n,endo,RMSE,MAE,MAPE,Ustat,CRPS_estimates,S1_estimates,S2_estimates,stringdates3,Fstartdate,Fcenddate,Fcperiods,fid);

% if forecast evaluation is activated but not possible, return a message to signal it
elseif Feval==1 && Fcomp==0

finfo1=['Forecast evaluation cannot be conducted.'];
fprintf('%s\n',finfo1);
fprintf(fid,'%s\n',finfo1);
finfo2=['Forecasts start in ' Fstartdate ', while observable data is available only until ' names{end,1} '.'];
fprintf('%s\n',finfo2);
fprintf(fid,'%s\n',finfo2);
finfo3=['To obtain forecast evaluation, the forecast start date must be anterior to the end of the data set.'];
fprintf('%s\n',finfo3);
fprintf(fid,'%s\n',finfo3);

% if forecast evaluation is not activated altogether, do not do anything
end

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

end

fclose(fid);

author: european-central-bank

repo: BEAR-toolbox


text:
function removed = removeIndicesInDimenion(input,dimension,indices)
% Remove the indices in a specified field in the given dimension
% USAGE:
%    removed = removeIndicesInDimenion(input, dimension, indices)
%
% INPUTS:
%
%    input:              The input matrix or array
%    dimension:          The dimension from which to remove the indices
%    indices:            The indices to remove
%
% OUTPUT:
%
%    removed:          The array/matrix with the given indices removed.
%
% NOTE:
%     Based on https://stackoverflow.com/questions/22537326/on-shape-agnostic-slicing-of-ndarrays
%
% .. Authors:
%                   - Thomas Pfau Sept 2017, adapted to merge all fields.

inputDimensions = ndims(input);
S.subs = repmat({':'},1,inputDimensions);
S.subs{dimension} = indices;
S.type = '()';
removed = subsref(input,S);

author: opencobra

repo: cobratoolbox


text:
%kvclose 'Close Boundaries of an Edge Image (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vclose.pane file
%
% Parameters:
% InputFile: i1 'Input Boundary Image', required: 'Input Boundary Image (result from vdrf, vgef, vsdef ...)'
% InputFile: i2 'Input Gradient Image', required: 'Input Gradient Image (result from vdrf, vsdef, sobel...)'
% Integer: grad 'Gradient Threshold ', default: 0: 'Gradient Threshold (used by the closing process)'
% Integer: len 'Minimum Gap Size ', default: 1: 'Minimum Gap Size In Pixels'
% OutputFile: o 'Output Image ', required: 'resulting output image or file'
%
% Example: o = kvclose({i1, i2}, {'i1','';'i2','';'grad',0;'len',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vclose - Close Boundaries of an Edge Image  (K1)
%
%  DESCRIPTION
% .I vclose
% uses a boundary image and a gradient image as input. Then it tries to close
% the boundaries by following the maximum gradient between two end points of
% two contours.
%
% The new closing part is kept if its length is at least
% "Minimum_Gap_Size" pixels. This parameter is used to prevent the result
% from being covered by small contours (edges of noisy areas...).
%
% The minimum gradient value "Gradient_Threshold" has to be set up to
% indicate to the closing process up to what value it can follow a maximum
% gradient curve to close a boundary.
%
% This Algorithm has been written by:
%
%
% Professor SERGE CASTAN, Prof CHEN and Dr ZHAO
%
%
% IRIT, CNRS, URA 1399
%
%
% 118, route de Narbonne  31062 Toulouse  FRANCE
%
%
%
%
%  EXAMPLES
% vclose -i1 moon_edge.xv -i2 moon_grad.xv -v 5 -f 0.5 -o closed_boundary.xv
%
% Will use moon_edge.xv and moon_grad.xv as the edges and respectively
% the gradient representation of an image. Then it will close the unclosed
% boundaries using a minimum gap size of 5 pixels and a gradient threshold of 0.5
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% .I vclose
% will only work with single band images. Then the first image has to be of data
% storage type VFF_TYP_1_BYTE. The second input can be anything from VFF_TYP_1_BYTE
% to VFF_TYP_FLOAT.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kvclose(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,..] = kvclose(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'grad', 0;'len', 1;'o', '__output'};
maxval={0,0,1,1,0};
minval={0,0,1,1,0};
istoggle=[0,0,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Integer','Integer','OutputFile'};
% 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 'vclose"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% check whether the sequences are of the same length.
% The sequences are aligned from begining. If some sequence is shorter, it
% will be padded with a big negative number, i.e. -1e10, in the first dimension.
%
function [data2, mask, variableLength] = ExtractVariableLengthTrajectory(data, mask)
if nargin<2
    mask = squeeze(data(1,:,:))== -1e10;
end

[D,T,N] = size(data);
if T==1 || N==1
    mask = mask';
end
if sum(sum(mask))==0    % all have equal length
    variableLength = 0;
else
    variableLength = 1;
end
for i=1:N
    data2{i} = data(:,mask(:,i)==0,i);
end
end

author: singaxiong

repo: SignalGraph


text:
classdef MimTimer < CoreBaseClass
    % MimTimer. Part of the internal framework of the TD MIM Toolkit.
    %
    %     You should not use this class within your own code. It is intended to
    %     be used internally within the framework of the TD MIM Toolkit.
    %
    %     MimTimer is used to time plugin execution
    %
    %     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 (Access = private)
        StartTicId       % ID from the tic function called at Start()
        PauseTicId       % ID from the tic function called at Resume()
        Reporting        % error/progress reporting
    end

    properties (SetAccess = private)
        SelfTime         % After Stop() is called, contains elapsed time in seconds, excluding when the timer was paused
        TotalElapsedTime % After Stop() is called, contains elapsed time in seconds, ignoring any Pause() calls
    end

    methods
        function obj = MimTimer(reporting)
            obj.Reporting = reporting;
        end

        function Start(obj)
            if ~isempty(obj.StartTicId)
                reporting.Error('MimTimer:TimerAlreadyStarted', 'The Start() method was called but the clock has already been started.');
            end
            obj.StartTicId = tic;
            obj.PauseTicId = obj.StartTicId;
        end

        function Stop(obj)
            if isempty(obj.StartTicId)
                reporting.Error('MimTimer:TimerNotStarted', 'The Stop() method was called but the clock has not been started.');
            end
            if ~isempty(obj.TotalElapsedTime)
                reporting.Error('MimTimer:ClockAlreadyStopped', 'The Stop() method was called but the clock has already been stopped.');
            end
            obj.TotalElapsedTime = toc(obj.StartTicId);
            elapsed_time_since_resume = toc(obj.PauseTicId);
            if isempty(obj.SelfTime)
                obj.SelfTime = elapsed_time_since_resume;
            else
                obj.SelfTime = obj.SelfTime + elapsed_time_since_resume;
            end
        end

        function Pause(obj)
            if isempty(obj.StartTicId)
                reporting.Error('MimTimer:TimerNotStarted', 'The Pause() method was called but the clock has not been started.');
            end
            if isempty(obj.PauseTicId)
                reporting.Error('MimTimer:TimerNotStarted', 'The Pause() method was called but the clock is already paused.');
            end
            if ~isempty(obj.TotalElapsedTime)
                reporting.Error('MimTimer:ClockAlreadyStopped', 'The Pause() method was called but the clock has already been stopped.');
            end
            elapsed_time_since_resume = toc(obj.PauseTicId);
            if isempty(obj.SelfTime)
                obj.SelfTime = elapsed_time_since_resume;
            else
                obj.SelfTime = obj.SelfTime + elapsed_time_since_resume;
            end
            obj.PauseTicId = [];
        end

        function Resume(obj)
            if isempty(obj.StartTicId)
                obj.Reporting.Error('MimTimer:TimerNotStarted', 'The Resume() method was called but the clock has not been started.');
            end
            if ~isempty(obj.TotalElapsedTime)
                obj.Reporting.Error('MimTimer:ClockAlreadyStopped', 'The Resume() method was called but the clock has already been stopped.');
            end
            if ~isempty(obj.PauseTicId)
                obj.Reporting.Error('MimTimer:NotPaused', 'The Resume() method was called but the clock has not been paused.');
            end
            obj.PauseTicId = tic;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [U,H,it] = qdwh(A,alpha,L,piv)
%QDWH   QR-based dynamically weighted Halley iteration for polar decomposition.
%   [U,H,it,res] = qdwh(A,alpha,L,PIV) computes the
%   polar decomposition A = U*H of a full rank M-by-N matrix A with
%   M >=  N.   Optional arguments: ALPHA: an estimate for norm(A,2),
%   L: a lower bound for the smallest singular value of A, and
%   PIV = 'rc' : column pivoting and row sorting,
%   PIV = 'c'           : column pivoting only,
%   PIV = ''    (default): no pivoting.
%   The third output argument IT is the number of iterations.

[m,n] = size(A);

tol1 = 10*eps/2; tol2 = 10*tol1; tol3 = tol1^(1/3);
if m == n && norm(A-A','fro')/norm(A,'fro') < tol2;
   symm = 1;
else
   symm = 0;
end

it = 0;

if m < n, error('m >= n is required.'), end

if nargin < 2 || isempty(alpha) % Estimate for largest singular value of A.
   alpha = normest(A,0.1);
end

% Scale original matrix to form X0.
U = A/alpha; Uprev = U;

if nargin < 3 || isempty(L) % Estimate for smallest singular value of U.
   Y = U; if m > n, [Q,Y] = qr(U,0); end
   smin_est =  norm(Y,1)/condest(Y);  % Actually an upper bound for smin.
   L = smin_est/sqrt(n);
end

if nargin < 4, piv = ''; end

col_piv = strfind(piv,'c');
row_sort = strfind(piv,'r');

if row_sort
   row_norms = sum(abs(U),2);
   [ignore,rind] = sort(row_norms,1,'descend');
   U = U(rind,:);
end

while norm(U-Uprev,'fro') > tol3 || it == 0 || abs(1-L) > tol1

      it = it + 1;
      Uprev = U;

      % Compute parameters L,a,b,c (second, equivalent way).
      L2 = L^2;
      dd = ( 4*(1-L2)/L2^2 )^(1/3);
      sqd = sqrt(1+dd);
      a = sqd + sqrt(8 - 4*dd + 8*(2-L2)/(L2*sqd))/2;
      a = real(a);
      b = (a-1)^2/4;
      c = a+b-1;
      % Update L.
      L = L*(a+b*L2)/(1+c*L2);

      if c > 100 % Use QR.
         B = [sqrt(c)*U; eye(n)];

          if col_piv
             [Q,R,E] = qr(B,0,'vector');
          else
             [Q,R] = qr(B,0); %E = 1:n;
          end

          Q1 = Q(1:m,:); Q2 = Q(m+1:end,:);
          U = b/c*U + (a-b/c)/sqrt(c)*Q1*Q2';

      else % Use Cholesky when U is well conditioned; faster.
          C = chol(c*(U'*U)+eye(n));
          % Utemp = (b/c)*U + (a-b/c)*(U/C)/C';
          % Next three lines are slightly faster.
          opts1.UT = true; opts1.TRANSA = true;
          opts2.UT = true; opts2.TRANSA = false;
          U = (b/c)*U + (a-b/c)*(linsolve(C,linsolve(C,U',opts1),opts2))';
    end
    if symm
       U = (U+U')/2;
    end
end
if row_sort
   U(rind,:) = U;
end

if nargout > 1
    H = U'*A; H = (H'+H)/2;
end

author: Sable

repo: mcbench-benchmarks


text:
tmptmptext = strvcat('                                  EEGLAB 1997-2010 (compiled version)', ...
                  '       ', ...
                  '       ', ...
                  'The EEGLAB electrophysiological data analysis software is an open source research', ...
                  'project developped at the Swartz Center for Computational Neuroscience, Institute', ...
                  'for Neural Computation at the University of California San Diego in collaboration', ...
                  'with The Salk Institute (La Jolla) and the CERCO laboratory (Toulouse) where part', ...
                  'of the core code was developped.', ...
                  ' ', ...
                  'EEGLAB incorporates several third party suites of functions, Fieldtrip for source ', ...
                  'localization, and BIOSIG for binary file import. A list of plugins is also included', ...
                  'in this compiled release. See http://sccn.ucsd.edu/eeglab/plugins.html for more', ...
                  'information.', ...
                  ' ', ...
                  'The EEGLAB software, its source code, and the source code of all additional third', ...
                  'party functions it includes is protected by the GNU GPL open source license.', ...
                  ' ');

warndlg2(tmptmptext);
clear tmptmptext;

author: PatternRecognition

repo: OpenBMI


text:
function test_suite=test_assert_greater_than()
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_assert_greater_than_exceptions
    assertExceptionThrown(@()assertGreaterThan([1 2],[1,2,3]),...
                          'assertGreaterThan:sizeNotEqual');
    assertExceptionThrown(@()assertGreaterThan([1,1],[0;0]),...
                          'assertGreaterThan:sizeNotEqual');
    assertExceptionThrown(@()assertGreaterThan([1],'a'), ...
                          'assertGreaterThan:classNotEqual');
    assertExceptionThrown(@()assertGreaterThan(1,1), ...
                          'assertGreaterThan:notGreaterThan');
    assertExceptionThrown(@()assertGreaterThan([2,3,3],[1,2,3]), ...
                          'assertGreaterThan:notGreaterThan');
    assertExceptionThrown(@()assertGreaterThan(nan,nan), ...
                          'assertGreaterThan:notGreaterThan');
    assertExceptionThrown(@()assertGreaterThan(nan,1), ...
                          'assertGreaterThan:notGreaterThan');
    assertExceptionThrown(@()assertGreaterThan(1,nan), ...
                          'assertGreaterThan:notGreaterThan');
    assertExceptionThrown(@()assertGreaterThan(ones(2,2,2,2),ones(2,2,2,2)*nan), ...
                          'assertGreaterThan:notGreaterThan');

    a = zeros(2,2,2,2);
    assertExceptionThrown(@()assertGreaterThan(a,ones(2,2,2,2)), ...
                          'assertGreaterThan:notGreaterThan');
    a(1,1,1,1) = 1;
    assertExceptionThrown(@()assertGreaterThan(a,ones(2,2,2,2)), ...
                          'assertGreaterThan:notGreaterThan');

function test_assert_greater_than_passes
    assertGreaterThan(1,0);
    assertGreaterThan([1,1],[0,0]);
    assertGreaterThan(ones(2,2,2,2),zeros(2,2,2,2));

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [v,psnrall] = gapdenoise( y , opt )
%GAPDENOISE Generalized alternating projection (GAP)-based denoising
%framework for compressive sensing reconstruction.
%   v=GAPDENOISE(y,opt) returns the reconstruction result v of the
%   measurements with CASSI or CACTI coding, where y is the measurement
%   matrix, opt is the parameters for the GAP-Denoise algorithm, typically
%   the denoiser applied in the framework.
%   Reference
%   [1] Y. Liu, X. Yuan, J. Suo, D.J. Brady, and Q. Dai, Rank Minimization
%       for Snapshot Compressive Imaging, IEEE Trans. Pattern Anal. Mach.
%       Intell. (TPAMI), DOI:10.1109/TPAMI.2018.2873587, 2018.
%   [2] X. Yuan, Generalized alternating projection based total variation
%       minimization for compressive sensing, in Proc. IEEE Int. Conf.
%       Image Process. (ICIP), pp. 2539-2543, 2016.
%   Code credit
%     Xin Yuan, Bell Labs, xyuan@bell-labs.com, initial version Jul 2,
%       2015.
%     Yang Liu, Tsinghua University, y-liu16@mails.tsinghua.edu.cn, last
%       update Jul 13, 2018.
%
%   See also TEST_GAPDENOISE.
if nargin<2
    opt = [];
end
% [0] default parameter configuration, to be specified
A  = @(x) M_func(x);
At = @(z) Mt_func(z);
if isfield(opt,'Mfunc'),  A  = @(x) opt.Mfunc(x);  end
if isfield(opt,'Mtfunc'), At = @(z) opt.Mtfunc(z); end
% Phisum   = ;
denoiser = 'vbm4d'; % video denoiser
% v0       = At(y); % start point (initialization of iteration)
lambda   = 0.2;     % correction coefficiency
maxiter  = 300;     % maximum number of iteration
acc      = 1;       % enable acceleration
tvweight = 0.07;    % weight for TV denoising
tviter   = 5;       % number of iteration for TV denoising
nosestim = true;    % enable noise estimation (if possible)
sigma    = 10/255;  % noise deviation
eta      = 0.5;     % coefficiency for noise estimation
flag_iqa = true;    % flag of showing image quality assessments (be sure
                    %  to turn it off for benchmark)
save_iter_image = false; % flag of saving all the images of the iteration
                         %  process (true only necessary)
                         %  require assigned directory opt.iter_image_dir

if isfield(opt,'Phisum'),     Phisum = opt.Phisum;   end
if isfield(opt,'denoiser'), denoiser = opt.denoiser; end
if isfield(opt,'v0'),             v0 = opt.v0;       end
if isfield(opt,'lambda'),   lambda   = opt.lambda;   end
if isfield(opt,'maxiter'),   maxiter = opt.maxiter;  end
if isfield(opt,'acc'),           acc = opt.acc;      end
if isfield(opt,'tvweight'), tvweight = opt.tvweight; end
if isfield(opt,'tviter'),     tviter = opt.tviter;   end
if isfield(opt,'nosestim'), nosestim = opt.nosestim; end
if isfield(opt,'sigma'),       sigma = opt.sigma;    end
if isfield(opt,'eta'),           eta = opt.eta;      end
if isfield(opt,'flag_iqa'), flag_iqa = opt.flag_iqa; end
if isfield(opt,'save_iter_image'), save_iter_image = opt.save_iter_image; end

if ~exist('v0','var') || isempty(v0)
    v0 = At(y); % start point (initialization of iteration)
end
y1 = zeros(size(y),'like',y);
psnrall = []; % return empty with no ground truth
% ssimall = []; % return empty with no ground truth
% [1] start iteration
v = v0; % initialization

k = 1; % current number of iteration
for isig = 1:length(maxiter) % extension for a series of noise levels
    nsigma = sigma(isig);
    opt.sigma = nsigma;
    for iter = 1:maxiter(isig)
        for ii = 1:1
            % [1.1] Euclidean projection
            yb = A(v);
            if acc % enable acceleration
                y1 = y1+(y-yb);
                v = v+lambda*(At((y1-yb)./Phisum)); % v=v+lambda*(At*A)^-1*At*dy
            else
                v = v+lambda*(At((y-yb)./Phisum));
            end
        end
        % [1.2] Denoising to match the video prior
        switch lower(denoiser)
            case 'tv' % TV denoising
                v = TV_denoising(v,tvweight,tviter);
            case 'vbm3d' % VBM3D denoising
                [~,v] = VBM3D(v,nsigma,0,0); % noise sigma scaled to [0,255]
            case 'vbm4d' % VBM4D denoising
                if nosestim % noise estimation enabled
                    v = vbm4d(v,-1,'lc',1,1,1,0); % -1 to enable noise estimation
                else % noise estimation disabled
                    v = vbm4d(v,nsigma,'lc',1,1,1,0); % -1 to enable noise estimation
                end
            case 'bm4d' % BM4D denoising
                if nosestim % noise estimation enabled
                    v = bm4d(v,'Gauss',0,'lc',1,0); % 0 to enable noise estimation
                else % noise estimation disabled
                    v = bm4d(v,'Gauss',nsigma,'lc',1,0); % 0 to enable noise estimation
                end
            case 'wnnm_c' % WNNM video denoising (earlier C-style strucuture version)
                v = wnnmvdenoiser(v,[],[],opt); % opt.sigma
            case 'wnnm' % WNNM video denoising (MATLAB-style matrix version)
                v = wnnm_vdenoise(v,[],opt); % opt.sigma
            otherwise
                error('Unsupported denoiser %s!',denoiser);
        end
        % % [1.3] update noise standard deviation
        % nsigma = eta*sqrt(abs(sigma^2-var(v(:)-v0(:))));
        % opt.sigma = nsigma;
        % [1.4] save and show intermediate results of psnr and ssim
        if flag_iqa && isfield(opt,'orig') && (~isempty(opt.orig))
            psnrall(k) = psnr(double(v),double(opt.orig)); % record all psnr
            % ssimall(k) = ssim(double(v),double(opt.orig)); % record all ssim
            if (mod(k,5)==0)
                fprintf('  GAP-%s iteration % 4d, sigma %.1f, PSNR %2.2f dB.\n',...
                    upper(opt.denoiser),k,nsigma*255,psnrall(k));
            end
            if save_iter_image % save all the images of the iteration process
                MAXB = opt.MAXB;
                % save each frame in the assigned directory
                nim = size(v,ndims(v)); % number of frames in the video
                for iim = 1:nim
                    im = v(:,:,iim);
                    impsnr = psnr(im,opt.orig(:,:,iim),1);
                    imssim = ssim(im,opt.orig(:,:,iim));
                    if k == 1 % mkdir for the first iteration
                        subdir = sprintf('%s/frame%02d',opt.iter_image_dir,iim); % frame-wise
                        if ~exist(subdir,'dir')
                            mkdir(subdir);
                        end
                    end
                    imwrite(uint8(im*MAXB),sprintf('%s/frame%02d/frame%02d_iter%03d_sigma%.1f_psnr%2.2f_ssim%.4f.png',...
                        opt.iter_image_dir,iim,iim,k,nsigma*MAXB,impsnr,imssim));
                end
            end
        end
        k = k+1;
    end % GAP loop [maxiter]
end % sigma loop [length(sigma)]

end

repo: DeSCI


text:
%INTIMAGE Compute integral image
%
% OUT = INTIMAGE(IM) is an integral image corresponding to IM.
%
% Integral images can be used for rapid computation of summations over
% rectangular regions.
%
% Examples::
% Create integral images for sum of pixels over rectangular regions
%        i = intimage(im);
%
% Create integral images for sum of pixel squared values over rectangular
% regions
%        i = intimage(im.^2);
%
% See also IISUM.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.

function ii = intgimage(I)

    ii = cumsum( cumsum(I)' )';

author: petercorke

repo: machinevision-toolbox-matlab