text:
% Struct Utility functions
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 7/14/2010
classdef StructUtil
    methods (Static=true)
        function A = catstruct(varargin)
            % CATSTRUCT - concatenate structures
            %
            %   X = CATSTRUCT(S1,S2,S3,...) concates the structures S1, S2, ... into one
            %   structure X.
            %
            %   Example:
            %     A.name = 'Me' ;
            %     B.income = 99999 ;
            %     X = catstruct(A,B)
            %     % -> X.name = 'Me' ;
            %     %    X.income = 99999 ;
            %
            %   CATSTRUCT(S1,S2,'sorted') will sort the fieldnames alphabetically.
            %
            %   If a fieldname occurs more than once in the argument list, only the last
            %   occurence is used, and the fields are alphabetically sorted.
            %
            %   To sort the fieldnames of a structure A use:
            %     A = CATSTRUCT(A,'sorted') ;
            %
            %   To concatenate two similar array of structs use simple concatenation:
            %     A = dir('*.mat') ; B = dir('*.m') ; C = [A ; B] ;
            %
            %   When there is nothing to concatenate, the result will be an empty
            %   struct (0x0 struct array with no fields).
            %
            %   See also CAT, STRUCT, FIELDNAMES, STRUCT2CELL

            % for Matlab R13 and up
            % version 2.2 (oct 2008)
            % (c) Jos van der Geest
            % email: jos@jasen.nl

            % History
            % Created:  2005
            % Revisions
            %   2.0 (sep 2007) removed bug when dealing with fields containing cell
            %                  arrays (Thanks to Rene Willemink)
            %   2.1 (sep 2008) added warning and error identifiers
            %   2.2 (oct 2008) fixed error when dealing with empty structs (Thanks to
            %                  Lars Barring)

            N = nargin ;

            error(nargchk(1,Inf,N)) ;

            if ~isstruct(varargin{end}),
                if isequal(varargin{end},'sorted'),
                    sorted = 1 ;
                    N = N-1 ;
                    if N < 1,
                        A = struct([]) ;
                        return
                    end
                else
                    error('catstruct:InvalidArgument','Last argument should be a structure, or the string "sorted".') ;
                end
            else
                sorted = 0 ;
            end

            FN = cell(N,1) ;
            VAL = cell(N,1) ;

            for ii=1:N,
                X = varargin{ii} ;
                if ~isstruct(X),
                    error('catstruct:InvalidArgument',['Argument #' num2str(ii) ' is not a structure.']) ;
                end
                if ~isempty(X),
                    % empty structs are ignored
                    FN{ii} = fieldnames(X) ;
                    VAL{ii} = struct2cell(X) ;
                end
            end

            FN = cat(1,FN{:}) ;
            VAL = cat(1,VAL{:}) ;
            [UFN,ind] = unique(FN) ;

            %BEGIN added by jkarr
            UVAL = {VAL{ind}}';
            for i=1:length(UFN)
                idxs = find(strcmp(FN,UFN{i}));
                if(length(idxs)==1)
                    continue;
                end

                recurse=true;
                for j=1:length(idxs)
                    if(~isstruct(VAL{idxs(j)}))
                        recurse=false;
                        break;
                    end
                end

                if(recurse)
                    UVAL{i} = edu.stanford.covert.util.StructUtil.catstruct(VAL{idxs});
                end
            end
            VAL=UVAL;
            FN=UFN;
            %END added by jkarr

            %BEGIN removed by jkarr
            %if numel(UFN) ~= numel(FN),
            %    warning('catstruct:DuplicatesFound','Duplicate fieldnames found. Last value is used and fields are sorted') ;
            %    sorted = 1 ;
            %end
            %
            %if sorted,
            %    VAL = VAL(ind) ;
            %    FN = FN(ind) ;
            %end
            %END removed by jkarr

            if ~isempty(FN),
                % This deals correctly with cell arrays
                A = cell2struct(VAL, FN);
            else
                A = struct([]) ;
            end
        end

        function diff(struct1, struct2, indent)
            if ~exist('indent', 'var')
                indent = 0;
            end

            fieldNames1 = fieldnames(struct1);
            fieldNames2 = fieldnames(struct2);

            fieldsNotIn2 = setdiff(fieldNames1, fieldNames2);
            fieldsNotIn1 = setdiff(fieldNames2, fieldNames1);

            if ~isempty(fieldsNotIn2)
                fprintf([repmat('\t', 1, indent) 'Fields in struct1 not in struct2:\n' repmat('- %s\n', 1, numel(fieldsNotIn2))], fieldsNotIn2{:});
            end
            if ~isempty(fieldsNotIn1)
                fprintf([repmat('\t', 1, indent) 'Fields in struct2 not in struct1:\n' repmat('- %s\n', 1, numel(fieldsNotIn1))], fieldsNotIn2{:});
            end

            commonFieldNames = intersect(fieldNames1, fieldNames2);
            for i = 1:numel(commonFieldNames)
                val1 = struct1.(commonFieldNames{i});
                val2 = struct2.(commonFieldNames{i});
                class1 = class(val1);
                class2 = class(val2);
                if ~isequal(class1, class2)
                    fprintf([repmat('\t', 1, indent) 'Types of %s''s fields differ (%s v. %s)\n'], commonFieldNames{i}, class1, class2);
                    continue;
                end

                size1 = size(val1);
                size2 = size(val2);
                if ~isequal(size1, size2)
                    fprintf([repmat('\t', 1, indent) 'Sizes of %s''s fields differ (%s v. %s)\n'], commonFieldNames{i}, ...
                        sprintf(['%d' repmat(' x %d', 1, numel(size1)-1)], size1), ...
                        sprintf(['%d' repmat(' x %d', 1, numel(size2)-1)], size2));
                    continue;
                end

                switch class1
                    case 'cell'
                        idxs = find(cellfun(@(x, y) ~isequal(x, y), val1, val2));
                        if isempty(idxs); continue; end;

                        fprintf([repmat('\t', 1, indent) 'Comparing %s fields\n'], commonFieldNames{i});
                        subs = cell(1, ndims(val1));
                        [subs{:}] = ind2sub(size1, idxs);
                        subs = cell2mat(subs);
                        val1 = val1(idxs);
                        for j = 1:numel(idxs)
                            fprintf([repmat('\t', 1, indent+1) '%s\t%s\n'], ...
                                sprintf(['(%d' repmat(', %d', 1, ndims(val1)-1)  ')'], subs(j, :)), ...
                                'differs');
                        end
                    case {'logical', 'int8', 'int16' ,'int32 ','int64', 'uint8', 'uint16', 'uint32', 'uint64', 'double', 'single', 'char'}
                        idxs = find(val1 ~= val2);
                        if isempty(idxs); continue; end;

                        fprintf([repmat('\t', 1, indent) 'Comparing %s fields\n'], commonFieldNames{i});
                        subs = cell(1, ndims(val1));
                        [subs{:}] = ind2sub(size1, idxs);
                        subs = cell2mat(subs);
                        val1 = val1(idxs);
                        val2 = val2(idxs);
                        for j = 1:numel(idxs)
                            fprintf([repmat('\t', 1, indent+1) '%s\t%e\t%e\n'], ...
                                sprintf(['(%d' repmat(', %d', 1, ndims(val1)-1)  ')'], subs(j, :)), ...
                                val1(j), val2(j));
                        end
                    case {'edu.stanford.covert.util.SparseMat', 'edu.stanford.covert.util.CircularSparseMat'}
                        subs = find(val1 ~= val2);
                        if isempty(subs); continue; end;

                        fprintf([repmat('\t', 1, indent) 'Comparing %s fields\n'], commonFieldNames{i});
                        val1 = val1(subs);
                        val2 = val2(subs);
                        for j = 1:numel(idxs)
                            fprintf([repmat('\t', 1, indent+1) '%s\t%e\t%e\n'], ...
                                sprintf(['(%d' repmat(', %d', 1, ndims(val1)-1)  ')'], subs(j, :)), ...
                                val1(j), val2(j));
                        end
                    otherwise
                        fprintf([repmat('\t', 1, indent) 'Comparing %s fields\n'], commonFieldNames{i});
                        edu.stanford.covert.util.StructUtil.diff(val1, val2, indent+1);
                end
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function X_rec = recoverData(Z, U, K)
%RECOVERDATA Recovers an approximation of the original data when using the
%projected data
%   X_rec = RECOVERDATA(Z, U, K) recovers an approximation the
%   original data that has been reduced to K dimensions. It returns the
%   approximate reconstruction in X_rec.
%

% You need to return the following variables correctly.
X_rec = zeros(size(Z, 1), size(U, 1));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the approximation of the data by projecting back
%               onto the original space using the top K eigenvectors in U.
%
%               For the i-th example Z(i,:), the (approximate)
%               recovered data for dimension j is given as follows:
%                    v = Z(i, :)';
%                    recovered_j = v' * U(j, 1:K)';
%
%               Notice that U(j, 1:K) is a row vector.
%
for i=1:size(Z, 1)
   v = Z(i,:)';

   X_rec(i,:) = v'*U(:,1:K)';

end

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

end

author: zzlyw

repo: machine-learning-exercises


text:
function output = calllinprog(interfacedata)

% Standard input interface
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
Q       = interfacedata.Q;
lb      = interfacedata.lb;
ub      = interfacedata.ub;
x0      = interfacedata.x0;

showprogress('Calling LINPROG',options.showprogress);

if isempty(F_struc)
    Aeq = [];
    beq = [];
    A = [];
    b = [];
else
    Aeq = -F_struc(1:1:K.f,2:end);
    beq = F_struc(1:1:K.f,1);
    A =-F_struc(K.f+1:end,2:end);
    b = F_struc(K.f+1:end,1);
end

switch options.verbose
case 0
    options.linprog.Display = 'off';
case 1
    options.linprog.Display = 'final';
otherwise
    options.linprog.Display = 'iter';
end

if isfield(options.linprog,'LargeScale')
    if ~isequal(options.linprog.LargeScale,'on')
        c = full(c);
        A = full(A);
        b = full(b);
        Aeq = full(Aeq);
        beq = full(beq);
    end
end

if ~options.warmstart
    x0 = [];
end

if options.savedebug
    ops = options.linprog;
    save linprogdebug c A b Aeq beq lb ub ops x0
end

solvertime = tic;
[x,fmin,flag,output,lambda] = linprog(c, A, b, Aeq, beq, lb, ub, x0,options.linprog);
solvertime = toc(solvertime);
problem = 0;

% Internal format for duals
if isempty(lambda)
    D_struc = [];
else
    D_struc = [lambda.eqlin;lambda.ineqlin];
end
if isempty(x)
    x = zeros(length(c),1);
end

% Check, currently not exhaustive...
if flag==0
    problem = 3;
elseif flag == -2 | flag==-5
    problem = 1;
elseif flag == -3
    problem = 2;
else
    if flag>0
        problem = 0;
    else
        if any((A*x-b)>sqrt(eps)) | any( abs(Aeq*x-beq)>sqrt(eps))
            problem = 1; % Likely to be infeasible
        else
            if c'*x<-1e10 % Likely unbounded
                problem = 2;
            else          % Probably convergence issues
                problem = 5;
            end
        end
    end
end

% Save all data sent to solver?
if options.savesolverinput
    solverinput.A = A;
    solverinput.b = b;
    solverinput.Aeq = Aeq;
    solverinput.beq = beq;
    solverinput.c = c;
    solverinput.options = options.linprog;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.fmin = fmin;
    solveroutput.flag = flag;
    solveroutput.output=output;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x(:),D_struc,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function a = r8cbb_add ( n1, n2, ml, mu, a, i, j, value )

%*****************************************************************************80
%
%% R8CBB_ADD adds a value to an entry of a R8CBB matrix.
%
%  Discussion:
%
%    The R8CBB storage format is for a compressed border banded matrix.
%    Such a matrix has the logical form:
%
%      A1 | A2
%      ---+---
%      A3 | A4
%
%    with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
%    are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
%    respectively.
%
%    The R8CBB format is the same as the DBB format, except that the banded
%    matrix A1 is stored in compressed band form rather than standard
%    banded form.  In other words, we do not include the extra room
%    set aside for fill in during pivoting.
%
%    A should be defined as a vector.  The user must then store
%    the entries of the four blocks of the matrix into the vector A.
%    Each block is stored by columns.
%
%    A1, the banded portion of the matrix, is stored in
%    the first (ML+MU+1)*N1 entries of A, using the obvious variant
%    of the LINPACK general band format.
%
%    The following formulas should be used to determine how to store
%    the entry corresponding to row I and column J in the original matrix:
%
%    Entries of A1:
%
%      1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
%
%      Store the I, J entry into location
%      (I-J+MU+1)+(J-1)*(ML+MU+1).
%
%    Entries of A2:
%
%      1 <= I <= N1, N1+1 <= J <= N1+N2.
%
%      Store the I, J entry into location
%      (ML+MU+1)*N1+(J-N1-1)*N1+I.
%
%    Entries of A3:
%
%      N1+1 <= I <= N1+N2, 1 <= J <= n1.
%
%      Store the I, J entry into location
%      (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%
%    Entries of A4:
%
%      N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
%
%      Store the I, J entry into location
%      (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%      (same formula used for A3).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N1, N2, the order of the banded and dense blocks.
%    N1 and N2 must be nonnegative, and at least one must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than N1-1.
%
%    Input, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the R8CBB matrix.
%
%    Input, integer I, J, the indices of the entry to be incremented.
%
%    Input, real VALUE, the value to be added to the (I,J) entry.
%
%    Output, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the modified R8CBB matrix.
%
  if ( value == 0.0 )
    return;
  end
%
%   Check for I or J out of bounds.
%
  if ( i <= 0 .or. n1+n2 < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8CBB_ADD - Fatal error!\n' );
    fprintf ( 1, '  Illegal input value of row index I = %d/n', i );
    error ( 'R8CBB_ADD - Fatal error!' );
  end

  if ( j <= 0 | n1+n2 < j )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8CBB_ADD - Fatal error!\n' );
    fprintf ( 1, '  Illegal input value of column index J = %d/n',j );
    error ( 'R8CBB_ADD - Fatal error!' );
  end
%
%  The A1 block of the matrix.
%
%  Check for out of band problems.
%
  if ( i <= n1 & j <= n1 )
    if ( mu < (j-i) | ml < (i-j) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8CBB_ADD - Warning!\n' );
      fprintf ( 1, '  Unable to add to entry (%d,%d).\n', i, j );
      return
    else
      ij = (i-j+mu+1)+(j-1)*(ml+mu+1);
    end
%
%  The A2 block of the matrix:
%
  elseif ( i <= n1 & n1 < j )
    ij = (ml+mu+1)*n1+(j-n1-1)*n1 + i;
%
%  The A3 and A4 blocks of the matrix.
%
  elseif ( n1 < i )
    ij = (ml+mu+1)*n1+n2*n1+(j-1)*n2 + (i-n1);
  end

  a(ij) = a(ij) + value;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pass = test_oldschool

%% Test the 'old-school' linop syntaxes for ATAP compatibility
% Toby Driscoll, 3 March 2014

d = domain(-1,1);

pref = cheboppref();
pref.discretization = @chebcolloc2;

%% Just test whether these will execute.
pass(1) = doesNotCrash( @() diff(d) );
pass(2) = doesNotCrash( @() eye(d) );
pass(3) = doesNotCrash( @() diag(chebfun('x',d([1 end])), d) );
pass(4) = doesNotCrash( @() cumsum(d) );
pass(5) = doesNotCrash( @() sum(d) );
pass(6) = doesNotCrash( @() feval(d, 0, 'left') );
pass(7) = doesNotCrash( @() zeros(d) );

%% Test the feval-style instantiation syntax.
D = diff(d);
warnState = warning('off', 'CHEBFUN:LINOP:feval:deprecated');
pass(8) = norm( D(8) - matrix(D,8,pref) ) < 2e-14;
warning(warnState);

%% Test boundary condition syntax
A = D^2;
[dummy1, dummy2, dummy3, A0] = matrix(A,10,pref);   % version with no BCs
[z,e,s,dif] = linop.primitiveFunctionals(d([1 end]));
A = addConstraint(A,e(d(1)),0);
A = addConstraint(A,e(d(end))*D,0);
A1 = matrix(A,10,pref);   % first two rows hold BCs
correct = cell2mat(A0);
correct([1 end],:) = A1(1:2,:);  % classic row replacement
warnState = warning('off', 'CHEBFUN:LINOP:feval:deprecated');
Aold = feval(A,12,'oldschool');
warning(warnState);
pass(9) = norm( correct - Aold ) < 2e-14;

end

function pass = doesNotCrash(fn)
try
    fn();
    pass = true;
catch ME
    pass = false;
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [elec] = read_polhemus_pos(filename)

% READ_POLHEMUS_POS reads electrode positions measured with the Polhemus tracker in
% one of the EEG labs at the DCCN. The software used with the Polhemus is from CTF.
%
% Use as:
%   [elec] = read_polhemus_pos(filename)
%
% This returns an electrode structure with
%   elec.label     cell-array with electrode labels (strings)
%   elec.pnt       position of each electrode

% Copyright (C) 2004-2019, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

fid = fopen_or_error(filename, 'rt');
line = fgetl(fid);
Nchan = str2double(line);

for i=1:Nchan
  line = fgetl(fid);
  [t, r] = strtok(line);
  elec.label{i} = char(t);
  elec.pnt(i,:) = sscanf(r, '%f')';
end
elec.label = elec.label(:);

try
  % read the fiducials
  line = fgetl(fid);
  [t, r] = strtok(line);
  fiducial.label{1} = char(t);
  fiducial.pnt(1,:) = sscanf(r, '%f')';
  line = fgetl(fid);
  [t, r] = strtok(line);
  fiducial.label{2} = char(t);
  fiducial.pnt(2,:) = sscanf(r, '%f')';
  line = fgetl(fid);
  [t, r] = strtok(line);
  fiducial.label{3} = char(t);
  fiducial.pnt(3,:) = sscanf(r, '%f')';
  % add the fiducials to the electrode array
  elec.label = cat(1, elec.label(:), fiducial.label(:));
  elec.pnt   = cat(1, elec.pnt, fiducial.pnt);
catch
  % do nothing
end

fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function demhint(nin, nhidden, nout)
%DEMHINT Demonstration of Hinton diagram for 2-layer feed-forward network.
%
%	Description
%
%	DEMHINT plots a Hinton diagram for a 2-layer feedforward network with
%	5 inputs, 4 hidden units and 3 outputs. The weight vector is chosen
%	from a Gaussian distribution as described under MLP.
%
%	DEMHINT(NIN, NHIDDEN, NOUT) allows the user to specify the number of
%	inputs, hidden units and outputs.
%
%	See also
%	HINTON, HINTMAT, MLP, MLPPAK, MLPUNPAK
%

%	Copyright (c) Ian T Nabney (1996-2001)

if nargin < 1 nin = 5; end
if nargin < 2 nhidden = 7; end
if nargin < 3 nout = 3; end

% Fix the seed for reproducible results
randn('state', 42);
clc
disp('This demonstration illustrates the plotting of Hinton diagrams')
disp('for Multi-Layer Perceptron networks.')
disp(' ')
disp('Press any key to continue.')
pause
net = mlp(nin, nhidden, nout, 'linear');

[h1, h2] = mlphint(net);
clc
disp('The MLP has been created with')
disp(['    ' int2str(nin) ' inputs'])
disp(['    ' int2str(nhidden) ' hidden units'])
disp(['    ' int2str(nout) ' outputs'])
disp(' ')
disp('One figure is produced for each layer of weights.')
disp('For each layer the fan-in weights are arranged in rows for each unit.')
disp('The bias weight is separated from the rest by a red vertical line.')
disp('The area of each box is proportional to the weight value: positive')
disp('values are white, and negative are black.')
disp(' ')
disp('Press any key to exit.');
pause;
delete(h1);
delete(h2);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef logistic_regression
% This file defines logistic regression (binary classifier) problem class
%
% Inputs:
%       x_train     train data matrix of x of size dxn.
%       y_train     train data vector of y of size 1xn.
%       x_test      test data matrix of x of size dxn.
%       y_test      test data vector of y of size 1xn.
%       lambda      l2-regularized parameter.
% Output:
%       Problem     problem instance.
%
%
% The problem of interest is defined as
%
%           min f(w) = 1/n * sum_i^n f_i(w),
%           where
%           f_i(w) = log(1 + exp(-y_i' .* (w'*x_i))) + lambda/2 * w^2.
%
% "w" is the model parameter of size d vector.
%
%
% This file is part of GDLibrary and SGDLibrary.
%
% Created by H.Kasai on Feb. 17, 2016
% Modified by H.Kasai on March 23, 2018

    properties
        name;
        dim;
        samples;
        lambda;
        classes;
        hessain_w_independent;
        d;
        n_train;
        n_test;
        x_train;
        y_train;
        x_test;
        y_test;
        x_norm;
        x;
    end

    methods
        function obj = logistic_regression(x_train, y_train, x_test, y_test, varargin)

            obj.x_train = x_train;
            obj.y_train = y_train;
            obj.x_test = x_test;
            obj.y_test = y_test;

            if nargin < 5
                obj.lambda = 0.1;
            else
                obj.lambda = varargin{1};
            end

            obj.d = size(obj.x_train, 1);
            obj.n_train = length(y_train);
            obj.n_test = length(y_test);
            obj.name = 'logistic_regression';
            obj.dim = obj.d;
            obj.samples = obj.n_train;
            obj.lambda = obj.lambda;
            obj.classes = 2;
            obj.hessain_w_independent = false;
            obj.x_norm = sum(obj.x_train.^2,1);
            obj.x = obj.x_train;
        end

        function f = cost(obj, w)
            %f_old = sum(log(1+exp(-y_train.*(w'*x_train)))/n_train,2)+ lambda * (w'*w) / 2; %is replaced below
            % becasuse log(sigmoid(a)) = log(1/(1+exp(-a))) = log1 - log(1+exp(-a)) = -log(1+exp(-a)).
            %f = -sum(log(sigmoid(y_train.*(w'*x_train)))/n_train,2) + lambda * (w'*w) / 2;

    %         Commented out below due to avoid '-Inf' values of g by HK on 2017/12/5
    %         f = -sum(log(sigmoid(y_train.*(w'*x_train))),2)/n_train + lambda * (w'*w) / 2;

    %         The above is replaced with below by HK on 2017/12/5
            sigmod_result = sigmoid(obj.y_train.*(w'*obj.x_train));
            sigmod_result = sigmod_result + (sigmod_result<eps).*eps;
            f = -sum(log(sigmod_result),2)/obj.n_train + obj.lambda * (w'*w) / 2;
        end

        function f = cost_batch(obj, w, indices)

            f = -sum(log(sigmoid(obj.y_train(indices).*(w'*obj.x_train(:,indices))))/obj.n_train,2) + obj.lambda * (w'*w) / 2;

        end

        function g = grad(obj, w, indices)

            % org code
            %g = -sum(((ones(obj.d,1) * (sigmoid(-obj.y_train(indices).*(w'*obj.x_train(:,indices))) .* obj.y_train(indices))) .* obj.x_train(:,indices))/length(indices),2) + obj.lambda * w;
            %
            % (log(1+exp(-y_train.*(w'*x_train)))' = -y_train.*x_train * (exp(-y_train.*(w'*x_train))/(1+exp(-y_train.*(w'*x_train))))
            %   = -y_train.*x_train * (1 - sigmoid(y_train.*(w'*x_train)))
            % or
            % (log(sigmoid(y_train.*(w'*x_train)))' = y_train.*x_train * 1/sigmoid(y_train.*(w'*x_train)) * (sigmoid(y_train.*(w'*x_train)))'
            %   = y_train.*x_train * 1/sigmoid(y_train.*(w'*x_train)) * (- sigmoid(y_train.*(w'*x_train)) * (1 - sigmoid(y_train.*(w'*x_train))))
            %   = -y_train.*x_train * (1 - sigmoid(y_train.*(w'*x_train)))

    %         Replace the commented-out lines with below (although it is a bit slower) due to avoid 'NAN' values of g by HK on 2017/12/5
            g = -sum(ones(obj.d,1) * obj.y_train(indices).*obj.x_train(:,indices) * (ones(1,length(indices))-sigmoid(obj.y_train(indices).*(w'*obj.x_train(:,indices))))',2)/length(indices)+ obj.lambda * w;

    %         Commented out below due to avoid 'NAN' values of g by HK on 2017/12/5
%             e = exp(-1*obj.y_train(indices)'.*(obj.x_train(:,indices)'*w));
%             s = e./(1+e);
%             g = -(1/length(indices))*((s.*obj.y_train(indices)')'*obj.x_train(:,indices)')';
%             g = full(g) + obj.lambda * w;

        end

        function g = full_grad(obj, w)

            %g = -sum(ones(d,1) * y_train.*x_train * (ones(1,n_train)-sigmoid(y_train.*(w'*x_train)))',2)/n_train+ lambda * w;
            g = grad(obj, w, 1:obj.n_train);
        end

        function g = ind_grad(obj, w, indices)

            g = -ones(obj.d,1) * obj.y_train(indices).*obj.x_train(:,indices) * diag(ones(1,length(indices))-sigmoid(obj.y_train(indices).*(w'*obj.x_train(:,indices))))+ obj.lambda* repmat(w, [1 length(indices)]);

        end

        function h = hess(obj, w, indices)

            %org code
            %temp = exp(-1*(y_train(indices)').*(x_train(:,indices)'*w));
            %b = temp ./ (1+temp);
            %h = 1/length(indices)*x_train(:,indices)*(diag(b-b.^2)*(x_train(:,indices)'))+lambda*eye(d);

            sigm_val = sigmoid(obj.y_train(indices).*(w'*obj.x_train(:,indices)));
            c = sigm_val .* (ones(1,length(indices))-sigm_val);
            h = 1/length(indices)* obj.x_train(:,indices) * diag(obj.y_train(indices).^2 .* c) * obj.x_train(:,indices)'+obj.lambda*eye(obj.d);

        end

        function h = full_hess(obj, w)

            h = hess(obj, w, 1:obj.n_train);

        end

        function hv = hess_vec(obj, w, v, indices)

            sigm_val = sigmoid(obj.y_train(indices).*(w'*obj.x_train(:,indices)));
            c = sigm_val .* (ones(1,length(indices))-sigm_val);
            hv = 1/length(indices)* obj.x_train(:,indices) * diag(obj.y_train(indices).^2 .* c) * (obj.x_train(:,indices)' * v) +obj.lambda*v;

        end

        function p = prediction(obj, w)

            p = sigmoid(w' * obj.x_test);

            class1_idx = p>0.5;
            class2_idx = p<=0.5;
            p(class1_idx) = 1;
            p(class2_idx) = -1;

        end

        function a = accuracy(obj, y_pred)

            a = sum(y_pred == obj.y_test) / obj.n_test;

        end

        function w_opt = calc_solution(obj, maxiter, method)

            if nargin < 3
                method = 'lbfgs';
            end

            options.max_iter = maxiter;
            options.verbose = true;
            options.tol_optgap = 1.0e-24;
            options.tol_gnorm = 1.0e-16;
            options.step_alg = 'backtracking';

            if strcmp(method, 'sd')
                [w_opt,~] = sd(obj, options);
            elseif strcmp(method, 'cg')
                [w_opt,~] = ncg(obj, options);
            elseif strcmp(method, 'newton')
                options.sub_mode = 'INEXACT';
                options.step_alg = 'non-backtracking';
                [w_opt,~] = newton(obj, options);
            else
                options.step_alg = 'backtracking';
                options.mem_size = 5;
                [w_opt,~] = lbfgs(obj, options);
            end
        end

        %% for NIM
        function [labels, samples] = get_partial_samples(obj, indices)
            samples = obj.x_train(:,indices);
            labels  = obj.y_train(indices);
        end

        function [s] = phi_prime(obj, w, indices)
            e = exp(-1.0 * obj.y_train(indices)' .* (obj.x_train(:,indices)'*w));
            s = e ./ (1.0+e);
        end

        function [ss] = phi_double_prime(obj, w, indices)
            e = exp(-1.0 * obj.y_train(indices)' .* (obj.x_train(:,indices)'*w));
            s = e ./ (1.0+e);
            ss = s .* (1.0 - s);
        end

        %% for Sub-sampled Newton
        function h = diag_based_hess(obj, w, indices, square_hess_diag)
            X = obj.x_train(:,indices)';
            h = X' * diag(square_hess_diag) * X / length(indices) + obj.lambda * eye(obj.d);
        end

        function square_hess_diag = calc_square_hess_diag(obj, w, indices)
            %hess_diag = 1./(1+exp(Y.*(X*w)))./(1+exp(-Y.*(X*w)));
            Xw = obj.x_train(:,indices)'*w;
            y = obj.y_train(indices)';
            yXw = y .* Xw;
            square_hess_diag = 1./(1+exp(yXw))./(1+exp(-yXw));
        end

    end
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function visual_scalpPlot(SMT, CNT, varargin)
% Description:
%   Draw  scalp topographies for all selected intervals,separately for each each class.
%   Scalp topographies of each classe are plotted in one row, and shared the same color map
%   scaling in each classes.
%
% Example Code:
%    visual_scalpPlot(SMT,CNT, {'Ival' , [start interval : time increase parameter: end intercal]});
%
% Input:
%   visual_scalpPlot(SMT,CNT, <OPT>);
%   SMT: Data structrue (ex) Epoched data structure
%   CNT: Continuous data structure
%
% Option:
%      .Ival - Selecting the interested time interval depending on time increase parameter
%                 (e.g. {'Ival' [ [-2000 : 1000: 4000]])
%
% Return:
%   Scalp topographies
%
% See also:
%    opt_getMontage, opt_cellToStruct
%
% Ji Hoon, Jeong
% jh_jeong@korea.ac.kr

%%
MNT = opt_getMontage(SMT);
opt = opt_cellToStruct(varargin{:});
opt.Interval = abs(opt.Ival(1)-opt.Ival(2));
for i = 1: size(SMT.class,1)
    for seg = 1: size(opt.Ival, 2)-1
        %     figure()
        subplot(size(SMT.class,1),size(opt.Ival,2)-1,((i-1)*(size(opt.Ival, 2)-1))+seg)

        center = [0 0];
        theta = linspace(0,2*pi,360);

        x = cos(theta)+center(1);
        y = sin(theta)+center(2);

        oldUnit = get(gcf,'units');
        set(gcf,'units','normalized');

        H = struct('ax', gca);
        set(gcf,'CurrentAxes',H.ax);

        % ----------------------------------------------------------------------
        % nose plot
        nose = [1 1.2 1];
        nosi = [83 90 97]+1;
        H.nose = plot(nose.*x(nosi), nose.*y(nosi), 'k', 'linewidth', 2.4 );

        hold on;

        % ----------------------------------------------------------------------
        % ears plot
        earw = .08; earh = .3;
        H.ears(1) = plot(x*earw-1-earw, y*earh, 'k',  'linewidth', 2);
        H.ears(2) = plot(x*earw+1+earw, y*earh, 'k',  'linewidth', 2);
        hold on;

        % ----------------------------------------------------------------------
        % main circle plot
        H.main = plot(x,y, 'k', 'linewidth', 2.2);
        set(H.ax, 'xTick',[], 'yTick',[]);
        axis('xy', 'tight', 'equal', 'tight');
        hold on;
        %         if (i == 1 && seg == 1)
        %             ylabel([SMT.class{i,2},' class'])
        %             hold on;
        %         elseif (i == 2 && seg == 1)
        %             ylabel([SMT.class{i,2},' class'])
        %             hold on;
        %         elseif (i == 3 && seg == 1)
        %             ylabel([SMT.class{i,2},' class'])
        %             hold on;
        %         elseif (i == 4 && seg == 1)
        %             ylabel([SMT.class{i,2},' class'])
        %             hold on;
        %         end
        % ----------------------------------------------------------------------
        % Rendering the contourf
        xe_org = MNT.x';
        ye_org = MNT.y';
        avgSMT= prep_average(SMT);

        SMTintervalstart = find(avgSMT.ival == opt.Ival(seg));
        SMTintervalEnd = SMTintervalstart+opt.Interval/10;
        UpdatedSMT{seg} = avgSMT.x(SMTintervalstart:SMTintervalEnd,:,:);
        w{seg,i} = UpdatedSMT{seg}(:,i,:);
        w{seg,i} = squeeze(w{seg,i});
        inputx{seg,i}  = mean(w{seg,i},1);

        % w = w(:);
        % resolution = 101;
        resolution = 500;

        maxrad = max(1,max(max(abs(MNT.x)),max(abs(MNT.y))));

        xx = linspace(-maxrad, maxrad, resolution);
        yy = linspace(-maxrad, maxrad, resolution)';

        xe_add = cos(linspace(0,2*pi,resolution))'*maxrad;
        ye_add = sin(linspace(0,2*pi,resolution))'*maxrad;
        w_add = ones(length(xe_add),1)*mean(inputx{seg,i});

        xe = [xe_org; xe_add];
        ye = [ye_org; ye_add];
        inputx{seg,i} = [inputx{seg,i}'; w_add];

        % xe_add = cos(linspace(0,2*pi,resolution))';
        % ye_add = sin(linspace(0,2*pi,resolution))';
        %
        % xe = [xe;xe_add];
        % ye = [ye;ye_add];
        % w = [w; zeros(length(xe_add),1)];

        [xg,yg,zg] = griddata(xe, ye, inputx{seg,i}, xx, yy);
        contourf(xg, yg, zg, 50, 'LineStyle','none'); hold on;

        % ----------------------------------------------------------------------
        % disp electrodes
        for j = 1:size(xe_org,1)
            plot(xe_org(j), ye_org(j), 'k*'); hold on;
            set(0,'defaultfigurecolor',[1 1 1])

        end

        axis off;
        title({[SMT.class{i,2},' class'];['[' , num2str(opt.Ival(seg)), ' ~ ' , num2str(opt.Ival(seg+1)) , '] ms']})
        %         colorbar('vert');
    end

end

end

author: PatternRecognition

repo: OpenBMI


text:
function [exp,waschanged] = exp_replacerepeated(exp,rules)
% Apply substitution rules to some expression until it no longer changes.
% [Out-Expression,Changed] = exp_replacerepeated(In-Expression,Rule-List)
%
% Rules may contain blanks and placeholders.
%
% In:
%   In-Expression : some expression in which substitution shall take place
%   Rule-List     : cell array of substutition rules; may also be a struct, where field names are
%                   taken as symbol names
%
% Out:
%   Out-Expression     : In-Expression with substitutions performed
%   Changed            : whether a replacement has taken place
%
% Notes:
%   Impure expressions degrade into pure expressions when substitutions are performed on them (i.e.
%   they lose their value).
%
% See also:
%   exp_replaceall, exp_match
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-19

if ~exp_beginfun('symbolic') return; end

[exp,waschanged] = utl_replacerepeated(exp,rules);

exp_endfun;

repo: REST


text:
function candoit = canGetApproxHessian(problem)
% Checks whether an approximate Hessian can be computed for this problem.
%
% function candoit = canGetApproxHessian(problem)
%
% Returns true if an approximate Hessian of the cost function is provided
% in the given problem description, false otherwise.
% If a Hessian is defined but no approximate Hessian is defined explicitly,
% returns false.
%
% Even if this returns false, calls to getApproxHessian may succeed, as
% they will be redirected to getHessianFD. The latter simply requires
% availability of gradients in problem, and vector transports in problem.M.
%
% See also: canGetHessian getHessianFD

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 8, 2015.
% Contributors:
% Change log:

    candoit = isfield(problem, 'approxhess');

end

author: NicolasBoumal

repo: manopt


text:
function err_shape = error_depth(param_dist,xcn,xpn,R,T,X_shape,ind);

X_new = depth_compute(xcn,xpn,[param_dist],R,T);

N_pt_calib = size(xcn,2);

% UnNormalized shape extraction:

X_shape2 = X_new;
X_shape2 = X_shape2 - (X_shape2(:,1)*ones(1,N_pt_calib));

% map the second vector at [1;0;0]:

omu = -cross([1;0;0],X_shape2(:,2));
omu = acos((dot([1;0;0],X_shape2(:,2)))/norm(X_shape2(:,2)))*(omu / norm(omu));
Ru = rodrigues(omu);

X_shape2 = Ru* X_shape2;

omu2 = -cross([0;1;0],[0;X_shape2(2:3,ind)]);
omu2 = acos((dot([0;1;0],[0;X_shape2(2:3,ind)]))/norm([0;X_shape2(2:3,ind)]))*(omu2 / norm(omu2));
Ru2 = rodrigues(omu2);

X_shape2 = Ru2* X_shape2;

% Error:

err_shape = X_shape2(:,2:end) - X_shape(:,2:end);

err_shape = err_shape(:);

%err_depth = Z_new - Z_ref;

repo: ekfmonoslam


text:
function [SIG,BAK,OVL]= comp_fwseg_variant(cleanFile, enhancedFile);

% ----------------------------------------------------------------------
%      Frequency-variant fwSNRseg Objective Speech Quality Measure
%
%   This function implements the frequency-variant fwSNRseg measure [1]
%   (see also Chap. 10, Eq. 10.24)
%
%
%   Usage:  [sig,bak,ovl]=comp_fwseg_variant(cleanFile.wav, enhancedFile.wav)
%
%         cleanFile.wav - clean input file in .wav format
%         enhancedFile  - enhanced output file in .wav format
%         sig           - predicted rating [1-5] of speech distortion
%         bak           - predicted rating [1-5] of noise distortion
%         ovl           - predicted rating [1-5] of overall quality
%
%
%  Example call:  [s,b,o] =comp_fwseg_variant('sp04.wav','enhanced.wav')
%
%
%  References:
%     [1] S. R. Quackenbush, T. P. Barnwell, and M. A. Clements,
%	    Objective Measures of Speech Quality.  Prentice Hall
%	    Advanced Reference Series, Englewood Cliffs, NJ, 1988,
%	    ISBN: 0-13-629056-6.
%
%   Author: Philipos C. Loizou
%  (critical-band filtering routines were written by Bryan Pellom & John Hansen)
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
% ----------------------------------------------------------------------

if nargin~=2
    fprintf('USAGE: [sig,bak,ovl]=comp_fwseg_variant(cleanFile.wav, enhancedFile.wav)\n');
    fprintf('For more help, type: help comp_fwseg_variant\n\n');
    return;
end

[data1, Srate1]= audioread(cleanFile);
[data2, Srate2]= audioread(enhancedFile);
% if ( Srate1~= Srate2) | ( Nbits1~= Nbits2)
%     error( 'The two files do not match!\n');
% end

len= min( length( data1), length( data2));
data1= data1( 1: len)+eps;
data2= data2( 1: len)+eps;

wss_dist_matrix= fwseg( data1, data2,Srate1);
wss_dist=mean(wss_dist_matrix);

% initialize  coefficients obtained from multiple linear
% regression analysis
%
b_sig=[0.021,-0.028,0.088,-0.031,0.048,-0.049,0.065,0.009,0.011,0.033,...
    -0.040,-0.002,0.041,-0.007,0.033,0.018,-0.007,0.044,-0.001,0.021,...
    -0.002,0.017,-0.03,0.073,0.043];
b_ovl=[-0.003,-0.026,0.066,-0.036,0.038,-0.023,0.037,0.022,0.014,0.009,...
    -0.03,0.004,0.044,-0.005,0.017,0.018,-0.001,0.051,0.009,0.011,...
    0.011,-0.002,-0.021,0.043,0.031];
b_bak=[-0.03,-0.022,0.03,-0.048,0.034,0.002,0.006,0.037,0.017,-0.016,-0.008,...
    0.019,0.024,-0.002,0.01,0.03,-0.018,0.046,0.022,0.005,0.03,-0.028,...
    -0.028,0.019,0.005];

SIG=0.567+sum(b_sig.*wss_dist);
SIG=max(1,SIG); SIG=min(5, SIG); % limit values to [1, 5]

BAK=1.013+sum(b_bak.*wss_dist);
BAK=max(1,BAK); BAK=min(5, BAK); % limit values to [1, 5]

OVL=0.446+sum(b_ovl.*wss_dist);
OVL=max(1,OVL); OVL=min(5, OVL); % limit values to [1, 5]

% ----------------------------------------------------------------------

function distortion = fwseg(clean_speech, processed_speech,sample_rate)

% ----------------------------------------------------------------------
% Check the length of the clean and processed speech.  Must be the same.
% ----------------------------------------------------------------------

clean_length      = length(clean_speech);
processed_length  = length(processed_speech);

if (clean_length ~= processed_length)
  disp('Error: Files  must have same length.');
  return
end

% ----------------------------------------------------------------------
% Global Variables
% ----------------------------------------------------------------------

winlength   = round(30*sample_rate/1000); 	   % window length in samples
skiprate    = floor(winlength/4);		   % window skip in samples
max_freq    = sample_rate/2;	   % maximum bandwidth
num_crit    = 25;		   % number of critical bands

n_fft       = 2^nextpow2(2*winlength);
n_fftby2    = n_fft/2;		   % FFT size/2

% ----------------------------------------------------------------------
% Critical Band Filter Definitions (Center Frequency and Bandwidths in Hz)
% ----------------------------------------------------------------------

cent_freq(1)  = 50.0000;   bandwidth(1)  = 70.0000;
cent_freq(2)  = 120.000;   bandwidth(2)  = 70.0000;
cent_freq(3)  = 190.000;   bandwidth(3)  = 70.0000;
cent_freq(4)  = 260.000;   bandwidth(4)  = 70.0000;
cent_freq(5)  = 330.000;   bandwidth(5)  = 70.0000;
cent_freq(6)  = 400.000;   bandwidth(6)  = 70.0000;
cent_freq(7)  = 470.000;   bandwidth(7)  = 70.0000;
cent_freq(8)  = 540.000;   bandwidth(8)  = 77.3724;
cent_freq(9)  = 617.372;   bandwidth(9)  = 86.0056;
cent_freq(10) = 703.378;   bandwidth(10) = 95.3398;
cent_freq(11) = 798.717;   bandwidth(11) = 105.411;
cent_freq(12) = 904.128;   bandwidth(12) = 116.256;
cent_freq(13) = 1020.38;   bandwidth(13) = 127.914;
cent_freq(14) = 1148.30;   bandwidth(14) = 140.423;
cent_freq(15) = 1288.72;   bandwidth(15) = 153.823;
cent_freq(16) = 1442.54;   bandwidth(16) = 168.154;
cent_freq(17) = 1610.70;   bandwidth(17) = 183.457;
cent_freq(18) = 1794.16;   bandwidth(18) = 199.776;
cent_freq(19) = 1993.93;   bandwidth(19) = 217.153;
cent_freq(20) = 2211.08;   bandwidth(20) = 235.631;
cent_freq(21) = 2446.71;   bandwidth(21) = 255.255;
cent_freq(22) = 2701.97;   bandwidth(22) = 276.072;
cent_freq(23) = 2978.04;   bandwidth(23) = 298.126;
cent_freq(24) = 3276.17;   bandwidth(24) = 321.465;
cent_freq(25) = 3597.63;   bandwidth(25) = 346.136;

bw_min      = bandwidth (1);	   % minimum critical bandwidth

% ----------------------------------------------------------------------
% Set up the critical band filters.  Note here that Gaussianly shaped
% filters are used.  Also, the sum of the filter weights are equivalent
% for each critical band filter.  Filter less than -30 dB and set to
% zero.
% ----------------------------------------------------------------------

min_factor = exp (-30.0 / (2.0 * 2.303));       % -30 dB point of filter

for i = 1:num_crit
  f0 = (cent_freq (i) / max_freq) * (n_fftby2);
  all_f0(i) = floor(f0);
  bw = (bandwidth (i) / max_freq) * (n_fftby2);
  norm_factor = log(bw_min) - log(bandwidth(i));
  j = 0:1:n_fftby2-1;
  crit_filter(i,:) = exp (-11 *(((j - floor(f0)) ./bw).^2) + norm_factor);
  crit_filter(i,:) = crit_filter(i,:).*(crit_filter(i,:) > min_factor);
end

% ----------------------------------------------------------------------
% For each frame of input speech, calculate the Weighted Spectral
% Slope Measure
% ----------------------------------------------------------------------

num_frames = floor(clean_length/skiprate-(winlength/skiprate)); % number of frames
start      = 1;					% starting sample
window     = 0.5*(1 - cos(2*pi*(1:winlength)'/(winlength+1)));

distortion=zeros(num_frames,num_crit);
for frame_count = 1:num_frames

   % ----------------------------------------------------------
   % (1) Get the Frames for the test and reference speech.
   %     Multiply by Hanning Window.
   % ----------------------------------------------------------

   clean_frame = clean_speech(start:start+winlength-1);
   processed_frame = processed_speech(start:start+winlength-1);
   clean_frame = clean_frame.*window;
   processed_frame = processed_frame.*window;

   % ----------------------------------------------------------
   % (2) Compute the magnitude Spectrum of Clean and Processed
   % ----------------------------------------------------------

       clean_spec     = abs(fft(clean_frame,n_fft));
       processed_spec = abs(fft(processed_frame,n_fft));

       % normalize so that spectra have unit area ----
        clean_spec=clean_spec/sum(clean_spec(1:n_fftby2));
        processed_spec=processed_spec/sum(processed_spec(1:n_fftby2));

   % ----------------------------------------------------------
   % (3) Compute Filterbank Output Energies (in dB scale)
   % ----------------------------------------------------------

   clean_energy=zeros(1,num_crit);
   processed_energy=zeros(1,num_crit);
   error_energy=zeros(1,num_crit);

   for i = 1:num_crit
      clean_energy(i) = sum(clean_spec(1:n_fftby2) ...
		            .*crit_filter(i,:)');
      processed_energy(i) = sum(processed_spec(1:n_fftby2) ...
			        .*crit_filter(i,:)');
      error_energy(i)=max((clean_energy(i)-processed_energy(i))^2,eps);
   end

   SNRlog=10*log10((clean_energy.^2)./error_energy);

   distortion(frame_count,:)=min(max(SNRlog,-10),35);

   start = start + skiprate;

end

author: anicolson

repo: DeepXi

save_path: github-repos/MATLAB/anicolson-DeepXi


text:
function showsolutionP2(node,elem,u)

[node,elem] = uniformrefine(node,elem);
showsolution(node,elem,u);

repo: Aero-matlab


text:
function TDOAshell

wave = wavread('sample.wav');
wave = wave(:,1);
scale = 0.8/max(wave);
wave = scale*wave;
Trials = 10;
Radius = 50;
N = 8;
Theta = linspace(0,2*pi,N+1);
X = Radius * cos(Theta(1:end-1));
Y = Radius * sin(Theta(1:end-1));
Z = [1:N];
Z = (-1).^Z;
Z = 5*Z+5;
Sen_position = [X.',Y.',Z.'];
Sen_position = [Sen_position];
True_position = zeros(Trials, 3);
Est_position = zeros(Trials,3);

% Generate position of source

for i=1:Trials
    r = rand(1)*50;
    t = rand(1)*2*pi;
    x = r*cos(t);
    y = r*sin(t);
    z = rand(1)*20;
    True_position(i,1) = x;
    True_position(i,2) = y;
    True_position(i,3) = z;
end

% Generate distances

Distances = zeros(Trials,8);
for i=1:Trials
    for j=1:8
        x1 = True_position(i,1);
        y1 = True_position(i,2);
        z1 = True_position(i,3);
        x2 = Sen_position(j,1);
        y2 = Sen_position(j,2);
        z2 = Sen_position(j,3);
        Distances(i,j) = sqrt((x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2);
    end
end

Distances;
TimeDelay = Distances./340.29;
Padding = TimeDelay*44100;

% Generate the signals
for i=1:Trials
   x = True_position(i,1);
   y = True_position(i,2);
   z = True_position(i,3);
   xstr = num2str(round(x));
   ystr = num2str(round(y));
   zstr = num2str(round(z));
   istr = num2str(i);
   name = strcat( 'Trial_', istr, '_', xstr, '_', ystr, '_', zstr, '_mdove.wav');
   mic1 = [zeros(round(Padding(i,1)),1) ; wave];
   mic2 = [zeros(round(Padding(i,2)),1) ; wave];
   mic3 = [zeros(round(Padding(i,3)),1) ; wave];
   mic4 = [zeros(round(Padding(i,4)),1) ; wave];
   mic5 = [zeros(round(Padding(i,5)),1) ; wave];
   mic6 = [zeros(round(Padding(i,6)),1) ; wave];
   mic7 = [zeros(round(Padding(i,7)),1) ; wave];
   mic8 = [zeros(round(Padding(i,8)),1) ; wave];
   l1 = length(mic1);
   l2 = length(mic2);
   l3 = length(mic3);
   l4 = length(mic4);
   l5 = length(mic5);
   l6 = length(mic6);
   l7 = length(mic7);
   l8 = length(mic8);
   lenvec = [l1 l2 l3 l4 l5 l6 l7 l8];
   m = max(lenvec);
   c = [m-l1, m-l2, m-l3, m-l4, m-l5,m-l6, m-l7, m-l8];
   mic1 = [mic1; zeros(c(1),1)];
   mic2 = [mic2; zeros(c(2),1)];
   mic3 = [mic3; zeros(c(3),1)];
   mic4 = [mic4; zeros(c(4),1)];
   mic5 = [mic5; zeros(c(5),1)];
   mic6 = [mic6; zeros(c(6),1)];
   mic7 = [mic7; zeros(c(7),1)];
   mic8 = [mic8; zeros(c(8),1)];

   mic1 = mic1./Distances(i,1);
   mic2 = mic2./Distances(i,2);
   mic3 = mic3./Distances(i,3);
   mic4 = mic4./Distances(i,4);
   mic5 = mic5./Distances(i,5);
   mic6 = mic6./Distances(i,6);
   mic7 = mic7./Distances(i,7);
   mic8 = mic8./Distances(i,8);

multitrack = [mic1, mic2, mic3, mic4, mic5, mic6, mic7, mic8];
%   wavwrite(multitrack, 44100, name);

[x y z] = Locate(Sen_position, multitrack);
Est_position(i,1) = x;
Est_position(i,2) = y;
Est_position(i,3) = z;

end

plot(True_position(:,1),True_position(:,2),'bd',Est_position(:,1),Est_position(:,2),'r+','LineWidth',2);
legend('True Position','Estimated Position');
xlabel('X coordinate of target');
ylabel('Y coordinate of target');
title('TDOA Hyperbolic Localization');
axis([-50 50 -50 50]);

end

function [x y z] = Locate(Sen_position, multitrack)
% sensor index shift of 1 occurrs here

s = size(Sen_position);
len = s(1);
timedelayvec = zeros(len,1);
for i=1:len
   timedelayvec(i) = timedelayfunc(multitrack(:,1),multitrack(:,i));
end

timedelayvec;

Amat = zeros(len,1);
Bmat = zeros(len,1);
Cmat = zeros(len,1);
Dmat = zeros(len,1);
for i=3:len
    x1 = Sen_position(1,1);
    y1 = Sen_position(1,2);
    z1 = Sen_position(1,3);
    x2 = Sen_position(2,1);
    y2 = Sen_position(2,2);
    z2 = Sen_position(2,3);
    xi = Sen_position(i,1);
    yi = Sen_position(i,2);
    zi = Sen_position(i,3);
   Amat(i) = (1/(340.29*timedelayvec(i)))*(-2*x1+2*xi) - (1/(340.29*timedelayvec(2)))*(-2*x1+2*x2);
   Bmat(i) = (1/(340.29*timedelayvec(i)))*(-2*y1+2*yi) - (1/(340.29*timedelayvec(2)))*(-2*y1+2*y2);
   Cmat(i) = (1/(340.29*timedelayvec(i)))*(-2*z1+2*zi) - (1/(340.29*timedelayvec(2)))*(-2*z1+2*z2);
   Sum1 = (x1^2)+(y1^2)+(z1^2)-(xi^2)-(yi^2)-(zi^2);
   Sum2 = (x1^2)+(y1^2)+(z1^2)-(x2^2)-(y2^2)-(z2^2);
   Dmat(i) = 340.29*(timedelayvec(i) - timedelayvec(2)) + (1/(340.29*timedelayvec(i)))*Sum1 - (1/(340.29*timedelayvec(2)))*Sum2;
end

M = zeros(len,3);
D = zeros(len,1);
for i=1:len
    M(i,1) = Amat(i);
    M(i,2) = Bmat(i);
    M(i,3) = Cmat(i);
    D(i) = Dmat(i);
end

M = M(3:len,:);
D = D(3:len);

D = D.*-1;

Minv = pinv(M);
T = Minv*(D);
x = T(1);
y = T(2);
z = T(3);

end

function out = timedelayfunc(x,y)
% suppose sampling rate is 44100
% Let Tx be transit time for x
% Let Ty be transit time for y
% out is Ty - Tx

c = xcorr(x, y);
[C I] = max(c);
out = ((length(c)+1)/2 - I)/44100;

end

author: StevenJL


text:
%
% A variational approach to SPCP (Aravkin et al. 2014)
%
% RPCA | SPCP-max-QN  | Stable PCP-max solved by Quasi-Newton (Aravkin et al. 2014)
% process_video('RPCA', 'flip-SPCP-max-QN', 'dataset/demo.avi', 'output/demo_flip-SPCP-max-QN.avi');

alg_path_aux = fullfile(lrs_conf.rpca_path,'SPGL1');
addpath(genpath(alg_path_aux));

nFrames     = size(M,2);
lambda      = 1/sqrt(max(size(M,1),size(M,2)));
L0          = repmat(median(M,2), 1, nFrames);
S0          = M - L0;
epsilon     = 5e-3*norm(M,'fro'); % tolerance for fidelity to data

% Flip-Flop version pf SPCP-max solved by Quasi-Newton
opts = struct('sum',false,'L0',L0,'S0',S0,'max',true,...
  'tau0',3e5,'SPGL1_tol',1e-1,'tol',1e-3);
[L,S] = solver_RPCA_SPGL1(M,lambda,epsilon,[],opts);

rmpath(genpath(alg_path_aux));

author: andrewssobral

repo: lrslibrary


text:
function [u,p,info] = mgstokes(A,B,f,g,u,p,elem,freeDof,option)
%% MGSTOKES: multigrid methods for Stokes equations
%
% Created by Ming Wang and Jie Zhou on Mar,3,2013 for P2 elements. Add
% other elements by Jie and Long on Mar 24,2013.

t = cputime;

%% Size of systems
N = max(elem(:));
Ndof = length(f) + length(g);

%% Options
% Assign default values to unspecified parameters
if ~exist('option','var'),
    option = [];
end
option = mgoptions(option,Ndof);    % parameters
tol = option.tol; maxIt = option.solvermaxit;
printlevel = option.printlevel;
solver = option.solver;

%% Mesh on each level
N0 = 10;
HB = zeros(N,3);
level = 20;
NL(level+1) = N; % now NL(1:level) = 0;
elemi = cell(level,1);
elemi{level} = elem;
for j = level: -1 : 2
    [elemi{j-1},newHB] = uniformcoarsenred(elemi{j});  % try coasen red refinement
    if (isempty(newHB)) || (size(elemi{j-1},1)< 2*N0)
    % no nodes are removed or it reaches the coarsest level
        NL = NL(j:end);
        break;
    end
    NL(j) = NL(j+1) - size(newHB,1); % update NL(k)
    HB(NL(j)+1:NL(j+1),1:3) = newHB(:,1:3);
end
level = length(NL)-1;    % actual level
elemi = elemi(end-level+1:end);

%% Transfer operators between levels
Pro_u = cell(level,1);
Res_u = cell(level,1);
Pro_p = cell(level,1);
Res_p = cell(level,1);

switch(upper(option.elemType(1:end-2))) % transfer operator of velocity
    case 'P2'
        P2FreeDof = freeDof(1:end/2);   % again only one component
        for i = level:-1:2
            [ProU0,P2FreeDof] = nodeP2TansferUniform(elemi{i-1},elemi{i},P2FreeDof);
            Pro_u{i-1} = ProU0;
            Res_u{i}   = ProU0';
        end
    case 'CR'
        CRfreeDof = freeDof(1:end/2);
        for i=level:-1:2
            [ProU0,CRfreeDof] = transferCRred(elemi{i-1},elemi{i},CRfreeDof);
            Pro_u{i-1} = ProU0;
            Res_u{i}   = ProU0';
        end
    case 'P1B'
        for i=level:-1:2
            Pro_u{i-1} = speye(length(ufreeDof),length(freeEdge));
            Res_u{i} = Pro_u{i-1}';
        end
    case 'ISOP2'
        P2FreeDof = freeDof(1:end/2);  % again only one component
        for i = level:-1:2
            [ProU0,P2FreeDof] = nodeisoP2TansferUniform(elemi{i-1},elemi{i},P2FreeDof);
            Pro_u{i-1} = ProU0;
            Res_u{i}   = ProU0';
        end
end
switch upper(option.elemType(end-1:end))
    case  'P0'
        for k = level:-1:2
            Np(k-1)    = size(elemi{k-1},1);
            Pro_p{k-1} = repmat(speye(Np(k-1)),4,1);
            Res_p{k} = Pro_p{k-1}';
        end
    case  'P1'
       [Pro_p,Res_p] = transferoperator(HB,NL);
end

%% Matrices in each level
% size of u and p
Nu = zeros(level,1);
Np = zeros(level,1);
Nu(level) = size(A,1)/2; % one component of velocity
Np(level) = size(B,1);

Ai = cell(level,1);
Bi = cell(level,1);
bigAi = cell(level,1);
Ai{level} = A;
Bi{level} = B;
bigAi{level} = [Ai{level} Bi{level}'; Bi{level} sparse(Np(level),Np(level))];

for j = level:-1:2
    Ai{j-1} = blkdiag(Res_u{j},Res_u{j})*Ai{j}*blkdiag(Pro_u{j-1},Pro_u{j-1});  % Ac = Res*Af*Pro
    Bi{j-1} = Res_p{j}*Bi{j}*blkdiag(Pro_u{j-1},Pro_u{j-1});
    Nu(j-1) = size(Ai{j-1},1)/2;
    Np(j-1) = size(Bi{j-1},1);
    bigAi{j-1} = [Ai{j-1} Bi{j-1}'; Bi{j-1} sparse(Np(j-1),Np(j-1))];
end
Ndof = 2*Nu+Np;

%% matrix used for smoothing
if isfield(option,'smootherType')
smoother = option.smootherType;
else
smoother ='LSCDGS';
end
auxMat = cell(level,1);
switch(smoother)
    case 'LSCDGS'
        for k = 2:level
            Bt = (Bi{k})';
            BBt = Bi{k}*Bt;
            BABt = Bi{k}*Ai{k}*Bt;
            Su = tril(Ai{k});
            Sp = tril(BBt);
            Spt = triu(BBt);
            DSp = diag(BBt);
            auxMat{k} = struct('Bt',Bt,'BBt',BBt,'BABt',BABt,'Su',Su,...
                               'Spt',Spt,'Sp',Sp,'DSp',DSp);
        end
    case 'IUzawa'
        for k = 2:level
            Bt = (Bi{k})';
            DA = 2*diag(Ai{k});
            invDA = spdiags(1./DA,0,2*Nu(k),2*Nu(k));
            BinvDABt = Bi{k}*invDA*Bt;
            auxMat{k} = struct('Bt',Bt,'DA',DA,'BinvDABt',BinvDABt);
        end
end
Ai_IU = cell(level,1);
Si_IU = cell(level,1);
SSi_IU = cell(level,1);
Res_IU = cell(level,1);
Pro_IU = cell(level,1);
if strcmp(option.smootherbarSp,'VCYCLE')
%     option.smootherbarSp = 'VCYCLE';
    innerMGoption.solver = 'NO';
    innerMGoption.N0 = 10;
    for j = 2:level
        [tempvar,tempvar,Ai_IU{j},Si_IU{j},SSi_IU{j},Res_IU{j},Pro_IU{j}] = ...
                mg(auxMat{j}.BBt,ones(Np(j),1),elemi{j},innerMGoption);
    end
end

%% Solver
%  initial set up
bigF = [f; g-mean(g)];
bigu = [u; p];
bigr = bigF - bigAi{level}*bigu;
nb = norm(bigF);
err = zeros(maxIt,1);
err(1) = norm(bigr)/nb;
k = 1;

while (max(err(k)) > tol) && (k <= maxIt)
    k = k + 1;
    switch (solver)
        case 'VCYCLE'
            bigerru = vcycle(bigr);
        case 'WCYCLE'
            bigerru = wcycle(bigr);
    end
    bigu = bigu+bigerru;
    bigr = bigr - bigAi{level}*bigerru;
    % compute the relative error
    err(k) = norm(bigr)/nb;
    if printlevel >= 2
        fprintf('#dof: %8.0u, MG %8s iter: %2.0u, err = %8.4e\n',...
            Ndof(level), solver,k-1, err(k));
    end
end
err = err(1:k);
itStep = k-1;
u = bigu(1:2*Nu(level));
p = bigu(2*Nu(level)+1:end);

%% Output
if k > maxIt
    flag = 1;
else
    flag = 0;
end
time = cputime-t;
if printlevel >= 1
    fprintf('#dof: %6.0u,  #nnz: %6.0u, level: %2.0u  MG %6s iter: %2.0u,  err = %8.4e,  time = %4.2g s\n',...
             Ndof(level), nnz(bigAi{level}), level, solver, itStep, err(end), time);
end
if (flag == 1) && (printlevel>0)
   fprintf('NOTE: the iterative method does not converge! \n');
end
info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions vcycle, wcycle, smoothing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% subfunctions v-cycle
     function e = vcycle(r,J)
        if nargin < 2, J = level; end
        if J == 1 % solver in the coaresest grid
            e = zeros(size(r));
            e(1:end-1) = bigAi{J}(1:end-1,1:end-1)\r(1:end-1);
            e(2*Nu(J)+1:end) = e(2*Nu(J)+1:end)-mean(e(2*Nu(J)+1:end));
            return
        end
        ru = r(1:2*Nu(J)); rp = r(2*Nu(J)+1:end);

        % pre-smoothing
        [eu,ep] = smoothing(zeros(2*Nu(J),1),zeros(Np(J),1),ru,rp,J);

        % form residual and restrict onto coarse grid
        rru = ru-Ai{J}*eu-(Bi{J})'*ep;
        rrp = rp-Bi{J}*eu;
        ruc = reshape(Res_u{J}*reshape(rru,Nu(J),2),2*Nu(J-1),1);
        rpc = Res_p{J}*rrp;

        % coarse grid correction
        rc = [ruc;rpc];
        ec = vcycle(rc,J-1);

        % correction on the fine grid
        tempeu = reshape(Pro_u{J-1}*reshape(ec(1:2*Nu(J-1)),Nu(J-1),2),2*Nu(J),1);
        tempep = Pro_p{J-1}*ec(2*Nu(J-1)+1:end);
        eu = tempeu + eu;
        ep = tempep + ep;

        % post-smoothing
        [eu,ep] = smoothing(eu,ep,ru,rp,J);
        e = [eu;ep];

     end

 %% Wcycle MG
     function e = wcycle(r,J)
        if nargin < 2, J = level; end
        if J == 1 % solver in the coaresest grid
            e = zeros(size(r));
            e(1:end-1) = bigAi{J}(1:end-1,1:end-1)\r(1:end-1);
            e(2*Nu(J)+1:end) = e(2*Nu(J)+1:end)-mean(e(2*Nu(J)+1:end));
            return
        end
        ru = r(1:2*Nu(J));
        rp = r(2*Nu(J)+1:end);
%
        % pre-smoothing in the fine grid
        [eu,ep] = smoothing(zeros(2*Nu(J),1),zeros(Np(J),1),ru,rp,J);

        % form residual and restrict onto coarse grid
        rru = ru - Ai{J}*eu-(Bi{J})'*ep;
        rrp = rp - Bi{J}*eu;
        ruc = reshape(Res_u{J}*reshape(rru,Nu(J),2),2*Nu(J-1),1);
        rpc = Res_p{J}*rrp;

        % coarse grid correction
        rc = [ruc;rpc];
        ec = wcycle(rc,J-1);
        % once more for w-cycle
        ec = ec+ wcycle(rc - bigAi{J-1}*ec,J-1);

        % correction on the fine grid
        tempeu = reshape(Pro_u{J-1}*reshape(ec(1:2*Nu(J-1)),Nu(J-1),2),2*Nu(J),1);
        tempep = Pro_p{J-1}*ec(2*Nu(J-1)+1:end);
        eu = eu + tempeu;
        ep = ep + tempep;

        % post-smoothing in the fine grid
        [eu,ep] = smoothing(eu,ep,ru,rp,J);
        e = [eu;ep];
     end

 %% Smoothing
    function [u,p] = smoothing(u,p,f,g,J)
        if strcmp(smoother,'LSCDGS')
            [u,p] = StokesLSCDGS(u,p,f,g,Ai{J},Bi{J},auxMat{J},elemi{J},option,...
                                 Ai_IU{J},Si_IU{J},SSi_IU{J},Res_IU{J},Pro_IU{J});
        elseif strcmp(smoother,'BDDGS')
            % e = StokesBiDGStri(Ai{J}, Bi{J}, e, r, smoothStep, invDM{J},area{J});
        elseif strcmp(smoother,'iUzawa')
            [u,p] = StokesIUzawa(u,p,f,g,Ai{J},Bi{J},auxMat{J},elemi{J},option,...
                                 Ai_IU{J},Si_IU{J},SSi_IU{J},Res_IU{J},Pro_IU{J});
        end
    end
end


text:
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Relevant Vector Regression 1D Example %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%               1) Load 1D Regression Datasets               %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  (simple example) Generate Data from sine function
clear all; close all; clc;
nbSamples = 200;
epsilon   = 0.1;
x_limits  = [0, 100];

% Generate True function and data
X         = linspace(x_limits(1),x_limits(2),nbSamples)';
y_true    = sin(X*0.05);
y         = y_true + normrnd(0,epsilon,1,nbSamples)';

% Plot data
options             = [];
options.points_size = 10;
options.title       = 'noisy sinusoidal data';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1      = ml_plot_data([X(:),y(:)],options); hold on;

% Plot True function and Data
plot(X,y_true,'--k','LineWidth',2);
legend({'data','true function'})

%% (complex example) Generate Data from a sinc function
clear all; close all; clc;
% Set parameters for sinc function data
nbSamples = 200;
epsilon   = 0.2;
y_offset  = 0.5;
x_limits  = [-5, 5];

% Generate True function and data
X = linspace(x_limits(1),x_limits(2),nbSamples) ;
y_true = sinc(X) + y_offset ;
y = y_true + normrnd(0,epsilon,1,nbSamples);

% Plot data
options             = [];
options.points_size = 15;
options.title       = 'noisy sinc data';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1      = ml_plot_data([X(:),y(:)],options); hold on;

% Plot True function and Data
plot(X,y_true,'--k','LineWidth',2);
legend({'data','true function'})

% Transform Data for CV
X = X'; y = y';

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                    2)  RELEVANT VECTOR REGRESSION                     %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% RVR + RBF Kernel

clear rvr_options

%Set RVR OPTIONS%
rvr_options.useBias = true;
rvr_options.maxIts  = 100;

%Set Kernel OPTIONS%
rvr_options.kernel_ = 'gauss';
rvr_options.width   = 1;

% Train RVR Model
clear model
[~, model] = rvm_regressor(X,y,rvr_options,[]);

% Plot RVR function
ml_plot_rvr_function(X, y, model, rvr_options);

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   3) Do K-fold cross validation on hyper-parameters for RVR      %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% K-fold cross validation

Kfold = 10;
disp('Parameter grid search RVR');

%Set RVR OPTIONS%
rvr_options.useBias = true;
rvr_options.maxIts  = 100;

%Set Kernel OPTIONS%
rvr_options.kernel_ = 'gauss';

% Modify these according to your data!
rbf_vars = [0.1:0.1:1];

test  = cell(length(rbf_vars),1);
train = cell(length(rbf_vars),1);

for i=1:length(rbf_vars)
    disp(['[' num2str(i) '/' num2str(length(rbf_vars)) ']']);

    rvr_options.width       = rbf_vars(i);   %  radial basis function: exp(-gamma*|u-v|^2), gamma = 1/(2*sigma^2)

    f                       = @(X,y,model)rvm_regressor(X,y,rvr_options,model);
    [test_eval,train_eval]  = ml_kcv(X,y,Kfold,f,'regression');

    test{i}                 = test_eval;
    train{i}                = train_eval;
    disp(' ');
end

%% Get Statistics

[ stats ] = ml_get_cv_grid_states_regression(test,train);

% Plot Statistics

options             = [];
options.title       = 'RVR k-CV';
options.metrics     = {'nmse'};     % <- you can add many other metrics, see list in next cell box
options.para_name   = 'variance rbf';

[handle,handle_test,handle_train] = ml_plot_cv_grid_states_regression(stats,rbf_vars,options);

author: epfl-lasa


text:
function boxes1 = cascade_test(model, prec, testset, year, suffix)

% boxes1 = cascade_test(cls, model, testset, year, suffix)
% Compute bounding boxes in a test set.
% boxes1 are detection windows and scores.
%
% prec specifies the target level of precision that wish to achieve.
% if you let prec = [], a threshold tuned to reach the prec-equals-recall
% point will be selected.

% Now we also save the locations of each filter for rescoring
% parts1 gives the locations for the detections in boxes1
% (these are saved in the cache file, but not returned by the function)

setVOCyear = year;
globals;
pascal_init;

cls = model.class;
ids = textread(sprintf(VOCopts.imgsetpath, testset), '%s');

[prec, recall, thresh] = cascade_thresh(model, year, prec);
model.thresh = thresh;
pca = 5;
model = cascade_model(model, model.year, pca, model.thresh);

% run detector in each image
try
  load([cachedir cls '_boxes_' testset '_' suffix]);
catch
  times = zeros(length(ids), 2);
  % parfor gets confused if we use VOCopts
  opts = VOCopts;
  % parallel implementation disabled for single-threaded tests
  %parfor i = 1:length(ids);
  for i = 1:length(ids);
    if strcmp('inriaperson', cls)
      % INRIA uses a mixutre of PNGs and JPGs, so we need to use the annotation
      % to locate the image.  The annotation is not generally available for PASCAL
      % test data (e.g., 2009 test), so this method can fail for PASCAL.
      rec = PASreadrecord(sprintf(opts.annopath, ids{i}));
      im = imread([opts.datadir rec.imgname]);
    else
      im = imread(sprintf(opts.imgpath, ids{i}));
    end
    th = tic();
    pyra = featpyramid(im, model);
    time_feat = toc(th);

    th = tic();
    [dets, boxes] = cascade_detect(pyra, model, model.thresh);
    time_det = toc(th);

    if ~isempty(boxes)
      [dets boxes] = clipboxes(im, dets, boxes);
      I = nms(dets, 0.5);
      boxes1{i} = dets(I,[1:4 end]);
      parts1{i} = boxes(I,:);
    else
      boxes1{i} = [];
      parts1{i} = [];
    end
    times(i,:) = [time_det time_feat];
    fprintf('%s: testing: %s %s, %d/%d (avg time %.3f)\n', cls, testset, year, ...
            i, length(ids), mean(times(1:i,1)));
  end
  save([cachedir cls '_boxes_' testset '_' suffix], ...
       'boxes1', 'parts1', 'times');
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [x,varargout] = ZeroToOne(x,varargin)

%ZeroToOne - Normalize values in [0,1].
%
%  USAGE
%
%    [y,b0,b1...] = ZeroToOne(x,a0,a1,...)
%
%    x              array to normalize (vector OR matrix)
%    a0...          additional inputs to transform using the same
%                   scale as x

% Copyright (C) 2008-2011 by Michaël Zugaro, modified GG 2015
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

if nargin < 1,
	error('Incorrect number of parameters (type ''help ZeroToOne'' for details).');
end

if nargin ~= nargout,
	error('Different numbers of input and output parameters (type ''help ZeroToOne'' for details).');
end

if isvector(x)
    m = min(x);
    M = max(x);
    x = (x-m)/(M-m);
elseif ismatrix(x)
    m = min(min(x));
    M = max(max(x));
    x = (x-m)./(M-m);
end

for i = 1:nargin-1,
	varargout{i} = (varargin{i}-m)/(M-m);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [dat] = read_brainvision_eeg(filename, hdr, begsample, endsample);

% READ_BRAINVISION_EEG reads raw data from an EEG file
% and returns it as a Nchans x Nsamples matrix
%
% Use as
%   dat = read_brainvision_eeg(filename, hdr, begsample, endsample)
% where the header should be first read using read_brainvision_vhdr
%
% See also READ_BRAINVISION_VHDR, READ_BRAINVISION_VMRK

% Copyright (C) 2003, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_brainvision_eeg.m 945 2010-04-21 17:41:20Z roboos $

if strcmpi(hdr.DataFormat, 'binary') && strcmpi(hdr.DataOrientation, 'multiplexed') && strcmpi(hdr.BinaryFormat, 'int_16')
  fid = fopen(filename, 'rb', 'ieee-le');
  fseek(fid, hdr.NumberOfChannels*2*(begsample-1), 'cof');
  [dat, siz] = fread(fid, [hdr.NumberOfChannels, (endsample-begsample+1)], 'int16');
  fclose(fid);
  % compute real microvolts using the calibration factor (resolution)
  res = sparse(diag(hdr.resolution));
  dat = res * dat;

elseif strcmpi(hdr.DataFormat, 'binary') && strcmpi(hdr.DataOrientation, 'multiplexed') && strcmpi(hdr.BinaryFormat, 'int_32')
  fid = fopen(filename, 'rb', 'ieee-le');
  fseek(fid, hdr.NumberOfChannels*4*(begsample-1), 'cof');
  [dat, siz] = fread(fid, [hdr.NumberOfChannels, (endsample-begsample+1)], 'int32');
  fclose(fid);
  % compute real microvolts using the calibration factor (resolution)
  res = sparse(diag(hdr.resolution));
  dat = res * dat;

elseif strcmpi(hdr.DataFormat, 'binary') && strcmpi(hdr.DataOrientation, 'multiplexed') && strcmpi(hdr.BinaryFormat, 'ieee_float_32')
  fid = fopen(filename, 'rb', 'ieee-le');
  fseek(fid, hdr.NumberOfChannels*4*(begsample-1), 'cof');
  [dat, siz] = fread(fid, [hdr.NumberOfChannels, (endsample-begsample+1)], 'float32');
  fclose(fid);

elseif strcmpi(hdr.DataFormat, 'binary') && strcmpi(hdr.DataOrientation, 'vectorized') && strcmpi(hdr.BinaryFormat, 'ieee_float_32')
  fid = fopen(filename, 'rb', 'ieee-le');
  fseek(fid, 0, 'eof');
  hdr.nSamples = ftell(fid)/(4*hdr.NumberOfChannels);
  fseek(fid, 0, 'bof');
  numsamples = (endsample-begsample+1);
  for chan=1:hdr.NumberOfChannels
    fseek(fid, (begsample-1)*4, 'cof');                 % skip the first N samples
    [tmp, siz] = fread(fid, numsamples, 'float32');     % read these samples
    fseek(fid, (hdr.nSamples-endsample)*4, 'cof');      % skip the last M samples
    dat(chan,:) = tmp(:)';
  end
  fclose(fid);

elseif strcmpi(hdr.DataFormat, 'ascii') && strcmpi(hdr.DataOrientation, 'multiplexed')
  fid = fopen(filename, 'rt');
  for line=1:(begsample-1)
    % read first lines and discard the data in them
    str = fgets(fid);
  end
  dat = zeros(endsample-begsample+1, hdr.NumberOfChannels);
  for line=1:(endsample-begsample+1)
    str = fgets(fid);           % read a single line with Nchan samples
    str(find(str==',')) = '.';      % replace comma with point
    dat(line,:) = str2num(str);
  end
  fclose(fid);
  % transpose the data
  dat = dat';

elseif strcmpi(hdr.DataFormat, 'ascii') && strcmpi(hdr.DataOrientation, 'vectorized')
  % this is a very inefficient fileformat to read data from, since it requires to
  % read in all the samples of each channel and then select only the samples of interest
  fid = fopen(filename, 'rt');
  dat = zeros(hdr.NumberOfChannels, endsample-begsample+1);
  for chan=1:hdr.NumberOfChannels
    % this is very slow, so better give some feedback to indicate that something is happening
    fprintf('reading channel %d from ascii file to get data from sample %d to %d\n', chan, begsample, endsample);

    str = fgets(fid);             % read all samples of a single channel
    str(find(str==',')) = '.';      % replace comma with point

    if ~isempty(regexp(str(1:10), '[a-zA-Z]', 'once'))
      % the line starts with letters, not numbers: probably it is a channel label
      % find the first number and remove the preceding part
      sel   = regexp(str(1:10), ' [-0-9]');   % find the first number, or actually the last space before the first number
      label = str(1:(sel));                   % this includes the space
      str   = str((sel+1):end);               % keep only the numbers
    end

    % convert the string into numbers and copy the desired samples over
    % into the data matrix
    tmp = str2num(str);
    dat(chan,:) = tmp(begsample:endsample);
  end
  fclose(fid);

else
  error('unsupported sub-fileformat');
end


text:
function test_tutorial_connectivityextended

% WALLTIME 00:45:00
% MEM 3gb
% DEPENDENCY ft_connectivityanalysis ft_connectivitysimulation ft_freqanalysis ft_connectivityplot ft_mvaranalysis

% This is the first section of the connectivity tutorial, which
% starts with an MVAR model and then uses parametric and nonparametric
% spectral decomposition for coherence and granger

% See also test_tutorial_connectivity2 and test_tutorial_connectivity3

%% simulate data
cfg             = [];
cfg.ntrials     = 500;
cfg.triallength = 1;
cfg.fsample     = 200;
cfg.nsignal     = 3;
cfg.method      = 'ar';
cfg.params(:,:,1) = [ 0.8    0    0 ;
  0  0.9  0.5 ;
  0.4    0  0.5];

cfg.params(:,:,2) = [-0.5    0    0 ;
  0 -0.8    0 ;
  0    0 -0.2];

cfg.noisecov      = [ 0.3    0    0 ;
  0    1    0 ;
  0    0  0.2];
data              = ft_connectivitysimulation(cfg);

figure
plot(data.time{1}, data.trial{1})
legend(data.label)
xlabel('time (s)')

cfg = [];
cfg.viewmode = 'vertical';  % you can also specify 'butterfly'
ft_databrowser(cfg, data);

%% mvaranalysis
cfg         = [];
cfg.order   = 5;
cfg.toolbox = 'bsmart';
mdata       = ft_mvaranalysis(cfg, data);

%% freqanalysis 1
cfg        = [];
cfg.method = 'mvar';
mfreq      = ft_freqanalysis(cfg, mdata);

%% freqanalysis 2
cfg           = [];
cfg.method    = 'mtmfft';
cfg.taper     = 'dpss';
cfg.output    = 'fourier';
cfg.tapsmofrq = 2;
freq          = ft_freqanalysis(cfg, data);

%% connectivityanalysis
cfg           = [];
cfg.method    = 'coh';
coh           = ft_connectivityanalysis(cfg, freq);
cohm          = ft_connectivityanalysis(cfg, mfreq);

%% visualisation
cfg           = [];
cfg.parameter = 'cohspctrm';
ft_connectivityplot(cfg, coh, cohm);

%% do the same for granger
cfg           = [];
cfg.method    = 'granger';
granger       = ft_connectivityanalysis(cfg, mfreq);

cfg           = [];
cfg.parameter = 'grangerspctrm';
ft_connectivityplot(cfg, granger);

figure
for row=1:3
  for col=1:3
    subplot(3,3,(row-1)*3+col);
    plot(granger.freq, squeeze(granger.grangerspctrm(row,col,:)))
    ylim([0 1])
  end
end

%% do the virtual channel stuff
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer_extended/source_coh_lft.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer_extended/source_diff.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer_extended/data_cmb.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer_extended/sourcemodel.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer_extended/hdm.mat'));

[maxval, maxcohindx] = max(source_coh_lft.avg.coh);
source_coh_lft.pos(maxcohindx, :)
[maxval, maxpowindx] = max(source_diff.avg.pow);
source_diff.pos(maxpowindx, :)

cfg                   = [];
cfg.covariance        = 'yes';
cfg.channel           = 'MEG';
cfg.vartrllength      = 2;
cfg.covariancewindow  = 'all';
tlock                 = ft_timelockanalysis(cfg, data_cmb);

% this is old-style stuff. as of end 2020 there's a ft_virtualchannel
% function that does the virtualchannel creation
cfg              = [];
cfg.method       = 'lcmv';
cfg.headmodel    = hdm;
cfg.sourcemodel.pos     = sourcemodel.pos([maxcohindx maxpowindx], :);
cfg.sourcemodel.inside  = true(2,1);
cfg.sourcemodel.unit    = sourcemodel.unit;
cfg.lcmv.keepfilter = 'yes';
source_idx       = ft_sourceanalysis(cfg, tlock);

%% old style
beamformer_lft_coh = source_idx.avg.filter{1};
beamformer_gam_pow = source_idx.avg.filter{2};

chansel = ft_channelselection('MEG', data_cmb.label); % find MEG sensor names
chansel = match_str(data_cmb.label, chansel);         % find MEG sensor indices

coh_lft_data = [];
coh_lft_data.label = {'coh_lft_x', 'coh_lft_y', 'coh_lft_z'};
coh_lft_data.time = data_cmb.time;
gam_pow_data = [];
gam_pow_data.label = {'gam_pow_x', 'gam_pow_y', 'gam_pow_z'};
gam_pow_data.time = data_cmb.time;
for i=1:length(data_cmb.trial)
  coh_lft_data.trial{i} = beamformer_lft_coh * data_cmb.trial{i}(chansel,:);
  gam_pow_data.trial{i} = beamformer_gam_pow * data_cmb.trial{i}(chansel,:);
end

cfg = [];
cfg.viewmode = 'vertical';  % you can also specify 'butterfly'
%ft_databrowser(cfg, gam_pow_data);

visualTimeseries = cat(2, gam_pow_data.trial{:});
motorTimeseries = cat(2, coh_lft_data.trial{:});
[u1, s1, v1] = svd(visualTimeseries, 'econ');
[u2, s2, v2] = svd(motorTimeseries, 'econ');

virtualchanneldata = [];
virtualchanneldata.label = {'visual', 'motor'};
virtualchanneldata.time = data_cmb.time;

for k = 1:length(data_cmb.trial)
  virtualchanneldata.trial{k}(1,:) = u1(:,1)' * beamformer_gam_pow * data_cmb.trial{k}(chansel,:);
  virtualchanneldata.trial{k}(2,:) = u2(:,1)' * beamformer_lft_coh * data_cmb.trial{k}(chansel,:);
end
virtualchanneldata_old = virtualchanneldata;

%% new-style
cfg = [];
cfg.pos            = source_idx.pos;
cfg.method         = 'svd';
cfg.numcomponent   = 1;
virtualchanneldata = ft_virtualchannel(cfg, data_cmb, source_idx);
virtualchanneldata.label = {'motor';'visual'}; % note the order is reversed w.r.t. old-style

% do a sanity check on whether the old and new style vcs match more or less
c = corr([virtualchanneldata_old.trial{1}' virtualchanneldata.trial{1}']);
assert(c(2,3)>0.99);
assert(c(1,4)>0.99);

% select the two EMG channels
cfg = [];
cfg.channel = 'EMG';
emgdata = ft_selectdata(cfg, data_cmb);

% combine the virtual channel with the two EMG channels
cfg = [];
combineddata = ft_appenddata(cfg, virtualchanneldata, emgdata);

%% compute the spectral decomposition
cfg            = [];
cfg.output     = 'fourier';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.tapsmofrq  = 5;
cfg.keeptrials = 'yes';
cfg.channel    = {'visual' 'motor' 'EMGlft' 'EMGrgt'};
freq    = ft_freqanalysis(cfg, combineddata);

cfg = [];
cfg.method = 'coh';
coherence = ft_connectivityanalysis(cfg, freq);

cfg = [];
cfg.zlim = [0 0.25];
figure
ft_connectivityplot(cfg, coherence);


text:
function [hPlot] = plot_CDF(vDistribution, hAxes, sColor)
    % Plots the cumulative density function of a given distribution.
    %
    % [hPlot] = plot_CDF(vDistribution, hAxes, sColor)
    % ---------------------------------------------------------
    % NaN values are removed from the distribution.
    %
    % plot_CDF(vDistribution) opens a figure and plots the cumulative density function
    %   with default parameters.
    %
    % Input parameters:
    %   vDistribution   Distribution of values to be plotted
    %   hAxes           Handle of axes to plot the cumulative density function
    %   sColor          Color of plot (refer to Matlab 'plot')
    %
    % Output parameters:
    %   hPlot           Handle of the plot
    %
    % Danijel Schorlemmer
    % July 9, 2003

    report_this_filefun();

    % Define missing input parameters
    if ~exist('hAxes', 'var')
        figure;
        hAxes = newplot;
    end
    if ~exist('sColor', 'var')
        sColor = 'k';
    end

    % Remove NaN-values
    vSel = ~isnan(vDistribution);
    vPlotDist = vDistribution(vSel,:);

    % Activate given axes
    axes(hAxes);

    % Plot the cumulative density function
    nLen = length(vPlotDist);
    vIndices = [1:nLen]/nLen;
    vDist = sort(vPlotDist);
    hPlot = plot(vDist, vIndices, sColor);
end

author: CelsoReyes


text:
function varargout = pnopt_backtrack( x, d, t, f_x, h_x, dg_x, smoothF, ...
  nonsmoothF, desc_param, xtol, maxIter )
% pnopt_backtrack : Backtracking line search for step that satisfies a sufficient
%   descent condition.
%
%   $Revision: 0.8.0 $  $Date: 2012/12/01 $
%
% --------------------Initialize--------------------

  % Set line search parameters
  beta = 0.5;

  % Set termination flags
  FLAG_SUFFDESC = 1;
  FLAG_TOLX     = 2;
  FLAG_MAXFUNEV = 3;

  iter = 0;

  desc = dg_x + nonsmoothF( x + d ) - h_x;

  % --------------------Main Loop--------------------
  while 1
    iter = iter + 1;

    % Evaluate trial point and function value.
    y = x + d * t;
    if nargout > 6
      [ g_y, Dg_y, D2g_y ] = smoothF( y );
    else
      [ g_y, Dg_y ] = smoothF( y );
    end
    h_y = nonsmoothF( y );
    f_y = g_y + h_y;

    % Check termination criteria
    if f_y < f_x + desc_param * t * desc         % Sufficient descent condition satisfied
      flag = FLAG_SUFFDESC;
      break
    elseif t <= xtol            % Step length too small
      flag = FLAG_TOLX;
      break
    elseif iter >= maxIter      % Too many line search iterations
      flag = FLAG_MAXFUNEV;
      break
    end

    % Backtrack if objective value not well-defined
    if isnan( f_y ) || isinf( f_y ) || abs( f_y - f_x - t * dg_x ) <= 1e-9
      t = beta * t;
    % Safeguard quadratic interpolation
    else
      t_interp = - ( dg_x * t^2) / ( 2 * ( f_y - f_x - t * dg_x ) );
      if 0.01 <= t_interp && t_interp <= 0.99*t
        t = t_interp;
      else
        t = beta * t;
      end
    end

  end

  if nargout > 6
    varargout = { y, f_y, Dg_y, D2g_y, t, flag ,iter };
  else
    varargout = { y, f_y, Dg_y, t, flag ,iter };
  end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
clear all; close all; clc;

% pull=[3 5 7 9 11];
% start=1:20;
% q=abs(start-pull(1));
% [x,n2]=min(q);
%
% for j=1:4
%    q=abs(start-pull(j));
%    [x,n2]=min(q)
%    start=start([1:n2-1 n2+1:(20-j)]);
%    start'
% end
%
% break

%%
L=10; x3=-L:0.1:L; n=length(x3)-1; % define domain
x2=x3(1:n); k=(2*pi/(2*L))*[0:n/2-1 -n/2:-1]; % k-vector
ye=exp(-(x2.^2)); ye2=exp((x2.^2)/2); % define Gaussians
for j=0:9      % loop through 10 modes
  yd=real(ifft(((i*k).^j).*fft(ye))); % 2nd derivative
  mode=((-1)^(j))*(((2^j)*factorial(j)*sqrt(pi))^(-0.5))*ye2.*yd;
  y(:,j+1)=(mode).';  % store modes as columns
end

x=x2(n/2+1-40:n/2+1+40);  % keep only -4<x<4
yharm=y(n/2+1-40:n/2+1+40,:);

%%
n=length(x);
f=(exp(-(x-0.5).^2)+3*exp(-2*(x+1.5).^2))';

for j=1:10  % full reconstruction
  a(j,1)=trapz(x,f.*yharm(:,j));
end
f2=yharm*a;
Efull(1)=log(norm(f2-f)+1);  % reconstruction error

for j=1:10  % matrix M reconstruction
   for jj=1:j
       Area=trapz(x,yharm(:,j).*yharm(:,jj));
       M(j,jj)=Area;
       M(jj,j)=Area;
   end
end
Cfull=cond(M)   % get condition number

%%  Willcox alg.

% 23    52    37    79

% sensor 1
for jloop=1:81

s=zeros(n,1);
s(jloop)=1;

for j=1:10
   for jj=1:j
       Area=trapz(x,s.*(yharm(:,j).*yharm(:,jj)));
       M2(j,jj)=Area; M2(jj,j)=Area;
   end
end
%figure(3), subplot(2,2,1), pcolor(10:-1:1,1:10,(M2'));, colormap(hot)

for j=1:10  % reconstruction using gappy
  ftild(j,1)=trapz(x,s.*(f.*yharm(:,j)));
end
atild=M2\ftild;   % compute error
f1=yharm*atild;
Er(1)=norm(f1-f);

con1(jloop)=cond(M2);

end

[s1,n1]=min(con1)

% sensor 2
jlook=[1:n1-1 n1+1:81];
for jloop=jlook

s=zeros(n,1);
s(n1)=1;
s(jloop)=1;

for j=1:10
   for jj=1:j
       Area=trapz(x,s.*(yharm(:,j).*yharm(:,jj)));
       M2(j,jj)=Area; M2(jj,j)=Area;
   end
end
%figure(3), subplot(2,2,2), pcolor(10:-1:1,1:10,(M2'));, colormap(hot)

for j=1:10  % reconstruction using gappy
  ftild(j,1)=trapz(x,s.*(f.*yharm(:,j)));
end
atild=M2\ftild;   % compute error
f2=yharm*atild;
Er(2)=norm(f2-f);

con2(jloop)=cond(M2);

end
[s2,n2]=min(con2(jlook))

% sensor 3
if n2>n1
   jlook=[1:n1-1 n1+1:n2-1 n2+1:81];
else
   jlook=[1:n2-1 n2+1:n1-1 n1+1:81];
end

for jloop=jlook

s=zeros(n,1);
s(n1)=1; s(n2)=1;
s(jloop)=1;

for j=1:10
   for jj=1:j
       Area=trapz(x,s.*(yharm(:,j).*yharm(:,jj)));
       M2(j,jj)=Area; M2(jj,j)=Area;
   end
end
%figure(3), subplot(2,2,3), pcolor(10:-1:1,1:10,(M2'));, colormap(hot)

for j=1:10  % reconstruction using gappy
  ftild(j,1)=trapz(x,s.*(f.*yharm(:,j)));
end
atild=M2\ftild;   % compute error
f3=yharm*atild;
Er(3)=norm(f3-f);

con3(jloop)=cond(M2);

end
[s3,n3]=min(con3(jlook))

% sensor 4
jlook=[1:n1-1 n1+1:n3-1 n3+1:n2-1 n2+1:81];
for jloop=jlook

s=zeros(n,1);
s(n1)=1; s(n2)=1;
s(jloop)=1;

for j=1:10
   for jj=1:j
       Area=trapz(x,s.*(yharm(:,j).*yharm(:,jj)));
       M2(j,jj)=Area; M2(jj,j)=Area;
   end
end
%figure(3), subplot(2,2,4),pcolor(10:-1:1,1:10,(M2'));, colormap(hot)

for j=1:10  % reconstruction using gappy
  ftild(j,1)=trapz(x,s.*(f.*yharm(:,j)));
end
atild=M2\ftild;   % compute error
f4=yharm*atild;
Er(4)=norm(f4-f);

con4(jloop)=cond(M2);

end
[s4,n4]=min(con4(jlook))

figure(1)
subplot(4,1,1), bar(log(con1),'Facecolor',[0.7 0.7 0.7]), axis([1 81 0 100])
subplot(4,1,2), bar(log(con2),'Facecolor',[0.7 0.7 0.7]), axis([1 81 0 100])
subplot(4,1,3), bar(log(con3),'Facecolor',[0.7 0.7 0.7]), axis([1 81 0 100])
subplot(4,1,4), bar(log(con4),'Facecolor',[0.7 0.7 0.7]), axis([1 81 0 100])

temp=zeros(n,1), temp(n1)=s1; subplot(4,1,1),hold on, bar(log(temp),'r')
temp=zeros(n,1), temp(n2)=s2; subplot(4,1,2),hold on, bar(log(temp),'r')
temp=zeros(n,1), temp(n3)=s3; subplot(4,1,3),hold on, bar(log(temp),'r')
temp=zeros(n,1), temp(n4)=s4; subplot(4,1,4),hold on, bar(log(temp),'r')

subplot(4,1,1), set(gca,'Xtick',[0 40 80],'Ytick',[0 50 100],'Xticklabel',{'','',''},'Yticklabel',{'','',''})
subplot(4,1,2), set(gca,'Xtick',[0 40 80],'Ytick',[0 50 100],'Xticklabel',{'','',''},'Yticklabel',{'','',''})
subplot(4,1,3), set(gca,'Xtick',[0 40 80],'Ytick',[0 50 100],'Xticklabel',{'','',''},'Yticklabel',{'','',''})
subplot(4,1,4), set(gca,'Xtick',[0 40 80],'Ytick',[0 50 100],'Xticklabel',{'','',''},'Yticklabel',{'','',''})

break

figure(2)
subplot(4,1,1), plot(x,f1,'k',x,f,'r','Linewidth',[2])
subplot(4,1,2), plot(x,f2,'k',x,f,'r','Linewidth',[2])
subplot(4,1,3), plot(x,f3,'k',x,f,'r','Linewidth',[2])
subplot(4,1,4), plot(x,f4,'k',x,f,'r','Linewidth',[2])

%%  Willcox:  condition number

% pull=[3 5 7 9 11];
% start=1:20;
% q=abs(start-pull(1));
% [x,n2]=min(q);
%
% for j=1:4
%    q=abs(start-pull(j));
%    [x,n2]=min(q)
%    start=start([1:n2-1 n2+1:(20-j)]);
%    start'
% end
%
% break
clear f1
clear f2
clear ns
ns=[];            % vector to keep track of measurement locations
jlook=1:81;   % loop over 81 positions

count=1;
for jsense=1:20   %  loop for 20 sensors

for jloop=1:(82-jsense)   % loop from 81 to 61 locations
s=zeros(n,1); s(ns)=1;
s(jlook(jloop))=1;

for j=1:10
   for jj=1:j    % matrix M
       Area=trapz(x,s.*(yharm(:,j).*yharm(:,jj)));
       M2(j,jj)=Area; M2(jj,j)=Area;
   end
end
con(jloop)=cond(M2);  % compute condition number

end

[s1,n1]=min(con);   % find minimum condition number location
kond(jsense)=s1; clear con
ns=[ns n1];  % add sensor location

jlook=jlook([1:n1-1 n1+1:(81-count+1)]);
count=count+1;

% reconstruct
s=zeros(n,1);  s(ns)=1;

for j=1:10
   for jj=1:j
       Area=trapz(x,s.*(yharm(:,j).*yharm(:,jj)));
       M2(j,jj)=Area; M2(jj,j)=Area;
   end
end

for j=1:10  % reconstruction using gappy
  ftild(j,1)=trapz(x,s.*(f.*yharm(:,j)));
end
atild=M2\ftild;   % compute error
f1(:,jsense)=yharm*atild;  % iterative reconstruction
Errr(jsense)=norm(f1(:,jsense)-f);  % iterative error
scum(:,jsense)=s;     % track iterative sensing matrix

end

figure(4)
subplot(2,1,2), bar(log(Errr+1))
subplot(2,1,1), bar(log(kond))

subplot(2,1,2), set(gca,'Xlim',[0 21],'Xtick',[1 5 10 15 20],'Xticklabel',{'','','','',''},'Ylim',[0 12],'Ytick',[0 4 8 12],'Yticklabel',{'','','',''})
subplot(2,1,1), set(gca,'Xlim',[0 21],'Xtick',[1 5 10 15 20],'Xticklabel',{'','','','',''},'Ylim',[0 42],'Ytick',[0 20 40],'Yticklabel',{'','',''})

%figure(6), bar(s)

figure(7),
titer=[1:20 25]; f1=[f1 f]; titer2=[9:20 25], f2=[f1(:,9:20) f];
subplot(2,2,1), waterfall(x,titer,f1.'), colormap([0 0 0]), view(-150,50)
subplot(2,2,2), waterfall(x,titer2,f2.'), colormap([0 0 0]), view(-150,50)

subplot(2,2,1), set(gca,'Zlim',[-2 8],'Xlim',[-4 4],'Xtick',[-4 0 4],'Xticklabel',{'','',''}, ...
    'Ylim',[0 25],'Ytick',[1 10 20],'Yticklabel',{'','',''},'Ztick',[-2 0 4 8],'Zticklabel',{'','','',''})
subplot(2,2,2), set(gca,'Zlim',[-.2 3],'Xlim',[-4 4],'Xtick',[-4 0 4],'Xticklabel',{'','',''}, ...
    'Ylim',[0 25],'Ytick',[1 10 20],'Yticklabel',{'','',''},'Ztick',[0 1 2 3],'Zticklabel',{'','','',''})

figure(8)
%bar3(scum)
axes('position',[.05 .05 .9 .35])
scum2=[scum(:,[1:11 13:19])];
pcolor(-scum2.'), colormap(hot), axis off

%bar3(scum2)

%caxis([0 1]), axis off

author: dynamicslab


text:
% Test file for singfun/diff.m

function pass = test_diff(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
d = 2;
x = 2*(1-10^(-d)) * rand(100, 1) - (1-10^(-d));

% The order of the exponents:
a = 0.56;
b = -0.56;
c = 1.28;
d = -1.28;

% Check for empty cases
f = singfun();
pass(1) = isempty(diff(f));

% Return type should not be a SINGFUN when a smooth function is differentiated.
% [TODO]: The return type should be SMOOTHFUN
f = singfun(@(x) sin(x));
g = diff(f);
pass(2) = ~isa(g, 'singfun');
%%
% Spot-check derivatives for a couple of functions.

% fractional root at the left endpoint
data.exponents = [a 0];
data.singType = {'root', 'none'};
f = singfun(@(x) (1+x).^a.*exp(x), data, pref);
df = diff(f);
vals_df = feval(df, x);
df_exact = @(x) (1+x).^(a-1).*(a+1+x).*exp(x);
vals_exact = feval(df_exact, x);
err = vals_df - vals_exact;
pass(3) = (norm(err, inf) < 1e2*eps*norm(vals_exact, inf));

% fractional pole at the left endpoint
data.exponents = [d 0];
data.singType = {'sing', 'none'};
f = singfun(@(x) (1+x).^d.*sin(x), data, pref);
df = diff(f);
vals_df = feval(df, x);
df_exact = @(x) (1+x).^(d-1).*(d*sin(x)+(1+x).*cos(x));
vals_exact = feval(df_exact, x);
err = vals_df - vals_exact;
pass(4) = (norm(err, inf) < 1e2*eps*norm(vals_exact, inf));

% fractional root at the right endpoint
data.exponents = [0 c];
data.singType = {'none', 'root'};
f = singfun(@(x) (1-x).^c.*cos(x), data, pref);
df = diff(f);
vals_df = feval(df, x);
df_exact = @(x) -(1-x).^(c-1).*(c*cos(x)+(1-x).*sin(x));
vals_exact = feval(df_exact, x);
err = vals_df - vals_exact;
pass(5) = (norm(err, inf) < 1e2*eps*norm(vals_exact, inf));

% fractional pole at the right endpoint
data.exponents = [0 b];
data.singType = {'none', 'sing'};
f = singfun(@(x) (1-x).^b.*(x.^5), data, pref);
df = diff(f);
vals_df = feval(df, x);
df_exact = @(x) (1-x).^(b-1).*(5-5*x-b*x).*(x.^4);
vals_exact = feval(df_exact, x);
err = vals_df - vals_exact;
pass(6) = (norm(err, inf) < 1e2*eps*norm(vals_exact, inf));

% a combination of fractional pole and fractional root
data.exponents = [b c];
data.singType = {'sing', 'root'};
f = singfun(@(x) (1+x).^b.*sin(x).*(1-x).^c, data, pref);
df = diff(f);
vals_df = feval(df, x);
df_exact = @(x) cos(x).*(1 - x).^c.*(x + 1).^b +...
    b*sin(x).*(1 - x).^c.*(x + 1).^(b - 1)...
    - c*sin(x).*(1 - x).^(c - 1).*(x + 1).^b;
vals_exact = feval(df_exact, x);
err = vals_df - vals_exact;
pass(7) = (norm(err, inf) < 1e2*eps*norm(vals_exact, inf));

%%
% Verify that calling diff() gives the reasonably accurate answer as direct
% construction.

data.exponents = [b b];
data.singType = {'sing', 'sing'};
f = singfun(@(x) (1+x).^b.*sin(2*x).*(1-x).^b, data, pref);
df = diff(f);
vals_df = feval(df, x);
data.exponents = [(b - 1) (b - 1)];
data.singType = {'sing', 'sing'};
df_exact = singfun(@(x) -2*(1 - x).^(b-1).*(x + 1).^(b-1) ...
    .*(x.^2.*cos(2*x) - cos(2*x) + b*x.*sin(2*x)), data, pref);
vals_exact = feval(df_exact, x);
err = vals_df - vals_exact;
pass(8) = (norm(err, inf) < 20*eps*norm(vals_exact, inf));

%%
% Check higher-order derivatives.

data.exponents = [a b];
data.singType = {'root', 'sing'};
f = singfun(@(x) (1+x).^a.*sin(x).*(1-x).^b, data, pref);
df2 = diff(f, 2);
vals_df2 = feval(df2, x);
df2_exact = @(x) 2*a*cos(x).*(1-x).^b.*(x+1).^(a-1)-...
    sin(x).*(1-x).^b.*(x+1).^a-2*b*cos(x).*(1-x).^(b-1).*(x+1).^a+...
    a*sin(x).*(a-1).*(1-x).^b.*(x+1).^(a-2)-...
    2*a*b*sin(x).*(1-x).^(b-1).*(x+1).^(a-1)+...
    b*sin(x).*(b-1).*(1-x).^(b-2).*(x+1).^a;
vals_exact = feval(df2_exact, x);
err = vals_df2 - vals_exact;
pass(9) = (norm(err, inf) < 1e2*eps*norm(vals_exact, inf));

end


text:
% This script file loads input files defined in vector 'n'  and
% selects some eartquakes based on various criteria
% (e.g. latidude, longitude, magnitude, depth error)

% The selected EQ are stored in matrix 'a'.
%
%  Stefan Wiemer, June 1994

report_this_filefun(mfilename('fullpath'));

clear
a = []
n = [ 'calsplitaa'
    'calsplitab'
    'calsplitac'
    'calsplitad'
    'calsplitae'
    'calsplitaf'
    'calsplitag'
    'calsplitah'
    'calsplitai'
    'calsplitaj'
    'calsplitak'
    'calsplital'
    'calsplitam'
    'calsplitan'
    'calsplitao'
    'calsplitap'
    'calsplitaq'
    'calsplitar'
    ]

for i = 1:length(n)

    lofi = ['load ' n(i,1:10) ]
    eval(lofi)
    comm = [' s = ' n(i,1:10) ';']
    eval(comm)
    a2 = [ -s(:,1) s(:,2) s(:,3) s(:,4) s(:,5) s(:,6) s(:,7)/100 s(:,8)/100];
    l = a2(:,8) < 2.0 & a2(:,1) < -121.2 & a2(:,1) > -122.4 & a2(:,2) > 36.65  & a2(:,2) < 37.8;
    a2 = a2(l,:);

    a = [a ; a2];
    size(a)
    comm = [' clear a2 s ' n(i,1:10)]
    eval(comm)
end

a2 = a;
load /FM1/ramon/matlab/zmapv1.1/eq_data/landers_cata.mat
a = a2;
clear a2;
save /FM1/ramon/matlab/zmapv1.1/eq_data/centcal_cata.mat


text:
%%% Nowcasting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script produces a nowcast of US real GDP growth for 2016:Q4
% using the estimated parameters from a dynamic factor model.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Clear workspace and set paths.
close all; clear; clc;
addpath('functions');

%% User inputs.
series = 'GDPC1' ; % Nowcasting real GDP (GDPC1) <fred.stlouisfed.org/series/GDPC1>
period = '2016q4'; % Forecasting target quarter

%% Load model specification and first vintage of data.
% Load model specification structure `Spec`
Spec = load_spec('Spec_US_example.xls');

%% Load DFM estimation results structure `Res`.
Res = load('ResDFM'); % example_DFM.m used the first vintage of data for estimation

%% Update nowcast and decompose nowcast changes into news.

%%% Nowcast update from week of December 7 to week of December 16, 2016 %%%
vintage_old = '2016-12-16'; datafile_old = fullfile('data','US',[vintage_old '.xls']);
vintage_new = '2016-12-23'; datafile_new = fullfile('data','US',[vintage_new '.xls']);
% Load datasets for each vintage
[X_old,~   ] = load_data(datafile_old,Spec);
[X_new,Time] = load_data(datafile_new,Spec);

% check if spec used in estimation is consistent with the current spec
if isequal(Res.Spec,Spec)
    Res = Res.Res;
else
    threshold = 1e-4; % Set to 1e-5 for robust estimates
    Res = dfm(X_new,Spec,threshold);
    save('ResDFM','Res','Spec');
end

update_nowcast(X_old,X_new,Time,Spec,Res,series,period,vintage_old,vintage_new);

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

author: FRBNY-TimeSeriesAnalysis

repo: Nowcasting


text:
function dosefilenames = find_pinnacle_dose_filenames(dosedim,planPath)
%
%	dosefilenames = find_pinnacle_dose_filenames(dosedim)
%
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%APA
% planPath = 'C:\Projects\PinnacleReader\PinnacleWholePatientPlanExample\Plan_2\plan.Trial.binary.*';
files = dir(planPath);
%files = dir('plan.Trial.binary.*');
N = length(files);

dosefilenames = cell(0);

for k = 1:N
	if files(k).bytes == prod(dosedim)*4
		if isempty(dosefilenames)
			dosefilenames{1} = files(k).name;
		else
			dosefilenames{end+1,1} = files(k).name;
		end
	end
end

% The number of binary dose files should be equal to 2 * number of beams

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function path_length = iTreePathLength(x, treenode, current_path_length)
%UNTITLED7 Summary of this function goes here
%   Detailed explanation goes here

% if strcmp(class(treenode), 'iTreeLeaf')
%    path_length = current_path_length + adjustment(treenode.Size);
% else
%    if x(treenode.SplitAttribute) < treenode.SplitValue
%        path_length = iTreePathLength(x, treenode.Left, current_path_length+1);
%    else
%        path_length = iTreePathLength(x, treenode.Right, current_path_length+1);
%    end
% end

path_length = 0;
while ~strcmp(class(treenode), 'iTreeLeaf')
    if x(treenode.SplitAttribute) < treenode.SplitValue
        treenode = treenode.Left;
    else
        treenode = treenode.Right;
    end
    path_length = path_length + 1;
end
path_length = path_length + adjustment(treenode.Size);
end

function value=adjustment(n)
%The average path length of an unsuccessful BST search
    if n<=1
        value = 0;
    else
        value = 2 * (log(n-1)+0.5772156649) - 2*(n-1) / n;
    end
end

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox