text:
function write_brainvoyager(filename, data, dataformat, vmpversion)

% helper function to write volumetric data for brainvoyager.
% this is old code that moved from ft_volumewrite to clean up
% the high level function a bit. it is assumed that the orientation
% of the volume is correct.

data = double(data);
siz  = size(data);
maxval = max(data(:));

[p, f, e] = fileparts(filename);
filename  = fullfile(p, f);

switch dataformat
  case 'vmp'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % write in BrainVoyager VMP format
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fid = fopen(sprintf('%s.vmp', filename),'w');
    if fid < 0
      ft_error('Cannot write to file %s.vmp\n',filename);
    end

    switch vmpversion
      case 1
        % write the header
        fwrite(fid, 1, 'short');      % version
        fwrite(fid, 1, 'short');      % number of maps
        fwrite(fid, 1, 'short');      % map type
        fwrite(fid, 0, 'short');      % lag

        fwrite(fid, 0, 'short');      % cluster size
        fwrite(fid, 1, 'float');      % thresh min
        fwrite(fid, maxval, 'float'); % thresh max
        fwrite(fid, 0, 'short');      % df1
        fwrite(fid, 0, 'short');      % df2
        fwrite(fid, 0, 'char');       % name

        fwrite(fid, siz, 'short');    % size
        fwrite(fid, 0, 'short');
        fwrite(fid, siz(1)-1, 'short');
        fwrite(fid, 0, 'short');
        fwrite(fid, siz(2)-1, 'short');
        fwrite(fid, 0, 'short');
        fwrite(fid, siz(3)-1, 'short');
        fwrite(fid, 1, 'short');      % resolution

        % write the data
        fwrite(fid, data, 'float');
      case 2
        % determine relevant subvolume
        % FIXME, this is not functional at the moment, since earlier in this function all nans have been replaced by zeros
        minx = min(find(~isnan(max(max(data,[],3),[],2))));
        maxx = max(find(~isnan(max(max(data,[],3),[],2))));
        miny = min(find(~isnan(max(max(data,[],3),[],1))));
        maxy = max(find(~isnan(max(max(data,[],3),[],1))));
        minz = min(find(~isnan(max(max(data,[],1),[],2))));
        maxz = max(find(~isnan(max(max(data,[],1),[],2))));

        % write the header
        fwrite(fid, 2, 'short');      % version
        fwrite(fid, 1, 'int');        % number of maps
        fwrite(fid, 1, 'int');        % map type
        fwrite(fid, 0, 'int');        % lag

        fwrite(fid, 0, 'int');        % cluster size
        fwrite(fid, 0, 'char');       % cluster enable
        fwrite(fid, 1, 'float');      % thresh
        fwrite(fid, maxval, 'float'); % thresh
        fwrite(fid, 0, 'int');        % df1
        fwrite(fid, 0, 'int');        % df2
        fwrite(fid, 0, 'int');        % bonf
        fwrite(fid, [255,0,0], 'uchar');   % col1
        fwrite(fid, [255,255,0], 'uchar'); % col2
        fwrite(fid, 1, 'char');       % enable SMP
        fwrite(fid, 1, 'float');      % transparency
        fwrite(fid, 0, 'char');       % name

        fwrite(fid, siz, 'int');      % original size
        fwrite(fid, minx-1, 'int');
        fwrite(fid, maxx-1, 'int');
        fwrite(fid, miny-1, 'int');
        fwrite(fid, maxy-1, 'int');
        fwrite(fid, minz-1, 'int');
        fwrite(fid, maxz-1, 'int');
        fwrite(fid, 1, 'int');        % resolution

        % write the data
        fwrite(fid, data(minx:maxx,miny:maxy,minz:maxz), 'float');
    end
    fclose(fid);

  case 'vmr'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % write in BrainVoyager VMR format
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fid = fopen(sprintf('%s.vmr',filename),'w');
    if fid < 0
      ft_error('Cannot write to file %s.vmr\n',filename);
    end

    % data should be scaled between 0 and 225
    data = data - min(data(:));
    data = round(225*data./max(data(:)));

    % write the header
    fwrite(fid, siz, 'ushort');
    % write the data
    fwrite(fid, data, 'uint8');
    fclose(fid);
end


text:
function [ar,ma,ASAsellog,ASAcontrol]=armasel_rs(varargin)
%ARMASEL_RS ARMAsel model identification
%   [AR,MA,SELLOG] = ARMASEL_RS(ARH,N_OBS) estimates autoregressive, moving
%   average, and autoregressive moving average models from the high-order
%   AR model ARH and selects the model with optimal predictive qualities.
%   ARH has been estimated from N_OBS observations. The AR and MA parts of
%   the selected model, each possibly of order 0, are returned in the
%   parameter vectors AR and MA. The structure SELLOG provides additional
%   information on the selection process.
%
%   SELLOG contains the fields 'ar', 'ma' and 'arma', in which SELLOG
%   structures are nested, as returned by the functions ARH2AR, ARH2MA
%   and ARH2ARMA, invoked by ARMASEL. In the field 'armasel' a structure
%   is nested that reports information about the final stage of model
%   selection, where the preselected AR, MA and ARMA models are compared.
%
%   ARMASEL(ARH,N_OBS,CAND_AR_ORDER,CAND_MA_ORDER,CAND_ARMA_ORDER)
%   narrows the selection to candidate models with orders provided by the
%   rows CAND_AR_ORDER, CAND_MA_ORDER and CAND_ARMA_ORDER AR, MA and
%   ARMA(r,r-1) models independently.
%   For any of these arguments it is allowed to pass an
%   empty array. Alternatively, additional arguments may be omitted from
%   the input list. In both cases, default values are automatically
%   determined and substituted for the missing arguments. The functions
%   ARH2AR, ARH2MA and ARH2ARMA provide additional information on
%   defining candidate orders. Note in this respect, that the candidate
%   AR orders of the ARMA model are called CAND_ARMA_ORDER in this help
%   text, while in ARH2ARMA they are called CAND_AR_ORDER.
%
%   The selection of MA and ARMA models can be conditioned to the
%   selection of AR models from a specific set of candidate orders
%   CAND_AR_ORDER. See ASAGLOB_AR_COND for more information.
%
%   ARMASEL_RS is an ARMASA_RS main function.
%   ARMASEL_RS can also be used to selected models based on an ARH model
%   that has been estimated from segments of data.
%
%   See also: ARMASEL, ARH2AR, ARH2MA, ARH2ARMA, DATA_SEGMENTS.

%   Reference:  P. M. T. Broersen and S. de Waele, Selection of Order and
%               Type of Time Series Models Estimated from Reduced Statistics
%               Proceedings of SYSID 2002, May 2002.

%Header
%===================================================================================================

%Declaration of variables
%------------------------

%Declare and assign values to local variables
%according to the input argument pattern
[ar_rs,n_obs,cand_ar_order,cand_ma_order,cand_arma_order,arma_order_diff,ASAcontrol] = ASAarg(varargin, ...
{'ar_rs'       ;'n_obs'      ;'cand_ar_order';'cand_ma_order';'cand_arma_order';'arma_order_diff';'ASAcontrol'}, ...
{'isnumeric'   ;'isnumeric'  ;'isnumeric'    ;'isnumeric'    ;'isnumeric'      ;'isnumeric'      ;'isstruct'  }, ...
{'ar_rs'       ;'n_obs'                                                                                }, ...
{'ar_rs'       ;'n_obs'      ;'cand_ar_order'                                                                 }, ...
{'ar_rs'       ;'n_obs'      ;'cand_ar_order';'cand_ma_order'                                                 }, ...
{'ar_rs'       ;'n_obs'      ;'cand_ar_order';'cand_ma_order';'cand_arma_order'                               }, ...
{'ar_rs'       ;'n_obs'      ;'cand_ar_order';'cand_ma_order';'cand_arma_order';'arma_order_diff'             });

%Declare ASAglob variables
ASAglob = {'ASAglob_subtr_mean';'ASAglob_mean_adj';'ASAglob_rc';'ASAglob_ar';'ASAglob_final_f'; ...
      'ASAglob_final_b';'ASAglob_ar_cond'};

%Assign values to ASAglob variables by screening the
%caller workspace
for ASAcounter = 1:length(ASAglob)
   ASAvar = ASAglob{ASAcounter};
   eval(['global ' ASAvar]);
   if evalin('caller',['exist(''' ASAvar ''',''var'')'])
      eval([ASAvar '=evalin(''caller'',ASAvar);']);
   else
      eval([ASAvar '=[];']);
   end
end

%ARMASA-function version information
%-----------------------------------

%This ARMASA-function is characterized by
%its current version,
ASAcontrol.is_version = [2000 12 30 20 0 0];
%and its compatability with versions down to,
ASAcontrol.comp_version = [2000 12 30 20 0 0];

%This function calls other functions of the ARMASA
%toolbox. The versions of these other functions
%must be greater than or equal to:
ASAcontrol.req_version.sig2ar = [2000 12 30 20 0 0];
ASAcontrol.req_version.sig2ma = [2000 12 30 20 0 0];
ASAcontrol.req_version.sig2arma = [2000 12 30 20 0 0];

%Checks
%------

if ~isfield(ASAcontrol,'error_chk') | ASAcontrol.error_chk
      %Perform standard error checks
   %Input argument format checks
   ASAcontrol.error_chk = 1;
   if ~isnum(ar_rs)
      error(ASAerr(11,'ar_rs'))
   elseif ~isvector(ar_rs)
      error([ASAerr(14) ASAerr(15,'ar_rs')])
   end
   if ~isempty(cand_ar_order)
      if ~isnum(cand_ar_order) | ~isintvector(cand_ar_order) |...
            cand_ar_order(1)<0 | ~isascending(cand_ar_order)
         error(ASAerr(12,{'candidate';'cand_ar_order'}))
      elseif size(cand_ar_order,1)>1
         cand_ar_order = cand_ar_order';
         warning(ASAwarn(25,{'column';'cand_ar_order';'row'},ASAcontrol))
      end
   end
   if ~isempty(cand_ma_order)
      if ~isnum(cand_ma_order) | ~isintvector(cand_ma_order) |...
            cand_ma_order(1)<0 | ~isascending(cand_ma_order)
         error(ASAerr(12,{'candidate';'cand_ma_order'}))
      elseif size(cand_ma_order,1)>1
         cand_ma_order = cand_ma_order';
         warning(ASAwarn(25,{'column';'cand_ma_order';'row'},ASAcontrol))
      end
   end
   if ~isempty(cand_arma_order)
      if ~isnum(cand_arma_order) | ~isintvector(cand_arma_order) |...
            cand_arma_order(1)<0 | ~isascending(cand_arma_order)
         error(ASAerr(12,{'candidate';'cand_arma_order'}))
      elseif size(cand_arma_order,1)>1
         cand_arma_order = cand_arma_order';
         warning(ASAwarn(25,{'column';'cand_arma_order';'row'},ASAcontrol))
      end
   end
   if ~isempty(arma_order_diff) & ...
         (~isnum(arma_order_diff) | ...
         ~isintscalar(arma_order_diff) |...
         arma_order_diff<0)
      error(ASAerr(17,'arma_order_diff'))
   end

   %Input argument value checks
   if ~isreal(ar_rs)
      error(ASAerr(13))
   end
   if max(cand_ar_order) > length(ar_rs)-1
      error(ASAerr(37,'cand_ar_order'))
   end
   if max(cand_ma_order) > length(ar_rs)-1
      error(ASAerr(37,'cand_ma_order'))
   end
   if ~isempty(cand_arma_order) & ...
         ~isempty(arma_order_diff)
      if cand_arma_order(1)~=0 & ...
            (arma_order_diff < 1 | ...
            arma_order_diff > cand_arma_order(1))
         error(ASAerr(18,{'arma_order_diff';'1';...
               num2str(cand_arma_order(1))}))
      elseif length(cand_arma_order)>1 & ...
            (arma_order_diff < 1 | ...
            arma_order_diff > cand_arma_order(2))
         error(ASAerr(18,{'arma_order_diff';'1';...
               num2str(cand_arma_order(2))}))
      end
   end
end

if ~isfield(ASAcontrol,'version_chk') | ...
      ASAcontrol.version_chk %Perform version check
   ASAcontrol.version_chk = 1;

   %Make sure the requested version of this function
   %complies with its actual version
   ASAversionchk(ASAcontrol);

   %Make sure the requested versions of the called
   %functions comply with their actual versions
   arh2ar(ASAcontrol);
   arh2ma(ASAcontrol);
   arh2arma(ASAcontrol);
end

if ~isfield(ASAcontrol,'run') | ASAcontrol.run
   ASAcontrol.run = 1;
   ASAdate = now;
end

if ASAcontrol.run %Run the computational kernel
   ASAcontrol.version_chk = 0;
   ASAcontrol.error_chk = 0;

%Main
%=====================================================

%Initialization of variables
%---------------------------

if ASAglob_ar_cond
   ASAglob_ar_cond = 1;
else
   ASAglob_ar_cond = 0;
end

n_obs_tot = sum(n_obs);

%AR-, MA- and ARMA-model identification
%--------------------------------------

[ar_ar,ar_sellog] = arh2ar...
   (ar_rs,n_obs,cand_ar_order,ASAcontrol);
[ma_ma,ma_sellog] = arh2ma...
   (ar_rs,n_obs,cand_ma_order,ASAcontrol);
[arma_ar,arma_ma,arma_sellog] = arh2arma...
   (ar_rs,n_obs,cand_arma_order,arma_order_diff,ASAcontrol);

%Selection of the ARMAsel model
%------------------------------

%Asess the selected model orders
order = ...
   [length(ar_ar)-1;length(ma_ma)-1;length(arma_ar)-1];

%Asess the corresponding prediction error estimates
sel_location = ...
   [find(order(1) == ar_sellog.cand_order);...
    find(order(2) == ma_sellog.cand_order);...
    find(order(3) == arma_sellog.cand_ar_order)];
pe_est = [ar_sellog.pe_est(sel_location(1));...
      ma_sellog.pe_est(sel_location(2));...
      arma_sellog.pe_est(sel_location(3))];

%Select the model with the smallest prediction error
%estimate
[sel_pe_est,model] = min(pe_est);

%Arranging output arguments
%--------------------------

%Retrieve the parameters of the selected model
ar = 1;
ma = 1;
switch model
case 1 %The AR model has been selected
   ar = ar_ar;
case 2 %The MA model has been selected
   ma = ma_ma;
case 3 %The ARMA model has been selected
   ar = arma_ar;
   ma = arma_ma;
end

%Gernerate a structure variable ASAsellog to report
%the selection process
ASAsellog.funct_name = mfilename;
ASAsellog.funct_version = ASAcontrol.is_version;
ASAsellog.date_time = ...
   [datestr(ASAdate,8) 32 datestr(ASAdate,0)];
ASAsellog.comp_time = ar_sellog.comp_time+...
   ma_sellog.comp_time+...
   arma_sellog.comp_time;
ASAsellog.armasel.ar = ar;
ASAsellog.armasel.ma = ma;
ASAsellog.armasel.ar_pe_est = pe_est(1);
ASAsellog.armasel.ma_pe_est = pe_est(2);
ASAsellog.armasel.arma_pe_est = pe_est(3);
ASAsellog.armasel.ar_cond = ASAglob_ar_cond;
ASAsellog.armasel.mean_adj = ASAglob_mean_adj;
ASAsellog.ar = ar_sellog;
ASAsellog.ma = ma_sellog;
ASAsellog.arma = arma_sellog;

%Footer
%=====================================================

else %Skip the computational kernel
   %Return ASAcontrol as the first output argument
   if nargout>1
      warning(ASAwarn(9,mfilename,ASAcontrol))
   end
   ar = ASAcontrol;
   ASAcontrol = [];
end

%Program history
%======================================================================
%
% Version                Programmer(s)          E-mail address
% -------                -------------          --------------
% former versions        P.M.T. Broersen        broersen@tn.tudelft.nl
%                        S. de Waele            waele@tn.tudelft.nl
% [2000 11  1 12 0 0]    W. Wunderink           wwunderink01@freeler.nl
% [2000 12 30 20 0 0]         ,,                          ,,

author: Sable

repo: mcbench-benchmarks


text:
function [charAcc, wordAcc] = ScoreModel (words, imageModel, pairwiseModel, tripletList)
% This function runs the Markov network model end-to-end and computes the
% per-character and per-word accuracy on provided data.
%
% Input:
%   words: A cell array where words{i} is the struct array for the ith
%     word (this is the structure of the provided 'allWords' data).
%   imageModel: The provided image model struct.
%   pairwiseModel: A K-by-K matrix (K is the alphabet size) where pairwiseModel(i,j)
%     is the factor value for the pairwise factor of character i followed by
%     character j.
%   tripletModel: The array of character triplets we will consider (along
%     with their corresponding factor values).
%
% Output:
%   charAcc: The percentage of all characters (across all words) correctly
%     identified. (Between 0 and 1)
%   wordAcc: The percentage of the words in which every character is
%     correctly identified. (Between 0 and 1)
%
% Copyright (C) Daphne Koller, Stanford University, 2012

predictions = ComputeWordPredictions(words, imageModel, pairwiseModel, tripletList);
[charAcc, wordAcc] = ScorePredictions(words, predictions, true);

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function obj = saveExpression(obj, export_path, varargin)
    % save the symbolic expression of system dynamical equations to a MX
    % binary files
    %
    % Parameters:
    %  export_path: the path to export the file @type char
    %  varargin: variable input parameters @type varargin
    %   ForceExport: force the export @type logical

    input_categories = {'Control','ConstraintWrench','External'};
    for k=1:3
        category = input_categories{k};
        % export the input map
        gmap_funcs = fieldnames(obj.Gmap.(category));
        if ~isempty(gmap_funcs)
            for i=1:length(gmap_funcs)
                fun = gmap_funcs{i};
                if ~isempty(obj.Gmap.(category).(fun))
                    save(obj.Gmap.(category).(fun),export_path,varargin{:});
                end
            end
        end

        % export the input vector fields
        gvec_funcs = fieldnames(obj.Gvec.(category));
        if ~isempty(gvec_funcs)
            for i=1:length(gvec_funcs)
                fun = gvec_funcs{i};
                if ~isempty(obj.Gvec.(category).(fun))
                    save(obj.Gvec.(category).(fun),export_path,varargin{:});
                end
            end
        end
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function params = updatePath(dirRoot,params)

if nargin < 1
    dirRoot = pwd;
end

if nargin < 2
    params = defaultParams(dirRoot);
end

params.trainingImages = [dirRoot '/Training/Images/'];
params.trainingExamples = [dirRoot '/Training/Images/Examples/'];
params.data = [dirRoot '/Data/'];
params.yourData = [dirRoot '/Data/yourData/'];
params.tempdir = [dirRoot '/tmpdir/'];
params.SS.soft_dir =[dirRoot '/segment/'];

author: Cloud-CV

repo: object-proposals


text:
function p00_newton_test ( problem_num )

%*****************************************************************************80
%
%% P00_NEWTON_TEST applies Newton's method to the perturbed starting point.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM_NUM, the number of problems.
%
  for problem = 1 : problem_num

    option_num = p00_option_num ( problem );

    fprintf ( 1, '\n' );

    for option = 1 : option_num

      seed = 123456789;

      fprintf ( 1, '\n' );
      fprintf ( 1, 'P00_NEWTON_TEST\n' );
      fprintf ( 1, '  Problem number = %d\n', problem );
      fprintf ( 1, '  Using option OPTION = %d\n', option );
%
%  Get the title.
%
      title = p00_title ( problem, option );

      fprintf ( 1, '  %s\n', title );
%
%  Get the number of variables.
%
      nvar = p00_nvar ( problem, option );

      fprintf ( 1, '  Number of variables is %d\n', nvar );
%
%  Get the starting point.
%
      x0 = p00_start ( problem, option, nvar );
%
%  Perturb the starting point.
%
      for i = 1 : nvar
        [ r, seed ] = r8_uniform_01 ( seed );
        dx = 0.10 * r * ( x0(i) + r8_sign ( x0(i) ) );
        x1(i) = x0(i) + dx;
      end
%
%  Choose a continuation parameter index.
%
      par_index = p00_par_index ( problem, option, nvar, x1 );

      fprintf ( 1, '  Fixing variable X(%d) = %f\n', par_index, x1(par_index) );
%
%  Apply Newton's method.
%
      [ x2, status ] = p00_newton ( problem, option, nvar, x1, par_index );

      if ( status == -3 )
        fprintf ( 1, '  The convergence test was not satisfied.\n' );
      elseif ( status == -2 )
        fprintf ( 1, '  The iteration seemed to be diverging, and was halted.\n' );
      elseif ( status == -1 )
        fprintf ( 1, '  The jacobian was singular, and the iteration was halted.\n' );
      else
        fprintf ( 1, '  Convergence was achieved in %d steps.\n', status );
      end

      if ( nvar <= 10 )

        fprintf ( 1, '\n' );
        fprintf ( 1, '        X0               X1=X0+dX          X2\n' );
        fprintf ( 1, '\n' );
        for i = 1 : nvar
          fprintf ( 1, '  %14f  %14f  %14f\n', x0(i), x1(i), x2(i) );
        end

      else

        fprintf ( 1, '\n' );
        fprintf ( 1, '      ||X0||           ||X1=X0+dX||        ||X2||\n' );
        fprintf ( 1, '\n' );
        fprintf ( 1, '  %14f  %14f  %14f\n', norm ( x0, 2 ), norm ( x1, 2 ), norm ( x2, 2 ) );

      end
%
%  Compute the initial function value.
%
      fx0 = p00_fun ( problem, option, nvar, x0 );
      fx1 = p00_fun ( problem, option, nvar, x1 );
      fx2 = p00_fun ( problem, option, nvar, x2 );

      if ( nvar <= 10 )

        fprintf ( 1, '\n' );
        fprintf ( 1, '       F(X0)           F(X1=X0+dX)        F(X2)\n' );
        fprintf ( 1, '\n' );
        for i = 1 : nvar - 1
          fprintf ( 1, '  %14f  %14f  %14f\n', fx0(i), fx1(i), fx2(i) );
        end

      else

        fprintf ( 1, '\n' );
        fprintf ( 1, '     ||F(X0)||       ||F(X1=X0+dX)||        ||F(X2)||\n' );
        fprintf ( 1, '\n' );
        fprintf ( 1, '  %14f  %14f  %14f\n', norm ( fx0, 2 ), norm ( fx1, 2 ), norm ( fx2, 2 ) );

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [logPr] = calcLogPrMatrixNormalInvWishart( A, invSigma, PP )

[logPrSigma, cholInvSigma] = calcLogPrInvWishart( invSigma, PP );

logDetInvSigma = 2*sum( log( diag( cholInvSigma ) ) );
logDetSigma    = -logDetInvSigma; % det(S) = 1/det(inv(S))

CC = PP.invAScaleMat;
[D DR] = size( CC );

if isfield( PP, 'MeanMat') && ~isempty( PP.MeanMat )
    UU = A - PP.MeanMat;
else
    UU = A; % assume zero mean
end

logPrA = 0.5*D*log( det(CC) ) - 0.5*DR*D*log( 2*pi ) ...
          - 0.5*DR*logDetSigma ...
          - 0.5*trace( UU'*invSigma*UU*CC );

logPr = logPrA + logPrSigma;

author: michaelchughes

repo: NPBayesHMM


text:
function v = isabspath(filename)

if isempty(filename),
  v = false;
  return;
end

if isunix,
  v = filename(1) == '/';
else
  v = numel(filename) > 1 && filename(2) == ':';
end

author: kristinbranson

repo: JAABA


text:
function [gx] = g_vgo(x,P,u,in)
gx = VBA_sigmoid(exp(P(1))*x+P(2));

author: MBB-team

repo: VBA-toolbox


text:
function gen_laguerre_rule ( order, alpha, a, b, filename )

%*****************************************************************************80
%
%% GEN_LAGUERRE_RULE generates a Gauss-Laguerre rule.
%
%  Discussion:
%
%    This program computes a standard or exponentially weighted
%    generalized Gauss-Laguerre quadrature rule and writes it to a file.
%
%    The user specifies:
%    * the ORDER (number of points) in the rule;
%    * ALPHA, the exponent of |X|;
%    * A, the left endpoint of integration;
%    * B, the scale factor in the exponential;
%    * FILENAME, the root name of the output files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GEN_LAGUERRE_RULE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute a generalized Gauss-Laguerre rule for approximating\n' );
  fprintf ( 1, '    Integral ( a <= x < oo ) |x-a|^ALPHA exp(-B*(x-a)) f(x) dx\n' );
  fprintf ( 1, '  of order ORDER.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The user specifies ORDER, ALPHA, A, B, and FILENAME.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ORDER is the number of points.\n' );
  fprintf ( 1, '  ALPHA is the exponent of |X|.\n' );
  fprintf ( 1, '  A is the left endpoint (typically 0).\n' );
  fprintf ( 1, '  B is the exponential scale factor (typically 1).\n' );
  fprintf ( 1, '  FILENAME is used to generate 3 files:\n' );
  fprintf ( 1, '  * filename_w.txt - the weight file\n' );
  fprintf ( 1, '  * filename_x.txt - the abscissa file.\n' );
  fprintf ( 1, '  * filename_r.txt - the region file.\n' );
%
%  Initialize the parameters.
%
  beta = 0.0;
%
%  Get ORDER.
%
  if ( nargin < 1 )
    order = input ( '  Enter the rule order ORDER.' );
  elseif ( ischar ( order ) )
    order = str2num ( order );
  end
%
%  Get ALPHA.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  ALPHA is the exponent of |X| in the weighting function.\n' );
    fprintf ( 1, '  ALPHA is a real number strictly greater than -1.\n' );
    fprintf ( 1, '\n' );
    alpha = input ( '  Enter the value of ALPHA.' );
  elseif ( ischar ( alpha ) )
    alpha = str2num ( alpha );
  end
%
%  Get A.
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  A is the left endpoint, typically 0.\n' );
    fprintf ( 1, '\n' );
    a = input ( '  Enter the value of A.' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end
%
%  Get B.
%
  if (  nargin < 4 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  B is the exponential scale factor, typically 1.\n' );
    fprintf ( 1, '\n' );
    b = input ( '  Enter the value of B.' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end
%
%  Get FILENAME:
%
  if ( nargin < 5 )
    fprintf ( 1,  '\n' );
    fprintf ( 1,  '  FILENAME specifies the ''root name'' of the quadrature files).\n' );
    filename = input ( '  Enter FILENAME as a quoted string:' );
  end
%
%  Input summary.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ORDER = %d\n', order );
  fprintf ( 1, '  ALPHA = %f\n', alpha );
  fprintf ( 1, '  A = %f\n', a );
  fprintf ( 1, '  B = %f\n', b );
  fprintf ( 1, '  FILENAME = "%s".\n', filename );
%
%  Construct the rule.
%
  kind = 5;
  [ x, w ] = cgqf ( order, kind, alpha, beta, a, b );
%
%  Write the rule.
%
  r = zeros ( 2, 1 );
  r(1) = a;
  r(2) = r8_huge ( );
  rule_write ( order, filename, x, w, r );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GEN_LAGUERRE_RULE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ t, wts ] = cdgqf ( nt, kind, alpha, beta )

%*****************************************************************************80
%
%% CDGQF computes a Gauss quadrature formula with default A, B and simple knots.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with a classical weight function with default values for A and B,
%    and only simple knots.
%
%    There are no moments checks and no printing is done.
%
%    Use routine EIQFS to evaluate a quadrature computed by CGQFS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%
  parchk ( kind, 2 * nt, alpha, beta );
%
%  Get the Jacobi matrix and zero-th moment.
%
  [ aj, bj, zemu ] = class_matrix ( kind, nt, alpha, beta );
%
%  Compute the knots and weights.
%
  [ t, wts ] = sgqf ( nt, aj, bj, zemu );

  return
end
function [ t, wts ] = cgqf ( nt, kind, alpha, beta, a, b )

%*****************************************************************************80
%
%% CGQF computes knots and weights of a Gauss quadrature formula.
%
%  Discussion:
%
%    The user may specify the interval (A,B).
%
%    Only simple knots are produced.
%
%    The user may request that the routine print the knots and weights,
%    and perform a moment check.
%
%    Use routine EIQFS to evaluate this quadrature formula.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Compute the Gauss quadrature formula for default values of A and B.
%
  [ t, wts ] = cdgqf ( nt, kind, alpha, beta );
%
%  All knots have multiplicity = 1.
%
  mlt = zeros(nt,1);
  mlt(1:nt) = 1;
%
%  NDX(I) = I.
%
  ndx = ( 1 : nt );
%
%  Scale the quadrature rule.
%
  [ t, wts ] = scqf ( nt, t, mlt, wts, nt, ndx, kind, alpha, beta, a, b );

  return
end
function [ aj, bj, zemu ] = class_matrix ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% CLASS_MATRIX computes the Jacobi matrix for a quadrature rule.
%
%  Discussion:
%
%    This routine computes the diagonal AJ and subdiagonal BJ
%    elements of the order M tridiagonal symmetric Jacobi matrix
%    associated with the polynomials orthogonal with respect to
%    the weight function specified by KIND.
%
%    For weight functions 1-7, M elements are defined in BJ even
%    though only M-1 are needed.  For weight function 8, BJ(M) is
%    set to zero.
%
%    The zero-th moment of the weight function is returned in ZEMU.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, integer M, the order of the Jacobi matrix.
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real AJ(M), BJ(M), the diagonal and subdiagonal
%    of the Jacobi matrix.
%
%    Output, real ZEMU, the zero-th moment.
%
  temp = eps;

  parchk ( kind, 2 * m - 1, alpha, beta );

  temp2 = 0.5;

  if ( 500.0 * temp < abs ( ( gamma ( temp2 ) )^2 - pi ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CLASS - Fatal error!\n' );
    fprintf ( 1, '  Gamma function does not match machine parameters.\n' );
    error ( 'CLASS - Fatal error!' );
  end

  bj = zeros(m,1);
  aj = zeros(m,1);

  if ( kind == 1 )

    ab = 0.0;

    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod ( i, 2 );
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 2 )

    zemu = pi;

    aj(1:m) = 0.0;

    bj(1) =  sqrt ( 0.5 );
    bj(2:m) = 0.5;

  elseif ( kind == 3 )

    ab = alpha * 2.0;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 )^2 ...
      / gamma ( ab + 2.0 );

    aj(1:m) = 0.0;
    bj(1) = 1.0 / ( 2.0 * alpha + 3.0 );
    for i = 2 : m
      bj(i) = i * ( i + ab ) / ( 4.0 * ( i + alpha )^2 - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 4 )

    ab = alpha + beta;
    abi = 2.0 + ab;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 ) ...
      * gamma ( beta + 1.0 ) / gamma ( abi );
    aj(1) = ( beta - alpha ) / abi;
    bj(1) = 4.0 * ( 1.0 + alpha ) * ( 1.0 + beta ) ...
      / ( ( abi + 1.0 ) * abi * abi );
    a2b2 = beta * beta - alpha * alpha;

    for i = 2 : m
      abi = 2.0 * i + ab;
      aj(i) = a2b2 / ( ( abi - 2.0 ) * abi );
      abi = abi^2;
      bj(i) = 4.0 * i * ( i + alpha ) * ( i + beta ) * ( i + ab ) ...
        / ( ( abi - 1.0 ) * abi );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 5 )

    zemu = gamma ( alpha + 1.0 );

    for i = 1 : m
      aj(i) = 2.0 * i - 1.0 + alpha;
      bj(i) = i * ( i + alpha );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 6 )

    zemu = gamma ( ( alpha + 1.0 ) / 2.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      bj(i) = ( i + alpha * mod ( i, 2 ) ) / 2.0;
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 7 )

    ab = alpha;
    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod(i,2);
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 8 )

    ab = alpha + beta;
    zemu = gamma ( alpha + 1.0 ) * gamma ( - ( ab + 1.0 ) ) ...
      / gamma ( - beta );
    apone = alpha + 1.0;
    aba = ab * apone;
    aj(1) = - apone / ( ab + 2.0 );
    bj(1) = - aj(1) * ( beta + 1.0 ) / ( ab + 2.0 ) / ( ab + 3.0 );
    for i = 2 : m
      abti = ab + 2.0 * i;
      aj(i) = aba + 2.0 * ( ab + i ) * ( i - 1 );
      aj(i) = - aj(i) / abti / ( abti - 2.0 );
    end

    for i = 2 : m - 1
      abti = ab + 2.0 * i;
      bj(i) = i * ( alpha + i ) / ( abti - 1.0 ) * ( beta + i ) ...
        / ( abti^2 ) * ( ab + i ) / ( abti + 1.0 );
    end

    bj(m) = 0.0;
    bj(1:m) =  sqrt ( bj(1:m) );

  end

  return
end
function [ d, z ] = imtqlx ( n, d, e, z )

%*****************************************************************************80
%
%% IMTQLX diagonalizes a symmetric tridiagonal matrix.
%
%  Discussion:
%
%    This routine is a slightly modified version of the EISPACK routine to
%    perform the implicit QL algorithm on a symmetric tridiagonal matrix.
%
%    The authors thank the authors of EISPACK for permission to use this
%    routine.
%
%    It has been modified to produce the product Q' * Z, where Z is an input
%    vector and Q is the orthogonal matrix diagonalizing the input matrix.
%    The changes consist (essentialy) of applying the orthogonal transformations
%    directly to Z as they are generated.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%    Roger Martin, James Wilkinson,
%    The Implicit QL Algorithm,
%    Numerische Mathematik,
%    Volume 12, Number 5, December 1968, pages 377-383.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real D(N), the diagonal entries of the matrix.
%
%    Input, real E(N), the subdiagonal entries of the
%    matrix, in entries E(1) through E(N-1).
%
%    Input, real Z(N), a vector to be operated on.
%
%    Output, real D(N), the diagonal entries of the diagonalized matrix.
%
%    Output, real Z(N), the value of Q' * Z, where Q is the matrix that
%    diagonalizes the input symmetric tridiagonal matrix.
%
  itn = 30;

  prec = eps;

  if ( n == 1 )
    return
  end

  e(n) = 0.0;

  for l = 1 : n

    j = 0;

    while ( 1 )

      for m = l : n

        if ( m == n )
          break
        end

        if ( abs ( e(m) ) <= prec * ( abs ( d(m) ) + abs ( d(m+1) ) ) )
          break
        end

      end

      p = d(l);

      if ( m == l )
        break
      end

      if ( j == itn )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'IMTQLX - Fatal error!\n' );
        fprintf ( 1, '  Iteration limit exceeded.\n' );
        error ( 'IMTQLX - Fatal error!' );
      end

      j = j + 1;
      g = ( d(l+1) - p ) / ( 2.0 * e(l) );
      r =  sqrt ( g * g + 1.0 );
      g = d(m) - p + e(l) / ( g + r8_sign ( g ) * abs ( r ) );
      s = 1.0;
      c = 1.0;
      p = 0.0;
      mml = m - l;

      for ii = 1 : mml

        i = m - ii;
        f = s * e(i);
        b = c * e(i);

        if ( abs ( f ) >= abs ( g ) )
          c = g / f;
          r =  sqrt ( c * c + 1.0 );
          e(i+1) = f * r;
          s = 1.0 / r;
          c = c * s;
        else
          s = f / g;
          r =  sqrt ( s * s + 1.0 );
          e(i+1) = g * r;
          c = 1.0 / r;
          s = s * c;
        end

        g = d(i+1) - p;
        r = ( d(i) - g ) * s + 2.0 * c * b;
        p = s * r;
        d(i+1) = g + p;
        g = c * r - b;
        f = z(i+1);
        z(i+1) = s * z(i) + c * f;
        z(i) = c * z(i) - s * f;

      end

      d(l) = d(l) - p;
      e(l) = g;
      e(m) = 0.0;

    end

  end

  for ii = 2 : n

     i = ii - 1;
     k = i;
     p = d(i);

     for j = ii : n
       if ( d(j) < p )
         k = j;
         p = d(j);
       end
     end

     if ( k ~= i )
       d(k) = d(i);
       d(i) = p;
       p = z(i);
       z(i) = z(k);
       z(k) = p;
     end

  end

  return
end
function parchk ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% PARCHK checks parameters ALPHA and BETA for classical weight functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, integer M, the order of the highest moment to
%    be calculated.  This value is only needed when KIND = 8.
%
%    Input, real ALPHA, BETA, the parameters, if required
%    by the value of KIND.
%
  if ( kind <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  KIND <= 0.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check ALPHA for Gegenbauer, Jacobi, Laguerre, Hermite, Exponential.
%
  if ( 3 <= kind && alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  3 <= KIND and ALPHA <= -1.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check BETA for Jacobi.
%
  if ( kind == 4 && beta <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  KIND == 4 and BETA <= -1.0.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check ALPHA and BETA for rational.
%
  if ( kind == 8 )
    tmp = alpha + beta + m + 1.0;
    if ( 0.0 <= tmp || tmp <= beta )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PARCHK - Fatal error!\n' );
      fprintf ( 1, '  KIND == 8 but condition on ALPHA and BETA fails.\n' );
      error ( 'PARCHK - Fatal error!' );
    end
  end

  return
end
function value = r8_huge ( )

%*****************************************************************************80
%
%% R8_HUGE returns a "huge" real number.
%
%  Discussion:
%
%    The value returned by this function is NOT required to be the
%    maximum representable R8.  This value varies from machine to machine,
%    from compiler to compiler, and may cause problems when being printed.
%    We simply want a "very large" but non-infinite number.
%
%    MATLAB provides a built-in symbolic constant "inf" that can be used
%    if a huge number is really what you want!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, a huge number.
%
  value = 1.0E+30;

  return
end
function value = r8_sign ( x )

%*****************************************************************************80
%
%% R8_SIGN returns the sign of an R8.
%
%  Discussion:
%
%    The value is +1 if the number is positive or zero, and it is -1 otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose sign is desired.
%
%    Output, real VALUE, the sign of X.
%
  if ( 0 <= x )
    value = +1.0;
  else
    value = -1.0;
  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6e', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16e', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function rule_write ( order, filename, x, w, r )

%*****************************************************************************80
%
%% RULE_WRITE writes a quadrature rule to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, string FILENAME, specifies the output files.
%    write files 'filename_w.txt', 'filename_x.txt', 'filename_r.txt' defining
%    weights, abscissas, and region.
%
%    Input, real X(ORDER), the abscissas.
%
%    Input, real W(ORDER), the weights.
%
%    Input, real R(2), the region.
%
  filename_x = strcat ( filename, '_x.txt' );
  filename_w = strcat ( filename, '_w.txt' );
  filename_r = strcat ( filename, '_r.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1,'  Creating quadrature files.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  "Root" file name is   "%s".\n', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weight file will be   "%s".\n', filename_w );
  fprintf ( 1, '  Abscissa file will be "%s".\n', filename_x );
  fprintf ( 1, '  Region file will be   "%s".\n', filename_r );

  r8mat_write ( filename_w, 1, order, w' );
  r8mat_write ( filename_x, 1, order, x' );
  r8mat_write ( filename_r, 1, 2,     r' );

  return
end
function [ t, wts ] = scqf ( nt, t, mlt, wts, nwts, ndx, kind, alpha, ...
  beta, a, b )

%*****************************************************************************80
%
%% SCQF scales a quadrature formula to a nonstandard interval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real T(NT), the original knots.
%
%    Input, integer MLT(NT), the multiplicity of the knots.
%
%    Input, real WTS(NWTS), the weights.
%
%    Input, integer NWTS, the number of weights.
%
%    Input, integer NDX(NT), used to index the array WTS.
%    For more details see the comments in CAWIQ.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real T(NT), the scaled knots.
%
%    Output, real WTS(NWTS), the scaled weights.
%
  temp = eps;

  parchk ( kind, 1, alpha, beta )

  if ( kind == 1 )

    al = 0.0;
    be = 0.0;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 2 )

    al = -0.5;
    be = -0.5;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 3 )

    al = alpha;
    be = alpha;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 4 )

    al = alpha;
    be = beta;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 5 )

    if ( b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = 1.0 / b;
    al = alpha;
    be = 0.0;

  elseif ( kind == 6 )

    if ( b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = 1.0 / sqrt ( b );
    al = alpha;
    be = 0.0;

  elseif ( kind == 7 )

    al = alpha;
    be = 0.0;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 8 )

    if ( a + b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  A + B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = a + b;
    al = alpha;
    be = beta;

  elseif ( kind == 9 )

    al = 0.5;
    be = 0.5;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  end

  p = slp^( al + be + 1.0 );

  for k = 1 : nt

    t(k) = shft + slp * t(k);
    l = abs ( ndx(k) );

    if ( l ~= 0 )
      tmp = p;
      for i = l : l + mlt(k) - 1
        wts(i) = wts(i) * tmp;
        tmp = tmp * slp;
      end
    end

  end

  return
end
function [ t, wts ] = sgqf ( nt, aj, bj, zemu )

%*****************************************************************************80
%
%% SGQF computes knots and weights of a Gauss Quadrature formula.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with simple knots from the Jacobi matrix and the zero-th
%    moment of the weight function, using the Golub-Welsch technique.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real AJ(NT), the diagonal of the Jacobi matrix.
%
%    Input, real BJ(NT), the subdiagonal of the Jacobi
%    matrix, in entries 1 through NT-1.  On output, BJ has been overwritten.
%
%    Input, real ZEMU, the zero-th moment of the weight function.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Exit if the zero-th moment is not positive.
%
  if ( zemu <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SGQF - Fatal error!\n' );
    fprintf ( 1, '  ZEMU <= 0.\n' );
    error ( 'SGQF - Fatal error!' );
  end
%
%  Set up vectors for IMTQLX.
%
  wts = zeros ( nt, 1 );

  wts(1) = sqrt ( zemu );
  wts(2:nt) = 0.0;
%
%  Diagonalize the Jacobi matrix.
%
  [ t, wts ] = imtqlx ( nt, aj, bj, wts );

  wts(1:nt) = wts(1:nt).^2;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function [data,units] = compute_min_wing_length(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  data{i} = min(trx(fly).wing_lengthl_mm,trx(fly).wing_lengthr_mm);

end
units = parseunits('mm');


text:
% PCA Assignment
% clc
% clear all;
% close all;
function [recp,rectime,TDS,DS]= PCA_NEW(TrNum);

TsNum=10-TrNum; % Remained Number of Images for training, out of 400

[Tr,irow,icol,DS]=CDTr(TrNum);
[Ts,irow,icol,TDS]=CDTs(TrNum);

Tr=double(Tr);
Ts=double(Ts);
M=mean(Tr,2);
A=Tr-repmat(M,1,size(Tr,2));

 L=cov(A); % surogate Matrix of mean subtracted
 [U D]=eig(L);% Eigen Vectors of Surogate Matrix

 L_eig_vec = [];
 for i = 1 : size(U,2)
        if( D(i,i)>4500 )
         L_eig_vec = [L_eig_vec U(:,i)];
        end
 end

Eigenfaces = A * L_eig_vec;% Eigenvectors of Mean Subtracted Matrix

%%%%%%%%%%%%%%%%%%%%
ProjectedImages = [];
Train_Number = size(Tr,2);
for i = 1 : Train_Number
    temp = Eigenfaces'*A(:,i); % Projection of centered images into facespace
    ProjectedImages = [ProjectedImages temp];
end
ProjectedTestImages=[];
rec=[];

for j=1:size(Ts,2)
Difference = Ts(:,j)-M; % Centered test image
ProjectedTestImages = [ProjectedTestImages Eigenfaces'*Difference]; % Test image feature vector
end
%%%%%%%%%%%%%%%%%%%%%%%% Calculating Euclidean distances
% Euclidean distances between the projected test image and the projection
% of all centered training images are calculated. Test image is
% supposed to have minimum distance with its corresponding image in the
% training database.
%tic

tic
for j=1:size(Ts,2)
Euc_dist = [];
for i = 1 : Train_Number
    %q = ProjectedImages(:,i);
    temp = ( norm( ProjectedTestImages(:,j) - ProjectedImages(:,i) ) )^2;
    Euc_dist = [Euc_dist temp];
end
%toc
%stem(Euc_dist);

[Euc_dist_min , Recognized_index] = min(Euc_dist);
rec=[rec Recognized_index];
end
rectime=toc;
recd=ceil(rec/TrNum);
%%%

 outd=[];
 for i = 1:DS/TrNum
     for j=1:TsNum
         outd=[outd i];
     end
 end

recp=1-(size(find((recd-outd)>0),2)/size(Ts,2));
sprintf('Recognition Percentage %1.2f%%  and recognition time is %1.3f sec for %d out of %d number of images',recp*100,rectime,TDS,DS)
end


text:
function R = get(DS, props)
%GET gets the properties of the dataset
%
% $ Syntax $
%   - R = get(DS, propname)
%   - R = get(DS, propname cell array)
%
% $ Description $
%   - R = get(DS, propname) gets a property value with its name specified
%     by propname
%
%   - R = get(DS, propname cell array) gets a set of properties with the
%     names given in a cell array. Consequently, the property values are
%     returns through a cell array of the same size.
%
% $ Remarks $
%   - The
% $ History $
%   - Created by Dahua Lin on Jul 23rd, 2006
%

%% Skeleton

if ischar(props)
    S = struct(DS);
    R = getprop(S, props);

elseif iscell(props)
    S = struct(DS);
    R = cell(size(props));
    n = numel(R);
    for i = 1 : n
        if ~ischar(props{i})
            error('dsdml:invalidarg', ...
                'The property name should be a char string');
        end
        R{i} = getprop(S, props{i});
    end

else
    error('dsdml:invalidarg', ...
        'The props should be a property name or a cell of property names');
end

%% Property getting function

function R = getprop(S, propname)

if isfield(S, propname)
    R = S.(propname);
elseif ~isempty(S.attribs) && isfield(S.attribs, propname)
    R = S.attribs.(propname);
else
    switch(propname)
        case 'numsamples'
            R = getnumsamples(S);
        case 'numgroups'
            R = getnumgroups(S);
        case 'numunits'
            R = length(S.units);

        case 'samples'
            R = getsamples(S);
        case 'groups'
            R = getgroups(S);
        case 'class_ids'
            R = getclass_ids(S);
        case 'filenames'
            R = getfilenames(S);
        case 'nums'
            R = getnums(S);

        case 'numspergrp'
            R = getnumspergrp(S);
        case 'grpnums'
            R = getgrpnums(S);
        case 'grpclass_ids'
            R = getgrpclass_ids(S);

        otherwise
            error('dsdml:invalidarg', ...
                'Invalid property name of the dataset: %s', propname);
    end
end

%% Specific property getting functions

function R = getnumsamples(S)

switch S.unittype
    case 'Sample'
        R = length(S.units);
    case 'SampleGroup'
        R = 0;
        ng = length(S.units);
        for i = 1 : ng
            R = R + length(S.units(i).samples);
        end
    otherwise
        invalid_utype(S);
end

function R = getnumgroups(S)

switch S.unittype
    case {'Sample', 'SampleGroup'}
        R = length(S.units);
    otherwise
        invalid_utype(S);
end

function R = getsamples(S)

switch S.unittype
    case 'Sample'
        R = S.units;
    case 'SampleGroup'
        R = vertcat(S.units.samples);
    otherwise
        invalid_utype(S);
end

function R = getgroups(S)

switch S.unittype
    case 'Sample'
        n = length(S.units);
        U = S.units;
        R = struct(...
            'class_id', cell(n, 1), ...
            'samples', cell(n, 1), ...
            'attribs', cell(n, 1));
        for i = 1 : n
            R(i).class_id = U(i).class_id;
            R(i).samples = U(i);
        end
    case 'SampleGroup'
        R = S.units;
    otherwise
        invalid_utype(S);
end

function R = getclass_ids(S)

switch(S.unittype)
    case 'Sample'
        n = length(S.units);
        R = zeros(n, 1);
        for i = 1 : n
            R(i) = S.units(i).class_id;
        end
    case 'SampleGroup'
        n = getnumsamples(S);
        R = zeros(n, 1);
        c = 0;
        ng = length(S.units);
        for i = 1 : ng
            ns = length(S.units(i).samples);
            R(c+1:c+ns) = S.units(i).class_id;
            c = c + ns;
        end
    otherwise
        invalid_utype(S);
end

function R = getfilenames(S)

switch(S.unittype)
    case 'Sample'
        n = length(S.units);
        R = cell(n, 1);
        for i = 1 : n
            R{i} = S.units(i).filename;
        end
    case 'SampleGroup'
        n = getnumsamples(S);
        R = cell(n, 1);
        c = 0;
        ng = length(S.units);
        for i = 1 : ng
            ns = length(S.units(i).samples);
            [R{c+1:c+ns}] = deal(S.units(i).samples.filename);
            c = c + ns;
        end
    otherwise
        invalid_utype(S);
end

function R = getnums(S)

R = countnums(getclass_ids(S));

function R = getnumspergrp(S)

switch (S.unittype)
    case 'Sample'
        n = length(S.units);
        R = ones(n, 1);

    case 'SampleGroup'
        n = length(S.units);
        R = zeros(n, 1);
        for i = 1 : n
            R(i) = length(S.units(i).samples);
        end

    otherwise
        invalid_utype(S);
end

function R = getgrpnums(S)

R = countnums(getgrpclass_ids(S));

function R = getgrpclass_ids(S)

switch (S.unittype)
    case {'Sample', 'SampleGroup'}
        n = length(S.units);
        R = zeros(n, 1);
        for i = 1 : n
            R(i) = S.units(i).class_id;
        end

    otherwise
        invalid_utype(S);
end

%% Auxiliary functions

function invalid_utype(S)

error('dsdml:invalidunittype', ...
            'Invalid unit type: %s', S.unittype);

function nums = countnums(A)

A = A(:);
difs = diff(A);
n = length(A);
nums = [0; find(difs ~= 0); n];
nums = diff(nums);

author: lmthang


text:
function poly2 = clipPolygonHP(poly, line)
%CLIPPOLYGONHP Clip a polygon with a Half-plane defined by a directed line
%
%   POLY2 = clipPolygonHP(POLY, LINE)
%   POLY is a [Nx2] array of points, and LINE is given as [x0 y0 dx dy].
%   The result POLY2 is also an array of points, sometimes smaller than
%   poly, and that can be [0x2] (empty polygon).
%
%   See also:
%   polygons2d, clipPolygon
%
% ---------
% author : David Legland
% created the 31/07/2005.
% Copyright 2010 INRA - Cepia Software Platform.
%

%   HISTORY
%   15/08/2005 add test to avoid empty polygons
%   13/06/2007 deprecate
%   10/10/2008 'reprecate'

% avoid to process empty polygons
if size(poly, 1)<3
    poly2 = zeros([0 2]);
    return;
end

% ensure the last point is the same as the first one
if sum(poly(end, :)==poly(1,:))~=2
    poly = [poly; poly(1,:)];
end

N = size(poly, 1);
edges = [poly([N 1:N-1], :) poly];

b = isLeftOriented(poly, line);

% case of totally clipped polygon
if sum(b)==0
    poly2 = zeros(0, 2);
    return;
end

poly2 = zeros(0, 2);

i=1;
while i<=N

    if isLeftOriented(poly(i,:), line)
        % keep all points located on the right side of line
        poly2 = [poly2; poly(i,:)]; %#ok<AGROW>
    else
        % compute of preceeding edge with line
        if i>1
            poly2 = [poly2; intersectLineEdge(line, edges(i, :))]; %#ok<AGROW>
        end

        % go to the next point on the left side
        i=i+1;
        while i<=N

            % find the next point on the right side
            if isLeftOriented(poly(i,:), line)
                % add intersection of previous edge
                poly2 = [poly2; intersectLineEdge(line, edges(i, :))]; %#ok<AGROW>

                % add current point
                poly2 = [poly2; poly(i,:)]; %#ok<AGROW>

                % exit the second loop
                break;
            end

            i=i+1;
        end
    end

    i=i+1;
end

% remove last point if it is the same as the first one
if sum(poly2(end, :)==poly(1,:))==2
    poly2 = poly2(1:end-1, :);
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%ASSEMBLECONSTR Assemble integral constraints.
%
%   [ C, D ] = ASSEMBLECONSTR( PROB, ICUB ) Assemble integral constraints.
%
%   Assemble subdomain/boundary integral constraints specified in the
%   PROB.CONSTR field. Each constr entry have TYPE (intsubd/intbdr),
%   DVAR, INDEX (indices to subdomains/boundaries), and EXPR
%   (resulting integrand expression, must evaluate to a numeric
%   scalar) entries, for example
%
%       const(1).type = 'intsubd';   % Subdomain integral constraint.
%       const(1).dvar = 'u';         % Applied to dependent variable u.
%       const(1).index = [1,2];      % Applied to domains 1 and 2.
%       const(1).expr  = 2.3;        % Resulting integrand.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function fixed = p09_fixed_points ( m, fixed_num )

%*****************************************************************************80
%
%% P09_FIXED_POINTS returns the fixed points in problem 09.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer FIXED_NUM, the number of fixed points.
%
%    Output, real FIXED(M,FIXED_NUM), the fixed points.
%
  center1 = [ 0.50, 0.50 ];
  center2 = [ 0.25, 0.75 ];
  center3 = [ 0.60, 0.40 ];

  r1 = 0.5;
  r2 = 0.1;
  r3 = 0.1;

  fixed[1:2,1:4] = [ ...
    center1(1) - r1, center1(2) - r1; ...
    center1(1) + r1, center1(2) - r1; ...
    center1(1) + r1, center1(2) + r1; ...
    center1(1) - r1, center1(2) + r1 ]';

  for j = 1 : 6
    angle = ( ( j - 1 ) * 2 ) * pi / 6.0;
    fixed(1,4+j) = center2(1) + r2 * cos ( angle );
    fixed(2,4+j) = center2(2) + r2 * sin ( angle );
  end

  for j = 1 : 6
    angle = ( ( j - 1 ) * 2 ) * pi / 6.0;
    fixed(1,10+j) = center3(1) + r3 * cos ( angle );
    fixed(2,10+j) = center3(2) + r3 * sin ( angle );
  end

  return
end


text:
function[result] = strStartsWith(str, prefix)
% [result] = strStartsWith(str, prefix)
%
% Check whether a given string starts with a given prefix.
%
% Copyright by Holger Caesar, 2014

result = false;
assert(~isempty(prefix));

if iscell(str),
    result = cellfun(@(x) strStartsWith(x, prefix), str);
else
    if numel(str) >= numel(prefix) && strcmp(str(1:numel(prefix)), prefix),
        result = true;
    end;
end;

author: nightrome

repo: matconvnet-calvin


text:
function [Y,MX] = logfsgram(X, N, SR, WIN, NOV, FMIN, BPO)
% [Y,MX] = logfsgram(X, N, SR, WIN, NOV, FMIN, BPO)
%    Calculate a log-frequency spectrogram
%    X is input signal; N is parent FFT window; SR is the source samplerate.
%    WIN is actual window length within FFT, NOV is number of overlapping
%    points between successive windows.
%    Optional FMIN is the lowest frequency to display (80Hz);
%    BPO is the number of bins per octave (12).
%    MX returns the nlogbin x nfftbin mapping matrix;
%    sqrt(MX'*(Y.^2)) is an approximation to the original FFT
%    spectrogram that Y is based on, suitably blurred by going
%    through the log-F domain.
% 2004-03-30 dpwe@ee.columbia.edu $Header: /homes/dpwe/matlab/columbiafns/RCS/logfsgram.m,v 1.3 2004/04/01 22:39:02 dpwe Exp $

if nargin < 2
  N = 1024;
end
if nargin < 3
  SR = 8000;
end
if nargin < 4
  WIN = [];
end
if nargin < 5
  NOV = [];
end
if nargin < 6
  FMIN = 80;
end
if nargin < 7
  BPO = 12;
end

if isempty(WIN)
  WIN = N;
end
if isempty(NOV)
  NOV = WIN/2;
end

% Calculate underlying STFT
XX = specgram(X,N,SR,WIN,NOV);

% Construct mapping matrix

% Ratio between adjacent frequencies in log-f axis
fratio = 2^(1/BPO);

% How many bins in log-f axis
nbins = floor( log((SR/2)/FMIN) / log(fratio) );

% Freqs corresponding to each bin in FFT
fftfrqs = [0:(N/2)]*(SR/N);
nfftbins = N/2+1;

% Freqs corresponding to each bin in log F output
logffrqs = FMIN * exp(log(2)*[0:(nbins-1)]/BPO);

% Bandwidths of each bin in log F
logfbws = logffrqs * (fratio - 1);

% .. but bandwidth cannot be less than FFT binwidth
logfbws = max(logfbws, SR/N);

% Controls how much overlap there is between adjacent bands
ovfctr = 0.5475;   % Adjusted by hand to make sum(mx'*mx) close to 1.0

% Weighting matrix mapping energy in FFT bins to logF bins
% is a set of Gaussian profiles depending on the difference in
% frequencies, scaled by the bandwidth of that bin
freqdiff = ( repmat(logffrqs',1,nfftbins) - repmat(fftfrqs,nbins,1) )./repmat(ovfctr*logfbws',1,nfftbins);
mx = exp( -0.5*freqdiff.^2 );
% Normalize rows by sqrt(E), so multiplying by mx' gets approx orig spec back
mx = mx ./ repmat(sqrt(2*sum(mx.^2,2)), 1, nfftbins);

% Perform mapping in magnitude-squared (energy) domain
y = sqrt( mx * (abs(XX).^2) );

% so, we lost phase information...

if nargout < 1
  imagesc([0 length(X)/SR],[1 nbins],20*log10(y));
  axis xy
  xlabel('Time');
  ylabel('Frequency');
  yt = get(gca,'YTick');
  for i = 1:length(yt)
    ytl{i} = sprintf('%.0f',logffrqs(yt(i)));
  end
  set(gca,'YTickLabel',ytl);
else
  Y = y;
  MX = mx;
end

author: posenhuang

repo: deeplearningsourceseparation


text:
function [w, infos] = cd_lasso_elasticnet(problem, options)
% Coordinate descent algorithm for LASSO problem.
%
% Inputs:
%       problem     function (cost/grad/hess)
%       options     options
% Output:
%       w           solution of w
%       infos       information
%
% This file is part of GDLibrary and SGDLibrary.
%
% Created by H.Kasai on Apr. 18, 2017

    % set dimensions and samples
    d = problem.dim();
    n = problem.samples();
    A = problem.A();

    % extract options
    if ~isfield(options, 'tol_optgap')
        tol_optgap = 1.0e-12;
    else
        tol_optgap = options.tol_optgap;
    end

    if ~isfield(options, 'tol_gnorm')
        tol_gnorm = 1.0e-12;
    else
        tol_gnorm = options.tol_gnorm;
    end

    if ~isfield(options, 'max_iter')
        max_iter = 100;
    else
        max_iter = options.max_iter;
    end

    if ~isfield(options, 'verbose')
        verbose = false;
    else
        verbose = options.verbose;
    end

    if ~isfield(options, 'w_init')
        w = randn(d,1);
    else
        w = options.w_init;
    end

    if ~isfield(options, 'f_opt')
        f_opt = -Inf;
    else
        f_opt = options.f_opt;
    end

    if ~isfield(options, 'store_w')
        store_w = false;
    else
        store_w = options.store_w;
    end

    if ~isfield(options, 'sub_mode')
        sub_mode = 'lasso';
    else
        sub_mode = options.sub_mode;
    end

    % initialise
    iter = 0;
    if strcmp(sub_mode, 'lasso')
        AtA = problem.AtA();
        squred_norm_col = diag(AtA);
    else
        AtA_l2 = problem.AtA_l2();
        squred_norm_col = diag(AtA_l2);
    end
    prox_th = ones(d, 1)./squred_norm_col;

    % store first infos
    clear infos;
    infos.iter = iter;
    infos.time = 0;
    infos.grad_calc_count = 0;
    f_val = problem.cost(w);
    infos.cost = f_val;
    optgap = f_val - f_opt;
    infos.optgap = optgap;
    grad = problem.full_grad(w);
    gnorm = norm(grad);
    infos.gnorm = gnorm;
    if ismethod(problem, 'reg')
        infos.reg = problem.reg(w);
    end
    if store_w
        infos.w = w;
    end

    % set start time
    start_time = tic();

    % print info
    if verbose
        fprintf('CD (%s): Iter = %03d, cost = %.24e, gnorm = %.4e, optgap = %.4e\n', sub_mode, iter, f_val, gnorm, optgap);
    end

    % main loop
    while (optgap > tol_optgap) && (gnorm > tol_gnorm) && (iter < max_iter)

        % update i-th coordinate
        if strcmp(sub_mode, 'lasso')
            for i = 1:d
                w_except_i = w;
                w_except_i(i) = 0;
                residual = problem.residual(w_except_i);

                snc = squred_norm_col(i);
                w(i) = problem.prox(A(:, i)'*residual/snc, prox_th(i));
            end
        else
            for i = 1:d
                w_except_i = w;
                w_except_i(i) = 0;
                residual = problem.residual(w_except_i, i);

                snc = squred_norm_col(i);
                w(i) = problem.prox(residual/snc, prox_th(i));
            end

        end

        % calculate gradient
        grad = problem.full_grad(w);

        % update iter
        iter = iter + 1;
        % calculate error
        f_val = problem.cost(w);
        optgap = f_val - f_opt;
        % calculate norm of gradient
        gnorm = norm(grad);

        % measure elapsed time
        elapsed_time = toc(start_time);

        % store infoa
        infos.iter = [infos.iter iter];
        infos.time = [infos.time elapsed_time];
        infos.grad_calc_count = [infos.grad_calc_count iter*n];
        infos.optgap = [infos.optgap optgap];
        infos.cost = [infos.cost f_val];
        infos.gnorm = [infos.gnorm gnorm];
        if ismethod(problem, 'reg')
            reg = problem.reg(w);
            infos.reg = [infos.reg reg];
        end
        if store_w
            infos.w = [infos.w w];
        end

        % print info
        if verbose
            fprintf('CD (%s): Iter = %03d, cost = %.24e, gnorm = %.4e, optgap = %.4e\n', sub_mode, iter, f_val, gnorm, optgap);
        end
    end

    if gnorm < tol_gnorm
        fprintf('Gradient norm tolerance reached: tol_gnorm = %g\n', tol_gnorm);
    elseif optgap < tol_optgap
        fprintf('Optimality gap tolerance reached: tol_optgap = %g\n', tol_optgap);
    elseif iter == max_iter
        fprintf('Max iter reached: max_iter = %g\n', max_iter);
    end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function error_frobenius = c8mat_is_symmetric ( m, n, a )

%*****************************************************************************80
%
%% C8MAT_IS_SYMMETRIC checks a complex matrix for symmetry.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Input, complex A(M,N), the matrix.
%
%    Output, real ERROR_FROBENIUS, measures the Frobenius norm of ( A - A' ).
%
  if ( m ~= n )
    error_frobenius = r8_huge ( );
    return
  end

  error_frobenius = sqrt ( sum ( sum ( ( abs ( a - a' ) ).^2 ) ) );

  return
end


text:
function objPC = loadPC(objICP, i)
% LOADPC Load point cloud to workspace.
% ------------------------------------------------------------------------------
% DESCRIPTION/NOTES
% With this method a point cloud which was previously added to the globalICP
% object, can be loaded to the workspace.
% ------------------------------------------------------------------------------
% INPUT
% [i]
% Index of point cloud to load.
% ------------------------------------------------------------------------------
% EXAMPLES
% Call 'help globalICP.globalICP' for a minimal working example, which also
% includes this method.
% ------------------------------------------------------------------------------
% philipp.glira@gmail.com
% ------------------------------------------------------------------------------

objPC = pointCloud(objICP.PC{i});

end

author: pglira


text:
function y = logMvGamma( x, D)
% Compute logarithm multivariate Gamma function.
%   operates on each entry of input matrix/vector x
%INPUT
%   x : any scalar/vector/matrix
%   D : integer dimension of the mv Gamma function
%MATH DETAILS -------------------------------------------------------
% Gamma_D(x) = pi^(D(D-1)/4) prod_(j=1)^p Gamma(x+(1-j)/2)
% log Gamma_D(x) = D(D-1)/4 log pi + sum_(j=1)^p log Gamma(x+(1-j)/2)
% Credit: Michael Chen (sth4nth@gmail.com).

s = size(x);
ds = (1-(1:D)')/2;

% Force input x to be a row vector
X = reshape(x,1,prod(s));

% X(dd, dim) = input(dim) - ds(dd)
X = bsxfun(@plus, X, ds   );

y = D*(D-1)/4*log(pi) + sum(gammaln(X),1);

% Force output back to size of original input
y = reshape(y,s);


text:
classdef prtClassConstant < prtClass

    properties (SetAccess=private)

        name = 'Constant' % Fisher Linear Discriminant
        nameAbbreviation = 'Constant'            % FLD
        isNativeMary = true;  % False
    end

    properties
        constant
    end

    methods
        % Allow for string, value pairs
        function self = prtClassConstant(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end

    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,ds)
            self.constant = mean(ds.getY);
        end

        function ds = runAction(self,ds)
            ds.X = repmat(self.constant,ds.nObservations,1);
        end

    end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [y,sy,st] = spm_mci_sens_init (R,P,M,U)
% Compute sensitivity to initial state
% FORMAT [y,sy,st] = spm_mci_sens_init (R,P,M,U)
%
% R         Initial state
% P         Parameters
% M         Model structure
% U         Inputs  [Nin x N]
%
% y         Outputs     [N x Nout]
% sy        Output Sensitivity, dy/dP [N x Nout x Nparams]
% st        Status flag (0 for OK, -1 for problem)
%           ... evaluated at the N time points in M.t
%
% M.f       Flow function dx/dt=f(x,u,P,M)
% M.g       Observation function y=g(x,u,P,M)
%
% This function uses Matlab's ODE suite
%
% B. Sengupta, K. Friston and W. Penny (2014) Efficient Gradient
% Computation for Dynamical Models. Neuroimage,98, 521-527.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny and Biswa Sengupta
% $Id: spm_mci_sens_init.m 6697 2016-01-27 14:57:28Z spm $

y=[];sy=[];x=[];sx=[];
st=0;

% Tolerances for ode15s
try, tol.rel=M.reltol; catch, tol.rel=1e-2; end
try, tol.abs=M.abstol; catch, tol.abs=1e-4; end

if isempty(U)
    U=zeros(1,M.N);
end

init_t=0;
final_t=M.T;
init_state=R;

% parameters for the integrator
options = odeset('AbsTol',tol.abs,'RelTol',tol.rel);

% Allocate matrices
x=zeros(M.N,M.n);
y=zeros(M.N,M.l);
sx=zeros(M.N,M.n,M.n);
sy=zeros(M.N,M.l,M.n);

% initialise state-sensitivites
Nx=length(M.x0);
Np=length(P);
S0=eye(M.n);
init_sens = S0(:);

% Use Klopfenstein-Shampine integrator from Matlab's ODE suite
[T,V] = ode15s(@(t,v) int_states_init(t,v,U,P,M),[init_t final_t], [init_state; init_sens], options);

% T     Times
% V     States and sensitivities

% Extract states
xm=V(:,1:Nx);

% Interpolate states to times M.t
x=interp1q(T,xm,M.t);

% Interpolate state sensitivities to times M.t
sxm=interp1q(T,V(:,Nx+1:end),M.t);
%sx=reshape(sxm,length(M.t),Nx,Np);
% Updated Sep 13 2014
sx=reshape(sxm,length(M.t),Nx,Nx);

% When computing output sensitivities, assume dydx=L, ie not a
% function of x. Generalise later
[tmp,L]=feval(M.g,M.x0,U(:,1),P,M);

% Compute output and output sensitivity
for n=1:M.N,
    yout = feval (M.g,x(n,:)',U(:,n),P,M);
    y(n,:) = yout';
    sy(n,:,:)=L*squeeze(sx(n,:,:));
end

end

%--------------------------------------------------------------------------
function DvDt = int_states_init (t,v,U,P,M)
% Integrate states and sensitivities to initial conditions

Nx=length(M.x0);
Nt=Nx+Nx*Nx;

state_ind=1:Nx;
sens_ind=Nx+1:Nt;

% Find nearest time point for which we have pre-computed input
if isempty(U)
    ut=[];
else
    [tmp,ind]=min(abs(t-M.t));
    ut=U(:,ind);
end

% initialise state matrix
DvDt = zeros(Nt,1);

% Flow
DvDt(1:Nx)=feval(M.f,v(state_ind),ut,P,M);

if isfield(M,'dfdx')
    Fx = feval(M.dfdx,v(1:Nx),ut,P,M);
else
    Fx = spm_diff(M.f,v(1:Nx),ut,P,M,1);
end

Sx_old = v(sens_ind);
Sx_old = reshape(Sx_old,Nx,Nx);

% Compute change in sensitivity from old - note absence of Fp term
Sx = Fx*Sx_old;

% Sensitivities
DvDt(sens_ind)=Sx(:);

end

author: spm


text:
classdef AbstractMesh < handle

    properties (GetAccess = public, SetAccess = protected)
        coord
        connec

        nelem
        ndim

        type

        coordElem
        interpolation

        edges
    end

    properties (Access = private)
       xFE
    end

    methods (Access = public)

       function xV = computeBaricenter(obj)
            xV = obj.xFE.computeValueInCenterElement();
       end

       function xGauss = computeXgauss(obj,xV)
            xGauss = obj.xFE.interpolateFunction(xV);
       end

       function dvolume = computeDvolume(obj,quad)
            s.mesh = obj;
            g = Geometry.create(s);
            g.computeGeometry(quad,obj.interpolation);
            dvolume = g.dvolu;
            dvolume = dvolume';
       end

       function q = computeElementQuality(obj)
            quad = Quadrature.set(obj.type);
            quad.computeQuadrature('CONSTANT');
            volume = obj.computeDvolume(quad);
            L(1,:) = obj.computeSquarePerimeter();
            q = 4*sqrt(3)*volume./L;
       end

       function v = computeVolume(obj)
            quad = Quadrature.set(obj.type);
            quad.computeQuadrature('CONSTANT');
            v = obj.computeDvolume(quad);
            v = sum(v(:));
       end

        function computeEdges(obj)
            s.nodesByElem = obj.connec;
            edge = EdgesConnectivitiesComputer(s);
            edge.compute();
            obj.edges = edge;
        end

    end

    methods (Access = protected)

        function createInterpolation(obj)
            obj.interpolation = Interpolation.create(obj,'LINEAR');
        end

        function computeElementCoordinates(obj)
            obj.computeCoordFEfunction();
            obj.coordElem = obj.xFE.fValues;
        end

    end

    methods (Access = private)

        function computeCoordFEfunction(obj)
            s.mesh    = obj;
            s.fValues = obj.coord;
            obj.xFE = P1Function(s);
        end

        function L = computeSquarePerimeter(obj)
            obj.computeEdges();
            nElem = size(obj.connec,1);
            L = zeros(nElem,1);
            for iedge = 1:obj.edges.nEdgeByElem
                edge = obj.edges.edgesInElem(:,iedge);
                nodesEdge = obj.edges.nodesInEdges(edge,:);
                for idim = 1:obj.ndim
                    xA = obj.coord(nodesEdge(:,1),idim);
                    xB = obj.coord(nodesEdge(:,2),idim);
                    L = L + (xA - xB).^2;
                end
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [params, names] = dnetExtractParam(model)

% DNETEXTRACTPARAM Extract weights and biases from an DNET.
% FORMAT
% DESC returns a vector of all the parameters from a density network.
% ARG model : the model from which we wish to extract the weights
% and biases.
% RETURN params : vector of all the weights and biases returned by
% the model. The structure is governed by dnetpak.
% RETURN names : optional additional returned cell array of the
% names of the parameters.
%
% SEEALSO : dnetCreate, dnetExpandParam, modelExtractParam
%
% COPYRIGHT : Neil D. Lawrence, 2008

% MLTOOLS

if nargout > 1
  [params, names] = modelExtractParam(model.mapping);
  names{end+1} = 'Output Precision';
else
  params = modelExtractParam(model.mapping);
end
func = str2func([model.betaTransform 'Transform']);
params(end+1) = func(model.beta, 'xtoa');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ x, seed ] = maxwell_sample ( a, seed )

%*****************************************************************************80
%
%% MAXWELL_SAMPLE samples the Maxwell PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the parameter of the PDF.
%    0 < A.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  a2 = 3.0;
  [ x, seed ] = chi_square_sample ( a2, seed );

  x = a * sqrt ( x );

  return
end


text:
function [Y] = mci_ramsay_gen (P,M,U)
% Generate data from Ramsay model
% FORMAT [Y] = mci_ramsay_gen (P,M,U)
%
% P         Parameters
% M         Model structure
% U         Inputs
%
% Y         Data
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_ramsay_gen.m 6548 2015-09-11 12:39:47Z will $

[G,x] = spm_mci_fwd (P,M,U);
[N,l] = size(G);
e = randn(N,l)*sqrt(M.Ce);
Y = G + e;


text: PitchTraking()

author: LimingShi


text:
function varargout = fillPolygon3d(varargin)
%FILLPOLYGON3D Fill a 3D polygon specified by a list of vertex coords.
%
%   fillPolygon3d(COORD, COLOR)
%   packs coordinates in a single [N*3] array.
%   COORD can also be a cell array of polygon, in this case each polygon is
%   drawn using the same color.
%
%   fillPolygon3d(PX, PY, PZ, COLOR)
%   specifies coordinates in separate numeric vectors (either row or
%   columns)
%
%   fillPolygon3d(..., PARAM, VALUE)
%   allows to specify some drawing parameter/value pairs as for the plot
%   function.
%
%   H = fillPolygon3d(...)
%   also returns a handle to the list of created patch objects.
%
%   Example
%     t = linspace(0, 2*pi, 100)';
%     xt = 10 * cos(t);
%     yt = 5 * sin(t);
%     zt = zeros(100,1);
%     poly = [xt yt zt];
%     figure; hold on; axis equal; fillPolygon3d(poly, 'c');
%     drawPolygon3d(poly, 'linewidth', 2, 'color', 'k');
%
%   See also
%   polygons3d, drawPolygon3d, drawPolyline3d
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2007-01-05
% Copyright 2007-2022 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas

% Check if axes handle is specified
if isAxisHandle(varargin{1})
    hAx = varargin{1};
    varargin(1) = [];
else
    hAx = gca;
end

% check case we want to draw several curves, stored in a cell array
var1 = varargin{1};
if iscell(var1)
    hold on;
    h = [];
    for i = 1:length(var1(:))
        h = [h; fillPolygon3d(hAx, var1{i}, varargin{2:end})]; %#ok<AGROW>
    end
    if nargout>0
        varargout{1}=h;
    end
    return;
end

% extract vertex coordinates
if min(size(var1)) == 1
    % if first argument is a vector (either row or column), then assumes
    % first argument contains x coords, second argument contains y coords
    % and third one the z coords
    px = var1;
    if length(varargin) < 3
        error('geom3d:fillPolygon3d:Wrong number of arguments in fillPolygon3d');
    end
    py = varargin{2};
    pz = varargin{3};
    varargin = varargin(4:end);
else
    % first argument contains all three coordinates
    px = var1(:, 1);
    py = var1(:, 2);
    pz = var1(:, 3);
    varargin = varargin(2:end);
end

% extract color information
if isempty(varargin)
    color = 'c';
else
    color = varargin{1};
    varargin = varargin(2:end);
end

% fill the polygon
h = fill3(hAx, px, py, pz, color, varargin{:});

if nargout>0
    varargout{1}=h;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function test_suite=test_randperm
% tests for cosmo_randperm
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_randperm_no_seed
    helper_test_deterministic(false);

function test_randperm_with_seed
    helper_test_deterministic(true);

function helper_test_deterministic(with_seed)
    if with_seed
        seed=ceil(rand()*1e5);
        args={'seed',seed};
    else
        args={};
    end

    f=@(varargin)cosmo_randperm(varargin{:},args{:});

    % single input
    n=randint();
    x1=f(n);
    assertEqual(size(x1),[1 n]);
    assertEqual(sort(x1),1:n);

    x1a=f(n);
    if with_seed
        assertEqual(x1,x1a);
    else
        assert(~isequal(x1,x1a));
    end

    % two inputs, select all
    x1a=f(n,n);

    if with_seed
        assertEqual(x1,x1a);
    else
        assert(~isequal(x1,x1a));
    end

    % two inputs, k<n
    k=randint();
    n=k+randint();
    assert(k<n);
    x2=f(n,k);
    assertEqual(size(x2),[1 k]);
    msk=bsxfun(@eq,(1:n)',x2);
    assertEqual(sum(msk(:)==1),k);
    assertEqual(sum(sum(msk==1,1),2),k); % each column has one 1
    assertEqual(sum(sum(msk==1,2),1),k); % each row has one 1

    if with_seed
        % deterministic
        x2a=f(n,k);
        assertEqual(x2,x2a);
    else
        found=false;
        for attempt=1:10
            x2a=f(n,k);
            if ~isequal(x2,x2a)
                found=true;
                break;
            end
        end

        if ~found
            error('Different calls do not lead to different outputs');
        end
    end

    x3=f(0);
    assertTrue(isempty(x3));

    x4=f(1);
    assertEqual(x4,1);

    x5=f(1,1);
    assertEqual(x5,1);

function test_randperm_different_seeds
    count=10;
    seed=0;
    result=cell(count,1);
    for k=1:count
        seed=seed+randint();

        result{k}=cosmo_randperm(1000,'seed',seed);
        for j=1:(k-1)
            assertFalse(isequal(result{j},result{k}));
        end
    end

function test_randperm_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                        cosmo_randperm(varargin{:}),'');
    illegal_args={-1,.5,[1 2],'foo',struct,cell(0)};
    narg=numel(illegal_args);

    % illegal first and/or second arguments
    for k=1:narg
        for second_arg=[false,true]
            for j=1:narg
                if second_arg
                    args=illegal_args([k j]);
                elseif j>1
                    continue;
                else
                    % single arg if j==1
                    args=illegal_args(k);
                end
                aet(args{:});
            end
        end
    end

    % illegal seed arguments
    for k=1:narg
        arg=illegal_args(k);
        aet(4,2,'seed',arg{:});
    end

    % missing seed
    aet(4,2,'seed');

    % double seed
    aet(4,2,'seed',1,'seed',1);

    % 3 numeric arguments
    aet(4,2,1);

    % second argument greater than first
    aet(2,4);
    aet(2,4,'seed',1);

    % unknown keyword
    aet(2,'foo',1);
    aet(2,'foo');

function x=randint(n)
    x=ceil(10+rand()*50);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function linplus_test30 ( )

%*****************************************************************************80
%
%% TEST30 tests R8GE_FA, R8GE_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST30\n' );
  fprintf ( 1, '  For a matrix in general storage,\n' );
  fprintf ( 1, '  R8GE_FA computes the LU factors,\n' );
  fprintf ( 1, '  R8GE_SL solves a factored system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  [ a, seed ] = r8ge_random ( n, n, seed );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  b = r8ge_mxv ( n, n, a, x );
%
%  Factor the matrix.
%
  [ a_lu, pivot, info ] = r8ge_fa ( n, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST30 - Fatal error!\n' );
    fprintf ( 1, '  R8GE_FA declares the matrix is singular!\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end
%
%  Solve the linear system.
%
  job = 0;
  x = r8ge_sl ( n, a_lu, pivot, b, job );

  r8vec_print ( n, x, '  Solution:' );
%
%  Set the desired solution.
%
  x(1:n) = 1.0E+00;
%
%  Compute the corresponding right hand side.
%
  job = 0;
  b = r8ge_ml ( n, a_lu, pivot, x, job );
%
%  Solve the system
%
  job = 0;
  x = r8ge_sl ( n, a_lu, pivot, b, job );

  r8vec_print ( n, x, '  Solution:' );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  job = 1;
  b = r8ge_ml ( n, a_lu, pivot, x, job );
%
%  Solve the system
%
  job = 1;
  x = r8ge_sl ( n, a_lu, pivot, b, job );

  r8vec_print ( n, x, '  Solution of transposed system:' );

  return
end

author: johannesgerer