text:
%kvgef 'First Derivative Operator for Symmetric Exponential Filter (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vgef.pane file
%
% Parameters:
% InputFile: i 'Input Image', required: 'input image'
% Double: a0 ' a0 ', default: 0.45: 'Exponential Filter Parameter for GEF'
% Integer: t1 ' T1 ', default: 3: '1st Hysteresis Threshold'
% Integer: t2 ' T2 ', default: 4: '2nd Hysteresis Threshold'
% Integer: l 'Length', default: 5: 'Minimum Number of Pixels in a Segment'
% OutputFile: o 'Output Image', required: 'output image'
%
% Example: o = kvgef(i, {'i','';'a0',0.45;'t1',3;'t2',4;'l',5;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vgef - First Derivative Operator for Symmetric Exponential Filter  (K1)
%
%  DESCRIPTION
%
% This Algorithm has been written by the team of:
%
%
% Professor SERGE CASTAN
%
%
% IRIT, URA 1399
%
%
% 118, route de Narbonne  31062 Toulouse  FRANCE
%
% OPTIMAL FILTER(ISEF) FOR EDGE DETECTION
%
% PRINCIPLE :
% I. Introduction
% Edge detection is one of the most important subjects in image processing,
% which finds wide applications in the pattern recognition, the scene analysis and
% the 3-D vision, because the edges correspond in general to the important changes
% of physical or geometrical properties of objects in the scene and they are widely
% used as primitives in the pattern recognition, the image matching etc.
%
% The edges coincide, generally speaking, grey level transition, they can
% be detected by maxima of gradient or the zero-crossing of the second
% derivatives calculated by some differential operators. Because the
% differential operators are sensitive to noise, a preprocessing such as
% smoothing is often necessary to eliminate the noise. A well-known smoothing
% filter is Gaussian filter and the edges can therefore be detected by a
% Laplacian-Gaussian filter. But there is an essential difficulty of
% the Laplacian-Gaussian filter which is the contradiction between the
% smoothing effect and the precision of edge localization. To overcome this
% difficulty,  We proposed the optimal linear filter based on one step model
% (a step edge and the white noise) and the multi-edge model [9][10][11].
% This optimal smoothing filter is a symmetric exponential filter of an
% infinitely large window size and can be realized by very simple recursive
% algorithm. It  is proved that the band limited Laplacian of an input image
% filtered by this filter can be calculated from the Difference between the input
% and the output of this Recursive Filter (DRF). The edges detected by DRF
% method are less noisy and with a much better precision of localization.
%
% The maxima of gradient or zeros of the second directional derivative
% along the gradient are a natural definition of intensity edges. Zeros of the
% Laplacian are only extensively used for their computational convenience.
% However, we must stress here that the zeros crossing of the Laplacian are
% not always coincided with the maxima of gradient, for example, the zeros
% of the Laplacian are farther apart than the maxima of gradient for circularly
% symmetric patterns, this lack of localization by the Laplacian can also be
% seen in the fact that zeros of Laplacian "swing wide" of corners. Therefore,
% it had better to detect the edges by maxima of gradient or zeros of the second
% directional derivative along the gradient .
%
% So we propose two methods for edges detection, one uses maxima of Gradient
% (GEF), another uses the zeros crossing of Second directional Derivative along
% the gradient (SDEF).
%
% II. The First and the Second Directional Derivative Operators for Symmetric
% Exponential Filter
% A normalized symmetric exponential filter on 1-D can be written :
% fL(x)=C*a0*(1-a0)**|x|=f1(x)#f2(x)=C*(f1(x)+f2(x)-a0*d(x)) (1)
%
% where:
% C=1/(2-a0), # means the convolution, d(x) is dirac function.
% f1(x)={0 if x<0  and a0*(1-a0)**x   if x>= 0  (2)
% f2(x)={0 if x>0  and a0*(1-a0)**-x  if x<= 0  (2)
% we can write the first derivative operator of exponential filter:
%
% fL'(x)=f2(x)-f1(x)         (3)
% And we can obtain the normalized second derivative operator of
% exponential filter :
% fL"(x)=f1(x)+f2(x)-2*d(x)  (4)
%
% Because the exponential function is separable, we can write out 2-D
% exponential filter:
% f(x,y) = fL(x)*fL(y)                                            (5)
%
% From the equations (18),(19) and (20),  the first and the second
% directional derivative operators for symmetric exponential filter can be
% written like this:
%
%
% fx(x,y) = fL(y)*(f2(x)-f1(x))          (6)
%
%
% fy(x,y) = fL(x)*(f2(y)-f1(y))          (7)
%
%
% fxx(x,y) = fL(y)*(f1(x)+f2(x)-2*d(x))  (8)
%
%
% fyy(x,y) = fL(x)*(f1(y)+f2(y)-2*d(y))  (9)
%
% III. The Recursive Algorithm for realizing the These Directional Derivative
% Operators of Symmetric Exponential Filter
%
% The exponential filter is an IIR filter corresponding to an infinite
% window size, so we should realize the functions f1(x) and f2(x) (see formula
% (2)) by a recursive algorithm.
%
% Supposing I(x,y) is the input image,
% I1(x,y)=I(x,y)#f1(x) and I2(x,y)=I(x,y)#f2(x),
% we have the recursive algorithm :
% I1(x,y)=I1(x-1,y)+a0*(I(x,y)-I1(x-1,y))    (10)
% I2(x,y)=I2(x+1,y)+a0*(I(x,y)-I2(x+1,y))
%
% From the equations (1),(6),(7),(8) and (9), the band-limited first
% and second directional derivative of input image can be calculated by the
% recursive algorithm f1 and f2 as follows
% Ix(x,y)=I(x,y)#f1(y)#f2(y)#(f2(x)-f1(x))   (11)
% Ixx(x,y)=I(x,y)#f1(y)#f2(y)#(f2(x)+f1(x))-2*I(x,y)#f1(y)#f2(y)  (12)
% Iy(x,y) = I(x,y)#f1(x)#f2(x)#(f2(y)-f1(y))                        (13)
% Iyy(x,y) = I(x,y)#f1(x)#f2(x)#(f2(y)+f1(y))-2*I(x,y)#f1(x)#f2(x)  (14)
%
% With this algorithm, we can calculate at the same time the band-limited
% first and second directional derivative Ix and Ixx ( or Iy and Iyy ) of
% input image.
%
% IV. Edges Detection
% The band-limited first and second directional derivative of input image
% can be obtained by the algorithms as stated above. Using them, we can then
% realize the edge detection for an image.
%
% The maxima of gradient or zeros of the second directional derivative
% along the gradient are a natural way of characterizing and localizing intensity
% edges,  so we present here to detect the edges from the maxima of gradient or
% zeros of the second directional derivative along the gradient by using the
% differential operators of exponential filter.
%
%
% 1. Edges from the maxima of gradient
% Using the first directional derivative operator of the exponential filter,
% the two band-limited first directional derivatives Ix and Iy can be calculated,
% and the gradient vector can be therefore determined approximatively for every
% point in image. The gradient magnitude image is then non maxima suppressed in
% the gradient direction and thresholded with hysteresis, i.e. if the entire
% segment of the contour lies above a low threshold T1, and at least one of part of
% which is above a high threshold T2, that contour is output. The non maxima
% suppression scheme requires three points, one of which will be the current
% point, and the other two should be estimated of the gradient magnitude at
% points displaced from the current point by vector normal to the edge direction.
%
%
% 2. Edges from the zero-crossings of the second directional derivative along the
% gradient direction
% Because edges detected from local gradient maxima can not be a pixel
% width (less good precision of localization), we propose an other method which
% detect the edges from the zeros crossing of the second directional derivative
% along the gradient direction.
%
% We can calculate Ix,Iy,Ixx and Iyy by using the method shown in
% paragraph III, and therefore obtain approximatively the gradient vector
% and the second derivative in the gradient direction for every point in image.
% We extract at first the zero crossing of second derivative along the gradient
% direction on which the gradient magnitude must be above a low threshold, so an
% edges image is obtained. To this image, the entire segment of the contour will
% be kept, if the gradient magnitude on at least one part of this contour is above
% a high threshold.
% V. Comparison of Performance of the Filters
% Filtering is a problem of estimation from noisy signal, and edge
% detection is a problem of estimating the position of maximal local signal
% change. Up to now, many works are done for edges detection in image, and
% different filters are proposed, for example, Gaussian filter, Canny filter,
% exponential filter, Deriche filter etc...
%
% We appreciate the performance of the filters as follows :
%
%
% (1) Precision of edge localization
% According to our analysis [10], we can calculate the average
% localization error xe for Gaussian filter, Canny filter [15], Deriche
% filter [16] and the exponential filter :
%
%
% xeG = (4*(2*e*3.14)**0.5)/a
%
%
% xec = 0.81/a
%
%
% xeD = 4*exp(-1)/a = 1.47/ a
%
%
% xeE = 0
%
%
% i.e. xeG > xeD > xec > xeE = 0.
%
% So,  we can see that the exponential filter localizes edge points with
% the best precision.
%
%
% (2) Signal/Noise ratio on the edge point detected
% Because xe is the average estimation for the position of the edge point
% detected, we propose to calculate Signal/Noise ratio (Eq.(7)) at the point xe.
%
% And the signal/noise ratio for the Gaussian filter, Canny filter,
% Deriche filter and the exponential filter is :
%
%
% SNRG = 2*s*exp(-32*s)/ ((3.14)**0.5)
%
%
% SNRc = 0.39/ a
%
%
% SNRD = 0.64/ a
%
%
% SNRE =1/ a
%
%
% i.e. SNRE > SNRD > SNRc > SNRG.
%
% Then, we see that the exponential filter has the best noise eliminating
% effect among the above four.
%
%
% (3) Complexity of calculation
% For the complexity of calculation, we only tell the difference from
% exponential filter and Deriche filter [16], because they are implemented by
% recursive algorithms which have a simpler calculation.
%
% Because ISEF can be realized by first order recursive filter, the
% ISEF algorithms are much simpler than that of Deriche filter. Besides, the
% ISEF algorithms can be implemented independently to every line and every column,
% it can be easily realized by a parallel system.
%
% According to the analysis results above, the ISEF filter is superior
% to the others at the 3 principal aspects of the performance of the filter.
%
% VI. Conclusion
% The symmetric exponential filter of an infinite large window size is
% an optimal linear filter deduced from one step edge model and the multi-edge
% model, now we further prove that the symmetric exponential filter is the optimal
% edge detection filter in the criteria of the signal to noise ratio, the
% localization precision and unique maximum. Obviously, the real images will be
% still more complicated than these models, however DRF method has already provided
% good results for different type of images. The results obtained through the
% two new methods further show the superior performance of this filter.
% The theoretical analysis for the performance of the filters shows also that
% the exponential filter is superior to the other current filters.
%
% The first and second directional derivative operators can be realized by
% recursive algorithm and calculated at the same time. The new algorithms
% are therefore very simple as well as DRF algorithm, and they are also easy
% to implemented in a parallel way.
%
%
%
%  EXAMPLES
% vgef -i cross.xv -o output.xv
%
% will compute a edge extraction by using the first directional
% derivative operator for symmetric exponential filter on
% image cross.xv and write the result in output.xv
%
% For the GEF edge extraction the filter parameter value is generally
% a0 between 0.40 and 0.9 (0.45  [default] is a good value)
%
% The hysteresis Thresholds T1 and T2 can be set to 12 and 15.
%
% The minimum Length of segment is up to the user. It is set by default to 10.
%
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% works only with single band byte images.
%
% Note that this routine was converted directly from Khoros 1.0.  Due to
% the mixing of unsigned versus signed char, and the difference in the rules
% that resolve the mixing of these data types between ANSI C and the original
% K&R C that the code was written in, the results of running this program
% on an image will not be the same as with vdrf in Khoros 1.0.  There is
% really no way to resolve which version is more correct, but the ANSI C version
% should at least be more consistent across machine architectures.
%
%  REFERENCES
%
% [1] W.K. PRATT, Digital Image Processing, New York, 1978.
%
% [2] J. PREWITT, Object Enhancement and Extraction,Picture Processing and
% Psychopictories, Etd. by B. Lipkin and A. Rosenfeld, New York, pp.75-149, 1970.
%
% [3] M. HUECKEL, An Operator Which Locates Edges in Digitized Pictures.
% J.A.C.M., Vol. 18, pp 113-125, 1971.
%
% [4] R.O. DUDA and P.E. HART, Pattern Classification and Scene Analysis.
% Wiley, New York, 1973.
%
% [5] R. HARALICK,Edge and Region Analysis for Digital Image Data. C.G.I.P.,
% Vol. 12, pp 60-73, 1980.
%
% [6] R.HARALICK and L.WATSON, A Facet Model for Image Data. C.G.I.P., Vol. 15,
% pp 113-129, 1981.
%
% [8] D. MARR and E.C. HILDRETH, Theory of Edge Detection. Proc. R. Soc.
% Lond. B, Vol. 207, pp 187-217, 1980.
%
% [9] J. SHEN and S. CASTAN, Un nouvel algorithme de detection de contours,
% proceedings of 5th Conf. on P.R.&.A.I. (in French), Grenoble, 1985.
%
% [10] J. SHEN and S. CASTAN, An Optimal Linear Operator for Edge Detection.
% Proc. CVPR'86,  Miami,1986.
%
% [11] J. SHEN and S. CASTAN, Edge Detection Based on Multi-Edge Models.Proc.
% SPIE'87, Cannes, 1987.
%
% [12] J. SHEN and S. CASTAN, Further Results on DRF Method for Edge Detection.
% 9th I.C.P.R., ROME, 1988.
%
% [13] V.TORRE and T.A.POGGIO, On Edge Detection IEEE Transaction on Pattern
% Analysis and    Machine Intelligence, Vol. Pami-8, N 2, March 1986.
%
% [14] J.S.CHEN and G.MEDIONI, Detection, Localization, and Estimation of Edges.
% IEEE  Transaction on Pattern  Analysis and Machine Intelligence, Vol. 11, N 2,
% February  1989.
%
% [15] J.F.CANNY, Finding Edges And Lines in Images.
% MIT Technical Report N 720, 1983.
%
% [16] R. DERICHE, Optimal Edge Detection Using Recursive Filtering.
% In proc. First International Conference on Computer Vision,
% London, June 8-12 1987.
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kvgef(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,..] = kvgef(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'a0', 0.45;'t1', 3;'t2', 4;'l', 5;'o', '__output'};
maxval={0,1,255,255,100000,0};
minval={0,0,0,0,0,0};
istoggle=[0,1,1,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','Double','Integer','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 'vgef"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function []=panel4irfdisp_mean(N,n,Units,endo,irf_estimates_m,strshocks_estimates_m,IRFperiods,IRFt,stringdates1,T,decimaldates1,pref)

% IRFs
if pref.plot
% plot the figure
irf=figure('Tag','BEARresults');
set(irf,'Color',[0.9 0.9 0.9]);
   if IRFt==1
   set(irf,'name',['impulse response functions (no structural identification)']);
   elseif IRFt==2
   set(irf,'name',['impulse response functions (structural identification by Choleski ordering)']);
   elseif IRFt==3
   set(irf,'name',['impulse response functions (structural identification by triangular factorisation)']);
   elseif IRFt==4
   set(irf,'name',['impulse response functions (structural identification by sign restrictions)']);
   end
% initiate the count
count=0;
% loop over endogenous variables
for ii=1:n
   % loop over shocks
   for jj=1:n
   % increment count
   count=count+1;
   % then plot
   subplot(n,n,count)
   temp=irf_estimates_m{ii,jj};
   hold on
   Xpatch=[(1:IRFperiods) (IRFperiods:-1:1)];
   Ypatch=[temp(1,:) fliplr(temp(3,:))];
   IRFpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
   set(IRFpatch,'facealpha',0.5);
   set(IRFpatch,'edgecolor','none');
   plot(temp(2,:),'Color',[0.4 0.4 1],'LineWidth',2);
   plot([1,IRFperiods],[0 0],'k--');
   hold off
   minband=min(temp(1,:));
   maxband=max(temp(3,:));
   space=maxband-minband;
   Ymin=minband-0.2*space;
   Ymax=maxband+0.2*space;
   set(gca,'XLim',[1 IRFperiods],'YLim',[Ymin Ymax],'FontName','Times New Roman');
      % top labels
      if count<=n
      title(endo{count,1},'FontWeight','normal');
      end
      % side labels
      if jj==1
      ylabel(endo{ii,1},'FontWeight','normal');
      end
   end
end
% top supertitle
ax=axes('Units','Normal','Position',[.11 .075 .85 .88],'Visible','off');
set(get(ax,'Title'),'Visible','on')
title('Shock:','FontSize',11,'FontName','Times New Roman','FontWeight','normal');
% side supertitle
ylabel('Response of:','FontSize',12,'FontName','Times New Roman','FontWeight','normal');
set(get(ax,'Ylabel'),'Visible','on');
% save on Excel
% create the cell that will be saved on excel
irfcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},IRFperiods+3,1);
horzspace=repmat({''},3,5*n);
% loop over endogenous variables
for ii=1:n
% initiate the cell of results
endocell={};
   % loop over shocks
   for jj=1:n
   % create a header
   header=[{['response of ' endo{ii,1} ' to ' endo{jj,1} ' shocks']} {''} {''} {''};{''} {''} {''} {''};{''} {'lower bound'} {'median'} {'upper bound'}];
   % complete the cell
   tempcell=[[header;num2cell((1:IRFperiods)') num2cell((irf_estimates_m{ii,jj})')] vertspace];
   % concatenate to the previous parts of unitcell
   endocell=[endocell tempcell];
   end
% concatenate to the previous parts of irfcell
irfcell=[irfcell;horzspace;endocell];
end
% trim
irfcell=irfcell(4:end,1:end-1);
% write in excel
warning on MATLAB:xlswrite:AddSheet
xlswrite([datapath '\results.xlsx'],irfcell,'IRF','B2');
warning off MATLAB:xlswrite:AddSheet

% structural shocks

% generate plot, if there is any structural decomposition
if IRFt~=1
% plot
strshocks=figure('Tag','BEARresults');
set(strshocks,'Color',[0.9 0.9 0.9]);
set(strshocks,'name','structural shocks')
% initiate the count
count=0;
   % loop over units
   for ii=1:N
      % loop over endogenous variables
      for jj=1:n
      % increment count
      count=count+1;
      % then plot
      subplot(N,n,count)
      hold on
      Xpatch=[decimaldates1' fliplr(decimaldates1')];
      Ypatch=[strshocks_estimates_m{jj,1,ii}(1,:) fliplr(strshocks_estimates_m{jj,1,ii}(3,:))];
      Fpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
      set(Fpatch,'facealpha',0.6);
      set(Fpatch,'edgecolor','none');
      strs=plot(decimaldates1,strshocks_estimates_m{jj,1,ii}(2,:),'Color',[0.4 0.4 1],'LineWidth',2);
      plot([decimaldates1(1,1),decimaldates1(end,1)],[0 0],'k--');
      hold off
      minband=min(strshocks_estimates_m{jj,1,ii}(1,:));
      maxband=max(strshocks_estimates_m{jj,1,ii}(3,:));
      space=maxband-minband;
      Ymin=minband-0.2*space;
      Ymax=maxband+0.2*space;
      set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'YLim',[Ymin,Ymax],'FontName','Times New Roman');
         % top labels
         if count<=n
         title(endo{count,1},'FontWeight','normal');
         end
         % side labels
         if jj==1
         ylabel(Units{ii,1},'FontWeight','normal');
         end
      end
   end
end
end
   % save on Excel
% create the cell that will be saved on excel
strshockcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);
horzspace=repmat({''},3,5*n);
   % loop over units
   for ii=1:N
   % initiate the cell of results
   unitcell={};
      % loop over endogenous variables (horizontal dimension)
      for jj=1:n
      % create a header
      header=[{[Units{ii,1} ': ' endo{jj,1}]} {''} {''} {''};{''} {''} {''} {''};{''} {'lower bound'} {'median'} {'upper bound'}];
      % complete the cell
      endocell=[[header;stringdates1 num2cell(strshocks_estimates_m{jj,:,ii}')] vertspace];
      % concatenate to the previous parts of unitcell
      unitcell=[unitcell endocell];
      end
   % concatenate to the previous parts of afcell
   strshockcell=[strshockcell;horzspace;unitcell];
   end

   % trim
strshockcell=strshockcell(4:end,1:end-1);
% write in excel
if pref.results==1
bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),strshockcell,'shocks','B2');
end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [b1, b2] = isManifoldMesh(varargin)
%ISMANIFOLDMESH Check whether the input mesh may be considered as manifold.
%
%   B = isManifoldMesh(V, F)
%   B = isManifoldMesh(V, E, F)
%   Checks if the specified mesh is a manifold. When mesh is a manifold,
%   all edges are connected to either 2 or 1 faces.
%
%   [B, HASBORDER] = isManifoldMesh(V, E, F)
%   Also checks whether the mesh contains border faces. Border faces
%   contains at least one edge which is ajacent to only one face.
%
%   Example
%     [V, F] = createOctahedron;
%     isManifoldMesh(V, F)
%     ans =
%       logical
%        1
%
%   See also
%     meshes3d, ensureManifoldMesh, trimMesh

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2019-01-31,    using Matlab 9.5.0.944444 (R2018b)
% Copyright 2019 INRA - Cepia Software Platform.

vertices = varargin{1};
faces = varargin{2};

% compute edge to vertex array
if nargin == 3
    edges = faces;
    faces = varargin{3};
else
    % compute edge to vertex array
    edges = meshEdges(faces);
end

% compute face to edge indices array
% as a nFaces-by-3 array (each face connected to exactly three edges)
faceEdgeInds = meshFaceEdges(vertices, edges, faces);

% compute number of faces incident each edge
edgeFaces = trimeshEdgeFaces(faces);
edgeFaceDegrees = sum(edgeFaces > 0, 2);

% for each face, concatenate the face degree of each edge
faceEdgeDegrees = zeros(size(faces, 1), 3);
for iFace = 1:size(faces, 1)
    edgeInds = faceEdgeInds{iFace};
    faceEdgeDegrees(iFace, :) = edgeFaceDegrees(edgeInds);
end

regFaces = sum(ismember(faceEdgeDegrees, [1 2]), 2) == 3;
innerFaces = sum(faceEdgeDegrees == 2, 2) == 3;
borderFaces = regFaces & ~innerFaces;

% check if mesh is manifold: all faces are either regular or border
b1 = all(regFaces);

% check if some faces are border
b2 = any(borderFaces);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Test file for chebtech/mrdivide.m

function pass = test_mrdivide(pref)

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

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

% Random number to use as a scalar constant.
alpha = -0.194758928283640 + 0.075474485412665i;

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    %%
    % Check division of a chebtech object by a numeric array.

    f_op = @(x) [sin(x) cos(x)];
    f = testclass.make(f_op, [], pref);
    pass(n, 1) = isnan(f / 0);

    g = f / alpha;
    g_exact = @(x) f_op(x) ./ alpha;
    pass(n, 2) = norm(feval(g, x) - g_exact(x), inf) < ...
        10*max(vscale(g)*eps);

    % A "least-squares" case where the solution is obvious.
    I = eye(2);
    g = f / I;
    err = g*I - f;
    pass(n, 3) = max(max(abs(feval(err, x)))) < 10*max(vscale(g)*eps);

    % A less trivial least-squares case for which we still know the answer.
    A = [1 1];
    g = f / A;
    g_exact = @(x) (sin(x) + cos(x))/2;
    pass(n, 4) = norm(feval(g, x) - g_exact(x), inf) < ...
        10*max(vscale(g)*eps);

    %%
    % Check division of a numeric array by a chebtech object.

    f = testclass.make(@(x) sin(x));
    g = alpha / f;
    pass(n, 5) = abs(innerProduct(f, g) - alpha) < 10*max(vscale(g)*eps);

    f = testclass.make(@(x) [sin(2*pi*x) cos(2*pi*x)]);
    g = [1 1]/f;
    g_exact = @(x) (sin(2*pi*x) + cos(2*pi*x));
    pass(n, 6) = norm(feval(g, x) - g_exact(x), inf) < ...
        1e2*max(vscale(g)*eps);

    %%
    % Check error conditions.

    % Catch dimension mismatch errors.
    try
        g = f / [1 2 3];
        pass(n, 7) = false;
    catch ME
        pass(n, 7) = strcmp(ME.identifier, 'CHEBFUN:CHEBTECH:mrdivide:size');
    end

    % Can't do f/g if both f and g are chebtech objects.
    try
        f = testclass.make(@(x) sin(x));
        g = testclass.make(@(x) cos(x));
        h = f / g;
        pass(n, 8) = false;
    catch ME
        pass(n, 8) = strcmp(ME.identifier, ...
            'CHEBFUN:CHEBTECH:mrdivide:chebtechDivChebtech');
    end

    % Can't call mrdivide on a chebtech and a non-chebtech or non-double
    % object.
    try
        f = testclass.make(@(x) sin(x));
        g = f / true;
        pass(n, 9) = false;
    catch ME
        pass(n, 9) = strcmp(ME.identifier, 'CHEBFUN:CHEBTECH:mrdivide:badArg');
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function catOBJs(input_files, output_file)
  % Concatenate meshes obj files into one mesh
  %
  % catOBJs(input_files, output_file)
  %
  % Inputs:
  %   input_files  cell of paths to input files
  % Output:
  %   output_file  path to output OBJ
  %
  [V,F] = load_mesh(input_files{1});
  for i = 2:numel(input_files)
    [Vi,Fi] = load_mesh(input_files{i});
    [V, F] = cat_meshes(V,F,Vi,Fi);
  end
  writeOBJ(output_file,V,F);
end

author: alecjacobson

repo: gptoolbox


text:
function g = real(f)
%REAL Real part of a BALLFUN.
%   REAL(F) is the real part of the BALLFUN F.
%
% See also IMAG.

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

g = compose( f, @real );

end


text:
function y = cot(x)
%COT          Implements  cot(x)  for intervals
%
%   y = cot(x)
%
%interval standard function implementation
%

% written  12/30/98     S.M. Rump
% modified 09/13/99     S.M. Rump  complex allowed, sparse input, NaN input,
%                                  major revision, improved accuracy
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 12/04/05     S.M. Rump  extreme values for approximate part
% modified 09/06/07     S.M. Rump  approximate std fcts removed
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/20/08     S.M. Rump  check for zero omitted
% modified 10/18/08     S.M. Rump  StdFctsException ignore/NaN
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if x.complex
    if issparse(x.mid)
      x.mid = full(x.mid);
      x.rad = full(x.rad);
    end
    y = cos(x)./sin(x);
    if rndold
      setround(rndold)
    end
    return
  end

  INTLAB_STDFCTS_PI = getappdata(0,'INTLAB_STDFCTS_PI');
  if issparse(x.inf)
    x.inf = full(x.inf);
    x.sup = full(x.sup);
  end

  y = x;

  % transform x.inf and x.sup mod pi/2
  [ xinfinf , xinfsup , Sinf ] = modpi2(x.inf);
  [ xsupinf , xsupsup , Ssup ] = modpi2(x.sup);

  % indices with result +/- inf
  setround(1)
  delta = x.sup-x.inf;
  Tinf = Sinf + 2;
  Tinf(Tinf>7) = Tinf(Tinf>7) - 8;
  Tsup = Ssup + 2;
  Tsup(Tsup>7) = Tsup(Tsup>7) - 8;
  indexinf = ( delta >= INTLAB_STDFCTS_PI.PIINF ) | ...
    ( floor(Tinf/4) ~= floor(Tsup/4) ) | ...
    ( ( x.inf<=0 ) & ( x.sup>=0 ) );
  Sinf(Sinf>3) = Sinf(Sinf>3) - 4;
  Ssup(Ssup>3) = Ssup(Ssup>3) - 4;

  % transformation of input arguments by modpi2:
  %   [ xinf,xsup,s ] = modpi2(y)  ==>  0 <= s <= 7, 0 <= x <= pi/4 and
  %   x = -pi/2 + y + s*pi/4 + 2k*pi        for s even
  %   x =       - y + (s-1)*pi/4 + 2k*pi    for s odd

  y.inf(:) = -inf;
  y.sup(:) = inf;

  % treat non-infinity intervals
  Sinf(indexinf) = -1;
  Ssup(indexinf) = -1;

  % save warning status
  wng = warning;
  warning off

  % treat infimum
  index = ( Ssup==0 );
  if any(index(:))
    y.inf(index) = - tan_pos(xsupsup(index),1);
  end
  index = ( Ssup==1 );
  if any(index(:))
    y.inf(index) = 1 ./ ( - tan_pos(xsupinf(index),-1) );
  end
  index = ( Ssup==2 );
  if any(index(:))
    res = tan_pos(xsupsup(index),1);
    setround(-1)
    y.inf(index) = 1 ./ res;
  end
  index = ( Ssup==3 );
  if any(index(:))
    y.inf(index) = tan_pos(xsupinf(index),-1);
  end

  % treat supremum
  index = ( Sinf==0 );
  if any(index(:))
    y.sup(index) = - tan_pos(xinfinf(index),-1);
  end
  index = ( Sinf==1 );
  if any(index(:))
    y.sup(index) = 1 ./ ( - tan_pos(xinfsup(index),1) );
  end
  index = ( Sinf==2 );
  if any(index(:))
    res = tan_pos(xinfinf(index),-1);
    setround(1)
    y.sup(index) = 1 ./ res;
  end
  index = ( Sinf==3 );
  if any(index(:))
    y.sup(index) = tan_pos(xinfsup(index),1);
  end

  % restore warning status
  warning(wng);
  setround(0)                        % set rounding to nearest

  index = isnan(x.inf);
  if any(index(:))
    y.inf(index) = NaN;
    y.sup(index) = NaN;
  end
  INTLAB_STDFCTS_EXCPTN = getappdata(0,'INTLAB_STDFCTS_EXCPTN');
  if INTLAB_STDFCTS_EXCPTN==3      % ignore input out of range (ignore-mode)
    index = ( x.inf==0 ) & ( x.sup==0 );
    if ~isempty(find(index))              % completely exceptional arguments to NaN
      setappdata(0,'INTLAB_STDFCTS_EXCPTN_',1);
      y.inf(index) = NaN;
      y.sup(index) = NaN;
    end
  elseif INTLAB_STDFCTS_EXCPTN==2  % any input out of range to NaN (NaN-mode)
    index = ( x.inf<=0 ) & ( 0<=x.sup );
    if ~isempty(find(index))              % exceptional arguments to NaN
      setappdata(0,'INTLAB_STDFCTS_EXCPTN_',1);
      y.inf(index) = NaN;
      y.sup(index) = NaN;
    end
  end

  setround(rndold)

repo: OpenMAS


text:
function [mi3] = m32mi3(m3)
% Convert volume from cubic meters to cubic miles.
% Chad Greene 2012
mi3 = m3*2.3991275858e-10;

author: Sable

repo: mcbench-benchmarks


text:
function runcilta() % autogenerated function wrapper
    %% This script runcias find first a maximum of all z-vlaues over all
    % frames and afterward produced a movie and start the movieviever
    % Last edit: Stefan Wiemer 11/94
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun(mfilename('fullpath'));
    think

    j = 0;
    it = 20
    iwl = iwl3/days(ZG.bin_days);
    [len, ncu] = size(cumuall);
    len = len -2;
    step = len/nustep;

    ma = [];
    mi = [];

    wai = waitbar(0,'Please wait...')
    set(wai,'Color',[0.8 0.8 0.8],'NumberTitle','off',...
        'Name','fin_max -Percent done','pos',[ZG.welcome_pos 300 80]);

    pause(0.1)
    % find maximu for scaling
    %
    for it = 1:step:len-iwl;
        fin_malt
        waitbar(it/len)
    end   % for i
    close(wai)

    ZG.maxc = max(max(ma));
    ZG.minc = min(min(mi));

    %set up movie axes
    %
    figure
    tmp = gcf

    cin_lta
    axes(has)
    fs_m = get(gcf,'pos');

    m = moviein(length(1:step:len-iwl));

    wai = waitbar(0,'Please wait...')
    set(wai,'Color',[0.8 0.8 0.8],'NumberTitle','off',...
        'Name','Movie -Percent done','pos',[ZG.welcome_pos 300 80]);
    pause(0.1)

    for it = 1:step:len-iwl;
        j = j+1;
        cin_lta
        axes(has)
        m(:,j) = getframe(has);
        figure(wai);
        waitbar(it/len)
    end   % for i

    close(tmp)
    close(wai)

    % save movie
    %
    clear newmatfile

    [newmatfile, newpath] = uiputfile('*.mat', 'Save As');

    if length(newpath)  > 1
        save([newpath newmatfile file1])
        showmovi
    else
        showmovi
    end

end

author: CelsoReyes


text:
function [ g , mu ] = gsp_design_translates(G, g0,N )
%GSP_DESIGN_TRANSLATES Create a filterbank by uniformly translating a window
%   Usage: g = gsp_design_translates( G, g0, Ntrans );
%
%   Inputs parameters:
%       G       : Graph structure
%       g0      : Mother window (anonymous function)
%       N       : Number of translate
%
%   Outputs parameters:
%       g       : filterbank
%       mu      : Centers of the filters
%
%   This function construct a filter bank of *N* uniformly translated
%   filter from the mother filter *g0*.
%

% Author : Nathanael Perraudin
% Date: 6 January 2016

if isstruct(G)
    if ~isfield(G,'lmax')
        if param.verbose
            fprintf('GSP_DESIGN_TRANSLATE has to compute lmax \n')
        end
        G = gsp_estimate_lmax(G);
    end
   lmax = G.lmax;
else
   lmax = G;
end

mu = linspace(0,lmax,N);

g = cell(length(mu),1);

for ii = 1:length(mu)
    g{ii} = @(x) g0(x-mu(ii));
end

end

repo: gspbox


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function [rho,drho] = rhoSpline(Tc,Rc,minT,maxT,nT,minR,maxR,nR,doDerivative)
%
% Parzen-Window Based density estimator for gray values of Tc and Rc, spline based.
%
% Input:
%  Tc, Rc          template and reference
%  minT, maxT, nT  discreted points in gray value range where splines are located
%  minR, maxR, nR  discreted points in gray value range where splines are located
%  doSerivative         FLAG for derivative compuatation
%
% Output:
%  rho          joint density estimator
%  drho         derivative of joint density estimator
%
% see also MIspline for an example.
%==============================================================================

function [rho,drho] = rhoSpline(Tc,Rc,minT,maxT,nT,minR,maxR,nR,doDerivative)

if nargin == 0,
  help(mfilename);
  return;
end;

if ~exist('doDerivative','var'), doDerivative = 1;       end;
doDerivative = doDerivative & (nargout>1);

% prepare the output and organize input
rho  = []; drho = []; Tc = reshape(Tc,[],1); Rc = reshape(Rc,[],1);

widthT = (maxT-minT)/nT;
widthR = (maxR-minR)/nR;

Tt   = linspace(minT,maxT,nT);
Rt   = linspace(minR,maxR,nR);
rho  = zeros(length(Tt),length(Rt));
drho = spalloc(numel(rho),length(Tc),5*length(Tc));

for j=1:length(Tc), % run over all samples
  IT = find(abs(Tt-Tc(j))<widthT); % find locations of interest in Tc
  IR = find(abs(Rt-Rc(j))<widthR); % find locations of interest in Rc
  if ~isempty(IT) && ~isempty(IR),
    [KT,dKT] = spline1D(Tc(j) - Tt(IT),widthT,doDerivative);
    KR       = spline1D(Rc(j) - Rt(IR),widthR,0);
    rho(IT,IR) = rho(IT,IR) +KT*KR';
    if doDerivative,
      drhoj = zeros(size(rho));
      drhoj(IT,IR) = dKT*KR';
      drho(:,j) = drho(:,j) + sparse(drhoj(:));
    end;
  end;
end;

% normalize rho
fac = (Tt(2)-Tt(1))*(Rt(2)-Rt(1))/length(Tc);
rho  = fac*rho(:);
drho = fac*drho;

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

% The Parzen Window Kernel is a spline function
%  function [s,ds] = spline1D(x,sigma,doDerivative);
% (c) Jan Modersitzki 2008/02/12, see FAIR.
% This evaluates a spline function s
% s(x,sigma) = 0 for x\notin[-sigma,sigma], \int s(x,sigma) dx = 1

function [s,ds] = spline1D(x,sigma,doDerivative);

if ~exist('doDerivative','var'), doDerivative = 1;       end;
doDerivative = doDerivative & (nargout>1);

% map [-sigm2,sigma] \to [-2,2]
x = 4/sigma*x;

if ~doDerivative,
  s  = 0*reshape(x,[],1);
  ds = [];
  J = find(x>=-2 & x<-1);   s(J)   = (x(J)+2).^3;
  J = find(x>=-1 & x< 0);   s(J)   = -x(J).^3 - 2*(x(J)+1).^3+6*(x(J)+1);
  J = find(x>= 0 & x< 1);   s(J)   = -2*(1-x(J)).^3 + x(J).^3 - 6*x(J) + 6;
  J = find(x>= 0 & x< 1);   s(J)   = -2*(1-x(J)).^3 + x(J).^3 - 6*x(J) + 6;
  J = find(x>= 1 & x< 2);   s(J)   = (2-x(J)).^3;
else
  s  = 0*reshape(x,[],1);
  ds = s;
  J = find(x>=-2 & x<-1);   s(J)   = (x(J)+2).^3;
                            ds(J)  = 3*(x(J)+2).^2;
  J = find(x>=-1 & x< 0);   s(J)   = -x(J).^3 - 2*(x(J)+1).^3+6*(x(J)+1);
                            ds(J)  = -3*x(J).^2 - 6*(x(J)+1).^2+6;
  J = find(x>= 0 & x< 1);   s(J)   = -2*(1-x(J)).^3 + x(J).^3 - 6*x(J) + 6;
                            ds(J)  = 6*(1-x(J)).^2 + 3*x(J).^2 - 6;
  J = find(x>= 1 & x< 2);   s(J)   = (2-x(J)).^3;
                            ds(J)  = -3*(2-x(J)).^2;
end;
% scale to make integral 1
s  = 2/(3*sigma)*s;
ds = 8/(3*sigma^2)*ds;
%==============================================================================


text:
classdef MONRP < PROBLEM
% <multi> <binary> <large/none>
% The multi-objective next release problem
% m --- 100 --- Number of customers

%------------------------------- Reference --------------------------------
% Y. Zhang, M. Harman, and S. A. Mansouri, The multi-objective next release
% problem, Proceedings of the Annual Conference on Genetic and Evolutionary
% Computation, 2007, 1129-1137.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties(Access = private)
        Cost;   % Cost of each requirement
        Value;  % Value of each customer on each requirement
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            % Parameter setting
            m = obj.ParameterSet(100);
            obj.M = 2;
            if isempty(obj.D); obj.D = 100; end
            obj.encoding = 4 + zeros(1,obj.D);
            % Randomly generate costs and values
            n    = obj.D;
            file = sprintf('MONRP-n%d-m%d.mat',n,m);
            file = fullfile(fileparts(mfilename('fullpath')),file);
            if exist(file,'file') == 2
                load(file,'Cost','Value');
            else
                Cost  = randi(9,1,n);
                Value = randi([0 5],n,m);
                save(file,'Cost','Value');
            end
            obj.Cost  = Cost;
            obj.Value = Value;
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj(:,1) = sum(repmat(obj.Cost,size(PopDec,1),1).*PopDec,2);
            PopObj(:,2) = sum(obj.Value(:)) - sum(PopDec*obj.Value,2);
        end
        %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [sum(obj.Cost),sum(obj.Value(:))];
        end
        %% Display a population in the objective space
        function DrawObj(obj,Population)
            PopObj(:,1) = sum(repmat(obj.Cost,length(Population),1).*Population.decs,2);
            PopObj(:,2) = sum(Population.decs*obj.Value,2);
            Draw(PopObj,{'Cost','Satisfaction score',[]});
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [centers, labels, info] = slkmeansex(X, n, estfunctor, clsfunctor, varargin)
%SLKMEANSEX Performs Generalized K-means
%
% $ Syntax $
%   - [centers, labels] = slkmeansex(X, n, estfunctor, clsfunctor, ...)
%   - [centers, labels, info] = slkmeansex(X, n, estfunctor, clsfunctor, ...)
%
% $ Arguments $
%   - X:            the samples to be clustered
%   - n:            the number of samples
%   - estfunctor:   the functor to estimate means(centers), as follows:
%                   centers = estfunc(centers, X, K, weights, labels, ...)
%                   when input centers is empty, it performs initial
%                   estimation, otherwise, it performs updating.
%                   In addition, it should ignore the samples with
%                   labels being zeros or negative numbers.
%   - clsfunctor:   the functor to classify samples
%                   labels = clsfunc(centers, X, n, ...)
%                   it should produce 1 x n row vector.
%   - centers:      the clustered centers
%   - labels:       the labels indicating which sample belong to which center
%                   a 1 x n row vector.
%   - info:         the information on iteration process
%
% $ Description $
%   - [centers, labels] = slkmeansex(X, n, estfunctor, clsfunctor, ...)
%     is a generalized version of K-means. It actually implements an
%     iterative process to estimate centers from clustered samples and
%     re-clustered the samples according to centers.
%     You can specify the following properties:
%       - 'K':              the number of initial number of clusters
%                           (default = 3)
%       - 'init_centers':   the initial centers.
%       - 'maxiter':        the maximum number of iterations
%                           (default = 100);
%       - 'annthres':       the threshold of annealing
%                           when the sum of sample weights for a center
%                           is below annthres * the total weight, the
%                           center will be discarded. (default = 0)
%       - 'annfunc':        the function to discard a set of centers
%                           centers = annfunc(centers, inds_discard);
%       - 'weights':        the weights of the samples (default = [])
%       - 'verbose':        whether to show progress information
%                           (default = true)
%
% $ Remarks $
%   - The X and centers can be in any form that conform to the specified
%     functors.
%
%   - If init_centers is specified, K should be exactly the number of
%     initial centers.
%
%   - If annthres is 0, then no centers will be discarded even some centers
%     have no support samples in the process. The estfunctor should keep
%     those centers unchanged.
%
% $ History $
%   - Created by Dahua Lin, on Aug 28, 2006
%   - Modified by Dahua Lin, on Aug 30, 2006
%       - utilize slevalfunctor and slsharedisp
%   - Modified by Dahua Lin, on Aug 31, 2006
%       - based on slreevallearn
%

%% parse and verify input arguments

if nargin < 4
    raise_lackinput('slkmeansex', 4);
end

opts.K = 3;
opts.init_centers = [];
opts.maxiter = 100;
opts.annthres = 0;
opts.annfunc = [];
opts.weights = [];
opts.verbose = true;
opts = slparseprops(opts, varargin{:});

if opts.K > n
    error('sltoolbox:rterror', ...
        'The initial K is larger than the number of samples');
end

if opts.annthres > 0
    if isempty(opts.annfunc)
        error('sltoolbox:invalidarg', ...
            'You should specify annfunc when annthres > 0');
    end
end

w = opts.weights;
if ~isempty(w)
    if ~isequal(w, [1 n])
        error('sltoolbox:sizmismatch', ...
            'The weights should be a 1 x n row vector');
    end
end

%% Initialization

slsharedisp_attach('slkmeansex', 'show', opts.verbose);

slsharedisp('Intialize K-Means');

if isempty(opts.init_centers)
    initcinds = randsample(n, opts.K);
    labels = zeros(1, n);
    labels(initcinds) = 1:opts.K;

    K = opts.K;
    centers = slevalfunctor(estfunctor, [], X, K, w, labels);
else
    K = opts.K;
    centers = opts.init_centers;
end

slsharedisp_incindent;
slsharedisp('initial K = %d', K);
slsharedisp_decindent;

labels = slevalfunctor(clsfunctor, centers, X, n);

%% Updating

slsharedisp('Update K-Means');
slsharedisp_incindent;

km_estfunctor = {@kmeansex_est, estfunctor, opts};
km_evalfunctor = {@kmeansex_eval, clsfunctor};
km_cmpfunctor = {@kmeansex_cmp};

models = {centers, K};
data = {X, n, w};
[models, labels, info] = slreevallearn(models, labels, data, ...
    km_estfunctor, km_evalfunctor, km_cmpfunctor, ...
    'iter', {'maxiter', opts.maxiter, 'titlebreak', false}, 'isrecorded', false);

centers = models{1};

slsharedisp_decindent;
slsharedisp_detach;

%% Core functions

% models = {centers, K}
% data = {X, n, w}

function models = kmeansex_est(models, data, labels, estfunctor, opts)

X = data{1};
w = data{3};
centers = models{1};
K = models{2};

if ~isempty(centers) && opts.annthres > 0
    if isempty(w)
        w = ones(1, length(labels));
    end
    cw = sllabeledsum(w, labels, 1:K);
    wthres = opts.annthres * sum(cw) / K;
    if any(cw < wthres)
        inds_ann = find(cw < wthres);
        centers = feval(opts.annfunc, centers, inds_ann);
        K = K - length(inds_ann);

        models = {centers, K};
        return;
    end
end

centers = slevalfunctor(estfunctor, centers, X, K, w, labels);
models = {centers, K};

function labels = kmeansex_eval(models, data, labels, clsfunctor)

X = data{1};
n = data{2};
centers = models{1};

slignorevars(labels);

labels = slevalfunctor(clsfunctor, centers, X, n);

function isconverged = kmeansex_cmp(models_prev, models, labels_prev, labels)

K_prev = models_prev{2};
K = models{2};
n = length(labels);

slsharedisp_attach('kmeansex_cmp');

isconverged = false;
if K == K_prev
    nchanged = sum(labels ~= labels_prev);
    slsharedisp('K = %d: %d / %d changed', K, nchanged, n);

    if nchanged == 0
        isconverged = true;
    end
else
    slsharedisp('K = %d ==> %d', K_prev, K);
end

slsharedisp_detach();

author: lmthang


text:
function varargout=tach(varargin)
%
% [hr]=tach(recordName,annotator,N,N0,ouputSize)
%
%    Wrapper to WFDB TACH:
%         http://www.physionet.org/physiotools/wag/tach-1.htm
%
%    Reads a WFDB record and Annotation file to return:
%
%
% hr
%       Nx1 vector of doubles representing a uniformly sampled and
%       smoothed instantaneous heart rate signal. The output are samples
%       of the instantaneous heart rate signal in units of beats per minute.
%
% Required Parameters:
%
% recorName
%       String specifying the name of the record in the WFDB path or
%       in the current directory.
%
% annotator  -
%       String specifying the name of the annotation file in the WFDB path or
%       in the current directory.
%
% Optional Parameters are:
%
% N
%       A 1x1 integer specifying the sample number at which to stop reading the
%       annotation file (default read all = N).
% N0
%       A 1x1 integer specifying the sample number at which to start reading the
%       annotation file (default 1 = begining of the record).
%
% outputSize
%
%    A 1x1 integer specifying the number of output samples (ie estimated
%    heart rate intervals) such that the output 'hr' is a vector of
%    size (outputSize-1) x 1.
%
%
% Written by Ikaro Silva, 2013
% Last Modified: January 24, 2014
% Version 1.1
%
% Since 0.0.1
%
% %Example 1- Read a signal and annotaion from PhysioNet's Remote server:
%[hr]=tach('challenge/2013/set-a/a01','fqrs');
%plot(hr);grid on;hold on

%endOfHelp

[javaWfdbExec,config]=getWfdbClass('tach');

%Set default pararamter values
inputs={'recordName','annotator','N','N0','ouputSize'};
outputs={'data(:,1)'};
N=[];
N0=1;
ouputSize=[];
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

N0=num2str(N0-1); %-1 is necessary because WFDB is 0 based indexed.
wfdb_argument={'-r',recordName,'-a',annotator,'-f',['s' N0]};

if(~isempty(N))
    wfdb_argument{end+1}='-t';
    %-1 is necessary because WFDB is 0 based indexed.
    wfdb_argument{end+1}=['s' num2str(N-1)];
end

if(~isempty(ouputSize))
    wfdb_argument{end+1}='-n';
    wfdb_argument{end+1}=[num2str(ouputSize)];
end

data=javaWfdbExec.execToDoubleArray(wfdb_argument);
data=wfdbjava2mat(data);
for n=1:nargout
        eval(['varargout{n}=' outputs{n} ';']);
end

author: ikarosilva

repo: wfdb-app-toolbox


text:
function box = mergeBoxes3d(box1, box2)
%MERGEBOXES3D Merge 3D boxes, by computing their greatest extent.
%
%   BOX = mergeBoxes3d(BOX1, BOX2);
%
%   Example
%   box1 = [5 20 5 30 10 50];
%   box2 = [0 15 0 15 0 20];
%   mergeBoxes3d(box1, box2)
%   ans =
%       0 20 0 30 0 50
%
%
%   See also
%   boxes3d, drawBox3d, intersectBoxes3d
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% unify sizes of data
if size(box1,1) == 1
    box1 = repmat(box1, size(box2,1), 1);
elseif size(box2, 1) == 1
    box2 = repmat(box2, size(box1,1), 1);
elseif size(box1,1) ~= size(box2,1)
    error('Bad size for inputs');
end

% compute extreme coords
mini = min(box1(:,1:2:end), box2(:,1:2:end));
maxi = max(box1(:,2:2:end), box2(:,2:2:end));

% concatenate result into a new box structure
box = [mini(:,1) maxi(:,1) mini(:,2) maxi(:,2) mini(:,3) maxi(:,3)];


text:
function selIdx = hierarchicalSelect(st)

selID = 567; % Cerebrum, default to start

[boxList, idList] = makeBoxList(st, selID);

ud.idList = idList; ud.st = st;

% make figure
f = figure; set(f, 'KeyPressFcn', @doOK);

% selector box
ud.hBox = uicontrol(f, 'Style', 'listbox', 'String', boxList, ...
    'Callback', @updateSel, 'Value', find(idList==selID),...
    'Units', 'normalized', 'Position', [0.1 0.2 0.8 0.7],...
    'KeyPressFcn', @doOK);

titleStr = boxList{idList==selID}; titleStr = titleStr(find(titleStr~=' ', 1):end);
ud.hSelTitle = uicontrol(f, 'Style', 'text', ...
    'String', sprintf('Selected: %s', titleStr), ...
    'Units', 'normalized', 'Position', [0.1 0.9 0.8 0.1]);

ud.hCancel = uicontrol(f, 'Style', 'pushbutton', ...
    'String', 'Cancel', 'Callback', @doCancel, ...
    'Units', 'normalized', 'Position', [0.1 0.1 0.2 0.1]);

ud.hOK = uicontrol(f, 'Style', 'pushbutton', ...
    'String', 'OK', 'Callback', @doOK, ...
    'Units', 'normalized', 'Position', [0.3 0.1 0.2 0.1]);

set(f, 'UserData', ud);
drawnow;

% AP: internal matlab command didn't work for Mac users, removed
% try
%     c = matlab.ui.internal.dialog.DialogUtils.disableAllWindowsSafely();
%     uiwait(f);
%     delete(c);
% catch
%     if ishghandle(f)
%         delete(f)
%     end
% end
uiwait(f);

if ishghandle(f)
    ud = get(f, 'UserData');
    idList = ud.idList;

    if ud.hBox.Value>1
        selID = idList(get(ud.hBox, 'Value'));

        selIdx = find(st.id==selID);
    else
        selIdx = [];
    end
    delete(f)
    drawnow;
else
    selIdx = [];
end

end

function [boxList, idList] = makeBoxList(st, selID)

idList = selID;
while idList(end)~=997 % root
    idList(end+1) = st.parent_structure_id(st.id==idList(end));
end
idList = idList(end:-1:1); % this is the tree of parents down to the selected one

% make the parent string representation
for q = 1:numel(idList)
    boxList{q} = sprintf('%s%s (%s)', repmat('  ', 1, q-1), ...
        st.acronym{st.id==idList(q)}, ...
        st.safe_name{st.id==idList(q)});
end
np = numel(idList);

% now add children
idList = [idList st.id(st.parent_structure_id==selID)'];

% make the parent string representation
for q = np+1:numel(idList)
    boxList{q} = sprintf('%s%s (%s)', repmat('  ', 1, np), ...
        st.acronym{st.id==idList(q)}, ...
        st.safe_name{st.id==idList(q)});
end

end

function updateSel(src, ~)

f = get(src, 'Parent');
ud = get(f, 'UserData');
st = ud.st; idList = ud.idList;

selID = idList(get(src, 'Value'));

[boxList, idList] = makeBoxList(st, selID);

ud.idList = idList;
set(f, 'UserData', ud);
set(src, 'String', boxList, 'Value', find(idList==selID));

titleStr = boxList{idList==selID}; titleStr = titleStr(find(titleStr~=' ', 1):end);
set(ud.hSelTitle, 'String', sprintf('Selected: %s', titleStr));

end

% OK callback
function doOK(~, ~)
    uiresume(gcbf);
end

% Cancel callback
function doCancel(~, ~)
ud = get(gcbf, 'UserData');
ud.hBox.Value = 1;
uiresume(gcbf);
end

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function [Xsol, Ssol] = generalized_eigenvalue_computation(A, B, p)
% Returns orthonormal basis of the dominant invariant p-subspace of B^-1 A.
%
% function [Xsol, Ssol] = generalized_eigenvalue_computation(A, B, p)
%
% Input: A is a real, symmetric matrix of size nxn,
%        B is a symmetric positive definite matrix, same size as A, and
%        an integer p < n.
% Output: Xsol: a real, orthonormal matrix X of size nxp such that
%         trace(X'*A*X) is maximized, subject to X'*B*X being identity.
%         That is, the columns of X form an orthonormal basis of a dominant
%         subspace of dimension p of B^(-1)*A. These are thus eigenvectors
%         associated with the largest eigenvalues of B^(-1)*A  (in no
%         particular order). Sign is important: 2 is deemed a larger
%         eigenvalue than -5.
%         Ssol: the eigenvalues associated with the eigenvectors Xsol.
%
% We intend to solve the homogeneous system A*X = B*X*S,
% where S is a diagonal matrix of dominant eigenvalues of B^-1 A.
%
%
% The optimization is performed on the generalized Grassmann manifold,
% since only the space spanned by the columns of X matters.
%
% The optimization problem that we are solving here is
% maximize trace(X'*A*X) subject to X'*B*X = eye(p).
% Consequently, the solutions remain invariant to transformation
% X --> XO, where O is a p-by-p orthogonal matrix. The search space, in
% essence, is set of equivalence classes
% [X] = {XO : X'*B*X = I and O is orthogonal matrix}. This space is called
% the generalized Grassmann manifold.
%
% See also dominant_invariant_subspace nonlinear_eigenspace

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Bamdev Mishra, June 30, 2015.
% Contributors:
% Change log:

    % Generate some random data to test the function
    if ~exist('A', 'var') || isempty(A)
        n = 128;
        A = randn(n);
        A = (A+A')/2;
    end
    if ~exist('B', 'var') || isempty(B)
        n = size(A, 1);
        e = ones(n, 1);
        B = spdiags([-e 2*e -e], -1:1, n, n); % Symmetric positive definite
    end

    if ~exist('p', 'var') || isempty(p)
        p = 3;
    end

    % Make sure the input matrix is square and symmetric
    n = size(A, 1);
	assert(isreal(A), 'A must be real.')
    assert(size(A, 2) == n, 'A must be square.');
    assert(norm(A-A', 'fro') < n*eps, 'A must be symmetric.');
	assert(p <= n, 'p must be smaller than n.');

    % Define the cost and its derivatives on the generalized
    % Grassmann manifold, i.e., the column space of all X such that
    % X'*B*X is identity.
    gGr = grassmanngeneralizedfactory(n, p, B);

    problem.M = gGr;
    problem.cost  = @(X)    -trace(X'*A*X);
    problem.egrad = @(X)    -2*(A*X); % Only Euclidean gradient needed.
    problem.ehess = @(X, H) -2*(A*H); % Only Euclidean Hessian needed.

    % Execute some checks on the derivatives for early debugging.
    % These things can be commented out of course.
    % checkgradient(problem);
    % pause;
    % checkhessian(problem);
    % pause;

    % Issue a call to a solver. A random initial guess will be chosen and
    % default options are selected except for the ones we specify here.
    options.Delta_bar = 8*sqrt(p);
    [Xsol, costXsol, info] = trustregions(problem, [], options); %#ok<NASGU,ASGLU>

    % To extract the eigenvalues, solve the small p-by-p symmetric
    % eigenvalue problem.
    Ssol = eig(Xsol'*(A*Xsol));

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [H,F] = make_bank_DFT(p,N)

% make_bank_DFT     Generate DFT Filter Bank with N Subbands
%
% Arguments:
% p                 Prototype filter (impulse response)
% N                 Number of subbands
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

p = p(:);                            % Make column
flen = max(size(p));                 % Length of prototype filter
H = zeros(flen,N/2+1);               % Assume N is an even number
n = (0:flen-1)';

k = 0; H(:,k+1) = p;                 % i = 0
k = N/2; H(:,k+1) = p.*((-1).^n);    % i = N/2+1

for k = 1:(N/2-1)
    H(:,k+1) = p.*exp(j*2*pi/N*k*n); % Complex modulation
end

D = H(:,2:N/2);
H = [H, conj(fliplr(D))];            % Complex-conjugate pair

% Synthesis filters are the same as the analysis filters

F = H;

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function output = condenseFields(data, regexp, nullData)
% "condenseFields"
%   Search a datastructure (data) recursively for instances of a fieldname
%   expression (regexp).  If the fields matching the regexp are encountered,
%   concatenate them together, creating an array, and replace the parent of
%   this list with the array.  New empty fields are filled with nullData if
%   it is specified, else they are [].
%   See dissimilarInsert for details on the concatenation.
%
%   JRA 10/28/03
%
% Warning:  This function must be used carefully, it can result in dataloss
% if you do not know beforehand how your data is structured.
%
% For DICOM plans with Item_1, Item_2, ... use 'Item_\d' as the regexp.
%
% Usage:
%   function output = condenseFields(data, regexp, nullData)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('nullData')
    nullData = [];
end
dataClass = class(data);

switch dataClass

    case 'struct'
		myFields = fields(data);
		for i=1:length(myFields)
            for j=1:length(data)
                if isstruct(data(j).(myFields{i}))
                    data(j).(myFields{i}) = condenseFields(data(j).(myFields{i}), regexp, nullData);
                end
            end
		end
		matches = regexpi(myFields, regexp);
		if length(matches) == 1;
            matches = {matches};
		end
		index = 1;
		for i=1:length(matches)
            if matches{i} == 1
                if ~exist('temp')
                    temp(index) = data.(myFields{i});
                else
                    temp = dissimilarInsert(temp, data.(myFields{i}), index, nullData);
                end
                index = index+1;
            end
		end
		if exist('temp')
            output = temp;
		else
            output = data;
		end

    case 'cell'
        numCells = length(data);
        for i=1:numCells
            data{i} = condenseFields(data{i}, regexp, nullData);
        end
        output = data;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function grid_print ( element_order, element_num, element_node )

%*****************************************************************************80
%
%% GRID_PRINT prints the elements that form a grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), the nodes that form
%    each element.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  GRID_PRINT: Element -> Node table.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Element order = %d\n', element_order );
  fprintf ( 1, '  Number of elements = %d\n', element_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  #   ' );
  for i = 1 : element_order
    fprintf ( 1, '%3d', i );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );

  for element = 1 : element_num
    fprintf ( 1, '  %3d   ', element );
    for i = 1 : element_order
      fprintf ( 1, '%3d', element_node(i,element) );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%data=Ka^ind+M;
%data must start from zero (K+M, Ka+M,....)
%data boun is the first and final ind values to be used (can be 0:end)

%ctyp: the method used to compute expoenential
%ctyp==0 => log based least square (all values must be postivie)
%ctyp==1 => normal equation (can be used also for negative values)

function [dmar_par exp_par]=expo_fit(dboun, data, ctyp, show_res, cval)
[ndim ndat]=size(data);
if (ctyp==2)
    nbas=size(cval,2);
else
    nbas=1;
end
a=zeros(ndim,nbas);
M=zeros(ndim,1);
K=zeros(ndim,nbas);

for i=1:ndim
    lb=dboun(i,1);
    ub=dboun(i,2);

    if (ctyp==1)
        %%normal equations based computation
        sr_a=data(i,(lb+1):ub)*data(i,(lb+1):ub)';
        sr_b=data(i,(lb+1):ub)*data(i,(lb+2):ub+1)';
        a_apx=sr_b/sr_a;
        dind=lb:ub;
    elseif (ctyp==0)
        %%apply a log based least square to compute a
        ind=(lb:ub);
        pind=find(data(i,ind+1)>eps*100);
        dind=ind(pind);

        H=[ones(length(dind),1) dind'];
        Y=log(data(i,dind+1))';
        x=inv(H'*H)*H'*Y;
        K_apx0=exp(x(1));
        a_apx=exp(x(2));
    elseif (ctyp==2)    %exponential bases are provided
        a_apx=cval(i,:);
        dind=lb:ub;
    end

    %optimal K & M given suboptimal "a"
    H=zeros(length(dind),nbas+1);
    for k=1:nbas
        H(:,k)=a_apx(k).^dind;
    end
    H(:,nbas+1)=1;

    apx_sol=inv(H'*H)*H'*data(i,dind+1)';
    K_apx=apx_sol(1:nbas)';
    M_apx=apx_sol(nbas+1);

    a(i,:)=a_apx;
    K(i,:)=K_apx;
    M(i)=M_apx;

    %compare the results
    if (show_res==1)
        figure;
        %compare the results
        ind=0:(length(data(i,:))-1);
        plot(ind,data(i,:));grid;
        hold on
        val_tot=0;
        for k=1:nbas
            val_th=K(i,k)*a(i,k).^ind+M(i);
            val_tot=val_tot+val_th-M(i);
            plot(ind,val_th,'r');
        end
        plot(ind,val_tot+M(i),'r');
    end
end

%exponential param
exp_par=[a K M];

%discrete markov param x'=ax+bu
b=(K.*(ones(ndim,1)-a.^2)).^0.5;
dmar_par=[a b];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     %optimal solution given suboptimal "a"
%     vr_a=a_apx.^dind;
%     sr_a=sum(vr_a);
%     sr_b=data(i,dind+1)*vr_a';
%     sr_c=sum(vr_a.^2);
%     sr_d=sum(data(i,dind+1));
%     sr_e=length(dind);
%
%     sr_f=sr_a^2-sr_e*sr_c;
%     M_apx=(sr_a*sr_b-sr_c*sr_d)/sr_f;
%     K_apx=(sr_a*sr_d-sr_b*sr_e)/sr_f;

%     %apply curve fit
%     %fit function
%     myfun2=@(p,k) p(1)+p(2)*p(3).^k;
%     opt_param=lsqcurvefit(myfun2,[M_apx, K_apx, a_apx],ind,data(i,:));
%     a(i)=opt_param(3);
%     K(i)=opt_param(2);
%     M(i)=opt_param(1);

repo: ekfmonoslam


text:
function lvmResultsClick(modelType, dataSet, number, dataType, varargin)

% LVMRESULTSCLICK Load a results file and visualise them with clicks
% FORMAT
% DESC loads results of a latent variable model and visualises them.
% ARG modelType : the type of model ran on the data set.
% ARG dataSet : the name of the data set to load.
% ARG number : the number of the run used.
% ARG dataType : the type of data to visualise.
% ARG arg1, arg2, arg3 ... : additional arguments to be passed to the
% lvmVisualise command.
%
% SEEALSO : lvmLoadResult, lvmVisualise
%
% COPYRIGHT : Neil D. Lawrence, 2008, 2009

% MLTOOLS

[model, lbls] = lvmLoadResult(modelType, dataSet, number);

% Visualise the results
switch size(model.X, 2)
 case 2
  lvmClickVisualise(model, lbls, [dataType 'Visualise'], [dataType 'Modify'], ...
                 varargin{:});

 otherwise
  lvmClickVisualise(model, lbls, [dataType 'Visualise'], [dataType 'Modify'], ...
                 varargin{:});
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function c = tapas_logrt_linear_binary_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the linear log-reaction time response model according to as
% developed with Louise Marshall and Sven Bestmann
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The Gaussian noise observation model assumes that responses have a Gaussian distribution around
% the inferred mean of the relevant state. The only parameter of the model is the noise variance
% (NOT standard deviation) zeta.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2016 Christoph Mathys, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'Linear log-reaction time for binary models';

% Sufficient statistics of Gaussian parameter priors
%
% Beta_0
c.be0mu = log(500);
c.be0sa = 4;

% Beta_1
c.be1mu = 0;
c.be1sa = 4;

% Beta_2
c.be2mu = 0;
c.be2sa = 4;

% Beta_3
c.be3mu = 0;
c.be3sa = 4;

% Beta_4
c.be4mu = 0;
c.be4sa = 4;

% Zeta
c.logzemu = log(log(20));
c.logzesa = log(2);

% Gather prior settings in vectors
c.priormus = [
    c.be0mu,...
    c.be1mu,...
    c.be2mu,...
    c.be3mu,...
    c.be4mu,...
    c.logzemu,...
         ];

c.priorsas = [
    c.be0sa,...
    c.be1sa,...
    c.be2sa,...
    c.be3sa,...
    c.be4sa,...
    c.logzesa,...
         ];

% Model filehandle
c.obs_fun = @tapas_logrt_linear_binary;

% Handle to function that transforms observation parameters to their native space
% from the space they are estimated in
c.transp_obs_fun = @tapas_logrt_linear_binary_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function fgplvmPrintPlot(model, lbls, capName, experimentNo)

% FGPLVMPRINTPLOT Print latent space for learnt model.
% FORMAT
% DESC prints a latent space repsresentation for an FGPLVM model.
% ARG model : the model to use for plotting the latent space.
% ARG lbls : any lables that are available for plotting.
% ARG capName : the name of the saved plots.
% ARG experimentNo : the experiment number to assign to the files.
%
% SEEALSO : lvmScatterPlot
%
% COPYRIGHT : Neil D. Lawrence, 2006

% FGPLVM

lvmScatterPlot(model, lbls);
fileName = ['dem' capName num2str(experimentNo)];
print('-depsc', ['../tex/diagrams/' fileName])
print('-deps', ['../tex/diagrams/' fileName 'NoColour'])

% make smaller for PNG plot.
pos = get(gcf, 'paperposition')
origpos = pos;
pos(3) = pos(3)/2;
pos(4) = pos(4)/2;
set(gcf, 'paperposition', pos);
fontsize = get(gca, 'fontsize');
set(gca, 'fontsize', fontsize/2);
lineWidth = get(gca, 'lineWidth');
set(gca, 'lineWidth', lineWidth*2);
print('-dpng', ['../html/' fileName])
set(gcf, 'paperposition', origpos);

figure
clf
ax = axes('position', [0.05 0.05 0.9 0.9]);
hold on
if ~isempty(lbls)
  lvmTwoDPlot(model.X, lbls, getSymbols(size(lbls, 2)));
else
  lvmTwoDPlot(model.X, lbls);
end
xLim = [min(model.X(:, 1)) max(model.X(:, 1))]*1.1;
yLim = [min(model.X(:, 2)) max(model.X(:, 2))]*1.1;
set(ax, 'xLim', xLim);
set(ax, 'yLim', yLim);

set(ax, 'fontname', 'arial');
set(ax, 'fontsize', 20);
print('-depsc', ['../tex/diagrams/' fileName 'NoGray'])
print('-deps', ['../tex/diagrams/' fileName 'NoGrayNoColour'])


text:
function [bbox3M, planC] = crop_for_constrictor(planC,paramS,varargin)
% Custom crop function for segmentation model
% AI 10/04/19

%% Limits of bounding box using chewing structures
indexS = planC{end};
strC = {planC{indexS.structures}.structureName};
bboxIdx = [];

%% Return bounding box if it exists
if isfield(paramS,'saveStrToPlanCFlag') && paramS.saveStrToPlanCFlag
    if isfield(paramS,'outStrName')
        outStrName = paramS.outStrName;
    else
        outStrName = 'crop_for_constrictor';
    end
    bboxIdx = getMatchingIndex(outStrName,strC,'EXACT');
end
if ~isempty(bboxIdx)
    bbox3M = getStrMask(bboxIdx,planC);
else
    bboxName = paramS.structureName.cropStructure;
    idx1 = getMatchingIndex(bboxName,strC,'EXACT');
    scanNum = getStructureAssociatedScan(idx1,planC);
    [mask3M, planC]  = getStrMask(idx1,planC);
    if sum(mask3M(:))>0
        [minr,~,minc,maxc,mins,~] = compute_boundingbox(mask3M);
    else
        minr = 1;
        minc = 1;
        maxc = size(mask3M,2);
        mins = 1;
    end
    %% Limits of bbox around larynx
    larynxStrName = paramS.structureName.larynx;
    idx2 = getMatchingIndex(larynxStrName,strC,'EXACT');
    [mask3M, planC]  = getStrMask(idx2,planC);
    if sum(mask3M(:))>0
        [~,maxr,~,~,~,maxs] = compute_boundingbox(mask3M);
    else
        maxr = size(mask3M,1);
        maxs = size(mask3M,3);
    end

    %% Get bounding box for constrictors
    tol_r = 30;
    tol_s = 15;
    maxr = min(maxr+tol_r,size(mask3M,1));
    maxs = min(maxs+tol_s,size(mask3M,3));

    bbox3M = false(size(getScanArray(scanNum,planC)));
    bbox3M(minr:maxr,minc:maxc,mins:maxs) = true;
end

end


text:
function result = applyToNonNan(data,func)
% Apply the given function columnwise to all non NaN values in the given
% data.
%
% USAGE:
%
%    result = applyToNonNan(data,func)
%
% INPUTS:
%    data:          Matrix of data (individuals x variables)
%    func:          function handle that takes an vector of data and
%                   computes single value result.
%
% OUTPUTS:
%
%    result:        A Vector of with dimensions (size(data,2) x 1) where
%                   the supplied function was applied to all columns
%                   ignoring NaN values.

result = zeros(1,size(data,2));

for i = 1:size(data,2)
    result(i) = func(data(~isnan(data(:,i)),i));
end

author: opencobra

repo: cobratoolbox


text:
function varargout = mean2(varargin)
%MEAN2   Mean of a CHEBFUN2.
%   V = MEAN2(F) returns the mean of a CHEBFUN2:
%
%   V = 1/A*integral2( f )
%
% 	where the A is the area of the domain of F.
%
% See also MEAN, STD2.

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

[varargout{1:nargout}] = mean2@separableApprox(varargin{:});

end


text:
function y = psdscale(ud,x,K,transp)
% y = psdscale(ud,x,K [,transp])
%
% PSDSCALE  Computes length lenud (=sum(K.s.^2)) vector y.
%   !transp (default) then y[k] = vec(Ldk' * Xk * Ldk)
%   transp == 1 then y[k] = vec(Udk' * Xk * Udk)
%   Uses pivot ordering ud.perm if available and nonempty.
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also scaleK, factorK.

% This file is part of SeDuMi 1.3 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

%The function is quite dirty as the type and dimension of the inputs may
%change, that is why we have so many subcases.

Ks = K.s;
if isempty(Ks)
    y = [];
    return
end
if nargin < 4
    transp = false;
end
Kq = Ks .* Ks;
nr = K.rsdpN;
nc = length(Ks);
N  = sum(Kq) + sum(Kq(nr+1:end));
y  = zeros(N,1);
xi = length(x) - N;
yi = 0;
ui = 0;
if isstruct(ud)
    perm = ud.perm;
    if isempty(perm)
        prep = false;
        postp = false;
    else
        prep = ~transp;
        postp = transp;
        pi = 0;
    end
    ud = ud.u;
else
    prep  = false;
    postp = false;
end
for i = 1 : nc
    ki = Ks(i);
    qi = Kq(i);
    TT = ud(ui+1:ui+qi); ui=ui+qi;
    if i > nr
        TT = TT + 1i*ud(ui+1:ui+qi); ui=ui+qi;
    end
    TT = reshape(TT,ki,ki);
    if transp
        TT = triu(TT);
    else
        TT = tril(TT);
    end
    XX = x(xi+1:xi+qi); xi=xi+qi;
    if i > nr
        XX = XX + 1i*x(xi+1:xi+qi); xi=xi+qi;
    end
    XX = reshape(XX,ki,ki);
    if prep
        PP = perm(pi+1:pi+ki); pi=pi+ki;
        if any(diff(PP)~=1)
            XX = XX(PP,PP);
        end
    end
    if nnz(XX) < 0.1 * qi
        XX = sparse(XX);
    end
    XX = TT' * XX * TT;
    if postp
        PP = perm(pi+1:pi+ki); pi=pi+ki;
        if any(diff(PP)~=1)
            XX(PP,PP) = XX;
        end
    end
    y(yi+1:yi+qi) = real(XX);
    yi = yi+qi;
    if i > nr
        XX = imag(XX);
        % Needed, otherwise psdfactor() will sometimes fail.
        XX(1:ki+1:end) = 0;
        y(yi+1:yi+qi) = XX;
        yi = yi+qi;
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
classdef prtPreProcClass < prtPreProc
    % prtPreProcClass
    %
    % prtPreProcClass is an abstract class and cannot be instantiated.
    %
    %   See Also: prtPreProc, prtPreProcPca, prtPreProcPls,
    %   prtPreProcHistEq, prtPreProcZeroMeanColumns, prtPreProcLda,
    %   prtPreProcZeroMeanRows, prtPreProcLogDisc, prtPreProcZmuv,
    %   prtPreProcMinMaxRows

    methods

        function obj = prtPreProcClass()
            obj.classTrain = 'prtDataSetClass';

            obj.isSupervised = true; % Overload this here.
        end
    end

    methods (Access = protected, Hidden = true)
        function ClassObj = preTrainProcessing(ClassObj, DataSet)
            % Overload preTrainProcessing() so that we can determine mary
            % output status
            assert(DataSet.isLabeled & DataSet.nClasses > 1,'The prtDataSetClass input to the train() method of a prtPreProcClass must have non-empty targets and have more than one class.');

            ClassObj = preTrainProcessing@prtAction(ClassObj,DataSet);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
classdef MimMarkerLast < MimGuiPlugin
    % MimMarkerLast. Gui Plugin for switching transparency of zero
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimMarkerLast is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     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 = 'Last marker'
        SelectedText = 'Last marker'
        ToolTip = 'Go to slice containing the last marker'
        Category = 'Slice navigation'
        Visibility = 'Dataset'
        Mode = 'Markers'

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

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.GotoLastMarker;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && isequal(gui_app.GetCurrentModeName, MimModes.MarkerMode);
        end

        function is_selected = IsSelected(gui_app)
            is_selected = false;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function PlotPolarityAroundATimeByDensity(aedat, time, numTimePoints, numDensities, minProportionOfPixels, maxProportionOfPixels, contrast, transposeVar, flipVertical, flipHorizontal)

%{
Takes 'aedat' - a data structure containing an imported .aedat file,
as created by ImportAedat, and creates a series of green/red plots of
polarity data.
Plots in Y have different densities (i.e. proportions of pixels). These
densities are distributed exponentially between the min and max parameters
Plots in X have different centre times. The numTimePoints param must be
odd, and the centre plot in X is the one centred at the parameter 'time'.

For numPlots > 1, the proportionOfPixels used to include events around that
time point is varied logarithmically from minProp.. to maxProp.
Defaults are 0.1 and 1.
flipVertical is assumed true, so that y=0 is considered the top of the image.
%}

%% Parameters

if ~exist('time', 'var')
	time = (aedat.info.firstTimeStamp + aedat.info.lastTimeStamp) / 2;
else
    time = time * 1e6;
end

if ~exist('numDensities', 'var')
	numDensities = 3;
end

if ~exist('numTimePoints', 'var')
	numTimePoints = 3;
end
if mod(numTimePoints,2) == 0
	numTimePoints = numTimePoints + 1;
end

if ~exist('distributeBy', 'var')
	distributeBy = 'time';
end

if ~exist('minProportionOfPixels', 'var') ...
        || (exist('minProportionOIfPixels', 'var') ...
            && minProportionOfPixels == 0)
    minProportionOfPixels = 0.001;
end
if ~exist('maxProportionOfPixels', 'var') ...
        || (exist('maxProportionOIfPixels', 'var') ...
            && maxProportionOfPixels == 0)
    maxProportionOfPixels = 0.1;
end

% The 'contrast' for display of events, as used in jAER.
if ~exist('contrast', 'var')
    contrast = 3;
end

%% Unpack

% break out the timeStamps for readability
timeStamp   = aedat.data.polarity.timeStamp;
x           = aedat.data.polarity.x;
y           = aedat.data.polarity.y;
polarity    = aedat.data.polarity.polarity;
numEvents   = aedat.data.polarity.numEvents;

%% Produce plots

% Find eventIndex nearest to timePoint
eventIndex = find(timeStamp >= time, 1, 'first');
if isempty(eventIndex)
    eventIndex = numEvents;
end

numPlots = numDensities * numTimePoints;

if numDensities > 1
    % distribute the
    logMin = log(minProportionOfPixels);
    logMax = log(maxProportionOfPixels);
    logStep = (logMax - logMin) / (numPlots - 1);
    proportionsOfPixels = exp(logMin : logStep : logMax);
    numPixelsInArray = aedat.info.deviceAddressSpace(1) * aedat.info.deviceAddressSpace(2);
    numPixelsToSelectEachWay = ceil(numPixelsInArray * proportionsOfPixels / 2);

    figure
else
    numPixelsToSelectEachWay = minProportionOfPixels; % Arbitrrary choice to use the minimum
end
for densityIndex = 1 : numDensities
    % First do the central timePoint
    firstIndex = max(1, eventIndex - numPixelsToSelectEachWay(plotIndex));
    lastIndex = min(numEvents, eventIndex + numPixelsToSelectEachWay(plotIndex));
    selectedLogical = [false(firstIndex - 1, 1); ...
					true(lastIndex - firstIndex + 1, 1); ...
					false(numEvents - lastIndex, 1)];
    eventsForFrame = struct;
    eventsForFrame.x        = x(selectedLogical);
    eventsForFrame.y        = y(selectedLogical);
    eventsForFrame.polarity = polarity(selectedLogical);
    frame = FrameFromEvents(eventsForFrame, contrast, aedat.info.deviceAddressSpace);
	if exist('transpose', 'var') && transposeVar
        frame = frame';
    end
    image(frame - 1);
    colormap(redgreencmap(contrast * 2 + 1))
	axis equal tight
	if ~exist('flipVertical', 'var') || flipVertical
		set(gca, 'YDir', 'reverse')
    end
    if exist('flipHorizontal', 'var') && flipHorizontal
		set(gca, 'XDir', 'reverse')
    end
	title(['Proportion of pixels: ' num2str(proportionsOfPixels(plotIndex))])

    % Then work outwards from the central time point
    lastIndexUp = lastIndex;
    firstIndexDown = firstIndex;
    for timePointsIndex = (numTimePoints - 1) / 2;
        subplot(numPlotsY, numPlotsX, plotIndex);
        hold all

        firstIndex = max(1, eventIndex - numPixelsToSelectEachWay(plotIndex));
        lastIndex = min(numEvents, eventIndex + numPixelsToSelectEachWay(plotIndex));
        selectedLogical = [false(firstIndex - 1, 1); ...
                        true(lastIndex - firstIndex + 1, 1); ...
                        false(numEvents - lastIndex, 1)];
        eventsForFrame = struct;
        eventsForFrame.x        = x(selectedLogical);
        eventsForFrame.y        = y(selectedLogical);
        eventsForFrame.polarity = polarity(selectedLogical);
        frame = FrameFromEvents(eventsForFrame, contrast, aedat.info.deviceAddressSpace);
        if exist('transpose', 'var') && transposeVar
            frame = frame';
        end
        image(frame - 1);
        colormap(redgreencmap(contrast * 2 + 1))
        axis equal tight
        if ~exist('flipVertical', 'var') || flipVertical
            set(gca, 'YDir', 'reverse')
        end
        if exist('flipHorizontal', 'var') && flipHorizontal
            set(gca, 'XDir', 'reverse')
        end
        title(['Proportion of pixels: ' num2str(proportionsOfPixels(plotIndex))])
    end
end

author: panpanfei


text:
function KM = computeKernelMatrix(P1, P2, kType, kParam)
%COMPUTEKERNELMATRIX computes the kernel matrix between two sets of patterns
%
%   KM = computeKernelMatrix(P1, P2, kType, kParam) computes kernel matrix
%   KM between sets of patterns P1 and P2. kType is the type of the kernel
%   function (Gauss, Linear, Quickrbf, Poly, or Sigmoid). KPARAM is the
%   kernel parameter, for instance the width of the Gaussian kernel.
%
% Function inputs:
%   patterns1, patterns2   - Two matrixes of patterns (train and test, or train and train)
%   kernelType             - Kernel function choice. Can be  Gauss, Linear, Quickrbf, Poly, or Sigmoid.
%   kernelMatrix           - The kernel matrix
%
%   This file is part of ORCA: https://github.com/ayrna/orca
%   Original authors: Pedro Antonio Gutiérrez, María Pérez Ortiz, Javier Sánchez Monedero
%   Citation: If you use this code, please cite the associated paper http://www.uco.es/grupos/ayrna/orreview
%   Copyright:
%       This software is released under the The GNU General Public License v3.0 licence
%       available at http://www.gnu.org/licenses/gpl-3.0.html
%
Nf1             = size(P1, 2);
Nf2             = size(P2, 2);
KM   = zeros(Nf1, Nf2);

switch upper(kType)
    case {'GAUSS','GAUSSIAN','RBF'}

        for i = 1:Nf2
            KM(:,i)  = exp(-kParam*sum((P1-P2(:,i)*ones(1,Nf1)).^2)');

        end

    case {'QUICKRBF'}
        d = pdistalt(P1,P2,'euclidean');
        KM = exp(-kParam*d.^2);

    case {'POLYNOMIAL', 'POLY', 'LINEAR'}
        if strcmpi(kType, 'LINEAR')
            kParam = 1;
            bias = 0;
        else
            bias = 1;
        end
        KM = ((P1' * P2 + bias)/size(P1,1)).^kParam;

    case 'SIGMOID'

        if (length(kParam) ~= 2)
            error('This kernel needs two parameters to operate!')
        end

        KM   = zeros(Nf1, Nf2);
        for i = 1:Nf2
            KM(:,i)  = tanh(P1'*P2(:,i)*kParam(1)+kParam(2));
        end

    otherwise
        error('Unknown kernel. Can be Gauss, Linear, Quickrbf, Poly, or Sigmoid.')

end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
function [forecast_record]=forecasttv2(data_endo_a,data_exo_p,It,Bu,beta_gibbs,omega_gibbs,F_gibbs,phi_gibbs,L_gibbs,gamma,sbar,Fstartlocation,Fperiods,n,p,k,q,const)

% create first the cell storing the results
forecast_record=cell(n,1);

% other preliminary tasks: generate the matrix of predicted exogenous variables
% if the constant has been retained, augment the matrices of exogenous with a column of ones:
if const==1
data_exo_p=[ones(Fperiods,1) data_exo_p];
% if no constant was included, do nothing
else
end

% then start simulations
% repeat the process a number of times equal to the number of simulations retained from Gibbs sampling
for ii=1:It-Bu

% compute the reduced matrix Y
Y=data_endo_a(end-p+1:end,:);

% step 3: draw beta, omega and sigma and F from their posterior distributions
% draw beta
beta=beta_gibbs{Fstartlocation-1,1}(:,ii);
% draw omega
omega=omega_gibbs(:,ii);
% create a choleski of omega, the variance matrix for the law of motion
cholomega=sparse(diag(omega));
% draw F from its posterior distribution
F=sparse(F_gibbs(:,:,ii));
% step 4: draw phi from its posterior
phi=phi_gibbs(ii,:)';
% also, compute the pre-sample value of lambda, the stochastic volatility process
lambda=L_gibbs(Fstartlocation-1,:,ii)';

   % then generate forecasts recursively
   % for each iteration ii, repeat the process for periods T+1 to T+h
   for jj=1:Fperiods

   % update beta
   beta=beta+cholomega*randn(q,1);
   % reshape it to obtain B
   B=reshape(beta,k,n);

   % use the function lagx to obtain the matrix temp
   temp=bear.lagx(Y,p-1);
   % define the reduced regressor matrix X
   % if no exogenous variable is present at all in the model (neither constant nor other exogenous), define X only from the endogenous variables
      if isempty(data_exo_p)==1
      X=[temp(end,:)];
      % if there are exogenous vaiables, concatenate them next to the endogenous
      else
      X=[temp(end,:) data_exo_p(jj,:)];
      end

   % update lambda_t and obtain Lambda_t
   % loop over variables
      for kk=1:n
      lambda(kk,1)=gamma*lambda(kk,1)+phi(kk,1)^0.5*randn;
      end
   % obtain Lambda_t
   Lambda=sparse(diag(sbar.*exp(lambda)));

   % recover sigma_t and draw the residuals
   sigma=full(F*Lambda*F');
   % draw the vector of residuals
   res=bear.trns(chol(bear.nspd(sigma),'Lower')*randn(n,1));

   % obtain predicted value for T+jj
   yp=X*B+res;

   % concatenate the transpose of yp to the top of Y
   Y=[Y;yp];

   % step 8: repeat until values are obtained for T+h
   end

   % record the results from current iteration in the cell forecast_record
   % loop over variables
   for kk=1:n
   % consider column kk of matrix Y and trim the p initial values: what remains is the predicted values for the period T+1 to T+h, for variable kk
   temp1=Y(p+1:end,kk);
   % record these values in the corresponding matrix of forecast_record
   forecast_record{kk,1}(ii,:)=temp1';
   end

% step 9: repeat until It-Bu iterations are obtained
end

author: european-central-bank

repo: BEAR-toolbox


text:
function COFMEn = COFMEn(data, m, r,fs)
%******************************************************
% $ This function is usded for calcualting the cofficient of fuzzy measure entropy (COFMEn) for physiological signal time sequence.
% $ Using COFMEn aims to improve the performance of COSEn.
%
% $ Reference: Liu C Y, Li K, Zhao L N, Liu F, Zheng D C, Liu C C and Liu S
% T. Analysis of heart rate variability using fuzzy measure entropy.
% Computers in Biology and Medicine, 2013, 43(2): 100-108
%
% $ Variable declaration:
% data is RR time series
% m is embedding dimension (usually m=1)
% r is threshold value (usually r=30 ms)
% local threshold r_l=0.2, global threshold r_g=0.2,
% local weight of sequence segments' similarity n_l=2
% global weight of sequence segments' similarity n_g=2
%
% $ Author:  Chengyu Liu (bestlcy@sdu.edu.cn)
%           Institute of Biomedical Engineering,
%           Shandong University
% $Last updated:  2015.10.15
%
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information

% data=[108,108,109,109,108,108,110,107,109,109,109,109,108,109,108,108,108,108,107,108,108,108,108,107,109,107,108,107,108,108,108,108,108,108,108,109,108,109,108,108,109,108,108,108,109,107,108,107,107,108,108,108,109,107,108,108,108,108,108,108]*4;
% data=[109,108,108,109,109,109,108,109,108,108,109,108,108,109,108,108,108,109,108,109,109,109,108,108,109,108,108,109,109,109]*4;
% m=1;
% r=30;

N=length(data);
for i=1:N-m
    x1(i,:)=data(i:i+m-1);
    x2(i,:)=data(i:i+m);
end

ratio=0.4;
if N>20
%     Thr=ratio*N;
    Thr=5;
else
    Thr=5;
end

Min_numerator=0;
kk=0;
while Min_numerator<Thr
    [Min_numerator,r]=ReGet_min_numerator(x1,x2,N,m,r,fs);
    kk=kk+1;
end
if kk==1 || kk>20
   if N<20
        r=r-1;
    else
        r=r-1000/fs;
    end
end

if Min_numerator==N-m
    while Min_numerator>=Thr
        [Min_numerator,r]=ReGet_min_numerator2(x1,x2,N,m,r,fs);
    end
    if r<0
        r=r+2;
    else
        r=r+1;
    end
end

x=data;
r_l   = r;
r_g   = r;
n_l   = 2;
n_g   = 2;
%% m=2
D_l   = zeros(N-m,1); % initiate the mean local distance vector
D_g   = zeros(N-m,1); % initiate the mean global distance vector
for i = 1:N-m
    distance_l = zeros(N-m,1); % initiate the local distance vector for the ith vector
    distance_g = zeros(N-m,1); % initiate the global distance vector for the ith vector
    for j = 1:N-m
        if j==i
            d_l = 0;
            d_g = 0;
        else
        d_l = max(abs(x(i:i+m-1)-mean(x(i:i+m-1))-x(j:j+m-1)+mean(x(j:j+m-1))));
        d_g = max(abs(x(i:i+m-1)-x(j:j+m-1)));
        end
        distance_l(j) = exp(-(d_l.^n_l/r_l));
        distance_g(j) = exp(-(d_g.^n_g/r_g));
    end
    D_l(i) = sum(distance_l)/(N-m-1);
    D_g(i) = sum(distance_g)/(N-m-1);
end
Bm_l = mean(D_l);
Bm_g = mean(D_g);

%% m=m+1
m=m+1;
D_l   = zeros(N-m,1);
D_g   = zeros(N-m,1);
for i = 1:N-m
    distance_l = zeros(N-m,1);
    distance_g = zeros(N-m,1);
    for j = 1:N-m
        if j==i
            d_l = 0;
            d_g = 0;
        else
        d_l = max(abs(x(i:i+m-1)-mean(x(i:i+m-1))-x(j:j+m-1)+mean(x(j:j+m-1))));
        d_g = max(abs(x(i:i+m-1)-x(j:j+m-1)));
        end
        distance_l(j) = exp(-(d_l.^n_l/r_l));
        distance_g(j) = exp(-(d_g.^n_g/r_g));
    end
    D_l(i) = sum(distance_l)/(N-m-1);
    D_g(i) = sum(distance_g)/(N-m-1);
end
Am_l = mean(D_l);
Am_g = mean(D_g);

%% Calculate local and global fuzzy measure entropy
FuzzyLMEn = -log(Am_l/Bm_l);  % local fuzzy measure entropy
FuzzyGMEn = -log(Am_g/Bm_g);  % global fuzzy measure entropy
%% Generate fuzzy measure entropy
COFMEn  = FuzzyLMEn+FuzzyGMEn+2*log(2*r/1000)-2*log(mean(data)/1000);

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function assocScansV = getDoseAssociatedScan(dosesV, planC)
% "getDoseAssociatedScan"
%   Returns a vector with the corresponding associated scan for each
%   dose number passed in structsV.
%
%   If planC is not specified, the global planC is used.
%
% DK 07/19/2006
%
%Usage:
%   [assocScansV, relStructNum] = getDoseAssociatedScan(dosesV, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC')
    global planC
end
indexS = planC{end};

assocScansV = [];

for i = 1:length(dosesV)
    assocScanTmp = getAssociatedScan(planC{indexS.dose}(dosesV(i)).assocScanUID, planC);
    if ~isempty(assocScanTmp)
        assocScansV(i) = assocScanTmp;
    end
end