text:
function intNum = whichInterval(disc, location, direction)
%WHICHINTERVAL Index of the subinterval of a given point and direction.
%   INTNUM = WHICHINTERVAL(DISC, LOC, DIRN) returns which interval within
%   DISC.domain a given point LOC is located in. The DIRN is +1 or -1 and
%   indicates approach from the right or left (zero means don't care).

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

% [TODO]: How does this differ from chebfun.whichInterval?

% Quick shortcut if there are no breakpoints.
numIntervals = length(disc.domain);
if (  numIntervals == 2 )
    intNum = 1;
    return
end

% If there were no direction, this would be it.
intNum = find( location >= disc.domain, 1, 'last' );

% LINOP already screened to make sure the location is in the interval, so we can
% check for being at the right endpoint.
if ( intNum == numIntervals )
    if ( direction > 0 )
        error('CHEBFUN:OPDISCRETIZATION:whichInterval:undefined', ...
            'Evaluation direction is undefined at the location.')
    end
    direction = -1;  % this forces the adjustment below
end

% Need to decrement if at a breakpoint coming from the left, or if at
% the right endpoint.
len = disc.domain(end) - disc.domain(1);
if ( direction < 0 && abs( location - disc.domain(intNum) ) < 10*eps*len )
    intNum = intNum - 1;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [Ex Ey Ez] = MiePECBall(cart, k, R, N)
% The code caculates the time-harmonic scattered field value at point 'cart'
% due to a PEC ball with radius R centered at the origin.
% The incident wave is E=exp(-ikz).
% The formulation is based on H.C. van de Hulst, P123
% 'light scattering by small particles', Dover, 1981
% Input: cart -- cartiesian coordiante (1x3) of the point ouside the ball
%        k    -- wavenumber
%        R    -- radius of the PEC ball
%        N    -- number of terms used in Mie series
% Output: [Ex Ey Ez] -- the scattering electric field
% Author: Jiguang Sun, 02/13/2011, jsun@desu.edu
% Please report bugs to jsun@desu.edu
% Copyright (c) 2011, Jiguang Sun, rights reserved.
% THIS SOFTWARE IS PROVIDED "AS IS".
% Redistribution and use in source and binary forms, with or without
% modification, are permitted.
i = sqrt(-1);
[phi,th,r] = cart2sph(cart(1), cart(2), cart(3)); th = pi/2-th;
if r < R
    % disp('The point is inside the ball!')
    Ex = 0; Ey = 0; Ez = 0;
    return;
end

n = 1:N+2;
x = k*R;
[j(n), ierr(1,n)] = besselj(n - 1 + 1/2, x); j = j*sqrt(pi/(2*x));
[h(n), ierr(1,n)] = besselh(n - 1 + 1/2, 2, x); h = h*sqrt(pi/(2*x));
if any(any(ierr))
    disp('Accuracy error in a Bessel or Hankel function.');
    Ex = 0; Ey = 0; Ez = 0;
    return;
end
% scattering coefficients
an = j(2:N+1)./h(2:N+1);
bn = (j(2:N+1)+k*R*j(1:N)-k*R*j(3:N+2))./(h(2:N+1)+k*R*h(1:N)-k*R*h(3:N+2));
% various terms ralated to Racati-Bessel's functions
n = 1:N+4;
z = k*r; k2r = k^2*r;
[h(n), ierr(2,n)] = besselh(n - 2 + 1/2, 2, z); h = h*sqrt(pi/(2*z));
hn = h(3:N+2);
drh = (h(3:N+2)+k*r*h(2:N+1)-k*r*h(4:N+3))/2;
ddrh=(k2r*h(1:N)+2*k*h(2:N+1)-(1/r+2*k2r)*h(3:N+2)-2*k*h(4:N+3)+k2r*h(5:N+4))/4;
%
if any(any(ierr))
    disp('Accuracy error in a Bessel or Hankel function.')
    Ex = 0; Ey = 0; Ez = 0;
    return;
end
% various terms related to associated Legendre polynomials
Pn(1) = 1;
Pn(2) = cos(th);
for n=2:N-1
    Pn(n+1) = ((2*n-1)*cos(th)*Pn(n)-(n-1)*Pn(n-1))/n;
end
% Pn1 starts from n = 1.
Pn1(1) = -sin(th);
Pn1(2) = -3.0*sin(th)*cos(th);
for n=2:N-1
    Pn1(n+1) = (2.0*n+1)*cos(th)*Pn1(n)/n-(n+1)*Pn1(n-1)/n;
end
% Pn1S starts from n = 1;
Pn1S(1) = -1;
Pn1S(2) = -3.0*cos(th);
for n = 2:N-1
    Pn1S(n+1) = Pn1S(n-1)-(2*n+1)*Pn(n+1);
end

for n = 1:N
    if n == 1
        dPn1(n) = -cos(th);
    elseif n == 2
        dPn1(n) = 3*(2*sin(th)*sin(th)-1);
    else
        dPn1(n) = (2*n-1)/(n-1)*Pn1(n-1)*(-sin(th))+(2*n-1)/(n-1)*cos(th)*dPn1(n-1)-n/(n-1)*dPn1(n-2);
    end
end
% compute M and N
Mrho = 0; Mth = 0; Mphi = 0;
for n = 1:N
    c = (-i)^(n)*(2*n+1)/(n*(n+1));
    Mth = Mth + c*an(n)*hn(n)*Pn1S(n)*cos(phi);
    Mphi = Mphi-c*an(n)*hn(n)*dPn1(n)*sin(phi);
end
Nrho = 0; Nth = 0; Nphi = 0;
for n = 1:N
    c = (-i)^(n)*(2*n+1)/(n*(n+1));
    Nrho = Nrho + bn(n)*c*(ddrh(n)+k^2*r*hn(n))*Pn1(n)*cos(phi)/(k);
    Nth = Nth + bn(n)*c/r*drh(n)*dPn1(n)*cos(phi)/(k);
    Nphi = Nphi + bn(n)*c/r*drh(n)*Pn1S(n)*(-sin(phi))/(k);
end
% compute scattering field in Spherical coordinate
Erho = Mrho + i*Nrho;
Ephi = Mphi + i*Nphi;
Eth = Mth + i*Nth;
% change back to cartisian coordinate
Ex = Erho*sin(th)*cos(phi)+Eth*cos(th)*cos(phi)-Ephi*sin(phi);
Ey = Erho*sin(th)*sin(phi)+Eth*cos(th)*sin(phi)+Ephi*cos(phi);
Ez = Erho*cos(th)-Eth*sin(th);

author: Sable

repo: mcbench-benchmarks


text:
function toolbarShortcuts = GetToolbarShortcuts(varargin)
%GETTOOLBARSHORTCUTS is deprecated.  Use GetShortcuts instead.

% $Author: rcotton $	$Date: 2010/08/20 10:21:24 $	$Revision: 1.4 $
% Copyright: Health and Safety Laboratory 2010

toolbarShortcuts = GetShortcuts(varargin{:});

end


text:
function [mantissa,exponent,signBit]=decomposeFloat(x,decodeExp,addLeading1)
%%DECOMPOSEFLOAT Given a floating point number, extract the sign but, the
%          mantissa and the exponent as three separate integer quantities.
%          Unlike the function [F,E]=log2(x), this provies the raw bits of
%          the mantissa rather than another float as in the F returned by
%          log2. The floating point value is assumed to be in IEEE 754
%          format, which is documented in [1].
%
%INPUTS: x The real, scalar floating point value to decompose. This can be
%          of class 'single' or 'double'.
% decodeExp The raw exponent of the floating point number is biased by 1023
%          for a double and by 127 for a single. The bits are shifted (from
%          the original double) so that the first bit in the return value
%          exponent is the least significant bit in the exponent. If this
%          is true, then the bias is removed. The default if this parameter
%          is omitted or an empty matrix is passed is true. Note that when
%          decoded, the exponent for 1 is 0, but for 0 it is actually -1023
%          for a double and is -127 for a single.
% addLeading1 For mantissas of normalized numbers, if this is true (the
%          default), then the leading 1 that is implied but not expressly
%          stored in the floating point register will be added. No 1 will
%          be added for subnormal numbers, where it is not implied that
%          there is a leading 1. Note that "leading" means that it is added
%          past the last bit in the mantissa.
%
%OUTPUTS: mantissa The mantissa of the float as a uint64 if x is a double
%                  or uint32 if x is a single. If addLeading1=true, then
%                  there are 53 significant bits for double and 24
%                  significant bits for a single. Otherwise, there are
%                  respectively 52 and 23 significant bits for a double or
%                  for a single.
%         exponent The exponent value. If decodeExp is false, then this is
%                  a uint64 if x is a double or a uint32 if x is a single.
%                  If decodeExp is true, then this is respectively an int64
%                  or an int32 for a double or for a single.
%          signBit This is 1 if the sign bit in x is 1 (negative) and 0 if
%                  the sign bit in x is zero. This is a uint64 if x is a
%                  double or a uint32 if x is a single.
%
%The IEEE 754 standard for floating point arithmetic in [1] formats a
%double as a 64-bit floating point data type as
%[mantissa (52 bits)] [exponent (11 bits)] [sign (1 bit)]
%where bit numbering starts at 0 in the mantissa and goes through 63 with
%the sign bit. For a 32-bit floating point single, the bit ordering is
%[mantissa (23 bits)] [exponent (8 bits)] [sign (1 bit)]
%If exponent is all ones and the mantissa is zero, then the value Inf is
%represented with the sign given by the sign bit. If the exponent is all
%ones and the mantissa is any nonzero value, then a NaN is given.
%
%To interpret the mantissa, consider that there is an implied binary point
%after the highest bit and if exponent is nonzero, then an implied 1. If
%exponent is zero, then it is implied that a 0 leads the binary point.
%Thus, the mantissa represents a fraction >=0 and < 0.5. If one reverses
%the order of the bits in the mantissa (because we wrtie numbers in big-
%endian format), then the number represented by the float is of the form
%1.mantissa*2^(exponent-bias) for a nonzero exponent and
%0.mantissa*2^(-bias) for a zero exponent. The bias is 1023 for a double
%and is 127 for a single.
%
%EXAMPLE:
% [mantissa,exponent,signBit]=decomposeFloat(4568.321)
%One gets mantissa=5022922058913284, exponent=12, signBit=0
%
%EXAMPLE:
%When dealing with integers, the mantissa with the leading 1 added
%multiplied by a power of two (shifted) is the integer value. Consider
% [mantissa,exponent,signBit]=decomposeFloat(11)
% %One gets mantissa=6192449487634432, exponent=3, signBit=0
% %Now, shift until the first nonzero bit is at the start.
% mantissaShifted=bitshift(mantissa,-findPosOfMin1Bit(mantissa)+1)
%One will see that the shifted mantissa is the orignal value, 11. The
%actual amount one has to shift the mantissa also depends on the exponent.
%
%REFERENCES:
%[1] IEEE Standard for Floating Point Arithmetic, Institute for Electrical
%    and Electronics Engineers Std. IEEE Std 754-2008, 29 Aug. 2008.
%
%January 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(decodeExp))
    decodeExp=true;
end

if(nargin<3||isempty(addLeading1))
    addLeading1=true;
end

if(~isreal(x)||~isscalar(x))
   error('x must be a real, scalar, floating point value.')
end

if(isa(x,'double'))
    N=typecast(x,'uint64');

    %Get the sign bit.
    theMask=getBitMask(63,'uint64');
    signBit=bitshift(bitand(N,theMask),-63);

    %Select the bits that make up the exponent.
    expMask=getBitMask([52,62],'uint64');
    exponent=bitshift(bitand(N,expMask),-52);

    %Select the bits that make up the mantissa.
    theMask=getBitMask([0,51],'uint64');
    mantissa=bitand(N,theMask);
    if(addLeading1&&exponent~=0&&exponent~=expMask)
        %If a leading 1 is implied by the floating point number, then add
        %it. No leading 1 is added if theExp is the value indicating an INf
        %or a NaN.

        mantissa=bitor(mantissa,bitshift(uint64(1),52));
    end

    if(decodeExp)
        %Get rid of the bias. Note that two exceptions are that +/-0 maps
        %to -1023 instead of 0 and +/-Inf and NaN map to 1024.
        exponent=int64(exponent)-1023;
    end
elseif(isa(x,'single'))
    N=typecast(x,'uint32');

    %Get the sign bit.
    theMask=getBitMask(31,'uint32');
    signBit=bitshift(bitand(N,theMask),-31);

    %Select the bits that make up the exponent.
    expMask=getBitMask([23,30],'uint32');
    exponent=bitshift(bitand(N,expMask),-23);

    %Select the bits that make up the mantissa.
    theMask=getBitMask([0,22],'uint32');
    mantissa=bitand(N,theMask);
    if(addLeading1&&exponent~=0&&exponent~=expMask)
        %If a leading 1 is implied by the floating point number, then add
        %it. No leading 1 is added if theExp is the value indicating an INf
        %or a NaN.

        mantissa=bitor(mantissa,bitshift(uint32(1),23));
    end

    if(decodeExp)
        %Get rid of the bias. Note that two exceptions are that +/-0 maps
        %to -127 instead of 0 and +/-Inf and NaN map to 127.
        exponent=int32(exponent)-127;
    end
