text:
function [ output_data ] = normalize_complex_file( input_filename, varargin )
%NORMALIZE_COMPLEX_FILE Normalize complex SAR data to make it have consistent
%characteristics in frequency domain
%
% Normalized complex data will have
%    1) Centered and deskewed frequency support
%    2) Uniform weighting
%    3) FFT sign of -1
%
% This function only applies these transforms in one dimension since we
% normally only process one dimension at a time, and since in general,
% frequency support deskew can only be applied in one dimension.
%
%    output_data = normalize_complex_file(input_filename, 'PropertyName', PropertyValue, ...)
%
%       Property name     Description
%
%       output_filename   Output data is stored to a complex single-
%                            precision float SIO file.  Default is to not
%                            save to file.
%       block_size        Determines how many lines are processed in memory
%                            at once.  Set to Inf to process whole image at
%                            once.  Default is 500 lines.
%       dim               Dimension over which to perform normalization.
%                            Default = 1 (slow-time).
%       azlimits          Min and max samples in azimuth over which to
%                            compute.  Default = full image.
%       rnglimits         Min and max lines in range over which to compute.
%                            Default = full image.
%       framenumber       Dataset to process if INPUT_FILENAME is a
%                            multi-image file.  Default is 1.
%
% If an output argument is given to the function call, the normalized
% complex data will be returned in an in-memory MATLAB array to
% OUTPUT_DATA. WARNING: Make sure your computer has enough memory to hold
% resulting array, if you use this option.
%
% Author: Wade Schwartzkopf, NGA/R
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Parse input parameters
p = inputParser; % Extract parameter-value pairs
p.KeepUnmatched=true;
p.addParamValue('output_filename','');
p.addParamValue('block_size',500);
p.addParamValue('dim',1);
p.addParamValue('azlimits',[]);
p.addParamValue('rnglimits',[]);
p.addParamValue('framenumber',1);
p.FunctionName = 'normalize_complex_file';
p.parse(varargin{:});

save_to_sio=~any(strcmp(p.UsingDefaults,'output_filename'))||~isempty(p.Results.output_filename);
save_to_variable=(nargout>0);
if ~save_to_sio&&~save_to_variable % Nothing to do
    warning('NORMALIZE_COMPLEX_FILE:NO_OUTPUT','Neither output filename nor output variable was specified.  No ouput will be calculated.');
    return;
end

%% Open data file
reader_obj=open_reader(input_filename);
if iscell(reader_obj) % If file contains more than one image
    for i=1:length(reader_obj)
        if i~=p.Results.framenumber
            reader_obj{i}.close();
        end
    end
    reader_obj=reader_obj{p.Results.framenumber};
end
meta=reader_obj.get_meta();
% Calculate AOI
if any(strcmp(p.UsingDefaults,'azlimits'))
    azlimits=[1 meta.ImageData.NumCols];
else
    azlimits=p.Results.azlimits;
    if strcmp(azlimits,'full')
        azlimits=[1 meta.ImageData.NumCols];
    end
end
if any(strcmp(p.UsingDefaults,'rnglimits'))
    rnglimits=[1 meta.ImageData.NumRows];
else
    rnglimits=p.Results.rnglimits;
    if strcmp(rnglimits,'full')
        rnglimits=[1 meta.ImageData.NumRows];
    end
end

%% Calculate input parameters for normalize_complex_mem
if p.Results.dim==1
    sicd_grid_struct = meta.Grid.Col;
else
    sicd_grid_struct = meta.Grid.Row;
end
% Vectors describing range and azimuth distances from SCP (in meters) for rows and columns
[ DeltaKCOAPoly, az_coords_m, rg_coords_m ] = deskewparams( meta, p.Results.dim );
coords_m_args = {az_coords_m, rg_coords_m};
coords_ind = {azlimits(1):azlimits(end), rnglimits(1):rnglimits(end)};
other_dim = 3-p.Results.dim;
coords_m = coords_m_args{other_dim};
coords_m_args{p.Results.dim}=coords_m_args{p.Results.dim}(coords_ind{p.Results.dim});
% Zeropad
if isfield(sicd_grid_struct,'ImpRespBW')&&isfield(sicd_grid_struct,'SS')
    zeropad = max(1,1/(sicd_grid_struct.ImpRespBW*sicd_grid_struct.SS));
else
    zeropad = 1;
end
% Weighting function
try
    weight_fun = sicdweight2fun(sicd_grid_struct);
catch % Unable to determine weighting function from metadata.
    % Estimate weighting from complex data instead
    weight_fun = estimate_weighting_file(reader_obj, p.Results.dim);
end
% FFT sign
if isfield(sicd_grid_struct,'Sgn')
    fft_sign = sicd_grid_struct.Sgn;
else
    fft_sign = -1;
end

%% Setup loop for processing through image in blocks of azimuth lines
if save_to_variable
    output_data = zeros(diff(azlimits)+1,diff(rnglimits)+1);
    output_data_args = {1:size(output_data,1),1:size(output_data,2)};
end
if save_to_sio
    writer_object=open_sio_writer(p.Results.output_filename);
    if p.Results.dim==1
        writer_object.writeline=writer_object.write_row;
    else
        writer_object.writeline=writer_object.write_column;
    end
end
if p.Results.dim==1
    linelimits = rnglimits;
    read_lines_fun=@(lims) reader_obj.read_chip(azlimits,lims);
else
    linelimits = azlimits;
    read_lines_fun=@(lims) reader_obj.read_chip(lims,rnglimits);
end
wb_hand=waitbar(0,'Normalizing data...');
first_line_in_block = linelimits(1);
while(first_line_in_block<linelimits(2))
    line_indices=first_line_in_block:...
        min(linelimits(2), first_line_in_block + p.Results.block_size - 1);
    coords_m_args{other_dim} = coords_m(line_indices);
    normalized_block = normalize_complex_mem(...
        single(read_lines_fun(line_indices([1 end]))), DeltaKCOAPoly, ...
        coords_m_args{:}, weight_fun, zeropad, fft_sign, p.Results.dim);
    if save_to_variable
        output_data_args{other_dim} = line_indices - linelimits(1) + 1;
        output_data(output_data_args{:}) = normalized_block;
    end
    if save_to_sio
        writer_object.writeline(normalized_block);
    end
    first_line_in_block = line_indices(end) + 1;
    waitbar(double(first_line_in_block-linelimits(1))/double(diff(linelimits)+1),wb_hand);
end

%% Cleanup
close(wb_hand);
reader_obj.close();
if save_to_sio, writer_object.close(); end;

end

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

author: ngageoint


author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% dftfilt() - discrete Fourier filter
%
% Usage:
%   >> b = dftfilt(n,W,c,k,q)
%
% Inputs:
%   n - number of input samples
%   W - maximum angular freq. relative to n, 0 < W <= .5
%   c - cycles
%   k - oversampling
%   q - [0;1] 0->fft, 1->c cycles
%
% Authors: Sigurd Enghoff, Arnaud Delorme & Scott Makeig,
%          SCCN/INC/UCSD, La Jolla, 8/1/98

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

% 01-25-02 reformated help & license -ad

% future developments
% -------------------
% input into dftfilt:
% - lowfreq and maxfreq (of interest)
% - lowcycle and maxcyle (ex: 3 cycles at low freq and 10 cycles at maxfreq)
% - the delta in frequency: ex 0.5 Hz
% The function should: compute the number of points (len) automatically
% Warning with FFT compatibility
% Still, something has to be done about the masking so that it would be comaptible

function b = dftfilt(len,maxfreq,cycle,oversmp,wavfact)

count = 1;
for index = 1:1/oversmp:maxfreq*len/cycle % scan frequencies
	w(:,count) = j * index * cycle * linspace(-pi+2*pi/len, pi-2*pi/len, len)'; % exp(-w) is a sinus curve
	count = count+1; % -2*pi/len ensures that we really scan from -pi to pi without redundance (-pi=+pi)
end;
b = exp(-w);

%srate = 2*pi/len;						    % Angular increment.
%w = j * cycle * [0:srate:2*pi-srate/2]';	% Column.
%x = 1:1/oversmp:maxfreq*len/cycle;		    % Row.
%b = exp(-w*x);					            % Exponentiation of outer product.

for i = 1:size(b,2),
	m  = round(wavfact*len*(i-1)/(i+oversmp-1));	% Number of elements to discard.
	mu = round(m/2);				                % Number of upper elemnts.
	ml = m-round(m/2);				                % Number of lower elemnts.
	b(:,i) = b(:,i) .* [zeros(mu,1) ; hanning(len-m) ; zeros(ml,1)];
end

% syemtric hanning function
function w = hanning(n)
if ~rem(n,2)
   w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
   w = [w; w(end:-1:1)];
else
   w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
   w = [w; w(end-1:-1:1)];
end

author: PatternRecognition

repo: OpenBMI


text:
function V = hogDraw( H, w )
% Create visualization of hog descriptor.
%
% USAGE
%  V = hogDraw( H, [w] )
%
% INPUTS
%  H          - [m n oBin*4] computed hog features
%  w          - [15] width for each glyph
%
% OUTPUTS
%  V          - [m*w n*w] visualization of hog features
%
% EXAMPLE
%
% See also hog
%
% Piotr's Image&Video Toolbox      Version 2.41
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% fold 4 normalizations
nFold=4; s=size(H); s(3)=s(3)/nFold; w0=H; H=zeros(s);
for o=0:nFold-1, H=H+w0(:,:,(1:s(3))+o*s(3)); end;

% construct a "glyph" for each orientaion
if(nargin<2 || isempty(w)), w=15; end
bar=zeros(w,w); bar(:,round(.45*w):round(.55*w))=1;
bars=zeros([size(bar) s(3)]);
for o=1:s(3), bars(:,:,o)=imrotate(bar,-(o-1)*180/s(3),'crop'); end

% make pictures of positive weights by adding up weighted glyphs
H(H<0)=0; V=zeros(w*s(1:2));
for r=1:s(1), rs=(1:w)+(r-1)*w;
  for c=1:s(2), cs=(1:w)+(c-1)*w;
    for o=1:s(3), V(rs,cs)=V(rs,cs)+bars(:,:,o)*H(r,c,o); end
  end
end

author: zouchuhang

repo: LayoutNet


text:
function    g = lfmGradientSigmaH3(gamma1, gamma2, sigma2, t1, t2, preFactor, mode, term)

% LFMGRADIENTSIGMAH3 Gradient of the function h_i(z) with respect \sigma.
% FORMAT
% DESC Computes the gradient of the function h_i(z) with respect to the
% length-scale of the input "force", \sigma.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode: indicates in which way the vectors t1 and t2 must be transposed
% RETURN g : Gradient of the function with respect to \sigma.
%
% COPYRIGHT : David Luengo, 2007, 2008, Mauricio Alvarez, 2008
%
% MODIFICATIONS : Mauricio Alvarez, 2008

% SEEALSO : lfmKernGradient, lfmXlfmKernGradient, lfmGradientSigmaUpsilon

% KERN

% Gradient

if nargin<8
    term = [];
end

if ~mode
    if ~term
        g = preFactor*lfmGradientSigmaUpsilonMatrix(gamma1,sigma2, t1,t2);
    else
        gradupsilon = lfmGradientSigmaUpsilonMatrix(gamma1,sigma2, t1,t2);
        g = -preFactor(1)*gradupsilon + preFactor(2)*conj(gradupsilon);
    end
else
    g =  preFactor(1)*lfmGradientSigmaUpsilonMatrix(gamma1,sigma2,t1,t2) + ...
        preFactor(2)*lfmGradientSigmaUpsilonMatrix(gamma2,sigma2,t1,t2);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function optval = PW_Delta_CallPut(S,K,C,r,T)
% S = NSim x 1 matrix of simulated prices
% K = Strike price
% C = 1 -> Call; C = 0 -> Put

    if(C==1)
        Indicator = (S(:,end)>K);
    else
        Indicator = (S(:,end)<=K);
    end
    optval = mean(S(:,end)./S(1,1).*Indicator)*exp(-r*T);
end

author: Sable

repo: mcbench-benchmarks


text:
function varargout = process_snapshot( varargin )
% PROCESS_SNAPSHOT: Save snapshot.
%
% USAGE:     sProcess = process_snapshot('GetDescription')
%                       process_snapshot('Run', sProcess, sInputs)

% @=============================================================================
% 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, 2012-2022

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Save snapshot';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'File';
    sProcess.Index       = 982;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Scripting';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw', 'data', 'results', 'timefreq', 'matrix', 'dipoles', 'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.OutputTypes = {'raw', 'data', 'results', 'timefreq', 'matrix', 'dipoles', 'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Definition of the options
    % === TARGET
    sProcess.options.type.Comment = 'Snapshot: ';
    sProcess.options.type.Type    = 'combobox_label';
    sProcess.options.type.Value   = {1, {...
        'Sensors/MRI registration',              'registration'; ...    % 1
        'SSP projectors',                        'ssp'; ...             % 2
        'Noise covariance',                      'noiscov'; ...         % 3
        'Data covariance',                       'ndatacov'; ...        % 12
        'Headmodel spheres',                     'headmodel'; ...       % 4
        'Recordings time series',                'data'; ...            % 5
        'Recordings topography (one time)',      'topo'; ...            % 6
        'Recordings topography (contact sheet)', 'topo_contact'; ...    % 7
        'Sources (one time)',                    'sources'; ...         % 8
        'Sources (contact sheet)',               'sources_contact'; ... % 9
        'Sources (MRI viewer)',                  'mriviewer'; ...
        'Frequency spectrum',                    'spectrum'; ...        % 10
        'Time-frequency maps',                   'timefreq'; ...        % 14
        'Connectivity matrix',                   'connectimage'; ...    % 11
        'Connectivity graph',                    'connectgraph'; ...
        'Dipoles',                               'dipoles'; ...         % 13
        }'};
    % === SENSORS
    sProcess.options.modality.Comment    = 'Sensor type: ';
    sProcess.options.modality.Type       = 'combobox';
    sProcess.options.modality.Value      = {1, {'MEG (All)', 'MEG (Gradiometers)', 'MEG (Magnetometers)', 'EEG', 'ECOG', 'SEEG', 'NIRS', 'EOG', 'ECG', 'EMG'}};
    sProcess.options.modality.InputTypes = {'raw', 'data', 'timefreq', 'pdata', 'ptimefreq'};
    % === Orientation
    sProcess.options.orient.Comment    = 'Orientation: ';
    sProcess.options.orient.Type       = 'combobox';
    sProcess.options.orient.Value      = {1, {'left', 'right', 'top', 'bottom', 'front', 'back', 'left_intern', 'right_intern'}};
    sProcess.options.orient.InputTypes = {'raw', 'data', 'results', 'timefreq', 'dipoles', 'pdata', 'presults', 'ptimefreq'};
    % === TIME: Single view
    sProcess.options.time.Comment = 'Time (in seconds):';
    sProcess.options.time.Type    = 'value';
    sProcess.options.time.Value   = {0, 's', 4};
    % === TIME: Contact sheet
    sProcess.options.contact_time.Comment = 'Contact sheet (start time, stop time):';
    sProcess.options.contact_time.Type    = 'value';
    sProcess.options.contact_time.Value   = {[0,.1], 'list', 4};
    sProcess.options.contact_nimage.Comment = 'Contact sheet (number of images):';
    sProcess.options.contact_nimage.Type    = 'value';
    sProcess.options.contact_nimage.Value   = {12, '', 0};
    % === THRESHOLD
    sProcess.options.threshold.Comment = 'Amplitude threshold:';
    sProcess.options.threshold.Type    = 'value';
    sProcess.options.threshold.Value   = {30, '%', 0};
    sProcess.options.threshold.InputTypes = { 'results', 'timefreq', 'presults', 'ptimefreq'};
    % === SMOOTHING
    sProcess.options.surfsmooth.Comment = 'Surface smoothing:';
    sProcess.options.surfsmooth.Type    = 'value';
    sProcess.options.surfsmooth.Value   = {30, '%', 0};
    sProcess.options.surfsmooth.InputTypes = {'results', 'timefreq', 'presults', 'ptimefreq'};
    % === FREQUENCY
    sProcess.options.freq.Comment = 'Frequency:';
    sProcess.options.freq.Type    = 'value';
    sProcess.options.freq.Value   = {0, 'Hz', 2};
    sProcess.options.freq.InputTypes = {'timefreq', 'ptimefreq'};
    % === ROW NAMES
    sProcess.options.rowname.Comment    = 'Signal name (empty=all): ';
    sProcess.options.rowname.Type       = 'text';
    sProcess.options.rowname.Value      = '';
    sProcess.options.rowname.InputTypes = {'raw', 'data', 'timefreq', 'matrix', 'pdata', 'ptimefreq', 'pmatrix'};
    % === MNI coordinates
    sProcess.options.mni.Comment = 'MNI coordinates:';
    sProcess.options.mni.Type    = 'value';
    sProcess.options.mni.Value   = {[0,0,0], 'list', 3};
    sProcess.options.mni.InputTypes = {'results', 'timefreq', 'presults', 'ptimefreq'};
    % === COMMENT
    sProcess.options.Comment.Comment = 'Comment: ';
    sProcess.options.Comment.Type    = 'text';
    sProcess.options.Comment.Value   = '';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    iType = find(strcmpi(sProcess.options.type.Value{1}, sProcess.options.type.Value{2}(2,:)));
    Comment = ['Snapshot: ' sProcess.options.type.Value{2}{1,iType}];
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Returned files: same as input
    OutputFiles = {sInputs.FileName};
    % Get options
    if isfield(sProcess.options, 'target') && isfield(sProcess.options.target, 'Value') && ~isempty(sProcess.options.target.Value)
        switch (sProcess.options.target.Value)
            case 1,  SnapType = 'registration';
            case 2,  SnapType = 'ssp';
            case 3,  SnapType = 'noiscov';
            case 4,  SnapType = 'spheres';
            case 5,  SnapType = 'data';
            case 6,  SnapType = 'topo';
            case 7,  SnapType = 'topo_contact';
            case 8,  SnapType = 'sources';
            case 9,  SnapType = 'sources_contact';
            case 10, SnapType = 'spectrum';
            case 11, SnapType = 'connectimage';
            case 12, SnapType = 'ndatacov';
            case 13, SnapType = 'dipoles';
            case 14, SnapType = 'timefreq';
        end
    else
        SnapType = sProcess.options.type.Value{1};
    end
    if isfield(sProcess.options, 'modality') && isfield(sProcess.options.modality, 'Value') && ~isempty(sProcess.options.modality.Value) && iscell(sProcess.options.modality.Value)
        switch (sProcess.options.modality.Value{1})
            case 1,  Modality = 'MEG';
            case 2,  Modality = 'MEG GRAD';
            case 3,  Modality = 'MEG MAG';
            case 4,  Modality = 'EEG';
            case 5,  Modality = 'ECOG';
            case 6,  Modality = 'SEEG';
            case 7,  Modality = 'NIRS';
            case 8,  Modality = 'EOG';
            case 9,  Modality = 'ECG';
            case 10, Modality = 'EMG';
            otherwise, Modality = [];
        end
    else
        Modality = [];
    end
    Time           = sProcess.options.time.Value{1};
    contact_time   = sProcess.options.contact_time.Value{1};
    contact_nimage = sProcess.options.contact_nimage.Value{1};
    if isfield(sProcess.options, 'orient') && isfield(sProcess.options.orient, 'Value') && ~isempty(sProcess.options.orient.Value) && iscell(sProcess.options.orient.Value)
        Orient = sProcess.options.orient.Value{2}{sProcess.options.orient.Value{1}};
    else
        Orient = [];
    end
    if isfield(sProcess.options, 'freq') && isfield(sProcess.options.freq, 'Value') && ~isempty(sProcess.options.freq.Value) && iscell(sProcess.options.freq.Value) && ~isequal(sProcess.options.freq.Value{1}, 0)
        Freq = sProcess.options.freq.Value{1};
    else
        Freq = [];
    end
    if isfield(sProcess.options, 'mni') && isfield(sProcess.options.mni, 'Value') && iscell(sProcess.options.mni.Value) && (length(sProcess.options.mni.Value{1}) == 3)
        XYZmni = sProcess.options.mni.Value{1};
    else
        XYZmni = [];
    end
    % If using "comment" instead of "Comment" (common scripting error)
    Comment = sProcess.options.Comment.Value;
    if isempty(Comment) && isfield(sProcess.options, 'comment') && isfield(sProcess.options.comment, 'Value') && ~isempty(sProcess.options.comment.Value)
        Comment = sProcess.options.comment.Value;
    end
    % Amplitude threshold
    if isfield(sProcess.options, 'threshold') && isfield(sProcess.options.threshold, 'Value') && ~isempty(sProcess.options.threshold.Value) && iscell(sProcess.options.threshold.Value)
        Threshold = sProcess.options.threshold.Value{1} / 100;
    else
        Threshold = 0.3;
    end
    % Surface smoothing
    if isfield(sProcess.options, 'surfsmooth') && isfield(sProcess.options.surfsmooth, 'Value') && ~isempty(sProcess.options.surfsmooth.Value) && iscell(sProcess.options.surfsmooth.Value)
        SurfSmooth = sProcess.options.surfsmooth.Value{1} / 100;
    else
        SurfSmooth = 0.3;
    end
    % Row names
    if isfield(sProcess.options, 'rowname') && isfield(sProcess.options.rowname, 'Value') && ~isempty(sProcess.options.rowname.Value)
        RowName = strtrim(sProcess.options.rowname.Value);
        iComa = find(RowName == ',');
        if ~isempty(iComa)
            bst_report('warning', sProcess, [], 'The option "signal name" should be empty or include only one signal.');
            RowName = strtrim(RowName(1,iComa-1));
        end
    else
        RowName = [];
    end
    % Contact sheet
    Contact = [contact_time, contact_nimage];
    % Select only one input file per channel file
    if ismember(SnapType, {'registration', 'ssp', 'noiscov', 'ndatacov', 'headmodel'})
        [AllChannelFile, iAllChan] = unique({sInputs.ChannelFile});
        sInputs = sInputs(iAllChan);
    end
    % For each file, capture view for each file
    for iFile = 1:length(sInputs)
        FileName = sInputs(iFile).FileName;
        switch (SnapType)
            case 'registration'
                bst_report('Snapshot', 'registration', FileName, Comment, Modality, Orient);
            case 'ssp'
                bst_report('Snapshot', 'ssp', FileName, Comment);
            case 'noiscov'
                bst_report('Snapshot', 'noisecov', FileName, Comment);
            case 'ndatacov'
                bst_report('Snapshot', 'ndatacov', FileName, Comment);
            case 'headmodel'
                bst_report('Snapshot', 'headmodel', FileName, Comment);
            case 'data'
                bst_report('Snapshot', 'data', FileName, Comment, Modality, Time, RowName);
            case 'topo'
                bst_report('Snapshot', 'topo', FileName, Comment, Modality, Time, Freq);
            case 'topo_contact'
                if (length(Contact) ~= 3)
                    bst_report('Error', sProcess, [], 'Invalid contact sheet time values');
                    return;
                end
                bst_report('Snapshot', 'topo', FileName, Comment, Modality, Contact, Freq);
            case 'sources'
                bst_report('Snapshot', 'sources', FileName, Comment, Time, Threshold, Orient, SurfSmooth, Freq);
            case 'sources_contact'
                if (length(Contact) ~= 3)
                    bst_report('Error', sProcess, [], 'Invalid contact sheet time values');
                    return;
                end
                bst_report('Snapshot', 'sources', FileName, Comment, Contact, Threshold, Orient, SurfSmooth, Freq);
            case 'mriviewer'
                bst_report('Snapshot', 'mriviewer', FileName, Comment, Time, Threshold, Freq, XYZmni);
            case 'spectrum'
                bst_report('Snapshot', 'spectrum', FileName, Comment, RowName, Freq);
            case 'timefreq'
                bst_report('Snapshot', 'timefreq', FileName, Comment, RowName, Time, Freq);
            case 'connectimage'
                bst_report('Snapshot', 'connectimage', FileName, Comment, Time, Freq);
            case 'connectgraph'
                bst_report('Snapshot', 'connectgraph', FileName, Comment, Threshold, Time, Freq);
            case 'dipoles'
                bst_report('Snapshot', 'dipoles', FileName, Comment, Threshold, Orient);
        end
    end
end

author: brainstorm-tools


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

function newImages = transImagesFwd(oldImages,T,rows,cols)
%newImages = transImagesFwd(oldImages)
%   images are all cell array.

m = length(oldImages);
if m~=length(T),error('size(oldImages,1) must equal to length(T)!');end;

if (m<100) || (matlabpool('size')==0)
    newImages = arrayfun(@(i)imtransform(oldImages{i},T(i),'XData',...
        [1 cols],'YData',[1 rows]),1:m,'UniformOutput',false);
    if (m>=100)
        warning('Please launch matlabpool to speed up your program!');
    end
else
    newImages = cell(m,1);
    parfor i = 1:m
        newImages{i} = imtransform(oldImages{i},T(i),'XData',...
            [1 cols],'YData',[1 rows]);
    end
end
for i = 1:m, if size(newImages{i},1)~=rows || size(newImages{i},2)~=cols, newImages{i} = imresize(newImages{i},[rows,cols]);end;end;

end

author: zhusz


text:
function cc=lpczz2cc(zz,np)
%LPCZZ2CC Convert poles to "complex" cepstrum CC=(ZZ,NP)

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

[nf,p]=size(zz);
if (nargin < 2) np=p; end
cc=zeros(nf,np);
yy=zz.';
if p<2			% special case 'cos sum() is weird
  cc(:,1)=real(zz);
  for k=2:np
    yy=yy.*zz.';
    cc(:,k)=real(yy).'/k;
  end
else
  cc(:,1)=sum(real(yy)).';
  for k=2:np
    yy=yy.*zz.';
    cc(:,k)=sum(real(yy)).'/k;
  end
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function calcRigor( config )
%RIGOR Summary of this function goes here
%   Detailed explanation goes here

rigorconfig = config.rigor;
param = rigorconfig.params.param_set;

%Check if image location exists or not.

if(~exist(config.imageLocation, 'dir'))
	fprintf('Image Location does not exist. Please check path once again \n');
	return;
end

if(~exist(config.outputLocation, 'dir'))
	fprintf('Image Location does not exist. Please check path once again \n');
	return;
end

%Load All images in a particular folder
images = dir(config.imageLocation);
images = regexpi({images.name}, '.*jpg|.*jpeg|.*png|.*bmp', 'match');
images = [images{:}];

for i=1:length(images)
    imname = char(images(i));
    impath = fullfile(config.imageLocation, imname);
    whos impath
    im=imread(impath);
    fprintf('Calculating RIGOR for %s\n', imname);
    proposals =calcrigorForIm(impath,rigorconfig);
    saveFile=[imname '.mat'];
    save([config.outputLocation saveFile], 'proposals');
end
end

author: Cloud-CV

repo: object-proposals


text:
% Least Mean Squares algorithm
%
% From A. H. Sayed, "Fundamentals of adaptive filtering}", Wiley-IEEE
% Press, 2003, Chapter 5.
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef lms < linear_filter

    properties (GetAccess = 'public', SetAccess = 'private')
        mu = 0.001; % learning rate
    end

    properties (GetAccess = 'public', SetAccess = 'private')
        w = []; % filter coefficients
    end

    methods
        function obj = lms(parameters) % constructor
            if (nargin > 0) % copy valid parameters
                for fn = fieldnames(parameters)'
                    if ismember(fn,fieldnames(obj))
                        obj.(fn{1}) = parameters.(fn{1});
                    end
                end
            end
        end

        function y_est = evaluate(obj,x) % evaluate the algorithm
            if numel(obj.w)>0
                y_est = x*obj.w;
            else
                y_est = zeros(size(x,1),1);
            end
        end

        function train(obj,x,y) % train the algorithm
            if numel(obj.w)==0 % initialize
                obj.w = zeros(length(x),1);
            end

            % Algorithm 5.2.1 in reference
            err = y - x*obj.w; % instantaneous error
            obj.w = obj.w + obj.mu*x'*err; % update filter coefficients
        end

    end
end

repo: kafbox


text:
function [dat_sf,W] = proc_multiBandSpatialFilter(dat,method)
%PROC_MULTIBANDSPATIALFILTER - Apply spatial filtering method to multiple
%frequency bands
%
%Synopsis:
% [DAT_SF, W] = proc_multiBandSpatialFilter(DAT, METHOD);
%
%Arguments:
% DAT    - data structure of epoched and pre-filtered data
% METHOD - cell array containing in the first entry a function handle to
%          the spatial filtering method and (optionally) in the following
%          entries parameters for that function
%
%Returns:
% DAT_SF - updated data structure, containing the spatially filtered data
%          of all bands, appended as channels
% W      - cell array with each entry containing the spatial filters (in
%          the columns) of each band
%
%Description:
% The input data structure is assumed to contain pre-filtered channels as
% returned by the function proc_filterbank. proc_multiBandSpatialFilter can
% be used either alone or as a processing step in the crossvalidation
% function.
%
%See also: crossvalidation proc_filterbank proc_multiBandLinearDerivation

% 10-2015: schultze-kraft@tu-berlin.de

misc_checkType(dat,'STRUCT(x clab y)');
misc_checkType(method,'CELL');
procFunc = method{1};
misc_checkType(procFunc,'!FUNC');
if length(method)>1
    procPar = method(2:end);
else
    procPar = {};
end

% get number of frequency bands
band_ix = zeros(1,length(dat.clab));
flt_ix = cellfun(@(x) strfind(x,'flt'),dat.clab);
for ii = 1:length(dat.clab)
    band_ix(ii) = str2double(dat.clab{ii}(flt_ix(ii)+3:end));
end
n_bands = max(band_ix);

% band-wise apply spatial filtering method
W = cell(1,n_bands);
dat_sf = [];
for bi = 1:n_bands
    dat2 = proc_selectChannels(dat,sprintf('*flt%d',bi));
    [dat2,W{bi}] = procFunc(dat2,procPar{:});
    dat_sf = proc_appendChannels(dat_sf,dat2);
end

author: bbci


function [res] = plotBeamformer(M, f, d, bfFilter, theta)
    c = 340;
    w = bfFilter;       % beamformer
    p = zeros(length(theta), length(f));
    for  j=1:length(theta)                    %scan angles
        a=exp(-1j*[0:M-1]'*2*pi*f*cos(theta(j))*d/c);
        p(j, :) = sum(w.*a);
    end
   % p=p';
    res = p;

    if length(f) == 1
        subplot(1, 2, 1);
        plot(theta/pi*180, abs(p.')); grid on;
        xlabel('Degree')
        subplot(1, 2, 2);
        polar(theta,abs(p.'))
    else
        subplot(1, 2, 1);
        mesh(f, theta/pi*180, abs(p))
        xlabel('Frequency(Hz)')
        ylabel('Degree')
        title('Beamformer Directivity')
        subplot(1, 2, 2);
        imagesc(abs(p).');
    end

end


text:
%CONFMAT Construct confusion matrix
%
%  [C,NE,LABLIST1,LABLIST2] = CONFMAT(LAB1,LAB2,METHOD,FID)
%
% INPUT
%  LAB1        Set of labels
%  LAB2        Set of labels
%  METHOD      'count' (default) to count number of co-occurences in
%	             LAB1 and LAB2, 'disagreement' to count relative
%		           non-co-occurrence.
%  FID         Write text result to file
%
% OUTPUT
%  C           Confusion matrix
%  NE          Total number of errors (empty labels are neglected)
%  LABLIST1    Label list for LAB1
%  LABLIST2    Label list for LAB2
%
% DESCRIPTION
% Constructs a confusion matrix C between two sets of labels LAB1
% (corresponding to the rows in C) and LAB2 (the columns in C). The order of
% the rows and columns is returned in LABLIST. NE is the total number of
% errors (sum of non-diagonal elements in C).
%
% If METHOD = 'count' (default), co-occurences in LAB1 and LAB2 are counted
% and returned in C.
% For METHOD = 'disagreement', the relative disagreement is returned in NE,
% and is split over all combinations of labels in C (such that the rows sum
% to 1). (The total disagreement for a class equals one minus the
% sensitivity for that class as computed by TESTC).
% For METHOD = 'ids' a cell array C is returned in which C(i,j) contains
% the indices of the objects for which LAB1 equals LABLIST1(i,:) and LAB2
% equals LABLIST2(j,:).
%
%   [C,NE,LABLIST] = CONFMAT(D,METHOD)
%
% If D is a classification result D = A*W, the labels LAB1 and LAB2 are
% internally retrieved by CONFMAT before computing the confusion matrix.
%
%   C = CONFMAT(D)
%
% This call also applies in case in D = A*W the dataset A has soft labels
% W is trained by a soft labeld classifier.
%
% When no output argument is specified, or when FID is given, the
% confusion matrix is displayed or written a to a text file. It is assumed
% that LAB1 contains true labels and LAB2 stores estimated labels.
%
% EXAMPLE
% Typical use of CONFMAT is the comparison of true and and estimated labels
% of a testset A by application to a trained classifier W:
% LAB1 = GETLABELS(A); LAB2 = A*W*LABELD.
% More examples can be found in PREX_CONFMAT, PREX_MATCHLAB.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, GETLABELS, LABELD

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: confmat.m,v 1.8 2010/06/01 08:45:31 duin Exp $

function [CC,ne,lablist1,lablist2] = confmat (arg1,arg2,arg3,fid)

		% Check arguments.
  if nargin < 4, fid = 1; end
	if nargin < 3 | isempty(arg3)
		if isdataset(arg1)
			if islabtype(arg1,'crisp')
				lablist1 = getlablist(arg1);
				nlab1 = getnlab(arg1);
				lab2 = arg1*labeld;
				nlab2 = renumlab(lab2,lablist1); % try to fit on original lablist
				if any(nlab2==0)                 % doesn't fit:  contruct its own one
					[nlab2,lablist2] = renumlab(lab2);
				else
					lablist2 = lablist1;
				end
				if nargin < 2| isempty(arg2)
					method = 'count';
					prwarning(4,'no method supplied, assuming count');
				else
					method = arg2;
				end
			elseif islabtype(arg1,'soft')
				true_labs = gettargets(arg1);
				if any(true_labs > 1 | true_labs < 0)
					error('True soft labels should be between 0 and 1')
				end
				est_labs = +arg1;
				if any(est_labs > 1 | est_labs < 0)
					error('Estimated soft labels should be between 0 and 1')
				end
				C = true_labs'*est_labs;
				output(C,getlablist(arg1),getfeatlab(arg1),'soft',fid);
				return
			else
				error('Confusion matrix can only be computed for crisp and soft labeld datasets')
			end
		else
			method = 'count';
			prwarning(4,'no method supplied, assuming count');
			if (nargin < 2 | isempty(arg2))
				error('Second label list not supplied')
			end
			[nlab1,lablist1] = renumlab(arg1);
			[nlab2,lablist2] = renumlab(arg2);
		end
	else
		[nlab1,lablist1] = renumlab(arg1);
		[nlab2,lablist2] = renumlab(arg2);
		%lab1 = arg1;
		%lab2 = arg2;
		method = arg3;
	end
	if nargin < 2
		if ~isdataset(arg1)
			error('two labellists or one dataset should be supplied')
		end
	end

	% Renumber LAB1 and LAB2 and find number of unique labels.

	m = size(nlab1,1);
	if (m~=size(nlab2,1))
		error('LAB1 and LAB2 have to have the same lengths.');
	end
	n1 = size(lablist1,1);
	n2 = size(lablist2,1);
	n = max(n1,n2);

	% Construct matrix of co-occurences (confusion matrix).

  if strcmp(method,'ids')
    C = cell(n1+1,n2+1); % we need to store object id's
  else
    C = zeros(n1+1,n2+1);
  end
	for i = 0:n1
		K = find(nlab1==i);
		if (isempty(K))
      if strcmp(method,'ids')
        C(i+1,:) = repmat({[]},1,n2+1);
      else
        C(i+1,:) = zeros(1,n2+1);
      end
		else
			for j = 0:n2
        if strcmp(method,'ids')
          C(i+1,j+1) = {K(nlab2(K)==j)};
        else
          C(i+1,j+1) = sum(nlab2(K)==j);
        end
			end
		end
	end

	% position rejects and unlabeled object at the end of the matrix

	D = C;
	D(1:end-1,1:end-1) = C(2:end,2:end);
	D(end,:) = [C(1,2:end) C(1,1)];
	D(1:end-1,end) = C(2:end,1);
	C = D;
	D = D(1:end-1,1:end-1);
	DD = D(1:min(n1,n2),1:min(n1,n2));

	% Calculate number of errors ('count') or disagreement ('disagreement').
	% Neglect rejects

	switch (method)
		case {'count','ids'}
			J = find(nlab1~=0 & nlab2~=0);
			ne = nlabcmp(lablist1(nlab1(J),:),lablist2(nlab2(J),:));
		case 'disagreement'
			ne = (sum(sum(D)) - sum(diag(DD)))/m;  % Relative sum of off-diagonal
			ne = ne/m;														 % entries.
			E = repmat(sum(D,2),1,n2);             % Disagreement = 1 -
			D = ones(n1,n2)-D./E;                  % relative co-occurence.
			D = D / (n-1);
		otherwise
			error('unknown method');
  end

  if ~strcmp(method,'ids')
    %Distinguish 'rejects / no_labels' from 'non_rejects / fully labeled'
    if (any(C(:,end) ~= 0) | any(C(end,:)~=0)) & strcmp(method,'count')
      n1 = n1+1; n2 = n2+1;
      labch = char(strlab(lablist2),'reject');
      labcv = char(strlab(lablist1),'No');
    else
      labch = strlab(lablist2);
      labcv = strlab(lablist1);
      %labcv = labch;
      C = D;
    end
  end

    % If no output argument is specified, pretty-print C.

	if ((nargout == 0) || nargin == 4) && (~strcmp(method,'ids'))

		output(C,labcv,labch,method,fid)
	end

	if nargout > 0 || strcmp(method,'ids')
		CC = C;
	end

return

function output(C,labcv,labch,method,fid)

	n1 = size(labcv,1);
	n2 = size(labch,1);

		% Make sure labels are stored in LABC as matrix of characters,
    % max. 6 per label.

		if (size(labch,2) > 6)
			labch = labch(:,1:6);
			%labcv = labcv(:,1:6);
		end
		if (size(labch,2) < 5)
			labch = [labch repmat(' ',n2,ceil((5-size(labch,2))/2))];
			labcv = [labcv repmat(' ',n1,ceil((5-size(labcv,2))/2))];
		end

		%C = round(1000*C./repmat(sum(C,2),1,size(C,2)));

		nspace = max(size(labcv,2)-7,0);
		cspace = repmat(' ',1,nspace);
		%fprintf(fid,['\n' cspace '        | Estimated Labels']);
		fprintf(fid,['\n  True   ' cspace '| Estimated Labels']);
		fprintf(fid,['\n  Labels ' cspace '|']);
		for j = 1:n2, fprintf(fid,'%7s',labch(j,:)); end
		fprintf(fid,'|');
		fprintf(fid,' Totals');
		fprintf(fid,'\n ');
		fprintf(fid,repmat('-',1,8+nspace));
		fprintf(fid,'|%s',repmat('-',1,7*n2));
		fprintf(fid,'|-------');
		fprintf(fid,'\n ');

		for j = 1:n1
			fprintf(fid,' %-7s|',labcv(j,:));
			switch (method)
				case 'count'
					fprintf(fid,'%5i  ',C(j,:)');
					fprintf(fid,'|');
					fprintf(fid,'%5i',sum(C(j,:)));
				case 'disagreement'
					fprintf(fid,' %5.3f ',C(j,:)');
					fprintf(fid,'|');
					fprintf(fid,' %5.3f ',sum(C(j,:)));
				case 'soft'
					fprintf(fid,' %5.2f ',C(j,:)');
					fprintf(fid,'|');
					fprintf(fid,' %5.2f ',sum(C(j,:)));
			end
			fprintf(fid,'\n ');
		end

		fprintf(fid,repmat('-',1,8+nspace));
		fprintf(fid,'|%s',repmat('-',1,7*n2));
		fprintf(fid,'|-------');
		fprintf(fid,['\n  Totals ' cspace '|']);

		switch (method)
			case 'count'
				fprintf(fid,'%5i  ',sum(C));
				fprintf(fid,'|');
				fprintf(fid,'%5i',sum(C(:)));
			case 'disagreement'
				fprintf(fid,' %5.3f ',sum(C));
				fprintf(fid,'|');
			%	fprintf(fid,' %5.3f ',sum(C(:)));
			case 'soft'
				fprintf(fid,' %5.2f ',C(j,:)');
				fprintf(fid,'|');
				fprintf(fid,' %5.2f ',sum(C(:)));
		end
		fprintf(fid,'\n\n');
	return

author: marianux

repo: ecg-kit

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


text:
% Here the training data is adapted from UCI ML repository, 'housing' data
% Input variables: 12 continous, one binary
% Ouput variables: continous
% The testing result trace is in the end of this script, it is same to the graph in page 219 of
% Leo Brieman etc. 1984 book titled "Classification and regression trees".

dtreeCPD=tree_CPD;

% load data
fname = fullfile(BNT_HOME, 'examples', 'static', 'uci_data', 'housing', 'housing.data');
data=load(fname);
data=data';
data=transform_data_into_bnt_format(data,[1:3,5:14]);

% learn decision tree from data
ns=1*ones(1,14);
ns(4)=2;
dtreeCPD1=learn_params(dtreeCPD,1:14,data,ns,[1:3,5:14],'stop_cases',5,'min_gain',0.006);

% evaluate on data
[score,outputs]=evaluate_tree_performance(dtreeCPD1,1:14,data,ns,[1:3,5:14]);
fprintf('Mean square deviation (using regression tree to predict) in old training data %6.3f\n',score);

% show decision tree using graphpad
% It should be easy, but still not implemented

% >> test_housing
% Create node 1 split at 6 gain  38.2205 Th 6.939000e+000. Mean  22.5328 Cases 506
% Create node 2 split at 13 gain  14.4503 Th 1.437000e+001. Mean  19.9337 Cases 430
% Create node 3 split at 8 gain   4.9809 Th 1.358000e+000. Mean  23.3498 Cases 255
% Create node 4 split at 1 gain   0.7722 Th 1.023300e+001. Mean  45.5800 Cases 5
% Create leaf node(samevalue) 5. Mean  50.0000 Std   0.0000 Cases 4
% Add subtree node 5 to 4. #nodes 5
% Create leaf node(samevalue) 6. Mean  27.9000 Std   0.0000 Cases 1
% Add subtree node 6 to 4. #nodes 6
% Add subtree node 4 to 3. #nodes 6
% Create node 7 split at 6 gain   2.8497 Th 6.540000e+000. Mean  22.9052 Cases 250
% Create node 8 split at 13 gain   0.5970 Th 7.560000e+000. Mean  21.6297 Cases 195
% Create leaf node(nogain) 9. Mean  23.9698 Std   1.7568 Cases 43
% Add subtree node 9 to 8. #nodes 9
% Create leaf node(nogain) 10. Mean  20.9678 Std   2.8242 Cases 152
% Add subtree node 10 to 8. #nodes 10
% Add subtree node 8 to 7. #nodes 10
% Create leaf node(nogain) 11. Mean  27.4273 Std   3.4512 Cases 55
% Add subtree node 11 to 7. #nodes 11
% Add subtree node 7 to 3. #nodes 11
% Add subtree node 3 to 2. #nodes 11
% Create node 12 split at 1 gain   2.2467 Th 6.962150e+000. Mean  14.9560 Cases 175
% Create node 13 split at 5 gain   0.5172 Th 5.240000e-001. Mean  17.1376 Cases 101
% Create leaf node(nogain) 14. Mean  20.0208 Std   3.0672 Cases 24
% Add subtree node 14 to 13. #nodes 14
% Create leaf node(nogain) 15. Mean  16.2390 Std   2.9746 Cases 77
% Add subtree node 15 to 13. #nodes 15
% Add subtree node 13 to 12. #nodes 15
% Create node 16 split at 5 gain   0.6133 Th 6.050000e-001. Mean  11.9784 Cases 74
% Create leaf node(nogain) 17. Mean  16.6333 Std   4.5052 Cases 12
% Add subtree node 17 to 16. #nodes 17
% Create leaf node(nogain) 18. Mean  11.0774 Std   3.0090 Cases 62
% Add subtree node 18 to 16. #nodes 18
% Add subtree node 16 to 12. #nodes 18
% Add subtree node 12 to 2. #nodes 18
% Add subtree node 2 to 1. #nodes 18
% Create node 19 split at 6 gain   6.0493 Th 7.420000e+000. Mean  37.2382 Cases 76
% Create node 20 split at 1 gain   1.9900 Th 7.367110e+000. Mean  32.1130 Cases 46
% Create node 21 split at 8 gain   0.6273 Th 1.877300e+000. Mean  33.3488 Cases 43
% Create leaf node(samevalue) 22. Mean  45.6500 Std   6.1518 Cases 2
% Add subtree node 22 to 21. #nodes 22
% Create leaf node(nogain) 23. Mean  32.7488 Std   3.5690 Cases 41
% Add subtree node 23 to 21. #nodes 23
% Add subtree node 21 to 20. #nodes 23
% Create leaf node(samevalue) 24. Mean  14.4000 Std   3.7363 Cases 3
% Add subtree node 24 to 20. #nodes 24
% Add subtree node 20 to 19. #nodes 24
% Create node 25 split at 1 gain   1.1001 Th 2.733970e+000. Mean  45.0967 Cases 30
% Create leaf node(nogain) 26. Mean  45.8966 Std   4.4005 Cases 29
% Add subtree node 26 to 25. #nodes 26
% Create leaf node(samevalue) 27. Mean  21.9000 Std   0.0000 Cases 1
% Add subtree node 27 to 25. #nodes 27
% Add subtree node 25 to 19. #nodes 27
% Add subtree node 19 to 1. #nodes 27
% Mean square deviation (using regression tree to predict) in old training data  9.405
%

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function line_monte_carlo_test01 ( )

%*****************************************************************************80
%
%% LINE_MONTE_CARLO_TEST01 estimates integrals in 1D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_MONTE_CARLO_TEST01\n' );
  fprintf ( 1, '  Use LINE01_SAMPLE to estimate integrals\n' );
  fprintf ( 1, '  along the length of the unit line in 1D.\n' );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '         N' );
  fprintf ( 1, '        1' );
  fprintf ( 1, '               X' );
  fprintf ( 1, '              X^2' );
  fprintf ( 1, '             X^3' );
  fprintf ( 1, '             X^4' );
  fprintf ( 1, '             X^5' );
  fprintf ( 1, '           X^6\n' );
  fprintf ( 1, '\n' );

  n = 1;

  while ( n <= 65536 )

    [ x, seed ] = line01_sample ( n, seed );

    fprintf ( 1, '  %8d', n );

    for j = 1 : 7

      e = j - 1;

      value = monomial_value_1d ( n, e, x );

      result = line01_length ( ) * sum ( value(1:n) ) / n;

      fprintf ( 1, '  %14.6g', result );

    end

    fprintf ( 1, '\n' );

    n = 2 * n;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '     Exact' );

  for j = 1 : 7

    e = j - 1;

    result = line01_monomial_integral ( e );
    fprintf ( 1, '  %14.6g', result );

  end

  fprintf ( 1, '\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
%  Context-Aware Correlation Filters
%
%  Written by Joao F. Henriques, 2014
%  Revised by Yang Li, August, 2014
%  Adapted by Matthias Mueller, 2016
%
%  This function takes care of setting up parameters, loading video
%  information and computing precisions. For the actual tracking code,
%  check out the TRACKER function.
%
%  RUN_TRACKER
%    Without any parameters, will ask you to choose a video,
%    and show the results in an interactive figure.
%    Press 'Esc' to stop the tracker early. You can navigate the
%    video using the scrollbar at the bottom.
%
%  RUN_TRACKER VIDEO
%    Allows you to select a VIDEO by its name. 'all' will run all videos
%    and show average statistics. 'choose' will select one interactively.
%
%  RUN_TRACKER(VIDEO, SHOW_VISUALIZATION, SHOW_PLOTS)
%    Decide whether to show the scrollable figure, and the precision plot.

function [precision, fps] = run_tracker(video, show_visualization, show_plots)

	%path to the videos (you'll be able to choose one with the GUI).
	base_path = 'sequences/';

	%default settings
	if nargin < 1, video = 'Car1'; end
	if nargin < 2, show_visualization = ~strcmp(video, 'all'); end
	if nargin < 3, show_plots = ~strcmp(video, 'all'); end

    %default settings
    kernel.type =  'linear';

    padding = 2;  %extra area surrounding the target
    lambda1 = 1e-4;     %regularization
    lambda2 = 0.4;
    interp_factor = 0.005;  %linear interpolation factor for adaptation
    output_sigma_factor = 0.1; %spatial bandwidth (proportional to target)

    features.gray = false;
    features.hog = false;
    features.hogcolor = true;
    features.hog_orientations = 9;

    cell_size = 4;

	switch video
	case 'choose',
		%ask the user for the video, then call self with that video name.
		video = choose_video(base_path);
		if ~isempty(video),
			[precision, fps] = run_tracker(video, show_visualization, show_plots);

			if nargout == 0,  %don't output precision as an argument
				clear precision
			end
		end

	case 'all',
		%all videos, call self with each video name.

		%only keep valid directory names
		dirs = dir(base_path);
		videos = {dirs.name};
		videos(strcmp('.', videos) | strcmp('..', videos) | ...
			strcmp('anno', videos) | ~[dirs.isdir]) = [];

		%the 'Jogging' sequence has 2 targets, create one entry for each.
		%we could make this more general if multiple targets per video
		%becomes a common occurence.
		videos(strcmpi('Jogging', videos)) = [];
		videos(end+1:end+2) = {'Jogging.1', 'Jogging.2'};

		all_precisions = zeros(numel(videos),1);  %to compute averages
		all_fps = zeros(numel(videos),1);

		if ~exist('matlabpool', 'file'),
			%no parallel toolbox, use a simple 'for' to iterate
			for k = 1:numel(videos),
				[all_precisions(k), all_fps(k)] = run_tracker(videos{k}, show_visualization, show_plots);
			end
		else
			%evaluate trackers for all videos in parallel
			if matlabpool('size') == 0,
				matlabpool open;
			end
			parfor k = 1:numel(videos),
				[all_precisions(k), all_fps(k)] = run_tracker(videos{k}, show_visualization, show_plots);
			end
		end

		%compute average precision at 20px, and FPS
		mean_precision = mean(all_precisions);
		fps = mean(all_fps);
		fprintf('\nAverage precision (20px):% 1.3f, Average FPS:% 4.2f\n\n', mean_precision, fps)
		if nargout > 0,
			precision = mean_precision;
		end

	case 'benchmark',
		%running in benchmark mode - this is meant to interface easily
		%with the benchmark's code.

		%get information (image file names, initial position, etc) from
		%the benchmark's workspace variables
		seq = evalin('base', 'subS');
		target_sz = seq.init_rect(1,[4,3]);
		pos = seq.init_rect(1,[2,1]) + floor(target_sz/2);
		img_files = seq.s_frames;
		video_path = [];

		%call tracker function with all the relevant parameters
		[positions,rect_results,t]= tracker(video_path, img_files, pos, target_sz, ...
			padding, kernel, lambda1, lambda2, output_sigma_factor, interp_factor, ...
			cell_size, features, 0);

		%return results to benchmark, in a workspace variable
		rects =rect_results;
%         [positions(:,2) - target_sz(2)/2, positions(:,1) - target_sz(1)/2];
% 		rects(:,3) = target_sz(2);
% 		rects(:,4) = target_sz(1);
		res.type = 'rect';
		res.res = rects;
		assignin('base', 'res', res);

	otherwise
		%we were given the name of a single video to process.

		%get image file names, initial state, and ground truth for evaluation
		[img_files, pos, target_sz, ground_truth, video_path] = load_video_info(base_path, video);

		%call tracker function with all the relevant parameters
		[positions,~, time] = tracker(video_path, img_files, pos, target_sz, ...
			padding, kernel, lambda1, lambda2, output_sigma_factor, interp_factor, ...
			cell_size, features, show_visualization);

		%calculate and show precision plot, as well as frames-per-second
		precisions = precision_plot(positions, ground_truth, video, show_plots);
		fps = numel(img_files) / time;

		fprintf('%12s - Precision (20px):% 1.3f, FPS:% 4.2f\n', video, precisions(20), fps)

		if nargout > 0,
			%return precisions at a 20 pixels threshold
			precision = precisions(20);
		end

	end
end

repo: Context-Aware-CF-Tracking


text:
function dz = cartPoleDynamicsHumanReadable(z,u,p)
% dz = cartPoleDynamicsHumanReadable(z,u,p)
%
% This function computes the first-order dynamics of the cart-pole.
%
% INPUTS:
%   z = [4, n] = [x;q;dx;dq] = state of the system
%   u = [1, n] = horizontal force applied to the cart
%   p = parameter struct
%       .g = gravity
%       .m1 = cart mass
%       .m2 = pole mass
%       .l = pendulum length
% OUTPUTS:
%   dz = dz/dt = time derivative of state
%
%

% x = z(1,:);   %Cart position (Not used in dynamics)
q = z(2,:);   % pendulum (pole) angle, measure from gravity vector
dx = z(3,:);  % cart velocity
dq = z(4,:);  %pendulum angle rate

% Unpack the physical parameters
l = p.l;  %Pendulum length
m1 = p.m1; % cart mass
m2 = p.m2; % pole mass
g = p.g;  %Gravity acceleration

% The following expressions were derived using the symbolic math toolbox,
% in the script: Derive_cartPole.m
% Typically I use the automatically generated code, but I've written out
% the code in a more readable way here for tutorial purposes.

mEff = m1 + m2 - m2.*cos(q).^2;
centr = dq.^2.*l.*m2.*sin(q);
ddx = (u + centr + g.*m2.*cos(q).*sin(q))./mEff;
ddq = -(u.*cos(q) + (m1 + m2)*g*sin(q) + centr.*cos(q))./(l.*mEff);

dz = [dx;dq;ddx;ddq];

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%****************************************************
% Function to calculate sub-band energy(SBC) parameters
% from enframed signal
% modifided on 17 jan 2008
% *********sarikaya paper algorithm***********
% fs samapling rate 8000
%   frame size used 192
%s signal
%
%
%***************************************************
function feature= sbc_2(s,fs)

frames=enframe(s,hamming(192),192);
no_offrames=size(frames);
no_of_frames=no_offrames(:,1);
feature=zeros(no_of_frames,12);
energy=zeros(24,1);

fe_frame=0;  % This counter is used to avoid frames with no energy is a sub band

f=statusbar('Extracting SBC Feature');% Creates status bar

for i=1:no_of_frames
    energy=ones(24);

    f=statusbar((i/no_of_frames),f);% updates status bar
    [coef,len]=wavedec(frames(i,:),6,'db3');
    s_no=cumsum(len);

    %% next two gives the last node in wavelet pacet

    energy(1)=en(coef(1:s_no(1)));

    energy(2)=en(coef((s_no(1)+1):s_no(2)));
   % coef((s_no(1)+1):s_no(2))
            %refer diagram in sarikaya paper
            [coef1,len1]=wavedec(coef(s_no(2):s_no(3)),1,'db3');
            len1=cumsum(len1);
                energy(3)=en(coef1(1:len1(1)) );
                energy(4)=en(coef1((len1(1)+1):len1(2)));

                    [coef2,len2]=wavedec(coef((s_no(3)+1):s_no(4)),2,'db3');
                        len2=cumsum(len2);
                        energy(5)=en(coef2(1:len2(1)));

                        energy(6)=en(coef2( (len2(1)+1):len2(2) ) );
                            [coef3,len3]=wavedec(coef2((len2(2)+1):len2(3)),1,'db3');
                            len3=cumsum(len3);
                                    energy(7)=en(coef3(1:len3(1)));
                                    energy(8)=en(coef3((len3(1)+1):len3(2)));

                % all nodes at level six are over

           [coef4,len4]=wavedec(coef(s_no(4):s_no(5)),2,'db3');
            len4=cumsum(len4);
                        energy(9)=en(coef4(1:len4(1)));
                        energy(10)=en(coef4((len4(1)+1):len4(2)));

                        [coef5,len5]=wavedec(coef4((len4(2)+1):len4(3)),1,'db3');
                        len5=cumsum(len5);
                                    energy(11)=en(coef5(1:len5(1)));
                                    energy(12)=en(coef5((len5(1)+1):len5(2)));

            [coef7,len7]=wavedec(coef(s_no(5):s_no(6)),3,'db3');
            len7=cumsum(len7);
                       energy(13)=en(coef7(1:len7(1)));
                       energy(14)=en(coef7((len7(1)+1):len7(2)));

                        [coef8,len8]=wavedec(coef7((len7(2)+1):len7(3)),1,'db3');
                        len8=cumsum(len8);
                           energy(15)=en(coef8(1:len8(1)));
                           energy(16)=en(coef8((len8(1)+1):len8(2)));

             [coef9,len9]=wavedec(coef7((len7(3)+1):len7(4)),2,'db3');
              len9=cumsum(len9);
                        energy(17)=en(coef9(1:len9(1)));
                        energy(18)=en(coef9((len9(1)+1):len9(2)));

                 energy(19)=en(coef9((len9(2)+1):len9(3)));

  [coef10,len10]=wavedec(coef((s_no(6)+1):s_no(7)),3,'db3');
  len10=cumsum(len10);

           energy(20)=en(coef10(1:len10(1)));
           energy(21)=en(coef10((len10(1)+1):len10(2)));
           energy(22)=en(coef10((len10(2)+1):len10(3)));

     [coef11,len11]=wavedec(coef10((len(3)+1):len10(4)),1,'db3');
     len11=cumsum(len11);
            energy(23)=en(coef11(1:len11(1)));
            energy(24)=en(coef11((len11(1)+1):len11(2)));

% Taking 12 filter bank equalent
%
%f=rdct(feature); function not woring
% 19 jan DCT modified file in dessai cmtr
if all(energy>0)
  fe_frame=fe_frame+1;
    log_en=log(energy.*1E+06);

   for j=1:12
   for k=1:24

        feature(fe_frame,j)=feature(fe_frame,j)+log_en(k)*cos((j*(k-0.5)*pi)/24);
   end
   end
 % feature=abs(feature); this needs to be checked
end
end
 delete(statusbar)


text:
classdef (HandleCompatible) AndNode < Node
    % AndNode are an class that represents AND connections in a logical formula
    % For further documentation please have a look at the Node Class.
    % .. Authors
    %     - Thomas Pfau 2016
    properties
    end

    methods
        function res = evaluate(self,assignment, printLevel)
            if ~exist('printLevel','var')
                printLevel = 0;
            end
            res = true;
            for i=1:numel(self.children)
                child = self.children(i);
                if not(child.evaluate(assignment,printLevel))
                    res = false;
                end
            end
            if printLevel >= 1
                fprintf('%s : %i\n',self.toString(0),res);
            end
        end

        function tf = deleteLiteral(self,literalID, keepClauses)
            tf = true;
            if ~exist('keepClauses','var')
                keepClauses = true;
            end
            if ~keepClauses
                % reduce to properly delete elements
                self.reduce();
            end
            arrayfun(@(x) ~isa(x,'LiteralNode') && x.deleteLiteral(literalID, keepClauses), self.children);
            % originalNodeString = self.toString(1);
            literalMatches = arrayfun(@(x) (isa(x, 'LiteralNode') && x.contains(literalID) ), self.children);
            emptyChildren = arrayfun(@(x) (~isa(x,'LiteralNode') && numel(x.children) <= 1), self.children);
            if ~keepClauses
                % if we don't keep and clauses containing the literal
                if any(literalMatches)
                    % and the literal is a direct child of this clause
                    % we empty this node
                    self.children(:) = [];
                    return
                end
            end
            % otherwise, we only remove the child.
            % and check for one element entries
            mergeChildren = arrayfun(@(x) ~isa(x,'LiteralNode') && numel(x.children) == 1, self.children);
            toDelete = literalMatches|emptyChildren;
            if any(mergeChildren)
                childsToMerge = self.children(mergeChildren);
                childrenToAdd = AndNode();
                for i = 1:numel(childsToMerge)
                    cchild = childsToMerge(i);
                    childrenToAdd(i) = cchild.children;
                end
                % the following works only on 2017b or newer, but is more
                % efficient.
                %childrenToAdd = arrayfun(@(x) x.children, self.children(mergeChildren));
            end
            self.children(toDelete) = [];
            if exist('childrenToAdd','var')
                for child = 1:numel(childrenToAdd)
                    newChild = childrenToAdd(child);
                    self.children(end+1) = newChild;
                    newChild.parent = self;
                end
            end
            %fprintf('Removing Literal %s from the following node:\n%s\nLeads to the node:\n%s\n',literalID,originalNodeString,self.toString(1));
        end

        function res = toString(self,PipeAnd)
            if nargin < 2
                PipeAnd = 0;
            end
            res = '';
            for i=1:numel(self.children)
                child = self.children(i);
                if PipeAnd
                    res = [res child.toString(PipeAnd) ' & '];
                else
                    res = [res child.toString(PipeAnd) ' and '];
                end

            end
            if length(res) > 2
                if PipeAnd
                    res = res(1:end-3);
                else
                    res = res(1:end-5);
                end
            end
        end
        function cnfNode = convertToCNF(self)
            cnfNode = AndNode();
            for i = 1:numel(self.children)
                if isa(self.children(i),'LiteralNode')
                    CNFChild = self.children(i).copy();
                else
                    CNFChild = self.children(i).convertToCNF();
                end
                cnfNode.addChild(CNFChild);
            end
        end

        function dnfNode = convertToDNF(self)
            dnfNode = OrNode();
            childNodes = [];
            sizes = [];
            for c=1:numel(self.children)
                child = self.children(c);
                if isempty(childNodes)
                    childNodes = child.convertToDNF();
                else
                    childNodes(end+1) = child.convertToDNF();
                end
                convNode = childNodes(end);
                sizes(end+1) = numel(convNode.children);
            end
            %Now make and combinations of all items in the children
            step = ones(numel(sizes),1);
            while self.isValid(sizes,step)
                nextNode = AndNode();
                for i=1:numel(step)
                    convNode = childNodes(i);
                    if strcmp(class(convNode),'LiteralNode')
                        nextNode.addChild(convNode);
                    else
                        nextNode.addChild(convNode.children(step(i)));
                    end
                end
                dnfNode.addChild(nextNode);
                step = self.nextcombination(sizes,step);
            end
            %finally, remove all duplicate literal nodes from this node.

            for c = 1:numel(dnfNode.children)
                literals = {};
                childrenToRemove = [];
                childNode = dnfNode.children(c);
                for i = 1 : numel(childNode.children)
                    if isa(childNode.children(i),'LiteralNode')
                        if ~any(~cellfun(@isempty, strfind(literals,childNode.children(i).toString())))
                            literals{end+1} = childNode.children(i).toString();
                        else
                            childrenToRemove(end+1) = i;
                        end
                    end
                end
                childNode.children(childrenToRemove) = [];
            end
        end

        function res = isValid(self,sizes,step)
            % Check whether a given step is a valid possibility (no step
            % element larger than sizes
            % USAGE:
            %    res = Node.isValid(sizes,step)
            %
            % INPUTS:
            %    sizes:     An array of sizes
            %    step:      An array of suggested selections
            %
            % OUTPUTS:
            %    res:       ~any(step > sizes')
            %
            res = ~any(step > sizes');
        end

        function combination = nextcombination(self,sizes,step)
            % Get the next combination given the current combination
            % USAGE:
            %    combination = Node.nextcombination(sizes,step)
            %
            % INPUTS:
            %    sizes:     An array of maximal sizes
            %    step:      The current combination
            %
            % OUTPUTS:
            %    combination:   The next allowed element of step
            %                   incremented, and potentially others reset
            %                   to 1.
            %
            combination = step;
            combination(1) = combination(1) + 1;
            for i=1:numel(sizes)
                if combination(i) > sizes(i)
                    if i < numel(sizes)
                        combination(i) = 1;
                        combination(i+1) = combination(i+1)  + 1;
                    end
                else
                    break;
                end
            end
        end

        function reduce(self)
            %we can merge any children of and nodes directly.
            mergeNode.children = [];
            childrenChanged = false;
            for i = 1:numel(self.children)
                cchild = self.children(i);
                cchild.reduce()
                %Check if the child has exactly one child. I
                if numel(cchild.children) == 1
                    %If there is only one child, we can directly add the
                    %child to this node.
                    mergeNode.children = [mergeNode.children,cchild.children];
                    childrenChanged = true;
                elseif isa(cchild,'AndNode')
                    mergeNode.children = [mergeNode.children,cchild.children];
                    childrenChanged = true;
                else
                    mergeNode.children = [mergeNode.children,cchild];
                end
            end
            if childrenChanged
                self.children = mergeNode.children;
                for i = 1:numel(self.children)
                    self.children(i).parent = self;
                end
            end
        end
    end

end

author: opencobra

repo: cobratoolbox


text:
function [ obj ] = vb_init( obj )
% Initialize posterior and auxiliary variabels for VB inversion of HUGE.
% Requires obj.dcm and obj.prior to be intialized.
%
% This is a protected method of the tapas_Huge class. It cannot be called
% from outside the class.
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch)
% Copyright (C) 2019 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This software is provided "as is", without warranty of any kind, express
% or implied, including, but not limited to the warranties of
% merchantability, fitness for a particular purpose and non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

% method for gradient calculation
obj.options.fncBold = @bold_grad_cd;

%% prior
prior = struct();
prior.alpha_0 = obj.prior.alpha_0;
prior.m_0 = obj.prior.m_0; % prior mean
prior.tau_0 = obj.prior.tau_0;
prior.nu_0 = max(obj.prior.nu_0, 1.5^obj.idx.P_c);
prior.nu_0 = min(prior.nu_0, double(realmax('single')));
% scaling matrix of inverse-Wishart
prior.S_0 = obj.prior.S_0*(prior.nu_0 - obj.idx.P_c - 1);
% prior mean of hemodynamic parameters
prior.mu_h = obj.prior.mu_h;
% prior Covariance of hemodynamic parameters
prior.Sigma_h = obj.prior.Sigma_h;
% prior inverse scale of observation noise (b_0 in Figure 1 of REF [1])
prior.b_0 = obj.prior.mu_lambda./obj.prior.s2_lambda;
% prior shape parameter of observation noise (a_0 in Figure 1 of REF [1])
prior.a_0 = obj.prior.mu_lambda.^2./obj.prior.s2_lambda;
% prior mean and covariance over confound coefficients
prior.m_beta_0 = obj.prior.m_beta_0;
prior.S_beta_0 = obj.prior.S_beta_0;

obj.prior = prior;

%% create posterior struct
obj.posterior = struct();
% cluster weights
obj.posterior.alpha = obj.prior.alpha_0;
obj.posterior.alpha(1) = obj.posterior.alpha(1) + obj.N;
% cluster parameters
obj.posterior.m = obj.options.start.clusters;
obj.posterior.tau = repmat(obj.prior.tau_0, obj.K, 1);
obj.posterior.tau(1) = obj.prior.tau_0 + obj.N;
obj.posterior.nu = repmat(obj.prior.nu_0, obj.K, 1);
obj.posterior.nu(1) = obj.prior.nu_0 + obj.N;
obj.posterior.S = repmat(obj.prior.S_0, 1, 1, obj.K);

% assigments
obj.posterior.q_nk = [ones(obj.N, 1),zeros(obj.N, obj.K - 1)];
% DCM parameters
obj.posterior.mu_n = obj.options.start.subjects;
Sigma = zeros(obj.idx.P_c + obj.idx.P_h);
Sigma(1:obj.idx.P_c, 1:obj.idx.P_c) = obj.prior.S_0/(obj.prior.nu_0 - obj.idx.P_c - 1);
Sigma(1:obj.idx.P_c, 1:obj.idx.P_c) = obj.prior.S_0;
Sigma(obj.idx.P_c+1:end, obj.idx.P_c+1:end) = obj.prior.Sigma_h;
obj.posterior.Sigma_n = repmat(Sigma, 1, 1, obj.N);
% noise
obj.posterior.b = repmat(obj.prior.b_0, obj.N, obj.R);
obj.posterior.a = repmat(obj.prior.a_0, obj.N, obj.R);
% confounds
if obj.options.confVar
    if obj.options.confVar > 1
        obj.options.confVar = obj.K;
    end
    obj.posterior.m_beta = repmat(obj.prior.m_beta_0, 1, obj.idx.P_c, ...
        obj.options.confVar);
    obj.posterior.S_beta = repmat(obj.prior.S_beta_0, 1, 1, obj.idx.P_c, ...
        obj.options.confVar);
else
    obj.posterior.m_beta = [];
    obj.posterior.S_beta = [];
end

%% auxiliary variables
obj.aux.epsilon = cell(obj.N, 1);   % residual
obj.aux.G       = cell(obj.N, 1);   % jacobian
if obj.options.confVar    % confounds
    % center and normalize confound
    obj.aux.x_n = [obj.data(:).confounds]';
    obj.aux.x_n = bsxfun(@minus, obj.aux.x_n, mean(obj.aux.x_n));
    tmp = std(obj.aux.x_n);
    tmp(tmp < eps) = 1;
    obj.aux.x_n = bsxfun(@rdivide, obj.aux.x_n, tmp);
    % confound-related auxiliary variabels
    obj.aux.x_n_2 = zeros(obj.M, obj.M, obj.N);
    tmp = zeros(obj.idx.P_c, obj.idx.P_c, obj.N);
    for n = 1:obj.N
        obj.aux.x_n_2(:,:,n) = obj.aux.x_n(n,:)'*obj.aux.x_n(n,:);
        tmp(:,:,n) = obj.aux.x_n(n,:)*obj.prior.S_beta_0*obj.aux.x_n(n,:)'...
            *eye(obj.idx.P_c);
    end
    obj.aux.Sigma_beta = repmat(tmp, 1, 1, 1, obj.options.confVar);
    obj.aux.Pi_beta_0 = inv(obj.prior.S_beta_0);
    obj.aux.Pi_m_beta_0 = obj.aux.Pi_beta_0*obj.prior.m_beta_0;
    obj.aux.ldS_beta = repmat(tapas_huge_logdet(obj.prior.S_beta_0), ...
        obj.idx.P_c, obj.options.confVar);
    obj.aux.mu_beta = zeros(obj.N, obj.idx.P_c, obj.options.confVar);
    for k = 1:obj.options.confVar
        obj.aux.mu_beta(:,:,k) = obj.aux.x_n*obj.posterior.m_beta(:,:,k);
    end
end

obj.aux.q_k = sum(obj.posterior.q_nk, 1)' + realmin;
obj.aux.ldS = repmat(tapas_huge_logdet(obj.prior.S_0), obj.K, 1);
obj.aux.ldSigma = repmat(tapas_huge_logdet(Sigma), obj.N ,1);
obj.aux.mu_prime_h = obj.prior.mu_h/obj.prior.Sigma_h;
obj.aux.Pi_h = inv(obj.prior.Sigma_h);
% nu_k*inv(S_k)
obj.aux.nu_inv_S = zeros(size(obj.posterior.S));
for k = 1:obj.K
    obj.aux.nu_inv_S(:,:,k) = obj.posterior.nu(k).*...
        inv(obj.posterior.S(:,:,k));
end
obj.aux.b_prime = zeros(obj.N, obj.R);

%% check stability of starting values
obj.aux.q_r = zeros(obj.N, 1);
for n = 1:obj.N
    obj.aux.q_r(n) = size(obj.data(n).bold, 1);
    obj = obj.options.fncBold(obj, n);
    assert(~(any(isnan(obj.aux.G{n}(:))) || any(isinf(obj.aux.G{n}(:)))), ...
        'TAPAS:HUGE:Init:Stability',...
        'Starting values for subject %u lead to instable DCM.', n);

    tmp = sum((obj.aux.G{n}*obj.posterior.Sigma_n(:,:,n)).*obj.aux.G{n}, 2);
    obj.aux.b_prime(n,:) = sum(obj.aux.epsilon{n}.^2 + ...
        reshape(tmp, obj.aux.q_r(n), obj.R), 1);
    obj.posterior.b(n,:) = obj.prior.b_0 + obj.aux.b_prime(n,:)/2;
    obj.posterior.a(n,:) = obj.prior.a_0 + obj.aux.q_r(n)/2;
end

%% initialize auxiliary variables
obj.aux.lambda_bar = obj.posterior.a./obj.posterior.b;

obj.aux.mu_k_c(1,:) = mean(obj.posterior.mu_n(:, 1:obj.idx.P_c), 1);
obj.aux.Sigma_k_c(:,:,1) = ...
    sum(obj.posterior.Sigma_n(1:obj.idx.P_c, 1:obj.idx.P_c, :), 3);

% offline part of negative free energy
F_off = -obj.K*(obj.idx.P_c*(obj.idx.P_c - 1)/4*log(pi) + ...
    sum(gammaln((obj.prior.nu_0-obj.idx.P_c+1:obj.prior.nu_0)/2)));
F_off = F_off +.5*obj.K*obj.prior.nu_0*tapas_huge_logdet(obj.prior.S_0);
F_off = F_off +.5*obj.K*obj.idx.P_c*log(obj.prior.tau_0);
F_off = F_off + gammaln(sum(obj.prior.alpha_0));
F_off = F_off - sum(gammaln(obj.prior.alpha_0));
F_off = F_off - gammaln(obj.N + sum(obj.prior.alpha_0));
F_off = F_off - obj.N*obj.R*gammaln(obj.prior.a_0);
F_off = F_off + sum(gammaln(obj.posterior.a(:)));
F_off = F_off + obj.N*obj.R*obj.prior.a_0*log(obj.prior.b_0);
F_off = F_off -.5*obj.N*tapas_huge_logdet(obj.prior.Sigma_h);
F_off = F_off -.5*obj.idx.P_c*tapas_huge_logdet(obj.prior.S_beta_0);
F_off = F_off +.5*obj.N*obj.idx.P_c*log(2);
F_off = F_off -.5*sum(obj.aux.q_r)*obj.R*log(2*pi);
F_off = F_off +.5*obj.N*(obj.idx.P_c + obj.idx.P_h);
F_off = F_off +.5*obj.idx.P_c*(obj.K*(1 + obj.prior.nu_0) + obj.N);
if obj.options.confVar
    F_off = F_off +.5*obj.M*obj.idx.P_c*obj.options.confVar;
end
obj.aux.F_off = F_off;

% evaluate negative free energy (VB)
obj.posterior.nfe = obj.vb_nfe( );

%% initialize trace
obj.trace = struct();
obj.trace.nfe = NaN;
obj.trace.nDcmUpdate = zeros(obj.N, 1);
obj.trace.nRetract   = zeros(obj.N,1);

end

author: translationalneuromodeling

repo: tapas


text:
function [x, E, s, lambdamax] = tvdip(y, lambda, display, stoptol, maxiter)
% Performs discrete total variation denoising (TVD) using a primal-dual
% interior-point solver. It minimizes the following discrete functional:
%
%  E=(1/2)||y-x||_2^2+lambda*||Dx||_1,
%
% over the variable x, given the input signal y, according to each
% value of the regularization parameter lambda > 0. D is the first
% difference matrix. Uses hot-restarts from each value of lambda to speed
% up convergence for subsequent values: best use of this feature is made by
% ensuring that the chosen lambda values are close to each other.
%
% Usage:
% [x, E, s] = tvdip(y, lambda, display, stoptol, maxiter)
%
% Input arguments:
% - y          Original signal to denoise, size N x 1.
% - lambda     A vector of positive regularization parameters, size L x 1.
%              TVD will be applied to each value in the vector.
% - display    (Optional) Set to 0 to turn off progress display, 1 to turn
%              on. If not specifed, defaults to progress display on.
% - stoptol    (Optional) Precision as determined by duality gap tolerance,
%              if not specified, defaults to 1e-3.
% - maxiter    (Optional) Maximum interior-point iterations, if not
%              specified defaults to 60.
%
% Output arguments:
% - x          Denoised output signal for each value of lambda, size N x L.
% - E          Objective functional at minimum for each lambda, size L x 1.
% - s          Optimization result, 1 = solved, 0 = maximum iterations
%              exceeded before reaching duality gap tolerance, size L x 1.
% - lambdamax  Maximum value of lambda for the given y. If
%              lambda >= lambdamax, the output is the trivial constant
%              solution x = mean(y).
%
% (c) Max Little, 2010. Based around code originally written by
% S.J. Kim, K. Koh, S. Boyd and D. Gorinevsky. If you use this code for
% your research, please cite:
% M.A. Little, Nick S. Jones (2010)
% "Sparse Bayesian Step-Filtering for High-Throughput Analysis of Molecular
% Machine Dynamics", in 2010 IEEE International Conference on Acoustics,
% Speech and Signal Processing, 2010, ICASSP 2010 Proceedings.
%
% This code is released under the terms of GNU General Public License as
% published by the Free Software Foundation; version 2 or later.

error(nargchk(2,5,nargin));
if (nargin < 3)
    display = 1;
end
if (nargin < 4)
    stoptol = 1e-3;
end
if (nargin < 5)
    maxiter = 60;
end

y = y(:);

% Search tuning parameters
ALPHA     = 0.01;   % Backtracking linesearch parameter (0,0.5]
BETA      = 0.5;    % Backtracking linesearch parameter (0,1)
MAXLSITER = 20;     % Max iterations of backtracking linesearch
MU        = 2;      % t update

N = length(y);    % Length of input signal y
M = N-1;          % Size of Dx

% Construct sparse operator matrices
I1 = speye(M,M);
O1 = spalloc(M,1,M);
D = [I1 O1]-[O1 I1];

DDT = D*D';
Dy  = D*y;

% Find max value of lambda
lambdamax = max(abs(DDT\Dy));

if (display)
    fprintf('lambda_max=%5.2e\n', lambdamax);
end

L = length(lambda);
x = zeros(N, L);
s = zeros(L, 1);
E = zeros(L, 1);

% Optimization variables set up once at the start
z    = zeros(M,1);   % Dual variable
mu1  = ones(M,1);    % Dual of dual variable
mu2  = ones(M,1);    % Dual of dual variable

% Work through each value of lambda, with hot-restart on optimization
% variables
for l = 1:L

    t    =  1e-10;
    step =  Inf;
    f1   =  z-lambda(l);
    f2   = -z-lambda(l);

    % Main optimization loop
    s(l) = 1;

    if (display)
        fprintf('Solving for lambda=%5.2e, lambda/lambda_max=%5.2e\nIter# Primal    Dual      Gap\n', ...
            lambda(l), lambda(l)/lambdamax);
    end
    for iters = 0:maxiter

        DTz  = (z'*D)';
        DDTz = D*DTz;
        w    = Dy-(mu1-mu2);

        % Calculate objectives and primal-dual gap
        pobj1 = 0.5*w'*(DDT\w)+lambda(l)*sum(mu1+mu2);
        pobj2 = 0.5*DTz'*DTz+lambda(l)*sum(abs(Dy-DDTz));
        pobj = min(pobj1,pobj2);
        dobj = -0.5*DTz'*DTz+Dy'*z;
        gap  = pobj - dobj;

        if (display)
            fprintf('%5d %7.2e %7.2e %7.2e\n', iters, pobj, dobj, gap);
        end

        % Test duality gap stopping criterion
        if (gap <= stoptol)
            s(l) = 1;
            break;
        end

        if (step >= 0.2)
            t = max(2*M*MU/gap, 1.2*t);
        end

        % Do Newton step
        rz      =  DDTz - w;
        S       =  DDT-sparse(1:M,1:M,mu1./f1+mu2./f2);
        r       = -DDTz + Dy + (1/t)./f1 - (1/t)./f2;
        dz      =  S\r;
        dmu1    = -(mu1+((1/t)+dz.*mu1)./f1);
        dmu2    = -(mu2+((1/t)-dz.*mu2)./f2);

        resDual = rz;
        resCent = [-mu1.*f1-1/t; -mu2.*f2-1/t];
        residual= [resDual; resCent];

        % Perform backtracking linesearch
        negIdx1 = (dmu1 < 0);
        negIdx2 = (dmu2 < 0);
        step = 1;
        if (any(negIdx1))
            step = min( step, 0.99*min(-mu1(negIdx1)./dmu1(negIdx1)) );
        end
        if (any(negIdx2))
            step = min( step, 0.99*min(-mu2(negIdx2)./dmu2(negIdx2)) );
        end

        for liter = 1:MAXLSITER
            newz    =  z  + step*dz;
            newmu1  =  mu1 + step*dmu1;
            newmu2  =  mu2 + step*dmu2;
            newf1   =  newz - lambda(l);
            newf2   = -newz - lambda(l);

            % Update residuals
            newResDual  = DDT*newz - Dy + newmu1 - newmu2;
            newResCent  = [-newmu1.*newf1-1/t; -newmu2.*newf2-1/t];
            newResidual = [newResDual; newResCent];

            if ((max(max(newf1),max(newf2)) < 0) && ...
                    (norm(newResidual) <= (1-ALPHA*step)*norm(residual)))
                break;
            end
            step = BETA*step;
        end

        % Update primal and dual optimization parameters
        z = newz;
        mu1 = newmu1;
        mu2 = newmu2;
        f1 = newf1;
        f2 = newf2;
    end

    x(:,l) = y-D'*z;
    E(l) = 0.5*sum((y-x(:,l)).^2)+lambda(l)*sum(abs(D*x(:,l)));

    % We may have a close solution that does not satisfy the duality gap
    if (iters >= maxiter)
        s(l) = 0;
    end

    if (display)
        if (s(l))
            fprintf('Solved to precision of duality gap %5.2e\n', gap);
        else
            fprintf('Max iterations exceeded - solution may be inaccurate\n');
        end
    end

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
clc;
clear;
close all;

T0 = 500 ; % initial temperature
r = 0.997 ; % temperature damping rate
Ts = 1 ; % stop temperature
iter = 300;

model = initModel();
set(gcf,'unit','normalized','position',[0,0.35,1,0.7]);

flag = 0;

% initialization
while(1)
route = randomSol(model);
if(isFeasible(route,model))
    break;
end
end

cost = calculateCost(route,model);
T = T0;

cnt = 1;
minCost = cost;
minRoute = route;

maxIterate = 2100;
costArray = zeros(maxIterate,1);

% SA
while(T > Ts)
    for k = 1:iter
    mode = randi([1 3]);
    newRoute = createNeibor(route,model,mode);
    newCost = calculateCost(newRoute,model);
    delta = newCost - cost;

    if(delta < 0)
        cost = newCost;
        route = newRoute;
    else
        p=exp(-delta/T);
        if rand() <= p
             cost = newCost;
             route = newRoute;

        end
    end
    end

    costArray(cnt) = cost;
    if cost<minCost
        minCost = cost;
        minRoute = route;
        flag = 1;
    end

    T = T*r; %  annealing
    disp(['Iteration ' num2str(cnt) ': BestCost = ' num2str(minCost) ': CurrentCost = ' num2str(cost) ' T = ' num2str(T)]);
    cnt = cnt+1;

   figure(1);
   if(flag == 1)
   plotSolution(minRoute,model);
    flag = 0;
   end
%    figure(2);
   subplot(1,2,2)
   plot(costArray);

   pause(0.0001);
end

author: lzane

repo: VRP-using-SA-with-Matlab


text: find(strcmp(labels,'occluder'))

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function depthRefined = SiftFuv2warp(sequenceName,frameIDtarget, interval, renderK, renderW, renderH)

if ~exist('sequenceName','var')
    sequenceName = 'Shuran/bedroom_funiturestore/2014-05-26_14-41-22_260595134347';
end

if ~exist('frameIDtarget','var')
    frameIDtarget = 34; %200;
end

if ~exist('interval','var')
    interval = 150;
end

data = loadSUN3Dv2(sequenceName);

frameIDs = max(frameIDtarget-interval,1):min(data.imageTotalFrames,frameIDtarget+interval);

frames = getRGBDframe(data,frameIDs);

Rts = repmat([eye(3) zeros(3,1)],[1,1,length(frameIDs)]);

itarget = find(frameIDs==frameIDtarget);

%{
depthOriginal = get_depth(imread(data.depth{frameIDs(itarget)}));
depth = denoise(depthOriginal,data.camera.D);
PCtarget(1,:)=reshape(data.camera.D.X .* depth,1,[]);
PCtarget(2,:)=reshape(data.camera.D.Y .* depth,1,[]);
PCtarget(3,:)=reshape(depth .* (~isnan(data.camera.D.X) & ~isnan(data.camera.D.Y)),1,[]);
PCtarget = PCtarget(:, PCtarget(3,:)>0);
%}

matched = false(1,length(frameIDs));
for i = 1:length(frameIDs)
    if i~=itarget
        fprintf('matching frame %d and frame %d\n',frameIDs(itarget),frameIDs(i));
        Rts(:,:,i) = align2RGBD(frames(:,:,:,itarget), frames(:,:,:,i));
        matched(i) = true;
        %{
        depth = get_depth(imread(data.depth{frameIDs(i)}));
        depth = denoise(depth,data.camera.D);
        clear PCi
        PCi(1,:)=reshape(data.camera.D.X .* depth,1,[]);
        PCi(2,:)=reshape(data.camera.D.Y .* depth,1,[]);
        PCi(3,:)=reshape(depth .* (~isnan(data.camera.D.X) & ~isnan(data.camera.D.Y)),1,[]);
        PCi = PCi(:, PCi(3,:)>0);

        try
            Rts(:,:,i) = align2RGBD(frames(:,:,:,itarget), frames(:,:,:,i), PCtarget, PCi);
            matched(i) = true;
        catch
        end
        %}
    end
end

% old Kinect
%renderK = [519.1638 0 320; 0 519.1638 240; 0 0 1];
if ~exist('renderK','var')
    renderK = data.camera.RGB.K2render; %[519 0 320; 0 519 240; 0 0 1];
end
if ~exist('renderW','var')
    renderW = data.camera.RGB.width; %640;
end
if ~exist('renderH','var')
    renderH = data.camera.RGB.height; %480;
end

% new Kinect
%{
renderK = data.camera.D.KK;
renderW = data.camera.D.width;
renderH = data.camera.D.height;
%}

for i=1:length(frameIDs)
    if matched(i) || i==itarget
        frameID=frameIDs(i);

        XYZcamera = frames(:,:,[5 6 4],i);
        XYZcamera(:,:,4) = XYZcamera(:,:,3)>0;
        %{
        depth = get_depth(imread(data.depth{frameID}));
        XYZcamera(:,:,1)=data.camera.D.X .* depth;
        XYZcamera(:,:,2)=data.camera.D.Y .* depth;
        XYZcamera(:,:,3)=depth .* (~isnan(data.camera.D.X) & ~isnan(data.camera.D.Y));
        XYZcamera(:,:,4)=depth>0 & ~isnan(data.camera.D.X) & ~isnan(data.camera.D.Y);
        %}

        [~,undistortDepth(:,:,i)] = WarpDepthMatlab(XYZcamera, renderK, Rts(:,:,i), renderW, renderH);
    end
end

undistortDepth(undistortDepth(:)==0) = NaN;

%depthRefined = nanmean(undistortDepth,3);
depthMedian = nanmedian(undistortDepth,3);

% to preserve details: check the value and the range of all values. if it is 25%-75%, keep the original values
minV = prctile(undistortDepth,75,3);
maxV = prctile(undistortDepth,25,3);
depthRaw = undistortDepth(:,:,itarget);

selV = (depthRaw < minV) | (depthRaw > maxV) | isnan(depthRaw);
ind = find(selV);
depthRefined = depthRaw;
depthRefined(ind) = depthMedian(ind);

depthRefined(isnan(depthRefined(:)))=0;

%{
[x,y] = meshgrid(1:renderW, 1:renderH);
clear X
X(1,:) = reshape((x-renderK(1,3)).*depthRefined/renderK(1,1),1,[]);
X(2,:) = reshape((y-renderK(2,3)).*depthRefined/renderK(2,2),1,[]);
X(3,:) = reshape(depthRefined,1,[]);
X = X(:,~isnan(depthRefined(:)));

% turn

figure;
subplot(2,2,1);
imagesc(depthOriginal,[0 20]); axis equal; axis tight; axis off; title('raw TOF depth');

subplot(2,2,2);
imagesc(depthRefined,[0 20]); axis equal; axis tight; axis off; title('undistorted denoised depth');

points2ply('pt_denoise101.ply', X);
%}

author: jianxiongxiao

repo: ProfXkit


text:
% a script to reproduce table 1 of paper :
%
%   ``Rotation, Scaling and Deformation Invariant Scattering
%   for Texture Discrimination"
%   Laurent Sifre, Stephane Mallat
%   Proc. IEEE CVPR 2013 Portland, Oregon
%
% Scattering classification rates for UIUC databases
%
% NOTE : Computing the scattering for the whole database takes time.
% We provide PRECOMPUTED scattering in the files
%   precomputed/uiuc/trans_scatt.mat
%   precomputed/uiuc/roto_trans_scatt.mat
%   precomputed/uiuc/roto_trans_scatt_log.mat
%   precomputed/uiuc/roto_trans_scatt_log_scale_avg.mat
%   precomputed/uiuc/roto_trans_scatt_log_scale_avg_multiscal_train.mat
% You can COMPUTE THE SCATTERING YOURSELF by changing the following line :
%   use_precomputed_scattering = 0;
%
% DOWNLOAD : The UIUC database can be downloaded at
%   http://www.nada.kth.se/cvap/databases/uiuc/kth_tips_grey_200x200.tar

%% load the database
clear; close all;
% NOTE : the following line must be modified with the path to the
% uiuc database in YOUR system.
path_to_db = '/Users/laurentsifre/TooBigForDropbox/Databases/cvr_uiuc';
src = uiuc_src(path_to_db);
db_name = 'uiuc';

use_precomputed_scattering = 0; % change to 0 to skip computation of scattering

grid_train = [5, 10, 20]; % number of training for classification
nb_split = 10; % number of split for classification

%% ---------------------------------------------------
%% ----------------- trans_scatt ---------------------
%% ---------------------------------------------------

%% compute scattering of all images in the database
feature_name = 'trans_scatt';
precomputed_path = sprintf('./precomputed/%s/%s.mat', db_name, feature_name);
if (use_precomputed_scattering)
    load(precomputed_path);
else
    %configure scattering
    options.J = 5; % number of octaves
    options.Q = 1; % number of scales per octave
    options.M = 2; % scattering orders

    % build the wavelet transform operators for scattering
    Wop = wavelet_factory_2d_spatial(options, options);

    % a function handle that
    %   - read the image
    %   - resize it to 200x200
    %   - compute its scattering
    fun = @(filename)(scat(imresize_notoolbox(imreadBW(filename)...
        ,[200 200]), Wop));

    % compute all scattering
    % (800 seconds on a 2.4 Ghz Intel Core i7)
    trans_scatt_all = srcfun(fun, src);

    % a function handle that
    %   - format the scattering in a 3d matrix
    %   - remove margins
    %   - average accross position
    % (10 seconds on a 2.4 Ghz Intel Core i7)
    fun = @(Sx)(mean(mean(remove_margin(format_scat(Sx),0),2),3));
    trans_scatt = cellfun_monitor(fun ,trans_scatt_all);

    % save scattering
    save(precomputed_path, 'trans_scatt');
end
% format the database of feature
db = cellsrc2db(trans_scatt, src);

%% classification
rsds_classif(db, db_name, feature_name, grid_train, nb_split);

%% ---------------------------------------------------
%% --------------- roto_trans_scatt ------------------
%% ---------------------------------------------------

%% compute scattering of all images in the database
feature_name = 'roto_trans_scatt';
precomputed_path = sprintf('./precomputed/%s/%s.mat', db_name, feature_name);
if (use_precomputed_scattering)
    load(precomputed_path);
else
    % configure scattering
    options.J = 5; % number of octaves
    options.Q = 1; % number of scales per octave
    options.M = 2; % scattering orders

    % build the wavelet transform operators for scattering
    Wop = wavelet_factory_3d_spatial(options, options, options);

    % a function handle that
    %   - read the image
    %   - resize it to 200x200
    %   - compute its scattering
    fun = @(filename)(scat(imresize_notoolbox(imreadBW(filename),[200 200]), Wop));
    % (1000 seconds on a 2.4 Ghz Intel Core i7)
    roto_trans_scatt_all = srcfun(fun, src);

    % a function handle that
    %   - format the scattering in a 3d matrix
    %   - remove margins
    %   - average accross position
    fun = @(Sx)(mean(mean(remove_margin(format_scat(Sx),1,[2,3]),2),3));
    % (10 seconds on a 2.4 Ghz Intel Core i7)
    roto_trans_scatt = cellfun_monitor(fun ,roto_trans_scatt_all);

    % save scattering
    save(precomputed_path, 'roto_trans_scatt');
end
% format the database of feature
db = cellsrc2db(roto_trans_scatt, src);

%% classification
rsds_classif(db, db_name, feature_name, grid_train, nb_split);

%% ---------------------------------------------------
%% ------------- roto_trans_scatt_log ----------------
%% ---------------------------------------------------

feature_name = 'roto_trans_scatt_log';
precomputed_path = sprintf('./precomputed/%s/%s.mat', db_name, feature_name);
if (use_precomputed_scattering)
    load(precomputed_path);
else
    % a function handle that
    %   - format the scattering in a 3d matrix
    %   - take the logarithm
    %   - remove margins
    %   - average accross position
    fun = @(Sx)(mean(mean(log(remove_margin(format_scat(Sx),1,[2,3])),2),3));
    roto_trans_scatt_log = cellfun_monitor(fun ,roto_trans_scatt_all);

    %save scattering
    save(precomputed_path);
end
% format the database of feature
db = cellsrc2db(roto_trans_scatt_log, src);

%% classification
rsds_classif(db, db_name, feature_name, grid_train, nb_split);

%% ---------------------------------------------------
%% ---------- roto_trans_scatt_log_scale_avg ---------
%% ---------------------------------------------------

%% compute scattering of all images in the database

feature_name = 'roto_trans_scatt_log_scale_avg';
precomputed_path = sprintf('./precomputed/%s/%s.mat', db_name, feature_name);

if (use_precomputed_scattering)
    load(precomputed_path);
else
    % configure scattering
    options.J = 5; % number of octaves
    options.Q = 1; % number of scales per octave
    options.M = 2; % scattering orders

    % build the wavelet transform operators for scattering
    Wop = wavelet_factory_3d_spatial(options, options, options);

    % a function handle that compute scattering given an image
    fun = @(x)(scat(x, Wop));

    % another function handle that
    %   - read the image
    %   - resize it to 200x200
    %   - apply fun to all scaled version of the image
    multi_fun = @(filename)(fun_multiscale(fun, ...
        imresize_notoolbox(imreadBW(filename),[200 200]), sqrt(2), 4));

    % (2748 seconds on a 2.4 Ghz Intel Core i7)
    roto_trans_scatt_multiscale = srcfun(multi_fun, src);

    %% log + spatial average
    fun = @(Sx)(mean(mean(log(remove_margin(format_scat(Sx),1,[2,3])),2),3));
    multi_fun = @(x)(cellfun_monitor(fun, x));
    roto_trans_scatt_multiscale_log_sp_avg = cellfun_monitor(multi_fun, roto_trans_scatt_multiscale);

    %% scale average
    fun = @(x)(mean(cell2mat(x),2));
    roto_trans_scatt_log_scale_avg = cellfun_monitor(fun, roto_trans_scatt_multiscale_log_sp_avg);

    % save scattering
    save(precomputed_path, 'roto_trans_scatt_log_scale_avg');

    % format the database of feature
    db = cellsrc2db(roto_trans_scatt_log_scale_avg, src);
end

%% classification
rsds_classif(db, db_name, feature_name, grid_train, nb_split);

%% ---------------------------------------------------
%% - roto_trans_scatt_log_scale_avg_multiscal_train --
%% ---------------------------------------------------

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
% Algorithm for finding connected components in a graph
% Valid for undirected graphs only
% INPUTS: adj - adjacency matrix
% OUTPUTS: a list of the components comp{i}=[j1,j2,...jk}
% Other routines used: find_conn_compI.m (embedded), degrees.m, kneighbors.m
% GB, Last updated: October 2, 2009

function comp_mat = find_conn_comp(adj)

[deg,~,~]=degrees(adj);   % degrees
comp_mat={};                       % initialize components matrix

for i=1:length(deg)
    if deg(i)>0
        done=0;
        for x=1:length(comp_mat) % for every component check whether i is included
            if length(find(comp_mat{x}==i))>0   % i in comp_mat(x).mat
                done=1;
                break
            end
        end
        if not(done)
            comp=find_conn_compI(adj,i);
            comp_mat{length(comp_mat)+1}=comp;
        end

    elseif deg(i)==0
        comp_mat{length(comp_mat)+1}=[i];
    end
end

function comp=find_conn_compI(adj,i)
% heuristic for finding the conn component to which "i" belongs to
% works well in practice for large datasets
% INPUTS: adjacency matrix and index of the key node
% OUTPUTS: all node indices of the nodes to which "i" belongs to

neigh1=kneighbors(adj,i,1);
neigh1=unique([neigh1 i]); % add i to its own component

while 1
  len0=length(neigh1);
  for j=1:len0
    neigh2=kneighbors(adj,neigh1(j),1);
    neigh1=unique([neigh1, neigh2]);
  end
  if len0==length(neigh1)
    comp=neigh1;
    return
  end
end
% Compute the total degree, in-degree and out-degree of a graph based on
% the adjacency matrix; should produce weighted degrees, if the input matrix is weighted
% INPUTS: adjacency matrix
% OUTPUTS: degree, indegree and outdegree sequences
% GB, Last Updated: October 2, 2009

function [deg,indeg,outdeg]=degrees(adj)

indeg = sum(adj);
outdeg = sum(adj');

% if isdirected(adj)
  deg = indeg + outdeg; % total degree
%
% else   % undirected graph: indeg=outdeg
%   deg = indeg + diag(adj)';  % add self-loops twice, if any

% end

% Finds the number of k-neighbors (k links away) for every node
% INPUTS: adjacency matrix, node index, k - number of links
% OUTPUTS: vector of k-neighbors indices
% GB, May 3, 2006

function kneigh = kneighbors(adj,ind,k)

adjk = adj;
for i=1:k-1;
    adjk = adjk*adj;
end;

kneigh = find(adjk(ind,:)>0);

f_neigh = find(adjk(:,ind)>0);
kneigh = unique([kneigh, f_neigh']);

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut