text:
%Once the image is processed using the gridextractor method, this function
%tries to fit the grid on the extraced points in order to establish the
%correspondences
%
%Input:     stats   The processed blobs of this image
%           id0     The index of the main bar
%           id1     The index of the second main bar
%           searchPoints    The potential grid points
%           doshow  1/0 to display / not display the results
%
%Output:    success The flag indicating whether it worked or not
%           x_      the 2D points of the grid points
%           X_      the corresponding 3D points of the grid points
%           searchPoints    The points
%
%This file uses some functions from Peter Kovesi's Matlab functions
%http://www.csse.uwa.edu.au/~pk/Research/MatlabFns/index.html
%
%Christian Wengert
%Computer Vision Laboratory
%ETH Zurich
%Sternwartstrasse 7
%CH-8092 Zurich
%www.vision.ee.ethz.ch/cwengert
%wengert@vision.ee.ethz.ch
function [success, x_, X_, searchPoints, error] = fitgrid(stats, searchPoints, id0, id1, dX, dY, show)

    %Grid small
    xmin = -12*dX;
    xmax = 12*dX;
    ymin = -12*dY;
    ymax = 12*dY;
    %Define a delta distance
    delta = 7; %[pixel]

    %Init
    x0 = [stats(id0).Centroid(1);stats(id0).Centroid(2)];
    x1 = [stats(id1).Centroid(1);stats(id1).Centroid(2)];
    dir = (x1-x0)/(norm(x1-x0));
    %Get a 90?angle
    dir90_ = [-dir(2);dir(1)];
    %get the orientation of main bar
    orientation = stats(id0).Orientation;
    y = tand(orientation);
    dir90 = [-1;y];
    dir90 = dir90/norm(dir90);

    %Compare direction of orientation with dir90_ and dir90
    sd90 = sign(dir90); sd90_  =sign(dir90_);
%     [dir90,dir90_]
    if(abs(1-abs(dir90_(1)))<0.01 & sd90(1)~=sd90_(1))
        dir90 = -dir90;
    elseif(abs(1-abs(dir90_(2)))<0.01 & sd90(2)~=sd90_(2))
        dir90 = -dir90;
    elseif(~(sd90(1)==sd90_(1) |sd90(2)==sd90_(2)))
        dir90 = -dir90;
    end

    error = 0;
    %Draw the initial lines
    if(show)
        line([x0(1), x0(1)+dir(1)*100], [x0(2), x0(2)+dir(2)*100],'Color','r');
        line([x0(1), x0(1)-dir(1)*100], [x0(2), x0(2)-dir(2)*100],'Color','r');
        line([x0(1), x0(1)+dir90(1)*100], [x0(2), x0(2)+dir90(2)*100]);
        line([x0(1), x0(1)-dir90(1)*100], [x0(2), x0(2)-dir90(2)*100]);
    end
    deltaDist2Line = 10;
    try
        %Create search Points
        a1 = x0+dir90*stats(id0).MajorAxisLength/2;
        a2 = x0-dir90*stats(id0).MajorAxisLength/2;
        b1 = x0-dir*stats(id1).MajorAxisLength/2;
        b2 = x1+dir*stats(id1).MajorAxisLength/2;
        if(show)
            plot(a1(1),a1(2),'mo','MarkerSize',40,'LineWidth',2)
            plot(a2(1),a2(2),'bo','MarkerSize',40,'LineWidth',2)
            plot(b1(1),b1(2),'co','MarkerSize',40,'LineWidth',2)
            plot(b2(1),b2(2),'ro','MarkerSize',40,'LineWidth',2)
        end
        da1 = [];da2 = []; db1 = []; db2 = [];
        for i=1:length(searchPoints)
            da1 = [da1;norm(a1 - searchPoints(2:3,i))];
            da2 = [da2;norm(a2 - searchPoints(2:3,i))];
            db1 = [db1;norm(b1 - searchPoints(2:3,i))];
            db2 = [db2;norm(b2 - searchPoints(2:3,i))];
        end
        [da1,ia1] = sort(da1);
        [da2,ia2] = sort(da2);
        [db1,ib1] = sort(db1);
        [db2,ib2] = sort(db2);
        ix0 = ia1(1);
        ix1 = ia2(1);
        iy(1) = ib1(1);
        iy(2) = ib2(1);

        %Show the points
        if(show)
            plot(searchPoints(2,ix0(1)),searchPoints(3,ix0(1)),'mx','MarkerSize',25,'LineWidth',2)
            plot(searchPoints(2,ix1(1)),searchPoints(3,ix1(1)),'bx','MarkerSize',25,'LineWidth',2)
            plot(searchPoints(2,iy(1)),searchPoints(3,iy(1)),'cx','MarkerSize',25,'LineWidth',2)
            plot(searchPoints(2,iy(2)),searchPoints(3,iy(2)),'rx','MarkerSize',25,'LineWidth',2)
        end
        %Refine the line
        dir = searchPoints(2:3,ix0(1)) - searchPoints(2:3,ix1(1));
        dir90 = searchPoints(2:3,iy(1)) - searchPoints(2:3,iy(2));

        %Make the initial grid
        [x,y] = meshgrid(xmin:dX:xmax, ymin:dY:ymax);
        x = reshape(x, 1, numel(x));
        y = reshape(y, 1, numel(y));
        xy = [x;y];
        usedPtsIndex = ones(1,length(x));
        %Get the homography
        %the points we have are: [pixels]
        xy__ = [searchPoints(2:3,ix0(1)), searchPoints(2:3,ix1(1)),searchPoints(2:3,iy(1)),searchPoints(2:3,iy(2))];
        xy_  = [[0;-2*dX],[0;2*dX], [-dY;0],[3*dY;0]];
        %Compute H
        [H, err] = invpersp(xy__, xy_);
        xyp = homoTrans(H,[xy;ones(1,length(xy))]);

        %Distance check
        for i=1:length(searchPoints)
            for j=1:length(xyp)
                d = norm(searchPoints(2:3,i)-xyp(1:2,j));
                error = error+d;
                %Check whether its good
                if(d<=delta)

                    if(show)
                        text(searchPoints(2,i)+5,searchPoints(3,i)+5,[num2str(xy(1,j)) ',' num2str(xy(2,j))],'Color','b')
                        plot(searchPoints(2,i), searchPoints(3,i),'r+');
                    end
                    searchPoints(4:5,i) = xy(1:2,j);
                    %Update the point vectors that they have been assigned
                    searchPoints(1,i) = 0;
                    usedPtsIndex(j) = 0;
                end
            end
        end
        success  = 1;
        finalindex = find(searchPoints(1,:)==0);
        x_ = searchPoints(2:3,finalindex);
        X_ = [searchPoints(4:5,finalindex)];
        X_ = [X_;zeros(1,length(X_))];

        error = error/length(x_)/length(xyp);
    catch
        disp('fitgrid::Could not find extension points')
        success = 0;
        x_ = []; X_ = [];
    end

author: christianwengert


text:
%  blockwise_optimization optimization by blocks of parameters rather than the whole vector
%
%  ::
%
%
%  Args:
%
%  Returns:
%     :
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih


text:
function [out] = interflow_5(p1,S)
%interflow_5

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Linear interflow
% Constraints:  -
% @(Inputs):    p1   - time coefficient [d-1]
%               S    - current storage [mm]

out = p1.*S;

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
% SELFINTERSECT Given a triangle mesh (V,F) compute a new mesh (VV,FF) which is the same as
% (V,F) except that any self-intersecting triangles in (V,F) have been
% subdivided (new vertices and face created) so that the self-intersection
% contour lies exactly on edges in (VV,FF). New vertices will appear in
% original faces or on original edges. New vertices on edges are "merged" only
% across original faces sharing that edge. This means that if the input
% triangle mesh is a closed manifold the output will be too.
%
% [VV,FF,IF] = selfintersect(V,F,'ParameterName',ParameterValue, ...)
%
% Inputs:
%   V  #V by 3 list of vertex positions
%   F  #F by 3 list of triangle indices into V
%   Optional:
%     'DetectOnly'  followed by bool. Whether to only detect intersecting pairs
%       (sets IF but not VV or FF) {false}
%     'FirstOnly'  followed by bool whether to only detect the first
%       intersection. {false}
%     'StitchAll'  followed by whether to stitch all vertices in the output, if
%       true then IM will be 1:size(VV,1) {false}
% Outputs:
%   VV  #VV by 3 list of vertex positions
%   FF  #FF by 3 list of triangle indices into V
%   IF  #intersecting face pairs by 2  list of intersecting face pairs,
%     indexing F
%   J   #FF list of indices into F of birth parents
%   IM  #VV list of indices into VV of unique vertices
%
% Example:
%   [SV,SF,~,~,IM] = selfintersect(V,F);
%   FF = IM(SF);
%   [U,IM] = remove_unreferenced(SV,FF);
%   G = IM(FF);
%
%   % Self-Intersect (V,F)+(U,G) and separate
%   VU = [V;U];
%   FG = [F;size(V,1)+G];
%   [SVU,SFG,~,J] = selfintersect(VU,FG);
%   SF = SFG(J<=size(F,1),:);
%   SG = SFG(J>size(F,1),:);
%   [SV,IM] = remove_unreferenced(SVU,SF);
%   SF = IM(SF);
%   [SU,IM] = remove_unreferenced(SVU,SG);
%   SG = IM(SG);
%
%

% See selfintersect.h, selfintersect.cpp for mex implementation

author: alecjacobson

repo: gptoolbox


text:
function [F_N, Y_N] = leastSquaresSHT(N, F, dirs, basisType, weights)
%LEASTSQUARESSHT Spherical harmonic transform of F using least-squares
%
%   N:  maximum order of transform
%   F:  the spherical function evaluated at K directions 'dirs'
%       if F is a KxM matrix then the transform is applied to each column
%       of F, and returns the coefficients at each column of F_N
%       respectively
%   dirs:   [azimuth1 inclination1; ...; azimuthK inclinationK] angles in
%           rads for each evaluation point, where inclination is the polar
%           angle from zenith: inclination = pi/2-elevation
%   basisType:  'complex' or 'real' spherical harmonics
%   weights:    vector of K weights for each direction, for weighted
%               least-squares
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Archontis Politis, 10/10/2013, updated 7/2/2015
%   archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % compute the harmonic coefficients
    Y_N = getSH(N, dirs, basisType);
    Npoints = size(dirs,1);

    % non-weighted case for uniform arrangements
    if nargin<5 || ~exist('weights','var') || isempty(weights)

        % perform transform in the least squares sense
        F_N = pinv(Y_N)*F;
    else

        % perform weighted least-squares transform
        F_N = (Y_N'*diag(weights)*Y_N) \ (Y_N'*diag(weights) * F);
    end

end

author: polarch

repo: Spherical-Harmonic-Transform


text:
%% Homework 12 (a)
% Heat equation with variable coeficients using Point Jacobi Method:
%
% $\frac{\partial}{\partial x}\left (k(x,y)\frac{\partial{T}}{\partial{x}} \right )+\frac{\partial}{\partial x}\left (k(x,y)\frac{\partial{T}}{\partial{y}} \right)=0$
%
% where:
%
% $k(x,y)=1+exp\left[-\left( \frac{x-0.25-0.5y}{0.25} \right)^2 \right]$

%close all
clear all

%% Grid
h = 0.02;
lx = [0,1]; dx = h; x = lx(1):dx:lx(2); n = length(x);
ly = [0,1]; dy = h; y = ly(1):dy:ly(2); m = length(y);
T = zeros(m,n); T_next = zeros(m,n);

% Computing k(x,y) on the grid
k = zeros(m,n);
for j=1:m
    for i=1:n
        %k(j,i)=1+exp(-((x(i)-0.25-0.5*y(j))/0.25)^2);
        k(j,i)=1;
    end
end

%% Boundary Conditions
T(1,:) = 0; % Bottom (last row of the matrix)
T(m,:) = 1-cos(pi*x); % Top (first row of the matrix)
% Left and Right BCs will be updated on every iteration

%% Method of Solution
%
A = zeros(m,n); B = zeros(m,n); C = zeros(m,n); D = zeros(m,n);
r = 0.00001;
iter = 0;
r_iter = 1;
while r_iter >= r
%for s=1:5000
    for i=2:m-1
        for j=2:n-1
            A(i,j)=k(i,j)+k(i+1,j);
            B(i,j)=k(i,j)+k(i,j+1);
            C(i,j)=k(i,j)+k(i-1,j);
            D(i,j)=k(i,j)+k(i,j-1);
            if j == 2 % Matrix left Column, Neumann BC
                T_next(i,j)=...
                    1/(k(i+1,j)+k(i-1,j)+3*k(i,j)+k(i,j+1))*(...
                    A(i,j)*T(i+1,j)+...
                    B(i,j)*T(i,j+1)+...
                    C(i,j)*T(i-1,j));
            elseif j == n-1 % Matrix right Column, Neumann BC
                T_next(i,j)=...
                    1/(k(i+1,j)+k(i-1,j)+3*k(i,j)+k(i,j-1))*(...
                    A(i,j)*T(i+1,j)+...
                    C(i,j)*T(i-1,j)+...
                    D(i,j)*T(i,j-1));
            else
                T_next(i,j)=...
                    1/(k(i+1,j)+k(i-1,j)+4*k(i,j)+k(i,j+1)+k(i,j-1))*(...
                    A(i,j)*T(i+1,j)+...
                    B(i,j)*T(i,j+1)+...
                    C(i,j)*T(i-1,j)+...
                    D(i,j)*T(i,j-1));
            end
        end
    end
    T_next(1,:) = 0;            % Matrix bottom row, Direchlet BC
    T_next(m,:) = 1-cos(pi*x);  % Matrix top row, Direchlet BC
    T_next(2:m-1,1)= T_next(2:m-1,2); % Matrix Left column, Neumann BC
    T_next(2:m-1,n)= T_next(2:m-1,n-1); % Matrix right column, Neumann BC
    %r_iter = norm(T_next-T,2);
    r_iter = max(max(abs(T_next-T)));
    T = T_next;
    iter=iter+1;
end
fprintf('iterations: %4.1f \n',iter);

%% Plot Figures
figure
contourf(T)
title(['Point Jacobi Method, h =',num2str(h),', iterations: ',num2str(iter)])
xlabel('x points')
ylabel('y points')
colormap hot
colorbar('location','southoutside')

repo: Aero-matlab


text:
%ComputeJointDistribution Computes the joint distribution defined by a set
% of given factors
%
%   Joint = ComputeJointDistribution(F) computes the joint distribution
%   defined by a set of given factors
%
%   Joint is a factor that encapsulates the joint distribution given by F
%   F is a vector of factors (struct array) containing the factors
%     defining the distribution
%

function Joint = ComputeJointDistribution(F)

  % Check for empty factor list
  if (numel(F) == 0)
      warning('Error: empty factor list');
      Joint = struct('var', [], 'card', [], 'val', []);
      return;
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE:
% Compute the joint distribution defined by F
% You may assume that you are given legal CPDs so no input checking is required.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Joint = struct('var', [], 'card', [], 'val', []); % Returns empty factor. Change this.

for i = 1:length(F)
	Joint = FactorProduct(Joint,F(i));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function Dist = par_alldist(x1,x2,varargin)
% a parallelized version of vl_alldist2

opts.measure = 'L2';
opts.maxParts = 100;
opts.maxMem = [];
opts.numWorkers = 20;
opts.verbose = true;
opts = vl_argparse(opts,varargin);

if isempty(x2) || isequal(x1,x2),
    if isempty(x2), x2 = x1; end
    withSelf = true;
else
    withSelf = false;
end

if opts.numWorkers<=1,
    Dist = vl_alldist(x1,x2,opts.measure);
    return;
end

typeX1 = class(x1); typeX2 = class(x2);
if ~strcmp(typeX1,typeX2),
    error('Wrong input: x1 and x2 should be of same data type');
end

[D1, n1] = size(x1); [D2, n2] = size(x2);
if D1~=D2,
    error('Wrong input: x1 and x2 don''t match in 1st dimention');
end;
D = D1;

% see how much memory is available
tmp = x1(1);
info = whos('tmp');
numbytes = info.bytes;
if ~ispc
    [~,w] = unix('free | grep Mem');
    stats = str2double(regexp(w, '[0-9]*', 'match'));
    freemem = (stats(3)+stats(end))*1e3;
else
    stats = memory;
    freemem = stats.MemAvailableAllArrays;
end
if isempty(opts.maxMem),
    maxmem = freemem * 0.8;
else
    maxmem = min(opts.maxMem,freemem);
end
maxNum = maxmem / numbytes;

% decide window width w.r.t. constraints
w = ceil(max(sqrt(n1*n2/opts.maxParts),2*n1*n2*D/(maxNum-2*n1*n2)));
npar = ceil([n1 n2]/w);
while npar(1)*npar(2)>opts.maxParts,
    npar = max(1,npar - 1);
    w = ceil(max([n1 n2]./npar));
    npar = ceil([n1 n2]/w);
end
sz1 = [ones(1,npar(1)-1)*w n1-(npar(1)-1)*w];
sz2 = [ones(1,npar(2)-1)*w n2-(npar(2)-1)*w];

% partition data
t = npar(1)*npar(2);
x1cell = cell(1,t);
x2cell = cell(1,t);
distcell = cell(1,t);
if opts.verbose,
    ts = tic;
    fprintf('[1/3] Partitioning into %dx ((%d,%d),(%d,%d)) parts \n', ...
        t, D,w,D,w);
end
for i = 1:t,
    i1 = mod(i-1,npar(1))+1;
    i2 = floor((i-1)/npar(1))+1;
    if withSelf && i2>i1, % skip top-right triangle for self-dist
        x1cell{i} = [];
        x2cell{i} = [];
    else
        x1cell{i} = x1(:,(i1-1)*w+(1:sz1(i1)));
        x2cell{i} = x2(:,(i2-1)*w+(1:sz2(i2)));
    end
    if opts.verbose,
        if mod(i,10)==0, fprintf('.');end;
        if mod(i,200)==0, fprintf(' [%d/%d]\n',i,t); end;
    end
end
if opts.verbose, fprintf(' done! (%s)\n', timestr(toc(ts))); end;

% estimate speed
tmp = rand(1000,1000);
tt=tic;vl_alldist2(tmp,tmp,opts.measure);tc=toc(tt);
if withSelf,
    estTime = (D*(n1*n2/2-n1/2)/1e9)*tc;
else
    estTime = (D*n1*n2/1e9)*tc;
end

% real work
%{-
% comment this block for ealier MATLAB versions
pool = gcp('nocreate');
if isempty(pool) || pool.NumWorkers<opts.numWorkers,
    if ~isempty(pool), delete(pool); end
    pool = parpool(opts.numWorkers);
end
if opts.verbose,
    ts = tic;
    fprintf('[2/3] Computing distances using %d workers (~%s) ...', ...
    pool.NumWorkers, timestr(estTime/pool.NumWorkers));
end
%}
fprintf('\n');
parfor_progress(t);
parfor i=1:t,
    distcell{i} = vl_alldist2(x1cell{i},x2cell{i},opts.measure);
    parfor_progress();
end
parfor_progress(0);
if opts.verbose, fprintf(' done! (%s)\n', timestr(toc(ts))); end;

% assemble results
if opts.verbose,
    ts = tic;
    fprintf('[3/3] Assembling distance matrix \n');
end;
Dist = zeros(n1,n2,typeX1);
for i = 1:t,
    i1 = mod(i-1,npar(1))+1;
    i2 = floor((i-1)/npar(1))+1;
    if withSelf && i2>i1,
        dist_block = distcell{(i1-1)*npar(1)+i2}';
    else
        dist_block = distcell{i};
    end
    Dist((i1-1)*w+(1:sz1(i1)),(i2-1)*w+(1:sz2(i2))) = dist_block;
    if opts.verbose,
        if mod(i,10)==0, fprintf('.'); end
        if mod(i,200)==0, fprintf(' [%d/%d]\n',i,t); end;
    end
end
if opts.verbose, fprintf(' done! (%s)\n', timestr(toc(ts))); end;

end

function str=timestr(nsecs)

muls = [60 60];
% label = {'H','M','S'};
label = {'','',''};
tvec = [nsecs];

for i=1:numel(muls)
    tvec = [floor(tvec(1)/muls(end-i+1)) tvec];
    tvec(2) = mod(tvec(2),muls(end-i+1));
end

if tvec(end)~=floor(tvec(end))
    str = sprintf('%02.2f%s',tvec(end),label{end});
else
    str = sprintf('%02d%s',tvec(end),label{end});
end
for i=2:numel(tvec),
    str = sprintf('%02d%s:%s',tvec(end-i+1),label{end-i+1},str);
end

end

author: suhangpro

repo: mvcnn

save_path: github-repos/MATLAB/suhangpro-mvcnn


text:
function value = f_01_2d ( dim_num, x )

%*****************************************************************************80
%
%% F_01_2D is the 2D test function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994,
%    ISBN: 0198534728,
%    LC: QA311.S56
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real X(DIM_NUM), the point where the function
%    is to be evaluated.
%
%    Output, real VALUE, the value of the function at X.
%
  e = 2.718281828459045;
  value = x(2) * exp ( x(1) * x(2) ) / ( e - 2.0 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function H = decode_nifti1(blob)
% function H = decode_nifti1(blob)
%
% Decodes a NIFTI-1 header given as raw 348 bytes (uint8) into a Matlab structure
% that matches the C struct defined in nifti1.h, with the only difference that the
% variable length arrays "dim" and "pixdim" are cut off to the right size, e.g., the
% "dim" entry will only contain the relevant elements:
% dim[0..7]={3,64,64,18,x,x,x,x} in C would become dim=[64,64,18] in Matlab.
%
% WARNING: This function currently ignores endianness !!!

% (C) 2010 S.Klanke

if class(blob)~='uint8'
	error 'Bad type for blob'
end
if length(blob)~=348
	error 'Blob must be exactly 348 bytes long'
end

% see nift1.h for information on structure
H = [];

magic = char(blob(345:347));
if blob(348)~=0 | magic~='ni1' & magic~='n+1'
	error 'Not a NIFTI-1 header!';
end

H.sizeof_hdr = typecast(blob(1:4),'int32');
H.data_type  = cstr2matlab(blob(5:14));
H.db_name    = cstr2matlab(blob(15:32));
H.extents    = typecast(blob(33:36),'int32');
H.session_error = typecast(blob(37:38),'int16');
H.regular    = blob(39);
H.dim_info   = blob(40);

dim = typecast(blob(41:56),'int16');
H.dim = dim(2:dim(1)+1);
H.intent_p1   = typecast(blob(57:60),'single');
H.intent_p2   = typecast(blob(61:64),'single');
H.intent_p3   = typecast(blob(65:68),'single');
H.intent_code = typecast(blob(69:70),'int16');
H.datatype    = typecast(blob(71:72),'int16');
H.bitpix      = typecast(blob(73:74),'int16');
H.slice_start = typecast(blob(75:76),'int16');
pixdim = typecast(blob(77:108),'single');
H.qfac        = pixdim(1);
H.pixdim      = pixdim(2:dim(1)+1);
H.vox_offset  = typecast(blob(109:112),'single');
H.scl_scope   = typecast(blob(113:116),'single');
H.scl_inter   = typecast(blob(117:120),'single');
H.slice_end   = typecast(blob(121:122),'int16');
H.slice_code  = blob(123);
H.xyzt_units  = blob(124);
H.cal_max     = typecast(blob(125:128),'single');
H.cal_min     = typecast(blob(129:132),'single');
H.slice_duration = typecast(blob(133:136),'single');
H.toffset     = typecast(blob(137:140),'single');
H.glmax       = typecast(blob(141:144),'int32');
H.glmin       = typecast(blob(145:148),'int32');
H.descrip     = cstr2matlab(blob(149:228));
H.aux_file    = cstr2matlab(blob(229:252));
H.qform_code  = typecast(blob(253:254),'int16');
H.sform_code  = typecast(blob(255:256),'int16');
quats = typecast(blob(257:280),'single');
H.quatern_b = quats(1);
H.quatern_c = quats(2);
H.quatern_d = quats(3);
H.quatern_x = quats(4);
H.quatern_y = quats(5);
H.quatern_z = quats(6);
trafo = typecast(blob(281:328),'single');
H.srow_x = trafo(1:4);
H.srow_y = trafo(5:8);
H.srow_z = trafo(9:12);
%H.S = [H.srow_x; H.srow_y; H.srow_z; 0 0 0 1];
H.intent_name = cstr2matlab(blob(329:344));
H.magic = magic;

function ms = cstr2matlab(cs)
if cs(1)==0
  ms = '';
else
  ind = find(cs==0);
  if isempty(ind)
    ms = char(cs)';
  else
    ms = char(cs(1:ind(1)-1))';
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function bern = bpab ( n, x, a, b )

%*****************************************************************************80
%
%% BPAB evaluates at X the Bernstein polynomials based in [A,B].
%
%  Formula:
%
%    BERN(N,I)(X) = [N!/(I!*(N-I)!)] * (B-X)^(N-I) * (X-A)^I / (B-A)^N
%
%  First values:
%
%    B(0,0)(X) =   1
%
%    B(1,0)(X) = (      B-X              ) / (B-A)
%    B(1,1)(X) = (                X-A    ) / (B-A)
%
%    B(2,0)(X) = (     (B-X)^2           ) / (B-A)^2
%    B(2,1)(X) = ( 2 * (B-X)   * (X-A)   ) / (B-A)^2
%    B(2,2)(X) = (               (X-A)^2 ) / (B-A)^2
%
%    B(3,0)(X) = (     (B-X)^3           ) / (B-A)^3
%    B(3,1)(X) = ( 3 * (B-X)^2 * (X-A)   ) / (B-A)^3
%    B(3,2)(X) = ( 3 * (B-X)   * (X-A)^2 ) / (B-A)^3
%    B(3,3)(X) = (               (X-A)^3 ) / (B-A)^3
%
%    B(4,0)(X) = (     (B-X)^4           ) / (B-A)^4
%    B(4,1)(X) = ( 4 * (B-X)^3 * (X-A)   ) / (B-A)^4
%    B(4,2)(X) = ( 6 * (B-X)^2 * (X-A)^2 ) / (B-A)^4
%    B(4,3)(X) = ( 4 * (B-X)   * (X-A)^3 ) / (B-A)^4
%    B(4,4)(X) = (               (X-A)^4 ) / (B-A)^4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the degree of the Bernstein polynomials to be used.
%    For any N, there is a set of N+1 Bernstein polynomials, each of
%    degree N, which form a basis for polynomials on [A,B].
%
%    Input, real X, the point at which the polynomials are to be evaluated.
%
%    Input, real A, B, the endpoints of the interval on which the
%    polynomials are to be based.  A and B should not be equal.
%
%    Output, real BERN(1:N+1), the values of the N+1 Bernstein polynomials at X.
%
  if ( b == a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BPAB - Fatal error!\n' );
    fprintf ( 1, '  A = B = %f\n', a );
    error ( 'BPAB - Fatal error!' );
  end

  if ( n == 0 )

    bern(1) = 1.0;

  elseif ( 0 < n )

    bern(1) = ( b - x ) / ( b - a );
    bern(2) = ( x - a ) / ( b - a );

    for i = 2 : n
      bern(i+1) = ( x - a ) * bern(i) / ( b - a );
      for j = i-1 : -1 : 1
        bern(j+1) = ( ( b - x ) * bern(j+1) + ( x - a ) * bern(j) ) / ( b - a );
      end
      bern(1) = ( b - x ) * bern(1) / ( b - a );
    end

  end

  return
end


text:
function y = sfb3D(lo, hi, sf1, sf2, sf3)

% 3D Synthesis Filter Bank
%
% USAGE:
%   y = sfb3D(lo, hi, sf1, sf2, sf3);
% INPUT:
%   lo, hi - lowpass subbands
%   sfi - synthesis filters for dimension i
% OUPUT:
%   y - output array
% See afb3D
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

if nargin < 4
   sf2 = sf1;
   sf3 = sf1;
end

LLL = lo;
LLH = hi{1};
LHL = hi{2};
LHH = hi{3};
HLL = hi{4};
HLH = hi{5};
HHL = hi{6};
HHH = hi{7};

% filter along dimension 3
LL = sfb3D_A(LLL, LLH, sf3, 3);
LH = sfb3D_A(LHL, LHH, sf3, 3);
HL = sfb3D_A(HLL, HLH, sf3, 3);
HH = sfb3D_A(HHL, HHH, sf3, 3);

% filter along dimension 3
L = sfb3D_A(LL, LH, sf2, 2);
H = sfb3D_A(HL, HH, sf2, 2);

% filter along dimension 1
y = sfb3D_A(L, H, sf1, 1);

% LOCAL FUNCTION

function y = sfb3D_A(lo, hi, sf, d)

% 3D Synthesis Filter Bank
% (along single dimension only)
%
% y = sfb3D_A(lo, hi, sf, d);
% sf - synthesis filters
% d  - dimension of filtering
% see afb2D_A

lpf = sf(:, 1);     % lowpass filter
hpf = sf(:, 2);     % highpass filter

% permute dimensions of lo and hi so that dimension d is first.
p = mod(d-1+[0:2], 3) + 1;
lo = permute(lo, p);
hi = permute(hi, p);

[N1, N2, N3] = size(lo);
N = 2*N1;
L = length(sf);
y = zeros(N+L-2, N2, N3);

for k = 1:N3
   y(:, :, k) = upfirdn(lo(:, :, k), lpf, 2, 1) + upfirdn(hi(:, :, k), hpf, 2, 1);
end
y(1:L-2, :, :) = y(1:L-2, :, :) + y(N+[1:L-2], :, :);
y = y(1:N, :, :);
y = cshift3D(y, 1-L/2, 1);

% permute dimensions of y (inverse permutation)
y = ipermute(y, p);

author: thomaskuestner


text:
function disp( x, name )
    %DISP Display TT/MPS tensor.
    %   disp(X) displays the TT/MPS tensor X as a tensor network.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    if (nargin < 2 || ~ischar(name))
      name = inputname(1);
    end

    disp([name, ' is a TT/MPS tensor of order ', num2str(x.order), ...
                        ' with size (', num2str(x.size), ...
                        ') and ranks (', num2str(x.rank), ')']);
    disp('');

    row1 = '';
    row2 = '';
    row3 = '';

    for i=1:x.order
        row1 = [row1, sprintf('%3i--(U%2i)--', x.rank(i), i)];
        row2 = [row2, '       |    '];
        row3 = [row3, sprintf('     %3i    ', x.size(i))];
    end
    row1 = [row1, sprintf( '%3i', x.rank(end) )];
    disp(row1)
    disp(row2)
    disp(row3)

end

author: NicolasBoumal

repo: manopt


text:
function [PA, P, PS] = reduce(disc, A, S)
%REDUCE   Dimension reduction for operator matrix.
%   PA = REDUCE(DISC, A) reduces the row dimension of each block column in the
%   cell array A (which is typically a discretization of DISC.SOURCE) so that
%   the reduced discretization, PA, can be formed as a matrix. In particular, PA
%   will have sum(DISC.dimension) rows and sum(cellfun(@(a) size(A, 2), A(1,:))
%   columns.
%
%   [PA, P] = REDUCE(DISC, A) returns also the block-diagonal reduction matrix
%   P. For COLLOC discretizations, this blocks are BARYMAT projections.
%
%   [PA, P, PS] = REDUCE(DISC, A, S) is required for consistency with other
%   opDiscretization reductions. Here S is ignored and PS = P.

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

% Setup:
r = disc.projOrder;
dim = disc.dimension;
dimAdjust = disc.dimAdjust(1,:);
if ( numel(dimAdjust) == 1 )
    dimAdjust = repmat(dimAdjust, 1, size(A, 2));
end
PA = cell(1, size(A, 2));
P = cell(1, size(A, 2));

% Do reduction for each block column:
for k = 1:size(A, 2)
    [PA{k}, P{k}] = reduceOne(disc, A(:,k), r(k), dim + dimAdjust(k));
%     [PS{k}, ignored] = reduceOne(disc, S(:,k), r(k), dim + dimAdjust(k));
end

% Convert cell arrays to matrices:
P = blkdiag(P{:});
PA = cell2mat(PA);
% PS = cell2mat(PS);
PS = P;

end

function [PA, P] = reduceOne(disc, A, m, n)
%REDUCEONE   Reduce one block column.
%   [PA, P] = REDUCEONE(DISC, A, M, N) reduces entries of the column cell arrays
%   A from a sum(N)xsum(N) discretization to sum(N-M)xsum(N) version, PA, using
%   the block-projection operator P.

% Step by intervals in the domain.
domain = disc.domain;
numInt = disc.numIntervals;
P = cell(1, numInt);

% Loop through intervals
for k = 1:numInt
    disc.domain = domain(k:(k+1));
    disc.dimension = n(k)-m;
    [xOut, ignored, ignored, tOut] = equationPoints(disc);
    disc.dimension = n(k);
    [xIn, ignored, baryWt, tIn] = functionPoints(disc);
    % Store the kth projection matrix in the cell P
%     P{k} = barymat(xOut, xIn, baryWt);
    P{k} = barymat(xOut, xIn, baryWt, tOut, tIn, 1);
end
% Convert the projection matrices P into a blockdiagonal matrix.
P = blkdiag(P{:});

% Project each of the entries of A:
PA = cell(size(A));
for j = 1:numel(A)
    if ( size(P, 2) == size(A{j}, 1) )
        PA{j} = P*A{j};
    else
        PA{j} = A{j};
    end
end
PA = cell2mat(PA);

if ( (m == 0) && (size(A{1}, 2) < sum(n)) )
    % We don't want to project scalars.
    P = eye(size(A, 2));
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%MDL_OFFSET3 A minimalistic 3DOF robot arm with shoulder offset
%
% MDL_OFFSET3 is a script that creates the workspace variable off3 which
% describes the kinematic characteristics of a simple arm manipulator with
% a shoulder offset, using standard DH conventions.
%
% Somewhat like a Puma arm without the wrist.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%
% Notes::
% - Unlike most other mdl_xxx scripts this one is actually a function that
%   behaves like a script and writes to the global workspace.
%
% See also SerialLink, mdl_offset6, mdl_simple6, mdl_puma560.

% MODEL: generic, 3DOF, standard_DH

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function r = mdl_offset3()

    % robot length values (metres)
    L1 = 1;
    L2 = 1;
    O1 = 0;
    O2 = 0.2;
    O3 = 0;

    % and build a serial link manipulator

    robot = SerialLink([
        Revolute('alpha', -pi/2, 'a', O1, 'd', 0)
        Revolute('alpha', 0,     'a', L1, 'd', O2)
        Revolute('alpha', pi/2,  'a', L2, 'd', O3)
        ], ...
        'name', 'Offset3');

    % place the variables into the global workspace
    if nargout == 1
        r = robot;
    elseif nargout == 0
        assignin('base', 'off3', robot);
        assignin('base', 'qz', [0 0 0 ]); % zero angles, arm up
    end
end

author: Allopart

repo: rbpf-gmapping


text:
function Z = zeros(d)
%ZEROS     Zero operator.
%   This function is deprecated. Use OPERATORBLOCK.ZEROS.

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

Z = linop( operatorBlock.zeros(double(d)) );

end


text:
function [Y] = spm_sdot(X,x,k)
% Sparse multidimensional dot (inner) product
% FORMAT [Y] = spm_sdot(X,x,[DIM])
%
% X   - numeric array
% x   - cell array of numeric vectors
% DIM - dimension to omit (asumes ndims(X) = numel(x))
%
% Y  - inner product obtained by summing the products of X and x along DIM
%
% If DIM is not specified the leading dimensions of X are omitted. This
% routine assumes X is sparse
%
% See also: spm_dot, spm_cross
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_sdot.m 7300 2018-04-25 21:14:07Z karl $

% initialise dimensions
%--------------------------------------------------------------------------
DIMS   = size(X);
NDIM   = numel(DIMS);
J      = 1:NDIM;
J(k)   = [];
I      = cell(1,NDIM);
j      = find(X > exp(-16));
[I{:}] = ind2sub(DIMS,j);

% sum of products
%--------------------------------------------------------------------------
Y      = zeros(numel(x{k}),1);
for i  = 1:numel(j)
    p  = X(j(i));
    for d = J
        p = p*x{d}(I{d}(i));
    end
    Y(I{k}(i)) = Y(I{k}(i)) + p;
end

author: spm


text:
% SLAMTBSLAVE  A slave version of SLAMTB.
%
%   SLAMTBSLAVE is intended to be invoked by other scripts initializing
%   some of the parameters. By now it is NEESANALYSIS the only one to do
%   so.
%
%   See also NEESANALYSIS, NEESPLOTS.

%   Created and maintained by
%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.
%   Programmers:
%   Copyright David Marquez and Jean-Marie Codol @ LAAS-CNRS

%% OK we start here

% clear workspace and declare globals
%clear
global Map          %#ok<NUSED>

%% I. Specify user-defined options - EDIT USER DATA FILE userData.m

userDataPnt;        % user-defined data. SCRIPT.
% userDataLin;        % user-defined data. SCRIPT.

% These inputs come from neesAnalysis and overwrite some userData values:
Time.lastFrame            = numFrames;
Opt.init.initType         = lmkType;
Opt.correct.reparametrize = reparametrize;
SimOpt.random.newSeed     = false;
SimOpt.random.fixedSeed   = randSeeds(nRun);
FigOpt.rendPeriod         = rendPeriod;
Sensor{1}.pixErrorStd     = pixelNoise;

%% II. Initialize all data structures from user-defined data in userData.m
% SLAM data
[Rob,Sen,Raw,Lmk,Obs,Tim]     = createSlamStructures(...
    Robot,...
    Sensor,...      % all user data
    Time,...
    Opt);

% Simulation data
[SimRob,SimSen,SimLmk,SimOpt] = createSimStructures(...
    Robot,...
    Sensor,...      % all user data
    World,...
    SimOpt);

% Graphics handles
[MapFig,SenFig]               = createGraphicsStructures(...
    Rob, Sen, Lmk, Obs,...      % SLAM data
    SimRob, SimSen, SimLmk,...  % Simulator data
    FigOpt);                    % User-defined graphic options

%% III. Init data logging
% TODO: Create source and/or destination files and paths for data input and
% logs.
% TODO: do something here to collect data for post-processing or
% plotting. Think about collecting data in files using fopen, fwrite,
% etc., instead of creating large Matlab variables for data logging.

logFile = openLogFile(logFileName,'w');
writeLogHeader(logFile,logFileHeader);

% Clear user data - not needed anymore
clear Robot Sensor World Time   % clear all user data

%% IV. Main loop
for currentFrame = Tim.firstFrame : Tim.lastFrame

    % 1. SIMULATION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Simulate robots
    for rob = [SimRob.rob]

        % Robot motion
        SimRob(rob) = simMotion(SimRob(rob),Tim);

        % Simulate sensor observations
        for sen = SimRob(rob).sensors

            % Observe simulated landmarks
            Raw(sen) = simObservation(SimRob(rob), SimSen(sen), SimLmk, SimOpt) ;

        end % end process sensors

    end % end process robots

    % 2. ESTIMATION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Process robots
    for rob = [Rob.rob]

        % Robot motion
        % FIXME: see how to include noise in a clever way.
        Rob(rob).con.u = SimRob(rob).con.u + Rob(rob).con.uStd.*randn(6,1)*noiseFactor;

        Rob(rob) = motion(Rob(rob),Tim);

        % Process sensor observations
        for sen = Rob(rob).sensors

            %TODO: see how to pass only used Lmks and Obs.
            % Observe knowm landmarks
            [Rob(rob),Sen(sen),Lmk,Obs(sen,:)] = correctKnownLmks( ...
                Rob(rob),   ...
                Sen(sen),   ...
                Raw(sen),   ...
                Lmk,        ...
                Obs(sen,:), ...
                Opt) ;

            % Initialize new landmarks
            if currentFrame == Tim.firstFrame
                numInit = Opt.init.nbrInits(1);
            else
                numInit = Opt.init.nbrInits(2);
            end
            for ini = 1:numInit
                [Lmk,Obs(sen,:)] = initNewLmk(...
                    Rob(rob),   ...
                    Sen(sen),   ...
                    Raw(sen),   ...
                    Lmk,        ...
                    Obs(sen,:), ...
                    Opt) ;
            end

        end % end process sensors

    end % end process robots

    % 3. VISUALIZATION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if currentFrame == Tim.firstFrame ...
            || currentFrame == Tim.lastFrame ...
            || mod(currentFrame,FigOpt.rendPeriod) == 0

        % Figure of the Map:
        MapFig = drawMapFig(MapFig,  ...
            Rob, Sen, Lmk,  ...
            SimRob, SimSen, ...
            FigOpt);

        makeVideoFrame(MapFig, ...
            sprintf('map-%04d.png',currentFrame), ...
            FigOpt, ExpOpt);

        % Figures for all sensors
        for sen = [Sen.sen]
            SenFig(sen) = drawSenFig(SenFig(sen), ...
                Sen(sen), Raw(sen), Obs(sen,:), ...
                FigOpt);

            makeVideoFrame(SenFig(sen), ...
                sprintf('sen%02d-%04d.png', sen, currentFrame),...
                FigOpt, ExpOpt);

        end

        % Do draw all objects
        drawnow;
    end

    % 4. DATA LOGGING
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % TODO: do something here to collect data for post-processing or
    % plotting. Think about collecting data in files using fopen, fwrite,
    % etc., instead of creating large Matlab variables for data logging.

    % NEES for robot 1
    logData = errorAnalysis(Rob(1), SimRob(1), errorAnalysisFunction);
%     logData = robotNees(Rob(1),SimRob(1));
    writeLogData(logFile,[currentFrame,logData]);

end

%% V. Post-processing

closeLogFile(logFile);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function a = plu_inverse ( n, pivot )

%*****************************************************************************80
%
%% PLU_INVERSE returns the inverse of a PLU matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer PIVOT(N), the list of pivot rows.  PIVOT(I)
%    must be a value between I and N, reflecting the choice of
%    pivot row on the I-th step.  For no pivoting, set PIVOT(I) = I.
%
%    Output, real A(N,N), the inverse matrix.
%
  [ p, l, u ] = plu_plu ( n, pivot );

  p_inverse = p';

  l_inverse = tri_l1_inverse ( n, l );

  u_inverse = tri_u_inverse ( n, u );

  a = u_inverse * l_inverse * p_inverse;

  return
end


text:
%  FILE:   cluster_rects.m
%
%    This function derive canonical bounding box shapes by applying K-medoid
%    clustering on a set of bounding boxes.
%
%  INPUT:  imdb        (dataset)
%          N           (number of medoids)
%          minsz/maxsz (filter out boxes that are either too big or too small)
%          vis         (whether we visualize the clustering results)
%
%  OUTPUT: C           (N medoids)

function C  = cluster_rects(imdb, N, minsz, maxsz, vis)
if nargin < 5
  vis = 0;
end

%% cluster based on shapes of training examples
idx = find(imdb.images.set == 1);
rects = imdb.labels.rects(idx);
rects = vertcat(rects{:});

%% centralize
hs = rects(:,4) - rects(:,2) + 1;
ws = rects(:,3) - rects(:,1) + 1;
rects = [-(ws-1)/2, -(hs-1)/2, (ws-1)/2, (hs-1)/2];

%% ignore faces with size out of the range
idx = find(hs<=maxsz(1)&ws<=maxsz(2)&hs>=minsz(1)&ws>=minsz(2));
rects = rects(idx,:);
fprintf('Ignored faces smaller than %dx%d, %d bboxes left.\n',minsz(1),minsz(2),numel(idx));

%% subsample for faster clustering
rects = rects(randsample(size(rects,1), min(size(rects,1), 1e5)), :);
fprintf('Clustering on %d/%d face bounding boxes.\n', size(rects,1), numel(idx));

%% build kmedoids
[Cidx,C,sumd,D,midx] = kmedoids(rects, N, 'Options', statset('UseParallel', true), 'Distance', @rect_dist);

%% reorder clusters based on bounding box areas
[~,I] = sort(C(:,3).*C(:,4),'descend');
C = C(I,:);

if ~vis, return; end
subplot = @(m,n,k) subtightplot(m,n,k,[0.1,0.1]);
clf;
[SI,SJ] = factorize(N);
for i = 1:N
  subplot(SI,SJ,i);
  plotBoxes(C(i,1),C(i,2),C(i,3)-C(i,1)+1,C(i,4)-C(i,2)+1,rand(1,3),0.5);
  title(num2str(i));
  axis([-250,250,-250,250]);
end

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
% line_fewer_markers - line with controlled amount of markers and correct legend behaviour
%
%   LINE_FEWER_MARKERS(X,Y,NUM_MARKERS) adds the line in vectors X and Y to the current axes
%   with exactly NUM_MARKERS markers drawn.
%
%   LINE_FEWER_MARKERS(X,Y,NUM_MARKERS,'PropertyName',PropertyValue,...) plots the data
%   stored in the vectors X and Y.
%
%   LINE_FEWER_MARKERS returns handles to LINE/MARKER objects.
%
%   [H1,H2,H3] = LINE_FEWER_MARKERS(X,Y,NUM_MARKERS,'PropertyName',PropertyValue,...)
%   performs the actions as above and returns the handles of all the plotted lines/markers.
%   H1    = handle to the main marker(1 point); it may be put in array and used with legend
%   H2    = handle to the continuous line (as in H2=plot())
%   H3    = handle to all other markers
%
%   Property/Value pairs and descriptions:
%
%       Spacing    - 'x'      : ordinary uniform along x
%                  - 'curve'  : equal lengths along curve y(x)
%                  - 'logx'   : to be used with logarithmic x scale
%
%       LockOnMax  - 0        : first marker on 1st data point
%                  - 1        : offset all markers such that one marker on first max of y(x)
%
%       LegendLine - 'on'     : default, reproduce linestyle also in legend
%                  - 'off'    : shows only marker in legend
%
%       LineSpec: same as for LINE: LineStyle,LineWidth,Marker,MarkerSize,MarkerFaceColor...
%
%
%   Example: plot 3 curves with 9,9, and 15 markers each, using different input styles
%
%      figure; hold on;
%      t  = 0:0.005:pi;
%      line_fewer_markers(t*180/pi,cos(t) ,9,  '--bs','spacing','curve');
%      line_fewer_markers(t*180/pi,sin(t) ,9,  '-.ro','MarkerFaceColor','g', ...
%                                                     'markersize',6,'linewidth',2);
%      grey1 = [1 1 1]*0.5;
%      line_fewer_markers(t*180/pi,sin(t).*cos(t) ,15, ':','marker','h','color',grey1, ...
%                                    'markerfacecolor',grey1,'linewidth',2,'LockOnMax',1);
%      leg = legend('cos','sin','sin*cos','location','best');
%
% Inspired by Ioannis Filippidis's answer:
% http://www.mathworks.com/matlabcentral/answers/2165-too-many-markers
%
% rev.4, Massimo Ciacci, October 17, 2014
%
function [H1,H2,H3] = line_fewer_markers(x,y,num_Markers, varargin)

%% find marker spec in varargin and remove it; extract special params: LockOnMax,Spacing
if mod(length(varargin),2)
    if ischar(varargin{1})
      linspec   = varargin{1};
      extraArgs = varargin(2:end);
      [varargInNoMk,varargInNoMkNoLn,lm,ms,mfc,LockOnMax,Spacing,LegendLine] = parseargsLineSpec(linspec,extraArgs);
    else
      error('odd sized [param | val] list, missing one param ?');
    end
else
      [varargInNoMk,varargInNoMkNoLn,lm,ms,mfc,LockOnMax,Spacing,LegendLine] = parseargs(varargin{:});
end

%% input size check
if  isvector(x) &&  isvector(y)
    % make x,y row vectors
    if iscolumn(x), x = x.'; end
    if iscolumn(y), y = y.'; end
else
    error('line_fewer_markers: input arguments must be 1D vectors');
end

% How the method works: plots 3 times:
% a) once only the line with all points with the style                                        'r--' and invisible handle,
% b) last time the markers, using fewer points with style                                     'ro' and again invisible handle.
% c) once with a visible handle, only the first point, using the complete style you specified (e.g. 'r--o')

%% a) once only the line with all points with the style
H2 = line(x   ,y   ,varargInNoMk{:});                               %no markers here
hasbehavior(H2,'legend',0);                                         %prevent to appear in legends!

%% b) last time the markers, using fewer points with style
%sort the xvalues
[x_t,I] = sort(x);
y_t = y(I);
[x_t,I] = unique(x_t);
y_t = y_t(I);

if     (strcmp(Spacing,'x') || strcmp(Spacing,'X'))
    ti = round(linspace(1,length(x_t),num_Markers));
elseif (strcmp(Spacing,'logx') || strcmp(Spacing,'log'))
    xi = logspace(log10(x_t(2)),log10(x_t(end-1)),num_Markers);
    ti = floor(interp1(x_t,(1:length(x_t)),xi));
elseif (strcmp(Spacing,'curve') || strcmp(Spacing,'Curve'))
    scaleY     = 3/4; % 1/1 figure aspect ratio
    yNrm       = (y-min(y))./(max(y)-min(y))*scaleY;             %NORMALIZE y scale in [0 1], height of display is prop to max(abs(y))
    xNrm       = (x-min(x))./(max(x)-min(x));                    %NORMALIZE x scale in [0 1]

    if (sum(isinf(yNrm))>0) || sum(isinf(x))>0                   %spacing along curve not possible with infinites
        ti = round(linspace(1,length(x),num_Markers));
    else
        t        = 1:length(x);
        s        = [0 cumsum(sqrt(diff(xNrm).^2+diff(yNrm).^2))];%measures length along the curve
        si       = (0:num_Markers-1)*s(end)/(num_Markers-1);     %equally spaced lengths along the curve
        si(end)  = s(end);                                       %fix last point to be within the curve
        ti       = round(interp1(s,t,si));                       %find x index of markers
    end
else
    error('invalid spacing parameter');
end
%if LockOnMax
    %set one ti on max if found
%    [Mv,idx]   = max(y); idx=idx(1);
%    [mv,idxti] = min(abs(idx-ti));
%    deltati    = ti(idxti)-idx;
%    ti         = max(1,min(ti-deltati,length(y)));
%end
xi = x_t(ti);
yi = y_t(ti);
H3 = line(xi,yi,varargInNoMkNoLn{:},'Marker',lm,'MarkerSize',ms,'MarkerFaceColor',mfc,'LineStyle','none');  %plot markers only
hasbehavior(H3,'legend',0); %prevent to appear in legends!

%% c) once with a visible handle, only the first point, using the complete style you specified
if strcmp(LegendLine,'on')
    H1 = line(xi(1),yi(1),varargInNoMk{:},'Marker',lm,'MarkerSize',ms,'MarkerFaceColor',mfc);
else
    H1 = line(xi(1),yi(1),varargInNoMk{:},'linestyle','none','Marker',lm,'MarkerSize',ms,'MarkerFaceColor',mfc);
end

%-------------------------------------------------------------
% PARSE FUNCTIONS
%-------------------------------------------------------------
% varargInNoMk = list of property pairs, marker specs removed
% varargInNoMkNoLn = list of property pairs, marker specs and line specs removed
function [varargInNoMk,varargInNoMkNoLn,lm,ms,mfc,LockOnMax,Spacing,LegendLine] = parseargs(varargin)
lm =[]; ms =[]; mfc=[]; LockOnMax=[]; Spacing=[]; LegendLine=[];
varargInNoMk = {};
varargInNoMkNoLn = {};
arg_index = 1;
while arg_index <= length(varargin)
	arg = varargin{arg_index};
    % extract special params and marker specs from arg list
    if strcmp(arg,'marker') || strcmp(arg,'Marker') || strcmp(arg,'Mk')  || strcmp(arg,'mk')
        lm              = varargin{arg_index+1};
    elseif strcmp(arg,'MarkerSize') || strcmp(arg,'markersize') || strcmp(arg,'Mks')  || strcmp(arg,'mks')
        ms              = varargin{arg_index+1};
    elseif strcmp(arg,'MarkerFaceColor') || strcmp(arg,'markerfacecolor')||strcmp(arg,'MFC')||strcmp(arg,'mfc')
        mfc             = varargin{arg_index+1};
    elseif strcmp(arg,'LockOnMax') || strcmp(arg,'lockonmax')
        LockOnMax       = varargin{arg_index+1};
    elseif strcmp(arg,'Spacing')   || strcmp(arg,'spacing')
        Spacing         = varargin{arg_index+1};
    elseif strcmp(arg,'LegendLine')   || strcmp(arg,'legendline')
        LegendLine      = varargin{arg_index+1};
    else
        % keep other params in arg list for line command
        varargInNoMk    = {varargInNoMk{:},  varargin{arg_index},  varargin{arg_index+1}};
        if ~strcmp(arg,'LineStyle') && ~strcmp(arg,'linestyle')
           % exclude line params for marker only plot
           varargInNoMkNoLn = {varargInNoMkNoLn{:},  varargin{arg_index},  varargin{arg_index+1}};
        end
    end
    arg_index = arg_index + 2;
end
%EXTRA DEFAULTS ARE SET HERE
if isempty(lm),         lm          = 'o'   ; end
if isempty(ms),         ms          = 10    ; end
if isempty(mfc),        mfc         = 'none'; end
if isempty(LockOnMax),  LockOnMax   = 1     ; end
if isempty(Spacing),    Spacing     = 'x'   ; end %%'x' -> marker delta-x constant; 'curve' : spacing constant along the curve length
if isempty(LegendLine), LegendLine  = 'on'  ; end

%-------------------------------------------------------------
% Parse LineSpec string and other arguments
% varargInNoMk     = list of property pairs, marker specs removed
% varargInNoMkNoLn = list of property pairs, marker specs and line specs removed
function [varargInNoMk,varargInNoMkNoLn,lm,ms,mfc,LockOnMax,Spacing,LegendLine] = parseargsLineSpec(linspec, extraArgs)
%          b     blue          .     point              -     solid
%          g     green         o     circle             :     dotted
%          r     red           x     x-mark             -.    dashdot
%          c     cyan          +     plus               --    dashed
%          m     magenta       *     star             (none)  no line
%          y     yellow        s     square
%          k     black         d     diamond
%          w     white         v     triangle (down)
%                              ^     triangle (up)
%                              <     triangle (left)
%                              >     triangle (right)
%                              p     pentagram
%                              h     hexagram
varargInNoMk     = {};
varargInNoMkNoLn = {};

foundLine           = false;
stringSearch        = {'-.','--','-',':'};
for ii=1:4
    if strfind(linspec, stringSearch{ii})
        foundLine   = true;
        ls          = stringSearch{ii};
        linspec     = setdiff(linspec,ls);
        break
    end
end
if foundLine
    varargInNoMk    = {varargInNoMk{:},'lineStyle',ls};
else
    varargInNoMk    = {varargInNoMk{:},'lineStyle','-'};
end

if ~isempty(linspec)
    foundCol            = false;
    stringSearch        = {'b','g','r','c','m','y','k','w'};
    for ii=1:8
        if strfind(linspec, stringSearch{ii})
            foundCol    = true;
            colspec     = stringSearch{ii};
            linspec     = setdiff(linspec,colspec);
            break
        end
    end
    if foundCol
        varargInNoMk    = {varargInNoMk{:},'color',colspec};
        varargInNoMkNoLn    = {varargInNoMkNoLn{:},'color',colspec};
    end
end

if ~isempty(linspec)
    foundMk             = false;
    stringSearch        = {'.','o','x','+','*','s','d','v','^','<','>','p','h'};
    for ii=1:13
        if strfind(linspec, stringSearch{ii})
            foundMk     = true;
            mkspec      = stringSearch{ii};
            break
        end
    end
    if foundMk, lm = mkspec; else lm = 'none'; end
else
    lm = 'none';
end

[extraArgs1,unused,lm2,ms,mfc,LockOnMax,Spacing,LegendLine] = parseargs(extraArgs{:});
if strcmp(lm,'none') && ~strcmp(lm2,'none') %if other marker specified in Property Pairs take that one
    lm = lm2;
end
varargInNoMk       = {varargInNoMk{:},extraArgs1{:}};
varargInNoMkNoLn   = {varargInNoMkNoLn{:},extraArgs1{:}};

author: konstantinberlin

repo: malware-windows-audit-log-detection


text:
function [C] = Clustering(data, k, varargin)
% Perform hierarchical agglomerative clustering to group points in
% matrix data into k clusters

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Roman Denysiuk

% get size of data
n = size(data, 1);

% compute distances between points
distPoints = zeros(n);
for i = 1:n
    distPoints(i,:) = sqrt(sum(power(repmat(data(i,:),n,1)-data, 2), 2));
    distPoints(i,i) = inf;
end

% initialize distances between clusters
distClusters = distPoints;

% initially, each points belongs to a distinct cluster
C = cell(n, 1);
for i = 1:n;
    C{i} = i;
end

update = 0;
toRemove = zeros(n-k, 1);
for i = 1:n-k

    % update distances between clusters
    if update>0

        for jj = 1:n

            if isempty(C{jj}) || jj == update
                continue
            end

            % calculate distances between clusters
            distClusters(update, jj) = sum(sum(distPoints(C{update}, C{jj})))/(numel(C{update})*numel(C{jj}));

        end

    end

    % find two clusters with min distance
    [a, b] = find(distClusters==min(min(distClusters)));

    % merge these clusters
    C{a(1)} = horzcat(C{a(1)}, C{b(1)});

    % removed cluster
    toRemove(i) = b(1);
    C{toRemove(i)} = [];
    update = a(1);

    % remove from distance matrix
    distClusters(toRemove(i), :) = inf;
    distClusters(:, toRemove(i)) = inf;

end

% final clusters
C(toRemove) = [];
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ mainDirect, score, angle] = findMainDirectionEMA( lines )
%FINDMAINDIRECTION compute vp from set of lines
%   Detailed explanation goes here
fprintf('Computing vanishing point:\n');

% arcList = [];
% for i = 1:length(edge)
%     panoLst = edge(i).panoLst;
%     if size(panoLst,1) == 0
%         continue;
%     end
%     arcList = [arcList; panoLst];
% end

%% initial guess
segNormal = lines(:,1:3);
segLength = lines(:,7);
segScores = ones(size(lines,1),1);%lines(:,8);

shortSegValid = segLength < 5*pi/180;
segNormal = segNormal(~shortSegValid,:);
segLength = segLength(~shortSegValid);
segScores = segScores(~shortSegValid);

numLinesg = size(segNormal,1);
[candiSet, tri] = icosahedron2sphere(3);
ang = acos(dot(candiSet(tri(1,1),:), candiSet(tri(1,2),:), 2)) / pi * 180;
binRadius = ang/2;
[ initXYZ, score, angle] = sphereHoughVote( segNormal, segLength, segScores, 2*binRadius, 2, candiSet );

if isempty(initXYZ)
    fprintf('Initial Failed\n');
    mainDirect = [];
    return;
end

fprintf('Initial Computation: %d candidates, %d line segments\n', size(candiSet,1), numLinesg);
fprintf('direction 1: %f %f %f\ndirection 2: %f %f %f\ndirection 3: %f %f %f\n', ...
        initXYZ(1,1), initXYZ(1,2), initXYZ(1,3), ...
        initXYZ(2,1), initXYZ(2,2), initXYZ(2,3), ...
        initXYZ(3,1), initXYZ(3,2), initXYZ(3,3));
%% iterative refine
iter_max = 3;
[candiSet, tri] = icosahedron2sphere(5);
numCandi = size(candiSet,1);
angD = acos(dot(candiSet(tri(1,1),:), candiSet(tri(1,2),:), 2)) / pi * 180;
binRadiusD = angD/2;
curXYZ = initXYZ;
tol = linspace(4*binRadius, 4*binRadiusD, iter_max); % shrink down #ls and #candi
for iter = 1:iter_max
    dot1 = abs(dot( segNormal, repmat(curXYZ(1,:), [numLinesg 1]), 2));
    dot2 = abs(dot( segNormal, repmat(curXYZ(2,:), [numLinesg 1]), 2));
    dot3 = abs(dot( segNormal, repmat(curXYZ(3,:), [numLinesg 1]), 2));
    valid1 = dot1<cos((90-tol(iter))*pi/180);
    valid2 = dot2<cos((90-tol(iter))*pi/180);
    valid3 = dot3<cos((90-tol(iter))*pi/180);
    valid = valid1 | valid2 | valid3;

    if(sum(valid)==0)
        fprintf('ZERO line segment for voting\n');
        break;
    end

    subSegNormal = segNormal(valid,:);
    subSegLength = segLength(valid);
    subSegScores = segScores(valid);

    dot1 = abs(dot( candiSet, repmat(curXYZ(1,:), [numCandi 1]), 2));
    dot2 = abs(dot( candiSet, repmat(curXYZ(2,:), [numCandi 1]), 2));
    dot3 = abs(dot( candiSet, repmat(curXYZ(3,:), [numCandi 1]), 2));
    valid1 = dot1>cos(tol(iter)*pi/180);
    valid2 = dot2>cos(tol(iter)*pi/180);
    valid3 = dot3>cos(tol(iter)*pi/180);
    valid = valid1 | valid2 | valid3;

    if(sum(valid)==0)
        fprintf('ZERO candidate for voting\n');
        break;
    end

    subCandiSet = candiSet(valid,:);

    [ tcurXYZ ] = sphereHoughVote( subSegNormal, subSegLength, subSegScores, 2*binRadiusD, 2, subCandiSet );

    if(isempty(tcurXYZ))
        fprintf('NO answer found!\n');
        break;
    end
    curXYZ = tcurXYZ;

    fprintf('%d-th iteration: %d candidates, %d line segments\n', iter, size(subCandiSet,1), length(subSegScores));

end
fprintf('direction 1: %f %f %f\ndirection 2: %f %f %f\ndirection 3: %f %f %f\n', ...
        curXYZ(1,1), curXYZ(1,2), curXYZ(1,3), ...
        curXYZ(2,1), curXYZ(2,2), curXYZ(2,3), ...
        curXYZ(3,1), curXYZ(3,2), curXYZ(3,3));
mainDirect = curXYZ;

mainDirect(1,:) = mainDirect(1,:).*sign(mainDirect(1,3));
mainDirect(2,:) = mainDirect(2,:).*sign(mainDirect(2,3));
mainDirect(3,:) = mainDirect(3,:).*sign(mainDirect(3,3));

uv = xyz2uvN(mainDirect);
[~,I1] = max(uv(:,2));
J = setdiff(1:3, I1);
[~,I2] = min(abs(sin(uv(J,1))));
I2 = J(I2);
I3 = setdiff(1:3, [I1 I2]);
mainDirect = [mainDirect(I1,:); mainDirect(I2,:); mainDirect(I3,:)];

mainDirect(1,:) = mainDirect(1,:)*sign(mainDirect(1,3));
mainDirect(2,:) = mainDirect(2,:)*sign(mainDirect(2,2));
mainDirect(3,:) = mainDirect(3,:)*sign(mainDirect(3,1));

mainDirect = [mainDirect; -mainDirect];

% score = 0;

end

author: zouchuhang

repo: LayoutNet


text:
function lambda = bis_eigenvalues ( alpha, beta, n )

%*****************************************************************************80
%
%% BIS_EIGENVALUES returns the eigenvalues of the BIS matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, BETA, the scalars which define the
%    diagonal and first superdiagonal of the matrix.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real LAMBDA(N,1), the eigenvalues of the matrix.
%
  lambda(1:n,1) = alpha;

  return
end


text:
function y=triangle_wave(x)

p=2.*pi;
y=(4.*(abs((x./p)-floor((x./p)+0.5))))-1;

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [m,v]=v_psycdigit(proc,r,mode,p,q,xp,noise,fn,dfile,ofile)
%V_PSYCDIGIT measures psychometric function using TIDIGITS stimuli
%
% Usage:
%         (1)[m,v]=v_psycdigit([],[],'GMPWrn10',[],[],[],[],[],dfile);
%                       % measure SRT using addditive white noise, repetitions allowed, data in .WAV format
%         (2)[m,v]=v_psycdigit(@v_specsub,[],'GMPn10',[],[],[],[],[],dfile);
%                       % compare spectral subtraction with unprocessed noisy speech
%
% Inputs:
%         proc    processing function handle (e.g. @v_specsub) called as follows:
%                    y=proc(x,fs,i,r)  % process noisy waveform x through model i with parameters r
%                    y=proc(x,fs,snr,i,r) % process clean speech degraded to snr through model i with parameters r
%                    y=proc()   % return comment text string describing algorithm (if 'c' option specified)
%         r       parameters for proc (omitted from call if r=[])
%         mode    string containing options (see below for list)
%         p,q,xp  parameters passed on to v_psycest or v_psycestu
%         noise   noise waveform or wav file containing noise
%         fn      noise waveform sample frequency [16000]
%         dfile   path to TIdigits folder
%         ofile   output text file (see below for output file format)
%
% Outputs:
%          m(2,n,3) estimated srt and slope of all models
%                   m(i,n,k): i={1=srt, 2=slope}, n=model number, k={1=mean, 2=mode (MAP), 3=marginal mode}
%          v(3,n)   estimated covariance matrix entries:
%                   [var(srt) cov(srt,slope) var(slope)]' of n'th model
%
% List of options for "mode" input:
%         a       proc adds its own noise
%         b*      base figure number for plotting results [100]
%         c       y=proc() returns comment string
%         e/E*    plot evolving psychometric functions *=1,2,3 for mean, mode, marginal mode (F=after each trial)
%         f/F*    plot psychometric functions *=1,2,3 for mean, mode, marginal mode (F=after each trial)
%         g       prompt with number of digits
%         G       prompt with SNR
%         l*      min token length (in digits)
%         L*      max token length (if different from min)
%         m*      use * external models [default 1]
%         M       include an extra model with no processing
%         n*      *=average number of trials per model
%         p/P     plot pdf (P=after each trial)
%         r       allow repetitions
%         s       respond s to save the noisy stimulus as a .wav file
%         t/T*    taper noise */10 seconds at start and end
%         v/V*    plot srt/slope convergence (V= after each trial)
%                 *=1,2,3 for mean, mode, marginal mode
%         x*      add */10 seconds of noise to the front of the speech before processing
%         X*      truncate */10 seconds from the start of the processed speech
%         W       data is in microsoft WAV format
%         z       omit tokens containing "oh"
%         Z       omit tokens containing "zero"
%
% Output file format
%   Each line starts 'x ' where x is one of the following
%       %  Comment
%       V  File version type
%       O  mode options
%       P  details about proc
%       C  Comment returned by proc
%       M  measurement

% Future mode options:
%        [ d     score as single digits ]
%        [ i/I   plot SRT improvement (I=after each trial) ]
%        [ j*    scaling: 0=autoscale each token, 1=constant speech,2=const noise, 3=const noisy ]
%        [N*     ignore the first * trials ]
%        [o/O    write to output file, O write result of each probe]
%        [ S     save all stimuli as wav files ]
%        [ u     do not normalize the noise level ]
%        [ y*    type of noise ]
%
% Bugs/Suggestions
% (1) Add sounds to indicate error and end of test
% (2) Routine could return a label to replace "SNR" in necessary
% (3) Add an input option to "discard" this sample
% (4) output file should include mode argument + date + IP address + computer name + r arguments + p/q values
% (5) Quote filenames in output file or else never have anything else on the line
% (6) silence detection doesn't work

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

%     25+4     m*      use * external models [default 1]
%     26    M       include an extra model with no processing
%      1   a       proc adds its own noise
%      27+5   n*      *=average number of trials per model
%     23+2    l*      min token length (in digits)
%      24+3   L*      max token length (if different from min)
%      7   d       score as single digits
%      37   s       speech-shaped noise instead of white
%       19+1  j*      scaling: 0=autoscale each token, 1=constant speech, 2=const noise, 3=const noisy
%      35   r       allow repetitions
%     41    u       unimodal psychometric function
%     32    P       plot pdf
%     31   p       plot srt convergence
%     13    g       prompt with number of digits
%      14   G       prompt with SNR
%      47+6   x*      add */10 seconds of noise to the front of the speech
%      48+7   X*      truncate */10 from the start of the processed speech
% input parameters
persistent tok mtok tigflgp digitsp
if nargin<10
    ofile='';
    if nargin<9
        dfile='F:\home\dmb\data\old_speech\tidigits';
        %         dfile='Y:\Speech\TIDigits\Disc 1\TIDIGITS';
        if nargin<8
            fn=16000;
            if nargin<7
                noise=[];
                if nargin<6
                    xp=[];
                    if nargin<5
                        q=[];
                        if nargin<4
                            p=[];
                            if nargin<3
                                mode='';
                                if nargin<2
                                    r=[];
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

% parse the mode string

i=1;
pv=[0 3 3 1 25  0  0 1 1 1 1 100]; % pv default values if option is unspecified
px=[0 3 3 1 25 20 18 1 1 1 1 100]; % pv values if option is given without a value
pvs='jlLmnxXvVfFb';
pf=zeros(1,52);
mval=0;
lmode=length(mode);
while i<=lmode
    if i<lmode  % read a following integer if it exists
        [v,nv,e,ni]=sscanf(mode(i+1:end),'%d',1);
    else
        nv=0;
        ni=1;
    end
    j=find(mode(i)==pvs,1);
    k=1+2*(double(lower(mode(i)))-'a')+(mode(i)<'a');
    if k>=1 && k<=52
        pf(k)=1-pf(k);
        if ~isempty(j)
            if nv==0
                pv(j)=px(j);
            else
                pv(j)=v;
                mval=mval || j==4;  % m has a value specified
            end
        end
    end
    i=i+ni;
end
if isempty(proc)
    pv(4)=0; % not allowed any processed models if not process is specified
end
% derived input parameters

varthr=20;   % variance threshold for "silence"
nxevo=30; % size of evloving pdf

if pf(23)>pf(24) || pv(2)>pv(3)
    pv(3)=pv(2);                % Make max token length reasonable
end
if pf(24)>pf(23)
    pv(2)=1;                    % min =1 if only max is specified
end
zmodel=pv(4)+(pf(26)>0); % total number of models (including null model)
ntrial=zmodel*pv(5);  % total number of trials
if pf(12)>pf(11)    % 'F' specified but not 'f'
    pv(10)=pv(11);  % copy across the average type selector
end
if pf(11)+pf(12)>0 && (pv(10)<1 || pv(10)>3)
    error('Invalid average type for option f/F');
end

% sort out p argument to v_psycest or v_psycestu
unimode=pf(41);
if isempty(p)
    p=0.75;     % default recognition rate at threshold
end
if size(p,1)<3-unimode
    p(2-unimode,1)=0.04;     % miss probability
    p(3-unimode,1)=0;     % dummy entry for guess probability
end
if size(p,2)<(zmodel)
    p=p(:,1+mod(0:zmodel-1,size(p,2)));
end

% first sort out the digit samples

tigflg=[pv(2:3) pf(51:52)];
if isempty(tok) || any(tigflg~=tigflgp) || ~strcmp(dfile,digitsp)
    disp('Scanning TIDIGITS files ...');
    digitsp=dfile;
    tigflgp=tigflg;
    if any(dfile(end)=='/\')
        dfile(end)=[];  % remove a trailing separator
    end
    dirlist{1}='';
    ntok=0;
    mtok=zeros(1,pv(3));
    tok=cell(1,2);
    while ~isempty(dirlist)
        dd=dir([dfile dirlist{1}]);
        for i=1:length(dd)
            name=dd(i).name;
            if name(1)~='.'   % ignore directories starting with '.'
                if dd(i).isdir
                    dirlist{end+1}=[dirlist{1} '\' name];
                elseif length(name)>4 && strcmpi(name(end-3:end),'.wav')
                    digs=name(1:end-4);
                    digz=upper(digs)=='Z';
                    digo=upper(digs)=='O';
                    digs(digo | digz)='0';
                    digs=digs(digs>='0' & digs<='9');
                    ndigs=length(digs);
                    if ndigs>=tigflg(1) && ndigs<=tigflg(2) && ~(tigflg(3) && any(digo)) && ~(tigflg(4) && any(digz))
                        ntok=ntok+1;
                        mtok(ndigs)=mtok(ndigs)+1;
                        tok{ntok,1}=[dirlist{1} '\' name];
                        tok{ntok,2}=digs;
                    end
                end
            end
        end
        dirlist(1)=[];  % remove this directory from the list
    end
end
ntok=size(tok,1);
if pf(46)
    [s,fs]=v_readwav([dfile tok{1,1}]); % get the first speech token to set fs
else
    [s,fs]=v_readsph([dfile tok{1,1}]); % get the first speech token to set fs
end

% calculate guess probability assuming you get the correct number of digits

if any(mode=='d')
    p(3-unimode,:)=0.1;
else
    p(3-unimode,:)=0.1.^(1:pv(3))*mtok'/ntok;
end

% now initialize the models

if unimode
    [xx,ii,m,v]=v_psycestu(-zmodel,p,q,xp); % initialize all models
    [pact,qact]=v_psycestu(0);  % save the actual parameters
else
    [xx,ii,m,v]=v_psycest(-zmodel,p,q,xp); % initialize all models
    [pact,qact]=v_psycest(0);  % save the actual parameters
end
if pv(4) && pf(5)
    x=[];  % set x=[] to force a description output
    ii=1;   % for now, only do model 1
    if pf(1)   % if process adds its own noise
        if isempty(r)  % proc does not require any auxilliary parameters
            if mval
                procdesc=proc(x,fs,xx,ii);  % process the noisy speech
            else
                procdesc=proc(x,fs,xx);  % process the noisy speech
            end
        else
            if mval
                procdesc=proc(x,fs,xx,ii,r);  % process the noisy speech
            else
                procdesc=proc(x,fs,xx,r);  % process the noisy speech
            end
        end
    else   % if process does not add its own noise
        if isempty(r)  % proc does not require any auxilliary parameters
            if mval
                procdesc=proc(x,fs,ii);  % process the noisy speech
            else
                procdesc=proc(x,fs);  % process the noisy speech
            end
        else
            if mval
                procdesc=proc(x,fs,ii,r);  % process the noisy speech
            else
                procdesc=proc(x,fs,r);  % process the noisy speech
            end
        end
    end
else
    procdesc='';
end

% now initialize the output file

if pf(29) || pf(30)   % o/O output info
    nw=fix(datevec(now));
    of=sprintf('psy%4d%02d%02d%02d%02d.txt',nw(1:5)); % filename includes date and time to nearest minute
    ofid=fopen(of,'wt');
    if ~ofid
        error('Cannot write to %s',of);
    end
    fprintf(ofid,'%% %s evaluation on %s\n',mfilename,datestr(now));
    fmfnm=[mfilename('fullpath') '.m'];
    dd=dir(fmfnm);
    fprintf(ofid,'%% %s = %d bytes = %s\n',fmfnm,dd.bytes,dd.date);
    fprintf(ofid,'V %d\n',2); % print file format version number
    fmfnm=which(func2str(proc));
    dd=dir(fmfnm);
    fprintf(ofid,'P %s = %d bytes = %s\n',fmfnm,dd.bytes,dd.date);
    if numel(procdesc)
        fprintf(ofid,'C %s\n',procdesc);
    end
    fprintf(ofid,'O %s\n',mode);
end

% now start testing

disp(['Testing ' procdesc]);
% now do the main loop
mnt=zeros(2+2*unimode,zmodel,3,ntrial+1);
vnt=zeros(3+unimode,zmodel,ntrial+1);
mnt(:,:,:,1)=m;
vnt(:,:,1)=v;
i=0;
imax=0;
quitit=0;
while ~quitit
    i=i+1;
    isp=min(1+floor(rand(1)*ntok),ntok); % select a token
    if pf(46)
        [s,fs]=v_readwav([dfile tok{isp,1}]); % get the speech token
    else
        [s,fs]=v_readsph([dfile tok{isp,1}]); % get the speech token
    end
    s=[zeros(pv(6)*round(fs/10),1); v_activlev(s(:),fs,'n')]; % preappend zeros and normalize speech level
    if pf(1) && ii<=pv(4)  % if process adds its own noise
        if isempty(r)
            if mval
                y=proc(s,fs,xx,ii);  % process the noisy speech
            else
                y=proc(s,fs,xx);  % process the noisy speech
            end
        else
            if mval
                y=proc(s,fs,xx,ii,r);  % process the noisy speech
            else
                y=proc(s,fs,xx,r);  % process the noisy speech
            end
        end
    else
        nn=randn(length(s),1);
        x=nn+s*10^(xx/20);   % create the data
        if ii<=pv(4)
            if isempty(r)
                if mval
                    y=proc(x,fs,ii);  % process the noisy speech
                else
                    y=proc(x,fs);  % process the noisy speech
                end
            else
                if mval
                    y=proc(x,fs,ii,r);  % process the noisy speech
                else
                    y=proc(x,fs,r);  % process the noisy speech
                end
            end
        else
            y=x;            % unprocessed for last model
        end
    end
    y=y(1+pv(7)*round(fs/10):end);   % remove junk from the start
    if pf(13)
        prg=sprintf(' %d',length(tok{isp,2}));
    else
        prg='';
    end
    if pf(14)
        prG=sprintf('SNR=%d dB, ',round(xx));
    else
        prG='';
    end
    if pf(35)
        prr=sprintf(', ''r'' to repeat');
    else
        prr='';
    end
    prompt=[prG 'enter' prg ' digits (''q'' to quit' prr '): '];
    %     ansr=-(var(y)>varthr);  % meant to detect silences but doesn't work
    ansr=-1;
    say=1;
    while ansr<0
        if say
            tdel=0;
            tic;
            soundsc(y,fs);      % *** probably shouldn't be ...sc
            say=0;
        end
        rv=input(prompt,'s');
        tdel=toc;
        if ~isempty(rv)
            if lower(rv(1))=='q'
                quitit=1;
                ansr=2;
            elseif lower(rv(1))=='r' && pf(35)
                say=1;
            elseif lower(rv(1))=='s' && pf(37)   % save the token
                ofn=input('File name: ','s');
                if numel(ofn)
                    v_writewav(y,fs,ofn);
                end
            elseif all(rv>='0') && all(rv<='9') && ( ~pf(13) || length(rv)==length(tok{isp,2}))
                ansr=strcmp(rv,tok{isp,2});
            end
        end
    end
    quitit=quitit || i==ntrial;   % mark as quiting if we have done all the trials
    jj=ii;  % remember which model has just been updated
    xxold=xx; % and what the SNR was
    if ansr<2  % valid answer: update the pdfs
        if unimode
            [xx,ii,m,v]=v_psycestu(ii,xx,ansr);
        else
            [xx,ii,m,v]=v_psycest(ii,xx,ansr);
        end
        mnt(:,:,:,i+1)=m;
        vnt(:,:,i+1)=v;
        imax=i;
    end
    if pf(30) || quitit && pf(29)       % 'O/o': output
        if ansr>1
            rv=num2str(double(rv(1)));
        end
        % could add in token name and length
        fprintf(ofid,'M %d %d %.3g %s %s %d %.1f',i,ii,xxold,tok{isp,2},rv,ansr,tdel);
        fprintf(ofid,' %.3g',m(:,jj,:),v(:,jj));
        fprintf(ofid,'\n');
    end
    if pf(32)                           % 'P': plot PDF: figures 1:m
        figure(jj);
        if unimode
            v_psycestu(jj);
        else
            v_psycest(jj);
        end
    elseif quitit && pf(31)             % 'p': plot PDF: figures 1:m
        for jj=1:zmodel
            figure(jj);
            if unimode
                v_psycestu(jj);
            else
                v_psycest(jj);
            end
        end
    end
    if pf(12) || quitit && pf(11)       % 'F/f': plot Psychometric function on figure 101
        figure(pv(12)+1);
        if unimode
            qqu.pk=m(1,jj,pv(10));      % peak position
            qqu.w=m(2,1,pv(10));        % peak width
            qqu.ll=m(3,1,pv(10));       % peak slope on low side
            qqu.lh=m(4,1,pv(10));       % peak slope on high side
            qqu.gu=pact(2);             % guess rate
            qqu.la=pact(1);             % lapse rate
            sw=qqu.w*qqu.ll*qqu.lh/(qqu.ll+qqu.lh);   % normalized distance of inflections from peak
            xax=linspace(qqu.pk-(4+sw)/qqu.ll,qqu.pk+(4+sw)/qqu.lh,200);
            bs=(qqu.lh-qqu.ll)/2;
            bu=(qqu.lh+qqu.ll)/2;
            plot(xax,qqu.gu+(1-qqu.gu-qqu.la)*(1+2*exp(-sw)*cosh(bs*(xax-qqu.pk)+bu*abs(xax-qqu.pk))+exp(-2*sw)).^(-1));
        else
            sd=(pact(1,:)-pact(3,:)).*(1-pact(2,:)-pact(1,:))./(m(2,:,pv(10)).*(1-pact(3,:)-pact(2,:)));
            md=m(1,:,pv(10))-sd.*log((pact(1,:)-pact(3,:))./(1-pact(2,:)-pact(1,:)));
            xax=linspace(min(md-3*sd),max(md+3*sd),100);
            for jj=1:zmodel
                plot(xax,v_psychofunc('',[pact(1,jj); m(:,jj,pv(10)); pact(2:3,jj); qact(10)],xax));
                hold on
            end
            hold off
        end
        axis([xax(1) xax(end) 0 1]);
        xlabel('SNR (dB)');
        ylabel('Recognition Probability');
        title(sprintf('Intelligibility: %s',procdesc));
    end
    if pf(10) || quitit && pf(9)        % 'E/e': plot evolving Psychometric function on figure 103
        figure(pv(12)+3);
        psyevo=zeros(ntrial+1,nxevo);   % space for evolving pdf
        if unimode                      % unimodal psychometric function
            qqu.pk=m(1,jj,pv(10));      % peak position
            qqu.w=m(2,1,pv(10));        % peak width
            qqu.ll=m(3,1,pv(10));       % peak slope on low side
            qqu.lh=m(4,1,pv(10));       % peak slope on high side
            qqu.gu=pact(2);             % guess rate
            qqu.la=pact(1);             % lapse rate
            sw=qqu.w*qqu.ll*qqu.lh/(qqu.ll+qqu.lh);   % normalized distance of inflections from peak
            xax=linspace(qqu.pk-(4+sw)/qqu.ll,qqu.pk+(4+sw)/qqu.lh,nxevo);
            for iet=1:imax+1
                qqu.pk=mnt(1,1,pv(10),iet);     % peak position
                qqu.w=mnt(2,1,pv(10),iet);    	% peak width
                qqu.ll=mnt(3,1,pv(10),iet);     % peak slope on low side
                qqu.lh=mnt(4,1,pv(10),iet);     % peak slope on high side
                sw=qqu.w*qqu.ll*qqu.lh/(qqu.ll+qqu.lh);   % normalized distance of inflections from peak
                bs=(qqu.lh-qqu.ll)/2;
                bu=(qqu.lh+qqu.ll)/2;
                psyevo(iet,:)=qqu.gu+(1-qqu.gu-qqu.la)*(1+2*exp(-sw)*cosh(bs*(xax-qqu.pk)+bu*abs(xax-qqu.pk))+exp(-2*sw)).^(-1);
            end
            imagesc(xax,0:ntrial,psyevo);
            axis 'ij'  % put trial 0 at the top
        else                            % monotonic psychometric function (not implemented)
        end
        xlabel('SNR (dB)');
        ylabel('After trial');
        title(sprintf('Intelligibility: %s',procdesc));
    end
    if pf(44) || quitit && pf(43)   % 'V/v': plot convergence on figure 102
        figure(pv(12)+2);
        if unimode                      % unimodal psychometric function
            sw=mnt(2,1,pv(10),1:imax+1).*mnt(3,1,pv(10),1:imax+1).*mnt(4,1,pv(10),1:imax+1)./(mnt(3,1,pv(10),1:imax+1).*mnt(4,1,pv(10),1:imax+1));
            subplot(221)
            plot(0:imax,squeeze(mnt(1,1,pv(10),1:imax+1)));
            set(gca,'xlim',[0 ntrial]);
            xlabel('After trial');
            ylabel('Peak position (dB SNR)');
            subplot(222)
            plot(0:imax,squeeze(mnt(1,1,pv(10),1:imax+1)-sw.*mnt(3,1,pv(10),1:imax+1)),'-b',0:imax,squeeze(mnt(1,1,pv(10),1:imax+1)+sw.*mnt(4,1,pv(10),1:imax+1)),'-b');
            set(gca,'xlim',[0 ntrial]);
            xlabel('After trial');
            ylabel('Inflections (dB)');
            subplot(223)
            plot(0:imax,squeeze(mnt(3,1,pv(10),1:imax+1)));
            set(gca,'xlim',[0 ntrial]);
            xlabel('After trial');
            ylabel('Upwards slope (prob/dB)');
            subplot(224)
            plot(0:imax,squeeze(mnt(4,1,pv(10),1:imax+1)));
            set(gca,'xlim',[0 ntrial]);
            xlabel('After trial');
            ylabel('Downwards slope (prob/dB)');
        else                            % monotonic psychometric function
            subplot(211);
            for jj=1:zmodel
                plot(0:imax,squeeze(mnt(1,jj,pv(10),1:imax+1)));
                hold on
            end
            hold off
            set(gca,'xlim',[0 ntrial]);
            xlabel('After trial');
            ylabel('SRT (dB SNR)');
            subplot(212);
            for jj=1:zmodel
                plot(0:i,squeeze(mnt(2,jj,pv(10),1:imax+1)));
                hold on
            end
            hold off
            set(gca,'xlim',[0 ntrial]);
            xlabel('After trial');
            ylabel('Slope (Prob/dB)');
        end
    end
end   % main loop for each probe value
if pf(29) || pf(30)
    fclose(ofid);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
% Test file for singfun constructor.

function pass = test_singfun_constructor(pref)

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

%%
% Select some random points as sample points
% These random points are in [-1, 1]
seedRNG(7890)
x = -1 + 2*rand(1, 100);
x = sort(x);

% Some arbitrary values to use for exponents.
a = 0.338745372057174;
b = 0.561224728136042;
a_int = 2;
b_int = 5;

%% Test calling syntax when the user provides exponents

% Negative fractional exponents
fh = @(x) sin(x)./((1+x).^a.*(1-x).^b);
data = struct();
data.exponents = [-a, -b];
f = singfun(fh, data, pref);
data.exponents = [-a, -b];
data.singType = {'sing', 'sing'};
g = singfun(fh, data, pref);
pass(1) = isequal(f,g);
pass(2) = ~any(f.exponents + [a,b]);
pass(3) = ~any(g.exponents + [a,b]);
pass(4) = norm(feval(fh,x) - feval(f,x), inf) < 1e2*eps;

%%
% Positive fractional exponents
fh = @(x) sin(x).*(1+x).^a.*(1-x).^b;
data = struct();
data.exponents = [a, b];
f = singfun(fh, data, pref);
data.exponents = [a, b];
data.singType = {'root', 'root'};
g = singfun(fh, data, pref);
pass(5) = isequal(f,g);
pass(6) = ~any(f.exponents - [a,b]);
pass(7) = ~any(g.exponents - [a,b]);
pass(8) = norm(feval(fh,x) - feval(f,x), inf) < 1e1*eps;

%%
% Negative integer exponents
fh = @(x) exp(x)./((1+x).^a_int.*(1-x).^b_int);
data = struct();
data.exponents = [-a_int, -b_int];
f = singfun(fh, data, pref);
data.exponents = [-a_int, -b_int];
data.singType = {'pole', 'pole'};
g = singfun(fh, data, pref);
pass(9) = isequal(f,g);
pass(10) = ~any(f.exponents + [a_int, b_int]);
pass(11) = ~any(g.exponents + [a_int, b_int]);
% don't check near end-points
xx = x(20:80);
pass(12) = norm(feval(fh,xx) - feval(f,xx), inf) < 1e2*eps;

%% Test Syntax and construction when the user doesn't provide exponents
%
% Negative fractional exponents
fh = @(x) exp(sin(x))./((1+x).^a.*(1-x).^b);
f = singfun(fh);
pass(13) = norm(f.exponents + [a,b], inf) < pref.blowupPrefs.exponentTol;
pass(14) = norm(feval(fh,x) - feval(f,x), inf) < 1e5*eps;

%%
% Positive fractional exponents
fh = @(x) sin(exp(cos(x))).*(1+x).^a.*(1-x).^b;
f = singfun(fh);
pass(15) = norm(f.exponents - [a,b], inf) < pref.blowupPrefs.exponentTol;
pass(16) = norm(feval(fh,x) - feval(f,x), inf) < 1e4*eps;

%%
% Negative integer exponents
fh = @(x) exp(sin(x.^2))./((1+x).^a_int.*(1-x).^b_int);
f = singfun(fh);
pass(17) = norm(f.exponents + [a_int, b_int], inf) < pref.blowupPrefs.exponentTol;
xx = x(20:80);
pass(18) = norm(feval(fh,xx) - feval(f,xx), inf) < 1e3*eps;

%%
% Construction with smoothfuns:
f = smoothfun.constructor( @(x) sin(x));
s = singfun(f);
pass(19) = iszero(f - s.smoothPart);
data = struct();
data.exponents = [-1.5, -1];
data.singType = {'sing', 'sing'};
s = singfun(f, data, pref);
pass(20) = iszero(f - s.smoothPart);
pass(21) = norm(s.exponents - [-1.5, -1], inf) < pref.blowupPrefs.exponentTol;

%%
% Construction from double:
f = singfun(42);
pass(22) = iszero(f - 42);
data = struct();
data.exponents = [1.5, 1];
data.singType = {'sing', 'sing'};
f = singfun(42, data, pref);
g = singfun(@(x) 42+0*x);
g.exponents = [1.5, 1];
pass(23) = iszero(f - g);
pass(24) = norm(s.exponents - [-1.5, -1], inf) < pref.blowupPrefs.exponentTol;

end


text:
classdef linopConstraint
%LINOPCONSTRAINT Class for linop constraints.
%
%   A linop operates on a set of chebfuns and scalars (the "variables" of the
%   linop). An instance of this class may be assigned to the 'constraint'
%   property of a linop in order to impose a constraint on that variable (for
%   purposes of solving a linear system or eigenvalue problem).
%
%   Each LINOPCONSTRAINT object has a 'functional' property, which is a
%   chebmatrix of vertically concatenated functionalBlocks, and a 'values'
%   property, a vector of the same length. If u is the variable of the
%   linop, then the constraint imposed by object C is:
%
%        C.functional * u = C.values.
%
%   The preferred way to add contraints to a linop is through LINOP.ADDBC,
%   not through direct use of this class.
%
% See also LINOP, LINOP.ADDBC.

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Access = public )
        functional    % applied to the variable to get values
        values        % constraint on the result of the functional
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        % Input a functional and a value to create a constraint, or create
        % an empty constraint if no inputs.
        function C = linopConstraint(op, vals)
            if ( nargin == 0 )
                return
            end
            C.functional = op;
            C.values = vals;
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function n = length(C)
        %LENGTH    Number of constraints in the object.
            n = size(C.functional, 1);
        end

        function e = isempty(C)
        %ISEMPTY   True if no constraints in the object.
            e = isempty(C.functional);
        end

        function C = append(C, func, value)
        %APPEND    Insert an additional constraint.
        %   C = APPEND(C, FUNC, VAL) appends the constraint FUNC*u=VAL to the
        %   current list. If VAL is omitted, it defaults to zero.

            % Appending empty does nothing:
            if ( isempty(func) )
                return
            end

            % Check if func is of an allowed class.
            validateattributes(func, {'linBlock', 'chebmatrix'}, {})

            % Check the value.
            if ( nargin < 3 )
                value = 0;
            end
            validateattributes(value, {'double'}, {})

            C.functional = [ C.functional ; func ];
            C.values = [ C.values ; value ];
        end

        function C = uminus(C)
        %-   Negate the VALUES property of a LINOPCONSTRAINT.
        %   This is useful at the CHEBOP level, where we need different signs
        %   for the boundary conditions of a LINOP when we call LINOP(N) where N
        %   is a CHEBOP, compared to what we want when we call LINEARIZE() from
        %   within a Newton iteration. This is because when problems are solved
        %   with LINOP backslash, the solution to the problem is the output
        %   itself, while in a Newton iteration, we have to add the output of
        %   the LINOP solution to the current guess.
            C.values = -C.values;
        end

    end

end


text:
% PROBABILITY DISTRIBUTION FUNCTIONS (in the dist-folder):
%
% Priors
%  PRIOR_FIXED       Fix parameter to its current value
%  PRIOR_GAMMA       Gamma prior structure
%  PRIOR_GAUSSIAN    Gaussian prior structure
%  PRIOR_INVGAMMA    Inverse-gamma prior structure
%  PRIOR_INVT        Inverse Student-t prior structure
%  PRIOR_INVT        Inverse uniform prior structure
%  PRIOR_LAPLACE     Laplace (double exponential) prior structure
%  PRIOR_LOGGAUSSIAN Log-Gaussian prior structure
%  PRIOR_LOGLOGUNIF  Uniform prior structure for the log(log(parameter))
%  PRIOR_LOGT        Student-t prior structure for the logarithm of the parameter
%  PRIOR_LOGUNIF     Uniform prior structure for the logarithm of the parameter
%  PRIOR_SINVCHI2    Scaled inverse-chi-square prior structure
%  PRIOR_SQRTT       Student-t prior structure for the square root of the
%                    parameter
%  PRIOR_SQRTUNIF    Uniform prior structure for the square root of the
%                    parameter
%  PRIOR_UNIF        Uniform prior structure
%  PRIOR_T           Student-t prior structure
%  PRIOR_SQRTUNIF    Uniform prior structure for the square root of the
%                    parameter
%
% Probability density functions
%
%    BETA_LPDF     - Beta log-probability density function (lpdf).
%    BETA_PDF      - Beta probability density function (pdf).
%    DIR_LPDF      - Log probability density function of uniform Dirichlet
%                    distribution
%    DIR_PDF       - Probability density function of uniform Dirichlet
%                    distribution
%    GAM_CDF       - Cumulative of Gamma probability density function (cdf).
%    GAM_LPDF      - Log of Gamma probability density function (lpdf).
%    GAM_PDF       - Gamma probability density function (pdf).
%    GEO_LPDF      - Geometric log probability density function (lpdf).
%    INVGAM_LPDF   - Inverse-Gamma log probability density function.
%    INVGAM_PDF    - Inverse-Gamma probability density function.
%    LAPLACE_LPDF  - Laplace log-probability density function (lpdf).
%    LAPLACE_PDF   - Laplace probability density function (pdf).
%    LOGN_LPDF     - Log normal log-probability density function (lpdf)
%    LOGT_LPDF     - Log probability density function (lpdf) for log Student's T
%    MNORM_LPDF    - Multivariate-Normal log-probability density function (lpdf).
%    MNORM_PDF     - Multivariate-Normal log-probability density function (lpdf).
%    NORM_LPDF     - Normal log-probability density function (lpdf).
%    NORM_PDF      - Normal probability density function (pdf).
%    POISS_LPDF    - Poisson log-probability density function.
%    POISS_PDF     - Poisson probability density function.
%    SINVCHI2_LPDF - Scaled inverse-chi log-probability density function.
%    SINVCHI2_PDF  - Scaled inverse-chi probability density function.
%    T_LPDF        - Student's T log-probability density function (lpdf)
%    T_PDF         - Student's T probability density function (pdf)
%
% Random number generators
%
%    CATRAND       - Random matrices from categorical distribution.
%    DIRRAND       - Uniform dirichlet random vectors
%    EXPRAND       - Random matrices from exponential distribution.
%    GAMRAND       - Random matrices from gamma distribution.
%    INTRAND       - Random matrices from uniform integer distribution.
%    INVGAMRAND    - Random matrices from inverse gamma distribution
%    INVGAMRAND1   - Random matrices from inverse gamma distribution
%    INVWISHRND    - Random matrices from inverse Wishart distribution.
%    NORMLTRAND    - Random draws from a left-truncated normal
%                    distribution, with mean = mu, variance = sigma2
%    NORMRTRAND    - Random draws from a right-truncated normal
%                    distribution, with mean = mu, variance = sigma2
%    NORMTRAND     - Random draws from a normal truncated to interval
%    NORMTZRAND    - Random draws from a normal distribution truncated by zero
%    SINVCHI2RAND  - Random matrices from scaled inverse-chi distribution
%    TRAND         - Random numbers from Student's t-distribution
%    UNIFRAND      - Generate unifrom random numberm from interval [A,B]
%    WISHRND       - Random matrices from Wishart distribution.
%
% Others
%    KERNELP       - Kernel density estimator for one dimensional distribution.
%    HAMMERSLEY    - Hammersley quasi-random sequence
%

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function final_result = check_yale_faces(num_subjects, solver, solver_name, solver_params)
if nargin < 4
    solver_params = struct;
end

% Create the directory for storing results
spx.fs.ensure_dir('bin');
yf = spx.data.image.ChongYaleFaces();
num_total_subjects = yf.num_subjects;

% R = yf.num_experiments_of_n_subjects(num_subjects);
if num_subjects == 38
    R = 1;
else
    R = 20;
end
fprintf('Number of subject sets: %d\n', R);

% We will carry out one experiment for each signal density level
test_result.failed_examples = [];
trials = cell(1, R);

for r=1:R
    trial = struct();
    rng(r * 38 + num_subjects);
    subject_list = randperm(num_total_subjects, num_subjects) % select #nCluster subjects
    fprintf('Subject set (%d): ', r);
    fprintf('%d ', subject_list);
    fprintf('\n');
    % [Y, cluster_sizes, true_labels] = yf.experiment_data(num_subjects, r);
    [Y, cluster_sizes, true_labels] = yf.data_for_subjects(subject_list);
    % Normalize the columns
    Y = spx.norm.normalize_l2(Y);
    % Ambient space dimension and number of data points
    [trial.M, trial.S] = size(Y);
    % Number of subspaces
    trial.K = num_subjects;
    % maximum dimension for each subspace
    trial.D = 5;
    trial.cluster_sizes = cluster_sizes;
    % Solve the sparse subspace clustering problem
    tstart = tic;
    try
        clustering_result = solver(Y, trial.D, trial.K, solver_params);
    catch ME
        % we will move on to next one
        fprintf('Problem in processing this example. %s: %s\n', ME.identifier, ME.message);
        % add to the list of failed examples
        test_result.failed_examples = [test_result.failed_examples r];
        % move on to next example
        throw(ME)
        continue;
    end
    trial.elapsed_time = toc (tstart);
    % graph connectivity
    trial.connectivity = clustering_result.connectivity;
    % estimated number of clusters
    trial.estimated_num_subspaces = clustering_result.num_clusters;
    % Time to compare the clustering
    cluster_labels = clustering_result.labels;
    comparsion_result = spx.cluster.clustering_error_hungarian_mapping(cluster_labels, true_labels, trial.K);
    trial.clustering_error_perc = comparsion_result.error_perc;
    trial.clustering_acc_perc = 100 - comparsion_result.error_perc;
    % Compute the statistics related to subspace preservation
    spr_stats = spx.cluster.subspace.subspace_preservation_stats(clustering_result.Z, trial.cluster_sizes);
    trial.spr_error = spr_stats.spr_error;
    trial.spr_flag = spr_stats.spr_flag;
    trial.spr_perc = spr_stats.spr_perc;
    fprintf('\nclustering error: %0.2f %% , clustering accuracy: %0.2f %%\n, mean spr error: %0.4f preserving : %0.2f %%\n, connectivity: %0.2f, elapsed time: %0.2f sec', trial.clustering_error_perc, trial.clustering_acc_perc, spr_stats.spr_error, spr_stats.spr_perc, trial.connectivity, trial.elapsed_time);
    fprintf('\n\n');
    trials{r} = trial;
end
% clear some variables no more required.
clear trial;
clear yf;
clear clustering_result;
clear comparison_result;
%  save rest of variables in file.
filepath = sprintf('bin/yale_faces_%d_subjects_test_%s.mat', num_subjects, solver_name);
save(filepath);
final_result = merge_results(num_subjects, solver_name);
end

author: indigits

repo: sparse-plex