text:
function [xfinal,ffinal,exitflag,xstart] = rmsearch(fun,optname,x0,LB,UB,varargin)
% RMSEARCH: Randomly multistarted overlay to several optimizers
% usage: [xfinal,ffinal,exitflag,xstart] = RMSEARCH(fun,optname,x0,LB,UB)
% usage: [xfinal,ffinal,exitflag,xstart] = RMSEARCH(fun,optname,x0,LB,UB,prop1,val1,prop2,val2,...)
%
% This tool is useful when little information is available
% for good starting values for an optimization, or when many
% local solutions exist. It is also useful to map out the various
% basins of attaction of each local solution.
%
% RMSEARCH is normally used for bounded optimization problems,
% but with care can be used when some variables are unbounded.
% Whenever possible I'd still recommend the use of both lower
% and upper bounds if you can possibly do so. Restriction of
% the solution space is an important way to improve the
% performance of any optimization.
%
%
% arguments: (Input)
%  fun - name of a function, inline function, or a function handle
%
%  optname - character string - the name of the
%        optimizer to be used. Valid options are
%
%        'fminsearchbnd', 'fminsearchcon', 'fmincon'
%        'fminbnd', 'fzero', 'lsqnonlin', 'fminsearch'
%
%        No default is allowed, but the function name may
%        be shortened as long as the shortening is unambiguous.
%        Thus 'fminb' or 'fz' are acceptable, but 'fm' is not.
%
%        Note that fzero and fminbnd, since they work in 1-d
%        will search EVERY bracketed interval they find that
%        contains either a root or a minimizer (as appropriate.)
%
%        Remember to set the options structure properly for
%        fmincon, turning off the largescale optimizer as
%        appropriate. Otherwise you will see many warning
%        messages generated, one for each call to the optimizer.
%
%  x0 -  starting value. Used to determine the dimension and
%        shape of the parameter vector. Also used to determine
%        the distribution of the random starts when the bounds
%        are incompletely specified.
%
%        No default is allowed.
%
%        If both lower and upper bounds are supplied for a
%        parameter, then the distribution for the multi-starts
%        is chosen to be uniform between those bounds. If no
%        upper bound or no lower bound, then the distribution
%        follows a two parameter exponential distribution, with
%        mean at x0(i). If both upper and lower bounds are
%        missing (i.e., they are +/- inf), then that parameter
%        will be normally distribued with mean at x0(i).
%
%  LB -  Lower bounds for the parameters, as used by fmincon,
%        fminsearchbnd, lsqnonlin, etc.
%
%        LB must be the same size as x0, and LB(i)<=UB(i)
%
%        When no lower bound exists for a variable, use -inf
%
%  UB -  Upper bounds for the parameters, as used by fmincon,
%        fminsearchbnd, lsqnonlin, etc.
%
%        UB must be the same size as x0, and LB(i)<=UB(i)
%
%        When no upper bound exists for a variable, use +inf
%
%
% Additional inputs must be in the form of property/value pairs.
%   Properties are character strings. They may be shortened
%   to the extent that they are unambiguous. Properties are
%   not case sensitive. Valid property names are:
%
%    'Plot', 'Nonlcon', 'FractionUsed', 'InitialSample',
%    'OverSample', 'A', 'B' 'Options'
%
%   All properties have default values, chosen as intelligently
%   as I could manage. Values that are character strings may
%   also be unambiguously shortened. The legal values for each
%   property are:
%
%   'Plot' - specifies if a plot is generated {'on', 'off'}
%         A plot can only be generated for 1 or 2 dimensional
%         problems. A warning is generated for more than 2
%         variables if plot is turned on.
%
%         DEFAULT: 'off'
%
%   'Nonlcon' - Nonlinear constraint function or function handle
%         Only used for 'fmincon' or 'fminsearchcon'. Other
%         optimizers will cause an error if this is provided.
%
%         DEFAULT: ''
%
%   'FractionUsed' - Fractional amount of the total sample that
%         are used as multiple starts for the optimizer. It must
%         lie in the half open interval (0,1].
%
%         DEFAULT: 1  (100% of the samples are used as starting
%         points.
%
%         Note: FractionUsed is ignored for fminbnd and fzero.
%
%   'InitialSample' - Size of the total initial random sample
%         generated.
%
%         DEFAULT: 100  (100 points are generated)
%
%   'OverSample' - Over-sampling ratio. If linear inequality or
%         or nonlinear inequality constraints are involved, a
%         simple sampling strategy within the bounds may fail
%         to produce as many feasible points as are requested.
%
%         Thus if 100 points are requested in the initial sample,
%         a 20 to 1 oversampling ratio will generate a total of
%         2000 initial samples that are all within the simple bound
%         constraints. All of these points are compared to the
%         inequality constraints, rejecting all of those that fail.
%         The first InitialSample set of points are retained, only
%         these points are considered as starting points for an
%         optimization.
%
%         DEFAULT: 20  (a 20 to 1 over-sampling ratio)
%
%   'A', 'B' - Linear Inequality constraint matrix (A) and right
%         hand side vector (B), used only for 'fmincon' or
%         'fminsearchcon'. Other optimizers will cause an error
%         if these are provided.
%
%   'Options' - options structure for the specified optimizer.
%
%        default: optimset('optname')
%
%        For fminsearchbnd or fminsearchcon, the default is
%        optimset('fminsearch')
%
% arguments: (Output)
%  xfinal - final solutions found for each xstart point. Each
%         row of xfinal is one solution.
%
%  ffinal - objective function at xfinal.
%
%         For lsqnonlin or lsqcurvefit, ffinal is the sum
%         of squares.
%
%  exitflag - exitflags for each "solution" found
%
%  xstart - starting value used for each point in xfinal
%
%
% Example usage:
%  Find all zeros of besselj(2,x) betweeen x=0 and x=50. Note
%  that x0 is actually ignored in this case. Generate a plot of
%  the sample points and the solutions found.
%
%  fun = @(x) besselj(2,x);
%  [xfinal,ffinal,exitflag,xstart] = rmsearch(fun,'fzero',...
%       10,0,50,'initialsample',500,'plot','on')
%
% Example usage:
%  Find local minimizers of besselj(2,x) betweeen x=0 and x=200
%  Generate a plot of the sample points and the solutions found.
%
%  fun = @(x) besselj(2,x);
%  [xfinal,ffinal,exitflag,xstart] = rmsearch(fun,'fminbnd',...
%       10,0,200,'initialsample',500,'plot','on')
%
% Example usage:
%  Find local minimizers of the peaks function. Generate
%  a plot of the sample points and the solutions found. Use
%  fminsearchbnd as the optimizer.
%
%  fun = inline('peaks(x(1),x(2))','x');
%  [xfinal,ffinal,exitflag,xstart] = rmsearch(fun,'fminsearchbnd',...
%      [0 0],[-5 -5],[5 5],'initialsample',100,'plot','on')
%
% Example usage:
%  Find the local minimizers of the peaks function within a circle
%  of radius 2 around the origin. Generate a plot of the sample points
%  and the solutions found. Use fmincon as the optimizer, setting the
%  appropriate options for its use. Note the trick with deal for nonlcon
%  inside an anonymous function.
%
%  fun = inline('peaks(x(1),x(2))','x');
%  opts = optimset('fmincon');
%  opts.LargeScale = 'off';
%  opts.Display = 'none';
%  nonlcon = @(x) deal(norm(x) - 2,[]);
%  [xfinal,ffinal,exitflag,xstart] = rmsearch(fun, ...
%     'fmincon',[0 0],[-2 -2],[2 2],'initialsample',100, ...
%     'plot','on','options',opts,'nonlcon',nonlcon)
%
%
% See also: fminsearchbnd, fmincon, fminsearchcon, fzero, fminbnd, lsqnonlin
%
%
% Author: John D'Errico
% E-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 1/1/07

% set defaults
par.Plot = 'off';
par.Nonlcon = '';
par.FractionUsed = 1;
par.InitialSample = 100;
par.OverSample = 20;
par.A = [];
par.B = [];
par.OptName = optname;
par.Options = [];

na = length(varargin);
if (rem(na,2)==1)
  error 'Property/value pairs must come as PAIRS of arguments.'
elseif na>0
  par = parse_pv_pairs(par,varargin);
end

% shape of x0, dimension of the problem
if (nargin<3) || isempty(x0)
  error 'x0 must be provided.'
else
  par.Nx0 = size(x0);
  x0 = x0(:);
  par.Dim = prod(par.Nx0);
end

% checks for validity of parameters, any problems
par = check_params(par);

% Was fun a string, or an inline/anonymous function?
if (nargin<1)
  help rmsearch
  return
elseif isempty(fun)
  error 'fun was not supplied.'
elseif ischar(fun)
  % a character function name
  fun = str2func(fun);
end

% Check the bounds
if (nargin<4) || isempty(LB)
  LB = repmat(-inf,par.Dim,1);
else
  LB = LB(:);
end
if (nargin<5) || isempty(UB)
  UB = repmat(inf,par.Dim,1);
else
  UB = UB(:);
end
% are they the same size as X0?
if length(LB)~=par.Dim
  error 'Lower bound array was inconsistent in size with x0'
end
if length(UB)~=par.Dim
  error 'Upper bound array was inconsistent in size with x0'
end

% Were there any unusable inequalities provided?
if ~ismember(par.OptName,{'fmincon','fminsearchcon'})
  if par.UseLC
    error(['Linear inequalities provided are not usable by ',par,OptName])
  elseif par.UseNLC
    error(['Nonlinear inequalities provided are not usable by ',par,OptName])
  end
end

% are any of the bounds inconsistent?
if any(LB>UB)
  error 'It is required that: LB <= UB'
end

% Choose the distributions used to sample
% from for each variable.
LBinf=isinf(LB);
UBinf=isinf(UB);
for i = 1:par.Dim
  if LBinf(i) && UBinf(i)
    % both are unbounded, use a normal
    vardist(i).dist = 'normal';
    vardist(i).mean = x0(i);
    vardist(i).sd = 1;
  elseif ~LBinf(i) && UBinf(i)
    % only lower bound, use an exponential
    vardist(i).dist = 'lower';
    vardist(i).lambda = x0(i) - LB(i);
    if vardist(i).lambda<=0
      vardist(i).lambda = 0.5;
    else
      vardist(i).lambda = 0.5 ./(vardist(i).lambda);
    end
    vardist(i).shift = LB(i);
  elseif LBinf(i) && ~UBinf(i)
    % only upper bound, use an exponential
    vardist(i).dist = 'upper';
    vardist(i).lambda = UB(i) - x0(i);
    if vardist(i).lambda <= 0
      vardist(i).lambda = 0.5;
    else
      vardist(i).lambda = 0.5 ./(vardist(i).lambda);
    end
    vardist(i).shift = UB(i);
  else
    % dual bounds, use a uniform
    vardist(i).dist = 'uniform';
    vardist(i).a = LB(i);
    vardist(i).b = UB(i);
  end
end

% Generate the intial random samples. Choose enough
% of them in case there are constraints to worry about.
if par.UseLC || par.UseNLC
  % oversampling for rejection
  xinitial = rand(par.OverSample*par.InitialSample,par.Dim);
else
  % No need to oversample
  xinitial = rand(par.InitialSample,par.Dim);
end

% generate a sample with the specified distributions
for i = 1:par.Dim
  switch vardist(i).dist
    case 'uniform'
      % uniform between two bounds
      xinitial(:,i) = vardist(i).a + ...
        (vardist(i).b - vardist(i).a)*xinitial(:,i);
    case 'normal'
      % unbounded - use a normal
      xinitial(:,i) = vardist(i).mean + ...
        vardist(i).sd*randn(par.InitialSample,1);
    case 'lower'
      % lower bound, use an exponential
      xinitial(:,i) = vardist(i).shift - ...
        log(xinitial(:,i))/vardist(i).lambda;
    case 'upper'
      % upper bound, use an exponential
      xinitial(:,i) = vardist(i).shift + ...
        log(xinitial(:,i))/vardist(i).lambda;
  end
end

% if 1-d problem, fminbnd and fzero will both do best
% if any bounds are included in the sample
if (par.Dim == 1)
  if ~isinf(LB) && ~isinf(UB)
    xinitial = [LB;UB;xinitial];
  elseif  (~isinf(LB) && isinf(UB))
    xinitial = [LB;xinitial];
  elseif (isinf(LB) && ~isinf(UB))
    xinitial = [UB;xinitial];
  end
end

% Do we need to reject any samples?
if par.UseLC
  k = ((repmat(par.B,1,size(xinitial,1)) - par.A*xinitial') < 0);
  k = any(k,1);
  xinitial(k,:) = [];
end
if size(xinitial,1)<1
  error 'Linear inequalities too strict. No feasible points were identified.'
end
if par.UseNLC
  % test any nonlinear constraints
  nx = size(xinitial,1);
  k = true(nx,1);
  for i = 1:nx
    [cineq,ceq] = par.Nonlcon(xinitial(i,:)); %#ok
    k(i) = all(cineq>0);
  end
  % reject those that failed
  xinitial(k,:) = [];
end
if size(xinitial,1)<1
  error 'Nonlinear constraints too strict. No feasible points were identified.'
end

% Some points got through, but not as many as requested
xsize = size(xinitial,1);
if xsize < par.InitialSample
  warning('RMSEARCH:FeasiblePoints',['Only ',num2str(xsize),' feasible points were found, ', ...
    num2str(par.InitialSample),' were requested. Insufficiently Oversampled.'])
end
% keep as many samples as we need if oversampled
xinitial = xinitial(1:min(xsize,par.InitialSample),:);
xsize = size(xinitial,1);

% Evaluate each sample through our objective function.
finitial = zeros(xsize,1);
for i = 1:xsize
  if ~strcmp(par.OptName,'lsqnonlin')
    % A simple optimizer:
    % 'fminbnd', 'fzero', 'fmincon', 'fminsearchbnd',
    % 'fminsearchcon','fminsearch'
    finitial(i) = fun(reshape(xinitial(i,:),par.Nx0));
  else
    % lsqnonlin, so we need to form the sum of squares
    f = fun(reshape(xinitial(i,:),par.Nx0));
    finitial(i) = sum(f.^2);
  end
end

% At which of these sample points will we
% start the indicated optimizer? Was it a 1-d
% problem? If so, then both fzero and fminbnd
% use a bracketed interval.
xfinal = zeros(size(xinitial));
ffinal = zeros(size(finitial));
if (par.Dim==1) && ismember(par.OptName,{'fzero', 'fminbnd'})
  switch par.OptName
    case 'fzero'
      % fzero - sort on x first
      [xinitial,tags] = sort(xinitial);
      finitial = finitial(tags);

      % find zero crossings in f
      ind = 1:(xsize-1);
      k = find((sign(finitial(ind)) .* sign(finitial(ind+1))) <= 0);

      % any crossings found?
      if isempty(k)
        error 'No zero crossings found in this sample for fzero'
      end

      xfinal = zeros(length(k),1);
      ffinal = xfinal;
      exitflag = xfinal;
      xstart = xfinal;
      fstart = xfinal;
      % loop over the candidate intervals
      for i = 1:length(k)
        [xfinal(i),ffinal(i),exitflag(i)] = fzero(fun, ...
          [xinitial(k(i)),xinitial(k(i)+1)],par.Options);

        % store the better of the two points in each bracket
        if abs(finitial(k(i))) <= abs(finitial(k(i)+1))
          xstart(i) = xinitial(k(i));
          fstart(i) = finitial(k(i));
        else
          xstart(i) = xinitial(k(i)+1);
          fstart(i) = finitial(k(i)+1);
        end
      end

    case 'fminbnd'
      % fminbnd - sort on x first
      [xinitial,tags] = sort(xinitial);
      finitial = finitial(tags);

      % find local minima in f
      ind = 2:(xsize-1);
      k = 1 + find((finitial(ind) <= finitial(ind-1)) & ...
        (finitial(ind) <= finitial(ind+1)));

      % do we look in the first interval?
      if finitial(1) < finitial(2)
        k = union(2,k);
      end
      % or the last one?
      if finitial(end) < finitial(end-1)
        k = union(xsize-1,k);
      end

      % any local minima found?
      xfinal = zeros(length(k),1);
      ffinal = xfinal;
      exitflag = xfinal;
      xstart = xfinal;
      fstart = xfinal;
      % loop over the candidate intervals
      for i = 1:length(k)
        [xfinal(i),ffinal(i),exitflag(i)] = fminbnd(fun, ...
          xinitial(k(i)-1),xinitial(k(i)+1),par.Options);

        % store in xstart
        xstart(i) = xinitial(k(i));
        fstart(i) = finitial(k(i));
      end
  end
else
  % it must be a multivariable minimization, or only one
  % variable, but not with a bracketed optimizer.

  % sort on f first in increasing order
  [finitial,tags] = sort(finitial);
  xinitial = xinitial(tags,:);

  if isempty(par.FractionUsed) || (par.FractionUsed==1)
    % use all the points as starting values
    nk = xsize;
  else
    nk = min(xsize,ceil(xsize*par.FractionUsed));
  end

  % any local minima found?
  xfinal = zeros(nk,par.Dim);
  ffinal = zeros(nk,1);
  exitflag = ffinal;
  xstart = xfinal;
  fstart = ffinal;
  % loop over the candidates chosen
  for i = 1:nk
    % store the start points in xstart and
    % fstart for plotting later
    xstart(i,:) = xinitial(i,:);
    fstart(i) = finitial(i);

    % but which optimizer?
    switch par.OptName
      case 'lsqnonlin'
        [xfinal(i,:),ffinal(i),residual,exitflag(i)] = ...
          lsqnonlin(fun,xinitial(i,:),LB,UB,par.Options); %#ok
      case 'fminsearchbnd'
        [xfinal(i,:),ffinal(i),exitflag(i)] = ...
          fminsearchbnd(fun,xinitial(i,:),LB,UB,par.Options);
      case 'fminsearchcon'
        [xfinal(i,:),ffinal(i),exitflag(i)] = ...
          fminsearchcon(fun,xinitial(i,:),LB,UB, ...
          par.A,par.B,par.Nonlcon,par.Options);
      case 'fminsearch'
        [xfinal(i,:),ffinal(i),exitflag(i)] = ...
          fminsearch(fun,xinitial(i,:),par.Options);
      case 'fmincon'
        [xfinal(i,:),ffinal(i),exitflag(i)] = ...
          fmincon(fun,xinitial(i,:),par.A,par.B, ...
          [],[],LB,UB,par.Nonlcon,par.Options);
    end
  end
end

% do we plot the results?
if strcmp(par.Plot,'on') && par.Dim <= 2
  switch par.Dim
    case 1
      % its a 1-d problem
      figure
      plot([xstart';xfinal'],[fstart';ffinal'],'g-')
      hold on
      plot(xfinal,ffinal,'ro')
      plot(xinitial,finitial,'b+')
      plot([xinitial(1),xinitial(end)],[0 0],'k:')
      hold off
      xlabel 'x'
      ylabel 'fun(x)'
      title(['Random multi-started search: ',par.OptName])

    case 2
      % a 2-d problem
      figure
      C = finitial-min(finitial);
      C = C/max(C);
      colormap hsv
      scatter3(xinitial(:,1),xinitial(:,2),finitial,50,C,'+')
      hold on
      plot3([xstart(:,1)';xfinal(:,1)'], ...
        [xstart(:,2)';xfinal(:,2)'],[fstart';ffinal'],'g-')
      plot3(xfinal(:,1),xfinal(:,2),ffinal,'ro')
      hold off
      xlabel 'x'
      ylabel 'y'
      zlabel 'fun([x,y])'
      title(['Random multi-started search: ',par.OptName])

  end
end

% ============================================
%     end of mainline rmsearch
% ============================================

% ============================================
% subfunction - check_params
% ============================================
function par = check_params(par)
% check the parameters for acceptability
%
% Defaults
% par.Plot = 'off'
% par.Nonlcon = '';
% par.FractionUsed = 1;
% par.InitialSample = 100;
% par.OptName = optname;

% Nonlcon == '' by default
if ~isempty(par.Nonlcon)
  if ischar(par.Nonlcon)
    par.Nonlcon = str2func(par.Nonlcon);
  end
  par.UseNLC = true;
else
  par.UseNLC = false;
end

% FractionUsed == 1 by default
if isempty(par.FractionUsed)
  par.FractionUsed = 1;
elseif (length(par.FractionUsed)>1) || (par.FractionUsed<=0) || (par.FractionUsed>1)
  error 'FractionUsed must be scalar, in the interval (0,1]'
end

% InitialSample == 100 by default
if isempty(par.InitialSample)
  par.InitialSample = 100;
elseif (length(par.InitialSample)>1) || (par.InitialSample<=0)
  error 'InitialSample must be positive and scalar'
else
  par.InitialSample = ceil(par.InitialSample);
end

% OverSample == 20 by default
if isempty(par.OverSample)
  par.OverSample = 20;
elseif (length(par.OverSample)>1) || (par.OverSample<1)
  error 'OverSample must be positive and scalar, >= 1'
end

% OptName is char
valid = {'fmincon', 'fminsearchcon', 'fminsearchbnd', ...
  'fminbnd', 'fzero', 'lsqnonlin', 'fminsearch'};
if isempty(par.OptName) || ~ischar(par.OptName)
  error 'Invalid OptName: Must be non-empty character'
end
ind = strmatch(par.OptName,valid,'exact');
if (length(ind)==1)
  par.OptName = valid{ind};
else
  ind = strmatch(par.OptName,valid);
  if isempty(ind) || (length(ind)>1)
    error(['Invalid OptName: ',par.OptName])
  else
    par.OptName = valid{ind};
  end
end
% 2+ dimensions for a 1-d optimizer?
if (par.Dim>1) && ismember(par.OptName,{'fminbnd', 'fzero'})
  error '1-d optimizer specified, but more than 1 variable to optimize'
end

% Options is a struct
if isempty(par.Options)
  switch par.OptName
    case 'fmincon'
      par.options = optimset('fmincon');
    case 'fzero'
      par.options = optimset('fzero');
    case 'fminbnd'
      par.options = optimset('fminbnd');
    case 'lsqnonlin'
      par.options = optimset('lsqnonlin');
    case {'fminsearch', 'fminsearchbnd', 'fminsearchcon'}
      par.options = optimset('fminbnd');
  end
elseif ~isstruct(par.Options)
  error 'If provided, Options must be a struct generated by optimset'
end

% Plot is char
valid = {'on', 'off'};
if isempty(par.Plot)
  error 'Invalid Plot: Must be character or empty'
end
ind = find(strncmpi(par.Plot,valid,length(par.Plot)));
if (length(ind)==1)
  par.Plot = valid{ind};
else
  error(['Invalid Plot: ',par.Plot])
end
if strcmp(par.Plot,'on') && (par.Dim>3)
  % No plots available for 4 or more variables
  Warning 'Sorry, plots are not generated for 4 or more variables.'
  par.Plot = 'off';
end

% A & B are both empty by default
if ~isempty(par.A) && (size(par.A,2)~=par.Dim)
  error 'A must have the same number of columns as the # of variables'
end
if ~isempty(par.B) && (size(par.A,1)~=size(par.B,1))
  error 'B must have the same number of rows as A'
end
if isempty(par.A)
  par.UseLC = false;
else
  par.UseLC = true;
end

% ============================================
% Included subfunction - parse_pv_pairs
% ============================================
function params=parse_pv_pairs(params,pv_pairs)
% parse_pv_pairs: parses sets of property value pairs, allows defaults
% usage: params=parse_pv_pairs(default_params,pv_pairs)
%
% arguments: (input)
%  default_params - structure, with one field for every potential
%             property/value pair. Each field will contain the default
%             value for that property. If no default is supplied for a
%             given property, then that field must be empty.
%
%  pv_array - cell array of property/value pairs.
%             Case is ignored when comparing properties to the list
%             of field names. Also, any unambiguous shortening of a
%             field/property name is allowed.
%
% arguments: (output)
%  params   - parameter struct that reflects any updated property/value
%             pairs in the pv_array.
%
% Example usage:
% First, set default values for the parameters. Assume we
% have four parameters that we wish to use optionally in
% the function examplefun.
%
%  - 'viscosity', which will have a default value of 1
%  - 'volume', which will default to 1
%  - 'pie' - which will have default value 3.141592653589793
%  - 'description' - a text field, left empty by default
%
% The first argument to examplefun is one which will always be
% supplied.
%
%   function examplefun(dummyarg1,varargin)
%   params.Viscosity = 1;
%   params.Volume = 1;
%   params.Pie = 3.141592653589793
%
%   params.Description = '';
%   params=parse_pv_pairs(params,varargin);
%   params
%
% Use examplefun, overriding the defaults for 'pie', 'viscosity'
% and 'description'. The 'volume' parameter is left at its default.
%
%   examplefun(rand(10),'vis',10,'pie',3,'Description','Hello world')
%
% params =
%     Viscosity: 10
%        Volume: 1
%           Pie: 3
%   Description: 'Hello world'
%
% Note that capitalization was ignored, and the property 'viscosity'
% was truncated as supplied. Also note that the order the pairs were
% supplied was arbitrary.

npv = length(pv_pairs);
n = npv/2;

if n~=floor(n)
  error 'Property/value pairs must come in PAIRS.'
end
if n<=0
  % just return the defaults
  return
end

if ~isstruct(params)
  error 'No structure for defaults was supplied'
end

% there was at least one pv pair. process any supplied
propnames = fieldnames(params);
lpropnames = lower(propnames);
for i=1:n
  p_i = lower(pv_pairs{2*i-1});
  v_i = pv_pairs{2*i};

  ind = strmatch(p_i,lpropnames,'exact');
  if isempty(ind)
    ind = find(strncmp(p_i,lpropnames,length(p_i)));
    if isempty(ind)
      error(['No matching property found for: ',pv_pairs{2*i-1}])
    elseif length(ind)>1
      error(['Ambiguous property name: ',pv_pairs{2*i-1}])
    end
  end
  p_i = propnames{ind};

  % override the corresponding default in params
  params = setfield(params,p_i,v_i); %#ok

end

% end % parse_pv_pairs

author: Sable

repo: mcbench-benchmarks


text:
function [ less, equal, more ] = r8r8vec_index_search ( n, x, y, indx, ...
  xval, yval )

%*****************************************************************************80
%
%% R8R8VEC_INDEX_SEARCH searches for an R8R8 value in an indexed sorted list.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the size of the current list.
%
%    Input, real X(N), Y(N), the list.
%
%    Input, integer INDX(N), the sort index of the list.
%
%    Input, real XVAL, YVAL, the value to be sought.
%
%    Output, integer LESS, EQUAL, MORE, the indexes in INDX of the
%    entries of X that are just less than, equal to, and just greater
%    than XVAL.  If XVAL does not occur in X, then EQUAL is zero.
%    If XVAL is the minimum entry of X, then LESS is 0.  If XVAL
%    is the greatest entry of X, then MORE is N+1.
%
  if ( n <= 0 )
    less = 0;
    equal = 0;
    more = 0;
    return
  end

  lo = 1;
  hi = n;

  xlo = x(indx(lo));
  ylo = y(indx(lo));

  xhi = x(indx(hi));
  yhi = y(indx(hi));

  compare = r8r8_compare ( xval, yval, xlo, ylo );

  if ( compare == -1 )
    less = 0;
    equal = 0;
    more = 1;
    return
  elseif ( compare == 0 )
    less = 0;
    equal = 1;
    more = 2;
    return
  end

  compare = r8r8_compare ( xval, yval, xhi, yhi );

  if ( compare == 1 )
    less = n;
    equal = 0;
    more = n + 1;
    return
  elseif ( compare == 0 )
    less = n - 1;
    equal = n;
    more = n + 1;
    return
  end

  while ( 1 )

    if ( lo + 1 == hi )
      less = lo;
      equal = 0;
      more = hi;
      return
    end

    mid = round ( ( lo + hi ) / 2 );
    xmid = x(indx(mid));
    ymid = y(indx(mid));

    compare = r8r8_compare ( xval, yval, xmid, ymid );

    if ( compare == 0 )
      equal = mid;
      less = equal - 1;
      more = equal + 1;
      return
    elseif ( compare == -1 )
      hi = mid;
    elseif ( compare == +1 )
      lo = mid;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = emailFeatures(word_indices)
%EMAILFEATURES takes in a word_indices vector and produces a feature vector
%from the word indices
%   x = EMAILFEATURES(word_indices) takes in a word_indices vector and
%   produces a feature vector from the word indices.

% Total number of words in the dictionary
n = 1899;

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

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return a feature vector for the
%               given email (word_indices). To help make it easier to
%               process the emails, we have have already pre-processed each
%               email and converted each word in the email into an index in
%               a fixed dictionary (of 1899 words). The variable
%               word_indices contains the list of indices of the words
%               which occur in one email.
%
%               Concretely, if an email has the text:
%
%                  The quick brown fox jumped over the lazy dog.
%
%               Then, the word_indices vector for this text might look
%               like:
%
%                   60  100   33   44   10     53  60  58   5
%
%               where, we have mapped each word onto a number, for example:
%
%                   the   -- 60
%                   quick -- 100
%                   ...
%
%              (note: the above numbers are just an example and are not the
%               actual mappings).
%
%              Your task is take one such word_indices vector and construct
%              a binary feature vector that indicates whether a particular
%              word occurs in the email. That is, x(i) = 1 when word i
%              is present in the email. Concretely, if the word 'the' (say,
%              index 60) appears in the email, then x(60) = 1. The feature
%              vector should look like:
%
%              x = [ 0 0 0 0 1 0 0 0 ... 0 0 0 0 1 ... 0 0 0 1 0 ..];
%
%

for i = 1:size(word_indices,1),
  x(word_indices(i)) = 1;
end

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

end

author: anirudhjayaraman

repo: Machine-Learning


text:
function [ less, equal, more ] = r8r8r8vec_index_search ( n, x, y, z, ...
  indx, xval, yval, zval )

%*****************************************************************************80
%
%% R8R8R8VEC_INDEX_SEARCH searches for an R8R8R8 value in an indexed sorted list.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the size of the current list.
%
%    Input, real X(N), Y(N), Z(N), the list.
%
%    Input, integer INDX(N), the sort index of the list.
%
%    Input, real XVAL, YVAL, ZVAL, the value to be sought.
%
%    Output, integer LESS, EQUAL, MORE, the indexes in INDX of the
%    entries of X that are just less than, equal to, and just greater
%    than XVAL.  If XVAL does not occur in X, then EQUAL is zero.
%    If XVAL is the minimum entry of X, then LESS is 0.  If XVAL
%    is the greatest entry of X, then MORE is N+1.
%
  if ( n <= 0 )
    less = 0;
    equal = 0;
    more = 0;
    return
  end

  lo = 1;
  hi = n;

  xlo = x(indx(lo));
  ylo = y(indx(lo));
  zlo = z(indx(lo));

  xhi = x(indx(hi));
  yhi = y(indx(hi));
  zhi = z(indx(hi));

  compare = r8r8r8_compare ( xval, yval, zval, xlo, ylo, zlo );

  if ( compare == -1 )
    less = 0;
    equal = 0;
    more = 1;
    return
  elseif ( compare == 0 )
    less = 0;
    equal = 1;
    more = 2;
    return
  end

  compare = r8r8r8_compare ( xval, yval, zval, xhi, yhi, zhi );

  if ( compare == 1 )
    less = n;
    equal = 0;
    more = n + 1;
    return
  elseif ( compare == 0 )
    less = n - 1;
    equal = n;
    more = n + 1;
    return
  end

  while ( 1 )

    if ( lo + 1 == hi )
      less = lo;
      equal = 0;
      more = hi;
      return
    end

    mid = round ( ( lo + hi ) / 2 );
    xmid = x(indx(mid));
    ymid = y(indx(mid));
    zmid = z(indx(mid));

    compare = r8r8r8_compare ( xval, yval, zval, xmid, ymid, zmid );

    if ( compare == 0 )
      equal = mid;
      less = equal - 1;
      more = equal + 1;
      return
    elseif ( compare == -1 )
      hi = mid;
    elseif ( compare == +1 )
      lo = mid;
    end

  end

  return
end


text:
function [beta_gibbs F_gibbs L_gibbs phi_gibbs sigma_gibbs lambda_t_gibbs sigma_t_gibbs sbar]=stvol1gibbs(Xbart,yt,beta0,omega0,alpha0,delta0,f0,upsilon0,betahat,sigmahat,gamma,G,I_o,omega,T,n,q,It,Bu,pick,pickf)

% preliminary elements for the algorithm
% compute the product G'*I_gamma*G (to speed up computations of deltabar)
GIG=G'*I_o*G;
% compute alphabar
alphabar=T+alpha0;

% initiate the Gibbs sampler
% initiate the counting of iterations
count=1;
pickcount=1;
% initiate the record matrices and cells
beta_gibbs=[];
F_gibbs=[];
L_gibbs=[];
phi_gibbs=[];
sigma_gibbs=[];
lambda_t_gibbs={};
sigma_t_gibbs={};

% step 1: determine initial values for the algorithm

% initial value for beta
beta=betahat;
% initial value for f_2,...,f_n
% obtain the triangular factorisation of sigmahat
[Fhat Lambdahat]=bear.triangf(sigmahat);
% obtain the inverse of Fhat
[invFhat]=bear.invltod(Fhat,n);
% create the cell storing the different vectors of invF
Finv=cell(n,1);
% store the vectors
for ii=2:n
Finv{ii,1}=invFhat(ii,1:ii-1);
end
% initial values for L_1,...,L_n
L=zeros(T,n);
% initial values for phi_1,...,phi_n
phi=ones(1,n);

% step 2: determine the sbar values and Lambda
sbar=diag(Lambdahat);
Lambda=sparse(diag(sbar));

% step 3: recover the series of initial values for lambda_1,...,lambda_T and sigma_1,...,sigma_T
lambda_t=repmat(diag(sbar),[1 1 T]);
sigma_t=repmat(sigmahat,[1 1 T]);

hbar = bear.parfor_progressbar(It,'Progress of the Gibbs sampler');  %create the progress bar

% run the Gibbs sampler
while count<=It
% count
   hbar.iterate(1);   % update progress by one iteration

% step 4: draw beta from its conditional posterior
% first compute the summations required for omegabar and betabar
summ1=zeros(q,q);
summ2=zeros(q,1);
   % run the summation
   for jj=1:T
   prodt=Xbart{jj,1}'/sigma_t(:,:,jj);
   summ1=summ1+prodt*Xbart{jj,1};
   summ2=summ2+prodt*yt(:,:,jj);
   end
% then obtain the inverse of omega0
invomega0=diag(1./diag(omega0));
% obtain the inverse of omegabar
invomegabar=summ1+invomega0;
% recover omegabar
C=chol(bear.nspd(invomegabar),'Lower')';
invC=C\speye(q);
omegabar=invC*invC';
% recover betabar
betabar=omegabar*(summ2+invomega0*beta0);
% finally, draw beta from its posterior
beta=betabar+chol(bear.nspd(omegabar),'lower')*randn(q,1);

% step 5: draw the series f_2,...,f_n from their conditional posteriors
   % recover first the residuals
   for jj=1:T
   epst(:,:,jj)=yt(:,:,jj)-Xbart{jj,1}*beta;
   end
   % then draw the vectors in turn
   for jj=2:n
   % first compute the summations required for upsilonbar and fbar
   summ1=zeros(jj-1,jj-1);
   summ2=zeros(jj-1,1);
      % run the summation
      for kk=1:T
      prodt=epst(1:jj-1,1,kk)*exp(-L(kk,jj));
      summ1=summ1+prodt*epst(1:jj-1,1,kk)';
      summ2=summ2+prodt*epst(jj,1,kk)';
      end
   summ1=(1/sbar(jj,1))*summ1;
   summ2=(-1/sbar(jj,1))*summ2;
   % then obtain the inverse of upsilon0
   invupsilon0=diag(1./diag(upsilon0{jj,1}));
   % obtain upsilonbar
   invupsilonbar=summ1+invupsilon0;
   C=chol(bear.nspd(invupsilonbar));
   invC=C\speye(jj-1);
   upsilonbar=full(invC*invC');
   % recover fbar
   fbar=upsilonbar*(summ2+invupsilon0*f0{jj,1});
   % finally draw f_i^(-1)
   Finv{jj,1}=fbar+chol(bear.nspd(upsilonbar),'lower')*randn(jj-1,1);
   end
% recover the inverse of F
invF=eye(n);
   for jj=2:n
   invF(jj,1:jj-1)=Finv{jj,1};
   end
% eventually recover F
F=bear.invltod(invF,n);
% then update sigma
sigma=F*Lambda*F';

% step 6: draw the series phi_1,...,phi_n from their conditional posteriors
% draw the parameters in turn
   for jj=1:n
   % estimate deltabar
   deltabar=L(:,jj)'*GIG*L(:,jj)+delta0;
   % draw the value phi_i
   phi(1,jj)=bear.igrandn(alphabar/2,deltabar/2);
   end

% step 7: draw the series lambda_i,t from their conditional posteriors, i=1,...,n and t=1,...,T
   % consider variables in turn
   for jj=1:n
      % consider periods in turn
      for kk=1:T
      % a candidate value will be drawn from N(lambdabar,phibar)
      % the definitions of lambdabar and phibar varies with the period, thus define them first
         % if the period is the first period
         if kk==1
         lambdabar=(gamma*L(2,jj))/(1/omega+gamma^2);
         phibar=phi(1,jj)/(1/omega+gamma^2);
         % if the period is the final period
         elseif kk==T
         lambdabar=gamma*L(T-1,jj);
         phibar=phi(1,jj);
         % if the period is any period in-between
         else
         lambdabar=(gamma/(1+gamma^2))*(L(kk-1,jj)+L(kk+1,jj));
         phibar=phi(1,jj)/(1+gamma^2);
         end
      % now draw the candidate
      cand=lambdabar+phibar^0.5*randn;
      % compute the acceptance probability
      prob=bear.mhprob2(jj,cand,L(kk,jj),sbar(jj,1),epst(:,1,kk),Finv{jj,1});
      % draw a uniform random number
      draw=rand;
         % keep the candidate if the draw value is lower than the prob
         if draw<=prob
         L(kk,jj)=cand;
         % if not, just keep the former value
         end
      end
   end
% then recover the series of matrices lambda_t and sigma_t
for jj=1:T
lambda_t(:,:,jj)=diag(sbar).*diag(exp(L(jj,:)));
sigma_t(:,:,jj)=F*lambda_t(:,:,jj)*F';
end

% record phase
   % if the burn-in sample phase is not yet over
   if count<=Bu
   % simply add 1 to the iteration count
   count=count+1;
   % on the other hand, if the burn-in sample phase is over
   elseif count>Bu
   % adding one iteration to the count will depend on wether post-burn selection applies
      % if there is no post burn selection
      if pick==0
      % record the results
      beta_gibbs(:,count-Bu)=beta;
      F_gibbs(:,:,count-Bu)=F;
      L_gibbs(:,:,count-Bu)=L;
      phi_gibbs(count-Bu,:)=phi;
      sigma_gibbs(:,count-Bu)=sigma(:);
         for jj=1:T
         lambda_t_gibbs{jj,1}(:,:,count-Bu)=lambda_t(:,:,jj);
         sigma_t_gibbs{jj,1}(:,:,count-Bu)=sigma_t(:,:,jj);
         end
      % then add one to the count
      count=count+1;
      % if there is post burn selection, only one draw over 'fpick' draws will be retained
      elseif pick==1
         % if the iteration does not correspond to fpick, don't record the results, don't increase the regular count, but do increase pickcount by 1, and do record the acceptance rate of the Metropolis-Hastings step
         if pickcount~=pickf
         pickcount=pickcount+1;
         % on the other hand, if the iteration does correspond to fpick
         elseif pickcount==pickf
         % do record the results
         beta_gibbs(:,count-Bu)=beta;
         F_gibbs(:,:,count-Bu)=F;
         L_gibbs(:,:,count-Bu)=L;
         phi_gibbs(count-Bu,:)=phi;
         sigma_gibbs(:,count-Bu)=sigma(:);
            for jj=1:T
            lambda_t_gibbs{jj,1}(:,:,count-Bu)=lambda_t(:,:,jj);
            sigma_t_gibbs{jj,1}(:,:,count-Bu)=sigma_t(:,:,jj);
            end
         % then increase the regular count by 1 and re-initialise pickcount
         count=count+1;
         pickcount=1;
         end
      end
   end

end
close(hbar);   %close progress bar

author: european-central-bank

repo: BEAR-toolbox


text:
function K = processConstraintSparse(K, bound)
% 这个函数对K进行处理，使其满足约束
% 用于稀疏矩阵方法
% 输入：
%      整体劲度矩阵 K
%      边界条件，格式参考README.MD
% 输出：
%      处理后的 K

for i = 1 : size(bound, 1)
    % 对所有的约束进行处理
    element_num = bound(i, 1);
    x_or_y      = bound(i, 2);
    j           = element_num * 2 - 2 + x_or_y;
    K(j, :) = 0;
    K(:, j) = 0;
    K(j, j) = 1;
end

author: Meelfy

repo: FEM

save_path: github-repos/MATLAB/Meelfy-FEM


text:
function rgb = brain

% returns a predefined color as [red green blue] values
%
% skin_surface         = [255 213 119]/255;
% outer_skull_surface  = [140  85  85]/255;
% inner_skull_surface  = [202 100 100]/255;
% cortex = [255 213 119]/255;
% black  = [0   0   0  ]/255;
% white  = [255 255 255]/255;
% red    = [255 0   0  ]/255;
% green  = [0   192 0  ]/255;
% blue   = [0   0   255]/255;
% yellow = [255 255 0  ]/255;
% cortex_light = [199 194 169]/255;
% cortex_dark  = [100 97 85]/255;

rgb = [202 100 100]/255;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [rh,rg,h,g]=rh2rg(rh)

%RH2RG    Calculates all the filters from the synthesis lowpass
%	  in the orthogonal case.
%
%	  [RH,RG,H,G]=RH2RG(RH) begins with the synthesis lowpass
%	  filter (RH) and returns the synthesis highpass filter (RG),
%	  the analysis lowpass filter (H) and the analysis highpass
%	  filter (G).
%
%	  It is an auxiliary function for orthogonal filters design.

%--------------------------------------------------------
% Copyright (C) 1994, 1995, 1996, by Universidad de Vigo
%
%
% Uvi_Wave 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, or (at your option) any
% later version.
%
% Uvi_Wave 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 Uvi_Wave; see the file COPYING.  If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%
%       Author: Nuria Gonzalez Prelcic
%       e-mail: Uvi_Wave@tsc.uvigo.es
%--------------------------------------------------------

% Calculate rg from rh.

for i=1:length(rh)
	rg(i) = -(-1)^i*rh(length(rh)-i+1);
end

% Calculate h and g

h=rh(length(rh):-1:1);
g=rg(length(rg):-1:1);

author: alexandrebarachant


text:
function timer_etime_test01 ( )

%*****************************************************************************80
%
%% TEST01 times the RAND routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 December 2005
%
%  Author:
%
%    John Burkardt
%
  n_log_min = 10;
  n_log_max = 20;
  n_min = 2^n_log_min;
  n_max = 2^n_log_max;
  n_rep = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  Time the MATLAB RAND routine:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    x = rand(n,1);\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data vectors will be of minimum size %d\n', n_min );
  fprintf ( 1, '  Data vectors will be of maximum size %d\n', n_max );
  fprintf ( 1, '  Number of repetitions of the operation: %d\n', n_rep );

  for i_rep = 1 : n_rep

    for n_log = n_log_min : n_log_max

      n = 2^n_log;

      time1 = clock;

      x = rand(n,1);

      time2 = clock;

      delta(n_log,i_rep) = etime ( time2, time1 );

    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01 Results:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Vector Size  Rep #1        Rep #2        Rep #3        ' );
  fprintf ( 1, 'Rep #4        Rep #5\n' );
  fprintf ( 1, '\n' );

  for n_log = n_log_min : n_log_max
    n = 2^n_log;
    fprintf ( 1, '%10d', n );
    for j = 1 : n_rep
      fprintf ( '%14f', delta(n_log,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function importGen3Model
% Imports Kinova Gen3 Ultra lightweight robot model into a Rigid Body Tree
%
% Copyright 2019 The MathWorks, Inc.

    % Load the URDF file
    %
    % NOTE: This requires you to have the kortex_description folder on
    % your path, which you can download from
    % https://github.com/Kinovarobotics/ros_kortex.git
    %
    % Then, you have to convert the gen3.xacro file to a URDF file using
    % the following commands in a ROS enabled terminal:
    %  $ cd PATH/TO/ros_kortex/kortex_description/robots
    %  $ rosrun xacro xacro --inorder -o gen3.urdf gen3.xacro
    addpath(genpath('kortex_description'))
    gen3 = importrobot('gen3.urdf');

    % Add a "dummy" gripper link
    gripperLength = 0.1; % Gripper length in meters
    gripperBody = rigidBody('Gripper');
    gripperJoint = rigidBodyJoint('GripperLink','fixed');
    T = rotm2tform([0 1 0;0 0 1;1 0 0]) * trvec2tform([gripperLength 0 0]);
    setFixedTransform(gripperJoint,T); % Move and orient the gripper
    gripperBody.Joint = gripperJoint;
    addBody(gen3,gripperBody,'end_effector_link');
    % Add a "dummy" mesh
    addVisual(gen3.Bodies{9},'Mesh','cylinder.stl', ...
              trvec2tform([-0.1 0 0]) * axang2tform([0 1 0 pi/2]));

    % Configure the data format and show the robot in the home position
    gen3.DataFormat = 'row';
    load gen3positions
    show(gen3,jointAnglesHome');

    % Save the Rigid Body Tree to a file
    curDir = pwd;
    saveDir = fileparts(mfilename('fullpath'));
    cd(saveDir)
    save gen3 gen3
    cd(curDir)

end

author: mathworks-robotics

repo: trajectory-planning-robot-manipulators


text:
function roisArray=dtiCreateRoisFromNifti(niftiFile)

%Create an array of ROIs corresponding to unique values in niftiFile.
%niftiFile should contain integer values labeling the voxels from the list
%of ROIs.

%ER 11/2009 wrote it

roiNii = niftiRead(niftiFile);

        % Find all unique, non-zero indices, each distinct index will be a
        % mask for a distinct ROI
        roiNii.data(isnan(roiNii.data))=0;
        roisToMake=unique(uint8(roiNii.data));
        roisToMake=roisToMake(roisToMake~=0);

        % For each distinct index, create a distinct ROI
        for jj=1:length(roisToMake)

            thisRoi = find(roiNii.data==roisToMake(jj));
            [x1,y1,z1] = ind2sub(size(roiNii.data), thisRoi);
            roisArray(jj) = dtiNewRoi([roiNii.fname '_' num2str(roisToMake(jj))], rand(1, 3));
            roisArray(jj).coords = mrAnatXformCoords(roiNii.qto_xyz, [x1,y1,z1]);
            roisArray(jj).name=[roiNii.fname '_' num2str(jj)];

        end
        if ~exist('roisArray', 'var')
            roisArray=[];
        end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [VEC,C,V,E] = prtUtilPca(DATA,P)
%[VEC,C,V,E] = getPcaVec(DATA,P);
% xxx Need Help xxx
%

%Peter Torrione
C = cov(DATA);
[V,E] = eig(C);

%sort
[E,IND] = sort(abs(diag(E)));
%re-sorta
E = flipud(E(:));
IND = flipud(IND(:));
V = V(:,IND);

if P < 1
    NRG = cumsum(E)./sum(E);
    NRGind = find(NRG > P,1,'first');
    VEC = V(:,1:NRGind);
else
    VEC = V(:,1:P);
end

PLOTTING = 0;
if PLOTTING
    plot(DATA'./max(abs(DATA(:))));
    hold on;
    h = plot(VEC./max(abs(VEC(:))));
    set(h,'linewidth',3);
    pause
    close all;
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function t = getLinesearch(problem, x, d, storedb, key)
% Returns a hint for line-search algorithms.
%
% function t = getLinesearch(problem, x, d)
% function t = getLinesearch(problem, x, d, storedb)
% function t = getLinesearch(problem, x, d, storedb, key)
%
% For a line-search problem at x along the tangent direction d, computes
% and returns t such that retracting t*d at x yields a good point around
% where to look for a line-search solution. That is: t is a hint as to
% "how far to look" along the line.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% See also: canGetLinesearch

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, July 17, 2014.
% Contributors:
% Change log:
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    if isfield(problem, 'linesearch')
    %% Compute the line-search hint function using linesearch.

        % Check whether this function wants to deal with storedb or not.
        switch nargin(problem.linesearch)
            case 2
                t = problem.linesearch(x, d);
            case 3
                % Obtain, pass along, and save the store for x.
                store = storedb.getWithShared(key);
                [t, store] = problem.linesearch(x, d, store);
                storedb.setWithShared(store, key);
            case 4
                % Pass along the whole storedb (by reference), with key.
                t = problem.linesearch(x, d, storedb, key);
            otherwise
                up = MException('manopt:getLinesearch:badfun', ...
                    'linesearch should accept 2, 3 or 4 inputs.');
                throw(up);
        end

    else
    %% Abandon computing the line-search function.

        up = MException('manopt:getLinesearch:fail', ...
            ['The problem description is not explicit enough to ' ...
             'compute a line-search hint.']);
        throw(up);

    end

end

author: NicolasBoumal

repo: manopt


text:
function [varargout] = dcblock(varargin);
%
% DCBLOCK   determines the filter coefficient a for the dc blocking/high-
%           pass filter y(t) = x(t) - x(t-1) + a*y(t-1).
%
%           [a] = dcblock(Fc) returns the filter coefficient a for
%           normalized cut-on frequency Fc.
%
%           [a] = dcblock(fc,fs)returns the filter coefficient a for cut-on
%           frequency fc, and sampling frequency fs, when both are given in
%           Hz.
%
%           [aQ] = dcblock(fc,fs,B) returns the quantized filter
%           coefficient aQ when the quantization level B is given.
%
%           dcblock(...) plots the frequency and phase response of the
%           filter.
%
%           Note that if two parameters are passed to the function without
%           an output argument, the plot routine always interprets the
%           first as fc and the second as fs.
%
% SYNTAX    [a] = dcblock(Fc);
%           [a] = dcblock(fc,fs);
%           [aQ] = dcblock(fc,fs,B);
%           [Fc,fc] = dcblock(a,fs);
%           dcblock(Fc)
%           dcblock(fc,fs)
%           dcblock(fc,fs,B)
%
% INPUT     Fc      Normalized frequency cut-on of the filter.
%                   Fc = fc[Hz]/(fs[Hz]/2); fs/2 is the Nyquist frequency.
%                   Note that 0 < Fc < 1/3;
%
%           fc      filter cut-on frequency in Hz (or other reciprocal
%                   units of frequency e.g. 1/metre).
%
%           fs      sampling frequency in Hz (or other reciprocal
%                   units of frequency e.g. 1/metre).
%
%           B       quantization level for fixed point processing e.g. B is
%                   16 for 16-bit operation.
%
% OUTPUT    a       filter coefficient; 0 < a < 1 for floating point
%                   operation.
%
%           aQ      quantized filter coefficient a; 0 < aQ < 2^(B-1) - 1.
%                   The quantization level is set by B which is assumed to
%                   be signed integer.
%
% EXAMPLES  (1) Calculate the filter coefficient a when the normalized
%               frequency is 1E-5:
%
%           >> format long                      % set display format
%           >> a = dcblock(1E-5)                % Fc = 0.00001
%           >> a = 0.99994558749947             % filter coefficient
%
%           (2) Calculate the normalized and true cut-on frequency when the
%               coefficient a is 0.9987 and the sampling frequency is 1MHz:
%
%           >> [Fc,fc] = dcblock(0.9987,1E6)    % a = 0.9987; fs = 1MHz
%           >> Fc = 2.3906E-004                 % this output is normalized
%           >> fc = 119.5323                    % this output is in Hz
%
%           Note the difference in inputs between [Fc,fc] = dcblock(a,fs),
%           the function call [a] = dcblock(fc,fs), and its plot facility
%           dcblock(fc,fs). You can plot the filter response using
%
%           >> dcblock(119.5323,1E6)            % fc = 119.5323; fs = 1E6
%
%           (3) Carry out filtering operation using the dc blocking filter
%               with cut-on frequency of 70Hz and sampling rate of 10kHz.
%
%           t = -0.1:1/10000:0.1-1/10000;       % set up time base
%           s = sin(2*pi*2.5*t+pi/6)+0.01*randn(size(t))+1.25;
%           x = gauspuls(t,250,0.6) + s;        % noisy Gaussian test pulse
%           p = dcblock(70,10000);              % get filter coefficient
%           b = [1 -1];                         % set up differentiator
%           a = [1 -p];                         % set up integrator
%           y = filter(b,a,x);                  % filter test data
%           plot(t,x,t,y)                       % plot results to screen
%           legend('Original','DC Filtered','location','northwest');
%
%           The forward pass of the filter creates a distortion in the
%           Gaussian pulse. This distortion can be removed by using the
%           filtfilt form of the filter. If you now try
%
%           y2 = filtfilt(b,a,x);               % filter test data
%           plot(t,x,t,y,t,y2)                  % plot results to screen
%           legend('Original','DC Filtered','Undistorted','location','northwest');
%
%           you will notice that the symmetry of the Gaussian pulse has
%           been restored.
%
% See also FILTER, FILTFILT, SPTOOL, FDATOOL.

% Author:   J M De Freitas
%           QinetiQ Ltd, Winfrith Technology Centre
%           Winfrith, Dorchester DT2 8XJ, UK.
%
% Email     jdefreitas@qinetiq.com
%
% Revision  1.0
%
% Date      11 October 2006.
%

% Perform checks on input and output
Lin = nargin;
Lout = nargout;
if Lin > 3
    error('dcblock error: Too many input arguments.');
end
if Lin < 1
    errror('dcblock error: No input argument used.');
end
if Lout > 2
    error('dcblock error: Too many output arguments.');
end

% Calculate a when Fc is given and plot response if required
if Lin == 1
    Fc = varargin{1};
    if (Fc < 0)|(Fc > 1/3)
        error('dcblock error: Normalized cut-on frequency must be between 0 and 1/3.');
    end
    a = coef(Fc);
    if Lout == 1
        varargout{1} = a;
    elseif Lout == 0
        PlotResp(1,a);
    end
end

% Calculate a when fc and fs are given and plot response if required
% or Fc and fc when a and fs are given
if Lin == 2
    if (Lout == 1)|(Lout == 0)
        fc = varargin{1};
        fs = varargin{2};
        Fc = 2*fc/fs;
        if (Fc < 0)|(Fc > 1/3)
            error('dcblock error: Cut-on frequency must be between 0 and fs/6.');
        end
        a = coef(Fc);
        if Lout == 1
            varargout{1} = a;
        else
            PlotResp(fs/2,a);
        end
    elseif Lout == 2
        a = varargin{1};
        fs = varargin{2};
        if (a < 0)|(a > 1)
            error('dcblock error: Filter coefficient a must be between 0 and 1.');
        end
        Fc = atan(sqrt(3)*(1-a^2)/(1+4*a+a^2))/pi;
        varargout{1} = Fc;
        varargout{2} = Fc*fs/2;
    end
end

% Calculate a when B, fc and fs are given and plot if required
if Lin == 3
    fc = varargin{1};
    fs = varargin{2};
    B = floor(varargin{3});
    Fc = 2*fc/fs;
    if (Fc < 0)|(Fc > 1/3)
        error('dcblock error: Cut-on frequency fc must be between 0 and fs/6.');
        return
    end
    if (B <= 2)|(B > 300)
        error('dcblock error: Quantization level B must be between 2 and 300.');
        return
    end
    a = coef(Fc);
    if Lout == 1
        varargout{1} = floor(a*(2^(B-1)-1));
    elseif Lout == 0
        aQ = floor(a*(2^(B-1)-1))/(2^(B-1)-1);
        PlotResp(fs/2,aQ);
    end
end

function a = coef(Fc);
% get up coefficient
    a = (sqrt(3) - 2*sin(pi*Fc))/(sin(pi*Fc) + sqrt(3)*cos(pi*Fc));

function PlotResp(Fn,a);
% Plot frequency and phase response
    f = 0.000001:0.000001:1-0.000001;
    h = (1+a)*sin(pi*f/2)./sqrt(1+a^2-2*a*cos(pi*f));
    phi = atan((1-a)*(sin(pi*f)./(1-cos(pi*f)))/(1+a));
    indigo = [.48 .06 .89];
    subplot(2,1,1);
        loglog(Fn*f,h,'LineWidth',1.5,'Color',indigo);
        title('\bf Frequency Response: IIR DC Blocking Filter','Color','k');
        if Fn == 1
            xlabel('Normalized Frequency');
        else
            xlabel('Frequency (Hz)');
        end
        ylabel('Magnitude ');
        minY = min(get(gca,'YTick'));
        minX = min(get(gca,'XTick'));
        maxX = max(get(gca,'XTick'));
        if minY >= 0.1
            TickY = [0.1 1];
        else
            for j = 1:abs(log10(minY))+1 TickY(j) = 10^(log10(minY)+j-1); end
        end
        for j = 1:(log10(maxX)-log10(minX))+1 TickX(j) = 10^(log10(minX)+j-1); end
        ylim([TickY(1) 1.1]);
        set(gca,'XGrid','on','YGrid','on','XMinorGrid','on',...
        'YMinorGrid','on','GridLineStyle','-',...
        'MinorGridLineStyle','-','Xtick',TickX,'Ytick',TickY,...
        'Color',[0.835 0.918 0.918]);
    subplot(2,1,2);
        loglog(Fn*f,phi*180/pi,'LineWidth',1.5,'Color',indigo);
        title('\bf Phase Response: IIR DC Blocking Filter','Color','k');
        if Fn == 1
            xlabel('Normalized Frequency');
        else
            xlabel('Frequency (Hz)');
        end
        ylabel('Phase (degree) ');
        ylim([0.1 100]);
        minX = min(get(gca,'XTick'));
        maxX = max(get(gca,'XTick'));
        for j = 1:(log10(maxX)-log10(minX))+1
            TickX2(j) = minX*10^(j-1) ;
        end
        set(gca,'XGrid','on','YGrid','on','XMinorGrid','on',...
        'YMinorGrid','on','GridLineStyle','-',...
        'MinorGridLineStyle','-','Xtick',TickX2,'YTick',[0.1 1 10 100],...
        'Color',[0.835 0.918 0.918]);


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

%*****************************************************************************80
%
%% LOBATTO_COMPUTE computes a Lobatto quadrature rule.
%
%  Discussion:
%
%    The integral:
%
%      Integral ( -1 <= X <= 1 ) F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( X(I) )
%
%    The quadrature rule will integrate exactly all polynomials up to
%    X^(2*N-3).
%
%    The Lobatto rule is distinguished by the fact that both endpoints
%    (-1 and 1) are always abscissas of the rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2015
%
%  Author:
%
%    Original MATLAB code by Greg von Winckel.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Claudio Canuto, Yousuff Hussaini, Alfio Quarteroni, Thomas Zang,
%    Spectral Methods in Fluid Dynamics,
%    Springer, 1993,
%    ISNB13: 978-3540522058,
%    LC: QA377.S676.
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%    Daniel Zwillinger, editor,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996,
%    ISBN: 0-8493-2479-3.
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights,
%
 if ( n == 1 )
    x(0) = -1.0;
    w(0) = 2.0;
    return
  end

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

  tolerance = 100.0 * eps;
%
%  Initial estimate for the abscissas is the Chebyshev-Gauss-Lobatto nodes.
%
  x(1:n,1) = cos ( pi * ( 0 : n - 1 ) / ( n - 1 ) )';
  xold(1:n,1) = 2.0;

  while ( tolerance < max ( abs ( x(1:n,1) - xold(1:n,1) ) ) )

    xold(1:n,1) = x(1:n,1);

    p(1:n,1) = 1.0;
    p(1:n,2) = x(1:n,1);

    for j = 2 : n-1
      p(1:n,j+1) = ( ( 2 * j - 1 ) * x(1:n,1) .* p(1:n,j)     ...
                   + (   - j + 1 ) *             p(1:n,j-1) ) ...
                   / (     j     );
    end

    x(1:n,1) = xold(1:n,1) - ( x(1:n,1) .* p(1:n,n) - p(1:n,n-1) ) ...
             ./ ( n * p(1:n,n) );
  end

  x(1:n) = r8vec_reverse ( n, x );

  w(1:n,1) = 2.0 ./ ( ( n - 1 ) * n * p(1:n,n).^2 );

  return
end


text:
function [opts] = FDDL_SpaCoef(ipts,par)
% ========================================================================
% Coefficient updating of FDDL, Version 1.0
% Copyright(c) 2011  Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang
% All Rights Reserved.
%
% -----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is here
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
%
% This is an implementation of the algorithm for updating the
% Coefficient matrix of FDDL (fix the dictionary)
%
% Please refer to the following paper
%
% Meng Yang, Lei Zhang, Xiangchu Feng, and David Zhang,"Fisher Discrimination
% Dictionary Learning for Sparse Representation", In IEEE Int. Conf. on
% Computer Vision, 2011.
% L. Rosasco, A. Verri, M. Santoro, S. Mosci, and S. Villa. Iterative
% Projection Methods for Structured Sparsity Regularization. MIT Technical
% Reports, MIT-CSAIL-TR-2009-050,CBCL-282, 2009.
% J. Bioucas-Dias, M. Figueiredo, ?A new TwIST: two-step iterative shrinkage
% /thresholding  algorithms for image restoration?, IEEE Transactions on
% Image Processing, December 2007.
%----------------------------------------------------------------------
%
%  Inputs :   (1) ipts :    the structre of input data
%                    .D     the dictionary
%                    .X     the training data
%                    .A     the coefficient matrix in the last iteration
%                    .trls  the labels of training data
%             (2) par :     the struture of input parameters
%                    .tau   the parameter of sparse constraint of coef
%                    .lambda  the parameter of within-class scatter
%                    .dls     the labels of dictionary's columns
%                    .index   the label of the class being processed
%
% Outputs:    (1) opts :    the structure of output data
%                    .A     the coefficient matrix
%                    .ert   the total energy sequence
%
%---------------------------------------------------------------------

par.nIter    =     200;   % maximal iteration number
par.isshow   =     false;
par.citeT    =     1e-6;  % stop criterion
par.cT       =     1e+10; % stop criterion

m            =    size(ipts.D,2);
fish_tau3    =    1;  % parameter of ||A_i-D_iX_i^i||_F^2 in fidelity term
fish_tau2    =    1;  % parameter of ||D_jX_i^j||_F^2 in fidelity term Eq.(4)
drls         =    par.dls;
D            =    ipts.D;
X            =    ipts.X;
A            =    ipts.A;
tau          =    par.tau;
lambda1      =    par.tau;
lambda2      =    par.lambda2;
lambda3      =    par.lambda2;
lambda4      =    par.lambda2;
trls         =    ipts.trls;
classn       =    length(unique(trls));
nIter        =    par.nIter;
c            =    par.c;
sigma        =    c;
tau1         =    tau/2;
index        =    par.index;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TWIST parameter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for_ever           =         1;
IST_iters          =         0;
TwIST_iters        =         0;
sparse             =         1;
verbose            =         1;
enforceMonotone    =         1;
lam1               =         1e-4;   %default minimal eigenvalues
lamN               =         1;      %default maximal eigenvalues
rho0               =         (1-lam1/lamN)/(1+lam1/lamN);
alpha              =         2/(1+sqrt(1-rho0^2));        %default,user can set
beta               =         alpha*2/(lam1+lamN);         %default,user can set

%%%%%%%%%%%%%%%%%%%%%%%%
%preprocessing
%%%%%%%%%%%%%%%%%%%%%%%%
Ai                 =          X;  % the i-th training data
Xa                 =          A;  %
Xi                 =          A(:,trls==index);
Xt_now             =          A(:,trls==index);
% Xi                 =          zeros(size(A(:,trls==index)));
% Xt_now             =          zeros(size(A(:,trls==index)));

newpar.n_d          =   size(Ai,2);             % the sample number of i-th training data
n                   =   size(Xa,2);             % the total sample number of training data

for ci = 1:classn
    t_n_d = sum(trls==ci);
    t_b_line_i     =   ones(t_n_d,newpar.n_d)./t_n_d;
    t_c_j = ones(t_n_d,newpar.n_d)./n;
    t_b_angle_i           =   t_b_line_i-t_c_j;
    CJ(ci).M  = t_c_j;
    BAI(ci).M = t_b_angle_i;
end

newpar.B_line_i     =   ones(newpar.n_d,newpar.n_d)./newpar.n_d;
newpar.C_j          =   ones(newpar.n_d,newpar.n_d)./n;
newpar.CjCj         =   (newpar.C_j)*(newpar.C_j)';
newpar.C_line       =   ones(n,newpar.n_d)./n;
B_i                 =   eye(newpar.n_d,newpar.n_d)-newpar.B_line_i;
newpar.BiBi         =   B_i*(B_i)';
B_angle_i           =   newpar.B_line_i-newpar.C_j;
newpar.Bai          =   B_angle_i;
newpar.BaiBai       =   B_angle_i*(B_angle_i)';
Xo                  =   Xa;
Xo(:,trls==index)   =   0;
G_X_i               =   Xo*newpar.C_line;
newpar.BaiGxi       =   B_angle_i*(G_X_i)';
newpar.DD           =   D'*D;
newpar.DAi          =   D'*Ai;
Di0                 =   D;
Di0(:,drls~=index)  =   0;
newpar.Di0Di0       =   (Di0)'*Di0;
newpar.Di0Ai        =   (Di0)'*Ai;

newpar.DoiDoi       =   zeros(size(D,2));
for t_i  =  1:classn
    if t_i ~= index
    Doi                 =   D;
    Doi(:,drls~=t_i)      =   0;
    newpar.DoiDoi       =   newpar.DoiDoi+(Doi)'*Doi;
    end
end

newpar.m            =   m;                 % the number of dictionary column atoms

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Xa(:,trls==index)  =  Xi;
xm2       =      Xi;%A(:,trls==index);
xm1       =      Xi;%A(:,trls==index); % now

[gap] = FDDL_Class_Energy(Ai,D,xm1,Xa,drls,trls,index,...
        lambda1,lambda2,lambda3,lambda4,classn,fish_tau2,fish_tau3);
prev_f   =   gap;
ert(1) = gap;
for n_it = 2 : nIter;

   Xa(:,trls==index)  =  Xi;

   while for_ever
        % IPM estimate

        grad = FDDL_Gradient_Comp(xm1,Xa,classn,index,...
        lambda2,lambda3,lambda4,fish_tau2,fish_tau3,trls,drls,newpar,...
        BAI,CJ);

        v        =   xm1(:)-grad./(2*sigma);
        tem      =   soft(v,tau1/sigma);
        x_temp   =   reshape(tem,[size(D,2),size(xm1,2)]);

        if (IST_iters >= 2) | ( TwIST_iters ~= 0)
            % set to zero the past when the present is zero
            % suitable for sparse inducing priors
            if sparse
                mask    =   (x_temp ~= 0);
                xm1     =   xm1.* mask;
                xm2     =   xm2.* mask;
            end
            % two-step iteration
            xm2    =   (alpha-beta)*xm1 + (1-alpha)*xm2 + beta*x_temp;
            % compute residual
           [gap] = FDDL_Class_Energy(Ai,D,xm2,Xa,drls,trls,index,...
                lambda1,lambda2,lambda3,lambda4,classn,fish_tau2,fish_tau3);

           f   =   gap;

            if (f > prev_f) & (enforceMonotone)
                TwIST_iters   =  0;  % do a IST iteration if monotonocity fails
            else
                TwIST_iters =   TwIST_iters+1; % TwIST iterations
                IST_iters   =    0;
                x_temp      =   xm2;
                if mod(TwIST_iters,10000) ==0
                   c = 0.9*c;
                   sigma = c;
                end
                break;  % break loop while
            end
        else

        [gap] = FDDL_Class_Energy(Ai,D,x_temp,Xa,drls,trls,index,...
        lambda1,lambda2,lambda3,lambda4,classn,fish_tau2,fish_tau3);

        f   =   gap;

         if f > prev_f
                % if monotonicity  fails here  is  because
                % max eig (A'A) > 1. Thus, we increase our guess
                % of max_svs
                c         =    2*c;
                sigma     =    c;
                if verbose
%                     fprintf('Incrementing c=%2.2e\n',c);
                end
                if  c > par.cT
                    break;  % break loop while
                end
                IST_iters = 0;
                TwIST_iters = 0;
           else
                TwIST_iters = TwIST_iters + 1;
                break;  % break loop while
           end
        end

    end

    citerion      =   abs(f-prev_f)/abs(prev_f);
    if citerion < par.citeT | c > par.cT
%        fprintf('Stop!\n c=%2.2e\n citerion=%2.2e\n',c,citerion);
       break;
    end

    xm2           =   xm1;
    xm1           =   x_temp;
    Xt_now        =   x_temp;
    Xi            =   Xt_now;
    prev_f        =   f;
    ert(n_it)     =   f;
%     fprintf('Iteration:%f  Total gap:%f\n',n_it,ert(n_it-1));
   end

opts.A     =       Xt_now;
opts.ert   =       ert;

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function camCreateHARDIFitScript(scr_modelfit_filename,scr_after_modelfit_filename,scheme_filename,raw_filename,num_procs,bgthresh)
%
%
%
%
%

% Parameters

% Number parallel processes
if ieNotDefined('num_procs')
    num_procs=1;
end
% Reduced Encoding for spherical deconvolution
re = 16;
% Threshold on b0 for processing background voxels
if ieNotDefined('bgthresh')
    bgthresh = 200;
end
% Volumetric information
raw = niftiRead(raw_filename);
vol_dim = raw.dim;
xd = vol_dim(1);
yd = vol_dim(2);
zd = vol_dim(3);
gd = vol_dim(4);

fid = fopen(scr_modelfit_filename,'wt');
fprintf(fid,'#!/usr/bin/env bash\n\n');

fprintf(fid,'cam_run_mp_hardi_slice.sh %d %d %d %d %d %d %d %s \n', num_procs, re, xd, yd, zd, gd, bgthresh, scheme_filename);

fclose(fid);

fid = fopen(scr_after_modelfit_filename,'wt');
fprintf(fid,'#!/usr/bin/env bash\n\n');

fprintf(fid,'cam_run_after_mp_hardi_slice.sh %d %d %d %d %d %d %d %s \n', num_procs, re, xd, yd, zd, gd, bgthresh, scheme_filename);

fclose(fid);

return


text:
function [newVectors, whiteningMatrix, dewhiteningMatrix] = whitenv ...
    (vectors, E, D, s_verbose);
%WHITENV - Whitenv vectors.
%
% [newVectors, whiteningMatrix, dewhiteningMatrix] = ...
%                               whitenv(vectors, E, D, verbose);
%
% Whitens the data (row vectors) and reduces dimension. Returns
% the whitened vectors (row vectors), whitening and dewhitening matrices.
%
% ARGUMENTS
%
% vectors       Data in row vectors.
% E             Eigenvector matrix from function 'pcamat'
% D             Diagonal eigenvalue matrix from function 'pcamat'
% verbose       Optional. Default is 'on'
%
% EXAMPLE
%       [E, D] = pcamat(vectors);
%       [nv, wm, dwm] = whitenv(vectors, E, D);
%
%
% This function is needed by FASTICA and FASTICAG
%
%   See also PCAMAT

% @(#)$Id: whitenv.m,v 1.3 2003/10/12 09:04:43 jarmo Exp $

% ========================================================
% Default value for 'verbose'
if nargin < 4, s_verbose = 'on'; end

% Check the optional parameter verbose;
switch lower(s_verbose)
 case 'on'
  b_verbose = 1;
 case 'off'
  b_verbose = 0;
 otherwise
  error(sprintf('Illegal value [ %s ] for parameter: ''verbose''\n', s_verbose));
end

% ========================================================
% In some cases, rounding errors in Matlab cause negative
% eigenvalues (elements in the diagonal of D). Since it
% is difficult to know when this happens, it is difficult
% to correct it automatically. Therefore an error is
% signalled and the correction is left to the user.
if any (diag (D) < 0),
  error (sprintf (['[ %d ] negative eigenvalues computed from the' ...
		   ' covariance matrix.\nThese are due to rounding' ...
		   ' errors in Matlab (the correct eigenvalues are\n' ...
		   'probably very small).\nTo correct the situation,' ...
		   ' please reduce the number of dimensions in the' ...
		   ' data\nby using the ''lastEig'' argument in' ...
		   ' function FASTICA, or ''Reduce dim.'' button\nin' ...
		   ' the graphical user interface.'], ...
		  sum (diag (D) < 0)));
end

% ========================================================
% Calculate the whitening and dewhitening matrices (these handle
% dimensionality simultaneously).
whiteningMatrix = inv (sqrt (D)) * E';
dewhiteningMatrix = E * sqrt (D);

% Project to the eigenvectors of the covariance matrix.
% Whiten the samples and reduce dimension simultaneously.
if b_verbose, fprintf ('Whitening...\n'); end
newVectors =  whiteningMatrix * vectors;

% ========================================================
% Just some security...
if ~isreal(newVectors)
  error ('Whitened vectors have imaginary values.');
end

% Print some information to user
if b_verbose
  fprintf ('Check: covariance differs from identity by [ %g ].\n', ...
    max (max (abs (cov (newVectors', 1) - eye (size (newVectors, 1))))));
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function doTheCalcThingMonte(obj)
		%everything for the calculation itself goes in here
		import mapseis.projector.*;
		import mapseis.calc.declustering.*;

		%Calculation with an added MonteCarlo simulation
		%(the things done by Thomas Van Stiphout)
		%LATER, first to the normal version

		%What the monte-carlo reasenberg does is simple

		%get parameter into variables
		ShortCatalog = getShortCatalog(obj.Datastore,obj.SendedEvents);
		taumin=obj.CalcParameter.Tau_min;
		taumax=obj.CalcParameter.Tau_max;
		xk=obj.CalcParameter.XMagFactor;
		xmeff=obj.CalcParameter.XMagEff;
		P=obj.CalcParameter.ProbObs;
		rfact=obj.CalcParameter.RadiusFactor;
		err=obj.CalcParameter.EpiError;
		derr=obj.CalcParameter.DepthError;

		%run declustering
		disp('please wait....');

		tic;
		[clusterID,EventType,AlgoInfo] = ReasenbergDecluster(taumin,taumax,xk,xmeff,P,rfact,err,derr,ShortCat);
		obj.CalcTime=toc;
		disp('finished!');

		%store result
		obj.CalcRes.clusterID=clusterID;
		obj.EventType=EventType;
		obj.AlgoInfo=AlgoInfo;

		%write to datastore
		obj.Datastore.setDeclusterData(EventType,clusterID,obj.SendedEvents,[]);
		obj.Datastore.DeclusterMetaData=AlgoInfo;

		%correct for old datastore version, it makes sense to do this here, because it is needed only by the declustering
		NumberedUserData=getNumberedFields(obj.Datastore);
		setNumberedFields(obj.Datastore,union(NumberedUserData,{'Month','Day','Hour','Minute','Second','MilliSecond' ,'DecYear'}));

end

author: CelsoReyes


text:
function linplus_test035 ( )

%*****************************************************************************80
%
%% TEST035 tests R83_GS_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 100;
  maxit = 1000;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST035\n' );
  fprintf ( 1, '  For a real tridiagonal system,\n' );
  fprintf ( 1, '  R83_GS_SL solves a linear system using\n' );
  fprintf ( 1, '    Gauss-Seidel iteration\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
  fprintf ( 1, '  Iterations per call = %d\n', maxit );
  fprintf ( 1, '\n' );
%
%  Set the matrix values.
%
  a(1,1)     =  0.0E+00;
  a(1,2:n)   = -1.0E+00;
  a(2,1:n)   =  2.0E+00;
  a(3,1:n-1) = -1.0E+00;
  a(3,n)     =  0.0E+00;

  for job = 0 : 1

    if ( job == 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Solving A * x = b.\n' );
      fprintf ( 1, '\n' );
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Solving A'' * x = b.\n' );
      fprintf ( 1, '\n' );
    end
%
%  Set the desired solution.
%
    x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
    if ( job == 0 )
      b = r83_mxv ( n, a, x );
    else
      b = r83_vxm ( n, a, x );
    end
%
%  Set the starting solution.
%
    x(1:n) = 0.0E+00;
%
%  Solve the linear system.
%
    for i = 1 : 3

      x_new = r83_gs_sl ( n, a, b, x, maxit, job );

      r8vec_print_some ( n, x_new, 1, 10, '  Current solution estimate:' );

      x(1:n) = x_new(1:n);

    end

  end

  return
end


text:
classdef blockCoeff
%BLOCKCOEFF   Class to convert linear operator to derivative coefficents.
%   This class is not intended to be called directly by the end user.
%
% See also LINOP, CHEBOP, CHEBOPPREF.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Developer notes
%
% This class converts a linBlock object into a list of coefficients for the
% (descending) powers of the derivative in the operator.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Access = public )
        coeffs = [];
        domain
        pref = chebfunpref;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function A = blockCoeff(varargin)
            % When called with no arguments, the returned object causes the
            % block's stack to be evaluated with these methods to produce
            % coefficients.
            if ( isempty(varargin{1}) )
                pref = cheboppref;
                A.domain = pref.domain;
                return

            % Calling the constructor with a linBlock argument initiates the
            % process of evaluating the stack with a dummy object of this class.
            elseif ( isa(varargin{1}, 'linBlock') )
                L = varargin{1};
                dummy = blockCoeff([]);
                dummy.domain = L.domain;
                if ( nargin == 2 )
                    dummy.pref = varargin{2};
                end
                A = L.stack(dummy);

            % If the constructor is called with data, just make a regular object
            % out of it.
            else
                f = varargin{1};
                if ( ~iscell(f) )
                    f = {f};
                end
                A.coeffs = f;
                A.domain = varargin{2};
                if ( nargin == 3 )
                    A.pref = varargin{3};
                end
            end
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% OTHER BASIC CONSTRUCTORS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        % These are the basic constructors.

        function I = eye(A)
            I = blockCoeff(chebfun(1, A.domain, A.pref), A.domain);
        end

        function I = zeros(A)
            I = blockCoeff(chebfun(0, A.domain, A.pref), A.domain);
        end

        function F = mult(A, f)
            F = blockCoeff(f, A.domain);
        end

        function C = cumsum(A, m)
            error('CHEBFUN:BLOCKCOEFF:cumsum:notSupported', ...
                'Conversion of integration to coefficients is not supported.')
        end

        function D = diff(A, order)

            if ( nargin < 2 )
                order = 1;
            end

            if ( ~isempty(A.coeffs) )
                % This syntax is used by the mtimes method to differentiate
                % an existing operator.
                D = A;
            else
                % Initialize with an identity.
                D = eye(A);
            end

            % Differentiate the requested number of times.
            c = D.coeffs;
            for d = 1:order
                m = numel(c);
                c = c([1, 1:m]);
                for k = 2:m
                    c{k} = diff(c{k}) + c{k+1};
                end
                c{m+1} = diff(c{m+1});
            end
            D = blockCoeff(c, A.domain);

        end

        function S = sum(A)
            error('CHEBFUN:BLOCKCOEFF:sum:notSupported', ...
                'Conversion of integration to coefficients is not supported.')
        end

        function E = feval(A, location, direction)
            error('CHEBFUN:BLOCKCOEFF:feval:notSupported', ...
                'Conversion of evaluation to coefficients is not supported.')
        end

        function F = inner(f)
            error('CHEBFUN:BLOCKCOEFF:inner:notSupported', ...
                'Conversion of inner product to coefficients is not supported.')
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function C = mtimes(A, B)

            if ( isnumeric(A) )
                % Allow multiplying a BLOCKCOEFF with a scalar.
                c = B.coeffs;
                for k = 1:numel(B.coeffs)
                    c{k} = A*c{k};
                end

                % Create the result.
                C = blockCoeff(c, B.domain);
                return
            elseif ( isnumeric(B) )
                C = mtimes(B, A);
                return
            end

            if ( isempty(A.coeffs) || isempty(B.coeffs) )
                C = blockCoeff([]);
                return
            end

            % Initialize constant term of A times the coeffs of B.
            c = B.coeffs;
            for k = 1:numel(B.coeffs)
                c{k} = A.coeffs{end}.*c{k};
            end

            % Do a convolution-style multiplication for the rest of A.
            if ( numel(A.coeffs) > 1 )
                z = {0*c{1}};  % an appropriate zero chebfun
            end
            for j = 1:numel(A.coeffs)-1
                B = diff(B);  % differentiate this operator
                c = [z, c];
                for k = 1:numel(B.coeffs)
                    c{k} = c{k} + A.coeffs{end-j}.*B.coeffs{k};
                end
            end

            % Create the result.
            C = blockCoeff(c, A.domain);
        end

        function C = plus(A, B)
            sA = numel(A.coeffs);
            sB = numel(B.coeffs);

            % Empty operand returns empty result.
            if ( (sA == 0) || (sB == 0) )
                C = blockCoeff([]);
                return
            end

            % Get to same length.
            if ( sA < sB )
                z = {0*A.coeffs{1}};
                A.coeffs = [repmat(z, 1, sB-sA), A.coeffs];
                sA = sB;
            elseif ( sB < sA )
                z = {0*B.coeffs{1}};
                B.coeffs = [repmat(z, 1, sA-sB), B.coeffs];
            end

            % Do the work.
            c = cell(1, sA);
            for k = 1:sA
                c{k} = A.coeffs{k} + B.coeffs{k};
            end
            C = blockCoeff(c, A.domain);
        end

        function A = uminus(A)
            A.coeffs = cellfun(@uminus, A.coeffs, 'uniform', false);
         end

        function A = uplus(A)
        end

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
x=@(s,t) 3*cos(s);
y=@(s,t) 2*sin(s);
z=@(s,t) t;
ezmesh(x,y,z)

author: Eurus-Holmes


text:
%%*****************************************************************
%% validate: validate data
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function  [blk,At,C,b,dim,nnblk,parbarrier] = ...
              validate(blk,At,C,b,par,parbarrier);

   if (nargin >= 5)
      spdensity = par.spdensity;
   else
      spdensity = 0.4;
   end
%%
   if ~iscell(blk);
      error('validate: blk must be a cell array'); end;
   if (size(blk,2) < 2)
      error('validate: blk must be a cell array with at least 2 columns');
   end
   if ~iscell(At) | ~iscell(C);
      error('validate: At, C must be cell arrays');
   end
   if (size(At,1) ~= size(blk,1))
      if (size(At,2) == size(blk,1));
         At = At';
      else
         error('validate: size of At is not compatible with blk');
      end
   end
   if (size(C,1) ~= size(blk,1))
      if (size(C,2) == size(blk,1))
         C = C';
      else
         error('validate: size of C is not compatible with blk');
      end
   end
   if (min(size(b)) > 1); error('validate: b must be a vector'); end
   if (size(b,1) < size(b,2)); b = b'; end
   m = length(b);
%%
%%-----------------------------------------
%% validate blk, At, C
%%-----------------------------------------
%%
   for p = 1:size(blk,1)
      if (size(blk{p,2},1) > size(blk{p,2},2))
         blk{p,2} = blk{p,2}';
      end
      pblk = blk(p,:);
      n = sum(pblk{2});
      numblk = length(pblk{2});
      if strcmp(pblk{1},'s');
         m1 = size(At{p,1},2);
         n2 = sum(pblk{2}.*pblk{2});  n22 = sum(pblk{2}.*(pblk{2}+1))/2;
         ntotal(p) = n22;
         if ~all(size(C{p}) == n)
            error('validate: blk and C are not compatible');
         end
         if (norm(C{p}-C{p}',inf) > 1e-13*norm(C{p},inf));
            error('validate: C is not symmetric');
         end
         if (size(At{p,1}) == [m1, n22] & m1~=n22);
            At{p,1} = At{p,1}';
         end
         if (~isempty(At{p,1})) & (size(At{p,1},1) ~= n22)
            error('validate: blk and At not compatible');
         end
         if (nnz(At{p,1}) < spdensity*n22*m1)
            if ~issparse(At{p,1}); At{p,1} = sparse(At{p,1}); end
         end
         if (length(pblk) > 2) %% for low rank constraints
            len = sum(pblk{3});
            if (size(pblk{1,3},2) < size(pblk{1,3},1))
                blk{p,3} = blk{p,3}';
            end
            if any(size(At{p,2})- [n,len])
               error(' low rank structure specified in blk and At not compatible')
            end
            if (length(At(p,:)) > 2)  & ~isempty(At{p,3})
               if all(size(At{p,3},2)-[1,4])
                  error(' low rank structure in At{p,3} not specified correctly')
               end
               if (size(At{p,3},2) == 1)
                  if (size(At{p,3},1) < size(At{p,3},2)); At{p,3} = At{p,3}'; end
                  lenn = length(At{p,3});
                  constrnum = mexexpand(pblk{3},[1:length(pblk{3})]');
                  At{p,3} = [constrnum, [1:lenn]', [1:lenn]', At{p,3}];
               elseif (size(At{p,3},2) == 4)
                  dd = At{p,3};
                  [dummy,idxsort] = sort(dd(:,1));
                  dd = dd(idxsort,:);
                  lenn = size(dd,1);
                  idxD = [0; find(diff(dd(:,1))); lenn];
                  ii = zeros(lenn,1); jj = zeros(lenn,1);
                  ss = [0,cumsum(pblk{3})];
                  for k = 1:length(pblk{3})
                     idx = [idxD(k)+1 : idxD(k+1)];
                     ii(idx) = dd(idx,2)+ss(k); %% convert to cumulative indexing
                     jj(idx) = dd(idx,3)+ss(k);
                  end
                  At{p,3} = [dd(:,1),ii,jj,dd(:,4)];
               end
            else
               constrnum = mexexpand(pblk{3},[1:length(pblk{3})]');
               At{p,3} = [constrnum,  [1:len]', [1:len]', ones(len,1)];
            end
         end
         if (nnz(C{p}) < spdensity*n2) | (numblk > 1);
            if ~issparse(C{p}); C{p} = sparse(C{p}); end;
         else
            if issparse(C{p}); C{p} = full(C{p}); end;
         end
      elseif strcmp(pblk{1},'q') | strcmp(pblk{1},'l') | strcmp(pblk{1},'u');
         ntotal(p) = n;
         if (min(size(C{p})) ~= 1 | max(size(C{p})) ~= n);
            error(['validate: blk and C are not compatible']);
         end;
         if (size(C{p},1) < size(C{p},2)); C{p} = C{p}'; end
         if (size(At{p,1}) == [m, n] & m~=n);
            At{p,1} = At{p,1}';
         end
         if ~all(size(At{p,1}) == [n,m]);
            error('validate: blk and At not compatible');
         end
         if ~issparse(At{p,1});
            At{p,1} = sparse(At{p,1});
         end
         if (nnz(C{p}) < spdensity*n);
            if ~issparse(C{p}); C{p} = sparse(C{p}); end;
         else
            if issparse(C{p}); C{p} = full(C{p}); end;
         end;
      else
         error(' blk: some fields are not specified correctly');
      end
   end
   if (sum(ntotal) < m)
      error(' total dimension of C should be > length(b)');
   end
%%
%%-----------------------------------------
%% problem dimension
%%-----------------------------------------
%%
   dim = zeros(1,4);  nnblk = zeros(1,2);
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      if strcmp(pblk{1},'s')
         dim(1) = dim(1) + sum(pblk{2});
         nnblk(1) = nnblk(1) + length(pblk{2});
         nn(p) = sum(pblk{2});
      elseif strcmp(pblk{1},'q')
         dim(2) = dim(2) + sum(pblk{2});
         nnblk(2) = nnblk(2) + length(pblk{2});
         nn(p) = length(pblk{2});
      elseif strcmp(pblk{1},'l')
         dim(3) = dim(3) + sum(pblk{2});
         nn(p) = sum(pblk{2});
      elseif strcmp(pblk{1},'u')
         dim(4) = dim(4) + sum(pblk{2});
         nn(p) = sum(pblk{2});
      end
   end
%%
%%-----------------------------------------
%% validate parbarrier
%%-----------------------------------------
%%
    if (nargin == 6)
       if ~iscell(parbarrier);
	  if (length(parbarrier) == size(blk,1))
	     tmp = parbarrier;
             clear parbarrier;
	     parbarrier = cell(size(blk,1),1);
	     for p = 1:size(blk,1)
		parbarrier{p} = tmp(p);
             end
          end
       end
       if (size(parbarrier,2) > size(parbarrier,1))
          parbarrier = parbarrier';
       end
       for p = 1:size(blk,1)
          pblk = blk(p,:);
          if (size(parbarrier{p},1) > size(parbarrier{p},2))
             parbarrier{p} = parbarrier{p}';
          end
          len = length(parbarrier{p});
          if strcmp(pblk{1},'s') | strcmp(pblk{1},'q')
             if (len == 1)
                parbarrier{p} = parbarrier{p}*ones(1,length(pblk{2}));
	     elseif (len == 0)
                parbarrier{p} = zeros(1,length(pblk{2}));
	     elseif (len ~= length(pblk{2}))
                error('blk and parbarrier not compatible');
             end
          elseif strcmp(pblk{1},'l')
             if (len == 1)
                parbarrier{p} = parbarrier{p}*ones(1,sum(pblk{2}));
	     elseif (len == 0)
		parbarrier{p} = zeros(1,sum(pblk{2}));
             elseif (len ~= sum(pblk{2}))
                error('blk and parbarrier not compatible');
             end
          elseif strcmp(pblk{1},'u')
             parbarrier{p}= zeros(1,sum(pblk{2}));
          end
       end
    end
%%***********************************************************************

author: yu-jiang

repo: radpbook

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


text:
%% Basic Usage
% This document introduces the basic usage of MaxwellFDFD with example codes.

%% Example Code
% We start with a simple example code that solves a 2D problem, which is to
% examine the transmission of a plane wave through a narrow slit:
[E, H] = maxwell_run(...
	'OSC', 1e-9, 1550, ...
	'DOM', {'vacuum', 'none', 1.0}, [-1100 1100; -1100 2600; 0 10], 10, BC.p, [100 100 0],...
	'OBJ', {'CRC/Ag', 'k'}, ...
		Box([-1100 -80; 0 1000; 0 10]), ...
		Box([80 1100; 0 1000; 0 10]), ...
	'SRCJ', PlaneSrc(Axis.y, -500, Axis.x));

%%%
% The above is actually one line of code; the ellipsis |...| at the end of each
% line is MATLAB's way of indicating the continuation of a line.  Upon
% execution, this code calculates the solution _E_- and _H_-fields of the
% problem.  The solutions are stored in the variables |E| and |H|, which are the
% outputs of <../comp/maxwell_run.html |maxwell_run|>, MaxwellFDFD's core
% function, that takes a long list of arguments describing the problem.

%% Visualizing the Solution
% Before discussing the details of the arguments of |maxwell_run|, let's first
% visualize the solution to understand the problem we are solving better. This
% can be done by another simple line of code:
vis2d(E{Axis.x}, Axis.z, 5);

%%%
% which produces a figure that looks like:
%
% <<../img/basic_01.png>>
%
% Note that the PML regions are excluded from the plot by default.

%% Visualizing of Objects and Sources
% The above field plot would have been much more informative if it was overlaid
% with the objects and source placed in the simulation domain.  The arrays of
% the objects and sources can be obtained from |maxwell_run| by changing the
% output arguments as:
[E, H, obj_array, src_array] = maxwell_run({ARGUMENTS});

%%%
% These objects and sources (actually only one source for the present problem)
% are drawn on top of the field plot when they are supplied to |vis2d| as extra
% input arguments as
vis2d(E{Axis.x}, Axis.z, 5, obj_array, src_array);

%%%
% The modified code will produce an updated figure that looks like:
%
% <<../img/basic_02.png>>
%
% The black lines indicate the boundaies of the metal pieces, and green line
% indicates the location of the source plane.

%% Input Arguments of |maxwell_run|
% |maxwell_run| takes many input arguments that are grouped into several
% *parameter groups*.  Each parameter group is specified by |'NAME'| (such as
% |'OSC'|, |'DOM'|, |'OBJ'|, |'SRCJ'|) that is followed by the arguments in the
% parameter group.  Below, the meanings of the arguments used in the above
% example code are explained for each parameter group.
%
% *'OSC'* specifies the oscillation parameter group, which describes the
% wavelength of the source.
%
% * |1e-9|: the unit of wavelength is 1 nm (= 1e-9 m).
% * |1550|: the wavelength is 1550 nm.
%
% The unit of wavelength specified in this parameter group serves as the unit of
% all the length arguments of |maxwell_run|.
%
% *'DOM'* specifies the domain parameter group, which describes the simulation
% domain of the problem.
%
% * |{'vacuum', 'none', 1.0}|: the background material that fills the simulation
% domain is named |'vacuum'| (just a user-defined name), will be visualized with
% color |'none'| (i.e., the material will not be shown in the figure generated
% by |vis2d|), and its dielectric constant is 1.0.
% * |[-1100 1100; -1100 2600; 0 10]|: the simulation domain lies between -1100
% and 1100 nm in the _x_-direction, -1100 and 2600 nm in the _y_-direction, 0
% and 10 nm in the _z_-direction.
% * |10|: the coarsest grid size used to descritize the simulation domain is 10
% nm.
% * |BC.p|: periodic boundary conditions are used in all the _x_-, _y_-,
% _z_-directions.
% * |[100 100 0]|: the thicknesses of PML are 100 nm at the _x_- and _y_-normal
% boundaries, and no PML is used at the _z_-normal boundaries.  The simulation
% domain except PML lies between -1000 and 1000 nm in the _x_-direction, -1000
% and 2500 nm in the _y_-direction, 0 and 10 nm in the _z_-direction.
%
% *'OBJ'* specifies the object parameter group, which describes the objects
% placed in the simulation domain.
%
% * |{'CRC/Ag', 'k'}|: the material used to create objects has
% frequency-dependent dielectric constants defined in the file |Ag.mat| under
% the subdirectory |dielconst/CRC/| under the main MaxwellFDFD directory.  (In
% fact, the dielectric constants are taken from the silver data in the CRC
% handbook.)  When the objects made of this material are plotted, they will be
% colored black (whose MATLAB color code is |'k'|).
% * |Box([-1100 -80; 0 1000; 0 10])|: place a box made of the material.  The
% box lies between -1100 and -80 nm in the _x_-direction, 0 and 1000 nm in
% the _y_-direction, 0 and 10 nm in the _z_-direction.
% * |Box([80 1100; 0 1000; 0 10])|: place another box made of the material.
%
% *'SRCJ'* specifies the _J_ source parameter group, which describes the
% electric current sources placed in the simulation domain.
%
% * |PlaneSrc(Axis.y, -500, Axis.x))|: at _y_ = -500 nm, place a plane of
% electric dipoles oscillating in the _x_-direction.
%
% Note that the order of appearance of the parameter groups in |maxwell_run| is
% interchangeable; for example, you can put the source parameter group before
% the object parameter group, and vice versa.

%% Inspecting the Design before Solution
% The above example is simple, so it does not hurt to solve the problem without
% making sure that the system you described with the arguments of |maxwell_run|
% is indeed what you intended to examine.  However, for larger and more complex
% problems, the solution process can be time-consuming.  Therefore, in general
% it is a good practice to inspect the design before solving the problem.
%
% You can command |maxwell_run| to visualize the simulation domain without
% solving the problem by appending a logical argument (|true| or |false|) as:
[E, H, obj_array, src_array] = maxwell_run({PARAMETER GROUPS}, true);

%%%
% The above code does not return solution _E_- and _H_-fields in the output
% arguments |E| and |H|, i.e., |E| and |H| are empty, but it produces the
% following figure:
%
% <<../img/basic_03.png>>
%
% The above figure visualizes the objects and source in the simulation domain,
% including the PML regions.

%% Complete Code
% Below is the complete code that has everything discussed in this document.
% Note that |vis2d| works only if the last argument of |maxwell_run| is flipped
% to |false| to calculate the solution.
[E, H, obj_array, src_array] = maxwell_run(...
	'OSC', 1e-9, 1550, ...
	'DOM', {'vacuum', 'none', 1.0}, [-1100 1100; -1100 2600; 0 10], 10, BC.p, [100 100 0],...
	'OBJ', {'CRC/Ag', 'k'}, ...
		Box([-1100 -80; 0 1000; 0 10]), ...
		Box([80 1100; 0 1000; 0 10]), ...
	'SRCJ', PlaneSrc(Axis.y, -500, Axis.x), ...
	false);  % true to inspect arguments without solving equation

vis2d(E{Axis.x}, Axis.z, 5, obj_array, src_array);

%%% See Also
% <intermediate.html |Intermediate Usage|>, <../comp/maxwell_run.html |maxwell_run|>

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd