text:
function [cij,flag] = makerandCIJdegreesfixed(in,out)
%MAKERANDCIJDEGREESFIXED        Synthetic directed random network
%
%   CIJ = makerandCIJdegreesfixed(N,K);
%
%   This function generates a directed random network with a specified
%   in-degree and out-degree sequence. The function returns a flag,
%   denoting whether the algorithm succeeded or failed.
%
%   Inputs:     in,     indegree vector
%               out,    outdegree vector
%
%   Output:     CIJ,    binary directed connectivity matrix
%               flag,   flag=1 if the algorithm succeeded; flag=0 otherwise
%
%
%   Notes:  Necessary conditions include:
%               length(in) = length(out) = n
%               sum(in) = sum(out) = k
%               in(i), out(i) < n-1
%               in(i) + out(j) < n+2
%               in(i) + out(i) < n
%
%           No connections are placed on the main diagonal
%
%
% Aviad Rubinstein, Indiana University 2005/2007

% intialize
n = length(in);
k = sum(in);
inInv = zeros(k,1);
outInv = inInv;
iIn = 1; iOut = 1;

for i = 1:n
    inInv(iIn:iIn+in(i) - 1) = i;
    outInv(iOut:iOut+out(i) - 1) = i;
    iIn = iIn+in(i);
    iOut = iOut+out(i);
end

cij = eye(n);
edges = [outInv(1:k)'; inInv(randperm(k))'];

% create cij, and check for double edges and self-connections
for i = 1:k
    if cij(edges(1,i),edges(2,i)),
        warningCounter = 1;
        while (1)
            switchTo = ceil(k*rand);
            if ~(cij(edges(1,i),edges(2,switchTo)) || cij(edges(1,switchTo),edges(2,i))),
                cij(edges(1,i),edges(2,switchTo)) = 1;
                if switchTo < i,
                    cij(edges(1,switchTo),edges(2,switchTo)) = 0;
                    cij(edges(1,switchTo),edges(2,i)) = 1;
                end
                temp = edges(2,i);
                edges(2,i) = edges(2,switchTo);
                edges(2,switchTo) = temp;
                break
            end
            warningCounter = warningCounter+1;
            % If there is a legitimate subtitution, it has a probability of 1/k of being done.
            % Thus it is highly unlikely that it will not be done after 2*k^2 attempts.
            % This is an indication that the given indegree / outdegree
            % vectors may not be possible.
            if warningCounter == 2*k^2
                flag = 0;  % no valid solution found
                return;
            end
        end
    else
        cij(edges(1,i),edges(2,i)) = 1;
    end
end

cij = cij - eye(n);

% a valid solution was found
flag = 1;

author: fieldtrip

repo: fieldtrip


text:
function pass = test_plotting( pref )
% Check that the very basic plotting commands do not crash in diskfunv

F = grad(diskfun(@(x,y) x.*cos(y)));
try
    hold off
    quiver(F),         j = ishold;
    close all
    if ( j == 0 )
        pass(1) = 1;
    else
        pass(1) = 0;
    end
catch
    pass(1) = 0;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [region_feats, regions] = convFeat_to_poolFeat_multi_region(...
    multiple_region_params, conv_feats, boxes, random_scale)
% convFeat_to_poolFeat_multi_region given the convolutional features of an
% image, it adaptively pools fixed size region features for each bounding
% box and for multiple type of regions.
%
% INPUTS:
% 1) multiple_region_params: is a M x 1 vector of objects of type struct that
% specify the region pooling parameters for each of the M types of regions
% 2) conv_feats: is a K x 1 cell vector or vector of objects, with K >=1
% different type of convolutional features.
% 3) boxes: a N x 4 array with the bounding box coordinates in the form of
% [x0,y0,x1,y1] (where (x0,y0) is the top-left corner and (x1,y1) the
% bottom left corner)
% 4) random_scale: a boolean value that if set to true then each bounding
% box is projected to the convolutional features of a random scale of the
% image.
%
% OUTPUTS:
% 1) region_feats: is a M x 1 cell array where region_feats{i} is a N x F_i
% array with the region features of each of the N bounding boxes for the
% i-th type of region. F_i is the number of features of the i-th type of
% region.
% 2) regions: is a M x 1 cell array with the region coordinates of each
% bounding box and for each type of region. Specifically, regions{i} is a
% N x 8 array that contains the region coordinates of each of the N bounding
% boxes for the i-th type of region. Note that each region is represented
% by 8 values [xo0, yo0, xo1, yo1, xi0, yi0, xi1, yi1] that correspond to
% its outer rectangle [xo0, yo0, xo1, yo1] and its inner rectangle
% [xi0, yi0, xi1, yi1].
%
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
% Part of the code in this file comes from the SPP-Net code:
% https://github.com/ShaoqingRen/SPP_net
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

if ~exist('random_scale', 'var'), random_scale = false; end

num_regions  = length(multiple_region_params);
region_feats = cell(num_regions,1);
regions      = cell(num_regions,1);

for p = 1:num_regions, region_feats{p} = single([]); end
for p = 1:num_regions, regions{p}   = single([]); end

if isempty(boxes), return; end

% pool the region features the bounding boxes for each type of region
if length(conv_feats) == 1
    for p = 1:num_regions
        [region_feats{p}, regions{p}] = convFeat_to_poolFeat(...
            multiple_region_params(p), conv_feats, boxes, random_scale);
    end
else
    if iscell(conv_feats)
        for p = 1:num_regions
            [region_feats{p}, regions{p}]  = convFeat_to_poolFeat(...
                multiple_region_params(p), conv_feats{multiple_region_params(p).feat_id}, ...
                boxes, random_scale);
        end
    elseif isstruct(conv_feats)
        for p = 1:num_regions
            [region_feats{p}, regions{p}]  = convFeat_to_poolFeat(...
                multiple_region_params(p), conv_feats(multiple_region_params(p).feat_id), ...
                boxes, random_scale);
        end
    end
end
end

author: gidariss

repo: mrcnn-object-detection


text:
function g = restrict(f, s)
%RESTRICT   Restrict a BNDFUN to a subinterval.
%   RESTRICT(F, S) returns a BNDFUN that is restricted to the subinterval
%   [S(1), S(2)] of F.domain.
%
%   If LENGTH(S) > 2, i.e., S = [S1, S2, S3, ...], then RESTRICT(F, S) returns
%   a cell-array of BNDFUN objects, where the cells contain F restricted to each
%   of the subintervals defined by S.

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

% Deal with empty case:
if ( isempty(f) )
    g = f;
    return
end

% Check if s is actually a subinterval:
hs = diff(f.domain)*eps; % Some wiggle room.
if ( s(1) < f.domain(1) )
    if ( abs(s(1) - f.domain(1)) < hs )
        s(1) = f.domain(1);
    else
        error('CHEBFUN:BNDFUN:restrict:badInterval', 'Not a valid interval.')
    end
end
if ( s(end) > f.domain(2) )
    if ( abs(s(end) - f.domain(2)) < hs )
        s(end) = f.domain(2);
    else
        error('CHEBFUN:BNDFUN:restrict:badInterval', 'Not a valid interval.')
    end
end
if ( any(diff(s) <= 0) )
    error('CHEBFUN:BNDFUN:restrict:badInterval', 'Not a valid interval.')
end

if ( (numel(s) == 2) && all( s == f.domain) )
    % Nothing to do here!
    g = f;
    return
end

% Compute the breaks in [-1,1] space and restrict the onefun:
t = f.mapping.Inv(s);

% Restrict the ONEFUN field of f.
restrictedOnefuns = restrict(f.onefun, t);

% In case S = [S1, S2], i.e., we are only restricting to one subinterval,
% restrictedOnefuns will be a ONEFUN. In case S = [S1, S2, ..., SN], i.e., we
% are restricting to multiple subinterval, restrictedOnefuns will be a
% cell-array of ONEFUN objects. We need to treat the cases differently, as we
% also want to either return a BNDFUN or cell-array of BNDFUN objects depending
% on how many subintervals we are restricting to.
if ( length(s) == 2 )
    % Only restricting to one subinterval -- return a BNDFUN.

    % Create an empty BNDFUN, and assign fields directly. This is faster than
    % using the BNDFUN constructor.
    g = bndfun();
    g.onefun = restrictedOnefuns;
    g.domain = s;
    g.mapping = bndfun.createMap(s);
else
    % Restricting to multiple subintervals -- return a cell-array of BNDFUN
    % objects.

    % Create a cell to be returned.
    g = cell(1, numel(s) - 1);

    % Create an empty BNDFUN:
    emptyBndfun = bndfun();

    % Loop over each of the new subintervals, make a bndfun with new mapping,
    % and store in the cell returned:
    for k = 1:(numel(s) - 1)
        % Assign fields directly to an empty temporary BNDFUN. This is faster
        % than using the BNDFUN constructor.
        gTemp = emptyBndfun;
        gTemp.onefun = restrictedOnefuns{k};
        gTemp.domain = s(k:k+1);
        gTemp.mapping = bndfun.createMap(s(k:k+1));
        g{k} = gTemp;
    end
end

end


text:
% Wrapper for vl_nnselection block
% inputs{1} :   X     : 1 x 2 x N x b
% obj       :   idx   : 1 x n
% outputs{1}:   y     : 1 x 2 x n x b

classdef selection < dagnn.Layer

    properties
        idx = [];
    end

    methods
        function outputs = forward(obj, inputs, ~)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
                outputs{1} = gpuArray( vl_nnselection(gather(inputs{1}), obj.idx) );
            else
                outputs{1} = vl_nnselection(inputs{1}, obj.idx);
            end

        end

        function [derInputs, derParams] = backward(obj, inputs, ~, derOutputs)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
                derInputs{1} = gpuArray( vl_nnselection(gather(inputs{1}), obj.idx, gather(derOutputs{1})) );
            else
                derInputs{1} = vl_nnselection(inputs{1}, obj.idx, derOutputs{1});
            end

            derParams = {};
        end

        function outputSizes = getOutputSizes(obj, inputSizes)
            outputSizes = {1 inputSizes{1}(2) length(obj.idx) inputSizes{1}(4)} ;
        end

        function obj = selection(varargin)
            obj.load(varargin);
        end
    end
end

author: anilbas


text:
function GAconstrain
% Genetic Algorithm(real coding)
% By: Javad Ivakpour
% E-mail: javad7@gmail.com
% May 2006
% Goal: find maximum of function that introduced in fun00.m file in current
% directory and can be plotted in plot00
% this file is also include the random serach for comparision
tic
clc
figure(1)
clf
clear all
format long

%------------------------        parameters        ------------------------
% befor using this function you must specified your function in fun00.m
% file in current directory and then set the parameters
var=2;            % Number of variables (this item must be equal to the
                  %   number of variables that is used in the function in
                  %   fun00.m file)
n=100;            % Number of population

m0=30;            % Number of generations that max value remains constant
                  %   (use for termination criteria)
nmutationG=20;                  %number of mutation children(Gaussian)
nmutationR=20;                  %number of mutation children(random)
nelit=2;                        %number of elitism children
valuemin=ones(1,var)*-5*pi;     % min possible value of variables
valuemax=ones(1,var)*5*pi;      % max possible value of variables
A=[1 1;-0.5 +1;1 -1;-1 -1];     % constrains: Ax+B>0
B=[5*pi 0 5*pi +5*pi+10];
%-------------------------------------------------------------------------
[r,c]=size(A);
if c~=var
    disp('Number of columns in A must be equivalent to var')
    A
    var
    return
end
nmutation=nmutationG+nmutationR;
sigma=(valuemax-valuemin)/10;    %Parameter that related to Gaussian
                                 %   function and used in mutation step
max1=zeros(nelit,var);
parent=zeros(n,var);
cu=[valuemin(1) valuemax(1) valuemin(2) valuemax(2)];
for l=1:var
    p(:,l)=valuemin(l)+rand(n,1).*(valuemax(l)-valuemin(l));
end
initial=p;
m=m0;
maxvalue=ones(m,1)*-1e10;
maxvalue(m)=-1e5;
g=0;
meanvalue(m)=0;
%-------------   ****    termination criteria   ****-------------
while abs(maxvalue(m)-maxvalue(m-(m0-1)))>0.001*maxvalue(m) &...
        (abs(maxvalue(m))>1e-10 & abs(maxvalue(m-(m0-1)))>1e-10)...
        & m<10000 & abs(maxvalue(m)-meanvalue(m))>1e-5 | m<20
    sigma=sigma./(1.05);% reducing the sigma value
    %  ------  **** % reducing the number of mutation()random   **** ----
    g=g+1;
    if g>10 & nmutationR>0
        g=0;
        nmutationR=nmutationR-1;
        nmutation=nmutationG+nmutationR;
    end

    %-------------   ****    function evaluation   ****-------------
    for i=1:n
        y(i)=fun00(p(i,:));
        flag(i)=0;
    end
    if var==2
        figure(1)
        hold off
        plot00(cu)
        hold on
        plot3(p(:,1),p(:,2),y,'ro')
        title({'Genetic Algorithm: constrained problem'...
            ,'Performance of GA ( o : each individual)',...
            'blue o:possible solution & red o impossible solution ',...
            'black o: best solution'},'color','b')
    end
    s=sort(y);
    maxy=max(y);
    miny=min(y);
    for i=1:n
        k=length(B);
        for j=1:k
            x=p(i,:);
            if (A(j,:)*x'+B(j))<=0
                y(i)=miny;
                flag(i)=1;
            end
        end
    end
    if var==2
        for i=1:n
            if flag(i)==0
                plot3(p(i,1),p(i,2),y(i),'bo')
            end
        end
    end
    s=sort(y);
    maxvalue1(1:nelit)=s(n:-1:n-nelit+1);
    if nelit==0
        maxvalue1(1)=s(n);
        for i=1:n
            if y(i)==maxvalue1(1)
                max1(1,:)=p(i,:);
            end
        end
    end
    for k=1:nelit
        for i=1:n
            if y(i)==maxvalue1(k)
                max1(k,:)=p(i,:);
            end
        end
    end
    if var==2
        hold on
        plot3(max1(1,1),max1(1,2),maxvalue1(1),'kh')
    end
    y=y-min(y)*1.02;
    sumd=y./sum(y);
    meanvalue=y./(sum(y)/n);

    %-------------   ****   Selection: Rolette wheel   ****-------------
    for l=1:n
        sel=rand;
        sumds=0;
        j=1;
        while sumds<sel
            sumds=sumds+sumd(j);
            j=j+1;
        end
        parent(l,:)=p(j-1,:);
    end
    p=zeros(n,var);

    %-------------   ****    regeneration   ****-------------
    for l=1:var

        %-------------   ****    cross-over   ****-------------
        for j=1:ceil((n-nmutation-nelit)/2)
            t=rand*1.5-0.25;
            p(2*j-1,l)=t*parent(2*j-1,l)+(1-t)*parent(2*j,l);
            p(2*j,l)=t*parent(2*j,l)+(1-t)*parent(2*j-1,l);
        end

        %-------------   ****    elitism   ****-------------
        for k=1:nelit
            p((n-nmutation-k+1),l)=max1(k,l);
        end

        %-------------   ****    mutation   ****-------------
        for i=n-nmutation+1:n-nmutationR
            phi=1-2*rand;
            z=erfinv(phi)*(2^0.5);
            p(i,l)=z*sigma(l)+parent(i,l);

        end
        for i=n-nmutationR+1:n
            p(i,1:var)=valuemin(1:var)+rand(1,var).*(valuemax(1:var)...
                -valuemin(1:var));
        end
        for i=1:n
            for l=1:var
                if p(i,l)<valuemin(l)
                    p(i,l)=valuemin(l);
                elseif p(i,l)>valuemax(l)
                    p(i,l)=valuemax(l);
                end
            end
        end
    end
    p;
    m=m+1;
    max1;
    maxvalue(m)=maxvalue1(1);
    maxvalue00(m-m0)=maxvalue1(1);
    mean00(m-m0)=sum(s)/n;
    meanvalue(m)=mean00(m-m0);
    figure(1)
    if var~=2
        hold off
        plot(maxvalue00)
        hold on
        plot(mean00,'g')
        hold on
        title({'Performance of GA',...
            'best value GA:blue, mean value GA:green',''}...
            ,'color','b')
        xlabel('number of generations')
        ylabel('value')
    end
    pause(0.001)
end

clc
disp('     Genetic Algorithm(real coding)   ')
disp('          By: Javad Ivakpour          ')
disp('       E-mail: javad7@gmail.com       ')
disp('**************************************')
num_of_fun_evaluation=n*m
max_point_GA=max1(1,:)
maxvalue_GA=maxvalue00(m-m0)
if var==2
    figure(1)
    hold on
    plot3(max1(1,1),max1(1,2),maxvalue1,'yp')
    hold on
end
figure(2)
title('Performance of GA(best value)','color','b')
xlabel('number of generations')
ylabel('max value of best solution')
hold on
plot(maxvalue00)
hold on
toc

author: Sable

repo: mcbench-benchmarks


text:
function script_labeling_KAIST_dataset( pth, setNum, vidNum )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% KAIST Multi-spectral Recognition Dataset
%                                Yukyung Choi et al. 2016.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% When the 'pngreadc.mexw64' makes a problem,
%   do followings:
%       >> clear mex; d = fileparts(which(mfilename));
%       >> delete( fullfile(d, videos, private, *png*) );
%
% - {rootPath} for KAIST dataset
%   - Set00 (AM05)
%       - V000 (North)
%           - T8  (Splitter)
%           - RGB (Splitter)
%       - V001 (West)
%       - V002 (East)
%   - Set01
%   - Set02
%   - Set03
%   - Set04
%   - Set05
%   - Set06
%   - Set07

d = fileparts(which(mfilename));
addpath( genpath( d ) );

if nargin < 3, return; end
if ~exist( pth, 'dir' ), pth = '../data/kaist-rgbt-ped/videos'; end

name = sprintf( '%s/Set%02d/V%03d/', pth, setNum, vidNum );

targetObjs = { 'person', 'car', 'people', 'ignore' };
occTypes = {'No-occ', 'Partial-occ', 'Heavy-occ'};
seqNames = { [ name 'RGB' ], [ name 'T8' ] };

hFig = vbbLabeler2( targetObjs, occTypes, seqNames );
uiwait( hFig );

tmpFiles = dir( 'tmp*' );
arrayfun( @(x) delete( x.name ), tmpFiles, 'uniformoutput', false );

author: SoonminHwang

repo: rgbt-ped-detection


text:
%[x,y]=find(fixations);
%densityMap= make_gauss_masks(y,x,[video_res_y,video_res_x]);

function my_mask= make_gauss_masks(x,y,my_size,fixTime)
%This function drop Gaussians centered at (x,y) locations. my_size is the
%size of your original image (also the size of the output of this function).

% Author    : Umesh Rajashekar
% Date      : 4 Jan. 2008

if (nargin < 4)
    fixTime = ones(size(x));
end

W = 1; %Full width at half max in visual degrees of angle for the Gaussian
W = W*30; %For my experiments, 1 degree of visual angle is 60 pixels

%Convert full width at half max to variance
my_sigma = W/(2*sqrt(2*log(2)));

R = my_size(1);C= my_size(2);
big_R=2*R+1;big_C= 2*C+1; %Double your gauss size

[X Y] = meshgrid(0:big_C-1,0:big_R-1);
my_mask = zeros(R,C);

temp_num = (X-C).^2 + (Y-R).^2;
temp = exp(-temp_num/(2*my_sigma^2));
big_gauss = temp /(2*pi*my_sigma);

for i = 1: length(x)
    temp_x = floor(x(i));temp_y = floor(y(i));
    temp = big_gauss(R-temp_y:R-temp_y+R-1 ,C-temp_x:C-temp_x+C-1);
    my_mask = my_mask + temp*fixTime(i);   %Weight by time of fixation
end

my_mask = my_mask./max(my_mask(:));

author: wenguanwang


text:
%kkimportasc 'Import ASCII Data File to Data Segment'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kimportasc.pane file
%
% Parameters:
% InputFile: i1 'ASCII Input File', required: 'ASCII input data file containing new data segment'
% InputFile: i2 'Data Object Input', optional: 'data object input to insert the new segment into'
% Integer: so 'Start Offset ', default: 0: 'starting point in ASCII input file to start importing at'
% Integer: sf 'Skip Factor ', default: 0: 'number of points in ASCII file to skip between reads'
% Integer: nr 'Number of Reads ', default: 1: 'number of points in ASCII file to read between skips'
% Integer: dim 'Segment Dimension ', default: 5: 'dimension of the segment'
% OutputFile: o 'Output Object', required: 'output data object'
% Integer: wsize 'WIDTH ', default: 1: 'size of width dimension in index specified'
% Integer: hsize 'HEIGHT ', default: 1: 'size of height dimension in index specified'
% Integer: dsize 'DEPTH ', default: 1: 'size of depth dimension in index specified'
% Integer: tsize 'TIME ', default: 1: 'size of time dimension in index specified'
% Integer: esize 'ELEMENTS ', default: 1: 'size of elements dimension in index specified'
%
% Example: o = kkimportasc({i1, i2}, {'i1','';'i2','';'so',0;'sf',0;'nr',1;'dim',5;'o','';'wsize',1;'hsize',1;'dsize',1;'tsize',1;'esize',1})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kimportasc - Import ASCII Data File to Data Segment
%
%  DESCRIPTION
% This routine will import an ASCII data file and insert it into a
% specified data object.  The new segment will be inserted into this
% data object.  If a source object was specified, then all of its data
% segments and the global object attributes will be copied to the
% destination.  If no source object was specified, then the resulting
% output object will contain only the imported segment.
% There are four attributes of the new segment that can be specified:
% dimension, index order, size and data type.  The number of indices
% and sizes specified must be equal to the dimension.
% The parameter dimension must be greater than 0 and and less than or
% equal to 5.
% The segment size is an optional argument specified by five different
% size variables: 'wsize', 'hsize', 'dsize', 'tsize', and 'esize'.  Each size
% variable must be >= to 1.  Note that although the size variables are
% optional, if no size variables are specified, the output will be of
% size 1.
% The index order is also specified by five different variables:
% 'windex', 'hindex', 'dindex', 'tindex', and 'eindex'.  The value for each
% index variable must be unique and range from 1 to 5.  Note that the index
% order dictates how data is stored and retrieved.  For example if the
% index order is specified as windex 1, hindex 2, dindex 3, tindex 4,
% eindex 5 then data will be accessed along width first, followed by height,
% depth, time, and elements.  If, however the index order is specified as
% eindex 1, windex 2, hindex 3, dindex 4, tindex 5 then data will be accessed
% along elements first, followed by width, height, depth, and time.
% The start offset specifies the point in the ASCII data from which to
% start importing in the data into the specified segment.  All points before
% the start offset are ignored.
% The skip factor specifies the number of points to skip in the ASCII data
% between reads.
% The number of reads specifies the number of points to read before skipping.
% For example, if the number of reads is 2 and the skip factor is 3, two
% values will be read, three will be skipped, followed by two more reads,
% three skips, etc. until all the data has been read.
% If the segment_name is one of the Polymorphic segments, the correct
% dimensionality must be provided :
%
%          Segment Name               Dimension
%        ---------------------------------------------------
%              KDA_SEGMENT_VALUE         5
%              KDA_SEGMENT_MAP           5
%              KDA_SEGMENT_MASK          5
%              KDA_SEGMENT_LOCATION      4
%              KDA_SEGMENT_TIME          1
%
% When importing ASCII data to either the value, map, or mask segment, all
% five size variables (wsize, hsize, dsize, tsize, and esize) may be specified
% since these segments are five dimensional.  Likewise, a five dimensional
% index order may be provided.
% When importing ASCII data to the locations segment, only four size variables
% (wsize, hsize, dsize, and esize) may be specified since the location segment
% is four dimensional.  Likewise only a four dimensional index order (windex,
% hindex, dindex, eindex) may be provided.
% When importing ASCII data to the time segment, only one size variable
% (tsize) may be specified since the location segment is one dimensional.
% Likewise only a one dimensional index order (tindex) may be provided.
% Complex values in ASCII are represented by number pairs, surrounded
% by parenthesis and separated by a comma.  These numbered pairs correspond
% to the real and imaginary values respectively.  For example, (8, 3) represents
% the complex number 8+3i.
%
%
%
%  EXAMPLES
% Consider a simple ASCII input file which looks as follows:
%
% 	file:  "ascii_data"
% 	1 2 3 4
% 	1 2 3 4
% 	1 2 3 4
% 	1 2 3 4
% 	1 2 3 4
%
% This data can be imported into the 'value' segment in any number of ways which
% is illustrated in the following four examples.
% .I
% Value Segment Example 1
% % kimportasc -i1 ascii_data -o output -wsize 2 -hsize 5 -dsize 1 -tsize 1 -esize 1 -so 0 -sf 1 -nr 1 -windex 1 -hindex 2 -dindex 3 -tindex 4 -eindex 5 -segment value -type 6
% This will import the above ASCII data into the value segment of a polymorphic
% data object such that the size is "2x5x1x1x1", index order is "w,h,d,t,e",
% and data type is integer.  The output from kimportasc would conceptually
% look as follows:
%
% 	1 3
% 	1 3
% 	1 3
% 	1 3
% 	1 3
%
% .I
% Value Segment Example 2
% % kimportasc -i1 ascii_data -o output -wsize 1 -hsize 5 -dsize 1 -tsize 1 -esize 1 -so 1 -sf 3 -nr 1 -windex 1 -hindex 2 -dindex 3 -tindex 4 -eindex 5 -segment value -type 4
% This will import the above ASCII data into the value segment of a polymorphic
% data object such that the size is "1x5x1x1x1", index order is "w,h,d,t,e",
% and data type is short.  The output from kimportasc would conceptually
% look as follows:
%
% 	2
% 	2
% 	2
% 	2
% 	2
%
% .I
% Value Segment Example 3
% % kimportasc -i1 ascii_data -o output -wsize 3 -hsize 5 -dsize 1 -tsize 1 -esize 1 -so 1 -sf 1 -nr 3 -windex 1 -hindex 2 -dindex 3 -tindex 4 -eindex 5 -segment value -type 10
% This will import the above ASCII data into the value segment of a polymorphic
% data object such that the size is "3x5x1x1x1", index order is "w,h,d,t,e",
% and data type is float.  The output from kimportasc would conceptually
% look as follows:
%
% 	2 3 4
% 	2 3 4
% 	2 3 4
% 	2 3 4
% 	2 3 4
%
% .I
% Value Segment Example 4
% % kimportasc -i1 ascii_data -o output -wsize 5 -hsize 4 -dsize 1 -tsize 1 -esize
%  1 -so 0 -sf 0 -nr 20 -segment value -type 8
% This will import the above ASCII data into the value segment of a polymorphic
% data object such that the size is "4x5x1x1x1", index order is "w,h,d,t,e",
% and data type is long.  The output from kimportasc would conceptually look
% as follows:
%
% 	1 2 3 4 1
% 	2 3 4 1 2
% 	3 4 1 2 3
% 	4 1 2 3 4
%
% The same input data can be imported into the 'map' segment which is
% illustrated in the following example.
% .I
% Map Segment Example
% % kimportasc -i1 ascii_data -i2 data_without_map -o output -wsize 4 -hsize 5 -dsize 1 -tsize 1 -esize 1 -so 0 -sf 0 -nr 1 -windex 1 -hindex 2 -dindex 3 -tindex 4 -eindex 5 -segment map -type 6
% This will import the above ASCII data into the map segment of the polymorphic
% data object specified by the input file "data_without_map" such that the size
% of the map is "4x5x1x1x1", index order is "mw,mh,md,mt,me", and data type is
% integer, where "mw,mh,md,mt,me" coorespond to map width, map height, map depth,
% map time, and map elements respectively.  The map segment output from
% kimportasc would conceptually look like the input ASCII data.
%
% 	1 2 3 4
% 	1 2 3 4
% 	1 2 3 4
% 	1 2 3 4
% 	1 2 3 4
%
% .I
% Mask Segment Example
% The following input data can be imported into the 'mask' segment which is
% illustrated in the following example.
%
% 	file:  "map_ascii_data"
% 	1 0 0 1
% 	0 0 1 0
% 	1 1 0 0
% 	1 1 1 1
% 	1 0 1 0
%
% % kimportasc -i1 map_ascii_data -i2 data_without_mask -o output -wsize 4 -hsize 5 -dsize 1 -tsize 1 -esize 1 -so 0 -sf 0 -nr 1 -windex 1 -hindex 2 -dindex 3 -tindex 4 -eindex 5 -segment mask -type 6
% This will import the above ASCII data into the mask segment of the polymorphic
% data object specified by the input file "data_without_mask" such that the size
% of the mask is "4x5x1x1x1", index order is "w,h,d,t,e" and data type is
% integer.  The mask segment output from kimportasc would conceptually look
% like the ASCII input.
% .I
% Time Segment Example
% The following input data can be imported into the 'time' segment which is
% illustrated in the following example.
%
% 	file: "time_ascii_data"
% 	1 5 10 15 19 36 40 52
%
% % kimportasc -i1 time_ascii_data -i2 data_with_time_series -o output -tsize 8 -so 0 -sf 0 -nr 1 -tindex 1 -segment time -type 6 -dim 1
% This will import the above ASCII data into the time segment of the polymorphic
% data object specified by the input file "data_with_time_series" such that the
% size of the explicit time segment is "8", index order is "t", and data type is
% integer.  The time segment output from kimportasc would conceptually look
% like the ASCII input.
% .I
% Location Segment Example
% The following input data can be imported into the 'location'segment which
% is illustrated in the following example.
%
% 	file:  "location_ascii_data"
% 	1  3  5  7
% 	9  11 13 15
% 	17 19 21 23
% 	25 27 29 31
% 	33 35 37 39
%
% % kimportasc -i1 location_ascii_data -i2 data_without_locations -o output -wsize 4 -hsize 5 -dsize 1 -esize 1 -so 0 -sf 0 -nr 1 -windex 1 -hindex 2 -dindex 3 -eindex 4 -segment location -type 6 -dim 4
% This will import the above ASCII data into the location segment of the
% polymorphic data object specified by the input file "data_without_locations"
% such that the size of the explicit location segment is "4x5x1x1", index
% order is "w,h,d,e", and data type is integer.  The location segment output
% from kimportasc would conceptually look like the ASCII input.
%
%  "SEE ALSO"
% kasc2val, kasc2mask, kasc2loc, kasc2time, kasc2map
%
%  RESTRICTIONS
% If no size variables are specified (wsize, hsize, dsize, tsize, or esize),
% the output data object will be of size 1.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkimportasc(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkimportasc(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'so', 0;'sf', 0;'nr', 1;'dim', 5;'o', '__output';'wsize', 1;'hsize', 1;'dsize', 1;'tsize', 1;'esize', 1};
maxval={0,1,1,1,2,2,0,2,2,2,2,2};
minval={0,1,1,1,2,2,0,2,2,2,2,2};
istoggle=[0,1,1,1,1,1,0,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Integer','Integer','Integer','Integer','OutputFile','Integer','Integer','Integer','Integer','Integer'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kimportasc"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function score = computeIntegralImageScores(integralImage,windows)

windows = round(windows);
windows(windows == 0) = 1;
%windows = [xmin ymin xmax ymax]
%computes the score of the windows wrt the integralImage
height = size(integralImage,1);
index1 = height*windows(:,3) + (windows(:,4) + 1);
index2 = height*(windows(:,1) - 1) + windows(:,2);
index3 = height*(windows(:,1) - 1) + (windows(:,4) + 1);
index4 = height*windows(:,3) + windows(:,2);
score = integralImage(index1) + integralImage(index2) - integralImage(index3) - integralImage(index4);

author: Cloud-CV

repo: object-proposals


text:
function [out] = evap_10(p1,S,Smax,Ep,dt)
%evap_10

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

% Flux function
% ------------------
% Description:  Evaporation from bare soil scaled by relative storage
% Constraints:  Ea <= Ep
%               Ea <= S/dt
% @(Inputs):    p1   - fraction of area that is bare soil [-]
%               S    - current storage [mm]
%               Smax - maximum storage [mm]
%               Ep   - potential evapotranspiration rate [mm/d]
%               dt   - time step size [d]

out = max(min(p1.*S./Smax.*Ep,S/dt),0);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [lat2,lon2,derr,aerr]=vdistinv(lat1,lon1,dist,azim)
% VDISTINV - Invert VDIST function using numerical inversion
%
% Usage:
%
% [lat2,lon2] = vdistinv(lat1,lon1,dist,azim)
% [lat2,lon2,derr,aerr] = vdistinv(lat1,lon1,dist,azim)
%
% Variables:
%
% lat1, lon1 = coordinates of intial point in degrees
% dist       = geodesic distance in meters
% azim       = geodesic azimuth in degrees clockwise from north
% lat2,lon2  = destination coordinates in degrees
% derr       = optional output: error between input distance and the
%              calculated length of the path to the computed endpoint
%              (as computed by VDIST)
% aerr       = optional output: error between the provided azimuth and the
%              calculated azimuth (as computed by VDIST)
%
% Notes: (1) This "quick and dirty" approach was written in response to
%            a user request. The use of a numerical optimization to invert
%            VDIST is a relatively slow and crude approach. (It would be
%            better to implement an algorithm written specifically for that
%            purpose by Vincenty or others.)
%        (2) The distance between essentially antipodal points on an
%            ellipsoid is very sensitive to small deviations in azimuth,
%            so such points should be avoided. (A warning is given.)
%        (3) For other cases, precision is set to about one part in 10^12
%        (3) Tested but no warranty; use at your own risk.
%        (4) Written by Michael Kleder, April 2006
%
% Example:
% >> [dist,azim]=vdist(10,20,30,40)
% dist =    3035728.95690893
% azim =    40.3196402221127
% >> [lat2,lon2]=vdistinv(10,20,dist,azim)
% lat2 =    29.9999999999981
% lon2 =    39.9999999999979
% >>

% initial guess for path endpoint is computed using spherical earth trig:
t1=lat1*0.0174532925199433; % degrees to radians
n1=lon1*0.0174532925199433; % degrees to radians
a=azim*0.0174532925199433; % degrees to radians
d=dist*1.56961230576048e-007; % meters to radians
lat2 = asin(sin(t1)*cos(d)+cos(t1)*sin(d)*cos(a));
lon2 = n1+atan2(sin(d)*sin(a),cos(t1)*cos(d)-sin(t1)*sin(d)*cos(a));
X=[lat2;lon2]*57.2957795130823; % radians to dgrees
% other parameters (start point, arc length, and azimuth) are fixed:
params=[lat1;lon1;dist;azim];
% optimization control settings:
opt=optimset('MaxFunEvals',5000,'TolFun',1e-12);
% solve for accurate end point:
X=fminsearch(@tryone,X,opt,params);
% recover coordinates of endpoint:
lat2=X(1);
lon2=X(2);
if nargout > 2 % if error data is requested
    % compute distance and azimuth from startpoint to endpoint
    [d,a] = vdist(lat1,lon1,lat2,lon2);
    % error in distance:
    derr=abs(dist-d);
    % error in azimuth:
    azim = mod(azim,360);
    a = mod(a,360);
    aerr = abs(azim-a);
    aerr = min(aerr,abs(360-aerr));
end
return
function err=tryone(X,params)
lat2=X(1);
lon2=X(2);
lat1=params(1);
lon1=params(2);
dist=params(3);
azim=params(4);
% crude catch for out-of-bounds attempts:
if lat1<-90 || lat1 > 90 || lat2 < -90 || lat2 > 90
    err = 1e6;
    return
end
[trydist,tryazim]=vdist(lat1,lon1,lat2,lon2);
% compute overall error. First convert distance to approximate arc length
% in degrees so as to weight the terms reasonably closely. (Both move to
% zero in the optimization, but reasonably balancing the units can help.)
err = sqrt((9e-6*(dist-trydist))^2 + (azim-tryazim)^2);
return
function varargout = vdist(lat1,lon1,lat2,lon2)
% VDIST - Using the WGS-84 Earth ellipsoid, compute the distance between
%         two points within a few millimeters of accuracy, compute forward
%         azimuth, and compute backward azimuth, all using a vectorized
%         version of Vincenty's algorithm.
%
% s = vdist(lat1,lon1,lat2,lon2)
% [s,a12] = vdist(lat1,lon1,lat2,lon2)
% [s,a12,a21] = vdist(lat1,lon1,lat2,lon2)
%
% s = distance in meters (inputs may be scalars, vectors, or matrices)
% a12 = azimuth in degrees from first point to second point (forward)
% a21 = azimuth in degrees from second point to first point (backward)
%       (Azimuths are in degrees clockwise from north.)
% lat1 = GEODETIC latitude of first point (degrees)
% lon1 = longitude of first point (degrees)
% lat2, lon2 = second point (degrees)
%
%  Original algorithm source:
%  T. Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid
%  with Application of Nested Equations", Survey Review, vol. 23, no. 176,
%  April 1975, pp 88-93.
%  Available at: http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
%
% Notes: (1) lat1,lon1,lat2,lon2 can be any (identical) size/shape. Outputs
%            will have the same size and shape.
%        (2) Error correcting code, convergence failure traps, antipodal
%            corrections, polar error corrections, WGS84 ellipsoid
%            parameters, testing, and comments: Michael Kleder, 2004.
%        (3) Azimuth implementation (including quadrant abiguity
%            resolution) and code vectorization, Michael Kleder, Sep 2005.
%        (4) Vectorization is convergence sensitive; that is, quantities
%            which have already converged to within tolerance are not
%            recomputed during subsequent iterations (while other
%            quantities are still converging).
%        (5) Vincenty describes his distance algorithm as precise to within
%            0.01 millimeters, subject to the ellipsoidal model.
%        (6) For distance calculations, essentially antipodal points are
%            treated as exactly antipodal, potentially reducing accuracy
%            slightly.
%        (7) Distance failures for points exactly at the poles are
%            eliminated by moving the points by 0.6 millimeters.
%        (8) The Vincenty distance algorithm was transcribed verbatim by
%            Peter Cederholm, August 12, 2003. It was modified and
%            translated to English by Michael Kleder.
%            Mr. Cederholm's website is http://www.plan.aau.dk/~pce/
%        (9) Distances agree with the Mapping Toolbox, version 2.2 (R14SP3)
%            with a max relative difference of about 5e-9, except when the
%            two points are nearly antipodal, and except when one point is
%            near the equator and the two longitudes are nearly 180 degrees
%            apart. This function (vdist) is more accurate in such cases.
%            For example, note this difference (as of this writing):
%            >>vdist(0.2,305,15,125)
%            18322827.0131551
%            >>distance(0.2,305,15,125,[6378137 0.08181919])
%            0
%       (10) Azimuths FROM the north pole (either forward starting at the
%            north pole or backward when ending at the north pole) are set
%            to 180 degrees by convention. Azimuths FROM the south pole are
%            set to 0 degrees by convention.
%       (11) Azimuths agree with the Mapping Toolbox, version 2.2 (R14SP3)
%            to within about a hundred-thousandth of a degree, except when
%            traversing to or from a pole, where the convention for this
%            function is described in (10), and except in the cases noted
%            above in (9).
%       (12) No warranties; use at your own risk.

% reshape inputs
keepsize = size(lat1);
lat1=lat1(:);
lon1=lon1(:);
lat2=lat2(:);
lon2=lon2(:);
% Input check:
if any(abs(lat1)>90 | abs(lat2)>90)
    error('Input latitudes must be between -90 and 90 degrees, inclusive.')
end
% Supply WGS84 earth ellipsoid axis lengths in meters:
a = 6378137; % definitionally
b = 6356752.31424518; % computed from WGS84 earth flattening coefficient
% preserve true input latitudes:
lat1tr = lat1;
lat2tr = lat2;
% convert inputs in degrees to radians:
lat1 = lat1 * 0.0174532925199433;
lon1 = lon1 * 0.0174532925199433;
lat2 = lat2 * 0.0174532925199433;
lon2 = lon2 * 0.0174532925199433;
% correct for errors at exact poles by adjusting 0.6 millimeters:
kidx = abs(pi/2-abs(lat1)) < 1e-10;
if any(kidx);
    lat1(kidx) = sign(lat1(kidx))*(pi/2-(1e-10));
end
kidx = abs(pi/2-abs(lat2)) < 1e-10;
if any(kidx)
    lat2(kidx) = sign(lat2(kidx))*(pi/2-(1e-10));
end
f = (a-b)/a;
U1 = atan((1-f)*tan(lat1));
U2 = atan((1-f)*tan(lat2));
lon1 = mod(lon1,2*pi);
lon2 = mod(lon2,2*pi);
L = abs(lon2-lon1);
kidx = L > pi;
if any(kidx)
    L(kidx) = 2*pi - L(kidx);
end
lambda = L;
lambdaold = 0*lat1;
itercount = 0;
notdone = logical(1+0*lat1);
alpha = 0*lat1;
sigma = 0*lat1;
sinsigma=nan*lat1;
cossigma=nan*lat1;
cos2sigmam = 0*lat1;
C = 0*lat1;
warninggiven = false;
while any(notdone)  % force at least one execution
    %disp(['lambda(21752) = ' num2str(lambda(21752),20)]);
    itercount = itercount+1;
    if itercount > 50
        if ~warninggiven
            warning('VDIST:antipodal',['Essentially antipodal points ' ...
                'encountered. Precision may be reduced.']);
        end
        lambda(notdone) = pi;
        break
    end
    lambdaold(notdone) = lambda(notdone);
    sinsigma(notdone) = sqrt((cos(U2(notdone)).*sin(lambda(notdone)))...
        .^2+(cos(U1(notdone)).*sin(U2(notdone))-sin(U1(notdone)).*...
        cos(U2(notdone)).*cos(lambda(notdone))).^2);
    cossigma(notdone) = sin(U1(notdone)).*sin(U2(notdone))+...
        cos(U1(notdone)).*cos(U2(notdone)).*cos(lambda(notdone));
    % eliminate rare imaginary portions at limit of numerical precision:
    sinsigma(notdone)=real(sinsigma(notdone));
    cossigma(notdone)=real(cossigma(notdone));
    sigma(notdone) = atan2(sinsigma(notdone),cossigma(notdone));
    alpha(notdone) = asin(cos(U1(notdone)).*cos(U2(notdone)).*...
        sin(lambda(notdone))./sin(sigma(notdone)));
    cos2sigmam(notdone) = cos(sigma(notdone))-2*sin(U1(notdone)).*...
        sin(U2(notdone))./cos(alpha(notdone)).^2;
    C(notdone) = f/16*cos(alpha(notdone)).^2.*(4+f*(4-3*...
        cos(alpha(notdone)).^2));
    lambda(notdone) = L(notdone)+(1-C(notdone)).*f.*sin(alpha(notdone))...
        .*(sigma(notdone)+C(notdone).*sin(sigma(notdone)).*...
        (cos2sigmam(notdone)+C(notdone).*cos(sigma(notdone)).*...
        (-1+2.*cos2sigmam(notdone).^2)));
    %disp(['then, lambda(21752) = ' num2str(lambda(21752),20)]);
    % correct for convergence failure in the case of essentially antipodal
    % points
    if any(lambda(notdone) > pi)
        warning('VDIST:antipodal',['Essentially antipodal points ' ...
            'encountered. Precision may be reduced.']);
        warninggiven = true;
        lambdaold(lambda>pi) = pi;
        lambda(lambda>pi) = pi;
    end
    notdone = abs(lambda-lambdaold) > 1e-12;
end
u2 = cos(alpha).^2.*(a^2-b^2)/b^2;
A = 1+u2./16384.*(4096+u2.*(-768+u2.*(320-175.*u2)));
B = u2./1024.*(256+u2.*(-128+u2.*(74-47.*u2)));
deltasigma = B.*sin(sigma).*(cos2sigmam+B./4.*(cos(sigma).*(-1+2.*...
    cos2sigmam.^2)-B./6.*cos2sigmam.*(-3+4.*sin(sigma).^2).*(-3+4*...
    cos2sigmam.^2)));
varargout{1} = reshape(b.*A.*(sigma-deltasigma),keepsize);
if nargout > 1
    % From point #1 to point #2
    % correct sign of lambda for azimuth calcs:
    lambda = abs(lambda);
    kidx=sign(sin(lon2-lon1)) .* sign(sin(lambda)) < 0;
    lambda(kidx) = -lambda(kidx);
    numer = cos(U2).*sin(lambda);
    denom = cos(U1).*sin(U2)-sin(U1).*cos(U2).*cos(lambda);
    a12 = atan2(numer,denom);
    kidx = a12<0;
    a12(kidx)=a12(kidx)+2*pi;
    % from poles:
    a12(lat1tr <= -90) = 0;
    a12(lat1tr >= 90 ) = pi;
    varargout{2} = reshape(a12 * 57.2957795130823,keepsize); % to degrees
end
if nargout > 2
    a21=NaN*lat1; %#ok this variable won't be computed if not needed
    % From point #2 to point #1
    % correct sign of lambda for azimuth calcs:
    lambda = abs(lambda);
    kidx=sign(sin(lon1-lon2)) .* sign(sin(lambda)) < 0;
    lambda(kidx)=-lambda(kidx);
    numer = cos(U1).*sin(lambda);
    denom = sin(U1).*cos(U2)-cos(U1).*sin(U2).*cos(lambda);
    a21 = atan2(numer,denom);
    kidx=a21<0;
    a21(kidx)= a21(kidx)+2*pi;
    % backwards from poles:
    a21(lat2tr >= 90) = pi;
    a21(lat2tr <= -90) = 0;
    varargout{3} = reshape(a21 * 57.2957795130823,keepsize); % to degrees
end
return


text:
import matlab.unittest.TestRunner;
import matlab.unittest.TestSuite;
import matlab.unittest.plugins.TestReportPlugin;
import matlab.unittest.plugins.XMLPlugin
import matlab.unittest.plugins.DiagnosticsValidationPlugin

try
    suite = testsuite({'DemoTests'});

    runner = matlab.unittest.TestRunner.withTextOutput('OutputDetail',1);
    runner.addPlugin(DiagnosticsValidationPlugin)

    xmlFile = 'BSPDemoTests.xml';
    plugin = XMLPlugin.producingJUnitFormat(xmlFile);
    runner.addPlugin(plugin);

    results = runner.run(suite);

    t = table(results);
    disp(t);
    disp(repmat('#',1,80));
    for test = results
        if test.Failed
            disp(test.Name);
        end
    end
catch e
    disp(getReport(e,'extended'));
    bdclose('all');
    exit(1);
end

save(['BSPInstallerTest_',datestr(now,'dd_mm_yyyy-HH:MM:SS'),'.mat'],'t');
bdclose('all');
exit(any([results.Failed]));

author: analogdevicesinc


text:
function [c] = datasets_feature(dataset_names, train_lists, test_lists, feature, c)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

if(~exist('c', 'var'))
  c = conf();
end

if (~exist('OCTAVE_VERSION','builtin'))
    openPool(c.cores);
end

if(c.common_dictionary)
    c.feature_config.(feature).dictionary = build_dictionary(train_lists, feature, c);
end
cache_folder = c.cache;
idx = randperm(length(train_lists));

for j=1:length(train_lists)
  i = idx(j);
  vprintf(c.verbosity, 0, 'Dataset: %s\n', dataset_names{i});
  c.cache = [cache_folder '/' dataset_names{i} '/'];

  if(~c.common_dictionary)
      c.feature_config.(feature).dictionary = build_dictionary(train_lists{i}, feature, c);
  end

  batch_feature(train_lists{i}, 'train', feature, c);
  batch_feature(test_lists{i}, 'test', feature, c);
end

c.cache = cache_folder;

author: adikhosla

repo: feature-extraction


text:
classdef SMMOP6 < PROBLEM
% <multi/many> <real> <large/none> <multimodal> <sparse/none>
% Sparse multi-modal multi-objective optimization problem
% theta --- 0.1 --- Sparsity of the Pareto sets
% np    ---   4 --- Number of the Pareto sets

%------------------------------- Reference --------------------------------
% Y. Tian, R. Liu, X. Zhang, H. Ma, K. C. Tan, and Y. Jin, A
% multipopulation evolutionary algorithm for solving large-scale multimodal
% multiobjective optimization problems, IEEE Transactions on Evolutionary
% Computation, 2021, 25(3): 405-418.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties(Access = private)
        theta = 0.1;    % Sparsity of the Pareto sets
        np    = 4;    	% Number of the Pareto sets
        POS;            % Pareto optimal set for IGDX calculation
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            [obj.theta,obj.np] = obj.ParameterSet(0.1,4);
            if isempty(obj.M); obj.M = 2; end
            if isempty(obj.D); obj.D = 100; end
            obj.lower    = [zeros(1,obj.M-1)+0,zeros(1,obj.D-obj.M+1)-1];
            obj.upper    = [zeros(1,obj.M-1)+1,zeros(1,obj.D-obj.M+1)+2];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            [N,D] = size(X);
            M     = obj.M;
            S     = ceil(obj.theta*(D-M));
            g     = zeros(N,obj.np);
            for i = 1 : obj.np
                g(:,i) = sum(g2(X(:,M+(i-1)*S:M+i*S-1),pi/3),2)+sum(g4(X(:,[M:M+(i-1)*S-1,M+i*S:end]),0),2);
            end
            PopObj = repmat(1+min(g,[],2)/(D-M+1),1,M).*fliplr(cumprod([ones(N,1),1-cos(X(:,1:M-1)*pi/2)],2)).*[ones(N,1),1-sin(X(:,M-1:-1:1)*pi/2)];
        end
        %% Generate Pareto optimal solutions
        function R = GetOptimum(obj,N)
            % Generate points in Pareto optimal set
            A       = GetPS(obj.D-obj.M+1,obj.np,ceil(obj.theta*(obj.D-obj.M)));
            X       = UniformPoint(N/size(A,1),obj.M-1,'grid');
            obj.POS = [repmat(X,size(A,1),1),A(repmat(1:end,size(X,1),1),:)];
            % Generate points on Pareto front
            R = UniformPoint(N,obj.M);
            c = ones(size(R,1),obj.M);
            for i = 1 : size(R,1)
                for j = 2 : obj.M
                    temp = R(i,j)/R(i,1)*prod(1-c(i,obj.M-j+2:obj.M-1));
                    c(i,obj.M-j+1) = (temp^2-temp+sqrt(2*temp))/(temp^2+1);
                end
            end
            x = acos(c)*2/pi;
            R = fliplr(cumprod([ones(size(x,1),1),1-cos(x(:,1:end-1)*pi/2)],2)).*[ones(size(x,1),1),1-sin(x(:,end-1:-1:1)*pi/2)];
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                a = linspace(0,pi/2,100)';
                R = [1-cos(a),1-sin(a)];
            elseif obj.M == 3
                a = linspace(0,pi/2,10)';
                R = {(1-cos(a))*(1-cos(a')),(1-cos(a))*(1-sin(a')),(1-sin(a))*ones(size(a'))};
            else
                R = [];
            end
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case 'IGDX'
                    score = feval(metName,Population,obj.POS);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Display a population in the objective space
        function DrawObj(obj,Population)
            PopDec = Population.decs;
            A      = GetPS(obj.D-obj.M+1,obj.np,ceil(obj.theta*(obj.D-obj.M)));
            [~,Label]  = min(pdist2(PopDec(:,obj.M:end),A),[],2);
            tempStream = RandStream('mlfg6331_64','Seed',2);
            if obj.M == 2
                for i = 1 : size(A,1)
                    color = rand(tempStream,1,3);
                    Draw(Population(Label==i).objs+(i-1)*0.05,'o','MarkerSize',6,'Marker','o','Markerfacecolor',sqrt(color),'Markeredgecolor',color,{'\it f\rm_1','\it f\rm_2',[]});
                    Draw(obj.PF+(i-1)*0.05,'-','LineWidth',1,'Color',color);
                end
            elseif obj.M == 3
                for i = 1 : size(A,1)
                    color = rand(tempStream,1,3);
                    ax = Draw(Population(Label==i).objs+(i-1)*0.05,'o','MarkerSize',8,'Marker','o','Markerfacecolor',sqrt(color),'Markeredgecolor',color,{'\it f\rm_1','\it f\rm_2','\it f\rm_3'});
                    surf(ax,obj.PF{1}+(i-1)*0.05,obj.PF{2}+(i-1)*0.05,obj.PF{3}+(i-1)*0.05,'EdgeColor',color,'FaceColor','none');
                end
            else
                for i = 1 : size(A,1)
                    Draw(Population(Label==i).objs,'-','Color',rand(tempStream,1,3),'LineWidth',2);
                end
            end
        end
    end
end

function g = g2(x,t)
    g = 2*(x-t).^2 + sin(2*pi*(x-t)).^2;
end

function g = g4(x,t)
    g = sqrt((x-t).^2)+sin(2*pi*(x-t)).^2;
end

function PS = GetPS(D,np,S)
    PS = zeros(np,D);
    for i = 1 : np
        PS(i,(i-1)*S+1:i*S) = 1;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [model, llh] = hmmEm(x, init)
% EM algorithm to fit the parameters of HMM model (a.k.a Baum-Welch algorithm)
% Input:
%   x: 1 x n integer vector which is the sequence of observations
%   init: model or k
% Output:s
%   model: trained model structure
%   llh: loglikelihood
% Written by Mo Chen (sth4nth@gmail.com).
n = size(x,2);
X = sparse(x,1:n,1);
d = size(X,1);
if isstruct(init)   % init with a model
    A = init.A;
    E = init.E;
    s = init.s;
elseif numel(init) == 1  % random init with latent k
    k = init;
    s = normalize(rand(k,1),1);
    A = normalize(rand(k,k),2);
    E = normalize(rand(k,d),2);
end
tol = 1e-4;
maxIter = 1000;
llh = -inf(1,maxIter);
for iter = 2:maxIter
    M = E*X;
%     E-step
    [gamma,alpha,beta,c] = hmmSmoother(M,A,s);
    llh(iter) = mean(log(c));
    if abs(llh(iter)-llh(iter-1)) < tol*abs(llh(iter-1)); break; end   % check likelihood for convergence
%     M-step
    s = gamma(:,1);                                                                             % 13.18
    A = normalize(A.*(alpha(:,1:n-1)*(beta(:,2:n).*M(:,2:n)./c(2:n))'),2);      % 13.19 13.43 13.65
    E = (gamma*X')./sum(gamma,2);                            % 13.23
end
model.s = s;
model.A = A;
model.E = E;
llh = llh(2:iter);

function [gamma, alpha, beta, c] = hmmSmoother(M, A, s)
[K,T] = size(M);
At = A';
c = zeros(1,T);
alpha = zeros(K,T);
[alpha(:,1),c(1)] = normalize(s.*M(:,1),1);
for t = 2:T
    [alpha(:,t),c(t)] = normalize((At*alpha(:,t-1)).*M(:,t),1);  % 13.59
end
beta = ones(K,T);
for t = T-1:-1:1
    beta(:,t) = A*(beta(:,t+1).*M(:,t+1))/c(t+1);   % 13.62
end
gamma = alpha.*beta;                  % 13.64

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
%% Copyright (C) 2014, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym divergence (@var{F})
%% @defmethodx @@sym divergence (@var{F}, @var{x})
%% Symbolic divergence of symbolic expression.
%%
%% Consider a vector expression @var{F}:
%% @example
%% @group
%% syms f_1(x,y,z) f_2(x,y,z) f_3(x,y,z)
%% F = [f_1; f_2; f_3]
%%   @result{} F = (sym 3×1 matrix)
%%       ⎡f₁(x, y, z)⎤
%%       ⎢           ⎥
%%       ⎢f₂(x, y, z)⎥
%%       ⎢           ⎥
%%       ⎣f₃(x, y, z)⎦
%% @end group
%% @end example
%% The divergence of @var{F} is the scalar expression:
%% @example
%% @group
%% divergence(F)
%%   @result{} (sym)
%%       ∂                 ∂                 ∂
%%       ──(f₁(x, y, z)) + ──(f₂(x, y, z)) + ──(f₃(x, y, z))
%%       ∂x                ∂y                ∂z
%% @end group
%% @end example
%%
%% Examples:
%% @example
%% @group
%% syms x y
%% F = [x^2/2  y^2/2];
%% divergence(F)
%%   @result{} (sym) x + y
%% @end group
%%
%% @group
%% syms z
%% F = [y x x*y];
%% divergence(F, [x; y; z])
%%   @result{} (sym) 0
%% @end group
%% @end example
%%
%% Note: assumes @var{x} is a Cartesian coordinate system.
%%
%% @seealso{@@sym/gradient, @@sym/curl, @@sym/laplacian, @@sym/jacobian,
%%          @@sym/hessian}
%% @end defmethod

function g = divergence(f, x)

  assert (isvector(f), 'divergence: defined for vectors')

  if (nargin == 1)
    x = symvar(f);
  elseif (nargin == 2)
    % np-op
  else
    print_usage ();
  end

  assert (length(f) == length(x), 'divergence: num vars must match vec length')

  idx1.type='()';
  if (iscell(x))
    idx2.type='{}';
  else
    idx2.type='()';
  end
  g = sym(0);
  for i = 1:length(f)
    idx1.subs={i};
    idx2.subs={i};
    g = g + diff (subsref(f,idx1), subsref(x,idx2));
  end

end

%!shared x,y,z
%! syms x y z

%!test
%! % 1D
%! f = x^2;
%! assert (isequal (divergence(f), diff(f,x)))
%! assert (isequal (divergence(f,{x}), diff(f,x)))
%! assert (isequal (divergence(f,[x]), diff(f,x)))
%! assert (isequal (divergence(f,x), diff(f,x)))

%!test
%! % const
%! f = [sym(1); 2; exp(sym(3))];
%! assert (isequal (divergence(f,{x,y,z}), 0))
%! f = [sym(1); 2; exp(sym('c'))];
%! assert (isequal (divergence(f,{x,y,z}), 0))

%!test
%! % double const
%! f = [1 2];
%! g = sym(0);
%! assert (isequal (divergence(f, [x y]), g))
%! % should fail, calls @double: divergence(f, {x y}), g))

%!test
%! % 1D fcn in 2d/3d
%! f = [x y z];
%! assert (isequal (divergence(f), 3))
%! assert (isequal (divergence(f, {x,y,z}), 3))
%! assert (isequal (divergence(f, [x,y,z]), 3))

%!test
%! % 2d fcn in 2d/3d
%! f = sin(exp(x)*y+sinh(z));
%! g2 = [diff(f,x); diff(f,y)];
%! l2 = diff(g2(1),x) + diff(g2(2),y);
%! g3 = [diff(f,x); diff(f,y); diff(f,z)];
%! l3 = diff(g3(1),x) + diff(g3(2),y) + diff(g3(3),z);
%! assert (isequal (divergence(g2, {x,y}), l2))
%! assert (isequal (divergence(g3, {x,y,z}), l3))

%!error divergence ([1 2], [sym('x')])
%!error divergence ([1 2], sym('x'), 42)

repo: octsympy


text:
%% kronecker
% http://www.mathworks.com/matlabcentral/fileexchange/24499-kronecker
%
function K = kronecker(A,B)
%KRONECKER   Kronecker tensor product.
%   KRONECKER(X,Y) is the Kronecker tensor product of X and Y.
%   The result is a large matrix formed by taking all possible
%   products between the elements of X and those of Y. For
%   example, if X is 2 by 3, then KRONECKER(X,Y) is
%
%      [ X(1,1)*Y  X(1,2)*Y  X(1,3)*Y
%        X(2,1)*Y  X(2,2)*Y  X(2,3)*Y ]
%
%   If either X or Y is sparse, only nonzero elements are multiplied
%   in the computation, and the result is sparse.
%
%   Class support for inputs X,Y:
%      float: double, single
%
%   NOTE: This function does exactly what Matlab KRON does, but for large
%      full matrices, the engine uses BSXFUN to accelerate the calculation.
%      Another advantage is no intermediates large matrices is generated
%      (four temporary arrays in case of KRON)
%
%   Benchmark on Intel Core2 Duo T7250 @2GHz and 2Go RAM
%   Size A/B  Speed gain
%      10       1.17
%      20       3.48
%      30       3.78
%      40       3.73
%      50       3.68
%      60       4.22
%      70       3.81
%
%   Restriction: MATLAB 2007A or later is required
%
%   See also: KRON
%
%   Author: Bruno Luong <brunoluong@yahoo.com>
%   History:
%       Original 21-Jun-2009

if ~issparse(A) && ~issparse(B)
    if ndims(A) > 2 || ndims(B) > 2
        error('kronecker:TwoDInput','Inputs must be 2-D.');
    end
    % Both inputs are full, result is full. This is faster than
    % MATLAB stock kron (indexing based)
    [ma na] = size(A);
    [mb nb] = size(B);
    A = reshape(A,[1 ma 1 na]);
    B = reshape(B,[mb 1 nb 1]);
    K = bsxfun(@times,A,B);
    K = reshape(K,[ma*mb na*nb]);

else % One of the input matrix is sparse

    % Call MATLAB stock KRON
    K = kron(A,B);

end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [vertex,face,d,c] = read_ply(filename)

% read_ply - read data from PLY file.
%
%   [vertex,face] = read_ply(filename);
%
%   'vertex' is a 'nb.vert x 3' array specifying the position of the vertices.
%   'face' is a 'nb.face x 3' array specifying the connectivity of the mesh.
%
%
%   Copyright (c) 2003 Gabriel Peyr?

[d,c] = plyread(filename);
vi = d.face.vertex_indices;
nf = length(vi);
% http://stackoverflow.com/a/6210539/148668
% This is slow...
%lengths = cellfun(@(x)numel(x),vi);
lengths = cellfun('length',vi);
maxLength=max(lengths);
if all(maxLength == lengths)
  face = cell2mat(vi)+1;
else
  face = cell2mat( ...
    cellfun( ...
      @(x)cat(2,x,-ones(1,maxLength-length(x))), ...
      vi, ...
      'UniformOutput',false))+1;
end
vertex = [d.vertex.x, d.vertex.y, d.vertex.z];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Elements,varargout] = plyread(Path,Str)
%PLYREAD   Read a PLY 3D data file.
%   [DATA,COMMENTS] = PLYREAD(FILENAME) reads a version 1.0 PLY file
%   FILENAME and returns a structure DATA.  The fields in this structure
%   are defined by the PLY header; each element type is a field and each
%   element property is a subfield.  If the file contains any comments,
%   they are returned in a cell string array COMMENTS.
%
%   [TRI,PTS] = PLYREAD(FILENAME,'tri') or
%   [TRI,PTS,DATA,COMMENTS] = PLYREAD(FILENAME,'tri') converts vertex
%   and face data into triangular connectivity and vertex arrays.  The
%   mesh can then be displayed using the TRISURF command.
%
%   Note: This function is slow for large mesh files (+50K faces),
%   especially when reading data with list type properties.
%
%   Example:
%   [Tri,Pts] = PLYREAD('cow.ply','tri');
%   trisurf(Tri,Pts(:,1),Pts(:,2),Pts(:,3));
%   colormap(gray); axis equal;
%
%   See also: PLYWRITE

% Pascal Getreuer 2004

[fid,Msg] = fopen(Path,'rt');	% open file in read text mode

if fid == -1, error(Msg); end

Buf = fscanf(fid,'%s',1);
if ~strcmp(Buf,'ply')
   fclose(fid);
   error('Not a PLY file.');
end

%%% read header %%%

Position = ftell(fid);
Format = '';
NumComments = 0;
Comments = {};				% for storing any file comments
NumElements = 0;
NumProperties = 0;
Elements = [];				% structure for holding the element data
ElementCount = [];		% number of each type of element in file
PropertyTypes = [];		% corresponding structure recording property types
ElementNames = {};		% list of element names in the order they are stored in the file
PropertyNames = [];		% structure of lists of property names

while 1
   Buf = fgetl(fid);   								% read one line from file
   BufRem = Buf;
   Token = {};
   Count = 0;

   while ~isempty(BufRem)								% split line into tokens
      [tmp,BufRem] = strtok(BufRem);

      if ~isempty(tmp)
         Count = Count + 1;							% count tokens
         Token{Count} = tmp;
      end
   end

   if Count 		% parse line
      switch lower(Token{1})
      case 'format'		% read data format
         if Count >= 2
            Format = lower(Token{2});

            if Count == 3 & ~strcmp(Token{3},'1.0')
               fclose(fid);
               error('Only PLY format version 1.0 supported.');
            end
         end
      case 'comment'		% read file comment
         NumComments = NumComments + 1;
         Comments{NumComments} = '';
         for i = 2:Count
            Comments{NumComments} = [Comments{NumComments},Token{i},' '];
         end
      case 'element'		% element name
         if Count >= 3
            if isfield(Elements,Token{2})
               fclose(fid);
               error(['Duplicate element name, ''',Token{2},'''.']);
            end

            NumElements = NumElements + 1;
            NumProperties = 0;
   	      Elements = setfield(Elements,Token{2},[]);
            PropertyTypes = setfield(PropertyTypes,Token{2},[]);
            ElementNames{NumElements} = Token{2};
            PropertyNames = setfield(PropertyNames,Token{2},{});
            CurElement = Token{2};
            ElementCount(NumElements) = str2double(Token{3});

            if isnan(ElementCount(NumElements))
               fclose(fid);
               error(['Bad element definition: ',Buf]);
            end
         else
            error(['Bad element definition: ',Buf]);
         end
      case 'property'	% element property
         if ~isempty(CurElement) & Count >= 3
            NumProperties = NumProperties + 1;
            eval(['tmp=isfield(Elements.',CurElement,',Token{Count});'],...
               'fclose(fid);error([''Error reading property: '',Buf])');

            if tmp
               error(['Duplicate property name, ''',CurElement,'.',Token{2},'''.']);
            end

            % add property subfield to Elements
            eval(['Elements.',CurElement,'.',Token{Count},'=[];'], ...
               'fclose(fid);error([''Error reading property: '',Buf])');
            % add property subfield to PropertyTypes and save type
            eval(['PropertyTypes.',CurElement,'.',Token{Count},'={Token{2:Count-1}};'], ...
               'fclose(fid);error([''Error reading property: '',Buf])');
            % record property name order
            eval(['PropertyNames.',CurElement,'{NumProperties}=Token{Count};'], ...
               'fclose(fid);error([''Error reading property: '',Buf])');
         else
            fclose(fid);

            if isempty(CurElement)
               error(['Property definition without element definition: ',Buf]);
            else
               error(['Bad property definition: ',Buf]);
            end
         end
      case 'end_header'	% end of header, break from while loop
         break;
      end
   end
end

%%% set reading for specified data format %%%

if isempty(Format)
	warning('Data format unspecified, assuming ASCII.');
   Format = 'ascii';
end

switch Format
case 'ascii'
   Format = 0;
case 'binary_little_endian'
   Format = 1;
case 'binary_big_endian'
   Format = 2;
otherwise
   fclose(fid);
   error(['Data format ''',Format,''' not supported.']);
end

if ~Format
   Buf = fscanf(fid,'%f');		% read the rest of the file as ASCII data
   BufOff = 1;
else
   % reopen the file in read binary mode
   fclose(fid);

   if Format == 1
      fid = fopen(Path,'r','ieee-le.l64');		% little endian
   else
      fid = fopen(Path,'r','ieee-be.l64');		% big endian
   end

   % find the end of the header again (using ftell on the old handle doesn't give the correct position)
   BufSize = 8192;
   Buf = [blanks(10),char(fread(fid,BufSize,'uchar')')];
   i = [];
   tmp = -11;

   while isempty(i)
   	i = findstr(Buf,['end_header',13,10]);			% look for end_header + CR/LF
   	i = [i,findstr(Buf,['end_header',10])];		% look for end_header + LF

      if isempty(i)
         tmp = tmp + BufSize;
         Buf = [Buf(BufSize+1:BufSize+10),char(fread(fid,BufSize,'uchar')')];
      end
   end

   % seek to just after the line feed
   fseek(fid,i + tmp + 11 + (Buf(i + 10) == 13),-1);
end

%%% read element data %%%

% PLY and MATLAB data types (for fread)
PlyTypeNames = {'char','uchar','short','ushort','int','uint','float','double', ...
   'char8','uchar8','short16','ushort16','int32','uint32','float32','double64'};
MatlabTypeNames = {'schar','uchar','int16','uint16','int32','uint32','single','double'};
SizeOf = [1,1,2,2,4,4,4,8];	% size in bytes of each type

for i = 1:NumElements
   % get current element property information
   eval(['CurPropertyNames=PropertyNames.',ElementNames{i},';']);
   eval(['CurPropertyTypes=PropertyTypes.',ElementNames{i},';']);
   NumProperties = size(CurPropertyNames,2);

%   fprintf('Reading %s...\n',ElementNames{i});

   if ~Format	%%% read ASCII data %%%
      for j = 1:NumProperties
         Token = getfield(CurPropertyTypes,CurPropertyNames{j});

         if strcmpi(Token{1},'list')
            Type(j) = 1;
         else
            Type(j) = 0;
			end
      end

      % parse buffer
      if ~any(Type)
         % no list types
         Data = reshape(Buf(BufOff:BufOff+ElementCount(i)*NumProperties-1),NumProperties,ElementCount(i))';
         BufOff = BufOff + ElementCount(i)*NumProperties;
      else
         ListData = cell(NumProperties,1);

         for k = 1:NumProperties
            ListData{k} = cell(ElementCount(i),1);
         end

         % list type
		   for j = 1:ElementCount(i)
   	      for k = 1:NumProperties
      	      if ~Type(k)
         	      Data(j,k) = Buf(BufOff);
            	   BufOff = BufOff + 1;
	            else
   	            tmp = Buf(BufOff);
      	         ListData{k}{j} = Buf(BufOff+(1:tmp))';
         	      BufOff = BufOff + tmp + 1;
            	end
            end
         end
      end
   else		%%% read binary data %%%
      % translate PLY data type names to MATLAB data type names
      ListFlag = 0;		% = 1 if there is a list type
      SameFlag = 1;     % = 1 if all types are the same

      for j = 1:NumProperties
         Token = getfield(CurPropertyTypes,CurPropertyNames{j});

         if ~strcmp(Token{1},'list')			% non-list type
	         tmp = rem(strmatch(Token{1},PlyTypeNames,'exact')-1,8)+1;

            if ~isempty(tmp)
               TypeSize(j) = SizeOf(tmp);
               Type{j} = MatlabTypeNames{tmp};
               TypeSize2(j) = 0;
               Type2{j} = '';

               SameFlag = SameFlag & strcmp(Type{1},Type{j});
	         else
   	         fclose(fid);
               error(['Unknown property data type, ''',Token{1},''', in ', ...
                     ElementNames{i},'.',CurPropertyNames{j},'.']);
         	end
         else											% list type
            if length(Token) == 3
               ListFlag = 1;
               SameFlag = 0;
               tmp = rem(strmatch(Token{2},PlyTypeNames,'exact')-1,8)+1;
               tmp2 = rem(strmatch(Token{3},PlyTypeNames,'exact')-1,8)+1;

               if ~isempty(tmp) & ~isempty(tmp2)
                  TypeSize(j) = SizeOf(tmp);
                  Type{j} = MatlabTypeNames{tmp};
                  TypeSize2(j) = SizeOf(tmp2);
                  Type2{j} = MatlabTypeNames{tmp2};
	       else
   	   	 fclose(fid);
               	 error(['Unknown property data type, ''list ',Token{2},' ',Token{3},''', in ', ...
                   ElementNames{i},'.',CurPropertyNames{j},'.']);
               end
            else
               fclose(fid);
               error(['Invalid list syntax in ',ElementNames{i},'.',CurPropertyNames{j},'.']);
            end
         end
      end

      % read file
      if ~ListFlag
         if SameFlag
            % no list types, all the same type (fast)
            Data = fread(fid,[NumProperties,ElementCount(i)],Type{1})';
         else
            % no list types, mixed type
            Data = zeros(ElementCount(i),NumProperties);

         	for j = 1:ElementCount(i)
        			for k = 1:NumProperties
               	Data(j,k) = fread(fid,1,Type{k});
              	end
         	end
         end
      else
         ListData = cell(NumProperties,1);

         for k = 1:NumProperties
            ListData{k} = cell(ElementCount(i),1);
         end

         if NumProperties == 1
            BufSize = 512;
            SkipNum = 4;
            j = 0;

            % list type, one property (fast if lists are usually the same length)
            while j < ElementCount(i)
               Position = ftell(fid);
               % read in BufSize count values, assuming all counts = SkipNum
               [Buf,BufSize] = fread(fid,BufSize,Type{1},SkipNum*TypeSize2(1));
               Miss = find(Buf ~= SkipNum);					% find first count that is not SkipNum
               fseek(fid,Position + TypeSize(1),-1); 		% seek back to after first count

               if isempty(Miss)									% all counts are SkipNum
                  Buf = fread(fid,[SkipNum,BufSize],[int2str(SkipNum),'*',Type2{1}],TypeSize(1))';
                  fseek(fid,-TypeSize(1),0); 				% undo last skip

                  for k = 1:BufSize
                     ListData{1}{j+k} = Buf(k,:);
                  end

                  j = j + BufSize;
                  BufSize = floor(1.5*BufSize);
               else
                  if Miss(1) > 1									% some counts are SkipNum
                     [Buf2,cc] = fread(fid,[SkipNum,Miss(1)-1],[int2str(SkipNum),'*',Type2{1}],TypeSize(1));
                     Buf2 = Buf2';

                     for k = 1:Miss(1)-1
                        ListData{1}{j+k} = Buf2(k,:);
                     end

                     j = j + k;
                  end
                  % Alec: check if done and rewind one step
                  if j >= ElementCount(i)
                    fseek(fid,-1,0);
                    break;
                  end

                  % read in the list with the missed count
                  SkipNum = Buf(Miss(1));
                  j = j + 1;
                  ListData{1}{j} = fread(fid,[1,SkipNum],Type2{1});
                  BufSize = ceil(0.6*BufSize);
               end
            end
         else
            % list type(s), multiple properties (slow)
            Data = zeros(ElementCount(i),NumProperties);

            for j = 1:ElementCount(i)
         		for k = 1:NumProperties
            		if isempty(Type2{k})
               		Data(j,k) = fread(fid,1,Type{k});
            		else
               		tmp = fread(fid,1,Type{k});
               		ListData{k}{j} = fread(fid,[1,tmp],Type2{k});
		            end
      		   end
      		end
         end
      end
   end

   % put data into Elements structure
   for k = 1:NumProperties
   	if (~Format & ~Type(k)) | (Format & isempty(Type2{k}))
      	eval(['Elements.',ElementNames{i},'.',CurPropertyNames{k},'=Data(:,k);']);
      else
      	eval(['Elements.',ElementNames{i},'.',CurPropertyNames{k},'=ListData{k};']);
		end
   end
end

clear Data ListData;
fclose(fid);

if (nargin > 1 & strcmpi(Str,'Tri')) | nargout > 2
   % find vertex element field
   Name = {'vertex','Vertex','point','Point','pts','Pts'};
   Names = [];

   for i = 1:length(Name)
      if any(strcmp(ElementNames,Name{i}))
         Names = getfield(PropertyNames,Name{i});
         Name = Name{i};
         break;
      end
   end

   if any(strcmp(Names,'x')) & any(strcmp(Names,'y')) & any(strcmp(Names,'z'))
      eval(['varargout{1}=[Elements.',Name,'.x,Elements.',Name,'.y,Elements.',Name,'.z];']);
   else
      varargout{1} = zeros(1,3);
	end

   varargout{2} = Elements;
   varargout{3} = Comments;
   Elements = [];

   % find face element field
   Name = {'face','Face','poly','Poly','tri','Tri'};
   Names = [];

   for i = 1:length(Name)
      if any(strcmp(ElementNames,Name{i}))
         Names = getfield(PropertyNames,Name{i});
         Name = Name{i};
         break;
      end
   end

   if ~isempty(Names)
      % find vertex indices property subfield
	   PropertyName = {'vertex_indices','vertex_indexes','vertex_index','indices','indexes'};

   	for i = 1:length(PropertyName)
      	if any(strcmp(Names,PropertyName{i}))
         	PropertyName = PropertyName{i};
	         break;
   	   end
      end

      if ~iscell(PropertyName)
         % convert face index lists to triangular connectivity
         eval(['FaceIndices=varargout{2}.',Name,'.',PropertyName,';']);
  			N = length(FaceIndices);
   		Elements = zeros(N*2,3);
   		Extra = 0;

			for k = 1:N
   			Elements(k,:) = FaceIndices{k}(1:3);

   			for j = 4:length(FaceIndices{k})
      			Extra = Extra + 1;
	      		Elements(N + Extra,:) = [Elements(k,[1,j-1]),FaceIndices{k}(j)];
   			end
         end
         Elements = Elements(1:N+Extra,:) + 1;
      end
   end
else
   varargout{1} = Comments;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function tests = test_spm_ncTcdf
% Unit Tests for spm_ncTcdf
%__________________________________________________________________________
% Copyright (C) 2018 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_ncTcdf.m 7258 2018-02-14 13:09:46Z guillaume $

tests = functiontests(localfunctions);

function test_spm_ncTcdf_1(testCase)
exp = 0.5; % nctcdf(0,1,0)
act = spm_ncTcdf(0,1,0);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

function test_spm_ncTcdf_2(testCase)
exp = 0.5; % nctcdf(0,4,0)
act = spm_ncTcdf(0,4,0);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

function test_spm_ncTcdf_3(testCase)
exp = 0.022750131948179; % nctcdf(0,4,2)
act = spm_ncTcdf(0,4,2);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

function test_spm_ncTcdf_4(testCase)
exp = 4.218150123125319e-05; % nctcdf(-4,4,2)
act = spm_ncTcdf(-4,4,2);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

function test_spm_ncTcdf_5(testCase)
exp = 0.920822617476195; % nctcdf(5,4,2)
act = spm_ncTcdf(5,4,2);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

function test_spm_ncTcdf_6(testCase)
exp = 0.992600181884681; % nctcdf(10,4,2)
act = spm_ncTcdf(10,4,2);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

function test_spm_ncTcdf_7(testCase)
exp = [0.000000011096548  % nctcdf(-4:10,4,4)
   0.000000032879594
   0.000000140813044
   0.000001174173412
   0.000031671241833
   0.001992604459483
   0.042091256308006
   0.202851164067265
   0.429784541497594
   0.622224194923879
   0.754499466251885
   0.838999243587669
   0.892202382905940
   0.926042984950848
   0.948002154097108]';
act = spm_ncTcdf(-4:10,4,4);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

author: spm


text:
function fibonacci_recursive_test ( )

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

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FIBONACCI_RECURSIVE_TEST\n' );
  fprintf ( 1, '  FIBONACCI_RECURSIVE computes the Fibonacci sequence.\n' );
  fprintf ( 1, '\n' );

  f = fibonacci_recursive ( n );

  for i = 1 : n
    fprintf ( 1, '  %4d  %8d\n', i, f(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function L = statistical_filter(SpikeMat,sdd)
% This functions does statistical filtering on detected spikes by
% thresholding on waveform std and mean, and returns a logical matrix
% containing location of detected outliers (L).

Max_std = sdd.stat_filter.max_std;
Max_mean = sdd.stat_filter.max_mean;

spike_std = std(SpikeMat,0,2);
spike_mean = mean(SpikeMat,2);

% finding locations of threshold crossing points
L = (spike_std > Max_std) | (spike_mean > Max_mean) | (spike_mean < -Max_mean);

% check if perecent of detected outliers exceeds maximum allowed percent
if sum(L)/length(L) > sdd.stat_filter.max_percent_of_outliers
    % limiting number of detected outlier to maximum allowed percent by
    % sorting them based on criterion
    criterion = spike_std/std(spike_std) + abs(spike_mean)/std(spike_mean);
    criterion_s = sort(criterion,'descend');
    L = criterion > criterion_s(ceil(length(L)*sdd.stat_filter.max_percent_of_outliers));
end

author: ramintoosi

repo: ROSS

save_path: github-repos/MATLAB/ramintoosi-ROSS


text:
function [dataIdx,pilotIdx] = helperSubcarrierIndices(cfg,fieldType)
% helperSubcarrierIndices Return the subcarrier indices within the FFT size
%
%   [DATAIDX,PILOTIDX] = helperSubcarrierIndices(CFGFORMAT,FIELDTYPE)
%   returns the data and pilot indices within the FFT size for the
%   specified format configuration object, CFGFORMAT, and specified
%   FIELDTYPE.
%
%   DATAIDX is a column vector containing the indices of data subcarriers.
%
%   PILOTIDX is a column vector containing the indices of pilot
%   subcarriers.
%
%   CFGFORMAT is the format configuration object of type wlanVHTConfig,
%   wlanHTConfig, or wlanNonHTConfig, which specifies the parameters for
%   the VHT, HT-Mixed, and Non-HT formats, respectively.
%
%   FIELDTYPE is a string specifying the format of the field and must be
%   one of 'Legacy','HT' or 'VHT'.
%
%   [DATAIDX,PILOTIDX] = helperSubcarrierIndices(CHANBW) returns indices
%   for the specified channel bandwidth. CHANBW must be one of 'CBW5',
%   'CBW10', 'CBW20', 'CBW40', 'CBW80' or 'CBW160'.
%
%   Example: Return indices for a VHT format configuration.
%
%   cfgVHT = wlanVHTConfig;
%   [dataIdx,pilotIdx] = helperSubcarrierIndices(cfgVHT,'VHT')

% Copyright 2015 The MathWorks, Inc.

%#codegen

if ischar(cfg)
    coder.internal.errorIf( ~any(strcmp(cfg, {'CBW5', 'CBW10', ...
            'CBW20', 'CBW40', 'CBW80', 'CBW160'})), ...
        'wlan:helperSubcarrierIndices:InvalidChBandwidth');

    chanBW = cfg;
else
    validateattributes(cfg,{'wlanVHTConfig','wlanHTConfig','wlanNonHTConfig'}, ...
        {'scalar'},mfilename,'format configuration object');

    coder.internal.errorIf(isa(cfg,'wlanNonHTConfig') && ...
        ~strcmp(cfg.Modulation,'OFDM'), ...
            'wlan:helperSubcarrierIndices:InvalidNonHTModulation')

    chanBW = cfg.ChannelBandwidth;
end

FFTLen = helperFFTLength(chanBW);

fieldType = validatestring(fieldType,{'Legacy','HT','VHT'}, ...
    mfilename,'Field type');
if strcmp(fieldType,'Legacy')
    switch chanBW
      case 'CBW40'
        Num20MHzChan = 2;
      case {'CBW80', 'CBW80+80'}
        Num20MHzChan = 4;
      case 'CBW160'
        Num20MHzChan = 8;
      otherwise   % For 'CBW20', 'CBW10', 'CBW5'
        Num20MHzChan = 1;
    end

    numGuardBands = [6; 5];
    pilotIdx20MHz =  [12; 26; 40; 54];

    % Get non-data subcarrier indices per 20MHz channel bandwidth
    nonDataIdxPerGroup = [(1:numGuardBands(1))'; 33; ...
        (64-numGuardBands(2)+1:64)'; pilotIdx20MHz];
    % Get non-data subcarrier indices for the whole bandwidth
    nonDataIdxAll = bsxfun(@plus, nonDataIdxPerGroup, 64*(0:Num20MHzChan-1));
    dataIdx  = setdiff((1:FFTLen)', sort(nonDataIdxAll(:)));
    pilotIdx = reshape(bsxfun(@plus, pilotIdx20MHz, 64*(0:Num20MHzChan-1)), [], 1);
else % 'HT' & 'VHT'
    switch chanBW
      case 'CBW40'
        numGuardBands = [6; 5];
        customNullIdx = [-1; 1];
        pilotIdx = [-53; -25; -11; 11; 25; 53];
      case 'CBW80'
        numGuardBands = [6; 5];
        customNullIdx = [-1; 1];
        pilotIdx = [-103; -75; -39; -11; 11; 39; 75; 103];
      case 'CBW80+80' % Merge with 80MHz case, if same. Separate for now
        numGuardBands = [6; 5];
        customNullIdx = [-1; 1];
        pilotIdx = [-103; -75; -39; -11; 11; 39; 75; 103];
      case 'CBW160'
        numGuardBands = [6; 5];
        customNullIdx = [(-129:-127)'; (-5:-1)'; (1:5)'; (127:129)'];
        pilotIdx = [-231; -203; -167; -139; -117; -89; -53; -25; 25; 53; 89; 117; 139; 167; 203; 231];
      otherwise  % CBW20
        numGuardBands = [4; 3];
        customNullIdx = [];
        pilotIdx = [-21; -7; 7; 21];
    end

    pilotIdx = pilotIdx + FFTLen/2 + 1; % Convert to 1-based indexing
    customNullIdx = customNullIdx + FFTLen/2 + 1; % Convert to 1-based indexing
    % Get non-data subcarrier indices for the whole bandwidth
    nonDataIdx = [(1:numGuardBands(1))'; FFTLen/2+1; ...
        (FFTLen-numGuardBands(2)+1:FFTLen)'; pilotIdx; customNullIdx];
    dataIdx = setdiff((1:FFTLen)', sort(nonDataIdx));
end
end

author: seemoo-lab


text:
function bw2 = hough_bin_pixels(bw, theta, rho, bin)
%HOUGH_BIN_PIXELS Find pixels corresponding to Hough accumulator bin.
%   BW2 = hough_bin_pixels(BW, THETA, RHO, BIN) finds the white pixels in a
%   binary image that correspond to a particular Hough transform
%   accumulator bin.  BW is the original binary image.  THETA and RHO are
%   the Hough parameter vectors returned by the hough function.  BIN is a
%   two-element vector containing the row-column coordinates of the Hough
%   transform bin.  BW2 is a binary image containing only the white pixels
%   in BW that contributed to the specified bin.
%
%   Example
%   =======
%   I  = imread('circuit.tif');
%   BW = edge(I,'canny');
%   imshow(BW)
%   [H,theta,rho] = hough(BW);
%   P = houghpeaks(H, 1);
%   BW2 = hough_bin_pixels(BW, theta, rho, P);
%   figure, imshow(BW2)
%   title('Pixels corresponding to maximum Hough transform bin')
%
%   See also hough, houghlines, houghpeaks.

%   Steven L. Eddins
%   The MathWorks, Inc.

[y, x] = find(bw);
x = x - 1;
y = y - 1;

theta_c = theta(bin(2)) * pi / 180;
rho_xy = x*cos(theta_c) + y*sin(theta_c);
nrho = length(rho);
slope = (nrho - 1)/(rho(end) - rho(1));
rho_bin_index = round(slope*(rho_xy - rho(1)) + 1);

idx = find(rho_bin_index == bin(1));

r = y(idx) + 1;
c = x(idx) + 1;

bw2 = false(size(bw));

bw2(sub2ind(size(bw), r, c)) = true;


text:
%converts frequency to mel
%>
%> @param fInHz: frequency
%> @param cModel: 'Fant','Shaughnessy', or 'Umesh'
%>
%> @retval mel pitch value
% ======================================================================
function [mel] = ToolFreq2Mel(fInHz, cModel)

    if (nargin < 2)
        cModel = 'Fant';
    end

    % set function handle
    hPitchFunc = str2func (['aca' cModel '_I']);

    mel = hPitchFunc(fInHz);
end

% Fant
function [mel] = acaFant_I(f)
    mel = 1000 * log2(1 + f/1000);
end

% Shaughnessy
function [mel] = acaShaughnessy_I(f)
    mel = 2595 * log10(1 + f/700);
end

% Umesh
function [mel] = acaUmesh_I(f)
    mel = f./(2.4e-4*f + 0.741);
end

author: alexanderlerch

repo: ACA-Code


text:
function matRad_writeVTK(filepath,cube,metadata)
% matRad function to write vtk cubes
%
% call
%   matRad_writeVTK(filepath,cube,metadata)
%
% input
%   filepath:   full filename (with extension)
%   cube:       3D array to be written into file
%   metadata:   struct of metadata. Writer will wrap the existing metadata
%               to VTK standard-specific fields
%               Necessary fieldnames are:
%               - resolution: [x y z]
%               - datatype: numeric MATLAB-Datatype
%
% output
%   file will be written to disk
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sanity checks and restrictions
dimensions = size(cube);
if numel(dimensions) ~= 3
    error('Sorry! matRad only supports 3-dimensional VTK output');
end

fid = fopen(filepath, 'wb');
if fid <= 0
    error('Could not open VTK destination file!');
end
cleaner = onCleanup(@() fclose(fid));

%We perform
if isfield(metadata,'axisPermutation')
    cube = permute(cube,metadata.axisPermutation);
end

fprintf(fid, '# vtk DataFile Version 3.0\n');
fprintf(fid, 'vtk output\n');
fprintf(fid, 'BINARY\n');
fprintf(fid, 'DATASET STRUCTURED_POINTS\n');
fprintf(fid, 'DIMENSIONS    %d   %d   %d\n', dimensions(1),dimensions(2),dimensions(3));
fprintf(fid, 'SPACING    %f   %f   %f\n',metadata.resolution(1),metadata.resolution(2),metadata.resolution(3));
fprintf(fid, 'ORIGIN    %f   %f   %f\n',metadata.imageOrigin(1),metadata.imageOrigin(2),metadata.imageOrigin(3));
fprintf(fid, 'POINT_DATA   %d\n',prod(dimensions));
if isfield(metadata,'dataName')
    dataName = metadata.dataName;
else
    dataName = 'scalars';
end
fprintf(fid, 'SCALARS %s %s\n',dataName,metadata.datatype);
fprintf(fid, 'LOOKUP_TABLE default\n');
fwrite(fid,cube,metadata.datatype,'b');

end

repo: matRad


text:
% LQRSIM: Algebraic Riccati Equation solution in Simulink
% Version 2.1 (R14SP1)  1-Jan-2009
%
% Giampiero Campa, Jason Hall, Riccardo Bevilacqua
%
%
% EXAMPLE :
%
% sfcndemo_lqry     : main example.
%

author: Sable


text:
function boolout = iscphd( filename )
%ISCPHD Checks to see whether input file is in the CPHD format
%
% Written by: Wade Schwartzkopf, NGA/Research
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

fid = fopen(filename,'r','b','UTF-8');
[cphd_str, version_str] = strtok(fgets( fid, 20 ),'/'); % First line is just CPHD/version
vers_parts = str2double(regexp(version_str(2:end),'\.','split'));
boolout = strcmp(cphd_str,'CPHD') && ...  % CPHD"X"
    ((vers_parts(1)==0&&vers_parts(2)>=3)||vers_parts(1)>0);  % 0.3 and above
fclose(fid);

end

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

author: ngageoint


% example_2_min_new.m VAR inference with Minnesota Prior
% Author:   Filippo Ferroni and  Fabio Canova
% Date:     27/02/2020, revised  14/10/2020
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute  posterior estimate with
% 1) mixed calibrated/estimated Minnesota prior
% 2) optimally chosen  Minnesota prior.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
warning off; clear; close all; clc;

addpath ../../cmintools/
addpath ../../bvartools/

% load Euro data
load Data
y= [IPI HICP CORE Euribor1Y M3 EXRATE];% the variables for the VAR

%% Case 1: Estimation with calibrated Minnesota Prior

lags = 6;
options.max_minn_hyper  = 1;       % start the  optimization routine
BVAR0                    = bvar_(y,lags,options);
pause;

%% Case 2: Estimation with calibrated/ estimated Minnesota Prior
lags = 6;
options.max_minn_hyper  = 1;       % start the  optimization routine
options.minn_prior_tau  = 10;      % set tau
options.index_est       = [3 4];   % hyper-parameters to maximize
options.lb              = [0 0];   % the lower bounds
options.ub              = [20 20]; % the upper bounds
options.max_compute     = 7;       % optimization  by Matlab Simplex
BVAR1                   = bvar_(y,lags,options);
pause;

%% Case 3:  Optimization of one dimension of the Minnesota Prior
%           (without  posterior draws to speed up computations)
%%  3.1 maximization of  tau
clear options
lags = 6;
% setting the default values for the hyperparameters
hyperpara(1)    = 3;		  % tau
hyperpara(2)    = 0.5;		  % decay
hyperpara(3)    = 5;		  % lambda
hyperpara(4)    = 2;		  % mu
hyperpara(5)    = 2;		  % omega
% setting the options
options.index_est	   = 1:1;    % hyper-parameter over which maximize
options.max_compute    = 2;      % maximize  using Matlab fmincon function
options.lb             = 0.8;    % Lower bound
options.ub             = 10;     % Upper bound
[postmode,logmlike,~] = bvar_max_hyper(hyperpara,y,lags,options);
pause;

%% 3.2    Take  optimal value for hyperparameter(1) and compute optimal
%%         values for tau, decay, lambda (without  posterior  draws)
hyperpara(1)            = postmode(1); % use as starting value previous mode
options.index_est       = 1:3;         % set hyper-parameters over which maximize
options.lb              = [0.1 0.1 0.1]; % Lower bounds
options.ub              = [10 10 10];    % Upper bounds
[postmode1,log_dnsty,~] = bvar_max_hyper(hyperpara,y,lags,options);

%% 3.3  Take  optimal value for hyperparameter(1:3) and compute optimal
%%           values for tau, decay, lambda, mu (without  posterior  draws)

hyperpara(1:3)          = postmode1(1:3); % use as starting value previous mode
options.index_est       = 1:4;         % set hyper-parameters over which maximize
options.lb              = [0.1 0.1 0.1 0.1]; % Lower bounds
options.ub              = [10 10 10 10];    % Upper bounds
[postmode,log_dnsty1,~] = bvar_max_hyper(hyperpara,y,lags,options);
% then run BVAR with  optimal  parameters
BVAR2                   = bvar_(y,lags,options);

%Plotting cholesky responses  to  Monetary  policy  shocks
% Define the IRF of Interest
indx_sho              = 4;
% Change the order of the variables for the plot
indx_var              = [4, 1, 2, 3];

% % IRFs to PLOT: compare different max hyperparam
mltple_irfs_to_plot_all(:,:,1,:) = BVAR0.ir_draws(indx_var,:,indx_sho,:);
mltple_irfs_to_plot_all(:,:,2,:) = BVAR1.ir_draws(indx_var,:,indx_sho,:);
mltple_irfs_to_plot_all(:,:,3,:) = BVAR2.ir_draws(indx_var,:,indx_sho,:);

% Customize the plot
% variables names for the plots
options.varnames      = {'1 year rate','IP','HICP','CORE'};
% name of the directory where the figure is saved
options.saveas_dir    = './irfs_plt';
% name of the figure to save
options.saveas_strng  = 'diff_minn_opt';
% name of the shock
options.shocksnames   = {'MP opt1','MP opt2','MP opt3'};
% additional 90% HPD set
options.conf_sig_2    = 0.9;
% the plotting command
plot_all_irfs_(mltple_irfs_to_plot_all,options)