else
    error('x must be a single or double floating point value.')
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% IN ORDER TO SIMULATE THE PROGRAM:
%   A) FIRST, LOAD A ROBOT
%       robot = load_robot('abb','irb140');
%       OR
%       robot = load_robot('abb','irb52');
%       or any other robot, the target points may not be reachable,
%       depending on the links' lengths.
%
%   B) NEXT, LOAD SOME EQUIPMENT.
%       robot.equipment = load_robot('equipment','tables/table_two_areas');
%
%
%
% The example presents the robot achieving different positions and orientations
% In addition, the conf values, change accordingly.

function test_configurations2

global  TD_tool0 robot

%definition of the end effector.
TD_tool0=[1,[[0,0,0],[1,0,0,0]],[0,[0,0,0],[1,0,0,0],0,0,0]];

main
end

function main
global  TD_tool0 robot

MoveJ([[217.4,1.57,1072.81]/1000,[0.972982,-0.01947,0.224767,-0.049073],[0,1,-3,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]],'vmax','z50',TD_tool0);
MoveJ( [[86.96,1.57,1072.8]/1000,[0.972992,-0.01944,0.224739,-0.049009],[0,0,-1,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]],'vmax','z50',TD_tool0);
MoveJ( [[28.15,1.55,1072.8]/1000,[0.973021,-0.019611,0.224611,-0.048949],[1,-2,0,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]],'vmax','z50',TD_tool0);
MoveJ( [[14.42,1.55,1072.8]/1000,[0.973015,-0.019655,0.224611,-0.049041],[1,-2,-1,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]],'vmax','z50',TD_tool0);
MoveJ( [[0.11,1.57,1072.82]/1000,[0.973031,-0.019796,0.224522,-0.049091],[1,-2,-1,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]],'vmax','z50',TD_tool0);
MoveJ( [[-12.28,1.55,1072.81]/1000,[0.973022,-0.019733,0.224521,-0.049298],[1,-2,-1,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]],'vmax','z50',TD_tool0);
MoveJ( [[-98.02,1.54,1072.82]/1000,[0.973023,-0.019651,0.22448,-0.049486],[1,-2,-1,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]],'vmax','z50',TD_tool0);
MoveJ( [[-231.93,1.54,1072.81]/1000,[0.973047,-0.019566,0.224335,-0.049716],[1,-2,-1,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]],'vmax','z50',TD_tool0);

end

repo: ARTE


text:
function fields = splitString(string, delimiter)
% Splits a string Perl style
%
% USAGE:
%
%     fields = splitString(string, delimiter)
%
% INPUTS:
%    string:       Either a single string or a cell array of strings
%    delimiter:    Splitting delimiter
%
% OUTPUT:
%    fields:       Either a single cell array of fields or a cell array of cell
%                  arrays of fields
%
% Default delimiter is '\s' (whitespace)
% Delimiters are perl regular expression style, e.g. '|' has to be expressed
% as '\|'
% Results are returned in the cell array fields
%
% .. Authors:  Markus Herrgard 07/14/04

if nargin < 2
    delimiter = '\s';
end

% Check if this is a list of strings or just a single string
if iscell(string)
    stringList = string;
    for i = 1:length(stringList)
        fields{i} = splitOneString(stringList{i}, delimiter);
    end
else
    fields = splitOneString(string, delimiter);
end

fields = columnVector(fields);

function fields = splitOneString(string, delimiter)
% Internal function that splits one string

[startIndex, endIndex] = regexp(string, delimiter);

if ~isempty(startIndex)
    cnt = 0;
    for i = 1:length(startIndex) + 1
        if i == 1
            if endIndex(i) > 1
                cnt = cnt + 1;
                fields{cnt} = string(1:endIndex(i) - 1);
            end
        elseif i == length(startIndex) + 1
            if startIndex(i - 1) < length(string)
                cnt = cnt + 1;
                fields{cnt} = string(startIndex(i - 1) + 1:end);
            end
        else
            cnt = cnt + 1;
            fields{cnt} = string(startIndex(i - 1) + 1:endIndex(i) - 1);
        end
    end
else
    fields{1} = string;
end

fieldsOut = {};
cnt = 0;
for i = 1:length(fields)
    if ~isempty(fields{i})
        cnt = cnt + 1;
        fieldsOut{cnt} = fields{i};
    end
end
fields = fieldsOut;

author: opencobra

repo: cobratoolbox


text:
%% A Unit Test Class for the OPTI Object
classdef opti_tests < matlab.unittest.TestCase

    properties

    end

    % Unit Tests
    methods (Test)

        %-- Single Precision: Nonlinear --%
        function singleNonlinear(testCase)

            optiObj = opti_tests.makeHS71();
            defaultOpts = optiset(optiObj.opts, 'warnings', 'none', 'derivCheck', 'on');

            testCase.verifyError(@() opti(optiObj, 'fun', @(x) single(opti_tests.hs71_obj(x))), 'OPTI:NotDouble');
            testCase.verifyError(@() opti(optiObj, 'nlcon', @(x) single(opti_tests.hs71_con(x))), 'OPTI:NotDouble');
            testCase.verifyError(@() opti(optiObj, 'grad', @(x) single(opti_tests.hs71_grad(x)), 'opts', defaultOpts), 'OPTI:NotDouble');
            testCase.verifyError(@() opti(optiObj, 'nljac', @(x) single(opti_tests.hs71_jac(x)), 'opts', defaultOpts), 'OPTI:NotDouble');
            testCase.verifyError(@() opti(optiObj, 'hess', @(x,s,l) single(opti_tests.hs71_hess(x,s,l)), 'opts', defaultOpts), 'OPTI:NotDouble');
        end
    end

    methods (Static)
        %-- Hock & Schittkowski #71 --%
        function optiObj = makeHS71()
            cl = [25;40];
            cu = [Inf;40];
            lb = ones(4,1);
            ub = 5*ones(4,1);
            x0 = [1 5 5 1]';
            opts = optiset('warnings', 'none');
            optiObj = opti('fun',@opti_tests.hs71_obj,'grad',@opti_tests.hs71_grad,'hess',@opti_tests.hs71_hess,...
                            'nl',@opti_tests.hs71_con,cl,cu,'nljac',@opti_tests.hs71_jac,'nljacstr',@opti_tests.hs71_jacStr,...
                            'hessstr',@opti_tests.hs71_hessStr,'bounds',lb,ub, 'x0', x0, 'opts', opts);
        end

        function obj = hs71_obj(x)
            obj = x(1)*x(4)*sum(x(1:3)) + x(3);
        end

        function grad = hs71_grad(x)
            grad = [ x(1)*x(4) + x(4)*sum(x(1:3))
                    x(1)*x(4)
                    x(1)*x(4) + 1
                    x(1)*sum(x(1:3)) ]';
        end

        function con = hs71_con(x)
            con = [ prod(x); sum(x.^2) ];
        end

        function jac = hs71_jac(x)
            jac = [ prod(x)./x'; 2*x' ];
        end

        function jacStr = hs71_jacStr()
            jacStr = sparse(ones(2,4));
        end

        function hess = hs71_hess(x, sigma, lambda)
            hess = sigma*[ 2*x(4)             0      0   0;
                              x(4)               0      0   0;
                              x(4)               0      0   0;
                              2*x(1)+x(2)+x(3)  x(1)  x(1)  0 ] + ...
                   lambda(1)*[   0          0         0         0;
                              x(3)*x(4)     0         0         0;
                              x(2)*x(4) x(1)*x(4)     0         0;
                              x(2)*x(3) x(1)*x(3) x(1)*x(2)     0  ] + ...
         		   lambda(2)*diag([2 2 2 2]);
        end

        function hessStr = hs71_hessStr()
            hessStr = sparse(tril(ones(4)));
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% OPTIMIZE_FILTER Optimize filter representation
%
% Usage
%    filter = OPTIMIZE_FILTER(filter_f, lowpass, options)
%
% Input
%    filter_f (numeric): The Fourier transform of the filter.
%    lowpass (boolean): If true, filter_f contains a lowpass filter.
%    options (struct): Various options on how to optimize the filter:
%       options.filter_format (char): Specifies the type of optimization,
%          either 'fourier', 'fourier_multires' or 'fourier_truncated'. See
%          description for more details.
%       options.truncate_threshold (numeric): If options.filter_format is
%          'fourier_truncated', this indicates the threshold to be passed on
%          to TRUNCATE_FILTER. See the documentation of this function for more
%          details.
%
% Output
%    filter (struct or numeric): The optimized filter structure.
%
% Description
%    Depending on the value of options.filter_format, OPTIMIZE_FILTER calls
%    different functions to optimize the filter representation. If 'fourier',
%    the function retains the Fourier representation of the filter. If
%    'fourier_multires', the filter is periodized and stored at all resolu-
%    tions using PERIODIZE_FILTER. Finally, if it equals 'fourier_truncated',
%    TRUNCATE_FILTER is called on filter_f.
%
% See also
%    PERIODIZE_FILTER, TRUNCATE_FILTER

function filter = optimize_filter(filter_f, lowpass, options)
	options = fill_struct(options,'truncate_threshold',1e-3);
	options = fill_struct(options,'filter_format','fourier_multires');

	if strcmp(options.filter_format,'fourier')
		filter = filter_f;
	elseif strcmp(options.filter_format,'fourier_multires')
		filter = periodize_filter(filter_f);
	elseif strcmp(options.filter_format,'fourier_truncated')
		filter = truncate_filter(filter_f,options.truncate_threshold,lowpass);
	else
		error(sprintf('Unknown filter format ''%s''',options.filter_format));
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
%% PSWARM Install for OPTI Toolbox
% Copyright (C) 2012 Jonathan Currie (I2C2)

% This file will help you compile PSwarm for use with MATLAB.

% My build platform:
% - Windows 8 x64
% - Visual Studio 2012
% - Intel Math Kernel Library

% To recompile you will need to get / do the following:

% 1) Get PSwarm
% PSwarm is available from http://www.norg.uminho.pt/aivaz/pswarm/.
% Download PPSwarm_vxx.zip (C Version) and unzip to a suitable location.

% 2) Compile PSwarm
% The easiest way to compile PSwarm is to use the Visual Studio Project
% Builder included with OPTI. Use the following commands, substituting the
% required path on your computer:
%
% %% Visual Studio Builder Commands
% path = 'full path to PSwarm here'; %e.g. 'C:\Solvers\PSwarm'
% sdir = path;
% name = 'libpswarm';
% opts.exPP = {'LINEAR','_CRT_SECURE_NO_WARNINGS','_CRT_NONSTDC_NO_DEPRECATE'};
% opts.exclude = {'showcache.c','user.c','pswarm_main.c','pswarm_py.c','pswarm_r.c','cache.c'};
% VS_WriteProj(sdir,name,[],opts)
% %%
%
% Once complete, you will have a directory called PSwarm\libpswarm. Open the
% Visual Studio 2012 project file, then complete the following steps:
%   a) A couple code changes are required in order to interface with OPTI:
%       - In pswarm.h add "int solveriters;" to the Stats structure as the
%       last line.
%       - In pswarm.h change the define on line 62 to
%         "#define SYS_RANDOM 1"
%       - In pswarm.c on line 666 add "stats.solveriters = iter;"
%       - Change the structure definitions in pswarm.c at line 66 to extern:
%           "extern struct Stats stats;"
%           "extern struct Options opt;"
%       - Comment the default options structure in pswarm.c on line 70
%       (comment all lines for this structure).
%   b) Build a Win32 or x64 Release to compile the code.
%   c) Copy the generated .lib file to the following folder:
%
%   OPTI/Solvers/pswarm/Source/lib/win32 or win64
%
%   You will also need to copy pswarm.h to the following folder:
%
%   OPTI/Solvers/pswarm/Source/Include

% 3) PSwarm MEX Interface
% The PSwarm MEX Interface is a simple MEX interface I wrote to use PSwarm.

% 4) Compile the MEX File
% The code below will automatically include all required libraries and
% directories to build the PSwarm MEX file. Once you have completed all the
% above steps, simply run this file to compile PSwarm! You MUST BE in the
% base directory of OPTI!

clear pswarm

% Modify below function if it cannot find Intel MKL on your system.
mkl_link = opti_FindMKL();
% Get Arch Dependent Library Path
libdir = opti_GetLibPath();

fprintf('\n------------------------------------------------\n');
fprintf('PSwarm MEX FILE INSTALL\n\n');

%Get PSwarm Libraries
post = [' -IInclude -L' libdir ' -llibpswarm -llibut -output pswarm'];
%Get MKL Libraries (for BLAS)
post = [post mkl_link];

%CD to Source Directory
cdir = cd;
cd 'Solvers/pswarm/Source';

%Compile & Move
pre = 'mex -v -largeArrayDims pswarmmex.c';
try
    eval([pre post])
    movefile(['pswarm.' mexext],'../','f')
    fprintf('Done!\n');
catch ME
    cd(cdir);
    error('opti:pswarm','Error Compiling PSwarm!\n%s',ME.message);
end
cd(cdir);
fprintf('------------------------------------------------\n');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [Z, Wsep, Wmix, m, U, S, V] = zca(X, portion)

if nargin < 2
    portion = 1;
end

rndidx = randperm(size(X, 1));
X_orig = X;
X = X(rndidx(1:round(size(X,1) * portion)), :);

m = mean(X, 1);
Xc = bsxfun(@minus, X, m);

sigma = Xc' * Xc / size(Xc, 1);

[U, S, V] = svd(sigma, 0);

Wsep = V / sqrt(S+1e-12);
Wmix = V';
Z = bsxfun(@minus, X_orig, m) * Wsep * Wmix;

author: kyunghyuncho

repo: deepmat


text:
function [status, fieldsUsed] = spm_matlab_version_chk(chk,tbx)
% Check a version number against a Toolbox version.
%    FORMAT [STATUS FIELDSUSED] = SPM_MATLAB_VERSION(CHK,TBX)
%
%    This function checks if a user supplied version number is less than,
%    equal to or greater than the version number of specified toolbox. If no
%    toolbox is specified the function checks the version of Matlab. User
%    defined toolboxes can be checked but the Contents.m file must conform
%    to the specification defined in ver.m
%
%    This function assumes that the version number is really a text string
%    with fields major.minor.revision.build. Other formats are not supported.
%    Checking is done to the level specified by the input version. Thus an
%    input of '7' will be rated as the same version as 7.1, but 7.0 would be
%    rated as earlier than 7.1.
%
%    INPUTS
%       chk         Version number to be checked [string].
%       tbx         Name of toolbox to check. Defaults to 'Matlab'.
%
%    OUPUTS
%       This is based on the fields the user supplies in the input.
%
%       status      Defines the outcome of the comparison
%            -1       Toolbox version is earlier than the user supplied version.
%             0       Toolbox and user versions are the same
%             1       Toolbox version is later than the user supplied version.
%                   Think of it this way, the sign of status is determined by
%                   MATLAB_TOOLBOX_VERSION - USER_VERSION (i.e., THE VERSION YOU
%                   INPUT).
%       fieldsUsed  Returns the version fields that were actually compared.
%
%    EXAMPLES
%       If the Matlab version is 7.1.0.83, and the user supplied version is '7':
%       status = spm_matlab_version_chk('7');
%       returns status == 0   : major revision numbers are the same.
%
%       If the Matlab version is 7.1.0.0, and the user supplied version is '7.1':
%       status = spm_matlab_version_chk('7');
%       returns status == 0   : major and minor revision numbers are the same.
%
%       If the Matlab version is 7.1.0.83, and the user supplied version is '7.2':
%       status = spm_matlab_version_chk('7.2');
%       returns status == -1   : major + minor revision is earlier for Matlab.
%
%       If the Matlab version is 6.5.1, and the user supplied version is '6.5.0'.
%       status = spm_matlab_version_chk('6.5.0');
%       returns status == 1     : major + minor + release revision is later
%                                 for matlab
%       The statement ( spm_matlab_version_chk('6.5.0') > 0 ) is true for
%       all Matlab Toolbox versions after 6.5.0.
%
%    See also  VERSION, VER.

%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Darren Gitelman
% $Id$

% output variable
%--------------------------------------------------------------------------
status = [];
fieldsUsed = [];

%  Check inputs
%--------------------------------------------------------------------------
% If no inputs then error; if 1 then set tbx to Matlab
%-------------------------------------------------------
switch nargin
    case 0
        error('Please provide a version number to be checked.');
    case 1
        tbx = 'Matlab';
end

% If a number is supplied then convert to text
%-------------------------------------------------------
if isnumeric(chk)
    chk = num2str(chk);
end

% If too many fields in input then error
%-------------------------------------------------------
if numel(findstr(chk,'.')) > 3
    error('Input string has too many fields. Only major.minor.release.build fields are supported.');
end

% parse the input string into a cell array of strings
%--------------------------------------------------------------------------
vCHK = textscan(chk,'%s %s %s %s','delimiter','.');

% find empty cells. These will be set to 0 in both cell arrays so that the
% same revision fields are compared
%--------------------------------------------------------------------------
emptyCHK = cellfun('isempty',vCHK);

% Get the requested toolbox version
%--------------------------------------------------------------------------
if strcmpi(tbx,'Matlab')
    % Matlab is a special case. The ver command does not report the
    % entire version string.
    [ tbxVer tmp ] = strtok(version);
else
    tbxStruct = ver(tbx);
    if numel(tbxStruct) >1
        error('Too many toolboxes found for given toolbox name.')
    end
    tbxVer = tbxStruct.Version;
end

% parse the requested toolbox version string into a cell array of strings.
%--------------------------------------------------------------------------
vTBX = textscan(tbxVer,'%s %s %s %s','delimiter','.');

% find empty cells. These will be removed from the cell array so that the
% same revision fields are compared
%--------------------------------------------------------------------------
emptyTBX = cellfun('isempty',vTBX);

% combine empty cell indices. notEmptyCells tell us which version
% fields are compared.
%--------------------------------------------------------------------------
emptyCellLogical = emptyCHK | emptyTBX;
emptyCells       = find(emptyCellLogical);
notEmptyCells    = find(~emptyCellLogical);

% and remove them.
%--------------------------------------------------------------------------
 vCHK(emptyCells) = [];
 vTBX(emptyCells) = [];

% final version fields to compare are converted to numbers
%--------------------------------------------------------------------------
vCHKMat = str2num(strvcat([vCHK{:}]))';
vTBXMat = str2num(strvcat([vTBX{:}]))';

% compare versions
%--------------------------------------------------------------------------
% This array will be used to decide which version is the later one. The
% differences between versions are computed on each element of the version
% number: major.minor.release.build and then multiplied by the base2array
% and summed. This is equivalent to logically combining the base 2 place
% values and ensures that the fields are evaluated by their significance.
% The sign of the summed values gives the result. Positive values mean the
% toolbox version is later, negative values mean the user supplied version,
% is later, and 0 means they are equal.
%--------------------------------------------------------------------------
base2array             = [8 4 2 1];
base2array(emptyCells) = [];

% take the difference between version elements.
%--------------------------------------------------------------------------
vDiff = vTBXMat - vCHKMat;

% get the sign of the differences
%--------------------------------------------------------------------------
signDiff = sign(vDiff);

% multiply by the base2array
%--------------------------------------------------------------------------
vVal = signDiff.*base2array;

% If the sign of the sum is positive then toolbox version is later; if 0
% they are the same, and if negative the user supplied version is later.
% Think of it this way: MATLAB_TOOLBOX_VERSION - USER_VERSION.
% Remember the comparison is made with respect to the FIELDS SUPPLED IN THE
% INPUT
%--------------------------------------------------------------------------
switch sign(sum(vVal))
    case -1     % MATLB_TOOLBOX_VERSION is EARLIER than USER_VERSION
        status = -1;
    case 0      % MATLB_TOOLBOX_VERSION = EQUALS = USER_VERSION
        status =  0;
    case 1      % MATLB_TOOLBOX_VERSION is LATER than USER_VERSION
        status =  1;
end

% also return the fields that were compared.
%--------------------------------------------------------------------------
fieldsUsed             = {'major','minor','revision','build'}';
fieldsUsed(emptyCells) = [];
return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef WindowData

  methods (Static)

    function wd = windowdata(nCls)
      % Windowdata constructor
      % nCls: number of classifiers
      % wd: struct array, indexed by classifier
      %
      %
      % .labelidx_cur
      % .labelidx_new
      % .labelidx_old: These three arrays take the same values as
      % Labels.labelidx.vals. A 0 value indicates 'unlabeled'; Otherwise
      % the values are label indices 1:nlabels.

      % ALTODO see below, NUA='not used anymore',
      % ALTODO: correct-sized empties (pass windowdataVerify)

      wd = struct(...
        'X',single([]), ...
        'exp',cell(nCls,1), ...
        'flies',[], ...
        't',[], ...
        'labelidx_cur',[], ...
        'labelidx_new',[], ...
        'labelidx_old',[], ...
        'labelidx_imp',[], ...
        'featurenames',{{}}, ...
        'predicted',[], ... % NUA
        'predicted_probs',[], ... % NUA
        'isvalidprediction',[], ... % NUA
        'distNdx',[], ...
        'scores',[], ... % NUA
        'scoreNorm',[], ...
        'binVals',[], ...
        'scores_old',[], ... % NUA
        'scores_validated',[], ...
        'postprocessed',[],... % NUA?
        'bins',[]);
    end

    function windowdataVerify(wd)
      for iCls = 1:numel(wd)
        w = wd(iCls);
        [n,nftr] = size(w.X);

        assert(isequal([n 1],size(w.exp),size(w.flies),size(w.t),...
          size(w.labelidx_cur),size(w.labelidx_new),size(w.labelidx_old),...
          size(w.labelidx_imp)) || ...
          isempty(w.exp) && isempty(w.flies) && isempty(w.t) && ... % AL: alternate condition here for legacy compat
          isempty(w.labelidx_cur) && isempty(w.labelidx_new) && ...
          isempty(w.labelidx_old) && isempty(w.labelidx_imp));

        assert(numel(w.featurenames)==nftr);
        if ~isempty(w.binVals),
          assert(size(w.bins,1) == nftr && size(w.bins,2) == n);
        end

        % scoreNorm?
        % binVals?
      end
    end

    function wd = windowdataClear(wd)
      for iCls = 1:numel(wd)
        wd(iCls).X = single([]);
        wd(iCls).exp = [];
        wd(iCls).flies = [];
        wd(iCls).t = [];
        wd(iCls).labelidx_cur = [];
        wd(iCls).labelidx_new = [];
        wd(iCls).labelidx_imp = [];
        wd(iCls).labelidx_old = [];
        wd(iCls).featurenames = {{}};

        wd(iCls).predicted = [];
        wd(iCls).predicted_probs = [];
        wd(iCls).isvalidprediction = [];
        wd(iCls).distNdx = [];
        wd(iCls).scores = [];
        wd(iCls).scoreNorm = [];
        wd(iCls).binVals = [];
        wd(iCls).bins = uint8([]);
        wd(iCls).scores_old = [];
        wd(iCls).scores_validated = [];
        wd(iCls).postprocessed  = [];
      end
    end

    function wd = windowdataTrim(wd,predFcn)
      assert(isa(predFcn,'function_handle'));

%       sizeLimit = obj.cacheSize*1e6;
%       classSize = 4;
%       ratioLimit = 0.2;
%
%       numUnlabeled = nnz(obj.windowdata.labelidx_new==0);
%       numLabeled = nnz(obj.windowdata.labelidx_new);
%
%       if (nargin < 2 || ~doforce) && (numel(obj.windowdata.X)*classSize < sizeLimit || numUnlabeled/numLabeled<ratioLimit);
%         return;
%       end
%
%       idx2remove = obj.windowdata.labelidx_new==0 & ...
%         ~obj.FlyNdx(obj.expi,obj.flies);

      for i = 1:numel(wd)
        n = size(wd(i).X,1);
        % MAYANK_JAN16_2015: Don't trim if there is no windowdata.
        if n == 0, continue ; end
        tfRm = predFcn(wd(i));
        assert(islogical(tfRm)&&isvector(tfRm)&&numel(tfRm)==n);
        % MAYANK MAY9 2016  don't trim if there is nothing to trim.
        if nnz(tfRm)<1, continue ;end
        wd(i).X(tfRm,:) = [];
        wd(i).exp(tfRm,:) = [];
        wd(i).flies(tfRm,:) = [];
        wd(i).t(tfRm,:) = [];
        wd(i).bins(:,tfRm) = [];
        wd(i).labelidx_cur(tfRm,:) = [];
        wd(i).labelidx_new(tfRm,:) = [];
        wd(i).labelidx_imp(tfRm,:) = [];
        wd(i).labelidx_old(tfRm,:) = [];
        wd(i).distNdx = [];
%        wd(i).binVals = [];
        wd(i).scores_validated(tfRm,:) = [];

