text:
function n = ndims(t)
%NDIMS Return the number of dimensions for a ttensor.
%
%   NDIMS(T) returns the number of dimensions of tensor T.
%
%   See also TTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

n = numel(t.u);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function classifieroutput = fcnnclassify(fcnn,X,R)
%FCNNCLASSIFY Fully-connected neural network classifier.
%   CLASSIFIEROUTPUT = FCNNCLASSIFY(FCNN,X,R) uses a feedforward pass
%   through a fully-connected neural net, fcnn, to classify a set of
%   input pattern vectors (rows of X). If the pattern class membership
%   matrix, R, is provided, this function also outputs the percent of
%   patterns classified correctly.
%
%   Type
%
%   >> help fcninfo
%
%   at the prompt for detailed explanations of the components of the
%   fully-connected neural net.
%
%   classifieroutput is a structure with the following fields.
%
%   classifieroutput.Class
%    A vector whose number of elements is equal to the number of input
%    patterns. The kth element of this vector gives the number of the
%    class to which the kth vector was assigned (i.e., classified).
%
%   classifieroutput.ClassificationRate
%    A scalar that gives the percent of patterns classified correctly,
%    assuming that classifierinput.R was provided. If this is not the
%    case then classifieroutput.ClassificationRate = [].
%
%   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

% INPUTS
np = size(X,2);

% CLASSIFY PATTERNS
% Feedforward.
fcnn = fcnnff(fcnn,X);
% Find the maximum value of fcnn(L).A in each column. The location of a
% maximum in a column gives the class of the pattern corresponding to
% that column. Could also do with softmax to associate a probability
% number with each output, but the net calssificationresult would be the
% same because the largest values of each approach would still
% correspond to the same input patterns.
for j = 1:np
    idx = find(fcnn(end).A(:,j)== max(max((fcnn(end).A(:,j)))));
    classifieroutput.Class(j) = idx;
end

% IF A CLASS MEMBERSHIP MATRIX WAS PROVIDED, USE IT TO COMPUTE THE
% CORRECT CLASSIFICATION RATE.
if nargin == 3
   numErrors = 0;
   % Compute the correct classification rate.
   for j = 1:np
      idx = find(R(:,j) == 1); % Class of input pattern j.
      if ~isequal(classifieroutput.Class(j),idx)
          numErrors = numErrors + 1;
      end
   end
    classifieroutput.ClassificationRate = ((np - numErrors)/np)*100;
else
    classifieroutput.ClassificationRate = [];
end

author: dipum

repo: dipum-toolbox

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


text:
function d = DistanceToSO3(mat)
% *** CHAPTER 3: RIGID-BODY MOTIONS ***
% Takes mat: A 3x3 matrix.
% Returns the Frobenius norm to describe the distance of mat from the SO(3)
% manifold.
% Computes the distance from R to the SO(3) manifold using the following
% method:
% If det(mat) <= 0, return a large number.
% If det(mat) > 0, return norm(mat' * mat - I).
% Example Inputs:
%
% clear; clc;
% mat = [1.0, 0.0,   0.0;
%        0.0, 0.1, -0.95;
%        0.0, 1.0,   0.1];
% d = DistanceToSO3(mat)
%
% Output:
% d =
%     0.0884

if det(mat) > 0
	d = norm(mat' * mat - eye(3), 'fro');
else
    d = 1e+9;
end
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

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

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE0 - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dr = Dice_Ratio(SEG, GT)
    % SEG, GT are the binary segmentation and ground truth areas, respectively.
    % dice ratio
    dr = 2*double(sum(uint8(SEG(:) & GT(:)))) / double(sum(uint8(SEG(:))) + sum(uint8(GT(:))));
end

repo: cluster


text:
function test_bug2303

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_read_header read_eeglabheader

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2303/Dopa2_Quart2_Av.set');

hdr = ft_read_header(filename);

% hdr =
%
%              Fs: 500
%          nChans: 14
%        nSamples: 450
%     nSamplesPre: 100
%         nTrials: 1
%           label: {14x1 cell}
%            elec: [1x1 struct]
%            orig: [1x1 struct]
%        chantype: {14x1 cell}
%        chanunit: {14x1 cell}

assert(ft_filetype(filename, 'eeglab_set'))
assert(hdr.Fs==500);
assert(hdr.nChans==14);
assert(hdr.nSamplesPre==100);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function word_indices = processEmail(email_contents)
%PROCESSEMAIL preprocesses a the body of an email and
%returns a list of word_indices
%   word_indices = PROCESSEMAIL(email_contents) preprocesses
%   the body of an email and returns a list of indices of the
%   words contained in the email.
%

% Load Vocabulary
vocabList = getVocabList();

% Init return value
word_indices = [];

% ========================== Preprocess Email ===========================

% Find the Headers ( \n\n and remove )
% Uncomment the following lines if you are working with raw emails with the
% full headers

% hdrstart = strfind(email_contents, ([char(10) char(10)]));
% email_contents = email_contents(hdrstart(1):end);

% Lower case
email_contents = lower(email_contents);

% Strip all HTML
% Looks for any expression that starts with < and ends with > and replace
% and does not have any < or > in the tag it with a space
email_contents = regexprep(email_contents, '<[^<>]+>', ' ');

% Handle Numbers
% Look for one or more characters between 0-9
email_contents = regexprep(email_contents, '[0-9]+', 'number');

% Handle URLS
% Look for strings starting with http:// or https://
email_contents = regexprep(email_contents, ...
                           '(http|https)://[^\s]*', 'httpaddr');

% Handle Email Addresses
% Look for strings with @ in the middle
email_contents = regexprep(email_contents, '[^\s]+@[^\s]+', 'emailaddr');

% Handle $ sign
email_contents = regexprep(email_contents, '[$]+', 'dollar');

% ========================== Tokenize Email ===========================

% Output the email to screen as well
fprintf('\n==== Processed Email ====\n\n');

% Process file
l = 0;

while ~isempty(email_contents)

    % Tokenize and also get rid of any punctuation
    [str, email_contents] = ...
       strtok(email_contents, ...
              [' @$/#.-:&*+=[]?!(){},''">_<;%' char(10) char(13)]);

    % Remove any non alphanumeric characters
    str = regexprep(str, '[^a-zA-Z0-9]', '');

    % Stem the word
    % (the porterStemmer sometimes has issues, so we use a try catch block)
    try str = porterStemmer(strtrim(str));
    catch str = ''; continue;
    end;

    % Skip the word if it is too short
    if length(str) < 1
       continue;
    end

    % Look up the word in the dictionary and add to word_indices if
    % found
    % ====================== YOUR CODE HERE ======================
    % Instructions: Fill in this function to add the index of str to
    %               word_indices if it is in the vocabulary. At this point
    %               of the code, you have a stemmed word from the email in
    %               the variable str. You should look up str in the
    %               vocabulary list (vocabList). If a match exists, you
    %               should add the index of the word to the word_indices
    %               vector. Concretely, if str = 'action', then you should
    %               look up the vocabulary list to find where in vocabList
    %               'action' appears. For example, if vocabList{18} =
    %               'action', then, you should add 18 to the word_indices
    %               vector (e.g., word_indices = [word_indices ; 18]; ).
    %
    % Note: vocabList{idx} returns a the word with index idx in the
    %       vocabulary list.
    %
    % Note: You can use strcmp(str1, str2) to compare two strings (str1 and
    %       str2). It will return 1 only if the two strings are equivalent.
    %

    for i = 1:size(vocabList)
      if (strcmp(str, vocabList(i)))
        word_indices = [ word_indices; i];
      endif
    endfor

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

    % Print to screen, ensuring that the output lines are not too long
    if (l + length(str) + 1) > 78
        fprintf('\n');
        l = 0;
    end
    fprintf('%s ', str);
    l = l + length(str) + 1;

end

% Print footer
fprintf('\n\n=========================\n');

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
function res = java_call(obj, objMethod, methodProto, varargin)
% JAVA_CALL: Call a method on a Java object on the EDT.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2013

% Parse inputs
if (nargin < 3)
    methodProto = [];
    nParam = 0;
else
    nParam = length(varargin);
end
% Initialize output
res = [];
% Call the creation with the latest available function
if exist('javaMethodEDT', 'builtin')
    if (nParam == 0)
        if (nargout == 1)
            res = javaMethodEDT(objMethod, obj);
        else
            javaMethodEDT(objMethod, obj);
        end
    else
        if (nargout == 1)
            res = javaMethodEDT(objMethod, obj, varargin{:});
        else
            javaMethodEDT(objMethod, obj, varargin{:});
        end
    end
else
    if (nParam == 0)
        if (nargout == 1)
            res = awtinvoke(obj, objMethod);
        else
            awtinvoke(obj, objMethod);
        end
    else
        if (nargout == 1)
            res = awtinvoke(obj, [objMethod, '(', methodProto, ')'], varargin{:});
        else
            awtinvoke(obj, [objMethod, '(', methodProto, ')'], varargin{:});
        end
    end
end

author: brainstorm-tools


text:
function [structout] = maternfit(varargin)
%MATERNFIT  Parametric spectral fit to the Matern form. [with A. Sykulski]
%
%   MATERNFIT performs a parametric fit of the spectrum of a time series to
%   that expected for a Matern process plus optional spin.
%
%   The time series may either be real-valued or complex-valued.
%
%   The Matern process plus spin has a spectrum given by, see MATERNSPEC,
%
%        SPP(F) = SIGMA^2 / [(F-NU)^2/LAMBDA^2 + 1]^ALPHA
%                           / (LAMBDA * MATERNC(ALPHA)) + 2*pi * EPSILON^2
%
%   where SIGMA is the standard deviation, NU is a frequency shift, LAMBDA
%   is a damping coefficient, ALPHA is 1/2 of the spectral slope, and
%   EPSILON^2 is the variance of an optional additive noise component.
%
%   The coefficient LAMBDA * MATERNC(ALPHA) in the above form lets us
%   parameterize the spectrum in terms of the process variance SIGMA^2.
%
%   The optimal parameters are found using a frequency-domain maximum
%   likelihood method, accounting for both aliasing and spectral blurring.
%
%   For further details on the parameter inference method method, see
%
%     Sykulski, Olhede, Lilly, and Danioux (2016).  Lagrangian time series
%        models for ocean surface drifter trajectories. Journal of the
%        Royal Statisical Society, Series C. 65 (1): 29--50.
%
%     Sykulski, Olhede, Guillaumin, Lilly, and Early (2019). The de-biased
%        Whittle likelihood. Biometrika, 106 (2): 251--266.
%
%   For details on the Matern process and its spectrum, see:
%
%     Lilly, Sykulski, Early, and Olhede, (2017).  Fractional Brownian
%        motion, the Matern process, and stochastic modeling of turbulent
%        dispersion.  Nonlinear Processes in Geophysics, 24: 481--514.
%   __________________________________________________________________
%
%   Usage
%
%   FIT=MATERNFIT(DT,Z,FO) where Z is a times series oriented as column
%   vector, returns the result of a fitting the periodogram of Z to a
%   spectrum having a Matern form, fit over all frequencies. The range
%   of frequencies involved in the fit can be modified, as described below.
%
%   DT is the sample interval, which has units of days, while FO is a
%   signed reference frequency in units of radians per day.  FO is used in
%   determining search ranges and initial guesses.  In oceanographic
%   applications, one would normally choose FO as the Coriolis frequency.
%
%   The output argument FIT is a structure which will be described later.
%
%   The default search ranges and initial guess values are as follows:
%
%                     Low  Guess  High
%        SIGMA    =  [0      1    100] * STD(Z)
%        ALPHA    =  [1/2    1     10]
%        LAMBDA   =  [1e-3  2e-3   10] * ABS(FO)
%        EPSILON  =  [0      0      0] * STD(Z)
%        NU       =  [0      0      0] * ABS(FO)
%
%   Note that by default, neither the EPSILON nor NU parameters are used.
%
%   The search ranges and guesses can all be modified, as described below.
%
%   A common special case is the inclusion of an additive noise component.
%   Calling flag MATERNFIT(...'noisy'), which sets
%
%        EPSILON  =  [0     0.01     2] * STD(Z)
%
%   for the range of the EPSILON parameter.  According to the spectral
%   normalizations used here, a white noice time series with standard
%   deviation EPSILON will have a spectral
%   __________________________________________________________________
%
%   Output format
%
%   The following parameters are output as fields of the structure FIT.
%
%      SIGMA     Standard deviation of currents in cm/s
%      ALPHA     Slope parameter
%      LAMBDA    Damping parameter in rad / day
%      NU        Oscillation frequency in rad /day
%      RANGE     Ranges of fit search for each parameter
%      PARAMS    Sub-structure with fields described below
%
%   The associated spectra can be then created from SIGMA, ALPHA, LAMBDA,
%   and NU using MATERNSPEC.
%
%   RANGE is substructure with fields SIGMA, ALPHA, LAMBDA, and NU, such
%   that RANGE.SIGMA specifies the *dimensional* values of the associated
%   range and initial guess with format [MIN GUESS MAX], and so forth for
%   all the other parameters.
%
%   PARAMS is a substructure containing various parameter values
%   characterizing the fit itself:
%
%      PARAMS.DT      Sample rate, as input to MATERNFIT
%      PARAMS.FO      Reference frequency, as input to MATERNFIT
%      PARAMS.A       Index into first frequency F(A) used in the fit
%      PARAMS.B       Index into last frequency F(B) used in the fit
%      PARAMS.P       Number of free parameters used in the fit
%      LIKE           Negative of the log-likelihood
%      AICC           Akaike Information Criterion, corrected version
%      ERR            Normalized error of fit to log spectra
%      EXITFLAG       The exit flag from the optimization routine
%      ITER           The number of iterations in the optimization routine
%      PARAMS.SIDE    String specifying frequency side options
%      PARAMS.ALG     String specifying algorithm options
%      PARAMS.VER     String specifiying raw or difference option
%      PARAMS.CORES   String specifying series or parallel computation
%
%   These parameters are described in more detail below.
%   __________________________________________________________________
%
%   Options for multiple input time series
%
%   FIT=MATERNFIT(DT,Z,FO) may have Z being a matrix with N columns, or a
%   cell array of N different time series.  In both of these cases, DT and
%   FO may be scalars or length N arrays.
%
%   In these cases, the fields SIGMA, ALPHA, LAMBDA, and NU of FIT will
%   also be arrays with N elements, as all non-string fields of PARAMS. The
%   fields of RANGE will then all be N x 3 arrays.
%
%   Alternatively, MATERNFIT(...,'average') with Z a matrix or a cell
%   array of matrices causes the columns of each matrix to be interpreted
%   as members of an ensemble, averaging over columns to create an average
%   spectrum. One fit per matrix is returned, rather than one per column.
%   __________________________________________________________________
%
%   Specifiying frequencies
%
%   MATERNFIT(DT,Z,FO,RA,RB), where RA and RB are both real-valued scalars,
%   applies the fit to only frequencies in the range
%
%               ABS(FO)*RA < F < ABS(FO)*RB
%
%   with the default behavior corresponding to MATERNFIT(DT,Z,FO,0,INF).
%   Thus RA is the smallest permissible ratio of F to the reference
%   frequency, while RB is similarly the largest permissible ratio.
%
%   MATERNFIT(DT,Z,FO,RA,[RB,RN]) also works, where the fifth argument is
%   an array of length two. In this case the fit is applied to the range
%
%               ABS(FO)*RA < F < MIN( ABS(FO)*RB, PI/DT*RN )
%
%   RB is the largest permissible ratio of F to the reference frequency,
%   while RN is the largest permissible ratio of F to the Nyquist PI/DT.
%   The fit extends to the smaller of these two frequencies.
%
%   If RA and RB are imaginary numbers, rather than real numbers, then
%   the fit is only applied to frequencies in the range
%
%               IMAG(RA) < F < IMAG(RB)
%
%   that is, the range is found without scaling by the reference frequency.
%
%   MATERNFIT can create the fit by utilizing both positive and negative
%   frequency sides of the spectrum, the default behavior, or to only one
%   side (plus the zero frequency). This is modified as follows:
%
%      MATERNFIT(...,'both',...), the default, uses both sides.
%      MATERNFIT(...,'positive',...), uses the side where F/FO is positive.
%      MATERNFIT(...,'negative',...), uses the side where F/FO is negative.
%
%   Thus, changing the sign of the reference frequency also changes the
%   side of the spectrum to be fit using the 'postive' or 'negative' flags.
%   __________________________________________________________________
%
%   Parameter range specification
%
%   The default search ranges and guess values can be modified.  As an
%   example, to modify the default range and guess for the SIGMA parameter
%   corresponding to the background flow, use
%
%       MATERNFIT(...,'range.sigma',[MIN GUESS MAX],...)
%
%   and so forth for other parameters. SIGMA ranges input to MATERNFIT
%   represent *fractions* of the total signal standard deviation.  Thus
%
%       MATERNFIT(...,'range.sigma',[0 1 100],...)
%
%   corresponds to the default setting.  The upper limit is set very high
%   because occasionally an optimum spectral fit is found that has much
%   larger total variance than the signal.
%
%   The ranges of LAMBDA and NU are nondimensional values representing
%   fractions of the magnitude of the reference frequency ABS(FO).  Thus
%
%       MATERNFIT(...,'range.lambda',[1/1000 2/1000 1],...)
%
%   corresponds to the default range for LAMBDA.
%
%   The slope parameter ALPHA can take on a value of no less than 1/2, so
%   the lower range for ALPHA cannot be less than 1/2.
%
%   This approach can be used to omit parameters from the fit, by setting
%   the MIN, GUESS, and MAX values to be identical.  For example,
%
%      MATERNFIT(...,'range.alpha',[2 2 2])
%
%   sets the ALPHA parameter to a value of 2.  Fixed parameters are then
%   not included in the optimization, thus speeding up the fit.
%   __________________________________________________________________
%
%   Parameter value specification
%
%   Parameters can also be set to particular values for each time series.
%   This is done by setting the 'value' field as follows
%
%       MATERNFIT(...,'value.sigma',SIGMA,...)
%
%   and so forth for the other parameters.
%
%   Here SIGMA is an array of the same length as the number of time series
%   in Z.  Thus SIGMA is a scalar if Z is a single array, an array of
%   LENGTH(Z) is Z is a cell array, or of SIZE(Z,2) if Z is a matrix.
%
%   Note that the values specified in this way are actual dimensional
%   values, not nondimensional values as with setting the range.
%
%   This approach works by internally setting the dimension MIN, GUESS, and
%   MAX to the value specified for each time series, overriding the default
%   choices. This will be reflected in the output RANGE fields.
%
%   If both ranges and values are specified for the same parameter, the
%   value settings take precedence.
%   __________________________________________________________________
%
%   Numerical options
%
%   MATERNFIT has options for specifying numerical details of the fit and
%   the optimization algorithm.
%
%   MATERNFIT can use one of two different optimization algorithms.
%
%     MATERNFIT(...,'bnd',...), the default, uses FMINSEARCHBND by J.
%       D'Errico included with JLAB in accordance with its license terms.
%       This in turn calls Matlab's FMINSEARCH using Nelder-Mead.
%
%     MATERNFIT(...,'con',...) alternately uses FMINCON with the default
%       interior-point algorithm.  This requires Matlab's Optimization
%       Toolbox to be installed.  This is mainly for testing.
%
%     MATERNFIT(...,'nlo',...) uses the Nelder-Mead algorithm from the
%       NLopt toolbox at https://nlopt.readthedocs.io.  This requires NLopt
%       to be installed.  Again, this is mainly for testing at the moment.
%
%   In tests, FMINCON is generally faster, and most of the fits agree
%   closely with those using FMINSEARCHBND.  However, occasionally FMINCON
%   produces fits that are significantly worse than those obtained from
%   FMINSEARCHBND, which is why the latter is preferred by default.
%
%   MATERNFIT can employ two different versions of the fit.
%
%     MATERNFIT(...,'difference',...) estimates the Matern parameters by
%       fitting the first difference of the time series to the first
%       difference of a Matern process.  This amounts to a form of pre-
%       whitening, and is the default behavior when a taper is not input.
%     MATERNFIT(...,'raw',...) fits the time series directly to a Matern.
%       This is the default behavior when a taper *is* input.
%
%   These choices are reflected in the output fields PARAMS.ALG and
%   PARAMS.VER, respectively.
%   __________________________________________________________________
%
%   Tapering
%
%   The default behavior of fitting the first difference of the spectrum
%   is usually sufficient to account for spectral blurring.  For very steep
%   spectra or those with a very large dynamic range, this is no longer the
%   case, because leakage from high-energy portions of the spectra will
%   obscure the structure of low-energy portions.
%
%   To addess this, MATERNFIT can optionally perform the fit by fitting
%   to the tapered and aliased spectrum, correctly accounting for the
%   influence of tapering.  This is accomplished with
%
%      MATERNFIT(...,'taper',PSI,...)
%
%   where PSI is a data taper of the same length as Z.  If Z is a cell
%   array, then PSI is a cell array of data tapers having the same length
%   as the components of Z.  PSI is typically computed by SLEPTAP.
%
%   When a taper is input, the default behavior is *not* to difference
%   the time series, corresponding to the 'raw' option.  If both a taper
%   and the 'difference' flag are both input, then the taper length must
%   be one less than the data length.
%   __________________________________________________________________
%
%   Error
%
%   While the maximum likelikehood method is not about finding the best
%   fit in a least squares sense, a measure of the mean squared error
%   provides a useful measure of evaluating the degree of misfit.
%
%   The error ERR returned by MATERNFIT is the squared difference between
%   the *natural log* of the periodogram and that of the fit spectrum,
%   summed over all frequencies used in the fit, and divided by the sum of
%   the squared natural log of the periodogram over the same frequencies.
%
%   The error is computed over the inertial side, anti-inertial side, or
%   both sides, depending on which frequency range is used for the fit.
%
%   When the 'difference' behavior is employed, as is the default, ERR will
%   reflect the error between the periodogram of the first difference of
%   the time series and the spectrum of the first difference of the fit.
%   __________________________________________________________________
%
%   Parallelization
%
%   With Matlab's Parallel Computing toolbox installed, when Z is a cell
%   array or a matrix, MATERNFIT(...,'parallel') will loop over the
%   elements of Z using a parfor loop to speed things up.
%
%   This choice is reflected in the output field PARAMS.CORES, which
%   will take on the values 'series' (the default) or 'parallel'.
%   __________________________________________________________________
%
%   'maternfit --t' runs a test.
%
%   Usage: fit=maternfit(dt,z,fo);
%          fit=maternfit(dt,z,fo,a,b);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2014--2021 J.M. Lilly and A.M. Sykulski
%                                --- type 'help jlab_license' for details

%   __________________________________________________________________
%
%   Extensions
%
%   Two extensions of the basic Matern form are supported, discussed in
%   more detail in MATERNSPEC.
%
%   An oscillatory version, with nonzero MU, and a noisy version, with
%   nonzero EPSILON, are available as options with these two extensions.
%
%   Generalized Matern
%
%   MATERFIT(...,'general') fits to the spectrum of generalized Matern
%   process. A parameter GAMMA is used with the default range
%
%        GAMMA  =  [1/2     1    20]
%
%   such that the initial guess is the standard Matern process.
%
%   The difference between the generalized Matern and standard Matern is
%   limited to frequency band in the vicinity of the falloff frequency, so
%   unless you are averaging over many time series (decribed below), you
%   can expect a lot of variability in the fit value of GAMMA.
%
%   Extended Matern
%
%   MATERFIT(...,'extended') fits to the spectrum of an extended Matern
%   process. A parameter MU is used with the default range
%
%            MU  =  [1/100  10  100] * ABS(FO)
%
%   and in this case, the ALPHA parameter has the default range
%
%       ALPHA    =  [-1/2    1    10]
%
%   which has a different lower bound than for the standard Matern process.
%   __________________________________________________________________

%   Composite Matern
%
%   MATERFIT(...,'composite') fits to the spectrum of an composite Matern
%   process. A parameter MU is used with the default range

%Notes to self:
%Note, this is all set up to run the 5-parameter Matern process, I have
%just not checked the code.  Also make sure to comment, which I can grab
%from SPECFIT, and also make sure to assign bgtype in output params.
%
%This is basically just a stripped-down version of SPECFIT.  Use Mac's
%Filemerge tool to compare them.

if nargin>0
%    if strcmpi(varargin{1}, '--f')
 %       maternfit_fig,return
    if strcmpi(varargin{1}, '--t')
        maternfit_test,return
    end
end
%--------------------------------------------------------------------------
%Set search ranges for parameters
range=struct;

%Ranges for the background process
range.sigma=[0 1 100];            %Range of standard deviation as fraction of total
range.alpha=[1/2 1 10];           %Range of slope parameter
range.lambda=[1/1000 2/1000 10];  %Range of decay parameter as multiple of Coriolis
range.mu=[0 0 0];                 %Range of damping parameter as multiple of inverse Coriolis
range.nu=[0 0 0];                 %Frequency shift is fixed at NU=0
range.epsilon=[0 0.01 2];         %Noise level
%--------------------------------------------------------------------------
%Initialize value structure
names=fieldnames(range);
for i=1:length(names)
    eval(['value.' names{i} '=[];'])
end
%--------------------------------------------------------------------------
dt=double(varargin{1});
z=varargin{2};
varargin=varargin(3:end);
%--------------------------------------------------------------------------
%Sorting out input options

opts.side='both';            %Determines side: opposite, same, or both
opts.alg='bnd';              %Algorithm: bnd, con, or nlopt
opts.ver='difference';       %Fit version: difference or raw
opts.cores='serial';         %Serial or parallel computation
opts.bgtype='mat';      %Type of model for background
opts.noise='clean';          %Noisy or clean
opts.cols='fit';             %Fit or average

psi=[];                      %The default data taper is the empty taper
verstrin=false;              %Flag for whether or not version string is input
for i=1:30
    if ischar(varargin{end})
        if strcmpi(varargin{end}(1:3),'con')||strcmpi(varargin{end}(1:3),'bnd')||strcmpi(varargin{end}(1:3),'nlo')
            opts.alg=varargin{end};
            if strcmpi(opts.alg(1:3),'con')
                if exist('fmincon.m')~=2
                    error('Sorry, MATERNFIT with ALG=''con'' requires the Optimization Toolbox')
                end
            end
        elseif strcmpi(varargin{end}(1:3),'raw')||strcmpi(varargin{end}(1:3),'dif')||strcmpi(varargin{end}(1:3),'sec')
            opts.ver=varargin{end};
            verstrin=true;
        elseif strcmpi(varargin{end}(1:3),'bot')||strcmpi(varargin{end}(1:3),'pos')||strcmpi(varargin{end}(1:3),'neg')
            opts.side=varargin{end};
        elseif strcmpi(varargin{end}(1:3),'ser')||strcmpi(varargin{end}(1:3),'par')
            opts.cores=varargin{end};
        elseif strcmpi(varargin{end}(1:3),'sta')||strcmpi(varargin{end}(1:3),'com')
            opts.model=varargin{end};
        elseif strcmpi(varargin{end}(1:3),'mat')||strcmpi(varargin{end}(1:3),'exp')||strcmpi(varargin{end}(1:3),'ext')||strcmpi(varargin{end}(1:3),'gen')
            opts.bgtype=varargin{end};
        elseif strcmpi(varargin{end}(1:3),'noi')||strcmpi(varargin{end}(1:3),'cle')
            opts.noise=varargin{end};
        elseif strcmpi(varargin{end}(1:3),'fit')||strcmpi(varargin{end}(1:3),'ave')
            opts.cols=varargin{end};
        end
        varargin=varargin(1:end-1);
    elseif length(varargin)>1
        if ischar(varargin{end-1})&&~ischar(varargin{end})
            if strcmpi(varargin{end-1}(1:3),'tap')
                psi=varargin{end};
            elseif strcmpi(varargin{end-1}(1:3),'ran')
                name=varargin{end-1}(7:end);
                range=setfield(range,name,varargin{end});
            elseif strcmpi(varargin{end-1}(1:3),'val')
                name=varargin{end-1}(7:end);
                value=setfield(value,name,varargin{end});
            end
            varargin=varargin(1:end-2);
        end
    end
end
%--------------------------------------------------------------------------
%Make sure ranges and flags are set correctly for specified model type
if strcmpi(opts.bgtype(1:3),'exp')
    range.alpha=[-1/2 -1/2 -1/2];
    range.mu=[1/100 10 100];
    %range.lambda=[1/100000 2/1000 10];  %Range of decay parameter as multiple of Coriolis
elseif strcmpi(opts.bgtype(1:3),'ext')
    range.alpha=[-1/2 1 10];
    range.mu=[1/100 10 100];
elseif strcmpi(opts.bgtype(1:3),'gen')
    range.mu=[1/2 1 20];
elseif strcmp(opts.bgtype(1:3),'mat')
    if range.alpha(1)<=1/2
        %Correction for alpha being less than permitted value
        range.alpha(1)=1/2+1e-10;
    end
end
if strcmpi(opts.noise(1:3),'cle')
    range.epsilon=[0 0 0];
end
if isfield(range,'gamma')
    range.mu=range.gamma;
    range=rmfield(range,'gamma');
end
%--------------------------------------------------------------------------
%Convert value structure to a numeric array
if iscell(z)
    N=length(z);
else
    N=size(z,2);
end
values=nan*zeros(6,3);
values(1,:)=range.sigma;
values(2,:)=range.alpha;
values(3,:)=range.lambda;
values(4,:)=range.mu;
values(5,:)=range.nu;
values(6,:)=range.epsilon;
if ~isempty(value.sigma),   values(1,:)=value.sigma;end
if ~isempty(value.alpha),   values(1,:)=value.alpha;end
if ~isempty(value.lambda),  values(1,:)=value.lambda;end
if ~isempty(value.mu),      values(1,:)=value.mu;end
if ~isempty(value.nu),      values(1,:)=value.nu;end
if ~isempty(value.epsilon), values(1,:)=value.epsilon;end
%Have to convert this to numeric values because I can't pass a structure
%through to the workers, for some unknown reason
%--------------------------------------------------------------------------
%Frequency range
flow=0;
fhigh=inf;
fc=varargin{1};
varargin=varargin(2:end);

if length(varargin)==2
    flow=varargin{1};
    fhigh=varargin{2};
end
if length(fhigh)==1
    fhigh=[fhigh inf];
end
%--------------------------------------------------------------------------
%Re-sizing of dt and fc
dt=dt(:);

if length(dt)==1
    if iscell(z)
        dt=dt+zeros(length(z),1);
     elseif ~strcmpi(opts.cols(1:3),'ave')
        dt=dt+zeros(size(z,2),1);
    end
end

fc=fc(:);
if length(fc)==1
    if iscell(z)
        fc=fc+zeros(length(z),1);
    elseif ~strcmpi(opts.cols(1:3),'ave')
        fc=fc+zeros(size(z,2),1);
    end
end
%--------------------------------------------------------------------------
%Difference z if requested, and compute standard deviation

%Default to raw algorithm if taper is input
if ~isempty(psi)&&~verstrin
    opts.ver='raw';
end

if ~iscell(z)
    stdz=std(z);
    if strcmpi(opts.ver(1:3),'dif')
        z=diff(z);
    elseif strcmpi(opts.ver(1:3),'sec')
        z=diff(diff(z));
    end
else
    for i=1:length(z)
        stdz{i}=std(z{i});
        if strcmpi(opts.ver(1:3),'dif')
            z{i}=diff(z{i});
        elseif strcmpi(opts.ver(1:3),'sec')
            z{i}=diff(diff(z{i}));
        end
    end
end
%--------------------------------------------------------------------------
%Compute data window from taper or periodogram
if ~isempty(psi)
    if iscell(psi)
        for i=1:length(psi)
            if size(psi{i},1)~=size(z{i},1)
                error('Length of taper does not match length of data.')
            end
            win{i}=conv(psi{i}(:),psi{i}(:));
            win{i}=win{i}(end-length(z{i})+1:end);
        end
    else
        if size(psi,1)~=size(z,1)
            error('Length of taper does not match length of data.')
        end
        win=conv(psi(:),psi(:));
        win=win(end-size(z,1)+1:end);
    end
else
    if iscell(z)
        for i=1:length(z)
            N=length(z{i});
            win{i}=[N:-1:1]'./N;
            psi{i}=[];
        end
    else
        N=size(z,1);
        win=[N:-1:1]'./N;
        psi=[];
    end
end
%--------------------------------------------------------------------------
if ~iscell(z)&&( size(z,2)==1 || strcmpi(opts.cols(1:3),'ave'))
    %Single time series input, no loop
    [x,xa,xb,xo,like,aicc,P,a,b,err,exitflag,iters]=...
        maternfit_one(dt,z,flow,fhigh,opts,fc,psi,win,stdz,values);
    %The rest of this simply implements different types of loops
elseif ~iscell(z)
    %Loop over matrix columns
    N=size(z,2);[x,xa,xb,xo]=vzeros(N,6);
    [like,aicc,P,a,b,err,exitflag,iters]=vzeros(N,1);
    if strcmpi(opts.cores(1:3),'ser')
        %For series loop over matrix columns
        for i=1:size(z,2)
            disp(['MATERNFIT processing series ' int2str(i) ' of ' int2str(size(z,2)) '.'])
            [x(i,:),xa(i,:),xb(i,:),xo(i,:),like(i),aicc(i),P(i),a(i),b(i),err(i),exitflag(i),iters(i)]=...
                maternfit_one(dt(i),z(:,i),flow,fhigh,opts,fc(i),psi,win,stdz(i),values);
        end
    elseif strcmpi(opts.cores(1:3),'par')
        %For parallel loop over matrix columns
        disp('MATERNFIT using parallel for loop...')
        parfor i=1:size(z,2)
            %dt(i),flow,fhigh,opts,fc(i),fit
            disp(['MATERNFIT processing series ' int2str(i) ' of ' int2str(size(z,2)) '.'])
            [x(i,:),xa(i,:),xb(i,:),xo(i,:),like(i),aicc(i),P(i),a(i),b(i),err(i),exitflag(i),iters(i)]=...
                maternfit_one(dt(i),z(:,i),flow,fhigh,opts,fc(i),psi,win,stdz(i),values);
        end
    end
elseif iscell(z)
    %Loop over cell array
    N=length(z);[x,xa,xb,xo]=vzeros(N,6);
    [like,aicc,P,a,b,err,exitflag,iters]=vzeros(N,1);
    if strcmpi(opts.cores(1:3),'ser')
        %For series loop over cell array
        for i=1:length(z)
            disp(['MATERNFIT processing series ' int2str(i) ' of ' int2str(length(z)) '.'])
            [x(i,:),xa(i,:),xb(i,:),xo(i,:),like(i),aicc(i),P(i),a(i),b(i),err(i),exitflag(i),iters(i)]=...
                maternfit_one(dt(i),z{i},flow,fhigh,opts,fc(i),psi{i},win{i},stdz{i},values);
        end
    elseif strcmpi(opts.cores(1:3),'par')
        %For parallel loop over cell array
        disp('MATERNFIT using parallel for loop...')
        parfor i=1:length(z)
            disp(['MATERNFIT processing series ' int2str(i) ' of ' int2str(length(z)) '.'])
            [x(i,:),xa(i,:),xb(i,:),xo(i,:),like(i),aicc(i),P(i),a(i),b(i),err(i),exitflag(i),iters(i)]=...
                maternfit_one(dt(i),z{i},flow,fhigh,opts,fc(i),psi{i},win{i},stdz{i},values);
        end
    end
end

sigma=x(:,1);
alpha=x(:,2);
lambda=x(:,3);
mu=x(:,4);
nu=x(:,5);
epsilon=x(:,6);

clear range
range.sigma=  [xa(:,1)  xo(:,1)  xb(:,1)];
range.alpha=  [xa(:,2)  xo(:,2)  xb(:,2)];
range.lambda= [xa(:,3)  xo(:,3)  xb(:,3)];
if strcmpi(opts.bgtype(1:3),'gen')
    range.gamma=[xa(:,4)  xo(:,4)  xb(:,4)];
else
    range.mu=[xa(:,4)  xo(:,4)  xb(:,4)];
end
range.nu=     [xa(:,5)  xo(:,5)  xb(:,5)];
range.epsilon=[xa(:,6)  xo(:,6)  xb(:,6)];

structout=[];

use opts

make params dt fc a b P like aicc err exitflag iters side alg ver cores
%make params dt fc a b P like aicc err exitflag iters side alg ver bgtype cores

if strcmpi(opts.bgtype(1:3),'gen')
    gamma=mu;
    make structout sigma alpha lambda gamma nu epsilon range params
else
    make structout sigma alpha lambda mu nu epsilon range params
end

%remove some fields I don't use
if strcmpi(opts.bgtype(1:3),'mat')
    structout=rmfield(structout,'mu');
    structout.range=rmfield(structout.range,'mu');
end
if strcmpi(opts.noise(1:3),'cle')
    structout=rmfield(structout,'epsilon');
    structout.range=rmfield(structout.range,'epsilon');
end
if all(range.nu==0)
    structout=rmfield(structout,'nu');
    structout.range=rmfield(structout.range,'nu');
end

%--------------------------------------------------------------------------
function[x,exitflag,iters]=maternfit_optimize(dt,N,om,spp,snn,index,xa,xo,xb,opts,realflag)
%Optimization for all parameters using both sides
%Optimizations like to work with doubles... see comment below
xa=double(xa);
xo=double(xo);
xb=double(xb);
%xa,xo,xb

xanorm=(xa./xo);
xbnorm=(xb./xo);
xanorm(~isfinite(xanorm))=0;
xbnorm(~isfinite(xbnorm))=0;

guess=ones(size(xo));
guess((xanorm==0)&(xbnorm==0))=0;

tol=1e-6;%tol=1e-3;
if strcmpi(opts.alg(1:3),'bnd')
    options=optimset('MaxFunEvals',10000,'MaxIter',10000,'TolFun',tol,'TolX',tol);
    [xf,fval1,exitflag,struct]=fminsearchbnd(@(z) specmodel(dt,z.*xo,spp,snn,index,N,opts,realflag),...
        guess,xanorm,xbnorm,options);
    iters=struct.iterations;
elseif strcmpi(opts.alg(1:3),'con')
    options=optimoptions('fmincon','Algorithm','interior-point','MaxFunEvals',10000,'MaxIter',10000,'TolFun',tol,'TolX',tol,'Display','off');
    [xf,fval1,exitflag,struct]=fmincon(@(z) specmodel(dt,z.*xo,spp,snn,index,N,opts,realflag),...
        guess,[],[],[],[],xanorm,xbnorm,[],options);
    iters=struct.iterations;
elseif strcmpi(opts.alg(1:3),'nlo')
    %opt.algorithm = NLOPT_GN_DIRECT_L;
    opt.algorithm = NLOPT_LN_NELDERMEAD;  %Works, not too much slower than Matlab
    %opt.algorithm = NLOPT_LN_PRAXIS;   %Works,but slower than NM
    %opt.algorithm = NLOPT_GN_CRS2_LM;
    %opt.algorithm = NLOPT_GD_STOGO;
    %opt.algorithm = NLOPT_GN_ISRES;
    %opt.algorithm = NLOPT_GN_ESCH;
    %opt.algorithm = NLOPT_LN_COBYLA;
    %opt.algorithm =NLOPT_LN_SBPLX;  %Works,but slower than NM
    %opt.algorithm = NLOPT_LN_BOBYQA;  %Somewhat slower than NM, and different results
    %opt.algorithm = NLOPT_AUGLAG;    opt.lower_bounds = xanorm;
    opt.lower_bounds = xanorm;
    opt.upper_bounds = xbnorm;
    opt.min_objective = @(z) specmodel(dt,z.*xo,spp,snn,index,N,opts,realflag);
    opt.fc_tol = [tol tol tol tol];
    opt.xtol_rel = tol;
    [xf, fmin, exitflag] = nlopt_optimize(opt, guess);
    iters=nan;  %Not sure how to return this
end

x=xf.*xo;% scale back to correct units
[like,Spp,Snn,f]=specmodel(dt,x,spp,snn,index,N,opts,realflag);
%--------------------------------------------------------------------------
function [like,Spp,Snn,f,indexn]=specmodel(dt,x,spphat,snnhat,index,N,opts,realflag)
% Computes the value of the Whittle likelihood for parameter set X,
% given periodogram estimates SPPHAT and SNNHAT of positive and negative
% sides of the spectrum, respectively.  Only frequencies in the locations
% given by INDEX are used to compute the likelihood.  N is the data length.

%Optimizations like to work with doubles... this prevents a strange error
%where, inside the optimization loop, the spectrum takes on negative values
x=double(x);

if strcmpi(opts.ver(1:3),'dif')
    N=N+1;
elseif strcmpi(opts.ver(1:3),'sec')
    N=N+2;
end

%Compute the blurred spectrum by first computing the covariance.
R=maternfit_materncov(dt,N,x,realflag,opts.bgtype);
[f,Spp,Snn]=maternfit_blurspec(dt,R,opts.ver,opts.win);

%length(R)
%[f,Spp2,Snn2]=blurspec(1,R,opts.ver,'window',opts.win);aresame(Spp,Spp2)
%Same answers but slower

%length(index),N,length(f),length(Spp)
indexn=fixindex(N,index);        %This deals with zero and Nyquist for negative frequencies
if strcmpi(opts.side(1:3),'pos')
    %Include only positive frequencies
    like=sum(log(Spp(index))+spphat(index)./Spp(index));
elseif strcmpi(opts.side(1:3),'neg')
    %Include only negative frequencies
    like=sum(log(Snn(index))+snnhat(index)./Snn(index));
else
    % vsize(Spp,spphat,index)
    likepp=sum(log(Spp(index))+spphat(index)./Spp(index));
    likenn=sum(log(Snn(indexn))+snnhat(indexn)./Snn(indexn));
    like=likepp+likenn;
end

%FMINCON likes to work with doubles
like=double(like);

%hold on,plot(f,spphat),plot(-f, Spp)
%--------------------------------------------------------------------------
%function[x,xar,xbr,like,aicc,P,f,spp,snn,Spp,Snn,a,b,err]=maternfit_one(dt,z,flow,fhigh,side,alg,ver,range,fc)
function[x,xa,xb,xo,like,aicc,P,a,b,err,exitflag,iters]=maternfit_one(dt,z,flow,fhigh,opts,fc,psi,win,stdz,values)

%This is pretty dumb, but Matlab doesn't want me to pass a structure
%through parfor.  Thus I have to convert to a matrix and then back again.
sigma=values(1,:);
alpha=values(2,:);
lambda=values(3,:);
mu=values(4,:);
nu=values(5,:);
epsilon=values(6,:);

[x,xa,xb,xo]=vzeros(1,6,nan);
[like,aicc,P,a,b,err,exitflag,iters]=vzeros(1,1,nan);
if anyany(~isfinite(z))
    return
end

N=size(z,1);
realflag=isreal(z);

%Put window into options to simplify argument passing
opts.win=win;

if strcmpi(opts.cols(1:3),'ave')
    stdz=sqrt(vmean(squared(stdz),2));
end

%Multiply standard deviation ranges by sample standard deviation
sigma    = sigma*stdz;
epsilon  = epsilon*stdz;

%Multiply damping and frequency ranges by Coriolis frequency
lambda  = lambda*abs(fc);
if ~strcmpi(opts.bgtype(1:3),'gen')
    mu      = mu.*abs(fc);   %Is this a time or frequency scale?
end
xa=[sigma(1)  alpha(1)   lambda(1)  mu(1)  nu(1)  epsilon(1)];
xo=[sigma(2)  alpha(2)   lambda(2)  mu(2)  nu(2)  epsilon(2)];
xb=[sigma(3)  alpha(3)   lambda(3)  mu(3)  nu(3)  epsilon(3)];

%Periodogram from mspec  ... use dt = 1
%vsize(dt,z,psi),length(find(isfinite(z))),isreal(z)
[om,spp,snn]=mspec(dt,z,psi);
if isreal(z)
    snn=spp;
end

if strcmpi(opts.cols(1:3),'ave')
    spp=vmean(spp,2);
    snn=vmean(snn,2);
end

if ~isnan(fc)
    if isreal(flow(1))
        a=find(om>=flow*abs(fc),1,'first');   %Look up multiple of Corilois frequency
    else
        a=find(om>=imag(flow),1,'first');
    end
    if isreal(fhigh(1))
        b=find(om<min(fhigh(1)*abs(fc),fhigh(2)*pi/dt),1,'last');    %Look up multiple of Corilios frequency
    else
        b=find(om<min(imag(fhigh(1)),real(fhigh(2))*pi/dt),1,'last');
    end
else
    a=1;
    b=length(om);
end

%figure,plot(om,spp),hold on,plot(-om,snn),ylog,vlines([om(a) om(b) -om(a) -om(b)])

%[like,Spp,Snn,f,indexn]=specmodel(dt,x(:,4),spp,snn,a:b,N,opts);
%figure,subplot(1,2,1),plot(f,[Spp spp]),ylog,vlines(x(9,2),'r')
%subplot(1,2,2),plot(f,[Snn snn]),ylog,

%Fit to background on requested side
%xa,xo,xb
%vsize(dt,N,om,spp,snn,a:b,xa,xo,xb,opts,realflag)
%[xa;xo;xb]
[x,exitflag,iters]=maternfit_optimize(dt,N,om,spp,snn,a:b,xa,xo,xb,opts,realflag);

%Final value of likelihood and spectra
%vsize(dt,x,spp,snn,a:b,N,opts,realflag)
[like,Spp,Snn,f,indexn]=specmodel(dt,x,spp,snn,a:b,N,opts,realflag);

% Compute the number of free parameters
P=0;
for i=1:size(x,1)
    if ~isnan(xo(i,1))
        P=P+length(find(xa(i,:)~=xb(i,:)));
    end
end

% Compute the value of the AICC information criterion
if (strcmpi(opts.side(1:3),'opp')||strcmpi(opts.side(1:3),'sam'))
    M=length(a:b);
else
    M=length(a:b)+length(indexn);
end
aicc=2*like + frac(4*M*P,2*M-P-1);

%figure,plot(f,[Spp spp]),hold on,plot(-f,[Snn snn]),ylog
%vlines(f(a)),vlines(f(b))

%Compute an error measure
if strcmpi(opts.side(1:3),'sam')
    numer=sum(squared(log(spp(a:b))-log(Spp(a:b))));
    denom=sum(squared(log(spp(a:b))));
    err=frac(numer,denom);
elseif strcmpi(opts.side(1:3),'opp')
    numer=sum(squared(log(snn(a:b))-log(Snn(a:b))));
    denom=sum(squared(log(snn(a:b))));
    err=frac(numer,denom);
else
    numerp=sum(squared(log(spp(a:b))-log(Spp(a:b))));
    denomp=sum(squared(log(spp(a:b))));
    numern=sum(squared(log(snn(a:b))-log(Snn(a:b))));
    denomn=sum(squared(log(snn(a:b))));
    err=frac(numerp+numern,denomp+denomn);
end

vtranspose(x,xa,xb,xo);
vcolon(x,xa,xb,xo);
vtranspose(x,xa,xb,xo);
%--------------------------------------------------------------------------
function[index]=fixindex(N,index)
%This modifies the frequency index appropriate for negative frequencies,
%due to the fact that zero is repeated for both even and odd, while the
%Nyquist is also repeated for even. See comments at MSPEC.

if iseven(N)
    if index(end)==((N/2)+1)
        index=index(1:end-1);
    end
end
if index(1)==1
    index=index(2:end);
end
%--------------------------------------------------------------------------
function[omega,Spp,Snn]=maternfit_blurspec(dt,R,ver,win)
%This is a version of BLURSPEC, stripped down for speed.

%I don't take dt into account during difference
if strcmpi(ver(1:3),'dif')
    R=2*R(1:end-1,:)-R(2:end,:)-[conj(R(2,:));R(1:end-2,:)];
elseif strcmpi(ver(1:3),'sec')
    R=2*R(1:end-1,:)-R(2:end,:)-[conj(R(2,:));R(1:end-2,:)];
    R=2*R(1:end-1,:)-R(2:end,:)-[conj(R(2,:));R(1:end-2,:)];
end

R=R.*win;
R(1,:)=R(1,:)./2;    %Don't forget to divide first element by two
S=dt*2*real(fft(R)); %But I do take it into account in spectrum

S=abs(S);  %Sometimes there are small negative parts after blurring
N=size(R,1);
omega=frac(1,dt)*2*pi*(0:floor(N/2))'./N;
%omega=fourier(N);

Spp=S(1:length(omega),:);
Snn=[S(1,:);S(end:-1:end-length(omega)+2,:)];
%--------------------------------------------------------------------------
function[R]=maternfit_materncov(dt,N,x,realflag,model)
sigma=x(1);
alpha=x(2);
lambda=x(3);
mu=x(4);
nu=x(5);
epsilon=x(6);

tau=dt*[0:N-1]';

if isnan(sigma)
    R=[];
else
    %This is just copied from MATERNCOV, but it is faster to have it internal
    if strcmpi(model(1:3),'mat')
        fact=2*frac(1,gamma(alpha-1/2).*pow2(alpha-1/2));
        R=fact.*((lambda*tau).^(alpha-1/2)).*besselk(abs(alpha-1/2),lambda*tau);
        R(1)=1;%because of being undefiend there
    elseif strcmpi(model(1:3),'ext')
        if alpha==-1/2
            tnorm=sqrt(tau.^2+mu.^2);
            fact=besselk(1,mu.*lambda);
            R=frac(1,fact).*frac(mu,tnorm).*besselk(1,lambda.*tnorm);
        else
            tnorm=lambda.*sqrt(tau.^2+mu.^2);
            fact=1./(abs(mu.*lambda)).^(alpha-1/2)./besselk(abs(alpha-1/2),abs(mu.*lambda));
            R=fact*tnorm.^(alpha-1/2).*besselk(abs(alpha-1/2),tnorm);
        end
    elseif strcmpi(model(1:3),'gen')
        M=10;P=10;  %Specifying oversampling rates for numerical computations
        [f,Spp,Snn]=maternspec(dt,M*N*P,sigma,alpha,lambda/P,mu,'generalized');
        S=[flipud(Snn(2:end));Spp];%plot(Spp),hold on
        Ri=ifft(ifftshift(S))./dt;  %Make sure it's ifftshift not fftshift
        Ri=Ri(1:P:end);
        R=Ri(1:N);%plot(Ri),hold on
        R=R./R(1);
    elseif strcmpi(model(1:3),'com')
        M=10;P=10;  %Specifying oversampling rates for numerical computations
        [f,Spp,Snn]=maternspec(dt,M*N*P,sigma,alpha,lambda/P,mu/P,nu/P,'composite');
        S=[flipud(Snn(2:end));Spp];%figure,plot(S)
        Ri=ifft(ifftshift(S))./dt;  %Make sure it's ifftshift not fftshift
        Ri=Ri(1:P:end);
        R=Ri(1:N);
        R=R./R(1);
    end
    if (nu~=0)&&~strcmpi(model(1:3),'com')
        R=R.*exp(sqrt(-1)*tau*nu);
    end
    R=R.*sigma.^2;
    R(1)=R(1)+squared(epsilon);
end

if realflag
     R=real(R);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%End of function body; begin tests and figures
function[]=maternfit_test
%--------------------------------------------------------------------------
sigo=17;
alpha=1.5;
h=1/10;

rng(0);
dt=1;
z=maternoise(dt,1000,sigo,alpha,h);
fit=maternfit(dt,z,frac(1,2)*pi);
bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.06)&&aresame(fit.lambda,h,0.002);
reporttest('MATERNFIT recovers Matern parameters with unit sample rate',allall(bool))

fit2=maternfit(dt,z,frac(1,2)*pi,'range.sigma',[0 17 100],'range.alpha',[0 1.5 100],'range.lambda',[1/1000 1/10 10]);

clear bool
bool(1)=aresame(fit.sigma,fit2.sigma,1e-2);
bool(2)=aresame(fit.alpha,fit2.alpha,1e-4);
bool(3)=aresame(fit.lambda,fit2.lambda,1e-4);

reporttest('MATERNFIT is independent of initial guess',allall(bool))

rng(0);
epsilono=2;
zn=epsilono.*(randn(length(z),1)+1i.*randn(length(z),1))./sqrt(2);
fit=maternfit(dt,zn,frac(1,2)*pi,'noisy','value.sigma',0);
bool=(abs(fit.sigma)<1e-7)&&aresame(abs(fit.epsilon./epsilono-1),0,0.05);
reporttest('MATERNFIT recovers Matern parameters with unit sample rate, noise only',allall(bool))

rng(0);
fit=maternfit(dt,z+zn,frac(1,2)*pi,'noisy');
bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.5)&&aresame(fit.lambda,h,0.1)&&aresame(abs(fit.epsilon./epsilono-1),0,0.3);
reporttest('MATERNFIT recovers Matern parameters with unit sample rate, noisy version',allall(bool))
%--------------------------------------------------------------------------
sigo=17;
alpha=1.5;
h=1/10;

rng(0);
dt=1;
z=maternoise(dt,1000,sigo,alpha,h,'real');
fit=maternfit(dt,z,frac(1,2)*pi);
bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.06)&&aresame(fit.lambda,h,0.0025);
reporttest('MATERNFIT recovers Matern parameters with unit sample rate and real signal',allall(bool))

fit2=maternfit(dt,z,frac(1,2)*pi,'range.sigma',[0 17 100],'range.alpha',[0 1.5 100],'range.lambda',[1/1000 1/10 10]);

clear bool
bool(1)=aresame(fit.sigma,fit2.sigma,1e-2);
bool(2)=aresame(fit.alpha,fit2.alpha,1e-4);
bool(3)=aresame(fit.lambda,fit2.lambda,1e-4);

reporttest('MATERNFIT is independent of initial guess for real signal',allall(bool))

rng(0);
epsilono=2;
zn=epsilono.*randn(length(z),1);
fit=maternfit(dt,zn,frac(1,2)*pi,'noisy','value.sigma',0);
bool=(abs(fit.sigma)<1e-7)&&aresame(abs(fit.epsilon./epsilono-1),0,0.1);
reporttest('MATERNFIT recovers Matern parameters with unit sample rate and real signal, noise only',allall(bool))

rng(0);
fit=maternfit(dt,z+zn,frac(1,2)*pi,'noisy');
bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.5)&&aresame(fit.lambda,h,0.1)&&aresame(abs(fit.epsilon./epsilono-1),0,0.3);
reporttest('MATERNFIT recovers Matern parameters with unit sample rate and real signal, noisy version',allall(bool))
%--------------------------------------------------------------------------

% rng(0);
% dt=3600;
% z=maternoise(dt,1000,sigo,alpha,h./dt);
% fit=maternfit(dt,real(z),frac(1,2)*pi./dt);
% bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.06)&&aresame(fit.lambda,h./dt,0.002);
% reporttest('MATERNFIT recovers Matern parameters with non-unit sample rate',allall(bool))

rng(0);
dt=3600;
z=maternoise(dt,1000,sigo,alpha,h./dt);
fit=maternfit(dt,z,frac(1,2)*pi./dt);
bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.06)&&aresame(fit.lambda,h./dt,0.002);
reporttest('MATERNFIT recovers Matern parameters with non-unit sample rate',allall(bool))

rng(0);
epsilono=2;
zn=epsilono.*(randn(length(z),1)+1i.*randn(length(z),1))./sqrt(2);
fit=maternfit(dt,zn,frac(1,2)*pi,'noisy','value.sigma',0);
bool=(abs(fit.sigma)<1e-7)&&aresame(abs(fit.epsilon./epsilono-1),0,0.1);
reporttest('MATERNFIT recovers Matern parameters with non-unit sample rate, noise only',allall(bool))

h=1;
rng(0);
dt=3600;

z=maternoise(dt,1000,sigo,alpha,h./dt);
psi=sleptap(size(z,1),3,1);
fit=maternfit(dt,z,frac(1,2)*pi./dt,'tapered',psi);
bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.06)&&aresame(fit.lambda,h./dt,0.002);
reporttest('MATERNFIT recovers Matern parameters with non-unit sample rate, tapered version',allall(bool))

psi=sleptap(size(z,1)-1,3,1);
tic;
fit=maternfit(dt,z,frac(1,2)*pi./dt,'tapered',psi,'difference');
bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.06)&&aresame(fit.lambda,h./dt,0.002);
reporttest('MATERNFIT recovers Matern parameters with non-unit sample rate, differenced tapered version',allall(bool))
etime1=toc;

sigo=17;
alpha=1.5;
h=1/10;
rng(0);

dt=3600;
gamma=5;
N=1000;
% clf
% z=maternoise(dt,[N,2000],sigo,alpha,h./dt,1,'generalized');
% plot(vmean(abs(fft(z)),2)),xlog,ylog
% z=maternoise(dt,[N,2000],sigo,alpha,h./dt,'chol');
% hold on,plot(vmean(abs(fft(z)),2)),xlog,ylog
z=maternoise(dt,[N,5000],sigo,alpha,h./dt,5,'generalized');
zn=epsilono.*randn(size(z,1),5000);
%hold on,plot(vmean(abs(fft(z)),2))

% [tau,tpz1]=materncov(dt,N,sigo,alpha,h./dt);
% [tau,tpz2]=materncov(dt,N,sigo,alpha,h./dt,1,'generalized');
% [tau,tpz3]=materncov(dt,N,sigo,alpha,h./dt,2,'generalized');
% [tau,tpz4]=materncov(dt,N,sigo,alpha,h./dt,5,'generalized');
% figure,plot(tau,[tpz1 tpz2 tpz3 tpz4])

%[f,s]=maternspec(dt,N,sigo,alpha,h./dt,5,'generalized');plot(f,s./maxmax(s))
%[f,s]=maternspec(dt,N,sigo,alpha,h./dt,1,'generalized');hold on,plot(f,s./maxmax(s))
%[f,s]=maternspec(dt,N,sigo,alpha,h./dt,1/2,'generalized');hold on,plot(f,s./maxmax(s))
%[f,s]=maternspec(dt,N,sigo,alpha,h./dt,1/10,'generalized');hold on,plot(f,s./maxmax(s))

fit=maternfit(dt,z,frac(1,2)*pi./dt,'generalized','average');
%fit=maternfit(dt,z,frac(1,2)*pi./dt,'generalized');
bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.06)&&aresame(fit.lambda,h./dt,0.002)&&aresame(fit.gamma,5,0.25);
reporttest('MATERNFIT recovers generalized Matern parameters with non-unit sample rate',allall(bool))

fit=maternfit(dt,z+zn,frac(1,2)*pi./dt,'generalized','average','noisy');
%fit=maternfit(dt,z,frac(1,2)*pi./dt,'generalized');
bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.06)&&aresame(fit.lambda,h./dt,0.002)&&aresame(fit.gamma,5,0.25)&&aresame(abs(fit.epsilon./epsilono-1),0,0.1);
reporttest('MATERNFIT recovers generalized Matern parameters with non-unit sample rate, noisy version',allall(bool))

% if exist('nlopt_optimize')==3
%     psi=sleptap(size(z,1)-1,3,1);
%     tic;
%     fit=maternfit(dt,z,frac(1,2)*pi./dt,'tapered',psi,'difference','nlopt');
%     bool=aresame(abs(fit.sigma./sigo-1),0,0.2)&&aresame(fit.alpha,alpha,0.06)&&aresame(fit.lambda,h./dt,0.002);
%     reporttest('MATERNFIT using NLopt version of previous',allall(bool))
%     etime2=toc;
%     disp(['MATERNFIT NLopt version took ' num2str(etime2./etime1) ' as much time as FMINSEARCH.'])
% end

% sigo=17;
% alpha=1.5;
% h=1/10;
%
% rng(0);
% dt=1;
% z=maternoise(dt,1000,sigo,alpha,h);
% fit=maternfit(dt,[z z],frac(1,2)*pi);

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function net_inputs = rfcn_get_minibatch(conf, image_roidb)
% net_inputs = rfcn_get_minibatch(conf, image_roidb)
% --------------------------------------------------------
% D&T implementation
% Modified from MATLAB  R-FCN (https://github.com/daijifeng001/R-FCN/)
% and Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2017, Christoph Feichtenhofer
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------
    visualize = false;

    num_images = length(image_roidb);
    if num_images == 0
      net_inputs = NaN; return;
    end
    if conf.bbox_class_agnostic
        num_classes = 1;
    else
        % Infer number of classes from the number of columns in gt_overlaps
        num_classes = size(image_roidb(1).overlap, 2);
    end
    % Sample random scales to use for each image in this batch
    random_scale_inds = randi(length(conf.scales), num_images, 1);

    if conf.batch_size > 0
        assert(mod(conf.batch_size, num_images) == 0, ...
        sprintf('num_images %d must divide BATCH_SIZE %d', num_images, conf.batch_size));

        rois_per_image = conf.batch_size / num_images;
        fg_rois_per_image = round(rois_per_image * conf.fg_fraction);
    else
        rois_per_image = inf;
        fg_rois_per_image = inf;
    end
    if conf.use_flipped
      flip_im = rand(num_images,1) > .5;
      if isfield(conf,'sample_vid') && conf.sample_vid, flip_im = repmat(flip_im(1), size(flip_im)); end;
    else
      flip_im = zeros(num_images,1) ;
    end
    % Get the input image blob
    [im_blob, im_scales] = get_image_blob(conf, image_roidb, random_scale_inds, flip_im);
    if isnan(im_blob), net_inputs = NaN; return; end
    % build the region of interest and label blobs
    rois_blob = zeros(0, 5, 'single');
    labels_blob = zeros(0, 1, 'single');
    bbox_targets_blob = zeros(0, 4 * (num_classes+1), 'single');
    bbox_loss_blob = zeros(size(bbox_targets_blob), 'single');
    track_rois_blob = rois_blob; track_targets_blob = bbox_targets_blob; track_loss_blob = bbox_loss_blob;
    for i = 1:num_images
        [labels, ~, im_rois, bbox_targets, bbox_loss] = ...
            sample_rois(conf, image_roidb(i), fg_rois_per_image, rois_per_image);

        if flip_im(i)
            im_rois(:, [1, 3]) = image_roidb(i).im_size(2) + 1 - im_rois(:, [3, 1]);
            %flip horizontal regression target: targets_dx = (gt_ctr_x - ex_ctr_x) ./ (ex_widths+eps);
            bbox_targets(:, 5) = - bbox_targets(:, 5);
        end

        if visualize
          im = (imread(  image_roidb(i).image_path ));
          boxes_cell = {im_rois(:,:)};

          figure('windowstyle','docked');
          imd(im);
          drawBBoxes(im_rois(:,:), 'LineWidth', 1);
        end

        % Add to ROIs blob
        feat_rois = rfcn_map_im_rois_to_feat_rois(conf, im_rois, im_scales(i));
        batch_ind = i * ones(size(feat_rois, 1), 1);
        rois_blob_this_image = [batch_ind, feat_rois];
        rois_blob = [rois_blob; rois_blob_this_image];

        % Add to labels, bbox targets, and bbox loss blobs
        labels_blob = [labels_blob; labels];
        bbox_targets_blob = [bbox_targets_blob; bbox_targets];
        bbox_loss_blob = [bbox_loss_blob; bbox_loss];

        if isfield(image_roidb(i), 'track_targets') && ~isempty(image_roidb(i).track_targets)
          [track_targets, track_loss_weights] = get_bbox_regression_labels(conf,  image_roidb(i).track_targets, num_classes);

          if flip_im(i)
             image_roidb(i).track_rois(:, [1, 3])  = image_roidb(i).im_size(2) + 1 - image_roidb(i).track_rois(:, [3, 1]);
             track_targets(:, 5) = - track_targets(:, 5);
          end
            track_rois_blob = [track_rois_blob; [ones(size(image_roidb(i).track_rois, 1), 1), rfcn_map_im_rois_to_feat_rois(conf, image_roidb(i).track_rois, im_scales(i))]];
            track_targets_blob = [track_targets_blob; track_targets];
            track_loss_blob = [track_loss_blob; track_loss_weights];
        end
    end

    % permute data into caffe c++ memory, thus [num, channels, height, width]
    im_blob(:, :, 1:3, :) = im_blob(:, :, [3, 2, 1], :); % from rgb to brg
    if size(im_blob,3) > 3
      im_blob(:, :, 4:6, :) = im_blob(:, :, [6, 5, 4], :); % from rgb to brg
    end
    im_blob = single(permute(im_blob, [2, 1, 3, 4]));
    rois_blob = rois_blob - 1; % to c's index (start from 0)
    rois_blob = single(permute(rois_blob, [3, 4, 2, 1]));
    labels_blob = single(permute(labels_blob, [3, 4, 2, 1]));
    bbox_targets_blob = single(permute(bbox_targets_blob, [3, 4, 2, 1]));
    bbox_loss_blob = single(permute(bbox_loss_blob, [3, 4, 2, 1]));

    track_rois_blob = track_rois_blob -1;
    track_rois_blob = single(permute(track_rois_blob, [3, 4, 2, 1]));
    track_targets_blob = single(permute(track_targets_blob, [3, 4, 2, 1]));
    track_loss_blob = single(permute(track_loss_blob, [3, 4, 2, 1]));

    try
    assert(~isempty(im_blob));
    assert(~isempty(rois_blob));
    assert(~isempty(labels_blob));
    assert(~isempty(bbox_targets_blob));
    assert(~isempty(bbox_loss_blob));
    catch
      fprintf('no object in image %s\n', image_roidb(1).image_path );
       net_inputs = NaN; return;
    end
    if ~isempty(track_rois_blob)
      net_inputs = {im_blob, rois_blob, labels_blob, bbox_targets_blob, bbox_loss_blob, ...
        track_rois_blob, track_targets_blob, track_loss_blob};
    else
      net_inputs = {im_blob, rois_blob, labels_blob, bbox_targets_blob, bbox_loss_blob};
    end

end

%% Build an input blob from the images in the roidb at the specified scales.
function [im_blob, im_scales] = get_image_blob(conf, images, random_scale_inds, flip_im)

    num_images = length(images);
    processed_ims = cell(num_images, 1);
    im_scales = nan(num_images, 1);
    for i = 1:num_images
      try
        if ~ispc
          images(i).image_path = strrep(images(i).image_path,'\','/');
        end
        nFrames = 1;
       if isfield(conf, 'nFrames')
         nFrames = conf.nFrames;
       end
        if isfield(conf, 'input_modality') && ~strcmp(conf.input_modality, 'rgb')
          if strcmp(conf.input_modality, '2framergb') || strcmp(conf.input_modality, 'framebatch')
            time_stride = 10;
             if isfield(conf, 'time_stride')
               time_stride = conf.time_stride;
             end
            [vid, frame, ext] = fileparts(images(i).image_path);
            nextframe = fullfile(vid, [sprintf('%06d',str2num(frame)+time_stride) ext]) ;
            frames = str2num(frame)-nFrames/2*time_stride:time_stride:str2num(frame)+nFrames/2*time_stride;
            im = {};
            for frame=frames
              nextframe = fullfile(vid, [sprintf('%06d',frame) ext]) ;
              im{end+1} = imread(nextframe);
            end
            if strcmp(conf.input_modality, 'framebatch')
              im = cat(4,im{:});
            else
              im = cat(3,im{:});
            end

          else
            if strcmp(images(i).imdb_name(1:12), 'ilsvrc15_vid')
              img_file_u = strrep(images(i).image_path, ['Data' filesep 'VID'], ...
                ['Data' filesep 'VID' filesep 'tvl1_flow_600' filesep 'u'] ) ;
              img_file_v = strrep(images(i).image_path, ['Data' filesep 'VID'], ...
                ['Data' filesep 'VID' filesep 'tvl1_flow_600' filesep 'v'] ) ;
              try
                [path_u, ~, ~] = fileparts(img_file_u);
                [path_v, frame, ext] = fileparts(img_file_v);
                frames = str2num(frame)-nFrames+1:str2num(frame)+nFrames-1;
                frames(frames <0) = 0; frames(frames > imdb.num_frames(i)) = imdb.num_frames(i);
                im_u = []; im_v = [];
                for frame=frames
                  img_file_u = fullfile(path_u, [sprintf('%06d',frame), ext]) ;
                  img_file_v = fullfile(path_v, [sprintf('%06d',frame), ext]) ;
                  im_u{end+1} = imread(img_file_u); im_v{end+1} = imread(img_file_v);
                end
                im_u = cat(3,im_u{:});im_v = cat(3,im_v{:});
              end
            else
              im_u = 128*ones(images(i).im_size); im_v = 128*ones(images(i).im_size);
            end
            sz = size(im_u);
          end
          switch conf.input_modality
            case 'flow'
              flow = single(cat(3,im_u,im_v)) - 128;
              flow = bsxfun(@minus,flow, median(median(flow,1),2)) ;
              mag_flow = sqrt(sum(flow.^2,3)) - 128;
              im = (cat(3,flow,mag_flow));
              im = imresize( im, images(i).im_size);
            case 'grayflow'
              imrgb = imread(images(i).image_path); sz = size(imrgb);
              imgray = rgb2gray(imrgb);
              flow = single(cat(3,im_u,im_v)) ;
              im = cat(3,imresize( flow, sz(1:2)),imgray);
            case 'rgbflow'
            if isfield(conf, 'subflow_mean') && ~conf.subflow_mean
              flow = single(cat(3,im_u,im_v));
            else
              flow = single(cat(3,im_u,im_v)) - 128;
            end
            imrgb = single(imread(images(i).image_path));

            if size(imrgb,3) < 3
              imrgb = imrgb(:,:,[1 1 1]);
            end
            sz = size(imrgb);
            im = cat(3,imrgb,imresize( flow, sz(1:2)));
          end
        else
            im = imread(images(i).image_path);
        end
      catch
        fprintf('could not read %s\n', images(i).image_path);
        im_blob = NaN;
        return;
      end
        target_size = conf.scales(random_scale_inds(i));

        if flip_im(i)
            im = fliplr(im);
        end
        [im, im_scale] = prep_im_for_blob(im, conf.image_means, target_size, conf.max_size);
        if isfield(conf,'image_std')
          im = bsxfun(@rdivide, im, conf.image_std);
        end
        im_scales(i) = im_scale;
        processed_ims{i} = im;
    end

    im_blob = im_list_to_blob(processed_ims);
end

%% Generate a random sample of ROIs comprising foreground and background examples.
function [labels, overlaps, rois, bbox_targets, bbox_loss_weights] = sample_rois(conf, image_roidb, fg_rois_per_image, rois_per_image)

    [overlaps, labels] = max(image_roidb(1).overlap, [], 2);
    rois = image_roidb(1).boxes;

    % Select foreground ROIs as those with >= FG_THRESH overlap
    fg_inds = find(overlaps >= conf.fg_thresh);
    % Guard against the case when an image has fewer than fg_rois_per_image
    % foreground ROIs
    fg_rois_per_this_image = min(fg_rois_per_image, length(fg_inds));
    % Sample foreground regions without replacement
    if ~isempty(fg_inds)
       fg_inds = fg_inds(randperm(length(fg_inds), fg_rois_per_this_image));
    end

    % Select background ROIs as those within [BG_THRESH_LO, BG_THRESH_HI)
    bg_inds = find(overlaps < conf.bg_thresh_hi & overlaps >= conf.bg_thresh_lo);
    % Compute number of background ROIs to take from this image (guarding
    % against there being fewer than desired)
    bg_rois_per_this_image = rois_per_image - fg_rois_per_this_image;
    bg_rois_per_this_image = min(bg_rois_per_this_image, length(bg_inds));
    % Sample foreground regions without replacement
    if ~isempty(bg_inds)
       bg_inds = bg_inds(randperm(length(bg_inds), bg_rois_per_this_image));
    end
    % The indices that we're selecting (both fg and bg)
    keep_inds = [fg_inds; bg_inds];
    % Select sampled values from various arrays
    labels = labels(keep_inds);
    % Clamp labels for the background ROIs to 0
    labels((fg_rois_per_this_image+1):end) = 0;
    overlaps = overlaps(keep_inds);
    rois = rois(keep_inds, :);

    if conf.bbox_class_agnostic
        assert(all((labels>0) == image_roidb.bbox_targets(keep_inds, 1)));
    else
        assert(all(labels == image_roidb.bbox_targets(keep_inds, 1)));
    end

    % Infer number of classes from the number of columns in gt_overlaps
    num_classes = size(image_roidb(1).overlap, 2);

    [bbox_targets, bbox_loss_weights] = get_bbox_regression_labels(conf, ...
        image_roidb.bbox_targets(keep_inds, :), num_classes);

end

function [bbox_targets, bbox_loss_weights] = get_bbox_regression_labels(conf, bbox_target_data, num_classes)
%% Bounding-box regression targets are stored in a compact form in the roidb.
 % This function expands those targets into the 4-of-4*(num_classes+1) representation used
 % by the network (i.e. only one class has non-zero targets).
 % The loss weights are similarly expanded.
% Return (N, (num_classes+1) * 4, 1, 1) blob of regression targets
% Return (N, (num_classes+1 * 4, 1, 1) blob of loss weights
    if conf.bbox_class_agnostic
        num_classes = 1;
    end

    clss = bbox_target_data(:, 1);
    bbox_targets = zeros(length(clss), 4 * (num_classes+1), 'single');
    bbox_loss_weights = zeros(size(bbox_targets), 'single');
    inds = find(clss > 0);
    for i = 1:length(inds)
       ind = inds(i);
       cls = clss(ind);
       targets_inds = (1+cls*4):((cls+1)*4);
       bbox_targets(ind, targets_inds) = bbox_target_data(ind, 2:end);
       bbox_loss_weights(ind, targets_inds) = 1;
    end
end

author: feichtenhofer

repo: Detect-Track


text:
function [VV,TT,FF,TN,IFF] = cdt(varargin)
  % CDT compute the constrained delaunay triangulation of a given mesh using
  % tetgen.
  %
  % [VV,TT,FF,TN,IFF] = cdt(V,F,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   V  #V by dim list of mesh positions
  %   F  #F by dim list of face indices
  %   Optional:
  %     'TriangleFlags' followed by a string of additional triangle flags, e.g.
  %       '-q'
  %     'TetgenFlags' followed by a string of additional tetgen flags, e.g.
  %       '-J'
  %     'UseBoundingBox' followed by whether to use an explicit bounding box
  %       rather than tetgen's convex hull
  %     'BoundingBoxPush' followed by a push factor to enlargen the bounding
  %       box around its centroid {1.01}
  %     'BoundingBoxUpsample' followed by a factor f, to upsample the bounding
  %       box faces until min(doublearea(BV,BF))<f*max(doublearea(V,F)) {inf}
  %     'Quiet' followed by true or {false}
  % Outputs:
  %   VV  #VV by dim list of vertices over which TT is defined, should
  %     contain V as a prefix and should be exactly V unless F self-intersects
  %   TT  #TT by dim+1 list of tetrahedra indices
  %   FF  #FF by dim list of F potentially subdivided to be defined over VV,
  %     should be exactly F unless F self-intersects
  %   TN  #TT by dim+1 list of tetrahedra neighbors
  %   IFF  #FF by 1 list of original facets indices corresponding to each facet
  %     in FF
  %

  % parse input
  V = varargin{1};
  if nargin>=2
    F = varargin{2};
  else
    F = [];
  end
  % default values for optional input
  tetgen_flags = '';
  triangle_flags = '';
  use_bounding_box = 0;
  bb_push = 1.01;
  bb_ups = inf;
  quiet = 0;

  % parse any optional input
  ii = 3;
  while(ii <= nargin)
    switch varargin{ii}
    case 'TriangleFlags'
      ii = ii + 1;
      assert(ii<=nargin);
      triangle_flags = varargin{ii};
    case 'TetgenFlags'
      ii = ii + 1;
      assert(ii<=nargin);
      tetgen_flags = varargin{ii};
    case 'UseBoundingBox'
      ii = ii+1;
      assert(ii<=nargin);
      use_bounding_box = 1*(varargin{ii});
    case 'BoundingBoxPush'
      ii = ii+1;
      assert(ii<=nargin);
      bb_push = varargin{ii};
    case 'BoundingBoxUpsample'
      ii = ii+1;
      assert(ii<=nargin);
      bb_ups = varargin{ii};
    case 'Quiet'
      ii = ii+1;
      assert(ii<=nargin);
      quiet = varargin{ii};
      quiet = quiet * 1;
    otherwise
      error(['Unsupported parameter: ' varargin{ii}]);
    end
    ii = ii + 1;
  end

  if quiet
    fid = fopen('/dev/null','w');
  else
    fid = 1;
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Check for cached result, do NOT edit variables until cache is checked,
  % your function code comes later. See below
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  clear varargin;
  % get a list of current variables in this scope, this is the input "state"
  variables = who;
  % get a temporary file's name
  tmpf = [tempname('.') '.mat'];
  % save the "state" to file, so we can get a md5 checksum
  save(tmpf,'-regexp',sprintf('^%s$|',variables{:}),'-ascii');
  % get md5 checksum on input "state", we append .cache.mat to the check sum
  % because we'll use the checksum as the cache file name
  [s,cache_name] = system(['/sbin/md5 -r ' tmpf ' | awk ''{printf "."$1".cache.mat"}''']);
  % clean up
  delete(tmpf);
  clear s tmpf variables;
  % If the checksum cache file exists then we've seen this input "state"
  % before, load in cached output "state"
  if(exist(cache_name,'file'))
    fprintf_(fid,'Cache exists. Using cache...\n');
    % use cache
    load(cache_name);
  % Otherwise this is the first time we've seen this input "state", so we
  % execute the function as usual and save the output "state" to the cache
  else
    fprintf_(fid,'First time. Creating cache...\n');

    dim = size(V,2);
    assert(numel(F) ==0 || dim >= size(F,2), sprintf( ...
      'dimension of V (%d) < simplex size (%d)',dim,size(F,2)));
    switch dim
    case 2
      % we're dealing with edges
      E = F;
      if use_bounding_box
        [BV,BF] = bounding_box([max(V)+eps;min(V)-eps]);
        dir = [sign(bsxfun(@minus,BV,mean(BV)))];
        mag_max = max(sqrt(sum(bsxfun(@minus,BV,mean(BV)).^2,2)));
        BV = bsxfun(@plus,bb_push*mag_max*dir,mean(BV));
        BE = [];
        EmBE = E;
      else
        BV = [];
        BF = [];
        % explicitly add convex hull
        F = delaunay(V);
        BE = outline(F);
        EmBE = setdiff(sort(E,2),sort(BE,2),'rows');
      end
      [VV,TT] = triangle([V;BV],[BE;EmBE;size(V,1)+BF],[],'Flags',triangle_flags);

      % This is a bullshit hack way of determining which edges were subdivided
      % A better way is using boundary markers in triangle like in
      % diffcurves_test.m
      if nargout>=3
      % First determine which edges were not subdivided
      E_indices = 1:size(E,1);
      [maintained] = ismember(sort(E,2),sort(edges(TT),2),'rows');
      EE = E(maintained,:);
      IEE = E_indices(maintained);
      % and those that were
      ES = E(~maintained,:);
      IES = E_indices(~maintained);
      % hopefully #ES is small, let Alec know if it isn't, we should rewrite this
      if size(ES,1)>100
        warning('many (%d) edges to subdivide, this could be slow',size(ES,1));
      end
      % build edge-length adjacency matrix
      C = adjacency_edge_cost_matrix(VV,edges(TT));
      % average edge length for tolerances
      h = avgedge(VV,TT);
      % loop over edges that were subdivided
      for e = 1:size(ES,1)
        % get shortest path from source to dest
        [D,P] = dijk(C,ES(e,2),ES(e,1));
        % check that D is what a straight line path should be

        orig_len = sqrt(sum((V(ES(e,2),:)-V(ES(e,1),:)).^2,2));
        if abs(orig_len-D)>1e-13*h
          warning('path %d to %d may not be straight: |%g - %g| = %g', ...
            ES(e,2),ES(e,1),orig_len,D,abs(orig_len-D));
        end
        d = ES(e,1);
        NE = [];
        % trace back path
        while d~=ES(e,2)
          NE = [NE;d P(d)];
          d = P(d);
        end
        % append to edges
        EE = [EE;NE];
        IEE = [IEE(:);repmat(IES(e),size(NE,1),1)];
      end
      % rename to match output
      FF = EE;
      end
      if nargout >= 5
        IFF = IEE;
      end
      if nargout >= 4
        % Compute element neighbors
        TN = triangle_triangle_adjacency(TT);
      end
    case 3

      % It's really important that we have at least tetgen 1.5
      [status,result] = system([path_to_tetgen ' -help']);
      tetgen_version = ...
        str2double(regexprep(result,'.*Version ([0-9].[0-9]).*','$1'));
      if isnan(tetgen_version) || tetgen_version < 1.5
        error('%s version (%g) < 1.5',path_to_tetgen,tetgen_version);
      end

      FF = F;

      % Q: Should clean up be the caller's responsibility?
      % A: Yes, I think so. We can parse tetgen and give hints.

      %% number of input points
      %n = size(V,1);

      %% Deal with duplicate input points
      %[V,I,J] = remove_duplicate_vertices(V);
      %if n ~= size(V,1)
      %  warning('cdt: %d exactly duplicate points removed from input',n-size(V,1));
      %  n = size(V,1);
      %  % Remap faces
      %  FF = J(F);
      %end

      %if ~isempty(F)
      %  NF = F( (F(:,1) ~= F(:,2)) & (F(:,2) ~= F(:,3)) & (F(:,3) ~= F(:,1)),:);
      %  if size(NF,1) < size(F,1)
      %    warning('cdt: %d faces with combinatorally duplicate vertices removed\n', ...
      %      size(F,1)-size(NF,1));
      %    FF = NF;
      %  end
      %end

      % Try to mesh with all faces included directly
      BM = -ones(size(FF,1),1);

      if use_bounding_box
        [BV,BF] = bounding_box(V);
        % Subdivide
        if ~isempty(FF) && size(FF,2)==3
          max_dblA = max(doublearea(V,FF));
          while min(doublearea(BV,BF))> bb_ups*max_dblA
            [BV,BF] = upsample(BV,BF);
          end
        end
        % push a little a away
        BV = bsxfun(@plus,bb_push*bsxfun(@minus,BV,mean(BV)),mean(BV));
        BF = BF + size(V,1);
        V = [V;BV];
        FF = [FF;BF];
        % true boundary marker
        BM = [BM;ones(size(BF,1),1)];
      end

      prefix = tempname;
      poly_filename = [prefix '.poly'];
      writePOLY_tetgen(poly_filename,V,FF,[],'BoundaryMarkers',BM);

      % -p: mesh PLC
      % -C: check consistency of result
      % -c: retain convex hull, doesn't seem to play nicely with -q
      %% -O0: optimization level
      %% -Y: no steiner points on boundary
      %% -J: Do not jettison duplicate vertices (we should have already handled
      %%   these)
      % -n: output neighbor information
      %% -F: suppress output of face information
      %% -S0: specifies maximum number of added points (may still add schoenhardt
      %   points)
      % There seems to be a secret -o[val] flag that optimizes the max dihedral
      % angle
      delaunay_flags = '-pCn';
      if ~use_bounding_box
        % then retain convex hull. Note: this seems broken, hence the option to
        % use bounding box
        delaunay_flags = [delaunay_flags 'c'];
      end
      %delaunay_flags = '-pq0.75 -a0.0002 -O5YcJnFg';

      command = [path_to_tetgen ' ' delaunay_flags ' ' tetgen_flags ' ' poly_filename];
      fprintf_(fid,'%s\n',command);
      if quiet
        [status, result] = system(command);
      else
        [status, result] = system(command,'-echo');
      end
      fprintf_(fid,'%s',result);
      if status ~= 0
        fprintf_(fid,'%s',result);
      else
        % Print some warnings based on tetgen output
        if strfind(result,'Steiner points on boundary edges')
          fprintf_(fid,'%s',result);
          warning('Tetgen added steiner points on boundary');
        end
        if strfind(result,'Jettisoning redundants points.')
          fprintf_(fid,'%s',result);
          warning('Tetgen removed duplicate points');
        end

      end

      % See terminatetetgen() in tetgen.h
      switch(status)
      case 0
      case 3
        error('Self intersections. Try selfintersect() or remove self-intersections');
      %% Self-intersection
      %  warning('Throwing away self-intersecting facets');
      %  % Seems tetgen sets status=3 if self-intersecting
      %  % Rerun and detect self-intersecting faces
      %  % -p: mesh PLC, -d: detect all self intersections, -J: don't jettison
      %  % duplicate points (since we want self-intersecting faces defined over
      %  % original vertices)
      %  selfintersection_flags = '-pdJ ';
      %  command = [path_to_tetgen ' ' selfintersection_flags ' ' poly_filename];
      %  fprintf_(fid,command);
      %  [status, result] = system(command);
      %  if strfind(result,'No faces are intersecting.')
      %    fprintf_(fid,'%s',result);
      %    error('Tetgen unsure whether there are self-intersections');
      %  end
      %  assert(status==0);
      %  % tetgen always writes output to file:
      %  face_filename = [prefix '.1.face'];
      %  node_filename = [prefix '.1.node'];
      %  % faces participating in self-intersections
      %  % Faces will index new nodes, but they better be the input ones
      %  IF = readFACE(face_filename);
      %  % clean up temp files
      %  delete(face_filename);
      %  %VV = readNODE(node_filename);
      %  %assert(size(VV,1)==size(V,1));
      %  %delete(node_filename);

      %  % Segments aren't supported by tetgen the way I'd like them to be (i.e. I
      %  % get seg faults)
      %  %% Get edges of intersecting faces
      %  %IFE = edges(IF);
      %  %IFsegs = mat2cell(IFE,ones(size(IFE,1),1),[2]);

      %  % get orignal faces \ intersecting faces
      %  % Using sort like this assumes that tetgen ignores orientation
      %  FmIF = setdiff(sort(FF,2),sort(IF,2),'rows');
      %  Facets = [];
      %  Facets.facets = mat2cell(FmIF,ones(size(FmIF,1),1),[size(FmIF,2)]);
      %  Facets.boundary_marker = -ones(numel(Facets.facets),1);
      %  Facets.holes = cell(numel(Facets.facets),1);

      %  % Rewrite poly file
      %  writePOLY(poly_filename,V,Facets,[]);

      %  % shoot for the moon
      %  command = [path_to_tetgen ' ' delaunay_flags ' ' tetgen_flags ' ' poly_filename];
      %  fprintf_(fid,command);
      %  [status, result] = system(command);
      %  result
      %  status
      case 4
      % Small feature
        error('Try again with reduced T?');
      case 5
      % close input facets
        error('Try again with -Y');
      case 6
      % close input facets
        error('Input error. Check yo` input');
      otherwise
        error('Tetgen returned status %d != 0',status);
      end

      % tetgen always writes output to file:
      %   xxxx.1.ele  tetrahedra
      %   xxxx.1.node tetrahedra vertices
      ele_filename = [prefix '.1.ele'];
      face_filename = [prefix '.1.face'];
      node_filename = [prefix '.1.node'];
      neigh_filename = [prefix '.1.neigh'];

      %delete(poly_filename);

      VV = readNODE(node_filename);
      TT = readELE(ele_filename);
      TN = readNEIGH(neigh_filename);
      TN(TN>=0) = TN(TN>=0)-1;

      [FF,BF] = readFACE(face_filename);
      %FF = FF-1;
      FF = FF(BF==-1,:);

      %[VV1,TT1] = readMESH(mesh_filename);
      % TODO: handle duplicate input vertices in a reasonable way
      % TODO: Tetgen seems to add steiners when input has coplanar neighboring
      % triangles (aka triangulated planar quads). It seems its enough to do:
      % V = V+eps*rand(size(V));
      if ~quiet
        if size(VV,1) > size(V,1)
          warning('New mesh has extra %d vertices',size(VV,1)-size(V,1));
        elseif size(VV,1) < size(V,1)
          warning('New mesh has %d fewer vertices',size(V,1)-size(VV,1));
        elseif max(abs(VV-V))>eps
          warning('max(abs(VV-V)) = %g > %g',max(abs(VV-V)),eps);
          %% Remap to original vertex list
          %TT = I(TT);
          %VV = V;
        end
      end
      %%% Reverse order so that boundary_faces(TT) meets CCW convention
      %TT = TT(:,[3 2 1 4]);
      %TN = TN(:,[3 2 1 4]);

      %% remove bounding box facets
      %if use_bounding_box
      %  old_FF_count = size(FF,1);
      %  FF = ...
      %    FF(all(reshape(~ismember(VV(FF(:),:),BV,'rows'),size(FF,1),3),2),:);
      %  rm_count = old_FF_count-size(FF,1);
      %  switch rm_count
      %  case 0:11
      %    warning('Only removed %d boundary facets (instead of 12)');
      %  case 12
      %    % alles gut
      %  otherwise
      %    warning('Removed %d boundary facets (instead of 12)',rm_count);
      %  end
      %end
    otherwise
      error(sprintf('Unsupported dimension (%d)',dim));
    end

    if nargout>=3
    [IN] = in_elements(FF,TT);
    if ~all(IN)
       warning('%d faces not in elements',sum(~IN));
    end
    end

    % get list of variables present in this scope at finish of function code,
    % this is the output "state"
    variables = who;
    % save output "state" to file, using md5 checksum cache file name
    save(cache_name,'-regexp',sprintf('^%s$|',variables{:}));

  end
end

%% local fprint variant for cross-latform functionality
function fprintf_(fid,varargin)
    if fid < 0
        % on windows there is no dev/null
        return
    else
        fprintf(fid,varargin{:});
    end
end

author: alecjacobson

repo: gptoolbox


text:
function y = prtUtilNanVar(x,dim)
% prtUtilNanVar - Calculated the variance of data, X, ignoring nans
%
% y = prtUtilNanVar(x)
% y = prtUtilNanVar(x,dim)

if nargin < 2 || isempty(dim)
    dim = 1; % Default of sum
    if isvector(x)
        if size(x,1) > size(x,2)
            dim = 1;
        else
            dim = 2;
        end
    end
end

meanX = prtUtilNanMean(x,dim);
y = prtUtilNanMean(bsxfun(@minus,x,meanX).^2,dim)*size(x,dim)./(size(x,dim)-1); % Variance is normalized by n-1 not n;

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function models = esvm_initialize_exemplars(e_set, params, ...
                                            models_name)

% Initialize script which writes out initial model files for all
% exemplars in an exemplar stream e_set (see get_pascal_stream)
% NOTE: this function is parallelizable (and dalalizable!)
%
% INPUTS:
% params.dataset_params: the parameters of the current dataset
% e_set: the exemplar stream set which contains
%   e_set{i}.I, e_set{i}.cls, e_set{i}.objectid, e_set{i}.bbox
% models_name: model name
% init_params: a structure of initialization parameters
% init_params.init_function: a function which takes as input (I,bbox,params)
%   and returns a model structure [if not specified, then just dump
%   out names of resulting files]
%
% OUTPUTS:
% allfiles: The names of all outputs (which are .mat model files
%   containing the initialized exemplars)
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

% DTstring = '';
% if dalalmode == 1
%   %Find the best window size from taking statistics over all
%   %training instances of matching class
%   hg_size = get_hg_size(cls);
%   DTstring = '-dt';
% elseif dalalmode == 2
%   hg_size = [8 8];
%   DTstring = '-gt';
% end

%if (dalalmode == 1) || (dalalmode == 2)
%  %Do the dalal-triggs anisotropic warping initialization
%  model = initialize_model_dt(I,bbox,SBIN,hg_size);
%else

if isfield(params,'dataset_params') && ...
      isfield(params.dataset_params,'localdir') && ...
      length(params.dataset_params.localdir)>0
  CACHE_FILE = 1;
else
  CACHE_FILE = 0;
  params.dataset_params.localdir = '';
end

if ~exist('models_name','var')
  models_name = '';
end

cache_dir =  ...
    sprintf('%s/models/',params.dataset_params.localdir);

cache_file = ...
    sprintf('%s/%s.mat',cache_dir,models_name);

if CACHE_FILE ==1 && fileexists(cache_file)
  models = load(cache_file);
  models = models.models;
  return;
end

results_directory = ...
    sprintf('%s/models/%s/',params.dataset_params.localdir, ...
            models_name);

if CACHE_FILE==1 && ~exist(results_directory,'dir')
  fprintf(1,'Making directory %s\n',results_directory);
  mkdir(results_directory);
end

%fprintf(1,'Writing %d Exemplars to directory %s\n',length(e_set),...
%        results_directory);

%Randomize creation order
if CACHE_FILE == 1
  myRandomize;
  rrr = randperm(length(e_set));

  if (params.dataset_params.display == 1)
    rrr = 1:length(rrr);
  end
  e_set = e_set(rrr);
end

%Create an array of all final file names
allfiles = cell(length(e_set), 1);

for i = 1:length(e_set)
  cls = e_set{i}.cls;
  objectid = e_set{i}.objectid;
  bbox = e_set{i}.bbox;
  curid = e_set{i}.curid;

  %[tmp,curid,tmp] = fileparts(e_set{i}.I);
  %filer = sprintf('%s/%s.%d.%s.mat',...
  %                results_directory, curid, objectid, cls);

  filer = sprintf('%s/%s',results_directory, e_set{i}.filer);

  allfiles{i} = filer;
  if ~isfield(params,'init_params')
    error('Warning, cannot initialize without params.init_params\n');
  end
  filerlock = [filer '.lock'];

  if CACHE_FILE == 1
    if fileexists(filer) || (mymkdir_dist(filerlock)==0)
      continue
    end
  end
  gt_box = bbox;
  fprintf(1,'.');

  I = convert_to_I(e_set{i}.I);

  %Call the init function which is a mapping from (I,bbox) to (model)
  [model] = params.init_params.init_function(I, bbox, params.init_params);

  clear m
  m.model = model;

  %Save filename (or original image)
  m.I = e_set{i}.I;
  m.curid = curid;
  m.objectid = objectid;
  m.cls = cls;
  m.gt_box = gt_box;

  m.sizeI = size(I);
  m.models_name = models_name;
  m.name = sprintf('%s.%d.%s',m.curid,m.objectid,m.cls);

  if CACHE_FILE == 1
    save(filer,'m');
    if exist(filerlock,'dir')
      rmdir(filerlock);
    end
  else
    allfiles{i} = m;
  end

  % %Print the bounding box overlap between the initial window and
  % %the final window
  % finalos = getosmatrix_bb(m.gt_box, m.model.bb(1,:));
  % fprintf(1,'Final OS is %.3f\n', ...
  %         finalos);

  % fprintf(1,'Final hg_size is %d %d\n',...
  %         m.model.hg_size(1), m.model.hg_size(2));

  %Show the initialized exemplars
  if params.dataset_params.display == 1
    esvm_show_exemplar_frames({m}, 1, params.dataset_params);
    drawnow
    snapnow;
  end
end

if CACHE_FILE == 0
  models = allfiles;
  return;
end

%sort files so they are in alphabetical order
[allfiles, bb] = sort(allfiles);

%Load all of the initialized exemplars
CACHE_FILE = 1;

if length(models_name) == 0
  CACHE_FILE = 0;
end

STRIP_FILE = 0;
DELETE_INITIAL = 0;

models = esvm_load_models(params.dataset_params, models_name, allfiles, ...
                          CACHE_FILE, STRIP_FILE, DELETE_INITIAL);

fprintf(1,'\n   --- Done initializing %d exemplars\n',length(e_set));

author: quantombone

repo: exemplarsvm


text:
%% spatial
src = kthtips_src;
options.J = 4;
options.L = 6;
options.parallel = 0;
w = wavelet_factory_2d_spatial(options, options);
features{1} = @(x)(mean(mean(format_scat(scat(x,w)),2),3));
db = prepare_database(src, features, options);

%% spatial with renorm
src = kthtips_src;
options.J = 4;
options.L = 6;
options.parallel = 0;
w = wavelet_factory_2d_spatial(options, options);
features{1} = @(x)(mean(mean(format_scat(renorm_scat_spatial(scat(x,w))),2),3));
db = prepare_database(src, features, options);

%%
db2 = db_pca(db, 100);

%% without order 0
db2 = db;
db2.features = db2.features(2:end, :);

%% this takes about an hour on 2.4 Ghz core i7
options.parallel = 0;
db = prepare_database(src, features, options);
options.J = 5;
options.antialiasing = 0;
w = wavelet_factory_3d([480, 640], options);
features{1} = @(x)(sum(sum(format_scat(scat(x,w)),2),3));

%% this takes about two hour on 2.4 Ghz core i7
options.parallel = 0;
db2 = prepare_database(src, features, options);

%% not averaged
src = uiuc_src;
options.parallel = 0;
options.J = 5;
options.L = 6;
w = wavelet_factory_2d_spatial(options, options);
features{1} = @(x)(scat(x,w));
db = prepare_database(src, features, options);

%% joint with slog
db3 = db2;
db3.features = log(db2.features);

%% classif with 200 randomn partition and size 5 10 20
grid_train = [5,20,40];
n_fold = 10;
clear error_2d;
for i_fold = 1:n_fold
	for i_grid = 1:numel(grid_train)
		n_train = grid_train(i_grid);
		prop = n_train/81;
		[train_set, test_set] = create_partition(src, prop);
		train_opt.dim = n_train;
		model = affine_train(db, train_set, train_opt);
		labels = affine_test(db, model, test_set);
		error_2d(i_fold, i_grid) = classif_err(labels, test_set, src);
		fprintf('fold %d n_train %g acc %g \n',i_fold, n_train, 1-error_2d(i_fold, i_grid));
	end
end
%  0.5278    0.6775    0.8379
%%
grid_train = [5,20,40];
n_fold = 10;
clear error_2d;
for i_fold = 1:n_fold
	for i_grid = 1:numel(grid_train)
		n_train = grid_train(i_grid);
		prop = n_train/81;
		[train_set, test_set] = create_partition(src, prop);
		train_opt.dim = n_train;
		model = affine_train(db2, train_set, train_opt);
		labels = affine_test(db2, model, test_set);
		error_2d(i_fold, i_grid) = classif_err(labels, test_set, src);
		fprintf('fold %d n_train %g acc %g \n',i_fold, n_train, 1-error_2d(i_fold, i_grid));
	end
end
%  0.5278    0.6775    0.8379

%%
grid_train = [5,10,20];
n_fold = 1;
clear error_2d;

for i_fold = 1:n_fold
	for i_grid = 1:numel(grid_train)
		n_train = grid_train(i_grid);
		prop = n_train/81;
		[train_set, test_set] = create_partition(src, prop);
		train_opt.dim = n_train;
		model = affine_train(db2, train_set, train_opt);
		labels = affine_test(db2, model, test_set);
		error_2d(i_fold, i_grid) = classif_err(labels, test_set, src);
		fprintf('fold %d n_train %g acc %g \n',i_fold, n_train, 1-error_2d(i_fold, i_grid));
	end
end

%% joint scatt
for i = 1:100
	[train_set, test_set] = create_partition(src, 1/2);
	train_opt.dim = 20;
	model = affine_train(db2, ...
		train_set, train_opt);

	labels = affine_test(db2, model, ...
		test_set);
	err_3d(i) = classif_err(labels, test_set, src);
	fprintf('fold %d error %d \n',i, err_3d(i));
end

for i = 1:100
	[train_set, test_set] = create_partition(src, 1/8);
	train_opt.dim = 5;
	model = affine_train(db3, ...
		train_set, train_opt);

	labels = affine_test(db3, model, ...
		test_set);
	err_3d_log(i) = classif_err(labels, test_set, src);
	fprintf('fold %d correct %d \n',i,100*(1- err_3d_log(i)));
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet