text:
% DEMUNLABELLEDIVM2 Test IVM code on a toy crescent data.
%
% Recreates the toy crescent data example shown in the NIPS paper.

% IVM

randn('seed', 1e6)
rand('seed', 1e6)

% Generate a toy data-set
dataSetName = 'unlabelledOne';
experimentNo = 2;

labelledProb = 0.1;
[X, y] = mapLoadData(['semi:' dataSetName ':' num2str(labelledProb)], 1e6);

ind = find(isnan(y));
X(ind, :) = [];
y(ind, :) = [];

options = ivmOptions;
options.noise = 'probit';
options.kern = {'rbf', 'white'};
options.display = 2;
options.numActive = 100;
prior = 0;

% Initialise the model.
model = ivmCreate(size(X, 1), size(y, 2), X, y, options);

prior.type = 'gamma';
prior = priorParamInit(prior);
prior.a = 1;
prior.b = 1;
prior.index = 2;
model.kern.comp{1}.priors(1) = prior;
prior.index = 1;
model.kern.comp{2}.priors(1) = prior;
if options.display > 1
  ivm3dPlot(model, 'ncnmContour', i); %incnmTwoDPlot(model, i);
end
for i = 1:15

  % Plot the data.
  % Select the active set.
  model = ivmOptimiseIvm(model, options.display);
  if options.display > 1
    ivm3dPlot(model, 'ncnmContour', i); %incnmTwoDPlot(model, i);
  end
  % Optimise the kernel parameters.
  model = ivmOptimiseKernel(model, options.display, options.kernIters);
  ivmDisplay(model);

end
model = ivmOptimiseIvm(model, options.display);
if options.display > 1
  ivm3dPlot(model, 'ncnmContour', i); %incnmTwoDPlot(model, i);
end
model = ivmOptimiseIvm(model, options.display);
% Display the final model.
ivmDisplay(model);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
[kern, noise, ivmInfo] = ivmDeconstruct(model);
save(['dem' capName num2str(experimentNo) '.mat'], ...
     'kern', ...
     'noise', ...
     'ivmInfo');

if exist('printDiagram') & printDiagram
  ivmPrintPlot(model, 'ivmContour', [], [], [], capName, experimentNo);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function value = i4_log_10 ( i )

%*****************************************************************************80
%
%% I4_LOG_10 returns the integer part of the logarithm base 10 of ABS(X).
%
%  Example:
%
%        I  VALUE
%    -----  --------
%        0    0
%        1    0
%        2    0
%        9    0
%       10    1
%       11    1
%       99    1
%      100    2
%      101    2
%      999    2
%     1000    3
%     1001    3
%     9999    3
%    10000    4
%
%  Discussion:
%
%    I4_LOG_10 ( I ) + 1 is the number of decimal digits in I.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number whose logarithm base 10 is desired.
%
%    Output, integer VALUE, the integer part of the logarithm base 10 of
%    the absolute value of X.
%
  i = floor ( i );

  if ( i == 0 )

    value = 0;

  else

    value = 0;
    ten_pow = 10;

    i_abs = abs ( i );

    while ( ten_pow <= i_abs )
      value = value + 1;
      ten_pow = ten_pow * 10;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function intVal=monomialIntCubicalShell(alpha,K1,K2)
%%MONOMIALINTCUBALSHELL Evaluate the integral of prod_{i=1}^N x(i)^alpha(i)
%              taken over the region K1<=x(i)<=K2 with 0<K1<K2<Inf. That
%              is, this function gives the value of the desired monomial
%              integral over an N-dimensional hypercubic shell. The region
%              is designated as C_N^{shell} in [1].
%
%INPUTS: alpha An NX1 or 1XN vector of the integer exponents of the
%              monomial term. All elements must be >=0.
%       K1, K2 Two positive, finite, real values such that 0<K1<K2.
%
%OUTPUTS: intVal The value of the specified integral.
%
%The formula is taken from Chapter 7.3 of [1]. These types of explicit
%moment formulae are useful for testing cubature integration points.
%
%EXAMPLE:
%Here we verify that the moment value produced here equals that obtained
%using cubature points of an appropriate order.
% K1=0.5;
% K2=1.25;
% [xi,w]=thirdOrderNDimShellCubPoints(3,K1,K2);
% alpha=[0;2;0];
% theMoment=findMomentFromSamp(alpha,xi,w)
% intVal=monomialIntCubicalShell(alpha,K1,K2)
%One will find that both theMoment and intVal are about 8.0547.
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%
%February 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(all(mod(alpha,2)==0))
    alphaSum=sum(alpha);
    n=length(alpha);

    intVal=(2^n/prod(alpha+1))*(K2^(n+alphaSum)-K1^(n+alphaSum));
else
    intVal=0;
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [edges]=renumberEdges(edges)

% create mapping to points
idx = unique(edges(:));
for i = 1 : length(idx)
    map(idx(i),1) = i;
end

% renumber bnde to correspond with pfix
for i = 1 : size(edges,1)
    edges(i,1) = map(edges(i,1));
    edges(i,2) = map(edges(i,2));
end

end

author: CHLNDDEV


text:
function m = computem (L, w, W, K, Phi)

  % Get the number of documents.
  D = length(L);

  % Initialize the counts "m".
  m = zeros(W,K);

  % Repeat for every document, then for every word in the document.
  is = 0;
  for d = 1:D
    is = is(end) + (1:L(d));
    for i = is
      j      = w(i);
      m(j,:) = m(j,:) + Phi(:,i)';
    end
  end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function motzkin_test ( )

%*****************************************************************************80
%
%% MOTZKIN_TEST tests MOTZKIN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'MOTZKIN_TEST\n' );
  fprintf ( 1, '  MOTZKIN computes the Motzkin numbers A(0:N).\n' );
  fprintf ( 1, '  A(N) counts the paths from (0,0) to (N,0).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I   A(I)\n' );
  fprintf ( 1, '\n' );

  a = motzkin ( n );

  for i = 0 : n
    fprintf ( 1, '  %2d  %10d\n', i, a(i+1) );
  end

  return
end


text:
classdef poisson_naive_bayes_CL

% poisson_naive_bayes_CL is a classifier (CL) object that implements a
%   Poisson Naive Bayes classifier.  For each class, the expected number of occurances
%   (denoted lambda) is calculated separately for each feature/neuron, by taking the mean
%   values from the training data for each class.  To evaluate whether a given test point
%   belongs to class i, the log of the likelihood function is calculated using the lambda values as
%   parameters of Poisson distributions for each feature/neuron (i.e., a separate Poisson distribution
%   is calculated for each neuron), and the probability of observing a given number of spikes for
%   a particular neuron is calculated using the lambda value for that
%   neuron.  The overall likelihood value is calculated by multiplying the
%   probabilities for each neuron together (i.e., Naive Bayes classifiers assume that each
%   neuron is independent), or equivalently, adding the log of the probabilities for
%   each neuron together.  The class with the highest likelihood value is chosen as the
%   predicted label, and the decision values are the log likelihood values.
%
%
% Like all CL objects, there are two main methods, which are:
%
%  1.  cl = train(cl, XTr, YTr)
%       This method takes the training data (XTr, YTr) and learns a mean vector
%       (i.e., the lambda values) for each class.
%
%  2. [predicted_labels decision_values] = test(cl, XTe)
%       This method takes the test data an[d calculates the log of the likelihood
%       function for each class (which are the decision values), and returns the class
%       with the highest likelihood as the predicted_label.
%
%  Notes:
%
%  1. This classifier assumes that all feature values are integers.
%
%  2. If the estimate rate parameter (lambda) for any feature is 0 (i.e., if the mean of a feature in the training data is 0),
%       then this 0 lambda estimate will be replaced by the value 1/(num_training_points_in_class_k +1); i.e., we will assume
%       that there is one more training point in which an event occurred.  This prevents errors if an event occurred on a test point
%       but the estimated lambda was 0.
%
%
%   XTr and XTe are in the form [num_features x num_examples]
%   YTr is in the form [num_examples x 1]
%

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

%     This code is part of the Neural Decoding Toolbox.
%     Copyright (C) 2011 by Ethan Meyers (emeyers@mit.edu)
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

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

  properties
        lambdas = [];   % the expected number of occurances for each neuron for each class
        labels = [];  % the unique labels for the different classes
    end

    methods

        % constructor
        function cl = poisson_naive_bayes_CL
        end

        function cl = train(cl, XTr, YTr)

            % added sanity check
            if size(YTr, 2) ~= size(XTr, 2)  &&  size(YTr, 1) ~= size(XTr, 2)
               error('Number of columns in YTr, and XTr must be the same (i.e., there must be one and exactly one label for each data point)')
            end

            % sanity check to make sure the training data only contains integers
%            if sum(sum(abs(round(XTr) - XTr))) > 10.^-2  % (not having it be exactly zero to account for round off error)
%                error(['The training data must only contain integers to use this classifier.  Make sure the data was loaded as integers (e.g., using the appropriate flag in ' ...
%                    'basic_DS/generalization_DS datasources, and that only appropriate feature preprocessors were used (e.g., do not use the zscore_normalize_FP)'])
%            end

            cl.labels = unique(YTr);

            for iLabel = 1:length(cl.labels )

                lambdas(:, iLabel) = mean(XTr(:, (YTr == cl.labels(iLabel))), 2);

                % If the rate parameter (lambda) is zero, there can be problems if the test data has a value greater than zero
                % (i.e., will get zero probability, which creates problems when taking the log).
                % We will deal with this problem by assuming that there is one more data point with a 1 (when all training data has zero occurances)
                zero_lambda_replacement_value = 1/(length(find(YTr == cl.labels(iLabel))) + 1);   % if found zeros for all trials, assume that there is one more trial where a 1 was found...
                lambdas((lambdas(:, iLabel) == 0), iLabel) = zero_lambda_replacement_value;

            end

            cl.lambdas = lambdas;

        end

        function [predicted_labels decision_values] = test(cl, XTe)

            % sanity check to make sure the test data only contains integers
%            if sum(sum(abs(round(XTe) - XTe))) > 10.^-2   % (not having it be exactly zero to account for round off error)
%                error(['The test data must only contain integers to use this classifier.  Make sure the data was loaded as integers (e.g., using the appropriate flag in ' ...
%                    'basic_DS/generalization_DS datasources, and that only appropriate feature preprocessors were used (e.g., do not use the zscore_normalize_FP)'])
%            end

            % compute simultaneously all p-values for all features and test points (easy b/c everything is independent)
            curr_lambdas = repmat(cl.lambdas, [1, 1, size(XTe, 2)]);
            XTe_repmat_for_all_classes = permute(repmat(XTe, [1, 1, size(cl.lambdas, 2)]), [1 3 2]);
            %p_vals = exp(-curr_lambdas + XTe_repmat_for_all_classes  .* log(curr_lambdas) - gammaln(XTe_repmat_for_all_classes  + 1));
            %log_likelihoods = squeeze(sum(log(p_vals)));

            log_likelihoods = squeeze(sum(-curr_lambdas + XTe_repmat_for_all_classes  .* log(curr_lambdas) - gammaln(XTe_repmat_for_all_classes  + 1), 1)); % possibly even faster...
            % % %log_likelihoods = squeeze(sum(-curr_lambdas + XTe_repmat_for_all_classes  .* log(curr_lambdas))); % since the factorial term does not depend on the class, code might be faster
                                                                                % if this term is removed.  The decision values returned
                                                                                % will no longer be exactly the log likelihood values but rather the log likelihood
                                                                                % values minus a constant.  Using this will make AUROC values worse, so don't do it!!!!

           % prior to version 1.0.4 the following line was used.  I added a the sum over the first dimention (i.e., sum(XXX, 1)) so that the code will run even when only a single site is used
           %log_likelihoods = squeeze(sum(-curr_lambdas + XTe_repmat_for_all_classes  .* log(curr_lambdas) - gammaln(XTe_repmat_for_all_classes  + 1))); % possibly even faster...

            [vals inds] = randmax(log_likelihoods);
            predicted_labels  = cl.labels(inds);
            decision_values = log_likelihoods';

        end

    end  % end public methods

end   % end classdef

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [FV2]=refinepatch(FV)
% This function "refinepatch" refines a triangular mesh with
% a spline interpolated 4-split method.
%
%   [FV2] = refinepatch(FV,options)
%
% inputs,
%   FV : Structure containing a Patch, with
%        FV.vertices the mesh vertices
%        FV.face the mesh faces (triangles), rows with each 3 vertex indices
% outputs,
%   FV2 : Structure Containing the refined patch
%
%
% Reference:
%  The spline interpolation of the face edges is done by the
%  Opposite Edge Method, described in: "Construction of Smooth Curves
%  and Surfaces from Polyhedral Models" by Leon A. Shirman
%
% How it works:
%  The tangents (normals) and velocity on the edge points of all edges
%  are calculated. Which are  later used for b-spline interpolation when
%  splitting the edges.
%
%  A tangent on an 3D line or edge is under defined and can rotate along
%  the line, thus an (virtual) opposite vertex is used to fix the tangent and
%  make it more like a surface normal.
%
%  B-spline interpolate a half way vertices between all existing vertices
%  using the velocity and tangent from the edgepoints. After splitting a
%  new facelist is constructed
%
% Speed:
%  Compile the c-functions for more speed with:
%   mex vertex_neighbours_double.c -v;
%   mex edge_tangents_double.c -v;
%
% Example:
%
% X=[-0.5000;  0.5000;  0.0000;  0.0000];
% Y=[-0.2887; -0.2887;  0.5774;  0.0000];
% Z=[ 0.0000;  0.0000;  0.0000;  0.8165];
% FV.vertices=[X Y Z];
%
% FV.faces=[2 3 4; 4 3 1; 1 2 4; 3 2 1];
%
% figure, set(gcf, 'Renderer', 'opengl'); axis equal;
% for i=1:4
%   patch(FV,'facecolor',[1 0 0]);
%   pause(2);
%   [FV]=refinepatch(FV);
% end
%
% Function is written by D.Kroon University of Twente (February 2010)

% Get the neighbour vertices of each vertice from the face list.
Ne=vertex_neighbours(FV);

% Calculate the tangents (normals) and velocity of all edges. Which is
% later used for b-spline interpolation and split of the edges
%
% A tangent on an 3D line or edge is under defined and can rotate along
% the line, thus an (virtual) opposite vertex is used to fix the tangent and
% make it more like a surface normal.
V=FV.vertices; F=FV.faces;
[ET_table,EV_table,ETV_index]=edge_tangents(V,Ne);

% B-spline interpolate a half way vertices between all existing vertices
% using the velocity and tangent from above
[V,HT_index, HT_values]=make_halfway_vertices(EV_table,ET_table,ETV_index,V,Ne);

% Make new facelist
Fnew=makenewfacelist(F,HT_index,HT_values);

FV2.vertices=V;
FV2.faces=Fnew;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function C = cat(dim, varargin)
%CAT   Concatenation of CHEBMATRIX objects.
%   C = CAT(dim, A, B,...) concatenates chebmatrices along the indicated
%   dimension following Matlab conventions.
%
% See also CHEBMATRIX.HORZCAT, CHEBMATRIX.VERTCAT.

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

% Any singleton operator blocks must be encased in cells to match up with
% CHEBMATRIX blocks.

if ( dim == 1 )
    C = vertcat(varargin{:});
elseif ( dim == 2 )
    C = horzcat(varargin{:});
else
    error('CHEBFUN:CHEBMATRIX:cat:dim', ...
        'Must concatenate along dimension 1 or 2.')
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Example
% This script is used to initialize the hand guiding functionality on KUKA
% robot remotely

% Mohammad SAFEEA, 9th of June 2017

% TO use this example:
% 1- Start the server on KUKA
% 2- Start the client on Matlab
% 3- Press the white button and hand guide the robot
% 4- To save the coordinates of the robot, long click on the green button
% 5- Repeat three times
% 6- The robot gives you an interval of five seconds, clear the area around
% the robot directly, the robot will reproduce the motion tought.

HandTeachingExample();

repo: KST-Kuka-Sunrise-Toolbox


text:
function [ status ] = export_poscar( filename, geometry )
%EXPORT_POSCAR Export a geometry struct as a VASP POSCAR file.
%   status = EXPORT_POSCAR(filename,geometry) exports the geometry as a
%   VASP POSCAR file.
%
%   See also IMPORT_POSCAR.

    % write selective dynamics and chemical symbols

    fid = fopen(filename,'w');
    if fid==-1
        error(['Error opening ' filename]);
    end

    fprintf(fid,[geometry.comment '\n']);
    fprintf(fid,'1.0\n'); % scale factor
    fprintf(fid, '%19.16f %19.16f %19.16f\n', geometry.lattice'); % lattice vectors
    if ~isempty(geometry.symbols)
        cellfun(@(x) fprintf(fid, '%s ', x), geometry.symbols);
        fprintf(fid, '\n');
    end
    fprintf(fid, '%d ', geometry.atomcount); % number of each species
    fprintf(fid, '\nDirect\n');
    fprintf(fid, '%19.16f %19.16f %19.16f\n', geometry.coords');

    fclose(fid);

    status = 0;

end

author: Sable

repo: mcbench-benchmarks


text:
function [p,p_lp] = addSymmetryCuts(p,p_lp)

for j = 1:length(p.sdpsymmetry)
    if length(p.sdpsymmetry{j}.variables{1}) <= 4
        % We can enumerate easily infeasibles
        excludes = [];
        n = sqrt(size(p.sdpsymmetry{j}.dataBlock,1));
        vars = p.sdpsymmetry{j}.variables{1};
        % Now enumerate, for historical reasons
        % also accept {-1,0}-models so translate and scale
        L = p.lb(vars);
        U = p.ub(vars);
        m = length(vars);
        combs = dec2decbin(0:2^m-1,m)';
        combs = repmat(L,1,size(combs,2)) + repmat(U-L,1,size(combs,2)).*combs;
        combs = unique(combs','rows')';
        feasible = ones(1,size(combs,2));
        for i = 1:size(combs,2)
            if min(eig(reshape(p.sdpsymmetry{j}.dataBlock*[1;combs(:,i)],n,n))) < -abs(p_lp.options.bnb.feastol)
                excludes = [excludes combs(:,i)];
                feasible(i) = 0;
            end
        end
        % We can model the inconsistencies in various ways
        % 1. Binary variables sum less
        % 2. Binary variables sum larger
        % 3. Generic, just exclude
        if ~isempty(excludes)
            done = 0;
            sum_feas = sum(combs(:,find(feasible)),1);
            sum_infeas = sum(excludes,1);
            newF = [];
            if all(sum_infeas==sum_infeas(1)) && all(sum_feas < sum_infeas(1) | sum_feas > sum_infeas(1))
                % There is a forbidden sum among these
                for s = 1:length(p.sdpsymmetry{j}.variables)
                    p.forbiddencardinality.variables{end+1} = p.sdpsymmetry{j}.variables{s};
                    p.forbiddencardinality.value{end+1} = sum_infeas(1);
                end
            end
            for k = min(sum(combs,1)):max(sum(combs,1))
                if all(sum_feas<=k) && all(sum_infeas>k)
                    for s = 1:length(p.sdpsymmetry{j}.variables)
                        a = spalloc(1,length(p_lp.c),1);
                        a(p.sdpsymmetry{j}.variables{s}) = -1;
                        newF = [newF;k a];
                    end
                    done = 1;
                    break
                elseif  all(sum_feas>=k) && all(sum_feas<k)
                    for s = 1:length(p.sdpsymmetry{j}.variables)
                        a = spalloc(1,length(p_lp.c),1);
                        a(p.sdpsymmetry{j}.variables{s}) = 1;
                        newF = [newF;-k a];
                    end
                    done = 1;
                    break
                end
            end
            if ~done
                % No good structure, just add no-goods
                infeasible_combinations = unique(excludes','rows')';
                for k = 1:size(infeasible_combinations,2)
                    % Local cut for reduced set
                    if min(min(excludes))==-1
                        binary_type = -1;
                    else
                        binary_type = 1;
                    end
                    [b,atemp] = exclusionCut(infeasible_combinations(:,k),binary_type);
                    % Add that cut for every variable groups
                    for s = 1:length(p.sdpsymmetry{j}.variables)
                        a = spalloc(1,length(p_lp.c),1);
                        a(p.sdpsymmetry{j}.variables{s}) = atemp;
                        newF = [newF;b a];
                    end
                end
            end
            p_lp = addInequality(p_lp,newF);
            % Delete, won't need these in the future
            p.sdpsymmetry{j}.dataBlock = [];
            p.sdpsymmetry{j}.variables = [];
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%

Sigma = (X'*X) ./ m;
[U, S, V] = svd(Sigma);

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

end

author: SaveTheRbtz

repo: ml-class


text:
function [ schema_struct ] = parse_sicd_schema( filename )
% PARSE_SICD_SCHEMA Parse SICD/CPHD schema XSD into a MATLAB structure
% It is MUCH faster to traverse through a MATLAB structure than DOM nodes,
% so we want to convert the schema info into a MATLAB structure before
% traversing through the SICD/CPHD XML.
%
% Author: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Open file
schema_root = xmlread(filename);
schema_root = schema_root.getDocumentElement;
schema_struct.ns = char(schema_root.getAttribute('xmlns'));

% Separate type definitions and master structure definition
for i=1:schema_root.getLength
    child=schema_root.item(i-1);
    if child.getNodeType == child.ELEMENT_NODE
        switch char(child.getNodeName)
            case {'xs:simpleType','xs:complexType','xs:group'} % Type definitions
                schema_struct.types.(char(child.getAttribute('name'))) = recursfun_schema(child);
            case 'xs:element' % Master node (should be onle one)
                schema_struct.master = recursfun_schema(child);
            otherwise
                % Should not be in other types of nodes in top level
                error('SICDXML2STRUCT:UNEXPECTED_NODE_TYPE','Unexpected node type in top level of XSD.');
        end
    end
end

    function output_struct = recursfun_schema(current_node)
        output_struct = struct();
        for j=1:current_node.getLength
            current_child=current_node.item(j-1);
            if current_child.getNodeType == current_child.ELEMENT_NODE
                switch char(current_child.getNodeName)
                    case 'xs:element'
                        name_str = char(current_child.getAttribute('name'));
                        type_str = char(current_child.getAttribute('type'));
                        if ~isempty(type_str)
                            output_struct.(name_str).SCHEMA_type = type_str;
                        else % Element with empty type.  Should have a structure defined within it.
                            output_struct.(name_str) = recursfun_schema(current_child);
                        end
                    case {'xs:restriction','xs:extension'}
                        output_struct = recursfun_schema(current_child); % Adds any attributes
                        output_struct.SCHEMA_type = char(current_child.getAttribute('base'));
                    case {'xs:simpleType','xs:simpleContent',...
                            'xs:complexType','xs:complexContent','xs:group'}
                        output_struct = recursfun_schema(current_child);
                    case {'xs:sequence','xs:choice','xs:all'}
                        output_struct = setstructfields(output_struct, recursfun_schema(current_child));
                    case {'xs:attribute'}
                        if isfield(output_struct,'SCHEMA_attributes')
                            output_struct.SCHEMA_attributes = ...
                                [output_struct.SCHEMA_attributes, ...
                                char(current_child.getAttribute('name'))];
                        else
                            output_struct.SCHEMA_attributes = ...
                                {char(current_child.getAttribute('name'))};
                        end
                    case {'xs:minInclusive','xs:maxInclusive',...
                            'xs:minExclusive','xs:maxExclusive',...
                            'xs:enumeration','xs:pattern'}
                        % These fields are expected, but we don't use them
                        % for anything.
                    otherwise
                        error('SICDXML2STRUCT:UNRECOGNIZED_NODE_TYPE','Unrecognized node type in XSD.');
                end
            end
        end
    end
end

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

author: ngageoint


text:
function newDepth = SiftFu(sequenceName, frameIDs)

%{
Please cite the following paper if you use this code

Citation:
J. Xiao, A. Owens and A. Torralba
SUN3D Database: Semantic RGB-D Bundle Adjustment with Human in the Loop
Proceedings of 14th IEEE International Conference on Computer Vision (ICCV2013)
%}

addpath(genpath('SIFTransac'))
vl_setup;

global toVisualize;
toVisualize = true;

%% IO
if ~exist('sequenceName','var')
    % load demo sequence
    % look for all sequence name list at http://sun3d.csail.mit.edu/player/list.html
    %sequenceName = 'hotel_mr/scan1';
    %sequenceName = 'hotel_umd/maryland_hotel3';
    %sequenceName = 'brown_bm_1/brown_bm_1';
    sequenceName = 'mit_32_d428/bs4j179mmv';
end

% the root path of SUN3D
% change it to local if you downloaded the data
%SUN3Dpath = '/data/vision/torralba/sun3d/record/scene_final';
SUN3Dpath = 'http://sun3d.csail.mit.edu/data/';

% read intrinsic
global K;
K = reshape(readValuesFromTxt(fullfile(SUN3Dpath,sequenceName,'intrinsics.txt')),3,3)';

% file list
imageFiles = dirSmart(fullfile(SUN3Dpath,sequenceName,'image/'),'jpg');
depthFiles = dirSmart(fullfile(SUN3Dpath,sequenceName,'depth/'),'png');

% read time stamp
imageFrameID = zeros(1,length(imageFiles));
imageTimestamp = zeros(1,length(imageFiles));
for i=1:length(imageFiles)
    id_time = sscanf(imageFiles(i).name, '%d-%d.jpg');
    imageFrameID(i) = id_time(1);
    imageTimestamp(i) = id_time(2);
end
depthFrameID = zeros(1,length(depthFiles));
depthTimestamp = zeros(1,length(depthFiles));
for i=1:length(depthFiles)
    id_time = sscanf(depthFiles(i).name, '%d-%d.png');
    depthFrameID(i) = id_time(1);
    depthTimestamp(i) = id_time(2);
end

% synchronize: find a depth for each image
frameCount = length(imageFiles);
IDimage2depth = zeros(1,frameCount);
for i=1:frameCount
    [~, IDimage2depth(i)]=min(abs(double(depthTimestamp)-double(imageTimestamp(i))));
end

%plot(double(imageTimestamp)-double(depthTimestamp(IDimage2depth)))

if ~exist('frameIDs','var')
    frameIDs = 1:frameCount;
else
    frameIDs = frameIDs(frameIDs>=1 & frameIDs<=frameCount);
end

imageFiles = imageFiles(frameIDs);
depthFiles = depthFiles(IDimage2depth(frameIDs));

%% kinect fusion

% for optimization

global VMap;
global NMap;
global CMap;

global XYZcam;
global Ncam;
global IMGcam;

%% tsdf

global voxel;
global tsdf_value;
global tsdf_weight;
global tsdf_color;    tsdf_color = [];

voxel.unit = 0.01; % Kevin: 4mm = 0.004 meter. Kinect cannot go better than 3mm
voxel.mu_grid = 10; % used to be 4
voxel.size_grid = [512; 512; 1024]; % [512; 512; 512];

voxel.range(1,1) = - voxel.size_grid(1) * voxel.unit / 2;
voxel.range(1,2) = voxel.range(1,1) + (voxel.size_grid(1)-1) * voxel.unit;

voxel.range(2,1) = - voxel.size_grid(2) * voxel.unit / 2;
voxel.range(2,2) = voxel.range(2,1) + (voxel.size_grid(2)-1) * voxel.unit;

voxel.range(3,1) = -0.5; % - voxel.size_grid(3) * voxel.unit / 2;
voxel.range(3,2) = voxel.range(3,1) + (voxel.size_grid(3)-1) * voxel.unit;

voxel.mu = voxel.mu_grid * voxel.unit;

fprintf('memory = %f GB\n',  prod(voxel.size_grid) * 4 / (1024*1024*1024));
fprintf('space = %.2f m x %.2f m x %.2f m ', voxel.size_grid(1) * voxel.unit, voxel.size_grid(2) * voxel.unit, voxel.size_grid(3) * voxel.unit);
fprintf('= [%.2f,%.2f] x [%.2f,%.2f] x [%.2f,%.2f]\n',voxel.range(1,1),voxel.range(1,2),voxel.range(2,1),voxel.range(2,2),voxel.range(3,1),voxel.range(3,2));

%tsdf_value  = -ones([voxel.size_grid(1),voxel.size_grid(2), voxel.size_grid(3)],'single');
tsdf_value  =  ones([voxel.size_grid(1),voxel.size_grid(2), voxel.size_grid(3)],'single');
tsdf_weight = zeros([voxel.size_grid(1),voxel.size_grid(2), voxel.size_grid(3)],'single');

% comment out this line to avoid using color
% tsdf_color  = zeros([voxel.size_grid(1),voxel.size_grid(2), voxel.size_grid(3)], 'uint32');

f = K(1,1);

global ViewFrustumC;
ViewFrustumC = [...
    0 -320 -320  320  320;
    0 -240  240  240 -240;
    0    f    f    f    f];
ViewFrustumC = ViewFrustumC/f * 8; % 8 meter is the furthest depth of kinect

% precompute
global raycastingDirectionC;
[pX,pY]=meshgrid(1:640,1:480);

raycastingDirectionC = [pX(:)'-K(1,3); pY(:)'-K(2,3); f*ones(1,640*480)]; % clipping at 8 meter is the furthest depth of kinect
raycastingDirectionC = raycastingDirectionC ./ repmat(sqrt(sum(raycastingDirectionC.^2,1)),3,1);

%%
useSIFT = false;
if useSIFT
    MatchPairs = cell(1,length(imageFiles)-1);
end

cameraRtC2W = repmat([eye(3) zeros(3,1)], [1,1,length(imageFiles)]);

dispclim = [0 5];

for frameID = 1:length(imageFiles)

    fprintf('================================ Frame %d ================================\n',frameID);

    IMGcam = imageRead(fullfile(fullfile(SUN3Dpath,sequenceName,'image',imageFiles(frameID).name)));

    %subplot(3,4,4)
    %imagesc(IMGcam);  axis equal; axis tight;
    %drawnow;
    %title('input color');

    if frameID==1
        IMGcam1 = IMGcam;
    end

    % read the frame

    depth = depthRead(fullfile(fullfile(SUN3Dpath,sequenceName,'depth',depthFiles(frameID).name)));
    XYZcam = depth2XYZcamera(K, depth);

    if frameID==1
        XYZcam1 = XYZcam;
    end

    Ncam = vertex2normal(XYZcam);

    if toVisualize
        if frameID==1
            subplot(3,4,9)
        else
            subplot(3,4,1)
        end
        imagesc(XYZcam(:,:,3),dispclim); axis equal; axis tight
        title(sprintf('Frame %d: Input Depth',frameID));
        axis off

        if frameID==1
            subplot(3,4,10)
        else
            subplot(3,4,2)
        end
        imagesc((Ncam+1)/2); axis equal; axis tight
        title(sprintf('Frame %d: Input Normal',frameID));
        axis off

        if frameID==1
            subplot(3,4,11)
        else
            subplot(3,4,3)
        end
        raycastingDirectionW = transformRTdir(raycastingDirectionC,[eye(3) zeros(3,1)]);
        DotMap = reshape(max(0,sum(-reshape(Ncam,[480*640 3])' .* raycastingDirectionW,1)),[480 640]);
        imagesc(DotMap);
        colormap('gray'); axis equal; axis tight
        title(sprintf('Frame %d: Input Phong',frameID));
        axis off
    end

    if frameID==1
        camRtC2W = [eye(3) [0;0;0]];
    else
        MatchPairs{frameID-1} =  align2view(1,IMGcam1,XYZcam1,frameID,IMGcam,XYZcam);
        camRtC2W = MatchPairs{frameID-1}.Rt;

        if size(MatchPairs{frameID-1}.matches,2)<5
            disp('SIFT matching failed, ignoring this frame');
            continue;
        end
    end

    cameraRtC2W(:,:,frameID) = camRtC2W;

    %% update TSDF
    disp('update TSDF...');
    %tic
    updateTSDF(camRtC2W);
    %toc

    %{
    %% visualizing the voxel
    subplot(3,4,10)
    for i=1:voxel.size_grid(3)
        if min(min(min(tsdf_value(:,:,i)))) ~= max(max(max(tsdf_value(:,:,i))))
            imagesc(tsdf_weight(:,:,i)'); axis equal; axis tight; xlabel('x'); ylabel('y');
            title(['frame ' num2str(i) ' = depth ' num2str((i-1)*voxel.unit+voxel.range(3,1)) ' meter']);
            pause(0.05);
        end
    end

    subplot(3,4,10)
    for i=1:voxel.size_grid(3)
        if min(min(min(tsdf_value(:,:,i)))) ~= max(max(max(tsdf_value(:,:,i))))
            imagesc(tsdf_value(:,:,i)'); axis equal; axis tight; xlabel('x'); ylabel('y');
            title(['frame ' num2str(i) ' = depth ' num2str((i-1)*voxel.unit+voxel.range(3,1)) ' meter']);
            pause(0.05);
        end
    end
    %}

    %% ray casting for the result
    disp('ray casting');
    tic
    % ray casting result is in the world coordinate
    [VMap,NMap,tMap,CMap] = raycastingTSDFvectorized([eye(3) [0;0;0]], [0.4 8]);
    %[VMap,NMap,tMap] = raycastingTSDFdump(camRtC2WrayCasting, [0.4 8]);
    toc

    newDepth = reshape(VMap(3,:),480,640);

    if toVisualize
        subplot(3,4,5)
        imagesc(newDepth,dispclim); axis equal; axis tight; axis off
        title(sprintf('Fused Depth',frameID));
    end

    VMap = reshape(double(VMap'),480,640,3);
    VMap(:,:,4) = ~isnan(VMap(:,:,1));
    NMap = vertex2normal(VMap);
    % normalize normal map
    %NMap = NMap ./ repmat(sqrt(sum(NMap.^2,1)),3,1);
    if toVisualize
        subplot(3,4,6)
        imagesc((NMap+1)/2); axis equal; axis tight; axis off
        title(sprintf('Fused Noraml',frameID));
    end

    if toVisualize
        subplot(3,4,7)
        raycastingDirectionW = transformRTdir(raycastingDirectionC,[eye(3) zeros(3,1)]);
        DotMap = reshape(max(0,sum(-reshape(NMap,[480*640 3])' .* raycastingDirectionW,1)),[480 640]);
        imagesc(DotMap);
        colormap('gray'); axis equal; axis tight; axis off
        title(sprintf('Fused Phong',frameID));
    end

    if toVisualize
        subplot(3,4,12)
        imagesc(abs(newDepth - XYZcam1(:,:,3)),dispclim); axis equal; axis tight; axis off
        title(sprintf('Difference of Depth',frameID));

        drawnow;
    end

    %{

    subplot(3,4,6)
    imagesc(min(1,max(0,(reshape(NMap',[480 640 3])+1)/2))); axis equal; axis tight
    title(sprintf('Frame %d: rayCasting Normal Map',frameID));
    drawnow;

    %figure
    %visNMap = reshape(NMap',[480 640 3]);
    %out = visualizeNormals(-visNMap);
    %imagesc(out)
    %title('rayCasting Normal Map');

    subplot(3,4,7)
    raycastingDirectionW = transformRTdir(raycastingDirectionC,camRtC2W);
    DotMap = reshape(max(0,sum(-NMap .* raycastingDirectionW,1)),[480 640]);
    imagesc(DotMap);
    colormap('gray'); axis equal; axis tight
    title(sprintf('Frame %d: phong shading',frameID));
    drawnow;

    if ~isempty(tsdf_color)
        subplot(3,4,8)
        imagesc(reshape(CMap',[480,640,3])/255);  axis equal; axis tight;
        drawnow;
        title(sprintf('Frame %d: ray casting color',frameID));
    end

    %}

    %subplot(3,4,9)
    %imagesc(reshape(tMap,480,640)); axis equal; axis tight
    %title(sprintf('Frame %d: rayCasting tMap',frameID));
    %drawnow;

end

end

%% IO functions

function values = readValuesFromTxt(filename)
try
    values = textscan(urlread(filename),'%f');
catch
    fid = fopen(filename,'r');
    values = textscan(fid,'%f');
    fclose(fid);
end
values = values{1};
end

function XYZcamera = depth2XYZcamera(K, depth)
[x,y] = meshgrid(1:640, 1:480);
XYZcamera(:,:,1) = (x-K(1,3)).*depth/K(1,1);
XYZcamera(:,:,2) = (y-K(2,3)).*depth/K(2,2);
XYZcamera(:,:,3) = depth;
XYZcamera(:,:,4) = depth~=0;
end

function depth = depthRead(filename)
depth = imread(filename);
depth = bitor(bitshift(depth,-3), bitshift(depth,16-3));
depth = single(depth)/1000;
end
%{
    % test to make sure it is correct
    for i=0:65535
        depth = uint16(i);
        code =bitor(bitshift(depth,3),bitshift(depth,3-16));
        recoverDepth = bitor(bitshift(code,-3), bitshift(code,16-3));
        if (depth~=recoverDepth)
            fprintf('error + %d\n',i);
        end
    end
%}

function image = imageRead(filename)
image = imread(filename);
end

function files = dirSmart(page, tag)
[files, status] = urldir(page, tag);
if status == 0
    files = dir(fullfile(page, ['*.' tag]));
end
end

function [files, status] = urldir(page, tag)
if nargin == 1
    tag = '/';
else
    tag = lower(tag);
    if strcmp(tag, 'dir')
        tag = '/';
    end
    if strcmp(tag, 'img')
        tag = 'jpg';
    end
end
nl = length(tag);
nfiles = 0;
files = [];

% Read page
page = strrep(page, '\', '/');
[webpage, status] = urlread(page);

if status
    % Parse page
    j1 = findstr(lower(webpage), '<a href="');
    j2 = findstr(lower(webpage), '');
    Nelements = length(j1);
    if Nelements>0
        for f = 1:Nelements
            % get HREF element
            chain = webpage(j1(f):j2(f));
            jc = findstr(lower(chain), '">');
            chain = deblank(chain(10:jc(1)-1));

            % check if it is the right type
            if length(chain)>length(tag)-1
                if strcmp(chain(end-nl+1:end), tag)
                    nfiles = nfiles+1;
                    chain = strrep(chain, '%20', ' '); % replace space character
                    files(nfiles).name = chain;
                    files(nfiles).bytes = 1;
                end
            end
        end
    end
end
end

author: jianxiongxiao

repo: ProfXkit


text:
function zInterp=biquadraticInterp(xyPts,x0,y0,deltaX,deltaY,zGrid)
%%BIQUADRATICINTERP Biquadratic interpolation is interpolation of a
%       function that depends on two independent variables, by first
%       performing quadratic interpolation across one dimension, and then
%       performing quadratic interpolation across the next dimension. This
%       function assumes that in x and y generating the values in zGrid is
%       uniform.
%
%INPUTS: xyPts A 2XnumPts set of [x;y] pairs at which the value of z should
%              be interpolated.
%       x0, y0 These are values of the independent variables corresponding
%              to the entry in zGrid(1,1).
% deltaX, deltaY The uniform spacings between the independent variables for
%              gridded elements in zGrid. Thus, zGrid(i,j) holds the z
%              value associated with x=x0+deltaX*(i-1) and
%              y=y0+deltaY*(i-1).
%        zGrid A numXXnumY matrix of values of the independent variable
%              for different values of the dependent variables. The first
%              index selects the x value and the second index the y value.
%              numX>=3 and numY>=3.
%
%OUTPUTS: zInterp A 1XnumPts vector of interpolated values.
%
%Biquadratic interpolation is described in [1]. However, this function does
%not use the implementation that is given at the end of [1].
%
%EXAMPLE:
%A 2D function is evaluated and plotted. Then, the function is evaluated at
%a much smaller number of points and those are used to interpolate and plot
%the function again. The basis points are ploted in black. The
%small discontinuities between interpolation regions can be seen.
% numXPts=400;
% numYPts=401;
% plotXPts=linspace(-1,1,numXPts);
% plotYPts=linspace(-1,1,numYPts);
% z=@(xy)sum(sin(1-(16/15)*xy).^2-(1/50)*sin(4-(64/15)*xy)-sin(1-(16/15)*xy),1);
% [Y,X]=meshgrid(plotYPts,plotXPts);
% xy=[X(:).';Y(:).'];
% Z=reshape(z(xy),[numXPts,numYPts]);
%
% figure(1)
% clf
% hold on
% surface(X,Y,Z,'edgeColor','None')
%
% %Get a small number of points to use for interpolation and mark the
% %points on the plot.
% numXBasis=10;
% numYBasis=11;
% x0=-1;
% y0=-1;
% xBasis=linspace(x0,1,numXBasis);
% yBasis=linspace(y0,1,numYBasis);
% deltaX=xBasis(2)-xBasis(1);
% deltaY=yBasis(2)-yBasis(1);
% [YB,XB]=meshgrid(yBasis,xBasis);
% xyB=[XB(:).';YB(:).'];
% ZGrid=reshape(z(xyB),[numXBasis,numYBasis]);
% scatter3(XB(:),YB(:),ZGrid(:),200,'.k')
% view(25,50)
%
% %Now, interpolate the surface using the above basis points.
% ZInterp=biquadraticInterp(xy,x0,y0,deltaX,deltaY,ZGrid);
% ZInterp=reshape(ZInterp,[numXPts,numYPts]);
%
% %Plot the biquadratically interpolated surface and the control points.
% figure(2)
% clf
% hold on
% surface(X,Y,ZInterp,'edgeColor','None')
% scatter3(XB(:),YB(:),ZGrid(:),200,'.k')
% view(25,50)
%
%REFERENCES:
%[1] D. Smith, "NOAA technical memorandum NOS NGS 84: Biquadratic
%    interpolation," National Oceanic and Atmospheric Administration,
%    National Geodetic Survey, Silver Spring, MD, Tech. Rep., Sep. 2020.
%
%July 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numPts=size(xyPts,2);
zInterp=zeros(1,numPts);

numX=size(zGrid,1);
numY=size(zGrid,2);

for curPt=1:numPts
    xCur=xyPts(1,curPt);
    yCur=xyPts(2,curPt);

    idxX=round((xCur-x0)/deltaX)+1;
    idxY=round((yCur-y0)/deltaY)+1;

    %The points cannot be endpoints.
    idxX=idxX+(idxX==1)-(idxX==numX);
    idxY=idxY+(idxY==1)-(idxY==numY);

    %First, interpolate across X for each Y.
    x0Cur=x0+deltaX*(idxX-2);
    y0Cur=y0+deltaY*(idxY-2);

    zInterp1=quadraticInterpUnif3Pt(xCur,x0Cur,deltaX,zGrid(idxX-1,idxY-1),zGrid(idxX,idxY-1),zGrid(idxX+1,idxY-1));
    zInterp2=quadraticInterpUnif3Pt(xCur,x0Cur,deltaX,zGrid(idxX-1,idxY),zGrid(idxX,idxY),zGrid(idxX+1,idxY));
    zInterp3=quadraticInterpUnif3Pt(xCur,x0Cur,deltaX,zGrid(idxX-1,idxY+1),zGrid(idxX,idxY+1),zGrid(idxX+1,idxY+1));

    zInterp(curPt)=quadraticInterpUnif3Pt(yCur,y0Cur,deltaY,zInterp1,zInterp2,zInterp3);
end
end

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


text:
function img = imresizeAspect(img, maxSize)
% function img = imresizeAspect(img, maxSize)
% If image is larger than max size, reduce size, preserving aspect ratio of input.
%
% If size(input) = [y x] and maxSize = [yy xx],
% then size(output) is given by the following (where a=y/x)
%  if y<yy & x<xx  then [y x]
%  else if a<1 then [a*xx xx]
%  else [yy yy/a]

% For why we use bilinear,
% See http://nickyguides.digital-digest.com/bilinear-vs-bicubic.htm

if isempty(maxSize), return; end

[y x c] = size(img);
a= y/x;
yy = maxSize(1); xx = maxSize(2);
if y <= yy & x <= xx
  % no-op
else
  if a < 1
    img = imresize(img, ceil([a*xx xx]), 'bilinear');
  else
    img = imresize(img, ceil([yy yy/a]), 'bilinear');
  end
  fprintf('resizing from %dx%d to %dx%d\n', y, x, size(img,1), size(img,2));
end

%test
if 0
  maxSize = [240 320];
  %img = imread('C:\Images\Wearables\web_static_office\8.jpg');
  %img = imread('C:\Images\Wearables\web_static_office\billandkirsten.jpg');
  %img = imread('C:\Images\Wearables\web_static_street_april\p5.jpg');
  img = imread('C:\Images\Wearables\Database_static_street\art11.jpg');
  img2 = imresizeAspect(img, maxSize);
  figure(1); clf; imshow(img)
  figure(2); clf; imshow(img2)
  fprintf('%dx%d (%5.3f) to %dx%d (%5.3f)\n', ...
	  size(img,1), size(img,2), size(img,1)/size(img,2), ...
	  size(img2,1), size(img2,2), size(img2,1)/size(img2,2));
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [slope, intercept, averslope, averIC] = fdsurfft(im)
% FDSURFFT Compute fractal dimension (slope) of surface image im and draw rose plots
%          of slope and intercept
%     im: input array of surface image (grayvalue or range image)
%     slope: an array of size 24 which stores the average slopes in 24
%            directions
%     intercept: an array of size 24 which stores the average intercepts in 24
%                directions
%     averslope: average slope for all directions
%     averIC: average intercept for all directions

% This is a matlab version of John C. Russ 's program

% Written by Mr. Jianbo Zhang, J.Zhang@ewi.utwente.nl
% 2 Nov, 2004

tic
NUM_DIR = 24; % number of directions that the frequency space is uniformally divided
NUM_RAD = 30; % number of points that the radius is uniformally divided
if nargin < 1,
    error('Missed input argument which must be an array!')
end

[M N] = size(im);
xctr = 1 + bitshift(N, -1); % x coordinate of center point
yctr = 1 + bitshift(M, -1); % y coordinate of center point
imMean = mean(im(:));
fim = fftshift(fft2(double(im) - imMean));

% power spectrum
mag = log(fim .* conj(fim)+ 10 ^ (-6));
sumBrite = zeros(NUM_DIR, NUM_RAD); %accumulation magnitude for each direction and radius
nCount = zeros(NUM_DIR, NUM_RAD); %number of magnitude
radius = zeros(2 * NUM_RAD,1);%accumulation magnitude for all directions
radCount = zeros(2 * NUM_RAD,1);% number of magintude for all directions

%Compute phase image and phase histogram
phaseim = zeros(M,N);
phase = zeros(180);
for j = 1:M
    for i = 1:N
        realv = real(fim(j,i));
        imagv = imag(fim(j,i));
        if realv == 0
            value = pi/2;
        else
            value = atan((imagv / realv));
            phaseim(j, i) = value;
            ang = floor(180 * (pi / 2 + value) / pi);
        end
        if ang < 0
            ang = 0;
        end
        if ang > 179
            ang = 179;
        end
        phase(ang + 1) = phase(ang + 1) + 1;
    end
end
maxphase = max(phase);
figure; imshow(phaseim, []);
title('Phase image');
figure; plot(phase / maxphase);
title('Phase histogram (0...2 \pi)');%
axis off

%accumulation of magnitude for each dirction and radius
rmax = log(min(M,N)/2);% maximum radius
for j = 1:M
    if j ~= yctr
        yval = yctr - j;
        y2 = yval * yval;
        for i = 1:N
            if i ~= xctr
                xval = i - xctr;
                rho = log(sqrt(y2 + xval * xval));
                if rho > 0 & rho <= rmax
                    mval = mag(j,i);
                    temp = yval /xval;
                    theta = atan(temp);
                    if xval < 0
                        theta = theta + pi;
                    end
                    if theta < 0
                        theta = theta + 2 * pi;
                    end

                    ang = floor(NUM_DIR * theta /(2* pi));
                    if ang > NUM_DIR - 1 | ang < 0
                        ang = NUM_DIR - 1;
                    end
                    k = floor(2*NUM_RAD * rho / rmax);
                    h = floor(k/2);
                    if k > 2 * NUM_RAD - 1
                        h = NUM_RAD - 1;
                        k = 2 * NUM_RAD - 1;
                    end

                    if h >= 5
                        sumBrite(ang + 1, h + 1) = sumBrite(ang + 1, h + 1) + mval;
                        nCount(ang + 1, h + 1) = nCount(ang + 1, h + 1) + 1;
                    end
                    if k >= 5
                        radius(k + 1) = radius(k + 1) + mval;
                        radCount(k + 1) = radCount(k + 1) + 1;
                    end
                end
            end
        end
    end
end

%linear regression
for ang = 1:NUM_DIR
    sumx = 0;
    sumy = 0;
    sumx2 = 0;
    sumxy = 0;
    sumn = 0;
    for range = 6:NUM_RAD
        if nCount(ang, range) > 0
            yval = sumBrite(ang, range)/nCount(ang, range);
            xval = (range -1) * rmax / NUM_RAD;
            sumx = sumx + xval;
            sumy = sumy + yval;
            sumx2 = sumx2 + xval * xval;
            sumxy = sumxy + xval * yval;
            sumn = sumn + 1;
        end
    end
    slope(ang) = (sumn * sumxy - sumx * sumy) / (sumn * sumx2 - sumx * sumx);
    intercept(ang) = (sumy - slope(ang) * sumx) / sumn;
end

%compute average slope over all directions and scales
sumn = 0;
for k = 6:(2 * NUM_RAD)
    if radCount(k) > 0
        sumn = sumn + 1;
        yval(sumn) = radius(k) / radCount(k);
        tempr(sumn) = (k -1) * rmax / (2 * NUM_RAD);
    end
end
p = polyfit(tempr,yval,1);
averslope = p(1);
averIC = p(2);

fitln = polyval(p,tempr);
figure; plot(tempr,yval,tempr,fitln,'r-');
title('Log Log plot of Magn. vs Freq.');
ylabel('Log Magnitude');
xlabel('Log Frequency');

slope(NUM_DIR + 1) = slope(1);
intercept(NUM_DIR + 1) = intercept(1);

%draw rose plot of slope and intercept
ang = 1: (NUM_DIR + 1);
figure;
polar(pi/NUM_DIR + (ang -1) * 2* pi / NUM_DIR, intercept(ang), 'r-');
title('Rose plot of intercept');
figure;
polar(pi/NUM_DIR + (ang - 1)* 2 * pi / NUM_DIR, abs(slope(ang)), '-');
title('Rose plot of slope');
disp(['Elapsed time: ' num2str(toc)]);


text:
function s = snrdbadd(x,n,db)

ls = length(x);

%% signal and noise power

sigpower = 10*log10(sum(x.^2)/ls);

noisepower = 10*log10(sum(n.^2)/ls);

%%

npower = sigpower-noisepower-db;

%%

s = x+sqrt(10^(npower/10)).*n;


text:
function status = checkGamma(Gamma,T,train,subj)
% Check that the state time courses are well-defined (no NaNs)
% and that the HMM inference wasn't trapped in a trivial local minima

if nargin<4, subj = 0; end

if any(isnan(sum(Gamma,2)))
    t = find(isnan(sum(Gamma,2)),1);
    if train.order>0
        d = train.order;
        Td = T-d;
    elseif length(train.embeddedlags)>1
        d = length(train.embeddedlags)-1;
        Td = T;
    else
        d = 0;
        Td = T;
    end
    cT = cumsum(Td);
    i = find(cT>t,1);
    t = t + i*d;
    if length(train.embeddedlags)>1, t = t - max(train.embeddedlags); end
    if subj > 0
        str = [' in subject ' num2str(subj) ', '];
    else
        str = ', ';
    end
    disp(['The state time courses have NaN in t=' num2str(t) str 'where t indexes time in the data.'])
    disp('This is probably due to an artifact or an event too extreme to model around that time point.')
    disp('Please check your data arount that time point, and remove or smooth this event')
	error('Issue of precision')
end

K = size(Gamma,2);
if train.episodic % no checking
    status = 0;
else
    status = (all(max(Gamma)<0.6) && all(min(Gamma)>(1/K/2)));
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function Q = IsentropicVortexIC2D(x, y, time)

% function Q = IsentropicVortexIC2D(x, y)
% Purpose: compute flow configuration given by
%     Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159

% based flow parameters
xo = 5; yo = 0; beta = 5; gamma = 1.4;
rho = 1; u = 1; v = 0; p = 1;

xmut = x-u*time; ymvt = y-v*time;
r = sqrt((xmut-xo).^2 + (ymvt-yo).^2);

% perturbed density
u   = u - beta*exp(1-r.^2).*(ymvt-yo)/(2*pi);
v   = v + beta*exp(1-r.^2).*(xmut-xo)/(2*pi);
rho1 = (1 - ((gamma-1)*beta^2*exp(2*(1-r.^2))/(16*gamma*pi*pi))).^(1/(gamma-1));
p1   = rho1.^gamma;

Q(:,:,1) = rho1; Q(:,:,2) = rho1.*u; Q(:,:,3) = rho1.*v;
Q(:,:,4) = p1/(gamma-1) + 0.5*rho1.*(u.^2 + v.^2);
return;

repo: Aero-matlab


text:
function [V_sorted, S_sorted] = sorted_eig(X, direction)
%SORTED_EIG Compute the sorted eigenvalue decomposition of a square matrix
%
%   Inputs:
%       X:          matrix to decompose
%       direction:  'ascend','descend' for ordering of eigenvectors and
%           eigenvalues
%
%   Outputs:
%       V_sorted:   matrix of eigenvectors sorted according to the
%           eigenvalues
%       S_sorted:   diagonal matrix of sorted eigenvalues
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SORTED_EIG.M - 5/10/2016
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%

if nargin<2, direction = 'descend'; end

if size(X,1) ~= size(X,2)
    error('input matrix should be square')
end

[V, S] = eig(X);
[~, perm] = sort(diag(S), 1, direction);
S_sorted = S(perm, perm); V_sorted = V(:, perm);

end

author: polarch

repo: Spherical-Array-Processing


text:
function add_io(hRD,design,board)

% add AXI4 and AXI4-Lite slave interfaces
switch lower(board)
	case 'ccfmc_lvds'
		hRD.addAXI4SlaveInterface( ...
		    'InterfaceConnection', 'axi_cpu_interconnect/M11_AXI', ...
		    'BaseAddress',         '0x43C00000', ...
		    'MasterAddressSpace',  'sys_ps7/Data');
    case 'ccbob_cmos'
		hRD.addAXI4SlaveInterface( ...
		    'InterfaceConnection', 'axi_cpu_interconnect/M06_AXI', ...
		    'BaseAddress',         '0x43C00000', ...
		    'MasterAddressSpace',  'sys_ps7/Data');
    case 'ccbob_lvds'
		hRD.addAXI4SlaveInterface( ...
		    'InterfaceConnection', 'axi_cpu_interconnect/M06_AXI', ...
		    'BaseAddress',         '0x43C00000', ...
		    'MasterAddressSpace',  'sys_ps7/Data');
    case 'ccbox_lvds'
		hRD.addAXI4SlaveInterface( ...
		    'InterfaceConnection', 'axi_cpu_interconnect/M05_AXI', ...
		    'BaseAddress',         '0x43C00000', ...
		    'MasterAddressSpace',  'sys_ps7/Data');
end

if ~strcmp(design,'Tx')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Rx Reference design interfaces
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data Valid OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'dut_data_valid', ...
        'PortWidth',      1, ...
        'InterfaceConnection', 'util_ad9361_adc_pack/adc_valid_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 0 OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'dut_data_0', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 1 OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'dut_data_1', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_1', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 2 OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'dut_data_2', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_2', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 3 OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'dut_data_3', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_3', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 ADC Data I0', ...
        'InterfaceType',  'IN', ...
        'PortName',       'sys_wfifo_0_dma_wdata', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 ADC Data Q0', ...
        'InterfaceType',  'IN', ...
        'PortName',       'sys_wfifo_1_dma_wdata', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_1', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 ADC Data I1', ...
        'InterfaceType',  'IN', ...
        'PortName',       'sys_wfifo_2_dma_wdata', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_2', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 ADC Data Q1', ...
        'InterfaceType',  'IN', ...
        'PortName',       'sys_wfifo_3_dma_wdata', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_3', ...
        'IsRequired',     false);
end

if ~strcmp(design,'Rx')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Tx Reference design interfaces
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 DAC Data I0', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'axi_ad9361_dac_data_i0', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 DAC Data Q0', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'axi_ad9361_dac_data_q0', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_1', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 DAC Data I1', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'axi_ad9361_dac_data_i1', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_2', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 DAC Data Q1', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'axi_ad9361_dac_data_q1', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_3', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 0 IN', ...
        'InterfaceType',  'IN', ...
        'PortName',       'util_dac_unpack_dac_data_00', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 1 IN', ...
        'InterfaceType',  'IN', ...
        'PortName',       'util_dac_unpack_dac_data_01', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_1', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 2 IN', ...
        'InterfaceType',  'IN', ...
        'PortName',       'util_dac_unpack_dac_data_02', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_2', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 3 IN', ...
        'InterfaceType',  'IN', ...
        'PortName',       'util_dac_unpack_dac_data_03', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_3', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Load Tx Data OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'util_dac_unpack_dac_valid_00', ...
        'PortWidth',      1, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/dac_valid_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Valid Tx Data IN', ...
        'InterfaceType',  'IN', ...
        'PortName',       'util_dac_unpack_upack_valid_00', ...
        'PortWidth',      1, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/upack_valid_0', ...
        'IsRequired',     false);

end

author: analogdevicesinc


text:
function [ n, indx ] = r8vec_indexed_heap_d_insert ( n, a, indx, indx_insert )

%*****************************************************************************80
%
%% R8VEC_INDEXED_HEAP_D_INSERT: insert value into heap descending indexed R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    An indexed R8VEC is an R8VEC of data values, and an R8VEC of N indices,
%    each referencing an entry of the data vector.
%
%    Note that the argument N must be a variable, and will be incremented before
%    return, and that INDX must be able to hold one more entry on output than
%    it held on input.
%
%    This is one of three functions needed to model a priority queue.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Thomas Cormen, Charles Leiserson, Ronald Rivest,
%    Introduction to Algorithms,
%    MIT Press, 2001,
%    ISBN: 0262032937,
%    LC: QA76.C662.
%
%  Parameters:
%
%    Input, integer N, the number of items in the
%    index vector.
%
%    Input, real A(*), the data vector.
%
%    Input, integer INDX(N), the index vector.
%
%    Input, integer INDX_INSERT, the index in A of the value
%    to be inserted into the heap.
%
%    Output, integer N, the number of items in the revised index vector.
%
%    Output, integer INDX(N), the revised index vector.
%
  n = n + 1;
  i = n;

  while ( 1 < i )

    parent = floor ( i / 2 );

    if ( a(indx_insert) <= a(indx(parent)) )
      break
    end

    indx(i) = indx(parent);
    i = parent;

  end

  indx(i) = indx_insert;

  return
end


text:
function [fL] = bruteloglike(vValues,time_as)
    % bruteloglike calculates the log likelihood function for the modeled aftershock sequence and the maximum likelihood estimate for k, c and p
    %
    % [fL] = bruteloglike(vValues,time_as);
    % -------------------------------------------------------------
    % Reference: Ogata, Estimation of the parameters in the modified Omori formula
    % for aftershock sequences by  the maximum likelihood procedure, J. Phys. Earth, 1983
    % (Formula 6)
    %
    % J. Woessner
    % updated: 29.07.03

    p = vValues(1);
    c = vValues(2);
    k = vValues(3);

    % Setting start end end time
    fTstart = min(time_as);
    fTend = max(time_as);

    if p ~= 1
        fAcp = ((fTend+c).^(1-p)-(fTstart+c).^(1-p))./(1-p);
        %cumnr_model = k/(p-1)*(c^(1-p)-(c+time_as(i)).^(1-p)); % integrated form of MOL
    else
        fAcp = log(fTend+c)-log(fTstart+c);
        %cumnr_model = k*log(time_as(i)/c+1); % integrated form of MOL
    end
    % rms = (sum((i-cumnr_model).^2)/length(i))^0.5; % RMS between observed data and MOL
    % Log likelihood
    nNumEvents = length(time_as);
    fL = -(nNumEvents*log(k)-p*sum(log(time_as+c))-k*fAcp);

author: CelsoReyes


text:
function boolVal=allElementsAreInts(x)
%%ALLELEMENTSAREINTS Given a vector or a matrix, this function returns true
%            if all of the elements in the matrix are exactly integers
%            (including complex integers. This differs from the isinteger
%            function that is built into Matlab as this actually checks the
%            values of the data, whereas isinteger only tells whether the
%            data type is an integer. THis function can handle floating
%            point numbers.
%
%INPUTS: x A scalar or matrix/hypermatrix where one wants to determine
%          whether all elements are integers.
%
%OUTPUTS: boolVal This is true if all elements of x are integer values
%               (regardless of data type) and this is false otherwise.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

boolVal=all(x(:)==fix(x(:)));

end

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


text:
function pass = test_ivp(varargin)

% This test solves the Van der Pol ODE in CHEBFUN using ode15s and ode45. It
% checks the solution with Matlab's inbuilt ode15s and ode45 solvers.

% Rodrigo Platte Jan 2009, modified by Asgeir Birkisson, modified by Nick Hale.

%% Test using Van der Pol:

% Test ode15s Using default tolerances (RelTol = 1e-3)
[t, y] = chebfun.ode15s(@vdp1, [0, 5], [2 ; 0]); % CHEBFUN solution
[tm, ym] = ode15s(@vdp1, [0, 5], [2 ; 0]);       % Matlab's solution
pass(1) = max(max(abs(ym - feval(y, tm)))) < 2e-2;

% Test ode45 Using default tolerances (RelTol = 1e-3)
[t, y] = chebfun.ode45(@vdp1, [0, 5], [2 ; 0]);  % CHEBFUN solution
[tm, ym] = ode45(@vdp1, [0, 5], [2 ; 0]);        % Matlab's solution
pass(2) = max(max(abs(ym - feval(y, tm)))) < 1e-2;

% Test ode113 with different tolerance (RelTol = 1e-6)
opts = odeset('RelTol', 1e-6);
[t, y] = chebfun.ode113(@vdp1, [0, 20], [2 ; 0], opts); % CHEBFUN solution
[tm, ym] = ode113(@vdp1, [0, 20], [2 ; 0], opts);       % Matlab's solution
pass(3) = max(max(abs(ym - feval(y,tm)))) < 1e-5;

%% Test some trivial complex-valued IVPs:

f = @(x, u) 1i*u;
d = [0, 1];
soln = exp(1i);

% Test ode15s:
y = chebfun.ode15s(f, d, 1);
pass(4) = abs(y(1) - soln) < 2e-2;

% Test ode45:
y = chebfun.ode45(f, d, 1);
pass(5) = abs(y(1) - soln) < 2e-2;

% Test ode113:
y = chebfun.ode113(f, d, 1);
pass(6) = abs(y(1) - soln) < 2e-2;

end


text:
function [Geometry,Geom_res] =fbkp_grid(theta,N,interp)

% theta are imput as angles in degrees
% Define the x & y axes for the reconstructed image so that the origin
% (center) is in the spot which RADON would choose.

%   N is a scalar that specifies the number of rows and columns in the
%   reconstructed image. N is determined from the length of the projections: !!!!
%
%       N = 2*floor(size(P,1)/(2*sqrt(2)))
%
theta = pi*theta/180;

xax = (1:N)-ceil(N/2);
x = repmat(xax,N,1) ;
% x coordinates, the y coordinates are rot90(x)
y = rot90(x);
costheta = cos(theta);
sintheta = sin(theta);

% Backprojection - vectorized in (x,y), looping over theta
sa=size(x,1); sb=size(x,2); sc=length(theta);

if strcmp(interp, 'nearest neighbor')
        Geometry=repmat(int16(0),[sa,sb,sc]);
        Geom_res=[];
    for i=1:length(theta)
        t=[];
        t = round( x*costheta(i) + y*sintheta(i) );
        Geometry(:,:,i)=t;
    end

elseif strcmp(interp, 'linear')
    Geometry=repmat(int16(0),[sa,sb,sc]); % parte intera
    Geom_res=repmat(uint8(0),[sa,sb,sc]); % parte decimale

    for i=1:length(theta)
        t = x.*costheta(i) + y.*sintheta(i);
        a = floor(t);
        % img = img + (t-a).*proj(a+1+ctrIdx) + (a+1-t).*proj(a+ctrIdx);
        Geometry(:,:,i)=a ; % floor
        Geom_res(:,:,i)=uint8( abs(t-a)*100) ; % decimal i.e. double - floor

    end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function disp(F)
%DISP   Display a DISKFUNV.
%
% See also DISPLAY.

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

loose = strcmp(get(0,'FormatSpacing'),'loose');

% Compact version
if ( isempty( F ) )
    fprintf('empty diskfunv\n\n')
    return
end

% Title in display:
disp('diskfunv object containing')
if ( loose )
    fprintf('\n');
end

% Display the two DISKFUN components.
for j = 1:F.nComponents
    disp( F.components{j} );
end

end


text:
% Set axes limit.

clear all;
addpath('../lib');

%% lets plot 3 cycles of 50Hz AC voltage
f = 50;
Vm = 10;
phi = 0;

% generate the signal
t = [0:0.0001:3/f];
th = 2*pi*f*t;
v = Vm*sin(th+phi);

%% plot now
plt = Plot(t*1E3, v);

plt.XLabel = 'Time, t (ms)'; % xlabel
plt.YLabel = 'Voltage, V (V)'; %ylabel
plt.XLim = [0, 40]; % set x axis limit
plt.YLim = [-11, 11]; % set y axis limit

% Save? comment the following line if you do not want to save
plt.export('plotAxisLimit.jpg');

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub


text:
% WAVELET_FACTORY_3D Build roto-translation wavelet operators
%
% Usage
%   [Wop, filters, filters_rot] = WAVELET_FACTORY_3D(size_in, filt_opt, filt_rot_opt, scat_opt)
%
% Input
%    size_in (2x1 int): the size of the image to be transformed
%    filt_opt (struct): the filter options, same as for MORLET_FILTER_BANK_2D
%    filt_rot_opt (struct): the filter options for the
%		filters along angular parameter, same as for MORLET_FILTER_BANK_1D
%	 scat_opt (struct): the scattering and wavelet options, same as
%		WAVELET_LAYER_3D
%
% Output
%    Wop: A cell array of wavelet transforms needed for the scattering trans-
%       form.
%    filters: A cell array of the filters used in defining the wavelets.
%
% Description
%   This function builds the wavelet operators used to compute the
%   roto-translation scattering. The first operator is a 2d wavelet
%   transform obtained with WAVELET_LAYER_2D. The second operator and third
%   operators are roto-translation wavelet transform obtained with
%   WAVELET_LAYER_3D.
%
% See also
%   SCAT, WAVELET_2D, WAVELET_LAYER_2D, WAVELET_3D, WAVELET_3D_LAYER

function [ Wop, filters, filters_rot ] = ...
		wavelet_factory_3d(size_in, filt_opt, filt_rot_opt, scat_opt)

    % initialize the non-provided options with empty structure
    if (nargin < 4)
        scat_opt = struct();
    end
    if (nargin < 3)
        filt_rot_opt = struct();
    end
    if (nargin < 2)
        filt_opt = struct();
    end

	% filters along spatial variable
	filters = morlet_filter_bank_2d(size_in, filt_opt);

	% filters along angular variable
	sz = filters.meta.L * 2; % L orientations between 0 and pi
	% periodic convolutions along angles
	filt_rot_opt.boundary = 'per';
	filt_rot_opt.filter_format = 'fourier_multires';
	filt_rot_opt.J = 3;
	filt_rot_opt.P = 0;
	filters_rot = morlet_filter_bank_1d(sz, filt_rot_opt);

	% number of layer
	scat_opt = fill_struct(scat_opt, 'M', 2);
    wav_opt = rmfield(scat_opt, 'M');

	% first wavelet transform is an usual wavelet transform
	Wop{1} = @(x)(wavelet_layer_2d(x, filters, wav_opt));

	% other wavelet transform are roto-translation wavelet
	for m = 2:scat_opt.M+1
		Wop{m} = @(x)(wavelet_layer_3d(x, filters, filters_rot, wav_opt));
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function test_model = testNRSFM(data, train_model)

globals;
params = get_params();

fprintf('Testing NRSfM model\n');

% Generate test set
[~,test] = prep_data(data.test,0,0);

% Take full train labels and test on test set. Just reorder.
test = normalize_test(test, train_model);
nTest = size(test.points,1)/2;
nTestIters = params.nrsfm.test_em_iters;

% Fit NRSFM parameters parallely to instances
p =  TimedProgressBar( nTest, 30, ...
'Testing NRSfM ', ', completed ', 'Testing NRSfM concluded in ' );
parfor i=1:nTest
    test_model_arr(i) = em_single_instance(train_model,test,i,nTestIters);
    p.progress();
end
p.stop();

% Merge test models into a single test_model
test_model = merge_model_arr(test_model_arr);
if(isfield(test,'subtype'))
    test_model.subtype = test.subtype;
end
if(isfield(test,'rotP3d'))
    test_model.rotP3d = test.rotP3d;
end
if(isfield(test,'rotationPred'))
    test_model.rotationPred = test.rotationPred;
end
if(isfield(test,'subtypePred'))
    test_model.subtypePred = test.subtypePred;
end

end

repo: CategoryShapes


text:
function x = p41_start ( n )

%*****************************************************************************80
%
%% P41_START returns a starting point for optimization for problem 41.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = [ 0.5, 1.0, -0.5, -1.0 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%ISFIXED Test on fixed mapping
%
%    I = ISFIXED(W)
%    ISFIXED(W)
%
% True if the mapping type of W is 'fixed' (see HELP MAPPINGS). If called
% without an output argument ISFIXED generates an error if the mapping type
% of W is not 'fixed'.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
disp(' ');
disp('fdct_wrapping_demo_basic.m -- Displays a curvelet')
disp ('both in the spatial and frequency domains.');
disp(' ');
disp(['This is achieved by setting all the coefficients in the curvelet'])
disp(['domain to zero except that at the required location (which'])
disp(['is set to one). The curvelet is obtained by taking the'])
disp(['adjoint curvelet transform. Notice how the curvelet is sharply '])
disp(['localized in both space and frequency.']);
disp(' ');

% fdct_wrapping_demo_basic.m -- Displays a curvelet both in the spatial and frequency domains.

% m = 512;
% n = 512;
%
% X = zeros(m,n);

%forward curvelet transform
disp('Take curvelet transform: fdct_wrapping');
tic; C = fdct_wrapping(X,0); toc;

%specify one curvelet
s = 5;
w = 1;
[A,B] = size(C{s}{w});6
a = ceil((A+1)/2);
b = ceil((B+1)/2);
C{s}{w}(a,b) = 1;

%adjoint curvelet transform
disp('Take adjoint curvelet transform: ifdct_wrapping');
tic; Y = ifdct_wrapping(C,0); toc;

%display the curvelet
F = ifftshift(fft2(fftshift(Y)));
subplot(1,2,1); colormap gray; imagesc(real(Y)); axis('image'); ...
    title('a curvelet: spatial viewpoint');
subplot(1,2,2); colormap gray; imagesc(abs(F)); axis('image'); ...
    title('a curvelet: frequency viewpoint');

%get parameters
[SX,SY,FX,FY,NX,NY] = fdct_wrapping_param(C);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function adjList = meshFaceAdjacency(vertices, edges, faces)
%MESHFACEADJACENCY Compute adjacency list of face around each face.
%
%
%   Example
%     % Create a sample 3D mesh
%     [v, e, f] = createDodecahedron;
%     adjList = meshFaceAdjacency(v, e, f);
%     figure; hold on; axis equal; view([100 40]);
%     drawMesh(v, f);
%     % draw sample face in a different color
%     drawMesh(v, f(1, :), 'faceColor', 'b');
%     % draw the neighbors of a sample face
%     drawMesh(v, f(adjList{1}, :), 'faceColor', 'g')
%
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2010-10-04, using Matlab 7.9.0.529 (R2009b)
% Copyright 2010-2022 INRA - Cepia Software Platform

edgeFaceList = meshEdgeFaces(vertices, edges, faces);

% allocate memory for adjacency list
nFaces = max(edgeFaceList(:));
adjList = cell(1, nFaces);

% iterate over edges to populate adjacency list
for iEdge = 1:size(edgeFaceList)
    f1 = edgeFaceList(iEdge, 1);
    f2 = edgeFaceList(iEdge, 2);
    adjList{f1} = [adjList{f1} f2];
    adjList{f2} = [adjList{f2} f1];
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function pdepetrans
% 1D transport - modelling with extensions for decay and fast sorption
%    using MATLAB pdepe
%
%   $Ekkehard Holzbecher  $Date: 2006/03/16 $
%--------------------------------------------------------------------------
T = 1;                     % maximum time [s]
L = 1;                     % length [m]
D = 1;                     % diffusivity [m*m/s]
v = 1;                     % velocity [m/s]
lambda = 0.0;              % decay constant [1/s]
sorption = 2;              % sorption-model: no sorption (0), linear (1),
                           %                 Freundlich (2), Langmuir (3)
k1 = 0.0004;               % sorption parameter 1 (R=0 for linear isotherm with Kd, else k1=R)
k2 = 0.5;                  % sorption parameter 2 (Kd for linear isotherm with Kd)
rhob = 1300;               % porous medium bulk density [kg/m*m*m]
theta = 0.2;               % porosity [-]
c0 = 0.0;                  % initial concentration [kg/m*m*m]
cin = 1;                   % boundary concentration [kg/m*m*m]

M = 40;                    % number of timesteps
N = 40;                    % number of nodes
%-------------------------- output parameters
gplot = 0;                 % =1: breakthrough curves; =2: profiles
gsurf = 0;                 % surface
gcont = 0;                 % =1: contours; =2: filled contours
ganim = 2;                 % animation of profiles; =1: single line; =2: all lines

t = linspace (T/M,T,M);    % time discretization
x = linspace (0,L,N);      % space discretization

%----------------------execution-------------------------------------------
if sorption == 1 && k1 <=0
    k1 = 1+k2*rhob/theta;
else
    if sorption > 1 k1 = rhob*k1/theta; end
end
options = odeset; if (c0 == 0) c0 = 1.e-20; end
c = pdepe(0,@transfun,@ictransfun,@bctransfun,x,[0 t],options,D,v,lambda,sorption,k1,k2,c0,cin);

%---------------------- graphical output ----------------------------------
switch gplot
    case 1
        plot ([0 t],c)        % breakthrough curves
        xlabel ('time'); ylabel ('concentration');
    case 2
        plot (x,c','--')      % profiles
        xlabel ('space'); ylabel ('concentration');
end
if gsurf                      % surface plot
    figure; surf (x,[0 t],c);
    xlabel ('space'); ylabel ('time'); zlabel('concentration');
end
if gcont figure; end
switch gcont
    case 1
        contour (x,[0 t],c)   % contours
        grid on; xlabel ('space'); ylabel ('time');
    case 2
        contourf(x,[0 t],c)   % filled contours
        colorbar; xlabel ('space'); ylabel ('time');
end
if (ganim)
    [FileName,PathName] = uiputfile('*.mpg');
    figure; if (ganim > 1) hold on; end
    for j = 1:size(c,1)
        axis manual;  plot (x,c(j,:),'r','LineWidth',2);
        ylim ([min(c0,cin) max(c0,cin)]);
        legend (['t=' num2str(T*(j-1)/M)]);
        Anim(j) = getframe;
        plot (x,c(j,:),'b','LineWidth',2);
    end
    mpgwrite (Anim,colormap,[PathName '/' FileName]);     % mgwrite not standard MATLAB
    movie (Anim,0);   % play animation
end

%----------------------functions------------------------------
function [c,f,s] = transfun(x,t,u,DuDx,D,v,lambda,sorption,k1,k2,c0,cin)
switch sorption
    case 0
        R = 1;
    case 1
        R = k1;
    case 2
        R = 1+k1*k2*u^(k2-1);
    case 3
        R = 1+k1*k2*u/(k2+u)/(k2+u);
end
c = R;
f = D*DuDx;
s = -v*DuDx -lambda*R*u;
% --------------------------------------------------------------
function u0 = ictransfun(x,D,v,lambda,sorption,k1,k2,c0,cin)
u0 = c0;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = bctransfun(xl,ul,xr,ur,t,D,v,lambda,sorption,k1,k2,c0,cin)
pl = ul-cin;
ql = 0;
pr = 0;
qr = 1;

author: Sable


text:
function [ uout, vout ] = triasimp ( x, y )

%*****************************************************************************80
%
%% TRIASIMP maps a point from the reference triangle to the simplex.
%
%  Discussion:
%
%    Map the reference triangle with vertices
%      (-1,-1/sqrt(3)), (1,-1/sqrt(3)), (0,2/sqrt(3))
%    to the simplex with vertices
%      (0,0), (1,0), (0,1).
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    26 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, real X, Y, the coordinates of a point in the
%    reference triangle.
%
%    Output, real UOUT, VOUT, the coordinates of the corresponding
%    point in the simplex.
%
  scale = 1.0 / sqrt ( 3.0 );

  uout = 0.5 * ( x + 1.0 ) - 0.5 * scale * ( y + scale );

  vout = 1.0 * scale * ( y + scale );

  return
end


text:
%%                   testinitWeighted.m

% This test file implements the Weighted initializer. The code builds
% a synthetic formulation of the Phase Retrieval problem, b = |Ax| and
% computes an estimate to x. The code finally outputs the correlation
% achieved.

% PAPER TITLE:
%              Solving Almost all Systems of Random Quadratic Equations.

% ARXIV LINK:
%              https://arxiv.org/pdf/1705.10407.pdf

% This is a test script for the weighted initializer.

% 1.) Each test script for an initializer starts out by defining the length
% of the unknown signal, n and the number of measurements, m. These
% mesurements can be complex by setting the isComplex flag to be true.

% 2.) We then build the test problem by generating random gaussian
% measurements and using b0 = abs(Ax) where A is the measurement matrix.

% 3.) We run x0 = initSpectral(A,[],b0,n,isTruncated,isScaled) which runs
% the initialier and and recovers the test vector with high correlation.

% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% ---------------------------------START-----------------------------

clc
clear
close all

% Parameters
n = 500;           % number of unknowns
m = 20*n;          % number of measurements
isComplex = true;  % use complex matrices? or just stick to real?

%  Build the test problem
xt = randn(n,1)+isComplex*randn(n,1)*1i; % true solution
A = randn(m,n)+isComplex*randn(m,n)*1i;  % matrix
b0 = abs(A*xt);                          % data

% Invoke the truncated spectral initial method
x0 = initWeighted(A,[],b0,n);

% Calculate the correlation between the recovered signal and the true signal
correlation = abs(x0'*xt/norm(x0)/norm(xt));

fprintf('correlation: %f\n', correlation);

author: tomgoldstein

repo: phasepack-matlab


text:
function c = he_polynomial_coefficients ( n )

%*****************************************************************************80
%
%% HE_POLYNOMIAL_COEFFICIENTS: coefficients of He(i,x).
%
%  Discussion:
%
%    He(i,x) represents the probabilist's Hermite polynomial.
%
%  First terms:
%
%    N/K     0     1      2      3       4     5      6    7      8    9   10
%
%     0      1
%     1      0     1
%     2     -1     0      1
%     3      0    -3      0      1
%     4      3     0     -6      0       1
%     5      0    15      0    -10       0     1
%     6    -15     0     45      0     -15     0      1
%     7      0  -105      0    105       0   -21      0     1
%     8    105     0   -420      0     210     0    -28     0      1
%     9      0   945      0  -1260       0   378      0   -36      0   1
%    10   -945     0   4725      0   -3150     0    630     0    -45   0    1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N, the highest order polynomial to compute.
%    Note that polynomials 0 through N will be computed.
%
%    Output, real C(1:N+1,1:N+1), the coefficients of the Hermite polynomials.
%
  if ( n < 0 )
    c = [];
    return
  end

  c(1:n+1,1:n+1) = 0.0;

  c(1,1) = 1.0;

  if ( n == 0 )
    return
  end

  c(2,2) = 1.0;

  for i = 2 : n
    c(i+1,1)     =                             - ( i - 1 ) * c(i-1,1);
    c(i+1,2:i-1) =                c(i  ,1:i-2) - ( i - 1 ) * c(i-1,2:i-1);
    c(i+1,  i  ) =                c(i  ,  i-1);
    c(i+1,  i+1) =                c(i  ,  i  );
  end

  return
end


text:
% Script to reproduce the experiments leading to the results provided in the
% Table 2 of the paper "Deep Scattering Spectrum" by J. Andén and S. Mallat.

% M=2 scattering, frequency scattering, multiple Q1, multiple T

run_name = 'DSS_Table2_GTZAN_m2_freq_multQ1_multT';

src = gtzan_src('/path/to/gtzan');

N = 5*2^17;

filt1_opt.filter_type = {'gabor_1d','morlet_1d'};
filt1_opt.Q = [8 1];
filt1_opt.J = T_to_J(8192,filt1_opt);

sc1_opt.M = 2;

ffilt1_opt.filter_type = 'morlet_1d';
ffilt1_opt.J = 7;

fsc1_opt.M = 1;

Wop1 = wavelet_factory_1d(N, filt1_opt, sc1_opt);
fWop1 = wavelet_factory_1d(128, ffilt1_opt, fsc1_opt);

scatt_fun1 = @(x)(log_scat(renorm_scat(scat(x,Wop1))));
fscatt_fun1 = @(x)(func_output(@scat_freq,2,scatt_fun1(x),fWop1));
feature_fun1 = @(x)(format_scat(fscatt_fun1(x)));

filt2_opt = filt1_opt;
filt2_opt.Q = [1 1];
filt2_opt.J = T_to_J(8192,filt2_opt);

sc2_opt = sc1_opt;

ffilt2_opt = ffilt1_opt;
ffilt2_opt.J = 5;

fsc2_opt = fsc1_opt;

Wop2 = wavelet_factory_1d(N, filt2_opt, sc2_opt);
fWop2 = wavelet_factory_1d(32, ffilt2_opt, fsc2_opt);

scatt_fun2 = @(x)(log_scat(renorm_scat(scat(x,Wop2))));
fscatt_fun2 = @(x)(func_output(@scat_freq,2,scatt_fun2(x),fWop2));
feature_fun2 = @(x)(format_scat(fscatt_fun2(x)));

filt3_opt = filt1_opt;
filt3_opt.J = T_to_J(2*8192,filt3_opt);

sc3_opt = sc1_opt;

Wop3 = wavelet_factory_1d(N, filt3_opt, sc3_opt);

scatt_fun3 = @(x)(log_scat(renorm_scat(scat(x,Wop3))));
fscatt_fun3 = @(x)(func_output(@scat_freq,2,scatt_fun3(x),fWop1));
feature_fun3 = @(x)(format_scat(fscatt_fun3(x)));

filt4_opt = filt2_opt;
filt4_opt.J = T_to_J(2*8192,filt4_opt);

sc4_opt = sc2_opt;

Wop4 = wavelet_factory_1d(N, filt4_opt, sc4_opt);

scatt_fun4 = @(x)(log_scat(renorm_scat(scat(x,Wop4))));
fscatt_fun4 = @(x)(func_output(@scat_freq,2,scatt_fun4(x),fWop2));
feature_fun4 = @(x)(format_scat(fscatt_fun4(x)));

filt5_opt = filt1_opt;
filt5_opt.J = T_to_J(4*8192,filt5_opt);

sc5_opt = sc1_opt;

Wop5 = wavelet_factory_1d(N, filt5_opt, sc5_opt);

scatt_fun5 = @(x)(log_scat(renorm_scat(scat(x,Wop5))));
fscatt_fun5 = @(x)(func_output(@scat_freq,2,scatt_fun5(x),fWop1));
feature_fun5 = @(x)(format_scat(fscatt_fun5(x)));

filt6_opt = filt2_opt;
filt6_opt.J = T_to_J(4*8192,filt6_opt);

sc6_opt = sc2_opt;

Wop6 = wavelet_factory_1d(N, filt6_opt, sc6_opt);

scatt_fun6 = @(x)(log_scat(renorm_scat(scat(x,Wop6))));
fscatt_fun6 = @(x)(func_output(@scat_freq,2,scatt_fun6(x),fWop2));
feature_fun6 = @(x)(format_scat(fscatt_fun6(x)));

features = {feature_fun1, feature_fun2, feature_fun3, ...
	feature_fun4, feature_fun5, feature_fun6};

for k = 1:length(features)
    fprintf('testing feature #%d...',k);
    tic;
    sz = size(features{k}(randn(N,1)));
    aa = toc;
    fprintf('OK (%.2fs) (size [%d,%d])\n',aa,sz(1),sz(2));
end

db = prepare_database(src,features);
db.features = single(db.features);
db = svm_calc_kernel(db,'gaussian','square',1:2:size(db.features,2));

rs = RandStream.create('mt19937ar','Seed',floor(pi*1e9));
RandStream.setGlobalStream(rs);
[train_set{1}, test_set{1}] = create_partition([src.objects.class], 0.9);
for k = 2:10
	[train_set{k}, test_set{k}] = ...
		next_fold([src.objects.class], train_set{k-1}, test_set{k-1});
end

optt.kernel_type = 'gaussian';
optt.C = 2.^[0:4:8];
optt.gamma = 2.^[-16:4:-8];
optt.search_depth = 3;
optt.full_test_kernel = 0;

for k = 1:10
	[dev_err_grid,C_grid,gamma_grid] = ...
		svm_adaptive_param_search(db,train_set{k},[],optt);

	[dev_err(k),ind] = min(mean(dev_err_grid{end},2));
	C(k) = C_grid{end}(ind);
	gamma(k) = gamma_grid{end}(ind);

	optt1 = optt;
	optt1.C = C(k);
	optt1.gamma = gamma(k);

	model = svm_train(db,train_set{k},optt1);
	labels(:,k) = svm_test(db,model,test_set{k});
	err(k) = classif_err(labels(:,k),test_set{k},db.src);

	fprintf('dev err = %f, test err = %f\n',dev_err(k),err(k));

	save([run_name '.mat'],'labels','dev_err','err','C','gamma');
end