%         wd(i).isvalidprediction(idxcurr(1:numel(wd(i).isvalidprediction))) = [];
%         wd(i).predicted(idxcurr(1:numel(wd(i).predicted))) = [];
%         wd(i).predicted_probs(idxcurr(1:numel(wd(i).predicted_probs))) = [];
%         wd(i).scores(idxcurr(1:numel(wd(i).scores))) = [];
%         wd(i).scores_old(idxcurr(1:numel(wd(i).scores_old))) = [];
%         wd(i).postprocessed(idxcurr(1:numel(wd(i).postprocessed)),:) = [];
      end
    end

    function wd = windowdataSetFeaturenames(wd,featurenames)
      assert(iscell(featurenames)&&numel(featurenames)==numel(wd));
      for i = 1:numel(wd)
        nFtr = numel(featurenames{i});
        if isequal(wd(i).X,[])
          wd(i).X = zeros(0,nFtr);
        else
          nCol = size(wd(i).X,2);
          assert(nCol==nFtr,...
            'Number of featurenames (%d) does not match number of cols in feature matrix (%d).',...
            nFtr,nCol);
        end
        wd(i).featurenames = featurenames{i};
      end
    end

    function wd = windowdataRemapLabelIdxs(wd,oldIdx2NewIdx)
      % oldIdx2NewIdx: vector of length number-of-old-labels.
      % oldIdx2NewIdx(i) contains the new label index corresponding to i.

      LABEL_FIELDS = {'labelidx_cur' 'labelidx_new' 'labelidx_old'};
      for iCls = 1:numel(wd)
        w = wd(iCls);
        for fld = LABEL_FIELDS, fld=fld{1}; %#ok<FXSET>
          tfnz = w.(fld)~=0;
          newNZVals = oldIdx2NewIdx(w.(fld)(tfnz));
          assert(all(newNZVals>0));
          w.(fld)(tfnz) = newNZVals;
        end
        wd(iCls) = w;
      end
    end

    % MK: Add thresholded bins to window data. May 10 2016
    function [wd,wdmodified ] = modernize(wd)
      wdmodified = false;
      if ~isfield(wd,'bins'),
        if numel(wd)>0,
          for ndx = 1:numel(wd),
            if ~isempty(wd(ndx).binVals),
              wd(ndx).bins = findThresholdBins(wd(ndx).X,wd(ndx).binVals);
            else
              wd(ndx).bins = uint8([]);
            end
          end
        else
          tmp=cell(size(wd)); [wd(:).bins]=deal(tmp{:});
        end
        wdmodified = true;
      end

    end

  end

end

author: kristinbranson

repo: JAABA


text:
function [mCatDecluster, mCatAfter, vCluster, vCl, vMainCluster] = calc_decluster_depth(mCatalog,nMethod)
% function [mCatDecluster, mCatAfter, vCluster, vCl, vMainCluster] = calc_decluster_depth(mCatalog,nMethod)
% ---------------------------------------------------------------------------------------------------------
% Function to decluster earthquake catalog using a windowing technique in space and time by
% Knopoff & Gardner, GJR astr. Soc, 28, 311-313, 1972 and Gardner & Knopoff, BSSA, 64,5, 1363-1367, 1974
% using different windows. Spatially, the algortihm uses a sphere to determine events.
%
% Incoming variables
% mCatalog : Incoming earthquake catalog (ZMAP format)
% nMethod  : Window length for declustering (see calc_windows.m)
%            1: Gardener & Knopoff, 1974
%            2: Gruenthal pers. communication
%            3: Urhammer, 1986
%
% Outgoing variables:
% mCatDecluster : Declustered earthquake catalog
% mCatAfter     : Catalog of aftershocks (and foreshocks)
% vCluster      : Vector indicating only aftershocks/foreshocls in cluster using a cluster number
% vCl           : Vector indicating all events in clusters using a cluster number
% vMainCluster  : Vector indicating mainshocks of clusters using a cluster number
%
% J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 29.08.02

%% Added:
% 31.07.02 Correction for problem of mainshocks with a cluster number as aftershocks belong to two sequences
% 31.07.02 Corrected fMaxClusterMag(nMagCount) to fMaxClusterMag, since counting variable not needed
% 31.07.02 Improved resizing time window by adding time difference from initial event to bigger aftershock
% 13.08.02 Added waitbars
% 28.08.02 Changed distance determination using now distance and repmat
% 29.08.02 Cluster determination strategy change: Now selecting all aftershocks using the window of the first shock,
%          adding the events from the bigger aftershocks (later labelled mainshocks); calc_decluster_ver3.m keeps
%          resizing technique

%%% Remember: Improve zero length cluster problem which might appear

%% Initialize Vectors and Matrices
mCatDecluster = [];
mCatAfter = [];
vCluster = zeros(length(mCatalog),1); % Initialize all events as mainshock
vCl = zeros(length(mCatalog),1); % Initialize all events as mainshock
vSel = zeros(length(mCatalog),1); % Initialize all events as mainshock
vMainCluster = zeros(length(mCatalog),1); % Initialize

[nXSize, nYsize] = size(mCatalog);
if nXSize == 0
    disp('Load new catalog');
    return
end

vDecDate = mCatalog(:,3);
nCount = 0;    % Variable of cluster number

fMagThreshold = min(mCatalog(:,6)); % Set Threshold to minimum magnitude of catalog
hWaitbar1 = waitbar(0,'Identifying clusters...');
set(hWaitbar1,'Numbertitle','off','Name','Decluster percentage');
for nEvent=1:length(mCatalog(:,6))
    %nEvent
    %nCount
    if vCluster(nEvent) == 0
        fMagnitude(nEvent) = mCatalog(nEvent, 6);
        if fMagnitude(nEvent) >= fMagThreshold
            %% Define first aftershock zone and determine magnitude of strongest aftershock
            fMag = fMagnitude(nEvent);
            [fSpace, fTime] = calc_windows(fMagnitude(nEvent), nMethod);
            fSpaceDeg = km2deg(fSpace);
            %% This first if is for events with no location given
            if isnan(mCatalog(nEvent,1))
                %vSel = (vDecDate(:,1)-vDecDate(nEvent,1) >= 0) & (vDecDate(:,1)-vDecDate(nEvent,1) <= fTime  & vCluster(nEvent) == 0);
                vSel = (mCatalog(:,3) == mCatalog(nEvent,3));
            else
                mPos = [mCatalog(nEvent, 1) mCatalog(nEvent,2)];
                mPos = repmat(mPos,length(mCatalog(:,1)), 1);
                mDist = abs(distance(mCatalog(:,1), mCatalog(:,2), mPos(:,1), mPos(:,2)));
                vDepthDist = abs(mCatalog(nEvent,7)-mCatalog(:,7));
                vSel = ((mDist <= fSpaceDeg) & vDepthDist <= fSpace & (vDecDate(:,1)-vDecDate(nEvent,1) >= 0) &...
                    (vDecDate(:,1)-vDecDate(nEvent,1) <= fTime) & vCluster(nEvent) == 0);
            end% End of isnan(mCatalog)
            mTmp = mCatalog(vSel,:);
            if length(mTmp(:,1)) == 1  % Only one event thus no cluster; IF to determine cluster or not
                fMaxClusterMag = fMag;
            else
                fMaxClusterMag = max(mTmp(:,6));
                [nIndiceMaxMag] = find(mTmp(:,6) == fMaxClusterMag);
                fTimeMaxClusterMag = mTmp(max(nIndiceMaxMag),3);
                % Search for event with bigger magnitude in cluster and add to cluster
                while fMaxClusterMag-fMag > 0
                    [fSpace, fTime] = calc_windows(fMaxClusterMag, nMethod);
                    fSpaceDeg = km2deg(fSpace);
                    %% Adding aftershocks from bigger aftershock
                    mPos = [mTmp(min(nIndiceMaxMag),1) mTmp(min(nIndiceMaxMag),2)];
                    mPos = repmat(mPos,length(mCatalog(:,1)), 1);
                    mDist = abs(distance(mCatalog(:,1), mCatalog(:,2), mPos(:,1), mPos(:,2)));
                    vDepthDist = abs(mTmp(min(nIndiceMaxMag),7)-mCatalog(:,7));
                    vSel2 = ((mDist <= fSpaceDeg) & vDepthDist <= fSpace & (vDecDate(:,1)-mTmp(min(nIndiceMaxMag),3) >= 0) &...
                        (vDecDate(:,1)-mTmp(min(nIndiceMaxMag),3) <= fTime) & vCluster == 0);
                    mTmp = mCatalog(vSel2,:);
                    vSel = (vSel > 0 | vSel2 > 0); % Actual addition
                    if isempty(mTmp) % no events in aftershock zone
                        break;
                    end
                    fMag = fMaxClusterMag;
                    fMaxClusterMag = max(mTmp(:,6));
                    [nIndiceMaxMag] = find(mTmp(:,6) == fMaxClusterMag);
                    fTimeMaxClusterMag = mTmp(max(nIndiceMaxMag),3);
                    if fMaxClusterMag - fMag == 0 % no bigger event in aftershock zone
                        break;
                    end
                end  % End of while
                nCount = nCount + 1; % Set cluster number
            end % End of if length(mTmp)

            [vIndice]=find(vSel); % Vector of indices with Clusters
            vTmpCluster(vIndice,:) = nCount;
            %length(vTmpCluster(vIndice,:));
            nI=1; % Variable counting the length of the cluster
            % Select the right numbers for the cluster using the indice vector vIndice
            % First: Insert cluster number after check for length
            % Second: Check if it's a mainshock
            % Third: Keep the former cluster indice;
            while nI <= length(vIndice)
                if (~isempty(vTmpCluster(vIndice(nI))) & length(vTmpCluster(vIndice,:)) > 1 & vCluster(vIndice(nI)) == 0)
                    vCluster(vIndice(nI)) = vTmpCluster(vIndice(nI));
                    %vEventnr(vIndice,:) = nEvent;
                elseif  (~isempty(vTmpCluster(vIndice(nI)))  &&  length(vTmpCluster(vIndice,:)) == 1  &&  vCluster(vIndice(nI)) == 0)
                    vCluster(vIndice(nI)) = 0;
                else
                    vCluster(vIndice(nI)) = vCluster(vIndice(nI));
                end
                nI=nI+1;
            end %End of while nI
            %                 nCount = nCount + 1; % Set cluster number %% Watch
            %             end % End of if to determine cluster or not %% Watch
            %%% Check if the Cluster is not just one event which can happen in case of keeping the former
            %%% cluster number in preceeding while-Loop
            vSelSingle = (vCluster == nCount);
            [vIndiceSingle] = find(vSelSingle);
            %vTmpSingle(vIndiceSingle,:);
            if length(vIndiceSingle) == 1
                %nCount
                %vIndiceSingle
                vCluster(vIndiceSingle)=0; % Set the event as mainsock
                nCount = nCount-1; % Correct the cluster number down by one
            end
        end % End of if checking magnitude threshold fMagThreshold
    end % End of if checking if vCluster == 0
    if rem(nEvent,100) == 0
        waitbar(nEvent/length(mCatalog(:,6)))
    end % End updating waitbar
end % End of FOR over mCatalog
close(hWaitbar1);
%nCount
%% vCL Cluster vector with mainshocks in it; vCluster is now modified to get rid of mainshocks
vCl = vCluster;

%% Matrix with cluster indice, magnitude and time
mTmpCat = [vCluster mCatalog(:,6) mCatalog(:,3)];
%% Delete largest event from cluster series and add to mainshock catalog
hWaitbar2 = waitbar(0,'Identifying mainshocks in clusters...');
set(hWaitbar2,'Numbertitle','off','Name','Mainshock identification ');
for nCevent = 1:nCount
    %nCevent
    vSel4 = (mTmpCat(:,1) == nCevent); % Select cluster
    mTmpCat2 = mCatalog(vSel4,:);
    fTmpMaxMag = max(mTmpCat2(:,6)); % Select max magnitude of cluster
    vSelMag = (mTmpCat2(:,6) == fTmpMaxMag);
    [nMag] = find(vSelMag);
    if length(nMag) == 1
        vSel5 = (mTmpCat(:,1) == nCevent & mTmpCat(:,2) == fTmpMaxMag); % Select the event
        [vIndiceMag] = find(vSel5); % Find indice
        vCluster(vIndiceMag) = 0;  % Set cluster value to zero, so it is a mainshock
        vMainCluster(vIndiceMag) = nCevent; % Set mainshock vector to cluster number
    elseif isempty(nMag)
        disp('Nothing in ')
        nCevent
    else
        vSel = (mTmpCat(:,1) == nCevent & mTmpCat(:,2) == fTmpMaxMag);
        mTmpCat3 = mCatalog(vSel,:);
        [vIndiceMag] = min(find(vSel)); % Find minimum indice of event with max magnitude in cluster
        vCluster(vIndiceMag) = 0;  % Set cluster value to zero, so it is a mainshock
        vMainCluster(vIndiceMag) = nCevent;  % Set mainshock vector to cluster number
    end
    if rem(nCevent,20) == 0
        waitbar(nCevent/nCount)
    end % End updating waitbar
end % End of For nCevent
close(hWaitbar2);
%% Create a catalog of aftershocks (mCatAfter) and of declustered catalog (mCatDecluster)
vSel = (vCluster(:,1) > 0);
mCatDecluster=mCatalog(~vSel,:);
mCatAfter = mCatalog(vSel,:);

author: CelsoReyes


text:
function [VMap,NMap,tMap,CMap] = raycastingTSDF(camRtC2W, castingRange)

% VMap and NMap are in the world coordinate

% DDA based ray casting

global raycastingDirectionC;
global voxel;
global tsdf_value;
global tsdf_color;

camCenterW = getCameraCenter(camRtC2W);

raycastingDirectionW = transformRTdir(raycastingDirectionC,camRtC2W);

camCenterWgrid = (camCenterW - voxel.range(:,1)) / voxel.unit + 1;

%camCenterWcopy = repmat( (camCenterW - voxel.range(:,1) ) / voxel.unit + 1,1,640*480);
%startPoints = camCenterWcopy + raycastingDirectionW * (castingRange(1)/ voxel.unit);
%endPoints = camCenterWcopy + raycastingDirectionW * (castingRange(2)/ voxel.unit);

% http://stellar.mit.edu/S/course/6/fa10/6.837/courseMaterial/topics/topic1/lectureNotes/13_RayTracing-Acceleration/13_RayTracing-Acceleration.pdf

% we assume the camera must be inside. otherwise, it is an error.
% so we don't test the camera

raycastingDirectionWinv = raycastingDirectionW.^-1;

maxTx = max((voxel.size_grid(1)-2-camCenterWgrid(1))*raycastingDirectionWinv(1,:),(2-camCenterWgrid(1))*raycastingDirectionWinv(1,:));
maxTy = max((voxel.size_grid(2)-2-camCenterWgrid(2))*raycastingDirectionWinv(2,:),(2-camCenterWgrid(2))*raycastingDirectionWinv(2,:));
maxTz = max((voxel.size_grid(3)-2-camCenterWgrid(3))*raycastingDirectionWinv(3,:),(2-camCenterWgrid(3))*raycastingDirectionWinv(3,:));
maxT = min(maxTx, min(maxTy, maxTz));

castingRangeGrid = castingRange/voxel.unit;
maxT = min(maxT, castingRangeGrid(2));

% parallel setting
parallelBlock = matlabpool('size');
blockDim = (640*480)/parallelBlock;
if blockDim ~= round(blockDim)
    error('thread does not align');
end

backMove = voxel.mu_grid;

tMap = NaN(1,640*480);
NMap = NaN(3,640*480);
CMap = NaN(3,640*480);

initDis = 3/ voxel.unit;

% http://www.cse.yorku.ca/~amana/research/grid.pdf
prevT = initDis;

for i=1:(640*480)
    % Ray-Box Intersection to get the range

    tMin = castingRangeGrid(1);
    tMax = maxT(i);

    if tMax>tMin % max range > min range

        tStart = max(castingRangeGrid(1),min(tMax,prevT-backMove));

        rayDir = raycastingDirectionW(:,i);

        startPoint = camCenterWgrid + rayDir*tStart;

        % The initialization phase begins by identifying the voxel in which the ray origin, ?u, is found.
        % If the ray origin is outside the grid, we ?nd the point in which the ray enters the grid and take the adjacent voxel.
        % The integer variables X and Y are initialized to the starting voxel coordinates.
        X = round(startPoint(1));
        Y = round(startPoint(2));
        Z = round(startPoint(3));

        valCurrentVoxel = tsdf_value(X,Y,Z);

        tOptimal = NaN;

        if valCurrentVoxel == 0
            % lucky, you got the surface immediately!
            %tMap(i) = tStart;
            prevT = tStart;

            tOptimal = tStart;
        else
            if valCurrentVoxel>0
                localDir = rayDir;
                tMax = tMax-tStart;
                lookforSign = -1;
            else
                localDir = -rayDir;
                tMax = tStart-tMin;
                lookforSign = +1;
            end

            % In addition, the variables stepX and stepY are initialized to either 1 or -1 indicating
            % whether X and Y are incremented or decremented as the ray crosses voxel boundaries
            % (this is determined by the sign of the x and y components of ?v).
            stepX = sign(localDir(1));
            stepY = sign(localDir(2));
            stepZ = sign(localDir(3));

            % Next, we determine the value of t at which the ray crosses the ?rst vertical voxel boundary
            % and store it in variable tMaxX. We perform a similar computation in y and store the result in tMaxY.
            % The minimum of these two values will indicate how much we can travel along the ray and still remain in the current voxel.

            if localDir(1)>0
                tMaxX = (X+0.5 - startPoint(1))/localDir(1);
            elseif localDir(1)<0
                tMaxX = (X-0.5 - startPoint(1))/localDir(1);
            else
                tMaxX = Inf;
            end
            if localDir(2)>0
                tMaxY = (Y+0.5 - startPoint(2))/localDir(2);
            elseif localDir(2)<0
                tMaxY = (Y-0.5 - startPoint(2))/localDir(2);
            else
                tMaxY = Inf;
            end
            if localDir(3)>0
                tMaxZ = (Z+0.5 - startPoint(3))/localDir(3);
            elseif localDir(3)<0
                tMaxZ = (Z-0.5 - startPoint(3))/localDir(3);
            else
                tMaxZ = Inf;
            end

            if min(min(tMaxX,tMaxY),tMaxZ)<0
                error('tStart<0');
            end

            % Finally, we compute tDeltaX and tDeltaY.
            % tDeltaX indicates how far along the ray we must move (in units of t) for the horizontal component of such a movement to equal the width of a voxel.
            % Similarly,we store in tDeltaY the amount of movement along the ray which has a vertical component equal to the height of a voxel.
            tDeltaX = 1/abs(localDir(1));
            tDeltaY = 1/abs(localDir(2));
            tDeltaZ = 1/abs(localDir(3));

            t = 0;
            while t<tMax
                if tMaxX < tMaxY
                    if tMaxX < tMaxZ
                        X= X + stepX;
                        tMaxX= tMaxX + tDeltaX;
                    else
                        Z= Z + stepZ;
                        tMaxZ= tMaxZ + tDeltaZ;
                    end
                else
                    if tMaxY < tMaxZ
                        Y= Y + stepY;
                        tMaxY= tMaxY + tDeltaY;
                    else
                        Z= Z + stepZ;
                        tMaxZ= tMaxZ + tDeltaZ;
                    end
                end
                prevT = t;
                t = min(min(tMaxX,tMaxY),tMaxZ);

                if tsdf_value(X,Y,Z) * lookforSign > 0

                    if (tsdf_value(X,Y,Z)==-1 && valCurrentVoxel==+1) || (tsdf_value(X,Y,Z)==+1 && valCurrentVoxel==-1)
                        tOptimal = NaN;
                    else

                        % found the zero crossing points!

                        % simplest
                        % tOptimal = t;

                        % simple average
                        tOptimal = (prevT + t)/2;

                        %{
                        % buggy: Ftdt == Ft

                        if lookforSign>0
                            prevTswap = prevT;
                            prevT = t;
                            t= prevTswap;
                        end
                        XYZt   = startPoint + localDir*t;
                        XYZtdt = startPoint + localDir*prevT;

                        Ft = interpolateTrilineary(XYZt(1),XYZt(2),XYZt(3));
                        Ftdt = interpolateTrilineary(XYZtdt(1),XYZtdt(2),XYZtdt(3));

                        tOptimal = t - abs(t-prevT) * Ft / (Ftdt - Ft);
                        if tOptimal>1000
                            error('tOptimal>1000');
                        end
                        %}

                        %fprintf('i=%d: t=%f tOptimal=%f Ft=%f Ftdt=%f\n',i,t,t - abs(t-prevT) * Ft / (Ftdt - Ft), Ft, Ftdt);
                        tOptimal = tStart - tOptimal * lookforSign;

                        prevT = tOptimal;

                    end
                    break;
                end

                valCurrentVoxel = tsdf_value(X,Y,Z);
            end

        end
        if ~isnan(tOptimal)
            tMap(i) = tOptimal;

            % compute normal map

            XYZgrid = camCenterWgrid + rayDir*tOptimal;

            NMap(1,i) = interpolateTrilineary(XYZgrid(1)+1,XYZgrid(2),XYZgrid(3))-interpolateTrilineary(XYZgrid(1)-1,XYZgrid(2),XYZgrid(3));
            NMap(2,i) = interpolateTrilineary(XYZgrid(1),XYZgrid(2)+1,XYZgrid(3))-interpolateTrilineary(XYZgrid(1),XYZgrid(2)-1,XYZgrid(3));
            NMap(3,i) = interpolateTrilineary(XYZgrid(1),XYZgrid(2),XYZgrid(3)+1)-interpolateTrilineary(XYZgrid(1),XYZgrid(2),XYZgrid(3)-1);

            if ~isempty(tsdf_color)
                CMap(:,i) = interpolateTrilinearyColor(XYZgrid(1),XYZgrid(2),XYZgrid(3));
            end
        end
    end
end

% normalize normal map
NMap = NMap ./ repmat(sqrt(sum(NMap.^2,1)),3,1);

% computer vertex map
VMap = repmat(camCenterW,1,640*480) + raycastingDirectionW .* (repmat(tMap*voxel.unit,3,1));

%imagesc(reshape(tMap,480,640)); axis equal; axis tight

author: jianxiongxiao

repo: ProfXkit


text:
%% This script simulates the Lyapunov Optimization-based Dynamic Computation Offloading (LODCO) algorithm.
%  author: Hailiang Zhao
clc, clear
opt = optimset('Display', 'none');

%% basic parameter settings (had better not change those paras)
k = 1e-28;                % effective switched capacitance (a constant decided by the chip architecture)
tau = 0.002;              % the length of time slot (in second)
phi = 0.002;              % the cost of task dropping (in second)
omega = 1e6;              % the bandwidth of MEC server (in Hz)
sigma = 1e-13;            % the noise power of the receiver (in W)
p_tx_max = 1;             % the maximum transmit power of mobile device (in W)
f_max = 1.5e9;            % the maximum CPU-cycle frequency of mobile device (in Hz)
E_max = 0.002;            % the maximum amout of battery output energy (in J)
L = 1000;                 % the input size of the computation task (in bit)
X = 737.5;                % the number of CPU cycles needed on processing one bit of task
W = L * X;                % the number of CPU cycles needed on processing one task
E_H_max = 48e-6;          % the upper bound of the energy arrive at the mobile device (in J)
p_H = E_H_max / (2*tau);  % the average Energy Harvesting (EH) power (in W)
g0 = power(10, -4);       % the path-loss constant

%% parameter control
T = 50000;                % the number of time slot (a.k.a. the size of the time horizon)
tau_d = 0.002;            % execution deadline (in second)
d = 50;                   % the distance between the mobile device and the MEC server (in meter)
E_min = 0.02e-3;          % the minimum amout of battery output energy (in J)
V = 1e-5;                 % the weight of penalty (the control parameter introduced by Lyapunov Optimization)
rho = 0.6;                % the probability that the computation task is requested

% the lower bound of perturbation parameter
E_max_hat = min(max(k * W * (f_max)^2, p_tx_max * tau), E_max);
theta = E_max_hat + V * phi / E_min;

%% allocate storage for valuable results
B = zeros(T, 1);          % the battery energy level (in J)
B_hat = zeros(T, 1);      % the virtual battery energy level ($B_hat = B - theta$)
e = zeros(T, 1);          % the amout of the harvested and stored energy (in J)
chosen_mode = zeros(T, 1);% {1: local, 2: remote, 3: drop, 4: no task request}
f = zeros(T, 1);          % the CPU-cycle frequency of local execution (in Hz)
p = zeros(T, 1);          % the transmit power of computation offloading (in W)
cost = zeros(T, 3);       % execution delay for mobile execution, MEC server execution and final choice, respectively (in second)
E = zeros(T, 3);          % energy consumption for mobile execution, MEC server execution and final choice, respectively (in J)

%% simulation begin
t = 1;
while t <= T
    disp(['===> Time slot #', num2str(t), ' <==='])

    %% initialization
    % generate the task request
    zeta = binornd(1, rho);
    % generate the virtual battery energy level
    B_hat(t) = B(t) - theta;

    %% step 1: get the optimal energy harvesting no matter whether task is requested
    E_H_t = unifrnd(0, E_H_max);
    if B_hat(t) <= 0
        e(t) = E_H_t;
    end

    %% step 2: get the optimal computation offloading strategy (I_m, I_s, I_d, f(t), p(t))
    if zeta == 0
        % chosen mode has to be 4
        disp('no task request generated!')
        chosen_mode(t) = 4;
    else
        % chosen_mode is chosen from {1, 2, 3}
        disp('task request generated!')
        % task request exists, generate the channel power gain
        h = exprnd(g0 / power(d, 4));

        %% step 2.1: solve the optimization problem $\mathcal{P}_{ME}$ (f(t) > 0)
        % calculate f_L and f_U
        f_L = max(sqrt(E_min / (k * W)), W / tau_d);
        f_U = min(sqrt(E_max / (k * W)), f_max);
        if f_L <= f_U
            % the sub-problem is feasible
            disp('mobile execution ($\mathcal{P}_{ME}$) is feasible!')

            if B_hat(t) < 0
                f_0 = (V / (-2 * B_hat(t) * k))^(1/3);
            else
                % complex number may exist, which may lead to error
                f_0 = -(V / (2 * B_hat(t) * k))^(1/3);
            end

            if (f_0 > f_U && B_hat(t) < 0) || (B_hat(t) >= 0)
                f(t) = f_U;
            elseif f_0 >= f_L && f_0 <= f_U && B_hat(t) < 0
                f(t) = f_0;
            elseif f_0 < f_L && B_hat(t) < 0
                f(t) = f_L;
            end
            % check whether f(t) is zero
            if f(t) == 0
                disp('Something wrong! f is 0!')
            end

            % calculate the delay of mobile execution
            cost(t, 1) = W / f(t);
            % calculate the energy consumption of mobile execution
            E(t, 1) = k * W * (f(t)^2);
            % calculate the value of optimization goal
            J_m = -B_hat(t) * k * W * (f(t))^2 + V * W / f(t);
        else
            % the sub-problem is not fasible because (i) the limited
            % computation capacity or (ii) time cosumed out of deadline or
            % (iii) the energy consumed out of battery energy level
            % If it is not feasible, it just means that we cannot choose
            % 'I_m=1'. It dosen't mean that the task has to be dropped.
            disp('mobile execution ($\mathcal{P}_{ME}$) is not feasible!')
            f(t) = 0;
            cost(t, 1) = 0;
            E(t, 1) = 0;
            % 'I_m=1' can never be chosen if mobile execution goal is inf
            J_m = inf;
        % Attention! We do not check whether the energy cunsumed is larger than
        % battery energy level because the problem $\mathcal{J}_{CO}$ does
        % not have constraint (8).
        end

        %% step 2.2: solve the optimization problem $\mathcal{P}_{SE}$ (p(t) > 0)
        E_tmp = sigma * L * log(2) / (omega * h);
        p_L_taud = (power(2, L / (omega * tau_d)) - 1) * sigma / h;
        % calculate p_L
        if E_tmp >= E_min
            p_L = p_L_taud;
        else
            % calculate p_E_min (use inline function and fsolve)
            y = @(x) x * L - omega * log2(1 + h*x/sigma) * E_min;
            % accroding to the function figure, p_L_taud is a positive
            % number around 0.2
            p_E_min = fsolve(y, 0.2, opt);
            p_L = max(p_L_taud, p_E_min);
        end
        % calculate p_U
        if E_tmp >= E_max
            p_U = 0;
        else
            % caculate p_E_max (use inline function and fsolve)
            y = @(x) x * L - omega * log2(1 + h*x/sigma) * E_max;
            % accroding to the function figure, p_E_max is a large positive
            % number around 20
            p_E_max = fsolve(y, 100, opt);
            p_U = min(p_tx_max, p_E_max);
        end

        if p_L <= p_U
            % the sub-problem is feasible
            disp('MEC server execution ($\mathcal{P}_{SE}$) is feasible!')
            % calculate p_0
            virtual_battery = B_hat(t);
            y = @(x) virtual_battery * log2(1 + h*x/sigma) + ...
                h * (V - virtual_battery*x) / log(2) / (sigma + h*x);
            p_0 = fsolve(y, 0.5, opt);

            if (p_U < p_0 && B_hat(t) < 0) || B_hat(t) >= 0
                p(t) = p_U;
            elseif p_0 < p_L && B_hat(t) < 0
                p(t) = p_L;
            elseif p_0 >= p_L && p_0 <= p_U && B_hat(t) < 0
                p(t) = p_0;
            end
            % check whether p(t) is zero
            if p(t) == 0
                disp('Something wrong! p is 0!')
            end

            % calculate the delay of MEC server execution
            cost(t, 2) = L / (omega * log2(1 + h*p(t)/sigma));
            % calculate the energy consumption of MEC server execution
            E(t, 2) = p(t) * cost(t, 2);
            % calculate the value of optimization goal
            J_s = (-B_hat(t) * p(t) + V) * cost(t, 2);
        else
            % the sub-problem is not feasible because (i) the limited transmit
            % power or (ii) time cosumed out of deadline or (iii) the energy
            % consumed out of battery energy level
            % If it is not feasible, it just means that we cannot choose
            % 'I_s=1'. It dosen't mean that the task has to be dropped.
            disp('MEC server execution ($\mathcal{P}_{SE}$) is not feasible!')
            p(t) = 0;
            cost(t, 2) = 0;
            E(t, 2) = 0;
            % 'I_s=1' can never be chosen if MEC server execution goal is inf
            J_s = inf;
        % Similarly, we do not check whether the energy cunsumed is larger than
        % battery energy level because the problem $\mathcal{J}_{CO}$ does
        % not have constraint (8).
        end

        %% step 3: choose the best execution mode
        J_d = V * phi;
        disp(['J_m:', num2str(J_m)])
        disp(['J_s:', num2str(J_s)])
        [~, mode] = min([J_m, J_s, J_d]);
        chosen_mode(t) = mode;
    end

    %% step 4: according to the chosen execution mode, calculate the real dealy and energy consumption
    if chosen_mode(t) == 1
        % mobile execution is chosen
        cost(t, 3) = cost(t, 1);
        E(t, 3) = E(t, 1);
    elseif chosen_mode(t) == 2
        % MEC server execution is chosen
        cost(t, 3) = cost(t, 2);
        E(t, 3) = E(t, 2);
    elseif chosen_mode(t) == 3
        % task is dropped, the delay is the task dropping penalty and the
        % energy consumption is zero
        cost(t, 3) = phi;
        E(t, 3) = 0;
    else
        % no task is requested, the delay and the energy consumption are
        % both zero
        cost(t, 3) = 0;
        E(t, 3) = 0;
    end

    %% step 5: update the battery energy level and go to next time slot
    B(t + 1) = B(t) - E(t, 3) + e(t);
    t = t + 1;
end

%% step 6: evaluate the simulation results
% 1. the battery energy level vs. time slot
figure
plot(1:T, B(1:T));
hold on
plot(1:T, repmat(theta + E_H_max, [T, 1]), '-')
title('Envolution of battery energy level')
xlabel('time slot')
ylabel('battery energy level $B_t$', 'Interpreter','latex')

% 2. the average execution cost vs. time slot
accumulated = 0;
average_cost = zeros(T, 1);
request_num = 0;
for t = 1: T
    accumulated = accumulated + cost(t, 3);
    if cost(t, 3) ~= 0
        % there exists task request
        request_num = request_num + 1;
    end
    average_cost(t) = accumulated / request_num;
end
figure
plot(1:T, average_cost);
title('Envolution of average execution cost')
xlabel('time slot')
ylabel('average execution cost $\frac{1}{T} \sum_{t=0}^{T-1} cost^t$', 'Interpreter','latex')

% 3. the average ratio of each chosen mode vs. time slot
average_ratio = zeros(T, 3);
mobile_exe = 0; server_exe = 0; drop = 0;
request_num = 0;
for t = 1: T
    if cost(t, 3) == 0
        continue
    else
        request_num = request_num + 1;
        if chosen_mode(t) == 1
            mobile_exe = mobile_exe + 1;
        elseif chosen_mode(t) == 2
            server_exe = server_exe + 1;
        else
            drop = drop + 1;
        end
    end
    average_ratio(t, :) = [mobile_exe, server_exe, drop] / request_num;
end
figure
plot(1:T, average_ratio(:, 1));
hold on
plot(1:T, average_ratio(:, 2));
hold on
plot(1:T, average_ratio(:, 3));
legend('mobile execution', 'MEC server execution', 'drop')
title('Envolution of average ratio of chosen modes')
xlabel('time slot')
ylabel('average  ratio of chosen modes $\frac{1}{T} \sum_{t=0}^{T-1} \{I_m^t, I_s^t, I_d^t\}$', 'Interpreter','latex')

author: hliangzhao

repo: Edge-Computing-Codes


text:
function value = orth_symm_condition ( n )

%*****************************************************************************80
%
%% ORTH_SYMM_DETERMINANT returns the L1 condition of the ORTH_SYMM matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real VALUE, the L1 condition.
%
  a_norm = 0.0;
  j = 1;
  for i = 1 : n
    angle = 2.0 * i * j * pi / ( 2 * n + 1 );
    a_norm = a_norm + 2.0 * abs ( sin ( angle ) ) / sqrt ( 2 * n + 1 );
  end
  b_norm = a_norm;
  value = a_norm * b_norm;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function BS2 = fit(v,varargin)
% function to fit Bingham parameters
%
% Description
% confidence ellipse for the mean direction based on Tanaka
% (1999) https://doi.org/10.1186/BF03351601
%
% Syntax
%   BS2 = BinghamS2.fit(v)
%
% Input
%  v - vector3d
%
% Output
%  BS2 - @BinghamS2
%
%
% Options
%  ConfElli - confidence level p (default at 0.95)
%
%
% Example
%
%   % simulate some directions
%   odf = unimodalODF(quaternion.id,'halfwidth',10*degree);
%   N = 100;
%   v = odf.discreteSample(N) .* ...
%     rotation.byAxisAngle(vector3d.X,rand(N,1)*2*pi) * vector3d.Y;
%
%   % fit a Bingham distribution
%   S2F = BinghamS2.fit(v)
%
%   % visualization
%   plot(S2F)
%   mtexColorMap LaboTeX
%   hold on
%   plot(t,'Markercolor','k','MarkerSize',3)
%   hold off
%

[a,kappa] = eig3(v*v);
kappa = kappa./sum(kappa);
Z =estimateZ(kappa);
BS2 = BinghamS2(Z, a);
BS2.N = BS2.normalizationConst;

% add the estimate of confidence level, given as ellipse half
% axes e.g.
% plot(v)
% ellipse(rotation('matrix',BS2.a.xyz'),BS2.cEllipse(1),BS2.cEllipse(2))
p = get_option(varargin,'ConfElli',0.95);
J = sqrt(chi2inv(p,2))/2;
BS2.cEllipse = [J/(-Z(2)*(kappa(3)-kappa(2))), ...
  J/(-Z(1)*(kappa(3)-kappa(1)))];

  function Z = estimateZ(kappa)
    % adapted from https://github.com/libDirectional/libDirectional
    % Igor Gilitschenski, Gerhard Kurz, Simon J. Julier, Uwe D. Hanebeck,
    % Efficient Bingham Filtering based on Saddlepoint Approximations, 2014
    % Proceedings of the 2014 IEEE International Conference on Multisensor Fusion and Information Integration (MFI 2014), Beijing, China, September 2014.

    f = @(z) findZ(z, kappa);
    Z = fsolve(f, -ones(2,1), optimset('display', 'off', 'algorithm', 'levenberg-marquardt'));
    Z = [Z; 0];
    [Z, ~] = sort(Z,'ascend');

    function R = findZ(Z, rhs)  % needs esternal mex
      Z=[Z;0];
      d = size(Z,1);

      % normalization constant
      A = numericalSaddlepointWithDerivatives(sort(-Z)+1)*exp(1);
      A = A(3);

      % derivative of normalization constant
      B = zeros(1,3);
      dim = size(Z,1);
      for i=1:dim
        mZ = Z([1:i i i:dim]);
        T = numericalSaddlepointWithDerivatives(sort(-mZ)+1)*exp(1)/(2*pi);
        B(i) = T(3);
      end

      R = zeros(d-1,1);
      for i=1:(d-1)
        R(i) = B(i)/A - rhs(i);
      end
    end
  end

end

author: mtex-toolbox

repo: mtex

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


text:
function name = p24_name ( )

%*****************************************************************************80
%
%% P24_NAME returns the name of problem 24.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string NAME, the name of the problem.
%
  name = '(|4X-2|+c)/(1+c)';

  return
end


text:
function [tn, xn,lambda] = calcDiscreteMap(obj, t, x, varargin)
    % calculates the discrete map of the dynamical system that maps
    % xminus from xplus. Subclasses must implement this method by its
    % own.
    %
    % @note By default, this discrete map will be a Identity map.
    % For specific discrete map, please implement as a subclass of
    % this class
    %
    % Parameters:
    % t: the time instant
    % x: the pre-impact states @type colvec
    %
    % Return values:
    % tn: the time after reset @type double
    % xn: the post-impact states @type colvec
    % lambda: the impulsive wrench @type struct

    R = obj.R;
    cstr_name = fieldnames(obj.ImpactConstraints);

    nx = obj.numState;
    qminus = x(1:nx);
    dqminus = x(nx+1:end);

    if isempty(cstr_name)
        qplus = R*qminus;
        dqplus = R*dqminus;
    else
        q = R*qminus;
        dq = R*dqminus;

        %% impact constraints
        cstr = struct2array(obj.ImpactConstraints);
        n_cstr = length(cstr);
        % determine the total dimension of the holonomic constraints
        cdim = sum([cstr.Dimension]);
        % initialize the Jacobian matrix
        Je = zeros(cdim,nx);

        idx = 1;
        for i=1:n_cstr
            c_obj = cstr(i);
            cstr_indices = idx:idx+c_obj.Dimension-1;

            % calculate the Jacobian
            [Jh] = calcJacobian(c_obj,q,dq);
            Je(cstr_indices,:) = Jh;
            idx = idx + c_obj.Dimension;
        end

        % inertia matrix
        De = calcMassMatrix(obj,q);

        % % Compute Dynamically Consistent Contact Null-Space from Lagrange
        % % Multiplier Formulation
        %         DbarInv = Je * (De \ transpose(Je));
        %         I = eye(nx);
        %         Jbar = De \ transpose(transpose(DbarInv) \ Je );
        %         Nbar = I - Jbar * Je;
        % % Apply null-space projection
        %         dqplus = Nbar * dq;

        nImpConstr = size(Je,1);
        A = [De -Je'; Je zeros(nImpConstr)];
        b = [De*dq; zeros(nImpConstr,1)];
        y = A\b;

        ImpF = y((obj.numState+1):end);

        idx = 1;
        for i=1:n_cstr
            c_obj = cstr(i);
            cstr_indices = idx:idx+c_obj.Dimension-1;

            % calculate the Jacobian
            lambda.(c_obj.InputName) = ImpF(cstr_indices);
            idx = idx + c_obj.Dimension;
        end

        dqplus = y(1:obj.numState);

        qplus  = q;

    end

    tn = t;
    xn = [qplus; dqplus];
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [structOut, errors] = checkStructureDefaults(params, defaults)
% Check params input parameter values against defaults put in structOut
errors = cell(0);
fNames = fieldnames(defaults);
structOut = defaults;
for k = 1:length(fNames)
    try
       nextValue = getStructureParameters(params, fNames{k}, ...
                                          defaults.(fNames{k}).value);
       validateattributes(nextValue, defaults.(fNames{k}).classes, ...
                         defaults.(fNames{k}).attributes);
       structOut.(fNames{k}).value = nextValue;
    catch mex
        errors{end + 1} = [fNames{k} ' invalid: ' mex.message]; %#ok<AGROW>
    end
end

author: VisLab

repo: EEG-Clean-Tools


text:
function x = MW2cv(P,N)
%MW1CV   critical Mann-Whitney's U associated to a p-value.
%It obtain a Mann-Whitney's U of two random variables with continuous cumulative
%distribution associated to a p-value. This procedure is highly recommended for sample sizes
%7< nx & ny <=40. For nx & ny <=7 it is recommended to use the MW1cv function;
%otherwise, U-value may be a poor approximation.
%It works with a procedure to get the nearest cumulative distribution relative value to P.
%[Based on the Fortran77 algorithm AS 62 Appl. Statist. (1973)]
%
%   Syntax: function x = MW2cv(P,N)
%
%     Inputs:
%          P - cumulative probability value of interest.
%          N - 2-element vector of sample sizes for the two samples [].
% The input quantities should be scalars.
%     Outputs:
%          x - Mann-Whitney's U statistic.
%
%    Example: For two independent samples we are interested to get the
%             Mann-Whitney's statistic U with an associated cumulative
%             probability P = 0.95. Sample sizes are n1 = 36 and n2 = 14.
%
%                              P = 0.95; N = [36,14];
%
%     Calling on Matlab the function:
%             x = MW2cv(P,N)
%
%       Answer is:
%
%                 328
%

%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.mx
%
%  May 23, 2003.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A. and R. Hernandez-Walls. (2003). MW2cv: Critical Mann-Whitney's U
%    associated to a p-value: nx or ny >7. A MATLAB file. [WWW document]. URL http://
%    www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=3555&objectType=FILE
%
%  References:
%
%  Mann, H. B. and Whitney, D. R. (1947), On a test of whether one of two
%           random variables is stochastically larger than the other. Annals
%           of Mathematical Statistics, 18: 50-60.
%  Algorithm AS 62 (1973). Journal of Applied Statistics, 22(2):1-3.
%

if nargin <  2,
   error('Requires two input arguments.');
end

nmin = min(N);  %largest sample size.
nmax = max(N);  %smallest sample size.

if (nmin <= 7) & (nmax <= 7);
   fprintf('Warning: For nx and ny <= 7, the p-value may be a poor approximation.\n');
   fprintf('It is recommended to use the MW1cv function you can find on the\n');
   fprintf('Matlab>File Exchange Antonio Trujillo-Ortiz'' Author Page.\n');
   disp(' ');
   cont=input('Do you want to continue anyway (y/n):','s');

   if (cont=='y');
      disp('Here it goes.');

      mn1 = prod(N)+1;
      n1 = nmax+1;
      freq = [ones(n1,1); zeros(mn1-n1,1)];

      lwrk = floor((mn1+1)/2 + nmin);
      work = zeros(lwrk,1);

% Generate successively higher-order distributions
      in = nmax;
      for i = 2:nmin
         in = in+nmax;
         n1 = in+2;
         l = 1 + in/2;
         k = i;

% Generate complete distribution from outside inwards
         for j = 1:l
            k = k+1;
            n1 = n1-1;
            summ = freq(j) + work(j);
            freq(j) = summ;
            work(k) = summ - freq(n1);
            freq(n1) = summ;
         end;
      end;

      freq = freq/sum(freq);  % Make distribution relative

% Cumulative frequency distribution
      cumfreq = cumsum(freq);

%Location of the interested Mann-Whitney's U on all the possible U's for this 2-sample sizes.
%Here we are using a procedure to get the nearest fc value to P.
      cumfreq=cumfreq-P;
      u = find(abs(cumfreq)==min(abs(cumfreq(:))));

      UU = [0:length(freq)-1];  %vector of all the possible Mann-Whitney's U values.

%Association of the interested Mann-Whitney's U with its cumulative distribution.
      x = UU(u);
   else
   end
else

   mn1 = prod(N)+1;
   n1 = nmax+1;
   freq = [ones(n1,1); zeros(mn1-n1,1)];

   lwrk = floor((mn1+1)/2 + nmin);
   work = zeros(lwrk,1);

%Generate successively higher-order distributions
   in = nmax;
   for i = 2:nmin
      in = in+nmax;
      n1 = in+2;
      l = 1 + in/2;
      k = i;

%Generate complete distribution from outside inwards
      for j = 1:l
         k = k+1;
         n1 = n1-1;
         summ = freq(j) + work(j);
         freq(j) = summ;
         work(k) = summ - freq(n1);
         freq(n1) = summ;
      end;
   end;

   freq = freq/sum(freq);  % Make distribution relative

% Cumulative frequency distribution
   cumfreq = cumsum(freq);

%Location of the interested Mann-Whitney's U on all the possible U's for this 2-sample sizes.
%Here we are using a procedure to get the nearest fc value to P.
   cumfreq=cumfreq-P;
   u = find(abs(cumfreq)==min(abs(cumfreq(:))));

   UU = [0:length(freq)-1];  %vector of all the possible Mann-Whitney's U values.

%Association of the interested Mann-Whitney's U with its cumulative distribution.
   x = UU(u);
end


text:
% 3D molecular dynamics (GPU vectorized)
% by Brian Haidet for AlphaPhoenix
% published 12/15/2018
% CC non-commercial, attribution

% Don't copy this straight into a homework or something rediculous - I'm
% posting this so people can hopefully learn from it! It's not the cleanest
% code, but I added a few more comments in the areas that need
% explaining. and deleted all the extra bits of commented-out code that had
% built up. Feel free to ask more questions on youtube here!
% https://youtu.be/6DlRsPo-dxY

%print out to images
%close all;
v=VideoWriter(['outMovie' num2str(now()) '.mp4'],'MPEG-4');
v.FrameRate = 30;
open(v);

numFrames = length(dir('c:\MatlabOutput\CsClSimFrames\data'))

breakearly=35404-1
breakearly=3500

tic
iteration = 0;
while true
    if iteration==breakearly
        break;
    end
    iteration = iteration + 1;
    if exist(['c:\MatlabOutput\CsClSimFrames\data\data' num2str(iteration,'%05d') '.mat'],'file')
        load(['c:\MatlabOutput\CsClSimFrames\data\data' num2str(iteration,'%05d') '.mat'])
    else
        break;
    end
    clf;
    figure(1)
    set(gcf, 'Position', [0 0 1920, 1080]); %<- Set size
    pause(.02);
            plotter(iteration);
        drawnow()

%     plotter(iteration);
    writeVideo(v,getframe(gcf));
    disp(['Printing frame ' num2str(iteration) '/' num2str(numFrames) '   |   Est Completion Time: ' datestr(addtodate(now,floor((1000*toc/(iteration)*(numFrames-iteration))),'millisecond'))]);
end
close(v)

author: BrianHaidet

repo: AlphaPhoenix


text:
function c = r8_sign_char ( r8 )

%*****************************************************************************80
%
%% R8_SIGN_CHAR returns the sign of an R8 as a character.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R8, the value whose sign is of interest.
%
%    Output, character C, is '-', '0', or '+'.

  if ( r8 < 0.0 )
    c = '-';
  elseif ( r8 == 0.0 )
    c = '0';
  else
    c = '+';
  end

  return
end


text:
function [pc_vec]=pca_kpm(features,N, method);
% PCA_KPM Compute top N principal components using eigs or svd.
% [pc_vec]=pca_kpm(features,N)
%
% features(:,i) is the i'th example - each COLUMN is an observation
% pc_vec(:,j) is the j'th basis function onto which you should project the data
% using pc_vec' * features

[d ncases] = size(features);
fm=features-repmat(mean(features,2), 1, ncases);

if method==1 % d*d < d*ncases
  fprintf('pca_kpm eigs\n');
  options.disp = 0;
  C = cov(fm'); % d x d matrix
  [pc_vec, evals] = eigs(C, N, 'LM', options);
else
  % [U,D,V] = SVD(fm), U(:,i)=evec of fm fm', V(:,i) = evec of fm' fm
  fprintf('pca_kpm svds\n');
  [U,D,V] = svds(fm', N);
  pc_vec = V;
end

if 0
X = randn(5,3);
X = X-repmat(mean(X),5,1);
C = X'*X;
C2=cov(X)
[U,D,V]=svd(X);
[V2,D2]=eig(X)
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% This file create a simple regression based network for speech
% dereverberation or enhancement
%
% Created by Xiong Xiao, Temasek Laboratories, NTU, Singapore.
% Last Modified: 08 Feb 2017
%
function layer = genNetworkSeparation_MaskingMagnitude(para)
para.freqBin = (0:1/para.fft_len:0.5)*2*pi; % w = 2*pi*f, where f is the normalized frequency k/N is from 0 to 0.5.
nFreqBin = length(para.freqBin);

% Part 1: generate the BF weight predicting subnet

if para.useWav  % input is row waveform
    layer{1}.name = 'Input';
    layer{end}.inputIdx = 1;
    layer{end}.dim = [1 1]*para.nCh;

    layer{end+1}.name = 'stft';
    layer{end}.prev = -length(layer)+1;
    layer{end}.fft_len = para.fft_len;
    layer{end}.frame_len = para.frame_len;
    layer{end}.frame_shift = para.frame_shift;
    layer{end}.removeDC = para.removeDC;
    layer{end}.win_type = para.win_type;
    layer{end}.dim = [(para.fft_len/2+1)*para.nCh layer{length(layer)+layer{end}.prev}.dim(1)];
    layer{end}.skipBP = 1;  % skip backpropagation

    layer{end+1}.name = 'Affine';       % scaling the Fourier transform
    layer{end}.prev = -1;
    layer{end}.W = [];
    layer{end}.b = [];
    layer{end}.dim = [1 1] * layer{length(layer)+layer{end}.prev}.dim(1);
    layer{end}.update = 0;
    layer{end}.skipBP = 1;  % skip backpropagation

    % extract a subset of dimensions for prediction. Not used currently, so we
    % will use all the dimensions.
    layer{end+1}.name = 'ExtractDims';
    layer{end}.prev = -1;
    layer{end}.dimIndex = 1:nFreqBin*para.nCh;
    layer{end}.dim = [length(layer{end}.dimIndex) layer{length(layer)+layer{end}.prev}.dim(1)];
    layer{end}.skipBP = 1;

    % get the log power spectrum and perform CMN
    layer{end+1}.name = 'Power';
    layer{end}.prev = -1;
    layer{end}.dim = [1 1] * layer{end-1}.dim(1);
    layer{end}.skipBP = 1;

    layer{end+1}.name = 'log';
    layer{end}.const = 0.01;
    layer{end}.prev = -1;
    layer{end}.dim = [1 1]*layer{end-1}.dim(1);
    layer{end}.skipBP = 1;

else    % input is log spectrogram
    layer{1}.name = 'Input';
    layer{end}.inputIdx = 1;
    layer{end}.dim = [1 1]*(para.fft_len/2+1)*para.nCh;
end

outputDim = nFreqBin;

switch para.NetType
    case 'DNN'
        layer{end+1}.name = 'Splice';
        layer{end}.prev = -1;
        layer{end}.context = para.contextSize;
        layer{end}.dim = [layer{end}.context 1]*layer{length(layer)+layer{end}.prev}.dim(1);
        layer{end}.skipBP = 1;

        layer{end+1}.name = 'Affine';       % scaling the Fourier transform
        layer{end}.prev = -1;
        layer{end}.W = [];
        layer{end}.b = [];
        layer{end}.dim = [1 1] * layer{length(layer)+layer{end}.prev}.dim(1);
        layer{end}.update = 0;

        layerRegression = genNetworkFeedForward_v2(layer{end}.dim(1), para.hiddenLayerSize, outputDim, 'mse', 'sigmoid');
    case 'LSTM'
        layer{end+1}.name = 'Affine';       % scaling the Fourier transform
        layer{end}.prev = -1;
        layer{end}.W = [];
        layer{end}.b = [];
        layer{end}.dim = [1 1] * layer{length(layer)+layer{end}.prev}.dim(1);
        layer{end}.update = 0;

        tmpTopology.inputDim = layer{end}.dim(1);
        tmpTopology.hiddenLayerSizeLSTM = para.hiddenLayerSize;
        tmpTopology.usePastState = zeros(1,length(para.hiddenLayerSize)); % do not use peeping hole
        tmpTopology.hiddenLayerSizeFF = [];
        tmpTopology.outputDim = outputDim;
        tmpTopology.costFn = 'mse';
        tmpTopology.LastActivation4MSE = 'sigmoid';
        layerRegression = genNetworkLSTM(tmpTopology);
end
layer = [layer layerRegression(2:end-2)];

if ~isempty(para.hiddenLayerSizeFF)
    layerFF = genNetworkFeedForward_v2(para.hiddenLayerSize(end), para.hiddenLayerSizeFF, outputDim, 'mse', 'sigmoid');
    layer = [layer(1:end-1) layerFF(2:end-2)];
end

layer = [layer layer(end-1:end)];
layer{end-1}.prev = -3;
mask_idx = [-2 0] + length(layer);

% apply the masking
if para.useWav  % input is row waveform
    power_idx = ReturnLayerIdxByName(layer, 'power');

    layer{end+1}.name = 'sqrt';
    layer{end}.prev = power_idx-length(layer);
    layer{end}.dim = [1 1]*layer{power_idx}.dim(1);
    layer{end}.skipBP= 1;

    mag_idx = length(layer);
else
    input_idx = ReturnLayerIdxByName(layer, 'input');
    layer{end+1}.name = 'exp';
    layer{end}.prev = input_idx(1)-length(layer);
    layer{end}.dim = [1 1]*layer{input_idx(1)}.dim(1);
    power_idx = length(layer);
end

MaskingLayer{1}.name = 'hadamard';
MaskingLayer{end}.dim = [1 1]*layer{end}.dim(1);
MaskingLayer{end}.prev = [mask_idx(1) mag_idx(1)] - length(layer) -1;

layer = [layer MaskingLayer];
output_idx(1) = length(layer);

MaskingLayer2 = MaskingLayer;
MaskingLayer2{1}.prev = [mask_idx(2) mag_idx(1)] - length(layer) -1;

layer = [layer MaskingLayer2];
output_idx(2) = length(layer);

if para.useWav
    power_idx = ReturnLayerIdxByName(layer, 'power');
    CleanLayer = layer(1:power_idx(1));
    CleanLayer{1}.inputIdx = 2;
    CleanLayer{1}.dim(:) = 1;
    CleanLayer{2}.dim(:) = CleanLayer{2}.dim(:)/para.nCh;
    CleanLayer{3}.dim(:) = CleanLayer{3}.dim(:)/para.nCh;
    CleanLayer(4) = [];
else
    CleanLayer{1}.name = 'Input';
    CleanLayer{end}.inputIdx = 2;
    CleanLayer{end}.dim = [1 1]*(para.fft_len/2+1);
end
CleanLayer{end+1}.name = 'sqrt';
CleanLayer{end}.prev = -1;
CleanLayer{end}.dim = [1 1]*CleanLayer{end-1}.dim(1);

for i=1:length(CleanLayer)
    CleanLayer{i}.skipBP = 1;
end

layer = [layer CleanLayer];
clean1_idx = length(layer);
CleanLayer{1}.inputIdx = 3;
layer = [layer CleanLayer];
clean2_idx = length(layer);

layer{end+1}.name = 'mixture_MSE';
layer{end}.prev = [output_idx clean1_idx clean2_idx] - length(layer);
layer{end}.dim = [1 layer{length(layer)+layer{end}.prev(1)}.dim(1)];

layer = FinishLayer(layer);
end

author: singaxiong

repo: SignalGraph


text:
function WriteMhaFile(filename, img_size, resolution, data_type, offset)
% WRITEMHAFILE  Write the header part of a MetaImage file (.mha)
%
% WRITEMHAFILE(FILENAME, IMG_SIZE, RESOLUTION, DATA_TYPE, OFFSET)
%
%   FILENAME is the path and name of the file to be written, e.g.
%   'foo.mha'.
%
%   IMG_SIZE is a 3-vector with the size of the output volume.
%
%   RESOLUTION is a 3-vector with the voxel size in the 3 directions.
%
%   DATA_TYPE is a string with the data type as given my Matlab, e.g.
%   'uint8', 'short', 'uint16'.
%
%   OFFSET is a vector with the real world coordinates in metres (not index
%   coordinates) of the first voxel in the volume. For example,
%   OFFSET=[0.014552, 0.010486, 0.00142]. By default, OFFSET=[0 0 0].
%
%   See also: WriteRawFile, writemetaimagefile.

% Author(s): Ramon Casero <rcasero@gmail.com> and Vicente Grau
% Copyright © 2009-2012 University of Oxford
% Version: 0.1.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error(nargchk(4, 5, nargin, 'struct'));
error(nargoutchk(0, 0, nargout, 'struct'));

% defaults
if (nargin < 5 || isempty(offset))
    offset = [0.0 0.0 0.0];
end

[path, name] = fileparts(filename);

% open file for writing
fid=fopen(filename, 'w');
if(fid<=0)
    fprintf('Impossible to open file %s\n', filename);
end

ndims=numel(resolution);

if(ndims == 3)
    fprintf(fid, 'NDims = 3\n');

    fprintf(fid, 'DimSize = %d %d %d\n', img_size(2), img_size(1), img_size(3));

    if(strcmp(data_type, 'uint8'))
        fprintf(fid, 'ElementType = MET_UCHAR\n');
    elseif(strcmp(data_type, 'short'))
        fprintf(fid, 'ElementType = MET_SHORT\n');
    elseif (strcmp(data_type, 'uint16'))
        fprintf(fid, 'ElementType = MET_USHORT\n');
    else
        error('Not implemented data type')
    end

    fprintf(fid, 'Offset = %1.6f %1.6f %1.6f\n', ...
        offset(1), offset(2), offset(3));

    fprintf(fid, 'ElementSpacing = %1.12f %1.12f %1.12f\n', resolution(1), resolution(2), resolution(3));

elseif(ndims==4)
    fprintf(fid, 'NDims = 4\n');

    fprintf(fid, 'DimSize = %d %d %d %d\n', img_size(2), img_size(1), img_size(3), img_size(4));

    if(strcmp(data_type, 'uint8'))
        fprintf(fid, 'ElementType = MET_UCHAR\n');
    elseif(strcmp(data_type, 'short'))
        fprintf(fid, 'ElementType = MET_SHORT\n');
    elseif (strcmp(data_type, 'uint16'))
        fprintf(fid, 'ElementType = MET_USHORT\n');
    else
        error('Not implemented data type')
    end

    fprintf(fid, 'Offset = %1.6f %1.6f %1.6f\n', ...
        offset(1), offset(2), offset(3));

    fprintf(fid, 'ElementSpacing = %1.12f %1.12f %1.12f %1.12f\n', resolution(1), resolution(2), resolution(3), resolution(4));

end

fprintf(fid, 'ElementByteOrderMSB = False\n');

fprintf(fid, 'ElementDataFile = %s\n', [name, '.raw']);

fclose(fid);


%**************************************************************************
% SARSA Learning applied to Cart-Pole balancing problem.
% The environment of the learning system is a black box, from which it has
% several lines and a reinforcement line. Its task is to learn to give
% responses which maximize the scalar signals on its reinforcement line.
% The Q-value update is the following, if the system takes action a from
%   state s at time t, and arrives at state ss with feedback r at time t+1:
%
%   Q(t+1, s, a) = Q(t, s, a)
%                  + alpha (r + gamma *Q(t,ss, b) - Q(t, s, a))

% get_box:  Given the current state, returns a number from 1 to 162
%           designating the region of the state space encompassing the current
%           state.
%           Returns a value of -1 if a failure state is encountered.

% cart_pole2: The cart and pole dynamics; given the Force and
%             current state, estimates next state

%Code written by: Savinay Nagendra
%email id:        nagsavi17@gmail.com
%**************************************************************************
clc;
clear all;
close all;
% Initialization

NUM_BOXES = 163;
ALPHA = 0.5;             % Learning rate parameter
GAMMA = 0.999;           % Discount factor for future reinf
Q = zeros(NUM_BOXES,2);  % State-Action Values
action = [10 -10];
MAX_FAILURES = 2000;
MAX_STEPS = 150000;
epsilon = 0;
steps = 0;
failures = 0;
thetaPlot = 0;
xPlot = 0;
%Pendulum state initialization
theta = 0;
thetaDot = 0;
x = 0;
xDot = 0;
box = getBox4(theta,thetaDot,x,xDot);
if(rand>epsilon)       % exploit
    [~,actionMax] = max(Q(box,:));
    currentAction = action(actionMax);
else                   % explore
    currentAction = datasample(action,1);
end
actionIndex1 = find(action == currentAction);

while(steps<=MAX_STEPS && failures<+MAX_FAILURES)
    steps = steps + 1;
    [thetaNext,thetaDotNext,thetaacc,xNext,xDotNext] = cart_pole2(currentAction,theta,thetaDot,x,xDot);
    thetaPlot(end + 1) = thetaNext;
    xPlot(end + 1) = xNext;
    newBox = getBox4(thetaNext,thetaDotNext,xNext,xDotNext);
    theta = thetaNext;
    thetaDot = thetaDotNext;
    x = xNext;
    xDot = xDotNext;
    if(newBox == 163)
        r = -1;
        Q(newBox,:) = 0;
        figure(2);
        plot((1:length(thetaPlot)),thetaPlot,'-ob');
        figure(3);
        plot((1:length(xPlot)),xPlot,'-og');
        thetaPlot = 0;
        xPlot = 0;
        %Swing Up. Find the box.
        theta = 0;
        thetaDot = 0;
        x = 0;
        xDot = 0;

        newBox = getBox4(theta,thetaDot,x,xDot);
        failures = failures + 1;
        fprintf('Trial %d was %d steps. \n',failures,steps);
        figure(1);
        plot(failures,steps,'-or');
        hold on;
        steps = 0;
    else
        r = 0;
    end
    if(rand > epsilon)  %exploit
        [~,newActionMax] = max(Q(newBox,:));
        newAction = action(newActionMax);
    else                %explore
        newAction = datasample(action,1);
    end
    actionIndex2 = find(action == newAction);
    Q(box,actionIndex1) = Q(box,actionIndex1) + ALPHA*(r + GAMMA*Q(newBox,actionIndex2) - Q(box,actionIndex1));
    box = newBox;
    currentAction = newAction;
    actionIndex1 = actionIndex2;
end

if(failures == MAX_FAILURES)
    fprintf('Pole not balanced. Stopping after %d failures.',failures);
else
    fprintf('Pole balanced successfully for at least %d steps\n', steps);
    figure(1);
    plot(failures+1,steps,'-or');
    hold on;
    figure(2);
    plot((1:length(thetaPlot)),thetaPlot,'-ob');
    figure(3);
    plot((1:length(xPlot)),xPlot,'-og');
    figure(4);
    plot((1:301),thetaPlot(1:301),'-ob');
    hold on;
    figure(5);
    plot((1:301),xPlot(1:301),'-og');
    hold on;
end


text:
function X = kr(U,varargin)
%KR Khatri-Rao product.
%   kr(A,B) returns the Khatri-Rao product of two matrices A and B, of
%   dimensions I-by-K and J-by-K respectively. The result is an I*J-by-K
%   matrix formed by the matching columnwise Kronecker products, i.e.,
%   the k-th column of the Khatri-Rao product is defined as
%   kron(A(:,k),B(:,k)).
%
%   kr(A,B,C,...) and kr({A B C ...}) compute a string of Khatri-Rao
%   products A x B x C x ..., where x denotes the Khatri-Rao product.
%
%   See also kron.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)

if ~iscell(U), U = [{U} varargin]; end
[J,K] = size(U{end});
if any(cellfun('size',U,2) ~= K)
    error('kr:U','Input matrices should have the same number of columns.');
end

X = reshape(U{end},[J 1 K]);
for n = length(U)-1:-1:1
    I = size(U{n},1);
    A = reshape(U{n},[1 I K]);
    X = reshape(bsxfun(@times,A,X),[I*J 1 K]);
    J = I*J;
end
X = reshape(X,[size(X,1) K]);

author: andrewssobral

repo: lrslibrary


text:
function [verticesList, facesList] = getVerticesAndFacesInSphere(mesh_outer, iV, radius)

% find all the vertices included in a sphere
verticesList = [];

for vertex = 1:length(mesh_outer.vertices)
    if isVertexInRadius(mesh_outer.vertices(vertex,:), mesh_outer.vertices(iV,:), radius)
        verticesList = [verticesList,vertex];
    end
end

% find faces to which those vertices belong
facesList = [];
for vert = verticesList
    facesList = [facesList, mesh_outer.facesOfVertex(vert).faceList];
end

facesList = unique(facesList);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft