text:
function tfrqview(tfr,sig,t,method,p1,p2,p3,p4,p5);
%TFRQVIEW Quick visualization of time-frequency representations.
%       TFRQVIEW(TFR,SIG,T,METHOD,P1,P2,P3,P4,P5) allows a quick
%       visualization of a time-frequency representation.
%
%       TFR     : time-frequency representation (MxN).
%       SIG     : signal in time. If unavailable, put sig=[] as input
%                 parameter.                    (default : []).
%       T       : time instants                 (default : 1:N).
%       METHOD  : name of chosen representation (default : 'TYPE1').
%                 See the TFR* files for authorized names.
%                 TYPE1 : the representation TFR goes in normalized
%                       frequency from -0.5 to 0.5 ;
%                 TYPE2 : the representation TFR goes in normalized
%                       frequency from 0 to 0.5.
%       P1...P5 : optional parameters of the representation : run the
%                 file TFRPARAM(METHOD) to know the meaning of P1..P5
%                 for your method.
%
%	When you use the 'save' option in the main menu, you save all your
%	variables as well as two strings, TfrQView and TfrView, in a mat
%	file. If you load this file and do eval(TfrQView), you will restart
%	the display session under tfrqview ; if you do eval(TfrView), you
%	will obtain the exact layout of the screen you had when clicking on
%	the 'save' button.
%
%       Example :
%        sig=fmsin(128); tfr=tfrwv(sig);
%        tfrqview(tfr,sig,1:128,'tfrwv');
%
%       See also TFRVIEW, TFRSAVE, TFRPARAM.

%       F. Auger, September 1994, July 1995
%       O. Lemoine, Oct 1995, May-July 1996.
%       Copyright (c) 1996 by CNRS (France).
%
%  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 St, Fifth Floor, Boston, MA  02110-1301  USA

comp=computer;

% Tests on the input arguments
if nargin<1,
 error('At least one parameter required');
end
[tfrrow,tfrcol]=size(tfr);
if nargin==1,
 sig=[]; t=1:tfrcol; method='type1';
elseif nargin==2,
 t=1:tfrcol; method='type1';
elseif nargin==3,
 method='type1';
end
[trow,tcol] = size(t);
if (trow~=1),
 error('T must only have one row');
end;
if (tfrcol~=tcol),
 error('T must have as much elements as tfr has columns');
end;
method=upper(method);
[Nsig,Ncol]=size(sig);
if Ncol>2,
 error('SIG must have less than two columns');
end

% Computation of Nf2, the number of interresting points in frequency
method=upper(method);
if strcmp(method,'TFRWV'   ) | strcmp(method,'TFRPWV'  ) | ...
   strcmp(method,'TFRSPWV' ) | strcmp(method,'TFRCW'   ) | ...
   strcmp(method,'TFRZAM'  ) | strcmp(method,'TFRBJ'   ) | ...
   strcmp(method,'TFRBUD'  ) | strcmp(method,'TFRGRD'  ) | ...
   strcmp(method,'TFRRSPWV') | strcmp(method,'TFRRPWV' ) | ...
   strcmp(method,'TFRRIDB' ) | strcmp(method,'TFRRIDH' ) | ...
   strcmp(method,'TFRRIDT' ) | strcmp(method,'TFRASPW' ) | ...
   strcmp(method,'TFRDFLA' ) | strcmp(method,'TFRSPAW' ) | ...
   strcmp(method,'TFRRIDBN') | strcmp(method,'TFRUNTER') | ...
   strcmp(method,'TFRBERT' ) | strcmp(method,'TFRSCALO') | ...
   strcmp(method,'TYPE2'   ),
 Nf2=tfrrow;
elseif strcmp(method,'TFRPMH'  )| strcmp(method,'TFRRPMH' )| ...
       strcmp(method,'TFRSP'   )| strcmp(method,'TFRRSP'  )| ...
       strcmp(method,'TFRPPAGE')| strcmp(method,'TFRRPPAG')| ...
       strcmp(method,'TFRMHS'  )| strcmp(method,'TFRRGAB' )| ...
       strcmp(method,'TFRMH'   )| strcmp(method,'TFRMMCE' )| ...
       strcmp(method,'TFRRMSC' )| strcmp(method,'TFRPAGE' )| ...
       strcmp(method,'TFRGABOR')| strcmp(method,'TFRRI'   )| ...
       strcmp(method,'TFRMSC'  )| strcmp(method,'TYPE1'   )| ...
       strcmp(method,'TFRSTFT' ),
 Nf2=tfrrow/2;
else
 error('Unknown representation. Use type1 or type2');
end;

% Computation of freq (vector of frequency samples)
if strcmp(method,'TFRASPW' ) | strcmp(method,'TFRSCALO') | ...
   strcmp(method,'TFRDFLA' ) | strcmp(method,'TFRSPAW' ) | ...
   strcmp(method,'TFRUNTER') | strcmp(method,'TFRBERT' ),
 freq=eval(['p',num2str(nargin-4)]);   	% last input argument is freqs.
else
 freq=(0.5*(0:Nf2-1)/Nf2);
end

% Initialization of the current figure
zoom off; clf;
set(gcf,'Resize','On','NextPlot','Add');

% Initialization of the variables
if exist('options.mat'),
 load options
else
 threshold=5.0;
 linlogtfr=0;
 linlogspec=1;
 levelnumb=6;
 colmap=1;
 display=1;
 layout=1;
 iscolorbar=0;
 isgridsig=0;
 isgridspe=0;
 isgridtfr=0;
 issig=0;
 isspec=0;
end
choice=4;
Grid=1;
displayold=1;
fs=1.0;

if ((layout>=2 & layout<=4) & isempty(sig)),
 layout=1;
end

% Initialization of the axes
axcb   = axes('Units','normal','Visible','off','Box','On');
axsig  = axes('Units','normal','Visible','off','Box','On');
axspec = axes('Units','normal','Visible','off','Box','On');
axtfr  = axes('Units','normal','Visible','off','Box','On');
if comp(1:2)=='PC',
 set(axsig ,'fontsize',10);
 set(axspec,'fontsize',10);
 set(axtfr ,'fontsize',10);
end
set(gcf,'UserData',[axcb axsig axspec axtfr]);

% Test of analycity
if ~isempty(sig),
 for k=1:Ncol,
  spec(:,k)=abs(fft(sig(min(t):max(t),k))).^2; Nsp=length(spec);
  spec1=sum(spec(1:round(Nsp/2),k));
  spec2=sum(spec(round(Nsp/2)+1:Nsp,k));
  if spec2>spec1/10,
   disp('Be careful : the signal is not analytic!');
  end
 end
end

% Test of reality
if (Ncol==2 & ~isreal(tfr)),
 disp('Cross distribution. As the result is complex, we display the real part.');
 tfr=real(tfr);
end

while choice~=12,

 axes(axtfr);				% The current axis is the tfr axis

 if (display~=7 & choice~=2 & choice~=3 & choice<8),	% Call to tfrview
  linlog=linlogtfr+2*linlogspec;
  if choice==4,
   access=3;
  else
   access=1;
  end
  state=issig+2*iscolorbar;
  isgrid=isgridsig+2*isgridspe+4*isgridtfr;
  param = [display,linlog,threshold,levelnumb,Nf2,layout,access,state,fs,isgrid];
  map=colormap;
  if (nargin<=4),
   tfrview(tfr,sig,t,method,param,map);
  elseif (nargin==5),
   tfrview(tfr,sig,t,method,param,map,p1);
  elseif (nargin==6),
   tfrview(tfr,sig,t,method,param,map,p1,p2);
  elseif (nargin==7),
   tfrview(tfr,sig,t,method,param,map,p1,p2,p3);
  elseif (nargin==8),
   tfrview(tfr,sig,t,method,param,map,p1,p2,p3,p4);
  elseif (nargin==9),
   tfrview(tfr,sig,t,method,param,map,p1,p2,p3,p4,p5);
  end;
 elseif (display==7),			% keep in mind the last value of display
  display=displayold;
 end

 if (linlogtfr==0),			% Lin/log scale of the tfr
  linlogstr='Change to a logarithmic scale';
 else
  linlogstr='Change to a linear scale';
 end;

 % Main menu
 choice=menu ('TFRQVIEW MENU :',...
	      'Change the display mode',...
	      'Change the display layout',...
	      'Change the color map',...
	      'Change the sampling frequency',...
	      'Change the threshold',...
	      linlogstr,...
	      'Change the number of levels',...
	      'Grid',...
	      'Save',...
	      'Save options',...
	      'Print',...
	      'Close');

 if ((choice==9 | choice==11) & comp(1:2)~='PC'),
  if exist('workdir.mat'),
   load workdir;
  else
   Path=pwd;
   str1=[' The current directory is ',Path];
   disp(str1); disp('');
   continue=0;
   while continue==0,
    continue=1;
    str2=' Name of the directory (with full path) on which the files will be saved : ';
    disp(str2);
    PathWorkDir=input('       > ','s');
    if (exist(PathWorkDir)~=2),
     rep=input(' This directory doesn''t exist. Do you want to create it (y/n) ? ','s');
     if upper(rep)=='Y',
      eval(['!mkdir ',PathWorkDir]);
     else
      continue=0;
     end
    end
   end
   if PathWorkDir(length(PathWorkDir))~='/',
    PathWorkDir=[PathWorkDir,'/'];
   end
   save workdir PathWorkDir;
  end
 end

 if (choice==1), 			% Change the display mode

  display=menu('DISPLAY MODE :',...
	       'Contour',...
	       'Imagesc',...
	       'Pcolor',...
	       'Surf',...
	       'Mesh',...
	       'Cancel');

 elseif (choice==2), 			% Change the display layout

  layoutold=layout; indlo=1; layout=1;
  while layout~=6,			% While not close

   if ~iscolorbar,
    colorbarstr='Add a color bar';
   else
    colorbarstr='Remove the color bar';
   end
   if ~linlogspec,
    specstr='TFR + Spectrum (linear)';
    tfrsstr='TFR + Signal + Spectrum (linear)';
   else
    specstr='TFR + Spectrum (logarithmic)';
    tfrsstr='TFR + Signal + Spectrum (logarithmic)';
   end

   layout=menu('DISPLAY LAYOUT :',...
	       'Time-Frequency Representation',...
	       'TFR + Signal ',...
	       specstr,...
	       tfrsstr,...
	       colorbarstr,...
	       'Close');

   if (layout==3|layout==4),
    linlogspec=1-linlogspec;
   end
   if ((layout>=2 & layout<=4) & isempty(sig)),
    disp('Impossible action : the signal is unavailable');
   else
    access=2;
    linlog=linlogtfr+2*linlogspec;
    state=issig+2*iscolorbar;
    isgrid=isgridsig+2*isgridspe+4*isgridtfr;
    param = [display,linlog,threshold,levelnumb,Nf2,layout,access,state,fs,isgrid];
    map=colormap;
    if (nargin<=4),
     tfrview(tfr,sig,t,method,param,map);
    elseif (nargin==5),
     tfrview(tfr,sig,t,method,param,map,p1);
    elseif (nargin==6),
     tfrview(tfr,sig,t,method,param,map,p1,p2);
    elseif (nargin==7),
     tfrview(tfr,sig,t,method,param,map,p1,p2,p3);
    elseif (nargin==8),
     tfrview(tfr,sig,t,method,param,map,p1,p2,p3,p4);
    elseif (nargin==9),
     tfrview(tfr,sig,t,method,param,map,p1,p2,p3,p4,p5);
    end;

    axesh=get(gcf,'UserData');
    axcb=axesh(1); axsig=axesh(2); axspec=axesh(3); axtfr=axesh(4);

    if layout==1,			% Time-Frequency Representation
     axes(axsig);  cla reset; set(gca,'Visible','off');
     axes(axspec); cla reset; set(gca,'Visible','off');
     issig=0; isspec=0;
    elseif layout==2,			% TFR + Signal
     axes(axspec); cla reset; set(gca,'Visible','off');
     issig=1; isspec=0;
    elseif layout==3,			% TFR + spectrum
     axes(axsig);  cla reset; set(gca,'Visible','off');
     issig=0; isspec=1;
    elseif layout==4,			% TFR + signal + spectrum
      issig=1; isspec=1;
    elseif layout==5,			% Colorbar
     iscolorbar=1-iscolorbar;
    elseif ((layout==2|layout==3|layout==4) & isempty(sig)),
     disp('Unallowed action : the signal is anavailable.');
    end;

    lo(indlo)=layout;
    indlo=indlo+1;
   end
  end
  while (layout==6|layout==5),
   indlo=indlo-1;
   if indlo>0,
    layout=lo(indlo);
   else
    layout=layoutold;
   end
  end

 elseif (choice==3),			% Change the color map

  colmap=1;
  while colmap~=12,
   colmap=menu('COLOR MAP :',...
	       'hsv','jet','cool','bone','gray','hot','prism',...
	       'brighten','darken','permute','spin','Close');
   if colmap==1,
    colormap(hsv);
   elseif colmap==2,
    colormap(jet);
   elseif colmap==3,
    colormap(cool);
   elseif colmap==4,
    colormap(bone);
   elseif colmap==5,
    colormap(gray);
   elseif colmap==6,
    colormap(hot);
   elseif colmap==7,
    colormap(prism);
   elseif colmap==8,
    brighten(0.25);
   elseif colmap==9,
    brighten(-0.25);
   elseif colmap==10,
    c = colormap; colormap(flipud(c));
   elseif colmap==11,
    spinmap;
   end
  end

 elseif (choice==4), 			% Change the sampling frequency

  fprintf(' Old sampling frequency: %f',fs); fsold=fs;
  fs=input(' New sampling frequency: ');
  if isempty(fs), fs=fsold; end

 elseif (choice==5),			% Change the threshold

  fprintf(' Old threshold: %f', threshold); throld=threshold;
  threshold=input(' New threshold: ');
  if isempty(threshold), threshold=throld; end

 elseif (choice==6),			% Change the lin/log scale of tfr

  linlogtfr=1-linlogtfr;

 elseif (choice==7),			% Change the number of levels

  fprintf(' Old number of levels: %f',levelnumb); levelold=levelnumb;
  levelnumb=input(' New number of levels: ');
  if isempty(levelnumb), levelnumb=levelold; end

 elseif (choice==8),			% Grids

  if ~issig & ~isspec,			% No signal and no spectrum
   isgridtfr=1-isgridtfr; axes(axtfr); grid

  elseif issig & ~isspec,		% A signal, no spectrum
   Grid=1;
   while Grid~=3,
    if ~isgridsig,
     gridsigstr='Add a grid on the signal';
    else
     gridsigstr='Remove the grid on the signal';
    end
    if ~isgridtfr,
     gridtfrstr='Add a grid on the TFR';
    else
     gridtfrstr='Remove the grid on the TFR';
    end
    Grid=menu('GRID MENU :',gridsigstr,gridtfrstr,'Close');
    if Grid==1,
     isgridsig=1-isgridsig; axes(axsig); grid
    elseif Grid==2,
     isgridtfr=1-isgridtfr; axes(axtfr); grid
    end
   end

  elseif ~issig & isspec,		% No signal, a spectrum
   Grid=1;
   while Grid~=3,
    if ~isgridspe,
     gridspestr='Add a grid on the spectrum';
    else
     gridspestr='Remove the grid on the spectrum';
    end
    if ~isgridtfr,
     gridtfrstr='Add a grid on the TFR';
    else
     gridtfrstr='Remove the grid on the TFR';
    end
    Grid=menu('GRID MENU :',gridspestr,gridtfrstr,'Close');
    if Grid==1,
     isgridspe=1-isgridspe; axes(axspec); grid
    elseif Grid==2,
     isgridtfr=1-isgridtfr; axes(axtfr); grid
    end
   end

  else					% A signal and a spectrum
   Grid=1;
   while Grid~=4,
    if ~isgridsig,
     gridsigstr='Add a grid on the signal';
    else
     gridsigstr='Remove the grid on the signal';
    end
    if ~isgridspe,
     gridspestr='Add a grid on the spectrum';
    else
     gridspestr='Remove the grid on the spectrum';
    end
    if ~isgridtfr,
     gridtfrstr='Add a grid on the TFR';
    else
     gridtfrstr='Remove the grid on the TFR';
    end
    Grid=menu('GRID MENU :',gridsigstr,gridspestr,gridtfrstr,'Close');
    if Grid==1,
     isgridsig=1-isgridsig; axes(axsig); grid
    elseif Grid==2,
     isgridspe=1-isgridspe; axes(axspec); grid
    elseif Grid==3,
     isgridtfr=1-isgridtfr; axes(axtfr); grid
    end
   end
  end

 elseif (choice==9),			% Save the parameters

  f=freq*fs;
  Nmethod=length(method);
  if comp(1:2)=='PC',
   namedflt=[method(4:Nmethod),num2str(Nsig),'.mat'];
   [name,PathWorkDir] = uiputfile(namedflt, 'Save As');
  else
   namedflt=[method(4:Nmethod),num2str(Nsig)];
   nameStr=[' Name of the MAT file [',namedflt,'] : '];
   name=input(nameStr,'s');
   while (length(name)>8),
    disp(' The name must have less than 8 characters');
    name=input(nameStr,'s');
   end
   if isempty(name),
    name=namedflt;
   end
  end
  access=0;
  linlog=linlogtfr+2*linlogspec;
  state=issig+2*iscolorbar;
  isgrid=isgridsig+2*isgridspe+4*isgridtfr;
  param = [display,linlog,threshold,levelnumb,Nf2,layout,access,state,fs,isgrid];
  map=colormap;
  if (nargin<=4),
   TfrQView=['tfrqview(tfr,sig,t,method)'];
   TfrView =['clf;tfrview(tfr,sig,t,method,param,map)'];
   eval(['save ',PathWorkDir,name,...
' tfr sig t f fs method param map TfrView TfrQView']);
  elseif (nargin==5),
   TfrQView=['tfrqview(tfr,sig,t,method,p1)'];
   TfrView =['clf; tfrview(tfr,sig,t,method,param,map,p1)'];
   eval(['save ',PathWorkDir,name,...
' tfr sig t f fs method param map p1 TfrView TfrQView']);
  elseif (nargin==6),
   TfrQView=['tfrqview(tfr,sig,t,method,p1,p2)'];
   TfrView =['clf; tfrview(tfr,sig,t,method,param,map,p1,p2)'];
   eval(['save ',PathWorkDir,name,...
' tfr sig t f fs method param map p1 p2 TfrView TfrQView']);
  elseif (nargin==7),
   TfrQView=['tfrqview(tfr,sig,t,method,p1,p2,p3)'];
   TfrView =['clf; tfrview(tfr,sig,t,method,param,map,p1,p2,p3)'];
   eval(['save ',PathWorkDir,name,...
' tfr sig t f fs method param map p1 p2 p3 TfrView TfrQView']);
  elseif (nargin==8),
   TfrQView=['tfrqview(tfr,sig,t,method,p1,p2,p3,p4)'];
   TfrView =['clf; tfrview(tfr,sig,t,method,param,map,p1,p2,p3,p4)'];
   eval(['save ',PathWorkDir,name,...
' tfr sig t f fs method param map p1 p2 p3 p4 TfrView TfrQView']);
  elseif (nargin==9),
   TfrQView=['tfrqview(tfr,sig,t,method,p1,p2,p3,p4,p5)'];
   TfrView =['clf; tfrview(tfr,sig,t,method,param,map,p1,p2,p3,p4,p5)'];
   eval(['save ',PathWorkDir,name,...
' tfr sig t f fs method param map p1 p2 p3 p4 p5 TfrView TfrQView']);
  end;
  disp(' ');
  SaveTxt1=[' The file is saved in the directory ',PathWorkDir];
  SaveTxt2=['  under the name ',name,'.mat'];
  disp(SaveTxt1); disp(SaveTxt2);

  disp(' If you want to find again the exact layout of this screen, do');
  Txt1=['	load ',PathWorkDir,name,'; eval(TfrView);'];
  disp(Txt1);
  disp(' If you want to restart the display session under tfrqview, do');
  Txt2=['	load ',PathWorkDir,name,'; eval(TfrQView);'];
  disp(Txt2);

 elseif (choice==10),			% Save options

  save options fs threshold linlogtfr linlogspec levelnumb colmap ...
       display layout iscolorbar isgridsig ...
       isgridspe isgridtfr issig isspec;
  disp(' '); disp(' Options saved');

 elseif (choice==11),			% Print the current figure

  Nmethod=length(method);
  if comp(1:2)=='PC',
   namedflt=[method(4:Nmethod),num2str(Nsig),'.eps'];
   [name,PathWorkDir] = uiputfile(namedflt, 'Save As');
  else
   namedflt=[method(4:Nmethod),num2str(Nsig)];
   nameStr=[' Name of the EPS file [',namedflt,'] : '];
   name=input(nameStr,'s');
   while (length(name)>8),
    disp('The name must have less than 8 characters');
    name=input(nameStr,'s');
   end
   if isempty(name),
    name=namedflt;
   end
  end
  eval(['print -deps ',PathWorkDir,name]); disp(' ');
  PrintTxt1=[' The file is saved in the directory ',PathWorkDir];
  if comp(1:2)=='PC',
   PrintTxt2=['  under the name ',name];
  else
   PrintTxt2=['  under the name ',name,'.eps'];
  end
  disp(PrintTxt1); disp(PrintTxt2);

 end;

 if display~=7,
  displayold=display;
 end

end

if (comp(1:2)~='PC' & exist('workdir.mat')),
 !rm workdir.mat
end

disp(' ');

author: HeLiangHIT


text:
function r8mat_nonzeros_test ( )

%*****************************************************************************80
%
%% R8MAT_NONZEROS_TEST tests R8MAT_NONZEROS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 December 2014
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 4;
  a = zeros ( m, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_NONZEROS_TEST\n' );
  fprintf ( 1, '  R8MAT_NONZEROS counts nonzeros in an R8MAT.\n' );

  c1 = 0;
  for i = 1 : m
    for j = 1 : n
      if ( mod ( i, 2 ) == 0 && mod (  j, 2 ) == 0 )
        a(i,j) = 1.0;
        c1 = c1 + 1;
      else
        a(i,j) = 0.0;
      end
    end
  end

  r8mat_print ( m, n, a, '  Matrix A:' );

  c2 = r8mat_nonzeros ( m, n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Expected nonzeros = %d\n', c1 );
  fprintf ( 1, '  Computed nonzeros = %d\n', c2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function self = seisan(varargin)
    % READ_CATALOG.SEISAN Read a Seisan event database - consisting of
    % S-files. Gathers a list of S-files matching the request and then uses
    % the read_sfile function, to return a structure for each S-file.

    debug.printfunctionstack('>')

    % Process input arguments
    % Process input arguments
    p = inputParser;
    p.addParameter('dbpath', '', @isstr);
    p.addParameter('startTime', -Inf);
    p.addParameter('endTime', Inf);
    p.addParameter('minimumMagnitude', -Inf, @isnumeric);
    p.addParameter('maximumMagnitude', Inf, @isnumeric);
    p.addParameter('minimumLatitude', -90., @isnumeric);
    p.addParameter('maximumLatitude', 90., @isnumeric);
    p.addParameter('minimumLongitude', -180., @isnumeric);
    p.addParameter('maximumLongitude', 180., @isnumeric);
    p.addParameter('minimumDepth', -12, @isnumeric);
    p.addParameter('maximumDepth', 6400., @isnumeric);
    p.addParameter('minimumRadius', 0., @isnumeric);
    p.addParameter('maximumRadius', deg2km(180.0), @isnumeric);
    p.addParameter('boxcoordinates', [], @isnumeric);    %[minLat, maxLat, minLon, maxLon]   % use NaN as a wildcard
    p.addParameter('radialcoordinates', [], @isnumeric); % [Lat, Lon, MaxRadius, MinRadius]   % MinRadius is optional
    p.addParameter('addarrivals', false, @islogical);

    % CUSTOM PARAMETERS
    p.addParameter('subclass', '*', @ischar);

    p.parse(varargin{:});
    fields = fieldnames(p.Results);
    for i=1:length(fields)
        field=fields{i};
        % val = eval(sprintf('p.Results.%s;',field));
        val = p.Results.(field);
        eval(sprintf('%s = val;',field));
    end

    self.request.dbpath = dbpath;
    self.request.startTime = startTime;
    self.request.endTime = endTime;

%     % SCAFFOLD: might have deleted something here. These 2 lines might not
%     % be correct
%     self.request.minimumMagnitude = minimumMagnitude;
%     p.addParameter('minimumMagnitude', [], @isnumeric);
%
%     p.addParameter('maximumMagnitude', [], @isnumeric);
%     p.addParameter('minimumLatitude', [], @isnumeric);
%     p.addParameter('maximumLatitude', [], @isnumeric);
%     p.addParameter('minimumLongitude', [], @isnumeric);
%     p.addParameter('maximumLongitude', [], @isnumeric);
%     p.addParameter('minimumDepth', [], @isnumeric);
%     p.addParameter('maximumDepth', [], @isnumeric);
%     p.addParameter('minimumRadius', [], @isnumeric);
%     p.addParameter('maximumRadius', [], @isnumeric);
%     p.addParameter('boxcoordinates', @isnumeric);    %[minLat, maxLat, minLon, maxLon]   % use NaN as a wildcard
%     p.addParameter('radialcoordinates', @isnumeric); % [Lat, Lon, MaxRadius, MinRadius]   % MinRadius is optional
%     p.addParameter('addarrivals', false, @islogical);

    if numel(boxcoordinates)==4
        minimumLatitude = boxcoordinates(1);
        maximumLatitude = boxcoordinates(2);
        minimumLongitude = boxcoordinates(3);
        maximumLongitude = boxcoordinates(4);
    end

    if numel(radialcoordinates)==4
        centerLatitude = radialcoordinates(1);
        centerLongitude = radialcoordinates(2);
        maximumRadius = radialcoordinates(3);
        %minimumRadius = radialcoordinates(4);
    end

    % Check start & end times
    snum = Catalog.read_catalog.ensure_dateformat(startTime);
    enum = Catalog.read_catalog.ensure_dateformat(endTime);

    if ~exist(dbpath, 'dir')
        fprintf('Directory %s not found\n',dbpath);
        self = struct;
        return;
    end

    % get dir list of matching sfilesSfile.list_sfiles
    sfiles = Sfile.list_sfiles(dbpath, snum, enum);

    % loop over sfiles
    for i=1:length(sfiles)
        % read
        fprintf('Processing %s\n',fullfile(sfiles(i).dir, sfiles(i).name));
        thiss = Sfile(fileread(fullfile(sfiles(i).dir, sfiles(i).name)));
        try
            s(i)=thiss;
        catch
            s(i)
            thiss
            error('Wrong number of fields?')
        end

        % add to catalog
        dnum(i)  = s(i).otime;
        etype{i} = s(i).subclass;
        lat(i) = s(i).latitude;
        lon(i) = s(i).longitude;
        depth(i) = s(i).depth;

        % SCAFFOLD
        mag(i) = NaN;
        try
            sfile_mags = [s(i).magnitude.value];
            if ~isempty(sfile_mags)
                disp('**************** ********************')
                mag(i) = max(sfile_mags);
            end
        end

        % SCAFFOLD also use durations (bbdur) and ampengfft info
        % Compute a magnitude from amp & eng, but need to know where
        % stations are. I can save these as MA and ME, to distinguish from
        % Ml, Ms, Mb, Mw if those exist
    end
    magtype = {};

    % save request
    fields = fieldnames(p.Results);
    for i=1:length(fields)
        field=fields{i};
        eval(sprintf('request.%s = eval(field);',field));
    end

    % create Catalog object
    self = Catalog(dnum', lon', lat', depth', mag', magtype', etype', 'request', request);
    request.startTime = snum;
    request.endTime = enum;
    self.request = request;
    debug.printfunctionstack('<')
end

author: geoscience-community-codes

repo: GISMO


text:
function rtk=detslip_tdcp(rtk,cur_obs,old_obs,nav,f)

global glc
dt=timediff(cur_obs(1).time,old_obs(1).time);
if abs(dt)>5
    return;
end

cur_sv=satposs(cur_obs,nav,rtk.opt.sateph);
old_sv=satposs(old_obs,nav,rtk.opt.sateph);

cur_pos=rtk.ins.pos;     cur_rr=blh2xyz(cur_pos);
old_pos=rtk.ins.oldpos;  old_rr=blh2xyz(old_pos);

n_cur=size(cur_obs,1);
n_old=size(old_obs,1);

nv=0; vsat=zeros(n_cur,1); v=zeros(n_cur,1);

for i=1:n_cur

    cur_obsi=cur_obs(i); cur_svi=cur_sv(i);
    lam=nav.lam(cur_obsi.sat,:);
    [sys,~]=satsys(cur_obsi.sat);
    if rtk.mask(sys)==0,continue;end

    if lam(f)==0,continue;end

    old_idx=0;
    for j=1:n_old
        if cur_obsi.sat==old_obs(j).sat
            old_idx=j;break;
        end
    end
    if old_idx==0,continue;end

    old_obsi=old_obs(old_idx); old_svi=old_sv(old_idx);

    if cur_svi.svh~=0||old_svi.svh~=0||norm(cur_svi.pos)<=0||norm(old_svi.pos)<=0
        continue;
    end
    if cur_obsi.L(f)==0||old_obsi.L(f)==0,continue;end

    [cur_range,~]=geodist(cur_svi.pos,cur_rr);
    [old_range,~]=geodist(old_svi.pos,old_rr);

    tdcp_obs=cur_obsi.L(f)*lam(f)-old_obsi.L(f)*lam(f);
    tdcp_pre=cur_range-old_range;

    v(i,1)=tdcp_obs-tdcp_pre+glc.CLIGHT*(cur_svi.dts-old_svi.dts);

    nv=nv+1;

    vsat(i)=1;

end

dt_drift=rtk.x(rtk.ic+6);
v=v-dt_drift;

slip=zeros(n_cur,1);
for i=1:n_cur
    sat=cur_obs(i).sat;
    if vsat(i)==0
        rtk.sat(sat).slip(f)=bitor(rtk.sat(sat).slip(f),1);
        continue;
    end

    vi=abs(v(i));
    v_other=abs(v);v_other(i)=[];
    ave_v=sum(abs(v_other))/(nv-1);
    ave_distance=sum(abs(vi-v_other))/(nv-1);
    if ave_distance>3*ave_v&&abs(v(i))>0.15
        slip(i)=1;
    end
end

idx=find(slip==1);
ratio=size(idx,1)/n_cur;
if ratio>=0.5
    return;
end

for i=1:n_cur
    sat=cur_obs(i).sat;
    if slip(i)==1
        rtk.sat(sat).slip(f)=bitor(rtk.sat(sat).slip(f),slip(i));
    end
end

return

repo: GINav


text:
function [merged_timestamps, merged_signal] = ...
    merge_samples4(timestamps1, samples1, timestamps2, samples2, ...
                   timestamps3, samples3, timestamps4, samples4)
    % Merge samples given their timestamps

    [merged_timestamps, ix] = sort([timestamps1(:); timestamps2(:); ...
        timestamps3(:); timestamps4(:)]);
    merged_signal = [samples1(:); samples2(:); samples3(:); samples4(:)];
    labels = [ones(length(samples1), 1); ones(length(samples2), 1) * 2; ...
              ones(length(samples3), 1) * 3; ones(length(samples4), 1) * 4];
    merged_signal = merged_signal(ix);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function test_suite=test_moxunit_fieldtrip_util_parse_walltime
% #   For MOxUnit_fieldtrip's copyright information and license terms,   #
% #   see the COPYING file distributed with MOxUnit_fieldtrip.           #

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

function varargout=randint(mx)
    n=nargout;
    varargout=arrayfun(@(unused)ceil(rand()*mx),1:n,'UniformOutput',false);

function assert_output(input, expected_output)
    output=moxunit_fieldtrip_util_parse_walltime(input);
    assertEqual(output,expected_output);

function test_fieldtrip_util_parse_walltime
    % seconds
    many_ss=randint(1000);
    assert_output(sprintf('%d',many_ss),...
                                many_ss);

    % hours, minutes, seconds
    [hh,mm,ss]=randint(60);
    assert_output(sprintf('%02d:%02d:%02d',hh,mm,ss),...
                                (hh*60+mm)*60+ss);

    % not a number
    [hh,mm,ss]=randint(60);
    assert_output('foo',[]);

    % should deal with comments after
    assert_output('  3 % comment',3);

function test_fieldtrip_util_parse_walltime_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
            moxunit_fieldtrip_util_parse_walltime(varargin{:}),'');
    % inputs that are not string throw an exception
    aet(struct)
    aet([1 2])
    aet(1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [B3, B2, B1, B] = construct_B_matrices(measurements)
% function [B3, B2, B1, B] = construct_B_matrices( measurements )
%
% Given a measurement struct of the same format that the SE-Sync method
% expects, this function computes and returns the (sparse) system matrices
% B1, B2, B3, and B defined in the first section of the appendix
% "Reformulating the estimation problem" in the SE-Sync paper

d = length(measurements.t{1});  % Dimension of observations
n = max(max(measurements.edges));  % Number of poses
m = size(measurements.edges, 1);  % Number of measurements

% B3 matrix:
B3_nnz = (d^3 + d^2)*m;
B3_rows = zeros(1, B3_nnz);
B3_cols = zeros(1, B3_nnz);
B3_vals = zeros(1, B3_nnz);

for e = 1:m
    tail = measurements.edges(e,1);
    head = measurements.edges(e,2);

    sqkappa = sqrt(measurements.kappa{e});
    Rt = measurements.R{e}';

    % Block entries corresponding to the tail of this edge

    for r = 1:d
        for c = 1:d

            % Block representation of the -sqrt(kappa) *Rt(i,j) * I_d that
            % appears in the Kronecker product

            idxs = [(d^3 + d^2)*(e-1) + d^2*(r-1) + d*(c - 1) + 1 : (d^3 + d^2)*(e-1) + d^2*(r-1) + d*(c - 1) + d];

            B3_rows(idxs) = [d^2 * (e-1) + d*(r-1) + 1 : d^2 * (e-1) + d*(r-1) + d];
            B3_cols(idxs) = [d^2 * (tail - 1) + d*(c-1) + 1 : d^2 * (tail - 1) + d*(c-1) + d];
            B3_vals(idxs) = -sqkappa * Rt(r,c) * ones(1, d);
        end
    end

    % Large d x d diagonal block entry corresponding to the head of this
    % edge

    idxs = [d^3 * e + d^2*(e-1) + 1 : (d^3 + d^2)*e];

    B3_rows(idxs) = [d^2 *(e-1) + 1 : d^2 * e];
    B3_cols(idxs) = [d^2 *(head -1) + 1 : d^2 * head];
    B3_vals(idxs) = sqkappa * ones(1, d^2);
end

B3 = sparse(B3_rows, B3_cols, B3_vals, d^2*m, d^2*n);

if nargout > 1

    % B2 matrix:
    B2_nnz = m*d^2;
    B2_rows = zeros(1, B2_nnz);
    B2_cols = zeros(1, B2_nnz);
    B2_vals = zeros(1, B2_nnz);

    for e = 1:m
        tail = measurements.edges(e, 1);

        sqtau = sqrt(measurements.tau{e});
        tij = measurements.t{e};

        for idx = 1:d
            % Block representation of -sqrt(tau)*t_ij(i) * I_d in the Kronecker
            % product:

            B2_rows(d^2 * (e-1) + d*(idx-1) + 1 : d^2 * (e-1) + d*idx) = [d*(e-1) + 1 : d*e];
            B2_cols(d^2 * (e-1) + d*(idx-1) + 1 : d^2 * (e-1) + d*idx) = [d^2 * (tail-1) + d*(idx - 1) + 1 : d^2 * (tail-1) + d*idx];
            B2_vals(d^2 * (e-1) + d*(idx-1) + 1 : d^2 * (e-1) + d*idx) = -sqtau * tij(idx)*ones(1,d);
        end
    end

    B2 = sparse(B2_rows, B2_cols, B2_vals, d*m, d^2*n);

    % B1 matrix:

    B1_nnz = 2*d*m;
    B1_rows = zeros(1, B1_nnz);
    B1_cols = zeros(1, B1_nnz);
    B1_vals = zeros(1, B1_nnz);

    for e = 1:m
        tail = measurements.edges(e,1);
        head = measurements.edges(e,2);

        % Block entry corresponding to the tail of this edge
        B1_rows(2*d*(e-1) + 1 : 2*d*e - d) = [d*(e-1) + 1 : d*e];
        B1_cols(2*d*(e-1) + 1 : 2*d*e - d) = [d*(tail - 1) + 1 : d*tail];
        B1_vals(2*d*(e-1) + 1 : 2*d*e - d) = -sqrt(measurements.tau{e})*ones(1, d);

        % Block entry corresponding to the head of this edge
        B1_rows(2*d*e - d + 1 : 2*d*e) = [d*(e-1) + 1 : d*e];
        B1_cols(2*d*e - d + 1 : 2*d*e) = [d*(head - 1) + 1 : d*head];
        B1_vals(2*d*e - d + 1 : 2*d*e) = sqrt(measurements.tau{e})*ones(1, d);
    end

    B1 = sparse(B1_rows, B1_cols, B1_vals, d*m, d*n);

end

if nargout >= 4
    B = [B1, B2;
        sparse(size(B3, 1), size(B1, 2)), B3];
end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [img_r1,img_r5,img_r10,img_med, img_map, txt_r1,txt_r5,txt_r10,txt_med,txt_map] = evaluate( ff1, ff2 ,img_id, txt_id)

parfor i = 1:size(ff1,1)
    %disp(i);
    tmp = ff1(i,:);
    score = tmp*(ff2)';
    [s, index] = sort(score, 'descend');
    good_index = find(txt_id==img_id(i));
    junk_index = [];
    [ap(i), CMC(i, :)] = compute_AP_rerank(good_index, junk_index, index);
end
CMC = mean(CMC);
rank = find(CMC>0.5);
img_r1 = CMC(1); img_map = mean(ap); img_med=rank(1);
img_r5 = CMC(5); img_r10 = CMC(10);

ap = [];
CMC = [];
% txt query
parfor i = 1:size(ff2,1)
    %disp(i);
    tmp = ff2(i,:);
    score = tmp*(ff1)';
    [s, index] = sort(score, 'descend');
    good_index = find(img_id==txt_id(i));
    %query_title = test_caption(i);
    junk_index = [];
    [ap(i), CMC(i, :)] = compute_AP_rerank(good_index, junk_index, index);
end
CMC = mean(CMC);
rank = find(CMC>0.5);
txt_r1 = CMC(1); txt_map = mean(ap); txt_med=rank(1);
txt_r5 = CMC(5); txt_r10 = CMC(10);

end

author: layumi

repo: Image-Text-Embedding


text:
%infsup   computes inf-sup eigenvalue distribution
%   IFISS scriptfile: DJS; 7 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
fprintf('inf-sup eigenvalue estimation \n')
leig=default('eig/eigs? 0/1 (default is eig)',0);
stab=default('stable/stabilised? 0/1 (default is stable)',0);
[np,nu]=size(Bst);
S=Bst*(Ast\Bst');
if stab==1, S=S+beta*C; end
%
% compute discrete eigenvalues
S=0.5*(S+S');Q=0.5*(Q+Q');
if leig==1
   [V,D,flag] = eigs(S,Q,4,'sm');
   e=diag(D); e=sort(real(e));
   fprintf('inf-sup constant is %g \n',e(2))
else
   [V,D]=eig(full(S),full(Q));
e=diag(D); [e,k]=sort(real(e)); V=V(:,k); pp=length(e);
fprintf('inf-sup constant is %g \n',e(2))
fprintf('     upper bound is %g \n',e(pp))
end


text:
function process_scores_Tex(conf, scores, nr)

fprintf('Writing results to Tex summary...\n');

valid_id = [2 3 4 5 6 7 8 9 11];
fprintf('\n');

    for j = 1:length(conf.desc)
        if sum(valid_id==j)==0
            continue;
        end;
        fprintf('& \\multicolumn{2}{|c|}{%s}', conf.desc{j});
    end
    fprintf('\\\\\n');

fprintf('image');
    for j = 1:length(conf.desc)
        if sum(valid_id==j)==0
            continue;
        end;
        fprintf('& PSNR & Time');
    end
    fprintf('\\\\\n');

for i =1:nr
    [p, f, x] = fileparts(conf.filenames{i});
    fprintf('%s',f);
    for j = 1:length(conf.desc)
        if sum(valid_id==j)==0
            continue;
        end;

        fprintf(' & %.1f & %.1f', scores(i,j), conf.countedtime(j-1,i));
    end
    fprintf('\\\\\n');
end

fprintf('average');
for j = 1:11
    if sum(valid_id==j)==0
       continue;
    end;

    fprintf(' & %.2f & %.2f', mean(scores(:,j)), mean(conf.countedtime(j-1,:)));
end
fprintf('\\\\\n\n\n');

author: thomas-koehler

repo: SupER


text:
function [fcfac,fcfb,fcfacreal,fcfbreal]=flopcounts(a,M,L,Lg)
%
%   Compute the flop count for the factorization algorithm.

N=L/a;

[c,h_a,h_m]=gcd(a,M);
h_a=-h_a;
p=a/c;
q=M/c;
d=N/q;

fcfac     = L*8*q+4*L*(1+q/p)*log2(d)+4*M*N*log2(M);

fcfb      = 8*L*Lg/a+4*M*N*log2(M);

fcfacreal = L*4*q+2*L*(1+q/p)*log2(d)+2*M*N*log2(M);

fcfbreal  = 2*L*Lg/a+2*M*N*log2(M);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% The COBRAToolbox: testSurfNet.m
%
% Purpose:
%     - testSurfNet tests the functionality of surfNet
%       and compares it to expected data.
%
% Authors:
%     - Lemmer El Assal March 2017
%     - Siu Hung Joshua Chan July 2017
%

global CBTDIR

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testSurfNet'));
cd(fileDir);

model = getDistributedModel('ecoli_core_model.mat');

% generate flux data for testing
s = optimizeCbModel(model, 'max', 'one');
model2 = changeRxnBounds(model, {'EX_glc(e)'; 'EX_fru(e)'}, [0; -10], {'L'; 'L'});
s2 = optimizeCbModel(model2, 'max', 'one');
fluxMatrix = [s.x, s2.x];
rxnsDiff = model.rxns(abs(fluxMatrix(:, 1) - fluxMatrix(:, 2)) > 1e-6);
[m, n] = size(model2.S);

% remove old generated file
if exist('surfNet.txt', 'file')
    delete('surfNet.txt');
end

% check normal functionalities
diary('surfNet.txt');

% start with a metabolite
metrxn = '13dpg[c]';
surfNet(model, metrxn);
% continue with a reaction (through clicking only)
surfNet([], 'GAPD', 0, NaN, 0, 1, [], 0);
% continue with a metabolite (through clicking only)
surfNet([], 'g3p[c]', 0, NaN, 0, 1, [], 0);
% print objective reactions given no second input
surfNet(model);
% print a list of reactions without showing details of metabolites
surfNet(model, {'GAPD'; 'FBA'}, [], [], [], 0);
% print *.metNames in reaction formulae
surfNet(model, {'13dpg[c]'; 'GAPD'}, 1);
% show previous steps (through clicking only)
surfNet([], [], 1, NaN, 0, 1, [], 0, [], struct('showPrev', true));

% print with a fixed number of characters per line
surfNet(model, [model.mets(1:10)'; model.rxns(1:10)'], [], [], [], 0, [], 60);
% show previous steps (through clicking only)
surfNet([], [], 0, NaN, 0, 0, [], 60, [], struct('showPrev', true));

% print connected reactions and the corresponding fluxes in a flux vector
surfNet(model, 'pyr[c]', [], s.x, 0);
% print connected reactions with nonzero fluxes only in a flux vector
surfNet(model, 'pyr[c]', [], s.x);
% compare multiple flux vectors
surfNet(model, rxnsDiff, [], fluxMatrix, [], 0);
surfNet(model, rxnsDiff, [], fluxMatrix', [], 0);

% print customized fields
[model2.newRxnProp, model2.newMetProp] = deal(cell(1, n), cell(1, m));
model2.newRxnProp{findRxnIDs(model2, 'GAPD')} = {'a', 'b'};  % cellstr with cells
model2.newMetProp{findMetIDs(model2, '13dpg[c]')} = {'c', 'd'};
[model2.newRxnProp2, model2.newMetProp2] = deal(char('C' * ones(n, 1)), char('E' * ones(m, 1)));  % character arrays
surfNet(model2, {'13dpg[c]'; 'GAPD'}, [], [], [], [], {'metFormulas', 'subSystems', ...
    'grRules', 'b', 'c', 'ub', 'newRxnProp', 'newMetProp', 'newRxnProp2', 'newMetProp2'});
surfNet(model2, {'13dpg[c]'; 'GAPD'}, [], [], [], [], 'lb');

% starting with a gene
surfNet(model, 'b1241')
% continue with another gene shown in the reactions
surfNet([], 'b1478', 0, NaN, 0, 1, [], 0);

% test the search function
surfNet(model, 'glucose')
surfNet(model, 'glucose', 't', 0.6)  % with a more relaxed threshold
surfNet(model, 'ATP')  % a query term with multiple matches in the same object
diary off;

% load the text files
[text1, text2] = deal('');
f = fopen('refData_surfNet.txt', 'r');
l = fgets(f);
while ~isequal(l, -1)
    text1 = [text1, l];
    l = fgets(f);
end
fclose(f);

f = fopen('surfNet.txt', 'r');
l = fgets(f);
while ~isequal(l, -1)
    text2 = [text2, l];
    l = fgets(f);
end
fclose(f);
% locations of linebreaks not constant. Replace linebreaks and consecutive spaces with single space
text1 = regexprep(text1, '\s*', ' ');
text2 = regexprep(text2, '\s*', ' ');

% remove the generated file
delete('surfNet.txt');

% compare the similarity of text using a simple scheme
% (no need to solve a DP here, the two strings are supposed to be highly similar)
[j1, j2, match] = deal(1, 1, 0);
while j1 <= numel(text1) && j2 <= numel(text2)
    if ~strcmp(text1(j1), text2(j2))
        % find the next closest identical character
        [j1skip, j2skip] = deal(1);
        while j2 + j2skip <= numel(text2) && ~strcmp(text1(j1), text2(j2 + j2skip))
            j2skip = j2skip + 1;
        end
        while j1 + j1skip <= numel(text1) && ~strcmp(text1(j1 + j1skip), text2(j2))
           j1skip = j1skip + 1;
        end
        % take the closer identical character from the two strings
        [j1, j2] = deal(j1 + (j1skip <= j2skip) * j1skip, j2 + (j1skip > j2skip) * j2skip);
    end
    match = match + 1;
    [j1, j2] = deal(j1 + 1, j2 + 1);
end
% subtract the last false match if the comparison is terminated due to the
% end of the strings rather than match found
match = match - (j1 == numel(text1) + 2 | j2 == numel(text2) + 2);
score = ((match / numel(text1)) * (match / numel(text2))) ^ 0.5;

fprintf('Compare the printed with the expected results ...\n')
assert(score > 1 - 1e-3);  % some mismatches due to linebreaks and space
fprintf('\nSuccess. Finish testing normal functionalities of surfNet.\n')

% check warnings
diary('surfNet.txt');

% fields not printable
%surfNet(model, object, metNameFlag, flux, nonzeroFluxFlag, showMets, printFields, charPerLine, similarity)
surfNet(model2, '13dpg[c]', [], [], [], [], {'S'});
surfNet(model2, '13dpg[c]', [], [], [], [], {{}, {'rxnGeneMat'}});
% non-existing met/rxn or incorrect input type
surfNet(model2, 'NOTEXIST');
surfNet({{}});

diary off;

% load the text file
textSurfNet = '';
f = fopen('surfNet.txt', 'r');
l = fgets(f);
while ~isequal(l, -1)
    textSurfNet = [textSurfNet, l];
    l = fgets(f);
end
textSurfNet = regexprep(textSurfNet, '\s*', ' ');
fclose(f);

fprintf('Compare the printed warnings with the expected results ...\n')
assert(contains(textSurfNet, 'Warning: surfNet does not support showing S. Ignore.'))
assert(contains(textSurfNet, 'Warning: surfNet does not support showing rxnGeneMat. Ignore.'))
assert(contains(textSurfNet, 'Warning: No related mets, rxns or genes are found from the search. Please try other query terms.'))
assert(contains(textSurfNet, 'Warning: The query term must be either a string or an array of string.'))
assert(contains(textSurfNet, '''NOTEXIST'' is not a metabolite, reaction or gene of the model. Searching for related objects:'))

fprintf('\nSuccess. Finish testing warning output of surfNet.\n')

% remove the generated file
delete('surfNet.txt');

% print a random reaction when the 2nd input 'metrxn' is not given and
% no objective reactions exist.
fprintf('Test printing random reactions ...\n')
model2.c(:) = 0;
surfNet(model2)
surfNet(model2, [], [], fluxMatrix);

fprintf('Test error messages ...\n')
% error messages
% no initialized model exists
clear surfNet
try
    surfNet;
    assert(false)
catch ME
    assert(strcmp(ME.message, 'The persistent variable modelLocal in surfNet is empty. Please supply a COBRA model.'))
end

% incorrect flux matrix input
try
    surfNet(model2, 'fdp[c]', [], fluxMatrix(1:(end - 1), :));
    assert(false)
catch ME
    assert(strcmp(ME.message, 'Input flux vector has incorrect dimension.'))
end

% add two fields with incorrect size and unsupported datatype
[model2.rxnABC, model2.metABC] = deal(repmat(struct(), n - 1, 2), repmat(struct(), m - 1, 2));
incorrectFieldInput = {{'NOTEXIST'}; 'description'; {'rxnABC', 'metABC'}; 0};
% nonexistent fields
errMsg = {'The following field(s) is(are) not in the model: NOTEXIST'; ...
    ... % unrecognizable fields
    'The following field(s) cannot be recognized as met or rxn field(s): description'; ...
    ... % incorrect size and data type
    sprintf(['Incorrect size of the following met field(s): metABC\n', ...
    'The following met field(s) is(are) neither numeric nor cell array of characters: metABC\n', ...
    'Incorrect size of the following rxn field(s): rxnABC\n', ...
    'The following rxn field(s) is(are) neither numeric nor cell array of characters: rxnABC']); ...
    ... % incorrect input format
    ['The value of ''printFields'' is invalid. Must be (1) a cell array of two cells, ', ...
    '1st cell being a character array for met fields and 2nd for rxn fields, ', ...
    'or (2) a character array of field names recognizable from the field names or the sizes.']};
for j = 1:numel(incorrectFieldInput)
    try
        surfNet(model2, 'fdp[c]', [], [], [], [], incorrectFieldInput{j});
        assert(false)
    catch ME
        assert(strcmp(ME.message, errMsg{j}))
    end
end

fprintf('\nSuccess. Finish testing error messages of surfNet.\n')

% test name-value pair input
for j = 1:3
    if exist(['surfNet' num2str(j) '.txt'], 'file')
        delete(['surfNet' num2str(j) '.txt']);
    end
end
surfNet(model, metrxn);

for j = 1:3
    if exist(['surfNet' num2str(j) '.txt'], 'file')
        delete(['surfNet' num2str(j) '.txt'])
    end
end

% reference data
diary('surfNet1.txt')
surfNet(model, 'GAPD');
surfNet([], [], [], [], [], 0);
% print a list of reactions without showing details of metabolites
surfNet(model, {'GAPD'; 'FBA'}, [], [], [], 0);
% print *.metNames in reaction formulae
surfNet(model, {'13dpg[c]'; 'GAPD'}, 1);
% print with a fixed number of characters per line
surfNet(model, [model.mets(1:10)'; model.rxns(1:10)'], [], [], [], 0, [], 60);
surfNet(model, 'pyr[c]', [], s.x, 0);
surfNet(model, 'pyr[c]', [], s.x);
surfNet(model, rxnsDiff, [], fluxMatrix, [], 0);
surfNet(model2, '13dpg[c]', [], [], [], [], {'default', 'lb', 'subSystems'});
diary off

% name-value pair input
diary('surfNet2.txt')
surfNet(model, 'GAPD');
surfNet('showMets', 0);
surfNet(model, {'GAPD'; 'FBA'}, 'showMets', 0)
surfNet(model, {'13dpg[c]'; 'GAPD'}, 'metNameFlag', 1);
surfNet(model, [model.mets(1:10)'; model.rxns(1:10)'], 'showMets', 0, 'charPerLine', 60);
surfNet(model, 'pyr[c]', 'flux', s.x, 'nonzeroFluxFlag', 0);
surfNet(model, 'pyr[c]', 'flux', s.x);
surfNet(model, rxnsDiff, 'flux', fluxMatrix, 'showMets', 0);
surfNet(model2, '13dpg[c]', 'printFields', {'default', 'lb', 'subSystems'});
diary off

% name-value pair input with partial matching
diary('surfNet3.txt')
surfNet(model, 'GAPD');
surfNet('s', 0);
surfNet(model, {'GAPD'; 'FBA'}, 's', 0)
surfNet(model, {'13dpg[c]'; 'GAPD'}, 'm', 1);
surfNet(model, [model.mets(1:10)'; model.rxns(1:10)'], 's', 0, 'c', 60);
surfNet(model, 'pyr[c]', 'f', s.x, 'n', 0);
surfNet(model, 'pyr[c]', 'f', s.x);
surfNet(model, rxnsDiff, 'f', fluxMatrix, 's', 0);
surfNet(model2, '13dpg[c]', 'p', {'d', 'lb', 'sub'});
diary off

textSurfNet = repmat({''}, 3, 1);
for j = 1:3
    f = fopen(['surfNet' num2str(j) '.txt'], 'r');
    l = fgets(f);
    while ~isequal(l, -1)
        textSurfNet{j} = [textSurfNet{j}, l];
        l = fgets(f);
    end
    fclose(f);
    delete(['surfNet' num2str(j) '.txt']);
end
assert(isequal(textSurfNet{1}, textSurfNet{2}))
assert(isequal(textSurfNet{1}, textSurfNet{3}))

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function region = p07_region ( )

%*****************************************************************************80
%
%% P07_REGION returns the name of the integration region for problem 07.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'box';

  return
end


text:
function matrixOut = smooth2a(matrixIn,Nr,Nc)
% Smooths 2D array data.  Ignores NaN's.
%
%function matrixOut = smooth2a(matrixIn,Nr,Nc)
%
% This function smooths the data in matrixIn using a mean filter over a
% rectangle of size (2*Nr+1)-by-(2*Nc+1).  Basically, you end up replacing
% element "i" by the mean of the rectange centered on "i".  Any NaN
% elements are ignored in the averaging.  If element "i" is a NaN, then it
% will be preserved as NaN in the output.  At the edges of the matrix,
% where you cannot build a full rectangle, as much of the rectangle that
% fits on your matrix is used (similar to the default on Matlab's builtin
% function "smooth").
%
% "matrixIn": original matrix
% "Nr": number of points used to smooth rows
% "Nc": number of points to smooth columns.  If not specified, Nc = Nr.
%
% "matrixOut": smoothed version of original matrix
%
%
% 	Written by Greg Reeves, March 2009.
% 	Division of Biology
% 	Caltech
%
% 	Inspired by "smooth2", written by Kelly Hilands, October 2004
% 	Applied Research Laboratory
% 	Penn State University
%
% 	Developed from code written by Olof Liungman, 1997
% 	Dept. of Oceanography, Earth Sciences Centre
% 	Gďż˝teborg University, Sweden
% 	E-mail: olof.liungman@oce.gu.se

%
% Initial error statements and definitions
%
if nargin < 2, error('Not enough input arguments!'), end

N(1) = Nr;
if nargin < 3, N(2) = N(1); else N(2) = Nc; end

if length(N(1)) ~= 1, error('Nr must be a scalar!'), end
if length(N(2)) ~= 1, error('Nc must be a scalar!'), end

%
% Building matrices that will compute running sums.  The left-matrix, eL,
% smooths along the rows.  The right-matrix, eR, smooths along the
% columns.  You end up replacing element "i" by the mean of a (2*Nr+1)-by-
% (2*Nc+1) rectangle centered on element "i".
%
[row,col] = size(matrixIn);
eL = spdiags(ones(row,2*N(1)+1),(-N(1):N(1)),row,row);
eR = spdiags(ones(col,2*N(2)+1),(-N(2):N(2)),col,col);

%
% Setting all "NaN" elements of "matrixIn" to zero so that these will not
% affect the summation.  (If this isn't done, any sum that includes a NaN
% will also become NaN.)
%
A = isnan(matrixIn);
matrixIn(A) = 0;

%
% For each element, we have to count how many non-NaN elements went into
% the sums.  This is so we can divide by that number to get a mean.  We use
% the same matrices to do this (ie, "eL" and "eR").
%
nrmlize = eL*(~A)*eR;
nrmlize(A) = NaN;

%
% Actually taking the mean.
%
matrixOut = eL*matrixIn*eR;
matrixOut = matrixOut./nrmlize;

author: Sable

repo: mcbench-benchmarks


text:
function C = subframecost(params,ref,tgt,cons,varargin)
% function C = subframecost(params,ref,tgt,cons,{projprojmodell})
%
% Compute sub-frame cost of synchronization parameters

% read in variable arguments
projmodel	= 'affine';
	if (length(varargin) > 0), projmodel	= varargin{1}; end

% extract beta and dt from params/cons
dt = params(1);
if isempty(cons),	a = params(2);
else,							a = cons(1);
end

% initialize some variables
npoints = size(ref,2) / 2;
nref		= size(ref,1);
ntgt		= size(tgt,1);

% compute corresponding target frames and mask out those out of range
fref	= 1:nref;
ftgt	= (a*fref) + dt;
mask	= (ftgt >= 1) & (ftgt <= size(tgt,1));

% if all frames are out of range then return infinite cost
if ~any(mask)
	C = inf; return
end

fref	= fref(mask);
ftgt	= ftgt(mask);

% compute sub-frame match cost over reference sequence
Cvec		= zeros(length(fref),1);

onevec	= ones(1,npoints);
zerovec	= zeros(1,3);
A				= [];

for f = 1:length(fref)
	refvec 	= reshape(ref(fref(f),:)',2,npoints);

	% interpolate frames of target sequence
	beta	 	= rem(ftgt(f),1);
	tgtvec 	= (1-beta)*reshape(tgt(floor(ftgt(f)),:),2,npoints) + ...
							(beta)*reshape(tgt(ceil(ftgt(f)),:),2,npoints);

	switch (projmodel)
		case{'affine'}
			A 	= [	refvec;
							tgtvec ];
			s 	= svd(A,0);	Cvec(f)	= s(4)*s(4);

		case{'fund','perspective'}
			x1	= [refvec; onevec];
			x2	= [tgtvec; onevec];
			A		= [	x1(1,:).*x2(1,:); x1(1,:).*x2(2,:); x1(1,:).*x2(3,:);
							x1(2,:).*x2(1,:); x1(2,:).*x2(2,:); x1(2,:).*x2(3,:);
							x1(3,:).*x2(1,:); x1(3,:).*x2(2,:); x1(3,:).*x2(3,:) ]';
			s		= svd(A,0);	Cvec(f)	= s(9)*s(9);

		case{'homog','homography'}
			x1	= [refvec; onevec];
			x2	= [tgtvec; onevec];
			row = 1;
			for p = 1:npoints
				a = x1(:,p)*x2(:,p)';
				A(row:row+1,:)	= [	zerovec		-a(3,:)		a(2,:);
														a(3,:)		zerovec		-a(1,:)];
				row = row + 2;
			end
			s		= svd(A,0);	Cvec(f)	= s(9)*s(9);
	end
end

% use mean to allow number of frames considered to vary
C = mean(Cvec);


text:
% [mu,S] = mean_cov(x)
% x is n x d
function [mu,S] = mean_cov(x)

%[n,d] = size(x);

z = size(x,1);
mu = sum(x,1)/z;

diffs = bsxfun(@minus,x,mu);
S = (diffs'*diffs)/z;

author: kristinbranson

repo: JAABA


function BER = get_ber(V_equal, W_equal)

global Ns Nsym Nr Vn H hMod hDemod;
%one channel realization for Nsym times data streams transmission
% 2 is for real and imaginary
data = randi([0 1],Nsym*Ns*2,1);

%QPSK modulation generate all original signals s in Nsym times
s = reshape(step(hMod,data),Ns, Nsym);

%generate noise vector u
u = sqrt(Vn/2).*(randn(Nr,Nsym)+1i*randn(Nr,Nsym));

%get the receive vector
%colloct Nsym receive vector in the r matrix, where nth column is a
%receive vector at n time
r = W_equal' * H * V_equal * s + W_equal' * u;

%QPSK demodulation
de_data = step(hDemod, r(:));
%get the number of error bits
ber = biterr(data,de_data);
BER = ber/length(data);


text:
function checkForCharUnitsInGUI(hGUI)
%checkForCharUnitsInGUI Summary of this function goes here
%   Detailed explanation goes here
    child_handles = allchild(hGUI);
    for(i=1:length(child_handles))
        hChild = child_handles(i);
        try
            unit = get(hChild,'Units');
            if(~strcmpi(unit,'pixels') && ~strcompi(unit,'data'))
                tag = get(hChild,'Tag');
                disp(['Unit of ', tag, ' is set to ', unit, '!']);
            end
        catch ME

        end
        checkForCharUnitsInGUI(hChild);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef source < iosr.dsp.audio
%SOURCE Class of sound source separation source.
%
%   iosr.bss.source objects contain information about a sound source
%   (including its spatial location). iosr.bss.source objects are passed to
%   iosr.bss.mixture objects in order to create a mixture.
%
%   IOSR.BSS.SOURCE is a subclass of IOSR.DSP.AUDIO.
%
%   IOSR.BSS.SOURCE properties:
%       azimuth     - Azimuth of the source
%       elevation   - Elevation of the source
%       numchans    - Number of audio channels in the source
%       precomposed - Logical flag indicating whether the source should be
%                     spatialised
%
%   IOSR.BSS.SOURCE methods:
%       source      - Create a source object
%       copy        - Create an independent copy of the source (but not its
%                     audio file)
%       write       - Write the source audio file (and resample/downmix)
%
%   See also IOSR.DSP.AUDIO, IOSR.BSS.MIXTURE.

%   Copyright 2016 University of Surrey.

    properties (AbortSet)
        azimuth     % Azimuth of the source
        elevation   % Elevation of the source
        numchans    % Number of audio channels in the source
        precomposed % Logical flag indicating whether the source should be spatialised
    end

    properties (Dependent, SetAccess = protected)
        signal      % The sampled data (read-only)
    end

    properties (GetAccess = private, SetAccess = public)
         parent
    end

    methods

        % constructor
        function obj = source(filename,varargin)
        %SOURCE Create a source object
        %
        %   OBJ = IOSR.BSS.SOURCE(FILENAME) creates a source whose signal
        %   is contained in the audio file FILENAME. The source will have
        %   an azimuth and elevation of 0. The number of channels and
        %   sampling frequency will be determined by the audio file.
        %
        %   OBJ = IOSR.BSS.SOURCE(...,'PARAMETER',VALUE) allows additional
        %   options to be specified. The options are ({} indicate
        %   defaults):
        %
        %       'azimuth'       : {0} | scalar
        %           The azimuth of the source for rendering.
        %       'elevation'     : {0} | scalar
        %           The elevation of the source for rendering.
        %       'fs'            : {fs of FILENAME} | scalar
        %           The sampling frequency of the source. If the sampling
        %           frequency of the audio file FILENAME does not match,
        %           the audio file will be resampled each time the signal
        %           is requested.
        %       'precomposed'   : {false} | true
        %           Logical flag indicating whether the source should be
        %           spatialised. By default it is assumed that the source
        %           is a point source (irrespective of its channel count),
        %           and it will be spatialised. When set to true, the
        %           source will be summed directly with the spatial signal.
        %       'numchans'      : {numchans of FILENAME} | scalar
        %           The number of channels in the source. If this does not
        %           match the channel count of the audio file FILENAME, the
        %           audio file will be up-/down-mixed each time the signal
        %           is requested.
        %
        %   The FS and NUMCHANS parameters do not affect the underlying
        %   audio file, but are implemented on-the-fly each time the signal
        %   is requested. To render the changes to a new audio file, use
        %   the write() method.
        %
        %   OBJ = IOSR.BSS.SOURCE creates an empty source object.
        %
        %   Note that this is a handle class. Sources are hence passed by
        %   reference. Use the COPY() method to create an independent copy
        %   of the source.
        %
        %   See also IOSR.DSP.AUDIO.UP_DOWN_MIX.

            if nargin > 0

                propNames = {'azimuth','elevation','fs','precomposed','numchans'};

                obj.filename = filename;
                info = audioinfo(obj.filename);

                % defaults
                obj.azimuth = 0;
                obj.elevation = 0;
                obj.precomposed = false;
                obj.fs = info.SampleRate;
                obj.numchans = info.NumChannels;
                obj.rendered = false;

                % read parameter/value inputs
                if nargin > 1 % if parameters are specified
                    obj.set_properties(varargin,propNames)
                end

            end

        end

        function write(obj,filename)
        %WRITE Write the source audio file (and resample/downmix)
        %
        %   IOSR.BSS.SOURCE.WRITE() overwrites the source's audio file. If
        %   the source audio file's sampling rate and/or channel count do
        %   not match the object's properties, then the audio file will be
        %   resampled and/or up- or down-mixed accordingly.
        %
        %   IOSR.BSS.SOURCE.WRITE(FILENAME) writes the audio file to the
        %   specified file FILENAME and updates SOURCE.FILENAME.
        %
        %   See also IOSR.DSP.AUDIO.UP_DOWN_MIX.

            % overwrite filename if one is specified
            if exist('filename','var')==1
                obj.filename = filename;
            end

            assert(~isempty(obj.filename), 'iosr:source:invalidFile', 'SOURCE.FILENAME is empty.')

            % ensure path exists
            obj.ensure_path(obj.filename)

            % write new file and update filename
            audiowrite(obj.filename,obj.normalize(obj.signal),obj.fs);
            obj.rendered = true;

        end

        % validate properties

        % validate azimuth
        function set.azimuth(obj,val)
            assert(isscalar(val), 'iosr:source:invalidAzimuth', 'AZIMUTH must be a scalar')
            obj.azimuth = val;
            obj.property_changed('azimuth',val);
        end

        % validate elevation
        function set.elevation(obj,val)
            assert(isscalar(val), 'iosr:source:invalidElevation', 'ELEVATION must be a scalar')
            obj.elevation = val;
            obj.property_changed('elevation',val);
        end

        % validate numchans
        function set.numchans(obj,val)
            assert(isscalar(val), 'iosr:source:invalidNumchans', 'NUMCHANS must be a scalar')
            obj.numchans = val;
            obj.property_changed('numchans',val);
        end

        % validate precomposed
        function set.precomposed(obj,val)
            assert(islogical(val) && numel(val)==1, 'iosr:source:invalidPrecomposed', 'PRECOMPOSED property must be true or false')
            obj.precomposed = val;
            obj.property_changed('precomposed',val);
        end

        % dependent properties

        % get dependent property signal
        function signal = get.signal(obj)

            % read audio file
            [signal,fs2] = audioread(obj.filename);

            % do more stuff if differences not rendered
            if ~obj.rendered
                N = size(signal,2);

                if fs2~=obj.fs % resample
                    signal = resample(signal,obj.fs,fs2);
                end
                if N~=obj.numchans % downmix
                    signal = iosr.dsp.audio.up_down_mix(signal,obj.numchans);
                end
            end

        end

    end

    methods(Access = protected)

        function property_changed(obj,name,val)
        %PROPERTY_CHANGED handle property changes

            obj.rendered = false;

            switch lower(name)
                case 'fs'
                    if isa(obj.parent,'mixture') && obj.parent.fs~=val
                        obj.parent.fs = val;
                    end
            end

        end

        function cpObj = copyElement(obj)
        %COPYELEMENT Overload copy method with additional functionality

            % Make a shallow copy of all properties
            cpObj = copyElement@iosr.dsp.audio(obj);

        end

    end

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
function f = p03_f ( m, c, w, n, x )

%*****************************************************************************80
%
%% P03_F evaluates the function for problem p03.
%
%  Discussion:
%
%    f(x) = 1 / ( 1 + sum ( c(1:m) * x(1:m) ) ) ^ ( m + 1 )
%
%    Default values are:
%
%    c(1:m) = 1/m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alan Genz,
%    A Package for Testing Multiple Integration Subroutines,
%    in Numerical Integration: Recent Developments, Software
%    and Applications,
%    edited by Patrick Keast and Graeme Fairweather,
%    Reidel, 1987, pages 337-340,
%    ISBN: 9027725144,
%    LC: QA299.3.N38.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, real C(M,1), W(M,1), the problem parameters.
%
%    Input, integer N, the number of points.
%
%    Input, real X(M,N), the evaluation points.
%
%    Output, real F(N,1), the function values.
%
  f = ones ( n, 1 );

  for i = 1 : m
    f = f + c(i) .* x(i,:)';
  end
  f = 1.0 ./ f .^ ( m + 1 );

  return
end


text:
%computes the key of the input audio (super simple variant)
%>
%> @param x: time domain sample data, dimension samples X channels
%> @param f_s: sample rate of audio data
%> @param afWindow: FFT window of length iBlockLength (default: hann), can be [] empty
%> @param iBlockLength: internal block length (default: 4096 samples)
%> @param iHopLength: internal hop length (default: 2048 samples)
%>
%> @retval cKey key string
% ======================================================================
function [cKey] = ComputeKey (x, f_s, afWindow, iBlockLength, iHopLength)

    % set default parameters if necessary
    if (nargin < 5)
        iHopLength = 2048;
    end
    if (nargin < 4)
        iBlockLength = 4096;
    end

    if (nargin < 3 || isempty(afWindow))
        afWindow = hann(iBlockLength, 'periodic');
    end

    % key names
    cMajor  = char ('C Maj','C# Maj','D Maj','D# Maj','E Maj','F Maj',...
        'F# Maj','G Maj','G# Maj','A Maj','A# Maj','B Maj');
    cMinor  = char ('c min','c# min','d min','d# min','e min','f min',...
        'f# min','g min','g# min','a min','a# min','b min');

    % template pitch chroma (Krumhansl major/minor)
    t_pc = [6.35 2.23 3.48 2.33 4.38 4.09 2.52 5.19 2.39 3.66 2.29 2.88
            6.33 2.68 3.52 5.38 2.60 3.53 2.54 4.75 3.98 2.69 3.34 3.17];
    t_pc = diag(1 ./ sum(t_pc, 2)) * t_pc;

    % compute FFT window function
    if (length(afWindow) ~= iBlockLength)
        error('window length mismatch');
    end

    % extract audio pitch chroma
    [v_pc, t] = computeFeature("SpectralPitchChroma", x, f_s, afWindow, iBlockLength, iHopLength);

    % average pitch chroma
    v_pc = mean(v_pc, 2);

    % compute manhattan distances for major and minor
    d = zeros(2,12);
    for (i = 0:11)
        d(:,i+1)= sum(abs(repmat(v_pc', 2, 1)-circshift(t_pc, [0 i])), 2);
    end
    [dist,iKeyIdx] = min(d,[],2);
    if (dist(1) < dist(2))
        cKey = deblank(cMajor(iKeyIdx(1), :));
    else
        cKey = deblank(cMinor(iKeyIdx(2), :));
    end
end

author: alexanderlerch

repo: ACA-Code


text:
function [mColormap] = gui_Colormap_ReadPovRay(sFilename, nSize)

if ~exist('nSize', 'var')
  nSize = 256;
end

mData = [];

hFile = fopen(sFilename, 'r');
while ~(feof(hFile))
  sLine = fgetl(hFile);
  % Ignore unused line
  if IsUsed(sLine)
    [vDataLine] = ConvertLine(sLine);
    mData = [mData; vDataLine];
  end
end
fclose(hFile);

%[mData] = ReduceData(mData);

[mColormap] = ComposeColormap(mData, nSize);

% - Subfunction -

function [bUsed] = IsUsed(sLine)

if ((sLine(1) == '/') | (sLine(1) == 'c') | (sLine(1) == '}'))
  bUsed = false;
else
  bUsed = true;
end

% ---

function [vDataLine] = ConvertLine(sLine)

fPart  = str2double(sLine(3:10));
fRed   = str2double(sLine(24:31));
fGreen = str2double(sLine(34:41));
fBlue  = str2double(sLine(44:51));
vDataLine = [fPart fRed fGreen fBlue];

% ---

function [mData] = ReduceData(mData)

nLen = length(mData(:,1));
for nCnt = (nLen-1):-1:1
  if mData(nCnt,1) == mData(nCnt+1,1)
    mData(nCnt+1,:) = [];
  end
end

% ---

function [mColormap] = ComposeColormap(mData, nSize)

mColormap = [];
nLen = length(mData(:,1));

for nCnt = 1:(nLen-1)
  nSteps = floor(nSize * (mData(nCnt+1,1) - mData(nCnt,1)));
  if nSteps > 0
    mColormap = [mColormap; gui_Interpolate(mData(nCnt, 2:4), mData(nCnt+1,2:4), nSteps)];
  end
end

%/* color_map file created by the GIMP */
%/* http://www.gimp.org/               */
%color_map {
%	[0.000000 color rgbt <0.128028, 0.725490, 0.128028, 0.000000>]
%	[0.168005 color rgbt <0.564014, 0.862745, 0.564014, 0.000000>]
%	[0.333333 color rgbt <1.000000, 1.000000, 1.000000, 0.000000>]
%	[0.333333 color rgbt <1.000000, 1.000000, 1.000000, 0.000000>]
%	[0.500678 color rgbt <0.996078, 0.904645, 0.538908, 0.000000>]
%	[0.668022 color rgbt <0.992157, 0.809289, 0.077816, 0.000000>]
%	[0.668022 color rgbt <0.992157, 0.809289, 0.077816, 0.000000>]
%	[0.834011 color rgbt <0.970588, 0.449304, 0.174595, 0.000000>]
%	[1.000000 color rgbt <0.949020, 0.089319, 0.271374, 0.000000>]
%} /* color_map */

author: CelsoReyes


text:
function [] = showMultiSuperquadrics(x, c, arclength)

n = size(x, 1);
rng(5);
% color = rand(100, 3);

for i = 1 : n
%     c = color(i, :);
    showSuperquadrics(x(i, :), 'Color', c, ...
        'FaceAlpha', 1, 'ShowAxis', 1, ...
        'Light', false, 'Arclength', arclength);
    hold on
end
hold off
rng('shuffle')
end

author: ChirikjianLab

repo: Marching-Primitives


text:
function [xf,yf,modebmp] = imagemode(x,y,mode,dx,dy);

% Produces a properly scaled color plot of a two-dimensional
% mode.  This routine is especially useful when x and y are
% non-uniformly spaced vectors.  In this case, the mode is
% interpolated over a uniformly-spaced grid before producing
% an image plot.  The output can be directly saved to a file
% using the imwrite() function.
%
% USAGE:
%
% [xf,yf,modebmp] = imagemode(x,y,mode);
% [xf,yf,modebmp] = imagemode(x,y,mode,dx,dy);
%
% INPUT:
%
% x,y - vectors describing horizontal and vertical grid points
% mode - the mode or field component to be plotted
% dx, dy (optional) - fine grid spacing at which to oversample
%   (interpolate) the mode.  If left unspecified, this routine
%   will use the smallest value of diff(x) and diff(y).
%
% OUTPUT:
%
% xf,yf - points at which the mode was interpolated
% modebmp - 8-bit unsigned integer array representing the mode
%    image

x = real(x);
y = real(y);

if (size(mode) == [length(x)-1,length(y)-1])
    x = (x(1:end-1) + x(2:end))/2;
    y = (y(1:end-1) + y(2:end))/2;
end

if (nargin == 3)
  [dx,ix] = min(diff(x));
  [dy,iy] = min(diff(y));
  xf = (min(x):dx:max(x))';
  yf = (min(y):dy:max(y));
  % line up with finest portion of grid
  delta = dx*(interp1(xf,(1:length(xf)),x(ix+1)) - ...
              round(interp1(xf,(1:length(xf)),x(ix+1))));
  xf = xf + delta;
  delta = dy*(interp1(yf,(1:length(yf)),y(iy+1)) - ...
              round(interp1(yf,(1:length(yf)),y(iy+1))));
  yf = yf + delta;
  % eliminate points outside of range
  kv = find((min(x) < xf) & (xf < max(x)));
  xf = xf(kv);
  kv = find((min(y) < yf) & (yf < max(y)));
  yf = yf(kv);
else
  xf = (min(x):dx:max(x))';
  yf = (min(y):dy:max(y));
end

cmax = size(colormap,1)-1;

modebmp = uint8(transpose(interp2(y,x, ...
                abs(cmax*mode),yf,xf)));
image(xf,yf,modebmp);
set(gca,'YDir','normal');
v = [min(xf),max(xf),min(yf),max(yf)];
axis(v);
set(gca,'PlotBoxAspectRatio',[v(2)-v(1) v(4)-v(3) 1]);


text:
function ADEM_occlusion
% Slow pursuit and occlusion under active inference:
%__________________________________________________________________________
% This demo illustrates slow pursuit in the context of visual occlusion. We
% start with a simulation of canonical slow pursuit of a visual target
% with sine wave motion. Crucially, the generative model is equipped with
% a simple empirical prior encoding the hidden motion of the target (using
% a linear oscillator, whose frequency is determined by a hidden cause).
% We then examine failures of tracking and anticipation during occlusion
% and when the target re-emerges from behind the occluder. We look at a
% simulation in which the precision of the oscillator dynamics modelling
% long-term behaviour of the target is reduced (cf., neuromodulatory
% deficits in cortical areas encoding biological motion). This has a
% an effect of producing a failure of pursuit, resulting in a catch-up
% saccade on target reappearance. The suppression of prior precision can
% however have beneficial effects when motion is itself unpredicted
% (as shown with differential pursuit performance under a reversal of
% the trajectory towards the end of motion). Finally, we look at how prior
% beliefs are acquired during exposure to the target - in terms of
% cumulative inference on the hidden causes encoding the frequency of
% periodic motion. This can be regarded as a high order form of evidence
% accumulation. Importantly, this (experience-dependent) inference is
% markedly impaired by the simulated lesion to precision above. In other
% words, a single failure of inference in modelling the motion of hidden
% states can have secondary consequences - such as a failure to even
% register and remember regularities. All these simulations are based upon
% active inference; with the prior belief that the centre of gaze is
% attracted to the same point responsible for target motion.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: ADEM_occlusion.m 7679 2019-10-24 15:54:07Z spm $

% hidden causes and states
%==========================================================================
% x    - hidden states:
%   x.o(1) - oculomotor angle
%   x.o(2) - oculomotor velocity
%   x.x(1) - target angle - extrinsic coordinates
%
% v    - causal states: force on target
%
% g    - sensations:
%   g(1) - oculomotor angle (proprioception)
%   g(2) - oculomotor velocity
%   g(:) - visual input - intrinsic coordinates
%--------------------------------------------------------------------------

% Set-up
%==========================================================================
M(1).E.s = 1/2;                               % smoothness
M(1).E.n = 4;                                 % order of
M(1).E.d = 1;                                 % generalised motion

% angular frequency of target motion
%--------------------------------------------------------------------------
w  = 2*pi/32;

% sensory mappings with and without occlusion
%--------------------------------------------------------------------------
g  = '[x.o; exp(-([-8:8]'' - x.x + x.o(1)).^2)*(x.x < 1/2)]';
h  = '[x.o; exp(-([-8:8]'' - x.x + x.o(1)).^2)]';

% oculomotor latencies (sinusoidal movement)
%==========================================================================
% Endow the model with internal dynamics (a simple oscillator) so that is
% recognises and remembers the trajectory to anticipate jumps in rectified
% sinusoidal motion. First, demonstrate canonical pursuit under occlusion:

% slow pursuit following with (second order) generative model
%--------------------------------------------------------------------------
x.o = [0;0];                                  % motor angle & velocity
x.x = 0;                                      % target location

% level 1: Displacement dynamics and mapping to sensory/proprioception
%--------------------------------------------------------------------------
M(1).f = '[x.o(2); (v - x.o(1))/4 - x.o(2)/2; v - x.x]';
M(1).g = g;
M(1).x = x;                                   % hidden states
M(1).V = exp(4);                              % error precision
M(1).W = exp(4);                              % error precision

% level 2: With hidden (memory) states
%--------------------------------------------------------------------------
M(2).f  = '[x(2); -x(1)]*v/8';
M(2).g  = 'x(1)';
M(2).x  = [0; 0];                             % hidden states
M(2).V  = exp(4);                             % error precision
M(2).W  = exp(4);                             % error precision

% level 3: Encoding frequency of memory states (U)
%--------------------------------------------------------------------------
M(3).v = 0;
M(3).V = exp(4);

% generative model
%==========================================================================

% first level
%--------------------------------------------------------------------------
G(1).f = '[x.o(2); a/4 - x.o(2)/8; v - x.x]';
G(1).g = g;
G(1).x = x;                                  % hidden states
G(1).V = exp(16);                            % error precision (errors)
G(1).W = exp(16);                            % error precision (motion)
G(1).U = sparse(1,[1 2],[1 1],1,19)*exp(4);  % motor gain

% second level
%--------------------------------------------------------------------------
G(2).v = 0;                                  % exogenous force
G(2).a = 0;                                  % action force
G(2).V = exp(16);

% Sine wave cause
%--------------------------------------------------------------------------
N      = 64;                                 % length of data sequence
dt     = 16;                                 % time step (ms)
t      = (1:N)*dt;                           % PST

DEM.M  = M;
DEM.G  = G;
DEM.C  = sin((1:N)*w).*((1:N) > 16);         % sinusoidal target motion
DEM.U  = zeros(1,N) + w*8;                   % prior beliefs
DEM    = spm_ADEM(DEM);

spm_figure('GetWin','Figure 1');
spm_DEM_qU(DEM.qU,DEM.pU)
subplot(3,2,1), title({'Slow pursuit:', 'prediction and error'},'FontSize',16)
subplot(3,2,2), title({'Occluded motion:', 'hidden states'},'FontSize',16)

% create movie in extrinsic and intrinsic coordinates
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 2'); clf
spm_dem_occlusion_movie(DEM)

% repeat with simulated lesion to precision
%==========================================================================
SEM        = DEM;
SEM.M(2).W = exp(0);
SEM        = spm_ADEM(SEM);
spm_DEM_qU(SEM.qU,SEM.pU)

spm_figure('GetWin','Figure 3'); clf
spm_dem_occlusion_movie(SEM)
subplot(2,2,3), title({'Angular position:','reduced precision'},'FontSize',16)
subplot(2,2,4), title({'Angular velocity:','and catch-up saccade'},'FontSize',16)

% show improved tracking of unexpected trajectories under reduced precision
%==========================================================================

% remove occlusion and switch target trajectory after one cycle
%--------------------------------------------------------------------------
i          = 50;
DEM.M(1).g = h;
DEM.G(1).g = h;
DEM.C(i:N) = -DEM.C(i:N);

% reduce precisions and integrate
%--------------------------------------------------------------------------
SEM        = DEM;
SEM.M(2).W = exp(0);

DEM        = spm_ADEM(DEM);
SEM        = spm_ADEM(SEM);

spm_figure('GetWin','Figure 4'); clf
spm_dem_occlusion_movie(DEM)
subplot(2,2,3), hold on, subplot(2,2,4), hold on
spm_dem_occlusion_movie(SEM)
subplot(2,2,3), hold off, subplot(2,2,4), hold off

% illustrate inference on hidden cause (motion of target)
%==========================================================================

% allow for uncertainty about hidden causes (frequency of motion)
%--------------------------------------------------------------------------
DEM.M(3).V = exp(-4);

% remove occlusion
%--------------------------------------------------------------------------
DEM.M(1).g = h;
DEM.G(1).g = h;

% create a longer stimulus and reduce prior expectation
%--------------------------------------------------------------------------
N      = 128;
DEM.C  = sin((1:N)*w).*((1:N) > 16);
DEM.U  = zeros(1,N) + w/8;
DEM    = spm_ADEM(DEM);

spm_figure('GetWin','Figure 5');
spm_DEM_qU(DEM.qU,DEM.pU)
subplot(3,2,5), hold on, plot([1 N],[w w]*8,'-.k','LineWidth',4), hold off

% repeat with simulated lesion to precision
%==========================================================================
SEM        = DEM;
SEM.M(2).W = exp(0);
SEM        = spm_ADEM(SEM);
spm_DEM_qU(SEM.qU,SEM.pU)

spm_figure('GetWin','Figure 6'); clf
spm_DEM_qU(SEM.qU,SEM.pU)
subplot(3,2,4), title({'hidden states','with reduced precision'},'FontSize',16)
subplot(3,2,5), hold on, plot([1 N],[w w]*8,'-.k','LineWidth',4), hold off
axis([1 N -1/2 2])

author: spm


text:
function variance = normal_truncated_ab_variance ( mu, s, a, b )

%*****************************************************************************80
%
%% NORMAL_TRUNCATED_AB_VARIANCE returns the variance of the truncated Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real MU, S, the mean and standard deviation of the
%    parent Normal distribution.
%
%    Input, real A, B, the lower and upper truncation limits.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  alpha = ( a - mu ) / s;
  beta = ( b - mu ) / s;

  alpha_pdf = normal_01_pdf ( alpha );
  beta_pdf = normal_01_pdf ( beta );

  alpha_cdf = normal_01_cdf ( alpha );
  beta_cdf = normal_01_cdf ( beta );

  variance = s * s * ( 1.0 ...
    + ( alpha * alpha_pdf - beta * beta_pdf ) / ( beta_cdf - alpha_cdf ) ...
    - ( ( alpha_pdf - beta_pdf ) / ( beta_cdf - alpha_cdf ) ) ^ 2 );

  return
end

author: johannesgerer


text:
function stat = ft_statfun_bayesfactor(cfg, dat, design)

% FT_STATFUN_BAYESFACTOR computes the Bayes factor for a H0 of the data in two
% conditions having the same mean, versus H1 of the data having different means. This
% function supports both unpaired and paired designs and assumes flat priors.
%
% Lee and Wagenmakers (2013) provide these guidelines for its interpretation
%   IF B10 IS...    THEN YOU HAVE...
%     > 100           Extreme evidence for H1
%     30 – 100        Very strong evidence for H1
%     10 – 30         Strong evidence for H1
%     3 – 10          Moderate evidence for H1
%     1 – 3           Anecdotal evidence for H1
%     1               No evidence
%     1/3 – 1         Anecdotal evidence for H0
%     1/3 – 1/10      Moderate evidence for H0
%     1/10 – 1/30     Strong evidence for H0
%     1/30 – 1/100    Very strong evidence for H0
%     < 1/100         Extreme evidence for H0
%
% Use this function by calling one of the high-level statistics functions as
%   [stat] = ft_timelockstatistics(cfg, timelock1, timelock2, ...)
%   [stat] = ft_freqstatistics(cfg, freq1, freq2, ...)
%   [stat] = ft_sourcestatistics(cfg, source1, source2, ...)
% with the following configuration option:
%   cfg.statistic = 'ft_statfun_bayesfactor'
%
% The experimental design is specified as:
%   cfg.ivar  = independent variable, row number of the design that contains the labels of the conditions to be compared (default=1)
%   cfg.uvar  = optional, row number of design that contains the labels of the units-of-observation, i.e. subjects or trials (default=2)
%
% The labels for the independent variable should be specified as the number 1 and 2.
% The labels for the unit of observation should be integers ranging from 1 to the
% total number of observations (subjects or trials).
%
% The cfg.uvar option is only needed for paired data, you should leave it empty
% for non-paired data.
%
% See https://www.statisticshowto.datasciencecentral.com/bayes-factor-definition/ for some background.
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS

% Copyright (C) 2020, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% ensure that the required toolbox dependency is on the path
ft_hastoolbox('bayesFactor', 1);

% set the defaults
cfg.ivar = ft_getopt(cfg, 'ivar', 1);
cfg.uvar = ft_getopt(cfg, 'uvar', []); % default is empty, which means non-paired

if isempty(cfg.uvar)
  sel1 = find(design(cfg.ivar,:)==1); % select replications that belong to condition 1
  sel2 = find(design(cfg.ivar,:)==2); % select replications that belong to condition 2
  n1 = length(sel1);
  n2 = length(sel2);

  x1 = dat(:,sel1);
  x2 = dat(:,sel2);

  nchan = size(x1, 1);
  bf10  = nan(nchan, 1);
  prob  = nan(nchan, 1);
  ci_lo = nan(nchan, 1);
  ci_hi = nan(nchan, 1);
  tstat = nan(nchan, 1);
  diff  = mean(x1, 2) - mean(x2, 2); % this is the raw effect size

  for i=1:nchan
    [bf10(i), prob(i), CI, stats] = bf.ttest2(x1(i,:), x2(i,:));
    ci_lo(i) = CI(1);
    ci_hi(i) = CI(2);
    tstat(i) = stats.tstat;
  end

else
  subj = unique(design(cfg.uvar,:)); % it can also be paired over trials

  n = length(subj);
  sel1 = nan(size(subj));
  sel2 = nan(size(subj));
  for i=1:n
    sel1(i) = find(design(cfg.uvar,:)==subj(i) & design(cfg.ivar,:)==1);
    sel2(i) = find(design(cfg.uvar,:)==subj(i) & design(cfg.ivar,:)==2);
  end
  x1 = dat(:,sel1);
  x2 = dat(:,sel2);

  nchan = size(x1, 1);
  bf10  = nan(nchan, 1);
  prob  = nan(nchan, 1);
  ci_lo = nan(nchan, 1);
  ci_hi = nan(nchan, 1);
  tstat = nan(nchan, 1);
  diff  = mean(x1 - x2, 2); % this is the raw effect size

  for i=1:nchan
    [bf10(i), prob(i), CI, stats] = bf.ttest(x1(i,:) - x2(i,:));
    ci_lo(i) = CI(1);
    ci_hi(i) = CI(2);
    tstat(i) = stats.tstat;
  end

end % if uvar

% return the results for all channel-time-frequency points, for all voxels, or for all source locations
stat.bf10  = bf10(:);
stat.prob  = prob(:);
stat.diff  = diff(:);
stat.ci_lo = ci_lo(:);
stat.ci_hi = ci_hi(:);
stat.tstat = tstat(:);


text:
function res = adjSPIRiT(kData, kernel)

%  res = adjSPIRiT(kData, kernel)
%
%	backprojection of the SPIRiT operator on the data
%
%	input:
%		kData: k-space data (sx x sy x ncoils)
%		kernel : SPIRiT kernel
%
%	(c) Michael Lustig 2006

	kernel = conj(kernel(end:-1:1,end:-1:1,:));
	res = kData*0;

	numCoils = size(kernel,3);
	for n=1:numCoils
		res(:,:,n) = filter2(kernel(:,:,n), kData(:,:),'same');
	end

author: thomaskuestner


text:
function cvx_optpnt = complex_lorentz( sx, dim )

%COMPLEX_LORENTZ   Complex second-order cone.
%   COMPLEX_LORENTZ(N), where N is a positive integer, creates a column
%   variable of length N and a scalar variable, and constrains them
%   to lie in a second-order cone. That is, given the declaration
%       variable x(n) complex
%       variable y
%   the constraint
%       {x,y} == complex_lorentz(n)
%   is equivalent to
%       norm(x,2) <= y
%   The inequality form is more natural, and preferred in most cases. But
%   in fact, the COMPLEX_LORENTZ set form is used by CVX itself to convert
%   complex NORM()-based constraints to solvable form.
%
%   COMPLEX_LORENTZ(SX,DIM), where SX is a valid size vector and DIM is a
%   positive integer, creates an array variable of size SX and an array
%   variable of size SY (see below) and applies the second-order cone
%   constraint along dimension DIM. That is, given the declarations
%       sy = sx; sy(min(dim,length(sx)+1))=1;
%       variable x(sx) complex
%       variable y(sy)
%   the constraint
%       {x,y} == complex_lorentz(sx,dim)
%   is equivalent to
%       norms(x,2,dim) <= y
%   Again, the inequality form is preferred, but CVX uses the set form
%   internally. DIM is optional; if it is omitted, the first non-singleton
%   dimension is used.
%
%   LORENTZ(SX,DIM,CPLX) creates real second-order cones if CPLX is FALSE,
%   and complex second-order cones if CPLX is TRUE. The latter case is
%   equivalent to COMPLEX_LORENTZ(SX,DIM).
%
%   Disciplined convex programming information:
%       LORENTZ is a cvx set specification. See the user guide for
%       details on how to use sets.

error( nargchk( 1, 2, nargin ) );
if nargin == 1,
    cvx_optpnt = lorentz( sx, [], true );
else
    cvx_optpnt = lorentz( sx, dim, true );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [dat, volInfo, masked_vol] = iimg_mask(maskindx,dat,varargin)
% Masks index image vector or filenames in dat
% With img vector or filenames in maskindx
% * Checks dimensions to make sure origins and vox. sizes are same
%
% :Usage:
% ::
%
%     [masked_indx, volInfo, masked_vol] = iimg_mask(maskindx,dat,[volInfo], [outname])
%
% Fastest way to mask an image:
% ::
%
%    [masked_indx_data, volInfo] = iimg_mask(maskindx,dat)
%
%    % or
%
%    [masked_indx, xyz, masked_vol] = iimg_mask(mask filename,image filenames)
%
% even slower: reconstruct a 3-D volume and return in masked_vol

[volInfo, maskindx] = iimg_read_img(maskindx,1);
[volInfod, dat] = iimg_read_img(dat);

% input volInfo overrides; in case maskindx is data vector instead of
% filename
if length(varargin) > 0 && ~isempty(varargin{1}), volInfo = varargin{1}; end

% deal with multiple images
% --------------------------------------------
[nvox,nimgs] = size(dat);

% reduce dat, if maskindx is reduced to in-mask only
% if both are image files, this accomplishes the masking
% --------------------------------------------
if nvox == volInfo.nvox && size(maskindx,1) ~= volInfo.nvox
    dat = dat(volInfo.wh_inmask,:);
end

% check sizes
% --------------------------------------------
if size(maskindx,1) - size(dat,1)
    error('Sizes of mask and image data do not match.');
end

if isstruct(volInfo) && isstruct(volInfod)
    % voxel sizes and origins
    m = [diag(volInfo.mat(1:3,1:3))' volInfo.mat(1:3,4)'];
    m2 = [diag(volInfod.mat(1:3,1:3))' volInfod.mat(1:3,4)'];

    if any(m-m2)
        error('Voxel sizes or origins for mask and image data do not match.');
    end
end

% masking
% --------------------------------------------
dat(~maskindx) = 0;

if nargout > 2
    % slower, reconstruct mask
    %if nargin < 3, error('You must enter volume structure V to reconstruct 3-D vol');, end

    for i = 1:nimgs
        if length(varargin) > 1
            outname = varargin{2};
            masked_vol(:,:,:,i) = iimg_reconstruct_3dvol(dat(:,i),volInfo, 'outname', outname);

        else
            masked_vol(:,:,:,i) = iimg_reconstruct_3dvol(dat(:,i),volInfo);
        end

    end
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ a, rcond, info ] = cpoco ( a, lda, n )

%*****************************************************************************80
%
%% CPOCO factors a complex hermitian positive definite matrix.
%
%  Discussion:
%
%    The routine also estimates the condition of the matrix.
%
%    If RCOND is not needed, CPOFA is slightly faster.
%
%    To solve A*X = B, follow CPOCO by CPOSL.
%
%    To compute inverse(A)*C, follow CPOCO by CPOSL.
%
%    To compute determinant(A), follow CPOCO by CPODI.
%
%    To compute inverse(A), follow CPOCO by CPODI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex A(LDA,N), the hermitian matrix to be factored.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, complex A(LDA,N); an upper triangular matrix R so that
%      A = hermitian(R)*R
%    where hermitian(R) is the conjugate transpose.  The strict lower
%    triangle is unaltered.  If INFO ~= 0, the factorization is not complete.
%
%    Output, real RCOND, an estimate of the reciprocal condition of
%    the matrix.  For the system A*X = B, relative perturbations in A and B
%    of size EPSILON may cause relative perturbations in X of size
%    (EPSILON/RCOND).  If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate underflows.
%
%    Output, integer INFO.
%    0, for normal return.
%    K, signals an error condition.  The leading minor of order K is not
%    positive definite.
%
%  Local Parameters:
%
%    Local, complex Z(N), a work vector whose contents are usually
%    unimportant.  If A is close to a singular matrix, then Z is an
%    approximate null vector in the sense that
%      norm(A*Z) = RCOND * norm(A) * norm(Z).
%

%
%  Find norm of A using only upper half.
%
  for j = 1 : n

    z(j) = scasum ( j, a(1:j,j), 1 );

    for i = 1 : j - 1
      z(i) = real ( z(i) ) + cabs1 ( a(i,j) );
    end

  end

  anorm = 0.0;
  for j = 1 : n
    anorm = max ( anorm, real ( z(j) ) );
  end
%
%  Factor.
%
  [ a, info ] = cpofa ( a, lda, n );

  if ( info ~= 0 )
    return;
  end
%
%  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
%
%  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
%
%  The components of E are chosen to cause maximum local
%  growth in the elements of W where hermitian(R)*W = E.
%
%  The vectors are frequently rescaled to avoid overflow.
%
%  Solve hermitian(R)*W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  for k = 1 : n

    if ( cabs1 ( z(k) ) ~= 0.0 )
      ek = csign1 ( ek, -z(k) );
    end

    if ( real ( a(k,k) ) < cabs1 ( ek - z(k) ) )
      s = real ( a(k,k) ) / cabs1 ( ek - z(k) );
      z(1:n) = z(1:n) * s;
      ek = s * ek;
    end

    wk = ek - z(k);
    wkm = -ek - z(k);
    s = cabs1 ( wk );
    sm = cabs1 ( wkm );
    wk = wk / a(k,k);
    wkm = wkm / a(k,k);
    kp1 = k + 1;

    if ( kp1 <= n )

      for j = kp1 : n
        sm = sm + cabs1 ( z(j) + wkm * conj ( a(k,j) ) );
        z(j) = z(j) + wk * conj ( a(k,j) );
        s = s + cabs1 ( z(j) );
      end

      if ( s < sm )
        t = wkm - wk;
        wk = wkm;
        z(kp1:n) = z(kp1:n) + t * conj ( a(k,kp1:n) );
      end

    end

    z(k) = wk;

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
%
%  Solve R * Y = W.
%
  for k = n : -1 : 1

    if ( real ( a(k,k) ) < cabs1 ( z(k) ) )
      s = real ( a(k,k) ) / cabs1 ( z(k) );
      z(1:n) = z(1:n) * s;
    end

    z(k) = z(k) / a(k,k);
    t = -z(k);
    z(1:k-1) = z(1:k-1) + t * transpose ( a(1:k-1,k) );

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = 1.0;
%
%  Solve hermitian(R) * V = Y.
%
  for k = 1 : n

    z(k) = z(k) - z(1:k-1) * conj ( a(1:k-1,k) );

    if ( real ( a(k,k) ) < cabs1 ( z(k) ) )
      s = real ( a(k,k) ) / cabs1 ( z(k) );
      z(1:n) = z(1:n) * s;
      ynorm = s * ynorm;
    end

    z(k) = z(k) / a(k,k);

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = s * ynorm;
%
%  Solve R * Z = V.
%
  for k = n: -1 : 1

    if ( real ( a(k,k) ) < cabs1 ( z(k) ) )
      s = real ( a(k,k) ) / cabs1 ( z(k) );
      z(1:n) = z(1:n) * s;
      ynorm = s * ynorm;
    end

    z(k) = z(k) / a(k,k);
    t = -z(k);
    z(1:k-1) = z(1:k-1) + t * transpose ( a(1:k-1,k) );

  end
%
%  Make ZNORM = 1.
%
  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end

repo: jburkardt-m


text:
function triangle_node = sphere_grid_t3 ( lat_num, long_num )

%*****************************************************************************80
%
%% SPHERE_GRID_T3 produces a triangle grid on a sphere.
%
%  Discussion:
%
%    The point numbering system is the same used in SPHERE_GRIDPOINTS,
%    and that routine may be used to compute the coordinates of the points.
%
%    A sphere in 3D satisfies the equation:
%
%      sum ( ( P(1:DIM_NUM) - pc(1:DIM_NUM) )^2 ) = R^2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer LAT_NUM, LONG_NUM, the number of latitude
%    and longitude lines to draw.  The latitudes do not include the North
%    and South poles, which will be included automatically, so LAT_NUM = 5,
%    for instance, will result in points along 7 lines of latitude.
%
%    Output, integer TRIANGLE_NODE(3,(LAT_NUM+1)*LONG_NUM*2), the
%    triangle vertices.
%
  triangle_node = zeros ( 3, ( lat_num + 1 ) * long_num * 2 );

  triangle_num = 0;
%
%  The first row.
%
  n = 1;

  sw = 2;
  se = sw + 1;

  s_min = 2;
  s_max = long_num + 1;

  for j = 0 : long_num - 1

    triangle_num = triangle_num + 1;
    triangle_node(1:3,triangle_num) = [ sw, se, n ]';

    sw = se;

    if ( se == s_max )
      se = s_min;
    else
      se = se + 1;
    end

  end
%
%  The intermediate rows.
%
  for i = 1 : lat_num

    n_max = s_max;
    n_min = s_min;

    s_max = s_max + long_num;
    s_min = s_min + long_num;

    nw = n_min;
    ne = nw + 1;
    sw = s_min;
    se = sw + 1;

    for j = 0 : long_num - 1

      triangle_num = triangle_num + 1;
      triangle_node(1:3,triangle_num) = [ sw, se, nw ]';

      triangle_num = triangle_num + 1;
      triangle_node(1:3,triangle_num) = [ ne, nw, se ]';

      sw = se;
      nw = ne;

      if ( se == s_max )
        se = s_min;
      else
        se = se + 1;
      end

      if ( ne == n_max )
        ne = n_min;
      else
        ne = ne + 1;
      end

    end

  end
%
%  The last row.
%
  n_max = s_max;
  n_min = s_min;

  s = n_max + 1;

  nw = n_min;
  ne = nw + 1;

  for j = 0 : long_num - 1

    triangle_num = triangle_num + 1;
    triangle_node(1:3,triangle_num) = [ ne, nw, s ]';

    nw = ne;

    if ( ne == n_max )
      ne = n_min;
    else
      ne = ne + 1;
    end

  end

  return
end


text:
function [ o, x, w ] = epn_lag_02_xiu ( n )

%*****************************************************************************80
%
%% EPN_LAG_02_XIU implements the Xiu rule for region EPN_LAG.
%
%  Discussion:
%
%    The rule has order
%
%      O = N + 1.
%
%    The rule has precision P = 2.
%
%    EPN is the N-dimensional positive space [0,+oo)^N with exponential
%    or Laguerre weight function:
%
%      w(x(1:n)) = exp ( - sum ( x(1:n) ) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dongbin Xiu,
%    Numerical integration formulas of degree two,
%    Applied Numerical Mathematics,
%    Volume 58, 2008, pages 1515-1520.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Output, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  o = n + 1;

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  for j = 1 : o

    i = 0;
    for r = 1 : floor ( n / 2 )
      arg = 2 * r * ( j - 1 ) * pi / ( n + 1 );
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * cos ( arg );
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * sin ( arg );
    end

    if ( i < n )
      i = i + 1;
      x(i,j) = r8_mop ( j - 1 );
    end

  end

  gamma0 = - 1.0;
  delta0 = 1.0;
  c1 = - 1.0;

  x(1:n,1:o) = ( sqrt ( gamma0 * c1 ) * x(1:n,1:o) - delta0 ) / gamma0;

  expon = 0;
  volume_1d = ep1_lag_monomial_integral ( expon );
  volume = volume_1d ^ n;

  w(1:o) = volume / o;

  return
end


text:
function [H_x] = H_xfn(kMinus1State, imuMeasurement, deltaT )
%H_XFN Compute the 6x6 matrix H_x_k

psiVec = imuMeasurement.omega*deltaT;
psiMag = norm(psiVec);
d = imuMeasurement.v*deltaT;

Psi = cos(psiMag)*eye(3) + (1 - cos(psiMag))*(psiVec/psiMag)*(psiVec/psiMag)' - sin(psiMag)*crossMat(psiVec/psiMag);

%Construct the H_x matrix
H_x = zeros(6,6);
H_x(1:3, 1:3) = eye(3);
H_x(4:6, 4:6) = Psi;
H_x(1:3, 4:6) = -kMinus1State.C_vi'*crossMat(d);

end

author: utiasSTARS

repo: msckf-swf-comparison