text:
function db=lpccc2db(cc,np,nc)
%LPCCC2DBf Convert complex cepstrum to dB power spectrum DB=(CC,NP,NC)
% cc = complex cepstrum coefs
% np+2 = number of frequency values 0 to nyquist (default: np=size(cc,2))
% nc = number of cepstral coefs to use (default: nc=np);
% for high speed make np one less than a power of 2

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpccc2db.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,mc]=size(cc);
if nargin<2 np=mc+2; end
if nargin<3 nc=np; end
k=10/log(10);
if nc==mc
   db=k*rsfft([zeros(nf,1) cc].',2*np+2).';
else
   db=k*rsfft([zeros(nf,1) lpccc2cc(cc,nc)].',2*np+2).';
end
if nargout==0
   plot((0:np+1)/(2*np+2),db.');
   ylabel('dB');
end

author: decouples


text:
function sD = som_normalize(sD,method,comps)

%SOM_NORMALIZE (Re)normalize data or add new normalizations.
%
% sS = som_normalize(sS,[method],[comps])
%
%   sS = som_normalize(sD)
%   sS = som_normalize(sS,sNorm)
%    D = som_normalize(D,'var')
%   sS = som_normalize(sS,'histC',[1:3 10])
%
%  Input and output arguments ([]'s are optional):
%   sS                The data to which the normalization is applied.
%                     The modified and updated data is returned.
%            (struct) data or map struct
%            (matrix) data matrix (a matrix is also returned)
%   [method]          The normalization method(s) to add/use. If missing,
%                     or an empty variable ('') is given, the
%                     normalizations in sS are used.
%            (string) identifier for a normalization method to be added:
%                     'var', 'range', 'log', 'logistic', 'histD' or 'histC'.
%            (struct) Normalization struct, or an array of such.
%                     Alternatively, a map/data struct can be given
%                     in which case its '.comp_norm' field is used
%                     (see below).
%            (cell array) Of normalization structs. Typically, the
%                     '.comp_norm' field of a map/data struct. The
%                     length of the array must be equal to data dimension.
%            (cellstr array) norm and denorm operations in a cellstr array
%                     which are evaluated with EVAL command with variable
%                     name 'x' reserved for the variable.
%   [comps]  (vector) the components to which the normalization is
%                     applied, default is [1:dim] ie. all components
%
% For more help, try 'type som_normalize' or check out online documentation.
% See also SOM_DENORMALIZE, SOM_NORM_VARIABLE, SOM_INFO.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_normalize
%
% PURPOSE
%
% Add/apply/redo normalization on data structs/sets.
%
% SYNTAX
%
%  sS = som_normalize(sS)
%  sS = som_normalize(sS,method)
%   D = som_normalize(D,sNorm)
%  sS = som_normalize(sS,csNorm)
%  sS = som_normalize(...,comps)
%
% DESCRIPTION
%
% This function is used to (initialize and) add, redo and apply
% normalizations on data/map structs/sets. If a data/map struct is given,
% the specified normalizations are added to the '.comp_norm' field of the
% struct after ensuring that all normalizations specified therein have
% status 'done'. SOM_NORMALIZE actually uses function SOM_NORM_VARIABLE
% to handle the normalization operations, and only handles the data
% struct/set specific stuff itself.
%
% The different normalization methods are listed below. For more
% detailed descriptions, see SOM_NORM_VARIABLE.
%
%   method     description
%   'var'      Variance is normalized to one (linear operation).
%   'range'    Values are normalized between [0,1] (linear operation).
%   'log'      Natural logarithm is applied to the values:
%                xnew = log(x-m+1)
%              where m = min(x).
%   'logistic' Logistic or softmax trasformation which scales all
%              possible values between [0,1].
%   'histD'    Histogram equalization, values scaled between [0,1].
%   'histC'    Approximate histogram equalization with partially
%              linear operations. Values scaled between [0,1].
%   'eval'     freeform operations
%
% To enable undoing and applying the exactly same normalization to
% other data sets, normalization information is saved into a
% normalization struct, which has the fields:
%
%   .type   ; struct type, ='som_norm'
%   .method ; normalization method, a string
%   .params ; normalization parameters
%   .status ; string: 'uninit', 'undone' or 'done'
%
% Normalizations are always one-variable operations. In the data and map
% structs the normalization information for each component is saved in the
% '.comp_norm' field, which is a cell array of length dim. Each cell
% contains normalizations for one vector component in a struct array of
% normalization structs. Each component may have different amounts of
% different kinds of normalizations. Typically, all normalizations are
% either 'undone' or 'done', but in special situations this may not be the
% case. The easiest way to check out the status of the normalizations is to
% use function SOM_INFO, e.g. som_info(sS,3)
%
% REQUIRED INPUT ARGUMENTS
%
%   sS                The data to which the normalization is applied.
%            (struct) Data or map struct. Before adding any new
%                     normalizations, it is ensured that the
%                     normalizations for the specified components in the
%                     '.comp_norm' field have status 'done'.
%            (matrix) data matrix
%
% OPTIONAL INPUT ARGUMENTS
%
%   method            The normalization(s) to add/use. If missing,
%                     or an empty variable ('' or []) is given, the
%                     normalizations in the data struct are used.
%            (string) Identifier for a normalization method to be added:
%                     'var', 'range', 'log', 'logistic', 'histD' or 'histC'. The
%                     same method is applied to all specified components
%                     (given in comps). The normalizations are first
%                     initialized (for each component separately, of
%                     course) and then applied.
%            (struct) Normalization struct, or an array of structs, which
%                     is applied to all specified components. If the
%                     '.status' field of the struct(s) is 'uninit',
%                     the normalization(s) is initialized first.
%                     Alternatively, the struct may be map or data struct
%                     in which case its '.comp_norm' field is used
%                     (see the cell array option below).
%            (cell array) In practice, the '.comp_norm' field of
%                     a data/map struct. The length of the array
%                     must be equal to the dimension of the given
%                     data set (sS). Each cell contains the
%                     normalization(s) for one component. Only the
%                     normalizations listed in comps argument are
%                     applied though.
%            (cellstr array) norm and denorm operations in a cellstr array
%                     which are evaluated with EVAL command with variable
%                     name 'x' reserved for the variable.
%
%   comps    (vector) The components to which the normalization(s) is
%                     applied. Default is to apply to all components.
%
% OUTPUT ARGUMENTS
%
%   sS                Modified and/or updated data.
%            (struct) If a struct was given as input argument, the
%                     same struct is returned with normalized data and
%                     updated '.comp_norm' fields.
%            (matrix) If a matrix was given as input argument, the
%                     normalized data matrix is returned.
%
% EXAMPLES
%
%  To add (initialize and apply) a normalization to a data struct:
%
%    sS = som_normalize(sS,'var');
%
%  This uses 'var'-method to all components. To add a method only to
%  a few selected components, use the comps argument:
%
%    sS = som_normalize(sS,'log',[1 3:5]);
%
%  To ensure that all normalization operations have indeed been done:
%
%    sS = som_normalize(sS);
%
%  The same for only a few components:
%
%    sS = som_normalize(sS,'',[1 3:5]);
%
%  To apply the normalizations of a data struct sS to a new data set D:
%
%    D = som_normalize(D,sS);
%  or
%    D = som_normalize(D,sS.comp_norm);
%
%  To normalize a data set:
%
%    D = som_normalize(D,'histD');
%
%  Note that in this case the normalization information is lost.
%
%  To check out the status of normalization in a struct use SOM_INFO:
%
%    som_info(sS,3)
%
%
% SEE ALSO
%
%  som_denormalize    Undo normalizations of a data struct/set.
%  som_norm_variable  Normalization operations for a set of scalar values.
%  som_info           User-friendly information of SOM Toolbox structs.

% Copyright (c) 1998-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 151199 150500

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check arguments

error(nargchk(1, 3, nargin));  % check no. of input arguments is correct

% sD
struct_mode = isstruct(sD);
if struct_mode,
  switch sD.type
   case 'som_map', D = sD.codebook;
   case 'som_data', D = sD.data;
   otherwise, error('Illegal struct.')
  end
else
  D = sD;
end
[dlen dim] = size(D);

% comps
if nargin<3 || (ischar(comps) && strcmp(comps,'all')),
  comps = [1:dim];
end
if isempty(comps), return; end
if size(comps,1)>1, comps = comps'; end  % make it a row vector

% method
csNorm = cell(dim,1);
if nargin<2 || isempty(method),
  if ~struct_mode,
    warning('No normalization method given. Data left unchanged.');
    return;
  end
  method = '';
else
  % check out the given method
  % (and if necessary, copy it for each specified component)
  if ischar(method),
    switch method,
    case {'var','range','log','histD','histC','logistic'},
      sN = som_set('som_norm','method',method);
    otherwise,
      error(['Unrecognized method: ' method]);
    end
    for i=comps, csNorm{i} = sN; end
  elseif isstruct(method),
    switch method(1).type,
    case {'som_map','som_data'}, csNorm = method(1).comp_norm;
    case {'som_norm'}, for i=comps, csNorm{i} = method; end
    otherwise,
      error('Invalid struct given as normalization method.')
    end
  elseif iscellstr(method),
    [dummy,sN] = som_norm_variable(1,method,'init');
    for i=comps, csNorm{i} = sN; end
  elseif iscell(method),
    csNorm = method;
  else
    error('Illegal method argument.')
  end
  % check the size of csNorm is the same as data dimension
  if length(csNorm) ~= dim,
    error('Given number of normalizations does not match data dimension.')
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialize

% make sure all the current normalizations for current
% components have been done
if struct_mode,
  alldone = 1;
  for i = comps,
    for j=1:length(sD.comp_norm{i}),
      sN = sD.comp_norm{i}(j);
      if ~strcmp(sN.status,'done'),
	alldone = 0;
        [x,sN] = som_norm_variable(D(:,i), sN, 'do');
        D(:,i) = x;
        sD.comp_norm{i}(j) = sN;
      end
    end
  end
  if isempty(method),
    if alldone,
      warning('No ''undone'' normalizations found. Data left unchanged.');
    else
      fprintf(1,'Normalizations have been redone.\n');
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% action

% add the new normalizations to the old ones
for i = comps,
  if ~isempty(csNorm{i}),
    [x,sN] = som_norm_variable(D(:,i), csNorm{i}, 'do');
    D(:,i) = x;
    if struct_mode,
      if isempty(sD.comp_norm{i}), sD.comp_norm{i} = sN;
      else sD.comp_norm{i} = [sD.comp_norm{i}, sN]; end
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% output

if struct_mode,
  switch sD.type
   case 'som_map', sD.codebook = D;
   case 'som_data', sD.data = D;
   otherwise, error('Illegal struct.')
  end
else
  sD = D;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: ilarinieminen

repo: SOM-Toolbox


text:
function [tpr,tnr,info] = vl_roc(labels, scores, varargin)
%VL_ROC   ROC curve.
%   [TPR,TNR] = VL_ROC(LABELS, SCORES) computes the Receiver Operating
%   Characteristic (ROC) curve [1]. LABELS is a row vector of ground
%   truth labels, greater than zero for a positive sample and smaller
%   than zero for a negative one. SCORES is a row vector of
%   corresponding sample scores, usually obtained from a
%   classifier. The scores induce a ranking of the samples where
%   larger scores should correspond to positive labels.
%
%   Without output arguments, the function plots the ROC graph of the
%   specified data in the current graphical axis.
%
%   Otherwise, the function returns the true positive and true
%   negative rates TPR and TNR. These are vectors of the same size of
%   LABELS and SCORES and are computed as follows. Samples are ranked
%   by decreasing scores, starting from rank 1. TPR(K) and TNR(K) are
%   the true positive and true negative rates when samples of rank
%   smaller or equal to K-1 are predicted to be positive. So for
%   example TPR(3) is the true positive rate when the two samples with
%   largest score are predicted to be positive. Similarly, TPR(1) is
%   the true positive rate when no samples are predicted to be
%   positive, i.e. the constant 0.
%
%   Setting a label to zero ignores the corresponding sample in the
%   calculations, as if the sample was removed from the data. Setting
%   the score of a sample to -INF causes the function to assume that
%   that sample was never retrieved. If there are samples with -INF
%   score, the ROC curve is incomplete as the maximum recall is less
%   than 1.
%
%   [TPR,TNR,INFO] = VL_ROC(...) returns an additional structure INFO
%   with the following fields:
%
%   info.auc:: Area under the ROC curve (AUC).
%     This is the area under the ROC plot, the parametric curve
%     (FPR(S), TPR(S)). The PLOT option can be used to plot variants
%     of this curve, which affects the calculation of a corresponding
%     AUC.
%
%   info.eer:: Equal error rate (EER).
%     The equal error rate is the value of FPR (or FNR) when the ROC
%     curves intersects the line connecting (0,0) to (1,1).
%
%   info.eerThreshold:: EER threshold.
%     The value of the score for which the EER is attained.
%
%   VL_ROC() accepts the following options:
%
%   Plot:: []
%     Setting this option turns on plotting unconditionally. The
%     following plot variants are supported:
%
%     tntp:: Plot TPR against TNR (standard ROC plot).
%     tptn:: Plot TNR against TPR (recall on the horizontal axis).
%     fptp:: Plot TPR against FPR.
%     fpfn:: Plot FNR against FPR (similar to a DET curve).
%
%     Note that this option will affect the INFO.AUC value computation
%     too.
%
%   NumPositives:: []
%   NumNegatives:: []
%     If either of these parameters is set to a number, the function
%     pretends that LABELS contains the specified number of
%     positive/negative labels. NUMPOSITIVES/NUMNEGATIVES cannot be
%     smaller than the actual number of positive/negative entries in
%     LABELS. The additional positive/negative labels are appended to
%     the end of the sequence as if they had -INF scores (as explained
%     above, the function interprets such samples as `not
%     retrieved'). This feature can be used to evaluate the
%     performance of a large-scale retrieval experiment in which only
%     a subset of highly-scoring results are recorded for efficiency
%     reason.
%
%   Stable:: false
%     If set to true, TPR and TNR are returned in the same order
%     of LABELS and SCORES rather than being sorted by decreasing
%     score.
%
%   About the ROC curve::
%     Consider a classifier that predicts as positive all samples whose
%     score is not smaller than a threshold S. The ROC curve represents
%     the performance of such classifier as the threshold S is
%     changed. Formally, define
%
%       P = overall num. of positive samples,
%       N = overall num. of negative samples,
%
%     and for each threshold S
%
%       TP(S) = num. of samples that are correctly classified as positive,
%       TN(S) = num. of samples that are correctly classified as negative,
%       FP(S) = num. of samples that are incorrectly classified as positive,
%       FN(S) = num. of samples that are incorrectly classified as negative.
%
%     Consider also the rates:
%
%       TPR = TP(S) / P,      FNR = FN(S) / P,
%       TNR = TN(S) / N,      FPR = FP(S) / N,
%
%     and notice that, by definition,
%
%       P = TP(S) + FN(S) ,    N = TN(S) + FP(S),
%       1 = TPR(S) + FNR(S),   1 = TNR(S) + FPR(S).
%
%     The ROC curve is the parametric curve (FPR(S), TPR(S)) obtained
%     as the classifier threshold S is varied in the reals. The TPR is
%     the same as `recall' in a PR curve (see VL_PR()).
%
%     The ROC curve is contained in the square with vertices (0,0) The
%     (average) ROC curve of a random classifier is a line which
%     connects (1,0) and (0,1).
%
%     The ROC curve is independent of the prior probability of the
%     labels (i.e. of P/(P+N) and N/(P+N)).
%
%   REFERENCES:
%   [1] http://en.wikipedia.org/wiki/Receiver_operating_characteristic
%
%   See also: VL_PR(), VL_DET(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

[tp, fp, p, n, perm, varargin] = vl_tpfp(labels, scores, varargin{:}) ;
opts.plot = [] ;
opts.stable = false ;
opts = vl_argparse(opts,varargin) ;

% compute the rates
small = 1e-10 ;
tpr = tp / max(p, small) ;
fpr = fp / max(n, small) ;
fnr = 1 - tpr ;
tnr = 1 - fpr ;

do_plots = ~isempty(opts.plot) || nargout == 0 ;
if isempty(opts.plot), opts.plot = 'fptp' ; end

% --------------------------------------------------------------------
%                                                      Additional info
% --------------------------------------------------------------------

if nargout > 2 || do_plots
  % Area under the curve. Since the curve is a staircase (in the
  % sense that for each sample either tn is decremented by one
  % or tp is incremented by one but the other remains fixed),
  % the integral is particularly simple and exact.

  switch opts.plot
    case 'tntp', info.auc = -sum(tpr .* diff([0 tnr])) ;
    case 'fptp', info.auc = +sum(tpr .* diff([0 fpr])) ;
    case 'tptn', info.auc = +sum(tnr .* diff([0 tpr])) ;
    case 'fpfn', info.auc = +sum(fnr .* diff([0 fpr])) ;
    otherwise
      error('''%s'' is not a valid PLOT type.', opts.plot);
  end

  % Equal error rate. One must find the index S in correspondence of
  % which TNR(S) and TPR(s) cross. Note that TPR(S) is non-decreasing,
  % TNR(S) is non-increasing, and from rank S to rank S+1 only one of
  % the two quantities can change. Hence there are exactly two types
  % of crossing points:
  %
  %  1) TNR(S) = TNR(S+1) = EER and TPR(S) <= EER, TPR(S+1) > EER,
  %  2) TPR(S) = TPR(S+1) = EER and TNR(S) > EER, TNR(S+1) <= EER.
  %
  % Moreover, if the maximum TPR is smaller than 1, then it is
  % possible that neither of the two cases realizes. In the latter
  % case, we return EER=NaN.

  s = max(find(tnr > tpr)) ;
  if s == length(tpr)
    info.eer = NaN ;
    info.eerThreshold = 0 ;
  else
    if tpr(s) == tpr(s+1)
      info.eer = 1 - tpr(s) ;
    else
      info.eer = 1 - tnr(s) ;
    end
    info.eerThreshold = scores(perm(s)) ;
  end
end

% --------------------------------------------------------------------
%                                                                 Plot
% --------------------------------------------------------------------

if do_plots
  cla ; hold on ;
  switch lower(opts.plot)
    case 'tntp'
      hroc = plot(tnr, tpr, 'b', 'linewidth', 2) ;
      hrand = spline([0 1], [1 0], 'r--', 'linewidth', 2) ;
      spline([0 1], [0 1], 'k--', 'linewidth', 1) ;
      plot(1-info.eer, 1-info.eer, 'k*', 'linewidth', 1) ;
      xlabel('true negative rate') ;
      ylabel('true positive rate (recall)') ;
      loc = 'sw' ;

    case 'fptp'
      hroc = plot(fpr, tpr, 'b', 'linewidth', 2) ;
      hrand = spline([0 1], [0 1], 'r--', 'linewidth', 2) ;
      spline([1 0], [0 1], 'k--', 'linewidth', 1) ;
      plot(info.eer, 1-info.eer, 'k*', 'linewidth', 1) ;
      xlabel('false positive rate') ;
      ylabel('true positive rate (recall)') ;
      loc = 'se' ;

    case 'tptn'
      hroc = plot(tpr, tnr, 'b', 'linewidth', 2) ;
      hrand = spline([0 1], [1 0], 'r--', 'linewidth', 2) ;
      spline([0 1], [0 1], 'k--', 'linewidth', 1) ;
      plot(1-info.eer, 1-info.eer, 'k*', 'linewidth', 1) ;
      xlabel('true positive rate (recall)') ;
      ylabel('false positive rate') ;
      loc = 'sw' ;

    case 'fpfn'
      hroc = plot(fpr, fnr, 'b', 'linewidth', 2) ;
      hrand = spline([0 1], [1 0], 'r--', 'linewidth', 2) ;
      spline([0 1], [0 1], 'k--', 'linewidth', 1) ;
      plot(info.eer, info.eer, 'k*', 'linewidth', 1) ;
      xlabel('false positive (false alarm) rate') ;
      ylabel('false negative (miss) rate') ;
      loc = 'ne' ;

    otherwise
      error('''%s'' is not a valid PLOT type.', opts.plot);
  end

  grid on ;
  xlim([0 1]) ;
  ylim([0 1]) ;
  axis square ;
  title(sprintf('ROC (AUC: %.2f%%, EER: %.2f%%)', info.auc * 100, info.eer * 100), ...
        'interpreter', 'none') ;
  legend([hroc hrand], 'ROC', 'ROC rand.', 'location', loc) ;
end

% --------------------------------------------------------------------
%                                                        Stable output
% --------------------------------------------------------------------

if opts.stable
  tpr(1) = [] ;
  tnr(1) = [] ;
  tpr_ = tpr ;
  tnr_ = tnr ;
  tpr = NaN(size(tpr)) ;
  tnr = NaN(size(tnr)) ;
  tpr(perm) = tpr_ ;
  tnr(perm) = tnr_ ;
end

% --------------------------------------------------------------------
function h = spline(x,y,spec,varargin)
% --------------------------------------------------------------------
prop = vl_linespec2prop(spec) ;
h = line(x,y,prop{:},varargin{:}) ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function PVal=computeCoherenceSignificance(TH,N)
% PVal=computeCoherenceSignificance(TH,N)
%
% PURPOSE:
% Returns a rough estimate of the significance value of a
% coherence threshold. 'Rough' because it assumes
% gaussian noise and takes no account of the fact that voxels
% are spatially correlated with their neighbours
% and have non-delta temporal auto-correlation functions.
%
% From Bantettini el al, 1993,
% Processing Strategies for Time Course Data Sets
% MRM 30:161-173 (1993) pp 171
%
% ARW 01/14/03 Wrote it
% $Author: wade $
% $Date: 2003/01/14 23:07:47 $

if (~exist('TH','var') | ~exist('N','var'))
    error('You must pass both a coherence threshold "TH" and the number of TRs "N"');
end

% It's just the complimentary error function...
PVal=1-erf(TH*sqrt(N/2))+eps;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function z = r82_to_c8 ( x, y )

%*****************************************************************************80
%
%% R82_TO_C8 makes a complex number from two real values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the real and imaginary parts.
%
%    Output, complex Z, the value of ( A + B*i).
%
  z = x + y * i;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_bug2148

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_connectivitysimulation ft_freqanalysis ft_connectivityanalysis ft_connectivityplot ft_freqdescriptives ft_checkdata

% ft_checkdata (ie fixdimord) always wants to reverts dimord to 'chan'
% although 'labelcmb' is present. The only way to make e.g.
% ft_connectivitplot work atm is to use:
% ft_checkdata(data, 'cmbstyle', 'full');
% This will make dimord chan_chan_freq and convert labelcmb to label.

cfg             = [];
cfg.ntrials     = 500;
cfg.triallength = 1;
cfg.fsample     = 200;
cfg.nsignal     = 4;
cfg.method      = 'ar';
cfg.params(:,:,1) = [.8   0   0 .2;
  0  .9  .5  0;
  .4   0  .5 .3;
  0  .2   0 .7];
cfg.params(:,:,2) = [-.5   0   0  0;
  0 -.8   0  0;
  0   0 -.2  0;
  0   0   0 .1];
cfg.noisecov      = [.3 0  0  0;
  0 1  0  0;
  0 0 .2  0;
  0 0  0  .4];

data            = ft_connectivitysimulation(cfg);
data.label      = {'AFz', 'Cz', 'Fz', 'Oz'};

% freqanalysis
cfgf           = [];
cfgf.method    = 'mtmfft';
cfgf.output    = 'fourier';
cfgf.tapsmofrq = 2;
cfgf.keeptrials= 'yes';
freq           = ft_freqanalysis(cfgf, data);
ft_freqdescriptives([], freq);

% connectivityanalysis
cfgc           = [];
cfgc.channelcmb = {'all' freq.label{1}};
cfgc.method    = 'coh';
c1             = ft_connectivityanalysis(cfgc, freq);
ft_connectivityplot([], c1);

% plotting
cfgp = [];
cfgp.layout    = 'EEG1010.lay';
cfgp.directionality = 'inflow';
cfgp.parameter = 'cohspctrm';
cfgp.refchannel = 'AFz';
ft_topoplotTFR(cfgp, c1);

% do JM's juicy transforms
tmp = ft_checkdata(freq, 'cmbstyle', 'fullfast');
c1             = ft_connectivityanalysis(cfgc, freq);
try
  ft_freqdescriptives([], tmp);
catch err
  fprintf('error: %s @ %s @ line %i\n ', err.message, err.stack(1).name, err.stack(1).line);
end

ft_topoplotTFR(cfgp, c1);
ft_connectivityplot([], c1);

cfga = [];
cfga.parameter = 'fourierspctrm';
try
  ft_appendfreq(cfga, freq, tmp);
catch err
  fprintf('error: %s @ %s @ line %i\n ', err.message, err.stack(1).name, err.stack(1).line);
end
vol = [];
vol.o = [0 0 4];
vol.r = 12;
vol.unit = 'cm';
vol.cond = 0.5;
vol.type = 'singlesphere';
cfgs.headmodel = vol;
cfgs.frequency = 10;
cfgs.elecfile = 'standard_1020.elc';
ft_sourceanalysis(cfgs, tmp);

tmp = ft_checkdata(tmp, 'cmbstyle', 'sparse');
c1             = ft_connectivityanalysis(cfgc, freq);
try
  ft_freqdescriptives([], tmp);
catch err
  fprintf('error: %s @ %s @ line %i\n ', err.message, err.stack(1).name, err.stack(1).line);
end
ft_topoplotTFR(cfgp, c1);
ft_connectivityplot([], c1);
try
  ft_appendfreq(cfga, freq, tmp);
catch err
  fprintf('error: %s @ %s @ line %i\n ', err.message, err.stack(1).name, err.stack(1).line);
end

vol = [];
vol.o = [0 0 4];
vol.r = 12;
vol.cond = 0.5;
vol.unit = 'cm';
vol.type = 'singlesphere';
cfgs.headmodel = vol;
cfgs.frequency = 10;
cfgs.elecfile = 'standard_1020.elc';
ft_sourceanalysis(cfgs, tmp);

tmp = ft_checkdata(tmp, 'cmbstyle', 'sparsewithpow');
c1             = ft_connectivityanalysis(cfgc, freq);
ft_freqdescriptives([], tmp);
ft_topoplotTFR(cfgp, c1);
ft_connectivityplot([], c1);
try
  ft_appendfreq(cfga, freq, tmp);
catch err
  fprintf('error: %s @ %s @ line %i\n ', err.message, err.stack(1).name, err.stack(1).line);
end

vol = [];
vol.o = [0 0 4];
vol.r = 12;
vol.cond = 0.5;
vol.unit = 'cm';
vol.type = 'singlesphere';
cfgs.headmodel = vol;
cfgs.elecfile = 'standard_1020.elc';
cfgs.frequency = 10;
ft_sourceanalysis(cfgs, tmp);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ stats ] = ml_get_cv_grid_states(test,train)
%ML_GET_CV_GRID_STATS Compute statistics on K-fold Cross Validation (CV).
%
%   input -----------------------------------------------------------------
%
%       o test  : cell array (M x 1), M is the number of parameters for
%                                     which K-fold CV has been evaluated on.
%
%           test{i}.accuracy  : (1 x K), K number of
%           test{i}.precision : (P x K), P is the total number of classes
%           test{i}.recall    : (P x K)
%           test{i}.fmeasure  : (P x K)
%
%
%       o train : cell array (M x 1),
%
%               train{i}: structure is the same as for test
%
%   output ----------------------------------------------------------------
%
%       o stats : struct , returns statistics on train and test sets
%
%           stats.train  : (statistics on the training data set)
%
%               - stats.train.acc : (statistics on the accuracy for the training set)
%
%                   - stats.train.acc.mean
%                   - stats.train.acc.std
%
%       the structure is the same for the test statistics
%
%% Check if grid search was done on > 1 parameters

[P, N] = size(train);

if N > 1
    test = test(:);
    train = train(:);
end

%% Extract Statitistics
M = length(train);
K = length(train{1}.accuracy);

% For Training Data
stats.train.acc.mean = zeros(1,K);
stats.train.acc.std  = zeros(1,K);

stats.train.fmeasure.mean = zeros(1,K);
stats.train.fmeasure.std  = zeros(1,K);

stats.train.fpr.mean = zeros(1,K);
stats.train.fpr.std  = zeros(1,K);

stats.train.tnr.mean = zeros(1,K);
stats.train.tnr.std  = zeros(1,K);

% bla bli
% Model Statistics for SVM
stats.model.totSV.mean   = zeros(1,K);
stats.model.totSV.std    = zeros(1,K);

stats.model.ratioSV.mean = zeros(1,K);
stats.model.ratioSV.std  = zeros(1,K);

stats.model.posSV.mean   = zeros(1,K);
stats.model.posSV.std    = zeros(1,K);

stats.model.negSV.mean   = zeros(1,K);
stats.model.negSV.std    = zeros(1,K);

stats.model.boundSV.mean = zeros(1,K);
stats.model.boundSV.std  = zeros(1,K);

% For Testing Data
if ~isempty(test)
    stats.test.acc.mean = zeros(1,K);
    stats.test.acc.std  = zeros(1,K);

    stats.test.fmeasure.mean = zeros(1,K);
    stats.test.fmeasure.std  = zeros(1,K);

    stats.test.fpr.mean = zeros(1,K);
    stats.test.fpr.std  = zeros(1,K);

    stats.test.tnr.mean = zeros(1,K);
    stats.test.tnr.std  = zeros(1,K);
end

for m=1:M

    % For Testing Data
    if ~isempty(test)
        stats.test.acc.mean(m)      = mean(test{m}.accuracy);
        stats.test.acc.std(m)       = std(test{m}.accuracy);

        stats.test.fmeasure.mean(m) = mean(test{m}.fmeasure);
        stats.test.fmeasure.std(m)  = std(test{m}.fmeasure);

        stats.test.fpr.mean(m)      = mean(test{m}.fpr);
        stats.test.fpr.std(m)       = std(test{m}.fpr);

        stats.test.tnr.mean(m)      = mean(test{m}.tnr);
        stats.test.tnr.std(m)       = std(test{m}.tnr);
    end

    % For Training Data
    stats.train.acc.mean(m)      = mean(train{m}.accuracy);
    stats.train.acc.std(m)       = std(train{m}.accuracy);

    stats.train.fmeasure.mean(m) = mean(train{m}.fmeasure);
    stats.train.fmeasure.std(m)  = std(train{m}.fmeasure);

    stats.train.fpr.mean(m)      = mean(train{m}.fpr);
    stats.train.fpr.std(m)       = std(train{m}.fpr);

    stats.train.tnr.mean(m)      = mean(train{m}.tnr);
    stats.train.tnr.std(m)       = std(train{m}.tnr);

    % Model Statistics for SVM
    stats.model.totSV.mean(m)   = mean(train{m}.totSV);
    stats.model.totSV.std(m)    = std(train{m}.totSV);

    stats.model.ratioSV.mean(m) = mean(train{m}.ratioSV);
    stats.model.ratioSV.std(m)  = std(train{m}.ratioSV);

    stats.model.posSV.mean(m)   = mean(train{m}.posSV);
    stats.model.posSV.std(m)    = std(train{m}.posSV);

    stats.model.negSV.mean(m)   = mean(train{m}.negSV);
    stats.model.negSV.std(m)    = std(train{m}.negSV);

    stats.model.boundSV.mean(m) = mean(train{m}.boundSV);
    stats.model.boundSV.std(m)  = std(train{m}.boundSV);

end

%% Reshape stats if grid search on > 1 parameters
if N > 1
    % For Testing Data
    if ~isempty(test)
        stats.test.acc.mean = reshape(stats.test.acc.mean, P, N);
        stats.test.acc.std  = reshape(stats.test.acc.std, P, N);

        stats.test.fmeasure.mean = reshape(stats.test.fmeasure.mean, P, N);
        stats.test.fmeasure.std  = reshape(stats.test.fmeasure.std, P, N);

        stats.test.fpr.mean = reshape(stats.test.fpr.mean, P, N);
        stats.test.fpr.std  = reshape(stats.test.fpr.std, P, N);

        stats.test.tnr.mean = reshape(stats.test.tnr.mean, P, N);
        stats.test.tnr.std  = reshape(stats.test.tnr.std, P, N);

    end

    % For Training Data
    stats.train.acc.mean = reshape(stats.train.acc.mean, P, N);
    stats.train.acc.std  = reshape(stats.train.acc.std, P, N);

    stats.train.fmeasure.mean = reshape(stats.train.fmeasure.mean, P, N);
    stats.train.fmeasure.std  = reshape(stats.train.fmeasure.std, P, N);

    stats.train.fpr.mean = reshape(stats.train.fpr.mean, P, N);
    stats.train.fpr.std  = reshape(stats.train.fpr.std, P, N);

    stats.train.tnr.mean = reshape(stats.train.tnr.mean, P, N);
    stats.train.tnr.std  = reshape(stats.train.tnr.std, P, N);

    % For model Statistics
    stats.model.totSV.mean   = reshape(stats.model.totSV.mean, P, N);
    stats.model.totSV.std    = reshape(stats.model.totSV.std, P, N);

    stats.model.ratioSV.mean = reshape(stats.model.ratioSV.mean, P, N);
    stats.model.ratioSV.std  = reshape(stats.model.ratioSV.std, P, N);

    stats.model.posSV.mean   = reshape(stats.model.posSV.mean, P, N);
    stats.model.posSV.std    = reshape(stats.model.posSV.std, P, N);

    stats.model.negSV.mean   = reshape(stats.model.negSV.mean, P, N);
    stats.model.negSV.std    = reshape(stats.model.negSV.std, P, N);

    stats.model.boundSV.mean = reshape(stats.model.boundSV.mean, P, N);
    stats.model.boundSV.std  = reshape(stats.model.boundSV.std, P, N);

end

end

author: epfl-lasa


text:
% Author: Junaid Ahmed Ansari @ Robotics Research Center, IIIT Hyderabad
% Email : junaid.ansari@research.iiit.ac.in
function generateCandidateFor2D2Dmatching(scoreMat_3D2D, thres, fpath, fno, detectionsQ, detectionsT)

    fname = sprintf('%s%06d.txt',fpath,fno);
    fid   = fopen(fname,'w+');

    Q_dno = zeros(length(detectionsQ), 1);
    T_dno = zeros(length(detectionsT), 1);

    % save detection numbers into array for speed
    for i = 1:length(detectionsQ)
        Q_dno(i) = detectionsQ{i}.dno;
    end

    for i = 1:length(detectionsT)
        T_dno(i) = detectionsT{i}.dno;
    end

    % find candidate matches based on 3D2D for 2D2D matching
    for i = 1:length(detectionsQ)
        i
        cand_match = [];
        cand_match = T_dno(find(scoreMat_3D2D(i,:) > thres))
        cand_match = [detectionsQ{i}.dno; length(cand_match); cand_match]'
        fprintf(fid,'%d ',cand_match);
        fprintf(fid,'\n');
    end

    fclose(fid);
end

repo: MOTBeyondPixels


text:
% textgui() - make sliding vertical window. This window contain text
%             with optional function calls at each line.
%
% Usage:
%   >> textgui( commandnames, helparray, 'key', 'val' ...);
%
% Inputs:
%   commandnames - name of the commands. Either char array or cell
%                  array of char. All style are 'pushbuttons' exept
%                  for empty commands.
%   helparray    - cell array of commands to execute for each menu
%                  (default is empty)
%
% Optional inputs:
%  'title'        - window title
%  'fontweight'   - font weight. Single value or cell array of value for each line.
%  'fontsize'     - font size. Single value or cell array of value for each line.
%  'fontname'     - font name. Single value or cell array of value for each line.
%  'linesperpage' - number of line per page. Default is 20.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% Example:
%  textgui({ 'function1' 'function2' }, {'hthelp(''function1'')' []});
%  % this function will call a pop_up window with one button on
%  % 'function1' which will call the help of this function.
%

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function tmp = textgui( textmenu, helparray, varargin);

XBORDERS   = 0; % pixels
TOPORDINATE = 1.09;
TOPLINES    = 2;

if nargin < 1
	help textgui;
	return;
end;
if nargin <2
	helparray = cell(1, length(	textmenu ));
end;

% optional parameters
% -------------------
if nargin >2
    for i = 1:length(varargin)
        if iscell(varargin{i}), varargin(i) = { varargin(i) }; end;
    end;
    g=struct(varargin{:});
else
    g = [];
end;
try, g.title;    catch, g.title = ''; end;
try, g.fontname; catch, g.fontname = 'courier'; end;
try, g.fontsize; catch, g.fontsize = 12; end;
try, g.fontweight; catch, g.fontweight = 'normal'; end;
try, g.linesperpage; catch, g.linesperpage = 20; end;

if isempty( helparray )
	helparray = cell(1,200);
else
	tmpcell = cell(1,200);
	helparray = { helparray{:} tmpcell{:} };
end;

% number of elements
% ------------------
if iscell(textmenu)	nblines = length(textmenu);
else				nblines = size(textmenu,1);
end;

% generating the main figure
% --------------------------
fig = figure('position', [100 100 800 25*15], 'menubar', 'none', 'numbertitle', 'off', 'name', 'Help', 'color', 'w');
pos = get(gca,'position'); % plot relative to current axes
q = [pos(1) pos(2) 0 0];
s = [pos(3) pos(4) pos(3) pos(4)]./100;

% generating cell arrays
% ----------------------
if isstr(g.fontname),     tmp1(1:nblines) = {g.fontname}; g.fontname = tmp1; end;
if isnumeric(g.fontsize), tmp2(1:nblines) = {g.fontsize}; g.fontsize = tmp2; end;
if isstr(g.fontweight),   tmp3(1:nblines) = {g.fontweight}; g.fontweight = tmp3; end;
switch g.fontname{1}
    case 'courrier', CHAR_WIDTH = 11; % pixels
    case 'times', CHAR_WIDTH = 11; % pixels
    otherwise, CHAR_WIDTH = 11;
end;

topordi = TOPORDINATE;
if ~isempty(g.title)
    addlines = size(g.title,1)+1+TOPLINES;
    if nblines+addlines < g.linesperpage
        divider = g.linesperpage;
    else
        divider = nblines+addlines;
    end;
    ordinate      = topordi-topordi*TOPLINES/divider;
    currentheight = topordi/divider;

    for index=1:size(g.title,1)
      ordinate      = topordi-topordi*(TOPLINES+index-1)/divider;
      h = text( -0.1, ordinate, g.title(index,:), 'unit', 'normalized', 'horizontalalignment', 'left', ...
		 'fontname', g.fontname{1}, 'fontsize', g.fontsize{1},'fontweight', fastif(index ==1, 'bold', ...
		  'normal'), 'interpreter', 'none' );
    end;
    %h = uicontrol( gcf, 'unit', 'normalized', 'style', 'text', 'backgroundcolor', get(gcf, 'color'), 'horizontalalignment', 'left', ...
    %		'position', [-0.1 ordinate 1.1 currentheight].*s*100+q, 'string', g.title, ...
    %	       'fontsize', g.fontsize{1}+1, 'fontweight', 'bold', 'fontname', g.fontname{1});
 else
    addlines = TOPLINES;
    if nblines+addlines < g.linesperpage
        divider = g.linesperpage;
    else
        divider = nblines+addlines;
    end;
end;

maxlen = size(g.title,2);
for i=1:nblines
	if iscell(textmenu)	tmptext = textmenu{i};
	else			tmptext = textmenu(i,:);
	end;
	ordinate      = topordi-topordi*(i-1+addlines)/divider;
	currentheight = topordi/divider;

	if isempty(helparray{i})
		h = text( -0.1, ordinate, tmptext, 'unit', 'normalized', 'horizontalalignment', 'left', ...
		 'fontname', g.fontname{i}, 'fontsize', g.fontsize{i},'fontweight', g.fontweight{i}, 'interpreter', 'none' );
%		h = uicontrol( gcf, 'unit', 'normalized', 'style', 'text', 'backgroundcolor', get(gcf, 'color'), 'horizontalalignment', 'left', ...
%			'position', [-0.1 ordinate 1.1 currentheight].*s*100+q, 'string', tmptext, 'fontname', g.fontname{i}, 'fontsize', g.fontsize{i},'fontweight', g.fontweight{i} );
	else
		h = text( -0.1, ordinate, tmptext, 'unit', 'normalized', 'color', 'b', 'horizontalalignment', 'left', ...
		 'fontname', g.fontname{i}, 'buttondownfcn', helparray{i}, 'fontsize', g.fontsize{i},'fontweight', g.fontweight{i}, 'interpreter', 'none' );
%		h = uicontrol( gcf,  'callback', helparray{i}, 'unit','normalized', 'style', 'pushbutton',  'horizontalalignment', 'left', ...
%			'position', [-0.1 ordinate 1.1 currentheight].*s*100+q, 'string', tmptext, 'fontname', g.fontname{i}, 'fontsize', g.fontsize{i},'fontweight', g.fontweight{i});
	end;
	maxlen = max(maxlen, length(tmptext));
end;
pos = get(gcf, 'position');
set(gcf, 'position', [pos(1:2) maxlen*CHAR_WIDTH+2*XBORDERS 400]);
% = findobj('parent', gcf);
%set(h, 'fontname', 'courier');

if nblines+addlines < g.linesperpage
    zooming  = 1;
else
    zooming  = (nblines+addlines)/g.linesperpage;
end;

slider(gcf, 0, 1, 1, zooming, 0);
axis off;

author: PatternRecognition

repo: OpenBMI


text:
function [ a, seed ] = npart_rsf_lex_random ( n, npart, seed )

%*****************************************************************************80
%
%% NPART_RSF_LEX_RANDOM returns a random RSF NPART partition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the integer to be partitioned.
%    N must be positive.
%
%    Input, integer NPART, the number of parts of the partition.
%    1 <= NPART <= N.
%
%    Input/output, integer SEED, a seed for the random number
%    generator.
%
%    Output, integer A(NPART), contains the partition.
%    A(1) through A(NPART) contain the nonzero integers which
%    sum to N.
%
  npartitions = npart_enum ( n, npart );

  [ rank, seed ] = i4_uniform ( 1, npartitions, seed );

  a = npart_rsf_lex_unrank ( rank, n, npart );

  return
end


text:
function do_all(dbName)

%DO_ALL execute all codes in the ttimes package.
% DO_ALL(dbName) runs each of the codes in the ttimes package on the
% database specified by dbName. See the ttimes.dbload function for a
% description.
%
% See also ttimres.dbload

ttimes.map(dbName);
ttimes.depth_section(dbName);
ttimes.tt_curve(dbName);
ttimes.arrival_histogram(dbName);
ttimes.write_lotos(dbName);
ttimes.dbload(dbName)

author: geoscience-community-codes

repo: GISMO


text:
function [fnormhat,t]=instfreq(x,t,L,trace);
%INSTFREQ Instantaneous frequency estimation.
%	[FNORMHAT,T]=INSTFREQ(X,T,L,TRACE) computes the instantaneous
%	frequency of the analytic signal X at time instant(s) T, using the
%	trapezoidal integration rule.
%	The result FNORMHAT lies between 0.0 and 0.5.
%
%	X : Analytic signal to be analyzed.
%	T : Time instants	        (default : 2:length(X)-1).
%	L : If L=1, computes the (normalized) instantaneous frequency
%	    of the signal X defined as angle(X(T+1)*conj(X(T-1)) ;
%	    if L>1, computes a Maximum Likelihood estimation of the
%	    instantaneous frequency of the deterministic part of the signal
%	    blurried in a white gaussian noise.
%	    L must be an integer       	(default : 1).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                (default : 0).
%	FNORMHAT : Output (normalized) instantaneous frequency.
%	T : Time instants.
%
%	Examples :
%	 x=fmsin(70,0.05,0.35,25); [instf,t]=instfreq(x); plot(t,instf)
%	 N=64; SNR=10.0; L=4; t=L+1:N-L; x=fmsin(N,0.05,0.35,40);
%	 sig=sigmerge(x,hilbert(randn(N,1)),SNR);
%	 plotifl(t,[instfreq(sig,t,L),instfreq(x,t)]); grid;
%	 title ('theoretical and estimated instantaneous frequencies');
%
%	See also  KAYTTH, SGRPDLAY.

%	F. Auger, March 1994, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least one parameter required');
end;
[xrow,xcol] = size(x);
if (xcol~=1),
 error('X must have only one column');
end

if (nargin == 1),
 t=2:xrow-1; L=1; trace=0.0;
elseif (nargin == 2),
 L = 1; trace=0.0;
elseif (nargin == 3),
 trace=0.0;
end;

if L<1,
 error('L must be >=1');
end
[trow,tcol] = size(t);
if (trow~=1),
 error('T must have only one row');
end;

if (L==1),
 if any(t==1)|any(t==xrow),
  error('T can not be equal to 1 neither to the last element of X');
 else
  fnormhat=0.5*(angle(-x(t+1).*conj(x(t-1)))+pi)/(2*pi);
 end;
else
 H=kaytth(L);
 if any(t<=L)|any(t+L>xrow),
  error('The relation L<T<=length(X)-L must be satisfied');
 else
  for icol=1:tcol,
   if trace, disprog(icol,tcol,10); end;
   ti = t(icol); tau = 0:L;
   R = x(ti+tau).*conj(x(ti-tau));
   M4 = R(2:L+1).*conj(R(1:L));

   diff=2e-6;
   tetapred = H * (unwrap(angle(-M4))+pi);
   while tetapred<0.0 , tetapred=tetapred+(2*pi); end;
   while tetapred>2*pi, tetapred=tetapred-(2*pi); end;
   iter = 1;
   while (diff > 1e-6)&(iter<50),
    M4bis=M4 .* exp(-j*2.0*tetapred);
    teta = H * (unwrap(angle(M4bis))+2.0*tetapred);
    while teta<0.0 , teta=(2*pi)+teta; end;
    while teta>2*pi, teta=teta-(2*pi); end;
    diff=abs(teta-tetapred);
    tetapred=teta; iter=iter+1;
   end;
   fnormhat(icol,1)=teta/(2*pi);
  end;
 end;
end;

author: HeLiangHIT


text:
function results = tensor_rpca_adal2( data, params )
% Solve
%       min_{X,E} \sum_i ||X_(i)||_* + \lambda1*||E||_1
%       s.t. X + E = T

% data.T
% params
% X, V are cell arrays of tensors.
%
% Algorithm: ADAL

tic;
T = data.T;
N = length( size(params.X0) );
X = cell( 1, N );
U = cell( 1, N );
V = cell( 1, N );
for i = 1:N
    X{i} = params.X0;
    V{i} = params.V0{i};
end
E = params.E0;

lambda = params.lambda;
mu = params.mu1;
verbose = params.verbose;
tnorm = norm(T);

for iter = 1:params.max_iter
    % solve X_i's
    for i = 1:N
        [X{i}, junk, U{i}] = tensor_shrinkage( T+mu*V{i}-E, mu, i );
    end

    % solve E
    Ep = E;
    D = cell( 1, N );
    for i = 1:N
        D{i} = X{i} - T - mu*V{i};
    end
    P = -(1/N)*ten_sum_all(D);
    E = shrinkage_t( P, lambda*mu/N );

    % compute optimality stats
    pres = 0;
    tdiff = cell( 1, N );
    for i = 1:N
        tdiff{i} = X{i} + E - T;
    end
    pres = tensor_array_norm(tdiff);

    Ediff = E - Ep;
    pres = pres / (sqrt(N)*tnorm);
    dres = norm(Ediff) / norm(Ep);
    Y = ten_sum_all(X) / N;
    rel_err = norm(Y-data.X) / norm(data.X);

    % print
    if verbose
    fprintf('Iter: %d,   pinf: %3.2e,   dinf: %3.2e,    rel_err: %3.2e\n', iter, pres, dres, rel_err );
    end

    if max(pres, dres) < params.opt_tol
%     if pres < params.opt_tol
        break;
    end

    % update Lagrange multipliers
    for i = 1:N
        V{i} = V{i} - tdiff{i}/mu;
    end

%     if rem(iter,10) == 0
%         mu = max(mu*0.9, params.mu_min);
%     end
end

results.X = Y;
results.E = E;
results.V = V;
results.T = T;
results.U = U;
results.iter = iter;
results.cpu = toc;
results.mu = mu;
results.lambda = lambda;

end

author: andrewssobral

repo: lrslibrary


text:
function r8_atan_test ( )

%*****************************************************************************80
%
%% R8_ATAN_TEST tests R8_ATAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  ntest = 8;

  xtest(1:8) = [ 1.0, 1.0, 0.0, -1.0, -1.0, -1.0,  0.0,  1.0 ];
  ytest(1:8) = [ 0.0, 1.0, 1.0,  1.0,  0.0, -1.0, -1.0, -1.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_ATAN_TEST\n' );
  fprintf ( 1, '  R8_ATAN computes the arc-tangent given Y and X;\n' );
  fprintf ( 1, '  ATAN2 is the system version of this routine.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X     Y     ATAN2(Y,X)   R8_ATAN(Y,X)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest
    x = xtest(i);
    y = ytest(i);
    fprintf ( 1, '  %12f  %12f  %12f  %12f\n', ...
      x, y, atan2 ( y, x ), r8_atan ( y, x ) );
  end

  return
end


text:
function  res = partialDCT(n,m,J)

res.adjoint = 0;
res.n = n;
res.m = m;
res.J = J(1:m);

% Register this variable as a partialDCT class
res = class(res,'partialDCT');

author: thomaskuestner


text:
function example2 ( )

%*****************************************************************************80
%
%% EXAMPLE2 uses BVP4C to solve the EXAMPLE2 problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2013
%
%  Author:
%
%    Original MATLAB version by Shampine, Kierzenka, Reichelt.
%    This version by John Burkardt.
%
%  Reference:
%
%    Lawrence Shampine, Jacek Kierzenka, Mark Reichelt,
%    Solving boundary value problems for ordinary differential equations
%    in MATLAB with bvp4c.
%
  global p

  p = 0.00001;

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXAMPLE2:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use BVP4C to solve the following boundary value problem:\n' );
  fprintf ( 1, '  y" + 3 p y / ( p + t^2 )^2 = 0\n' );
  fprintf ( 1, '  y(-0.1) = -0.1 / sqrt ( p + 0.01 )\n' );
  fprintf ( 1, '  y(+0.1) = +0.1 / sqrt ( p + 0.01 )\n' );
  fprintf ( 1, '  Use initial guesses:\n' );
  fprintf ( 1, '  y(x) = 0\n' );
  fprintf ( 1, '  y''(x) = 10\n' );
  fprintf ( 1, '  Use parameter value p = %g\n', p );
%
%  Set SOLINIT, the structure defining the initial guess.
%
  x_init = linspace ( -0.1, 0.1, 11 );
  y_init = [ 0.0, 10.0 ];
  solinit = bvpinit ( x_init, y_init );
%
%  Because this problem is difficult to solve accurately, we must request
%  a more stringent tolerance.  This is done through the optional OPTIONS
%  input.
%
  options = bvpset ( 'RelTol', 1.0E-04 );
%
%  Have BVP4C solve the problem.
%
  sol = bvp4c ( @example2_ode, @example2_bc, solinit, options );
%
%  Use DEVAL to evaluate the solution.
%
  x = linspace ( -0.1, 0.1, 101 );
  sol = deval ( sol, x );

  exact = x ./ sqrt ( p + x.^2 );
%
%  Display a plot of Y and the exact solution.
%
  plot ( x, sol(1,:), 'b*', ...
         x, exact, 'r-', 'Linewidth', 2 );
  xlabel ( '<--- X --->', 'Fontsize', 16 );
  ylabel ( '<--- Y(X) --->', 'Fontsize', 16 );
  title ( 'EXAMPLE2: Exact solution (red), Estimate (blue)', 'Fontsize', 16 )
  grid on
  filename = 'example2.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot file as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXAMPLE2:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function dydx = example2_ode ( x, y )

%*****************************************************************************80
%
%% EXAMPLE2_ODE evaluates the right hand side of the ODE.
%
%  Discussion:
%
%    We assume that the differential equation has been rewritten as a
%    system of first order equations of the form
%
%      dydx = f(x,y)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the point at which the ODE is to be evaluated.
%
%    Input, real Y(M), the value of the solution at X.
%
%    Output, real DYDX(M), the value of the right hand side given X and Y.
%
  global p

  dydx(1) = y(2);
  dydx(2) = - 3.0 * p * y(1) / ( p + x^2 )^2;

  return
end
function bc = example2_bc ( ya, yb )

%*****************************************************************************80
%
%% EXAMPLE2_BC evaluates the boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real YA(M), YB(M), the solution value at the left and right endpoints.
%
%    Output, real BC(M), the value of the boundary conditions.
%
  global p

  bc(1) = ya(1) + 0.1 / sqrt ( p + 0.01 );
  bc(2) = yb(1) - 0.1 / sqrt ( p + 0.01 );

  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:
% LTFAT - Signal processing tools
%
%  Peter L. Søndergaard, 2007 - 2018.
%
%  General
%    RMS            -  Root Mean Square norm of signal.
%    SETNORM      -  Normalize signal by specified norm.
%    GAINDB         -  Scale input signal
%    CRESTFACTOR    -  Compute the crest factor of a signal.
%    UQUANT         -  Simulate uniform quantization.
%
%  Window functions
%    FIRWIN         -  FIR windows (Hanning,Hamming,Blackman,...).
%    FIRKAISER      -  FIR Kaiser-Bessel window.
%    FIR2LONG       -  Extend FIR window to LONG window.
%    LONG2FIR       -  Cut LONG window to FIR window.
%    FREQWIN        -  Freq responses (Gauss,Gammatone,Butterworth,...)
%
%  Wavelet functions
%    FREQWAVELET    -  Frequency responses of wavelets (Cauchy, Morse)
%
%  Filtering
%    FIRFILTER      -  Construct an FIR filter.
%    BLFILTER       -  Construct a band-limited filter.
%    WARPEDBLFILTER -  Warped, band-limited filter.
%    FREQFILTER     -  Construct a full length frequency side filter.
%    PFILT          -  Apply filter with periodic boundary conditions.
%    MAGRESP        -  Magnitude response plot.
%    TRANSFERFUNCTION - Compute the transfer function of a filter.
%    PGRPDELAY      -  Periodic Group Delay
%
%  Ramping
%    RAMPUP         -  Rising ramp.
%    RAMPDOWN       -  Falling ramp.
%    RAMPSIGNAL     -  Ramp a signal.
%
%  Thresholding methods
%    THRESH         -  Coefficient thresholding.
%    LARGESTR       -  Keep largest ratio of coefficients.
%    LARGESTN       -  Keep N largest coefficients.
%    DYNLIMIT       -  Limit the dynamical range.
%    GROUPTHRESH    -  Group thresholding.
%
%  Image processing
%    RGB2JPEG       -  Convert RGB values to the JPEG colour model
%    JPEG2RGB       -  Convert values from the JPEG colour model to RGB
%
%  Tools for OFDM
%    QAM4           -  Quadrature amplitude modulation, order 4
%    IQAM4          -  Inverse QAM of order 4
%
%  Tonal-transient separation
%    tfjigsawsep      - Tonal-transient-residual separation using the T-F jigsaw puzzle algorithm.
%    plottfjigsawsep  - Plot the separated layers.
%
%  For help, bug reports, suggestions etc. please visit
%  http://github.com/ltfat/ltfat/issues
%

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [varargout] = triangleRayIntersection (V_ori,R,V,F,optStruct)

% function [V_intersect,L_intersect,T] = triangleRayIntersection (V_ori,R,V,F,optStruct)
% ------------------------------------------------------------------------
%
% Ray/triangle intersection using the algorithm proposed by Muller and
% Trumbore (1997)
%
% Based on implementation by: Jarek Tuszynski (jaroslaw.w.tuszynski@saic.com)
%
%
% Algorithm:
%  Function solves
%        |t|
%    M * |u| = (o-v0)
%        |v|
%  for [t; u; v] where M = [-d, v1-v0, v2-v0]. u,v are barycentric coordinates
%  and t - the distance from the ray origin in |d| units
%  ray/triangle intersect if u>=0, v>=0 and u+v<=1
%
% Note:
%  The algorithm is able to solve several types of problems:
%  * many faces / single ray  intersection
%  * one  face  / many   rays intersection
%  * one  face  / one    ray  intersection
%  * many faces / many   rays intersection
%  In order to allow that to happen all imput arrays are expected in Nx3
%  format, where N is number of vertices or rays. In most cases number of
%  vertices is different than number of rays, so one of the imputs will
%  have to be cloned to have the right size. Use "repmat(A,size(B,1),1)".
%
% Input (all arrays in in Nx3 format, where N is number of vertices or rays):
%  * orig : ray's origin
%  * dir  : ray's direction
%  * vert0, vert1, vert2: vertices of the triangle
%  * options: aditional customization options
%    * options.triangle - 'one sided' or 'two sided' (default) - how to treat
%        triangles. In 'one sided' version only intersections in single
%        direction are counted and intersections with back facing
%           tringles are ignored
%    * options.ray - 'ray' (default) or 'segment' - how to treat ray as an
%        infinite line (ray) or as line segment defined by a vector
%    * option.border - controls border handling. If 'normal'(default)
%        border points are included, but can be easily lost due to
%        rounding errors. If option.border='inclusive' borders points are
%        included, with a margin of option.eps. If option.border='exclusive'
%        borders points are excluded, with margin of option.eps.
%    * options.epsilon (default = 1e-5)
%
% Output:
%   * Intersect - boolean array of length N
%   * t   - distance from the ray origin to the intersection point in |dir|
%   * u,v - barycentric coordinates of the intersection point units
%
% Based on:
%  *"Fast, minimum storage ray-triangle intersection". Tomas Muller and
%    Ben Trumbore. Journal of Graphics Tools, 2(1):21--28, 1997.
%    http://www.graphics.cornell.edu/pubs/1997/MT97.pdf
%  * http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/raytri/
%  * http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/raytri/raytri.c
%
% Change log:
% 2013/12/18 Changed for patch type input and Cartesian coordinate output
% 2019/07/03 Added varargout type output
% 2021/08/18 Listed as depricated in favour of triSurfRayTrace
% ------------------------------------------------------------------------

%%

warning('This function is depricated, it is recommended to update your code to use triSurfRayTrace instead');

%%

%Get triangle vertices
V0=V(F(:,1),:);
V1=V(F(:,2),:);
V2=V(F(:,3),:);

%% verify that inputs are in correct format

if (size(V_ori,1)==3 && size(V_ori,2)~=3)
    V_ori =V_ori';
end

if (size(R,1)==3 && size(R,2)~=3)
    R =R';
end

if (size(V0,1)==3 && size(V0,2)~=3)
    V0=V0';
end

if (size(V1,1)==3 && size(V1,2)~=3)
    V1=V1';
end

if (size(V2,1)==3 && size(V2,2)~=3)
    V2=V2';
end

if (any(size(V_ori)~=size(V0)) || ...
        any(size(V_ori)~=size(V1)) || ...
        any(size(V_ori)~=size(V2)) || ...
        any(size(V_ori)~=size(R  )) )
    error('All input vectors have to be of the same size.');
end

if (size(V_ori,2)~=3)
    error('All input vectors have to be in Nx3 format.');
end

%% Read options from optStruct

%Create defaults
if (nargin<5)
    optStruct.eps      = 1e-5;
    optStruct.triangle = 'two sided';
    optStruct.ray      = 'ray';
    optStruct.border   = 'normal';
end

%% initialize default output
L_intersect = false(size(V_ori,1),1);
T=nan(size(V_ori,1),1);
V_bar=T;
V_intersect=nan(size(V_ori));

%% Find faces parallel to the ray
E1 = V1-V0;          % find vectors for two edges sharing vert0
E2 = V2-V0;
Tvec  = V_ori -V0;          % distance from vert0 to ray origin
Pvec  = cross(R, E2,2);  % begin calculating determinant - also used to calculate U parameter

detMat   = sum(E1.*Pvec,2);   % determinant of the matrix M = dot(E1,Pvec)

logicParalel = (abs(detMat)<optStruct.eps);    % if determinant is near zero then ray lies in the plane of the triangle

if all(logicParalel) % if all parallel than no intersections
    varargout{1}=V_intersect;
    varargout{2}=L_intersect;
    varargout{3}=T;
    return
end

switch optStruct.border
    case 'normal'
        zeroLim=0.0;
    case 'inclusive'
        zeroLim=optStruct.eps;
    case 'exclusive'
        zeroLim=-optStruct.eps;
end

%% Different behavior depending on one or two sided triangles
if strcmpi(optStruct.triangle,'two sided')          % treats triangles as two sided

    detMat(logicParalel) = 1;                       % change to avoid division by zero

    U_bar  = sum(Tvec.*Pvec,2)./detMat;             % calculate U parameter used to test bounds
    L_ok = (~logicParalel & U_bar>=-zeroLim & U_bar<=1.0+zeroLim);% mask which allows performing next 2 operations only when needed

    if ~any(L_ok) % if all ray/plane intersections are outside the triangle than no intersections
        varargout{1}=V_intersect;
        varargout{2}=L_intersect;
        varargout{3}=T;
        return
    end

    Qvec = cross(Tvec(L_ok,:), E1(L_ok,:),2); % prepare to test V parameter
    V_bar(L_ok,:) = sum(R(L_ok,:).*Qvec,2) ./ detMat(L_ok,:);  % calculate V parameter used to test bounds
    L_intersect = (V_bar>=-zeroLim & U_bar+V_bar<=1.0+zeroLim & L_ok);

    T(L_ok,:) = sum(E2(L_ok,:).*Qvec,2)./detMat(L_ok,:);

    if ~(strcmpi(optStruct.ray,'ray'))
        L_intersect = (L_intersect & T>=-zeroLim & T<=1.0+zeroLim);
    end

else % treats triangles as one sided

    U_bar = sum(Tvec.*Pvec,2);                   % calculate U parameter used to test bounds
    L_ok = (detMat>optStruct.eps & U_bar>=0.0 & U_bar<=detMat);        % mask which allows performing next 2 operations only when needed

    if ~any(L_ok) % if all ray/plane intersections are outside the triangle than no intersections
        varargout{1}=V_intersect;
        varargout{2}=L_intersect;
        varargout{3}=T;
        return
    end

    Qvec = cross(Tvec(L_ok,:), E1(L_ok,:),2); % prepare to test V parameter
    V_bar(L_ok,:) = sum(R(L_ok,:).*Qvec,2);        % calculate V parameter used to test bounds
    L_intersect = (detMat>optStruct.eps & U_bar>=-zeroLim & V_bar>=-zeroLim & U_bar+V_bar<=detMat*(1+zeroLim));

    T(L_ok,:)  = sum(E2(L_ok,:).*Qvec,2);
    inv_det = zeros(size(detMat));
    inv_det(L_ok,:) = 1./detMat(L_ok,:);
    T = T.*inv_det;  % calculate t - distance from origin to the intersection in |d| units
    U_bar = U_bar.*inv_det;
    V_bar = V_bar.*inv_det;

    if ~(strcmpi(optStruct.ray,'ray'))
        L_intersect = (L_intersect & T>=-zeroLim & T<=1.0+zeroLim); % intersection between origin and destination
    end
end

V_intersect=((1-U_bar(:,ones(1,3))-V_bar(:,ones(1,3))).*V0) + (U_bar(:,ones(1,3)).*V1) + (V_bar(:,ones(1,3)).*V2);
V_intersect(~L_intersect,:)=NaN;
T(~L_intersect,:)=NaN;

%% Collect output
varargout{1}=V_intersect;
varargout{2}=L_intersect;
varargout{3}=T;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function data_num = p07_data_num ( )

%*****************************************************************************80
%
%% P07_DATA_NUM returns the dimension of the data vector for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DATA_NUM, the dimension of the data vector.
%
  data_num = 9;

  return
end


text:
function h = dtiImportFibers(h, filename, xform)
%
% handles = dtiImportFibers(handles, filename)
%
% Imports non-mrDiffusion fiber/path file formats. Currently supports:
%
% 1. MetroTrac (*.pdb) format
% 2. Camino Raw (*.Bfloat) format
%
%
% HISTORY:
% 2006.08.17 RFD: wrote it.
%

fg = mtrImportFibers(filename, xform);
if ~isempty(fg)
    h = dtiAddFG(fg, h);
end
return;


text:
function g = disk_grid_fibonacci ( n, r, c )

%*****************************************************************************80
%
%% DISK_GRID_FIBONACCI computes Fibonacci grid points inside a disk.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Swinbank, James Purser,
%    Fibonacci grids: A novel approach to global modelling,
%    Quarterly Journal of the Royal Meteorological Society,
%    Volume 132, Number 619, July 2006 Part B, pages 1769-1793.
%
%  Parameters:
%
%    Input, integer N, the number of points desired.
%
%    Input, real R, the radius of the disk.
%
%    Input, real C(2), the coordinates of the center of the disk.
%
%    Output, real G(2,N), the grid points.
%
  r0 = r / sqrt ( n - 0.5 );
  phi = ( 1.0 + sqrt ( 5.0 ) ) / 2.0;

  i = 1 : n;

  gr = r0 * sqrt ( i - 0.5 );
  gt = 2.0 * pi * i / phi;

  g = zeros ( 2, n );

  g(1,1:n) = c(1);
  g(1,1:n) = g(1,1:n) + gr .* cos ( gt );
  g(2,1:n) = c(2);
  g(2,1:n) = g(2,1:n) + gr .* sin ( gt );

  return
end


text:
function [xi,w]=fifthOrderHexagonCubPoints()
%%FIFTHORDERHEXAGONCUBPOINTS Generate fifth-order cubature points for
%               integration over a two-dimensional hexagon with vertices
%               (+/-1,0) and (+/-(1/2),+/-(1/2)*sqrt(3)).
%
%INPUTS: None
%
%OUTPUTS: xi A numDim X numCubaturePoints matrix containing the cubature
%            points (Each "point" is a vector).
%          w A numCubaturePoints X 1 vector of the weights associated with
%            the cubature points.
%
%This function implements Formula H_2 5-1 in [1], pg. 336, 7 points.
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

r=sqrt(14/25);
s=sqrt(14/100);
t=sqrt(42/100);

xi=[[0;0],PMCombos([r;0]),PMCombos([s;t])];

V=3*sqrt(3)/2;
A=(258/1008)*V;
B=(125/1008)*V;
w=[A;B*ones(6,1)];
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [Ax,Ay,b,K] = convexhullFromSampled(z,fz,xL,xU);

k1 = max((fz(2:end)-fz(1))./(z(2:end)-xL))+1e-12;
k2 = min((fz(2:end)-fz(1))./(z(2:end)-xL))-1e-12;
k3 = min((fz(1:end-1)-fz(end))./(z(1:end-1)-xU))+1e-12;
k4 = max((fz(1:end-1)-fz(end))./(z(1:end-1)-xU))-1e-12;
Ax = [-k1;k2;-k3;k4];
Ay = [1;-1;1;-1];
b =  [k1*(-z(1)) + fz(1);-(k2*(-z(1)) + fz(1));k3*(-z(end)) + fz(end);-(k4*(-z(end)) + fz(end))];
K.f = 0;
K.l = length(b);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [Hx,Hy,Ez,time] = MaxwellHNonCon2D(Hx, Hy, Ez, FinalTime)

% function [Hx,Hy,Ez] = MaxwellHNonCon2D(Hx, Hy, Ez, FinalTime)
% Purpose  : Integrate TM-mode Maxwell's until FinalTime starting with initial conditions Hx,Hy,Ez
%

Globals2D;
time = 0;

% Runge-Kutta residual storage
resHx = zeros(Np,K); resHy = zeros(Np,K); resEz = zeros(Np,K);

% compute time step size
rLGL = JacobiGQ(0,0,N); rmin = abs(rLGL(1)-rLGL(2));
dtscale = dtscale2D; dt = min(dtscale)*rmin*2/3;

% find non-conforming neighbors
neighbors = BuildHNonCon2D(N+1, 1e-6);

% outer time step loop
tstep = 1;
while (time<FinalTime)

  if(time+dt>FinalTime), dt = FinalTime-time; end

   for INTRK = 1:5
      % compute right hand side of TM-mode Maxwell's equations
      [rhsHx, rhsHy, rhsEz] = MaxwellHNonConRHS2D(Hx,Hy,Ez, neighbors);

      % initiate and increment Runge-Kutta residuals
      resHx = rk4a(INTRK)*resHx + dt*rhsHx;
      resHy = rk4a(INTRK)*resHy + dt*rhsHy;
      resEz = rk4a(INTRK)*resEz + dt*rhsEz;

      % update fields
      Hx = Hx+rk4b(INTRK)*resHx;
      Hy = Hy+rk4b(INTRK)*resHy;
      Ez = Ez+rk4b(INTRK)*resEz;
   end;

   % Increment time
   time = time+dt; tstep = tstep+1;

   if(~mod(tstep, 5))
     PlotField2D(N, x, y, Ez);
     hold on; ha = trimesh(EToV, VX, VY); set(ha, 'Color', 'black'); hold off
     axis([-1 1 -1 1 -1 1]); drawnow; pause(.02);
   end

end
return

repo: Aero-matlab


text:
function val=erfcxComplex(z)
%%ERFCXCOMPLEX Evaluate the scaled complementary error function in a manner
%             valid for complex arguments. The complementary error function
%             is one minus the error function (erf). The scaled
%             complementary error function is exp(z^2)*erfc(z). The
%             functions erfcx and erfc that are built into Matlab cannot
%             handle complex arguments.
%
%INPUTS: z A scalar, vector, or matrix of values at which one wishes to
%          evaluate the scaled complementary error function.
%
%OUTPUTS: val A matrix having the same dimensions as z in which the values
%            of the scaled complementary error function are computed for
%            the values in z.
%
%The scaled complementary error function is just Faddeeva(1i*z) based on
%the definition of the Faddeeva function.
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

val=Faddeeva(1i*z);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
function [g,NR,SI,TI] = regiongrow(f,S,T)
%REGIONGROW Image segmentation using region growing.
%   [G,NR,SI,TI] = REGIONGROW(F,S,T). S can be a matrix of the same size
%   as the input image F, with a 1 at the coordinates of every seed
%   point and 0s elsewhere. S can also be a single seed value. Similarly
%   T can be a matrix of the same size as F, containing a threshold
%   value for each pixel in F. T can also be a scalar, in which case it
%   becomes a global threshold. All values in S and T must be in the
%   range [0,1].
%
%   G is the segmented image resulting from region growing, with each
%   region labeled by a different integer, NR is the number of regions,
%   SI is the final seed image used by the algorithm, and TI is the
%   image consisting of the pixels in F that satisfied the threshold
%   test, but before they were processed for connectivity.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

f = tofloat(f);
% If S is a scalar, obtain the seed image.
if numel(S) == 1
   SI = f == S;
   S1 = S;
else
   % S is a matrix. Eliminate duplicate, connected seed locations to
   % reduce the number of loop executions in the following section of
   % code.
   SI = bwmorph(S,'shrink',Inf);
   S1 = f(SI); % Matrix of seed values.
end

TI = false(size(f));
for K = 1:length(S1)
   seedvalue = S1(K);
   S = abs(f - seedvalue) <= T; % Re-use variable S.
   TI = TI | S;
end

% Use function imreconstruct with SI as the marker image to obtain the
% regions corresponding to each seed in S. Function bwlabel assigns a
% different integer to each connected region.
[g,NR] = bwlabel(imreconstruct(SI,TI));

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
function [lp,dlp] = priorTransform(g,dg,ig,prior,x)

% Transformed univariate hyperparameter prior distribution.
% Compute log-likelihood and its derivative or draw a random sample.
% The prior distribution is parameterized as:
%
%  prior(y), where y = g(x),
%
% g is the transformation function (monotonically increasing and invertable),
% dg is its derivative, ig is its inverse, prior is the distribution to be
% transformed and x(1xN) contains query hyperparameters for prior evaluation.
% Note that p(x) = prior(g(x)) is not necessarily normalised w.r.t. x.
%
% For example, to put a nonnegative gamma prior on the lenthscale
% parameter ell of covSEiso, one can use
%   {@priorTransform,@exp,@exp,@log,{@priorGamma,k,t}}
% since ell is represented in the log domain.
%
% For more help on design of priors, try "help priorDistributions".
%
% Copyright (c) by Roman Garnett and Hannes Nickisch, 2014-09-08.
%
% See also PRIORDISTRIBUTIONS.M.

if nargin<4, error('g, dg, ig and prior parameters need to be provided'), end
if nargin<5, lp = ig(feval(prior{:})); return, end      % apply inverse sampling

if ~iscell(prior), prior = {prior}; end
[lp,dlp] = feval(prior{:},g(x)); dlp = dlp.*dg(x);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ROImask] = selectROI(IM, varargin)
%% function for selecting the ROI on speckle images in step2
% options:
%[] = selectROI(IM);
%[] = selectROI(IM,nROI);
%
% function for selecting the region of interest (ROI) for stereo DIC (2 views)
% requirements: GIBBON toolbox
% The idea is to draw a polygon around the ROI, which is the overlapping
% area that is visible in both views. The ROI is selected on the reference
% image, but in order to select it properly, it is advisable to look at
% both the reference and the "deformed" views.
%
% INPUT:
% * IMset - a 2nX1 cell array containing paths 2n grayscale images. The first n
% images are from camera A (the "reference" camera), and the last n images
% are from camera B (the "deformed" camera). The first image in the set is
% considered as the reference image, on which the reference grid of points
% is defined, and all the correlated points and consequent displacements
% and strains, are relative to this image.
% * nROI (optional, default=1) - number of ROIs (closed polygons). must be
% integer>=1; If nROI is empty [] or not given then it is set to the default(=1).

%%
nVars = length(varargin);

switch nVars
    case 0
        nROI=1;
    case 1
        nROI=varargin{1};
        if isempty(nROI), nROI=1; end
    otherwise
        error('Wrong number of input arguments');
end

polyOK='No';
while ~strcmp(polyOK,'Yes')
    hf=figure('Name','Select ROI');
    hf.Units='normalized'; hf.OuterPosition=[.05 .05 .9 .9]; hf.Units='pixels';

    gtitle(['Draw ' num2str(nROI) ' ROI polygons on the reference images. double-click on the polygon to finish'],20);

    imshow(IM); hold on; axis on;
%     title('Reference image');
    axis off

    for iROI=1:nROI
        hp(iROI) = impoly;
        wait(hp(iROI));
    end
    polyOK = questdlg('polygon mask OK? (if you select "No", this GUI will start over)', 'polygon mask OK?', 'Yes', 'No', 'Yes');
end

for iROI=1:nROI
    nodes = getPosition(hp(iROI));
    ROImask(:,:,iROI) = poly2mask(nodes(:,1),nodes(:,2),size(IM,1),size(IM,2));
end

ROImask=logical(sum(ROImask,3));

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function fftfilegui(varargin)
%FFTFILEGUI Read image, generate, and view data in the dispersed domain using MATLAB GUIs.
%
%    fftfile('PropertyName',PropertyValue,...)
%
% Function will query user for input file and AOI.  Only real parameter
% that can be passed in dimension for FFT to operate.  Default is 2D FFT.
%
%       Property name     Description
%       dims              dimensions along which to compute FFT (options:
%                            1, 2, or [1 2] (2-D FFT), default = [1 2])
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

demo_core(@fftfile, varargin{:});

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [bit] = MB2bit(MB)
% Convert computery things from megabytes to bits.
% Chad A. Greene 2012
bit = MB*8388608 ;

author: Sable

repo: mcbench-benchmarks


text:
clear all
clc

%load vgg19
net = load('./models/imagenet-vgg-verydeep-19.mat');
net = vl_simplenn_tidy(net);

% testing dataset
test_path_ir = './IV_images/IR/';
fileFolder_ir=fullfile(test_path_ir);
dirOutput_ir =dir(fullfile(fileFolder_ir,'*'));
num_ir = length(dirOutput_ir);

test_path_vis = replace(test_path_ir, '/IR/', '/VIS/');
fileFolder_vis=fullfile(test_path_vis);
dirOutput_vis =dir(fullfile(fileFolder_vis,'*'));

for i=3:num_ir
index = i;
disp(num2str(index));

path1 = [test_path_ir,dirOutput_ir(i).name]; % IR image
path2 = [test_path_vis,dirOutput_vis(i).name]; % VIS image
fused_path = ['./fused_iv/fused_vggml_',dirOutput_ir(i).name];

image1 = imread(path1);
image2 = imread(path2);
image1 = im2double(image1);
image2 = im2double(image2);

[h,w,d] = size(image1);
isRe = 0;
if d > 1
    image1 = rgb2gray(image1);
    image2 = rgb2gray(image2);
end
if h < 224
    isRe = 1;
    image1 = imresize(image1, [224, w]);
    image2 = imresize(image2, [224, w]);
end
if w < 224
    isRe = 1;
    image1 = imresize(image1, [h, 224]);
    image2 = imresize(image2, [h, 224]);
end

tic;
% Highpass filter test image
npd = 16;
fltlmbd = 5;
[I_lrr1, I_saliency1] = lowpass(image1, fltlmbd, npd);
[I_lrr2, I_saliency2] = lowpass(image2, fltlmbd, npd);

%% fuison lrr parts
F_lrr = (I_lrr1+I_lrr2)/2;
% figure;imshow(F_lrr);

%% fuison saliency parts use VGG19
% disp('VGG19-saliency');
saliency_a = make_3c(I_saliency1);
saliency_b = make_3c(I_saliency2);
saliency_a = single(saliency_a) ; % note: 255 range
saliency_b = single(saliency_b) ; % note: 255 range

res_a = vl_simplenn(net, saliency_a);
res_b = vl_simplenn(net, saliency_b);

%% relu1_1
% disp('relu1_1');
out_relu1_1_a = res_a(2).x;
out_relu1_1_b = res_b(2).x;
unit_relu1_1 = 1;

l1_featrues_relu1_a = extract_l1_feature(out_relu1_1_a);
l1_featrues_relu1_b = extract_l1_feature(out_relu1_1_b);
% average
[m,n,k] = size(out_relu1_1_a);

[F_saliency_relu1, l1_featrues_relu1_ave_a, l1_featrues_relu1_ave_b] = ...
            fusion_strategy(l1_featrues_relu1_a, l1_featrues_relu1_b, I_saliency1, I_saliency2, unit_relu1_1);

%% relu2_1
% disp('relu2_1');
out_relu2_1_a = res_a(7).x;
out_relu2_1_b = res_b(7).x;
unit_relu2_1 = 2;

l1_featrues_relu2_a = extract_l1_feature(out_relu2_1_a);
l1_featrues_relu2_b = extract_l1_feature(out_relu2_1_b);

[F_saliency_relu2, l1_featrues_relu2_ave_a, l1_featrues_relu2_ave_b] = ...
            fusion_strategy(l1_featrues_relu2_a, l1_featrues_relu2_b, I_saliency1, I_saliency2, unit_relu2_1);

%% relu3_1
% disp('relu3_1');
out_relu3_1_a = res_a(12).x;
out_relu3_1_b = res_b(12).x;
unit_relu3_1 = 4;

l1_featrues_relu3_a = extract_l1_feature(out_relu3_1_a);
l1_featrues_relu3_b = extract_l1_feature(out_relu3_1_b);

[F_saliency_relu3, l1_featrues_relu3_ave_a, l1_featrues_relu3_ave_b] = ...
            fusion_strategy(l1_featrues_relu3_a, l1_featrues_relu3_b, I_saliency1, I_saliency2, unit_relu3_1);

%% relu4_1
% disp('relu4_1');
out_relu4_1_a = res_a(21).x;
out_relu4_1_b = res_b(21).x;
unit_relu4_1 = 8;

l1_featrues_relu4_a = extract_l1_feature(out_relu4_1_a);
l1_featrues_relu4_b = extract_l1_feature(out_relu4_1_b);

[F_saliency_relu4, l1_featrues_relu4_ave_a, l1_featrues_relu4_ave_b] = ...
            fusion_strategy(l1_featrues_relu4_a, l1_featrues_relu4_b, I_saliency1, I_saliency2, unit_relu4_1);

%% fusion strategy
% disp('choose max');
F_saliency = max(F_saliency_relu1, F_saliency_relu2);
F_saliency = max(F_saliency, F_saliency_relu3);
F_saliency = max(F_saliency, F_saliency_relu4);
% figure;imshow(F_saliency);

% imwrite(F_saliency,'./decomposition/F_s.png','png');

fusion_im = F_lrr + F_saliency;
toc;
% figure;imshow(fusion_im);
if isRe==1
    fusion_im = imresize(fusion_im, [h, w]);
end
imwrite(fusion_im,fused_path,'png');
end

author: thfylsty


text:
% SYNTAX:
% [stim, tRange] = hmrR_StimRejection(data, stim, tIncAuto, tIncMan, tRange)
%
% UI NAME:
% Stim_Exclude
%
% DESCRIPTION:
% Excludes stims that fall within the time points identified as
% motion artifacts from HRF calculation.
%
% INPUT:
% data:     SNIRF data object
% stim:     SNIRF stim object
% tIncAuto: Cell array of length equal to the # of time bases in data. Each
%           cell element is time points (#time points x 1) identified as motion
%           artifacts by processing stream.
% tIncMan:  Cell array of length equal to the # of time bases in data. Each
%           cell element is time points (#time points x 1) identified as motion
%           artifacts by user.
% tRange:   an array of 2 numbers [t1 t2] specifying how many
%           seconds surrounding motion artifacts, tIncMan and tIncAuto,
%           to consider as excluded data and therefore exclude any stimuli
%           which fall within those buffers.
%           Typically values are t1=-2 and t2 equal to the stimulus
%           duration.
%
% OUTPUT:
% stim:     SNIRF stim object
% tRange:   same tRange array as in the input
%
% USAGE OPTIONS:
% Stim_Exclude: [stim,tRange] = hmrR_StimRejection(dod,stim,tIncAuto,tIncMan,tRange)
%
% PARAMETERS:
% tRange: [-5.0, 10.0]
%
% PREREQUISITES:
% Motion_Artifact: tIncAuto = hmrR_MotionArtifact(dod, probe, mlActMan, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh), Motion_Artifact_By_Channel: [tIncAuto, tIncAutoCh] = hmrR_MotionArtifactByChannel(dod, probe, mlActMan, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)
%

function [stim, tRange] = hmrR_StimRejection(data, stim, tIncAuto, tIncMan, tRange)

if isempty(tIncAuto)
    tIncAuto = cell(length(data),1);
end
if isempty(tIncMan)
    tIncMan = cell(length(data),1);
end

snirf = SnirfClass(data, stim);
for iBlk=1:length(snirf.data)
    t = snirf.data(iBlk).GetTime();
    s = snirf.GetStims(t);

    dt = (t(end)-t(1))/length(t);
    tRangeIdx = [floor(tRange(1)/dt):ceil(tRange(2)/dt)];

    smax = max(s,[],2);
    lstS = find(smax==1);
    for iS = 1:length(lstS)
        lst = round(min(max(lstS(iS) + tRangeIdx,1),length(t)));
        if ~isempty(tIncAuto{iBlk}) && min(tIncAuto{iBlk}(lst))==0
            s(lstS(iS),:) = -2*abs(s(lstS(iS),:));
        end
        if ~isempty(tIncMan{iBlk}) && min(tIncMan{iBlk}(lst))==0
            s(lstS(iS),:) = -1*abs(s(lstS(iS),:));
        end
    end
    snirf.SetStims_MatInput(s, t);
end
stim = snirf.stim;

author: BUNPC


text:
function env = envelopeDetection(x)
%ENVELOPEDETECTION  Extract signal envelope using the Hilbert Transform.
%
% DESCRIPTION:
%       envelopeDetection applies the Hilbert transform to extract the
%       envelope from an input vector x. If x is a matrix, the envelope
%       along each row is returned.
%
%       Example:
%           x = toneBurst(10e6, 0.5e6, 10);
%           plot(0:length(x)-1, x, 'k-', 0:length(x)-1, envelopeDetection(x), 'r-');
%           legend('Input Signal', 'Envelope');
%
% USAGE:
%       env = envelopeDetection(x)
%
% INPUTS:
%       x           - input function
%
% OUTPUTS:
%       env         - envelope of input function
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 20th December 2010
%       last update - 14th November 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% if the input is one-dimensional, force to be 1 x Nx
if numDim(x) == 1
    x = reshape(x, 1, []);
end

% compute the FFT of the input function (use zero padding to prevent
% effects of wrapping at the beginning of the envelope), if x is a matrix
% the fft's are computed across each row
X = fft(x, length(x(1, :))*2, 2);

% multiply the fft by -1i
X = -1i*X;

% set the DC frequency to zero
X(1) = 0;

% calculate where the negative frequencies start in the FFT
neg_f_index = ceil(length(X(1, :))/2) + 1;

% multiply the negative frequency components by -1
X(:, neg_f_index:end) = -1 * X(:, neg_f_index:end);

% compute the Hilbert transform using the inverse fft
z = ifft(X, [], 2);

% extract the envelope
env = abs(x + 1i*z(:, 1:length(x)));


text:
% POP_RUNSCRIPT - Run Matlab script
%
% Usage: >> pop_runscript;
%        >> pop_runscript( filename );
%
% Input:
%   filename - [string] name of the file.
%
% Author: Arnaud Delorme, SCCN / INC / UCSD, August 2009

% Copyright (C) Arnaud Delorme, August 2009
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function com = pop_runscript(filename);

com = [];
if nargin <1
    [filename filepath] = uigetfile('*.*', 'Please select input script -- pop_runscript()');

    if filename(1) == 0, return; end

    filename = fullfile(filepath, filename);
end

str = readtxtfile(filename);
try
    evalin('base', str);
catch
    lasterr
end
com = sprintf('pop_runscript(''%s'');', filename);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
##% Script to run a MARRMoT model through the BMI interface
##
##%% 1. Create config file
##
##% Get some data
##load MARRMoT_example_data.mat
##
##% Catchment settings
##% -------------------------------------------------------------------------
##
##% 'data_origin' is a 2x1 array with [lat,lon]
##data_origin = [35.29,87.49];
##
##% 'forcing' is a structure with fields 'precip, 'temp', 'pet',
##% 'delta_t_days', 'time_unit'
##forcing.precip       = data_MARRMoT_examples.precipitation;
##forcing.temp         = data_MARRMoT_examples.temperature;
##forcing.pet          = data_MARRMoT_examples.potential_evapotranspiration;
##forcing.delta_t_days = 1;   % 1 [d]
##forcing.time_unit    = 'day';
##
##% 'time_start' & 'time_end' are a vector with starting/end date
##time_start = datevec(data_MARRMoT_examples.dates_as_datenum(1));
##time_end   = datevec(data_MARRMoT_examples.dates_as_datenum(end));
##
##% Model settings
##% -------------------------------------------------------------------------
##
##%'model_name' is a string with the name of the model function
##model_name = 'm_01_collie1_1p_1s';
##
##% 'parameters' is a vector of length n (number of parameters)
##parameters = 10; % Smax, [mm]
##
##% 'store_ini' is a vector of length m (number of stores)
##store_ini = 5;
##
##% 'solver' is a structure with fields 'name', 'resnorm_tolerance',
##% 'resnorm_maxiter'
##solver.name = 'createOdeApprox_IE';
##solver.resnorm_tolerance = 0.1;
##solver.resnorm_maxiter = 6;
##
##% Save as a config file
##% -------------------------------------------------------------------------
##clear data_MARRMoT_examples
##save BMI_testcase_m01_BuffaloRiver_TN_USA
##
##
##%% 2. Run the model through BMI
##clear all % start with a clean slate

addpath(genpath('../..'))

% Specify config path
filepath = '/Users/rwhut/Documents/eWaterCycle/repos/MARRMoT/BMI/Config/BMI_testcase_m01_BuffaloRiver_TN_USA';

% Create a model object
obj_marrmot_m01 = marrmotBMI_oct();

% Initialize the model
obj_marrmot_m01.initialize(filepath)

% Get the number of time steps
ts = datenum(obj_marrmot_m01.startTime);
te = datenum(obj_marrmot_m01.endTime);
tn = te-ts;

% Run a loop
qsim = NaN.*zeros(tn,1);

while obj_marrmot_m01.time <= tn

    % Display time step
    disp(obj_marrmot_m01.time)

    % Advance 1 time step
    obj_marrmot_m01.update();

    % Get simulated flow
    tmp = obj_marrmot_m01.get_value('flux_out'); % gives a structure
    qsim(obj_marrmot_m01.time-1) = tmp.Q;

end

%% 3. Do some visualisation

% get observations
load '/Users/rwhut/Documents/eWaterCycle/repos/MARRMoT/BMI/Config/MARRMoT_example_data.mat'

% Figure
figure('color','w');
hold on;
plot(data_MARRMoT_examples.streamflow);
plot(qsim)
title([obj_marrmot_m01.model_name, ' at ', num2str(obj_marrmot_m01.lat), ...
    ' lat, ',num2str(obj_marrmot_m01.lon), ' lon'],'interpreter','none')
xlabel('time [d]')
ylabel('discharge [mm/d]')

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [J, grad] = costFunctionReg(theta, X, y, lambda)
%COSTFUNCTIONREG Compute cost and gradient for logistic regression with regularization
%   J = COSTFUNCTIONREG(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
h=sigmoid(X*theta);
theta(1,1)=0;
J=sum(-y'*log(h)-(1-y)'*log(1-h))/m+lambda/2/m*sum(power(theta,2));
grad=((h-y)'*X)/m+lambda/m*theta';

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

end

author: loserChen

repo: Coursera-MachineLearning


text:
function vl_help(topic, varargin)
% VL_HELP   VLFeat toolbox builtin help
%   VL_HELP('doc') displays the HTML documentation in MATLAB
%   web-browser. VL_HELP('vl_function') displays the HTML
%   documentation relative to the specified funtion. VL_HELP('topic')
%   displays the HTML documentation relative to a particular 'topic'.
%
%   VL_HELP(..., 'ONLINE', true) uses the online documentation instead of
%   the local copy.
%
%   Example::
%     To get information about the library, a topic such as Fisher vectors,
%     or a function such as VL_FISHER(), use:
%        vl_help doc
%        vl_help fisher
%        vl_help vl_fisher
%
%   See also the core VLFeat
%   documentation.

% Author: Andrea Vedaldi

% Copyright (C) 2013 Andrea Vedaldi
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if nargin == 0
  help vl_help ;
  return ;
end

opts.online = false ;
opts = vl_argparse(opts, varargin) ;

local = fullfile(vl_root, 'doc') ;
if opts.online
  prefix = 'http://www.vlfeat.org/doc' ;
else
  prefix = local ;
end

if strcmp(topic, 'doc')
  topic = 'home' ;
end

[ids,urls] = textread(fullfile(local, 'index.txt'),'%s%s','delimiter','|') ;
k = find(strcmp(topic,ids)) ;
if numel(k) > 0
  web(fullfile(prefix, urls{k(1)}))
else
  error('''%s'' is not a valid function or topic.', topic) ;
end

author: yihui-he

repo: panorama


text:
function neqn = p09_neqn ( )

%*****************************************************************************80
%
%% P09_NEQN returns the number of equations for problem P09.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 3;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function X = mrdivide(A, B)
%/   Right matrix divide for SINGFUN objects.
%
%   Since SINGFUN does not support array-valued objects, right matrix divide is
%   exactly same as right divide.
%
% See also RDIVIDE, LDIVIDE, MLDIVIDE.

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

X = rdivide(A, B);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun