text:
classdef TestNet
    %TestNet

    properties (Constant)
        im1 = fullfile(mexopencv.root(), 'test', 'cat.jpg');
        im2 = fullfile(mexopencv.root(), 'test', 'space_shuttle.jpg');
    end

    methods (Static)
        function test_blobs
            imgs = cell(1,4);
            for i=1:numel(imgs)
                imgs{i} = rand(50, 50, 3, 'single');
            end
            blob = cv.Net.blobFromImages(imgs, 'SwapRB',false, 'Crop',false);
            out = cv.Net.imagesFromBlob(blob);
            assert(isequal(out, imgs));
        end

        function test_caffe_googlenet
            % load net and images
            net = load_bvlc_googlenet();
            img1 = cv.imread(TestNet.im1, 'FlipChannels',false);
            img2 = cv.imread(TestNet.im2, 'FlipChannels',false);

            % forward pass (blob from one image)
            blob = cv.Net.blobFromImages(img1, 'Size',[224 224]);
            net.setInput(blob, 'data');
            prob = net.forward('prob');
            validateattributes(prob, {'numeric'}, ...
                {'vector', 'real', 'numel',1000, '>=',0, '<=',1});
            [p,idx] = max(prob,[],2);
            %assert(abs(sum(prob) - 1) < 1e-4);

            % forward pass (blob from multiple images)
            blob = cv.Net.blobFromImages({img1,img2}, 'Size',[224 224]);
            net.setInput(blob);
            prob = net.forward();
            validateattributes(prob, {'numeric'}, ...
                {'size',[2 1000], 'real', '>=',0, '<=',1});
            [p,idx] = max(prob,[],2);
            %assert(all(abs(sum(prob,2) - 1) < 1e-4));
        end

        function test_intermediate_blobs
            net = load_bvlc_googlenet();
            img1 = cv.imread(TestNet.im1, 'FlipChannels',false);
            img2 = cv.imread(TestNet.im2, 'FlipChannels',false);

            blob = cv.Net.blobFromImages(img1, 'Size',[224 224]);
            net.setInput(blob, 'data');
            blobs = net.forward({
                'conv1/7x7_s2'
                'conv1/relu_7x7'
                'inception_4c/1x1'
                'inception_4c/relu_1x1'
            });

            blob = cv.Net.blobFromImages({img1,img2}, 'Size',[224 224]);
            net.setInput(blob, 'data');
            blobs = net.forwardAndRetrieve('conv1/7x7_s2');
        end

        function test_shrink_caffe_fp16
            model = fullfile(mexopencv.root(), 'test', 'dnn', 'GoogLeNet', ...
                'bvlc_googlenet.caffemodel');
            model_fp16 = fullfile(tempdir(), 'bvlc_googlenet_fp16.caffemodel');
            if exist(model, 'file') ~= 2
                error('mexopencv:testskip', 'missing data');
            end

            cObj = onCleanup(@() delete(model_fp16));
            cv.Net.shrinkCaffeModel(model, model_fp16);
            assert(exist(model_fp16, 'file') == 2);
        end

        function test_layers
            % Convolution layer
            % (64 filters each of size 7x7, 3x3 padding, 2x2 stride)
            % output size of convolution is: (H + 2*pad - kernel_size)/stride + 1
            lp = struct();
            lp.name = 'conv/7x7_s2';
            lp.type = 'Convolution';
            lp.dict.num_output = int32(64);
            lp.dict.kernel_size = int32(7);
            lp.dict.stride = int32(2);
            lp.dict.pad = int32(3);
            lp.dict.bias_term = true;
            lp.blobs = {
                rand(64,3,7,7,'single')*2-1,  % weights
                rand(1,1,1,64,'single')*2-1   % biases
            };

            % add layer
            net = cv.Net();
            lp_id = net.addLayer(lp.name, lp.type, lp);
            validateattributes(lp_id, {'int32'}, {'scalar', 'integer'});
            assert(lp_id ~= -1);

            % connect input layer to conv layer
            net.setInputsNames('data');
            if true
                net.connect(0, 0, lp_id, 0);
            elseif true
                %net.connect('_input.0', lp.name);
                net.connect('_input', lp.name);
            else
                net.connect('_input.data', lp.name);
            end

            % feed net batch of two 224x224 BGR images, and run forward pass
            net.setInput(rand(2,3,224,224,'single')*255, 'data');
            out = net.forward(lp.name);
            assert(isequal(size(out), [2 64 112 112]))

            b0 = net.getParam(lp_id, 0);
            validateattributes(b0, {'numeric'}, {'real'});
            %assert(isequal(b0, lp.blobs{1}));

            b1 = net.getParam(lp.name, 1);
            validateattributes(b1, {'numeric'}, {'real'});
            %assert(isequal(b1, lp.blobs{2}));

            net.setParam(lp_id, 0, lp.blobs{1});
            net.setParam(lp.name, 1, lp.blobs{2});

            net.setPreferableBackend('Default') % 'Halide'
            net.setPreferableTarget('CPU')      % 'OpenCL'
            net.enableFusion(true);

            id = net.getUnconnectedOutLayers();
            validateattributes(id, {'numeric'}, {'vector', 'integer'});
            %assert(isequal(id, lp_id))  % output layer

            id = net.getLayerId('_input');
            assert(id == 0);  % special input layer with id=0

            types = net.getLayerTypes();
            assert(iscellstr(types));

            counts = cellfun(@(t) net.getLayersCount(t), types);
            validateattributes(counts, {'numeric'}, {'vector', 'integer'});

            names = net.getLayerNames();
            assert(iscellstr(names));
            %assert(isequal(names{1}, lp.name))

            id = net.getLayerId(names{1});
            validateattributes(id, {'numeric'}, {'scalar', 'integer'});
            assert(lp_id ~= -1);
            %assert(isequal(id, lp_id))

            if true
                layer = net.getLayer(id);
            else
                layer = net.getLayer(names{1});
            end
            assert(isstruct(layer) && isscalar(layer));
            validateattributes(layer.blobs, {'cell'}, {'vector', 'numel',2});
            %assert(isequal(layer.name, lp.name))
            %assert(isequal(layer.type, lp.type))
            %assert(isequal(layer.blobs(:), lp.blobs(:)))

            layers = net.getLayerInputs(id);
            assert(isstruct(layers));

            %if true
            %    net.deleteLayer(lp_id);
            %else
            %    net.deleteLayer(lp.name);
            %end
        end
    end

end

function net = load_bvlc_googlenet()
    %TODO: download files (or run the sample "caffe_googlenet_demo.m")
    rootdir = fullfile(mexopencv.root(), 'test', 'dnn', 'GoogLeNet');
    modelTxt = fullfile(rootdir, 'deploy.prototxt');
    modelBin = fullfile(rootdir, 'bvlc_googlenet.caffemodel');
    if exist(modelTxt, 'file') ~= 2 || exist(modelBin, 'file') ~= 2
        error('mexopencv:testskip', 'missing data');
    end

    net = cv.Net('Caffe', modelTxt, modelBin);
    assert(~net.empty());
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Script to run a MARRMoT model through the BMI interface

%% 1. Create config file

% Get some data
load MARRMoT_example_data.mat

% Catchment settings
% -------------------------------------------------------------------------

% 'data_origin' is a 2x1 array with [lat,lon]
data_origin = [35.29,87.49];

% 'forcing' is a structure with fields 'precip, 'temp', 'pet',
% 'delta_t_days', 'time_unit'
forcing.precip       = data_MARRMoT_examples.precipitation;
forcing.temp         = data_MARRMoT_examples.temperature;
forcing.pet          = data_MARRMoT_examples.potential_evapotranspiration;
forcing.delta_t      = 1;   % 1 [d]
forcing.time_unit    = 'day';

% 'time_start' & 'time_end' are a vector with starting/end date
time_start = datevec(data_MARRMoT_examples.dates_as_datenum(1));
time_end   = datevec(data_MARRMoT_examples.dates_as_datenum(end));

% Model settings
% -------------------------------------------------------------------------

%'model_name' is a string with the name of the model function
model_name = 'm_01_collie1_1p_1s';

% 'parameters' is a vector of length n (number of parameters)
parameters = 10; % Smax, [mm]

% 'store_ini' is a vector of length m (number of stores)
store_ini = 5;

% 'solver' is a structure with fields 'name', 'resnorm_tolerance',
% 'resnorm_maxiter'
solver.name = 'createOdeApprox_IE';
solver.resnorm_tolerance = 0.1;
solver.resnorm_maxiter = 6;

% Save as a config file
% -------------------------------------------------------------------------
clear data_MARRMoT_examples
save BMI_testcase_m01_BuffaloRiver_TN_USA

%% 2. Run the model through BMI
clear % start with a clean slate

% Specify config path
filepath = '../Config/BMI_testcase_m01_BuffaloRiver_TN_USA';

% Create a model object
obj_marrmot_m01 = marrmotBMI();

% Initialize the model
obj_marrmot_m01.initialize(filepath)

% Get the number of time steps
ts = datenum(obj_marrmot_m01.startTime);
te = datenum(obj_marrmot_m01.endTime);
tn = te-ts;

% Run a loop
qsim = NaN.*zeros(tn,1);

while obj_marrmot_m01.time <= tn

    % Advance 1 time step
    obj_marrmot_m01.update();

    % Get simulated flow
    tmp = obj_marrmot_m01.get_value('flux_out'); % gives a structure
    qsim(obj_marrmot_m01.time-1) = tmp.Q;

end

%% 3. Do some visualisation

% get observations
load '../Config/MARRMoT_example_data.mat'

% Figure
figure('color','w');
hold on;
plot(data_MARRMoT_examples.streamflow);
plot(qsim)
title([obj_marrmot_m01.model_name, ' at ', num2str(obj_marrmot_m01.lat), ...
    ' lat, ',num2str(obj_marrmot_m01.lon), ' lon'],'interpreter','none')
xlabel('time [d]')
ylabel('discharge [mm/d]')

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function test_suite = test_spatial2

%   Run specific demo and save values for comparison.
%
%   See also
%     TEST_ALL, DEMO_SPATIAL2

initTestSuite;

function testDemo
% Set random number stream so that failing isn't because randomness. Run
% demo & save test values.
prevstream=setrandstream(0);

disp('Running: demo_spatial2')
demo_spatial2
Ef = Ef(1:100);
Varf = Varf(1:100);
C = C(1:50, 1:50);
path = which('test_spatial2.m');
path = strrep(path,'test_spatial2.m', 'testValues');
if ~(exist(path, 'dir') == 7)
    mkdir(path)
end
path = strcat(path, '/testSpatial2');
save(path, 'Ef', 'Varf', 'C');

% Set back initial random stream
setrandstream(prevstream);
drawnow;clear;close all

% Compare test values to real values.

function testPredictionsEP
values.real = load('realValuesSpatial2.mat', 'Ef', 'Varf');
values.test = load(strrep(which('test_spatial2.m'), 'test_spatial2.m', 'testValues/testSpatial2.mat'), 'Ef', 'Varf');
assertElementsAlmostEqual(mean(values.test.Ef), mean(values.real.Ef), 'relative', 0.1);
assertElementsAlmostEqual(mean(values.test.Varf), mean(values.real.Varf), 'relative', 0.1);

function testCovarianceMatrix
values.real = load('realValuesSpatial2.mat', 'C');
values.test = load(strrep(which('test_spatial2.m'), 'test_spatial2.m', 'testValues/testSpatial2.mat'), 'C');
assertElementsAlmostEqual(mean(values.real.C), mean(values.test.C), 'relative', 0.1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function g = fgplvmSequenceGradient(xvec, model, Y, varargin)

% FGPLVMSEQUENCEGRADIENT Wrapper function for gradient of a latent sequence.
% FORMAT
% DESC is a wrapper function for the gradient of the log probability
% of a sequence under the posterior distribution induced by the
% training data with respect to the latent positions.  The GP-LVM
% model is one that is already trained with a specific data set.
% ARG vecx : the positions in the latent space that are being
% optimised as a row vector using the matlab X(:)' notation.
% ARG model : the trained GP-LVM model that is being optimised.
% ARG Y : the position in data space for which the latent sequence is
% being optimised.
% ARG P1, P2, P3 ... : optional additional arguments to be passed
% to the model sequence log likelihood gradient.
% RETURN g : the gradient of the log likelihood with respect to the
% latent positions.
%
% SEEALSO : fgplvmSequenceLogLikeGradient, fgplvmOptimiseSequence
%
% COPYRIGHT Neil D. Lawrence, 2006

% FGPLVM

X = reshape(xvec, size(Y, 1), model.q);
g = - fgplvmSequenceLogLikeGradient(model, X, Y, varargin{:});

g = g(:)';

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_pull434

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY

% this script tests the appending of sensor structures
% and specifically the handling of overlapping channels and/or electrodes
% see https://github.com/fieldtrip/fieldtrip/pull/434
% and the preceding https://github.com/fieldtrip/fieldtrip/pull/428

%% construct some electrode structures

elec1.label = {'1'; '2'; '3'};
elec1.elecpos = [1 1 1; 2 2 2; 3 3 3];
elec1.chanpos = [1 1 1; 2 2 2; 3 3 3];
elec1.tra = eye(3);

elec2.label = {'3'; '4'};
elec2.elecpos = [3 3 3; 4 4 4];
elec2.chanpos = [3 3 3; 4 4 4];
elec2.tra = eye(2);

elec3.label = {'4'; '5'; '6'};
elec3.elecpos = [4 4 4; 5 5 5; 6 6 6];
elec3.chanpos = [4 4 4; 5 5 5; 6 6 6];
elec3.tra = eye(3);

%% append them in different ways

cfg = [];

append1   = ft_appendsens(cfg, elec1);
append11  = ft_appendsens(cfg, elec1, elec1);
append111 = ft_appendsens(cfg, elec1, elec1, elec1);

assert(isequal(append1.label,   elec1.label));
assert(isequal(append1.chanpos,   elec1.chanpos));
assert(isequal(append1.tra,   elec1.tra));
assert(isequal(append11.label,   elec1.label));
assert(isequal(append11.chanpos,   elec1.chanpos));
assert(isequal(append11.tra,   elec1.tra));
assert(isequal(append111.label,   elec1.label));
assert(isequal(append111.chanpos,   elec1.chanpos));
assert(isequal(append111.tra,   elec1.tra));

%%

append13  = ft_appendsens(cfg, elec1, elec3);
append12  = ft_appendsens(cfg, elec1, elec2);
append23  = ft_appendsens(cfg, elec2, elec3);

assert(numel(append13.label)==6);
assert(numel(append12.label)==4);
assert(numel(append23.label)==4);

assert(size(append13.chanpos,1)==6);
assert(size(append12.chanpos,1)==4);
assert(size(append23.chanpos,1)==4);

assert(size(append13.elecpos,1)==6);
assert(size(append12.elecpos,1)==4);
assert(size(append23.elecpos,1)==4);

assert(isequal(size(append13.tra), [6 6]));
assert(isequal(size(append12.tra), [4 4]));
assert(isequal(size(append23.tra), [4 4]));

%%

append123 = ft_appendsens(cfg, elec1, elec2, elec3);

assert(numel(append123.label)==6);
assert(size(append123.chanpos,1)==6);
assert(size(append123.elecpos,1)==6);
assert(isequal(size(append123.tra), [6 6]));

%%

elec1b.label = {'1b', '2b', '3b'};
elec1b.elecpos = [1 1 1; 2 2 2; 3 3 3];      % same electrodes
elec1b.chanpos = [1 1 1; 2 2 2; 3 3 3] + 10; % different channels
elec1b.tra = eye(3);

elec1c.label = {'1c', '2c', '3c'};
elec1c.elecpos = [1 1 1; 2 2 2; 3 3 3] + 10; % different electrodes
elec1c.chanpos = [1 1 1; 2 2 2; 3 3 3];      % same channels
elec1c.tra = eye(3);

%%

append11b = ft_appendsens(cfg, elec1, elec1b); % same elec, diff chan

assert(numel(append11b.label)==6);
assert(size(append11b.chanpos,1)==6);
assert(size(append11b.elecpos,1)==3);
assert(isequal(size(append11b.tra), [6 3]));

%%

if false
append11c = ft_appendsens(cfg, elec1, elec1c); % same chan, diff elec
% this should crash because there are 6 labels and only 3 chanpos
end


text:
function [fTafseq] = calc_Afseqlengtwomol(fPval1,fCval1,fKval1,fBgrate,nMod,fTsecaf,fPval2,fCval2,fKval2)
    % [fTafseq] = calc_Afseqlengtwomol(fPval1,fCval1,fKval1,fBgrate,nMod,fTsecaf,fPval2,fCval2,fKval2)
    % -------------------------------------------------------------------------
    % Determines the length of an aftershock sequence by determining the
    % intersection with the DAILY backgraund rate; MOL and nested MOL are allowed
    % (only one secondary sequence)
    %
    % Incoming:
    % fPval1   : Modified Omori law p-value
    % fCval1   : Modified Omori law c-value [days]
    % fKval1   : Modified Omori law k-value
    % fBgrate  : Daily background rate of events above Mc
    % nMod     : Aftershock sequence model choice
    %               1. Modified Omori law (MOL)
    %               2. MOL with one secondary sequence
    % fTsecaf
    % fPval2   : Modified Omori law p-value
    % fCval2   : Modified Omori law c-value [days]
    % fKval2   : Modified Omori law k-value

    % Outgoing:
    % fTafseq  : Length of aftershock sequence
    %
    % last update: 08.07.04
    % jochen.woessner@sed.ethz.ch

    % Time vector in days
    vT = [0:0.1:10000];
    vT = vT';

    % Check input
    if nargin < 4; disp('Not enough input parameters'); return; end
    if nargin == 4; nMod=1; end

    % Switch between MOL and nested MOL
    switch nMod
        case 1 % MOL
            vRate = abs(fKval1.*(vT+fCval1).^-fPval1-fBgrate);
            vSel = (min(vRate) == vRate);
            fTafseq = vT(vSel);
        case 2 % Nested Models
            vSelT = (vT >= fTsecaf);
            vRate1 = abs(fKval1.*(vT(~vSel)+fCval1).^-fPval1-fBgrate);
            vRate2 = abs(fKval1.*(vT(vSel)+fCval1).^-fPval1+fKval2.*(vT(vSel,:)-fTsecaf+fCval2).^-fPval2-fBgrate);
            vRate = [vRate1; vRate2];
            vSel = (min(vRate) == vRate);
            fTafseq = vT(vSel);
        otherwise
            disp('Not a valid model!');
    end
    % % Possible plot
    % figure
    % loglog(vT,vRate)
    % xlabel('Time [days after mainshock]')
    % ylabel('Daily rate of earthquakes')

author: CelsoReyes


text:
function cache = create_loc_cc_fftn_cache(metric, metric_param_pix, volmov, volfix, internal_dtype, maskfix, deps, loc_cc_approximate)
%     deps = 1e-5;
    Nd = 3;
    if size(volmov, 3) == 1
        Nd = 2;
    end
    szvolfix = [size(volfix, 1), size(volfix, 2), size(volfix, 3)];

    cache = [];
    if numel(metric_param_pix) == 1
        metric_param_pix = metric_param_pix * ones(1,Nd);
    end
    metric_param_pix = metric_param_pix(1:Nd);

    metr_threshold = 0.4;
    metr_threshold = 0.8;
    metric_param_pix(metric_param_pix < metr_threshold) = metr_threshold;
    sgm = metric_param_pix;
    hsz = max(2*ceil(2*sgm) + 1, 5);
    filterRadius = (hsz-1)/2;

    pad_size = ceil(filterRadius);
    start = pad_size+1;

    stop = start + size(volfix) - 1;
    vol_fix_p = padarray(volfix, pad_size, 'symmetric');

    if Nd == 2
        [n1, n2] = ndgrid(1:size(vol_fix_p,1), 1:size(vol_fix_p,2));
    elseif Nd == 3
        [n1, n2, n3] = ndgrid(1:size(vol_fix_p,1), 1:size(vol_fix_p,2), 1:size(vol_fix_p,3));
    end
    c = ceil(size(n1)/2)+1;
    if Nd == 2
        g = exp(- (n1-c(1)).^2/(2*sgm(1)^2) - (n2-c(2)).^2/(2*sgm(2)^2));
    elseif Nd == 3
        g = exp(- (n1-c(1)).^2/(2*sgm(1)^2) - (n2-c(2)).^2/(2*sgm(2)^2) - (n3-c(3)).^2/(2*sgm(3)^2));
    end

    g = g/sum(g(:));
    g = fftshift(g);
    fg = fftn(g, size(vol_fix_p));
    convop = @(fx) (ifftn(fftn(fx) .* fg, 'symmetric'));

    mean_fix = convop(vol_fix_p);
    sgm_fix = convop(vol_fix_p.^2) - mean_fix.^2 + deps;
    sgm_fix = sqrt(sgm_fix);

%     maskfix_pad = [];
    if ~isempty(maskfix)
        maskfix_pad = padarray(maskfix, pad_size, 'symmetric');
        sgm_fix = sgm_fix ./ maskfix_pad;
        sgm_fix_inv = maskfix_pad ./ sgm_fix;
    else
        sgm_fix_inv = 1 ./ sgm_fix;
    end

    cache.fg = fg;
    cache.vol_fix_p = vol_fix_p;
    cache.mean_fix = mean_fix;
    cache.sgm_fix = sgm_fix;
    cache.sgm_fix_inv = sgm_fix_inv;
    cache.start = start;
    cache.stop = stop;

    cache.pad_size = pad_size;
    cache.metric_name = metric;
    cache.sigma = metric_param_pix;
    cache.internal_dtype = internal_dtype;
    cache.deps = deps;
%     cache.maskfix = maskfix;
%     cache.maskfix_pad = maskfix_pad;
    cache.maskfix_pad = [];
    cache.loc_cc_approximate = loc_cc_approximate;

    if strcmp(metric, 'loc_cc_fftn_gpu')
        cache.fg = gpuArray((cache.fg));
        cache.vol_fix_p = gpuArray((cache.vol_fix_p));
        cache.mean_fix = gpuArray((cache.mean_fix));
        cache.sgm_fix = gpuArray((cache.sgm_fix));
        cache.deps = gpuArray(cache.deps);
        cache.sgm_fix_inv = gpuArray(sgm_fix_inv);
%         cache.maskfix = gpuArray(cache.maskfix);
%         cache.maskfix_pad = gpuArray(cache.maskfix_pad);
    elseif strcmp(metric, 'loc_cc_fftn_gpu_single')
        cache.fg = gpuArray(single(cache.fg));
        cache.vol_fix_p = gpuArray(single(cache.vol_fix_p));
        cache.mean_fix = gpuArray(single(cache.mean_fix));
        cache.sgm_fix = gpuArray(single(cache.sgm_fix));
        cache.deps = gpuArray(single(cache.deps));
        cache.sgm_fix_inv = gpuArray(single(sgm_fix_inv));
%         cache.maskfix = gpuArray(single(cache.maskfix));
%         cache.maskfix_pad = gpuArray(single(cache.maskfix_pad));
    elseif strcmp(metric, 'loc_cc_fftn_single')
        cache.fg = (single(cache.fg));
        cache.vol_fix_p = (single(cache.vol_fix_p));
        cache.mean_fix = (single(cache.mean_fix));
        cache.sgm_fix = (single(cache.sgm_fix));
        cache.deps = (single(cache.deps));
        cache.sgm_fix_inv = (single(sgm_fix_inv));
    end
end

repo: pTVreg


text:
close all;
clearvars;
clc;
dataset_file = fullfile(spx.data_dir, 'clustering', ...
    'self_tuning_paper_clustering_data');
data = load(dataset_file);
datasets = data.XX;
raw_data = datasets{1};
num_clusters = data.group_num(1);
X = raw_data(:, 1);
Y = raw_data(:, 2);
figure;
axis equal;
plot(X, Y, '.', 'MarkerSize',16);
export_fig images/demo_sc_1_unscaled.png -r120 -nocrop;
% Scale the raw_data
raw_data = raw_data - repmat(mean(raw_data),size(raw_data,1),1);
raw_data = raw_data/max(max(abs(raw_data)));
X = raw_data(:, 1);
Y = raw_data(:, 2);
figure;
axis equal;
plot(X, Y, '.', 'MarkerSize',16);
export_fig images/demo_sc_1_scaled.png -r120 -nocrop;

% Compute the distance matrix
sqrt_dist_mat = spx.commons.distance.sqrd_l2_distances_rw(raw_data);
%% Scale to use for standard spectral clustering
scale = 0.04;
% Compute the similarity matrix
sim_mat = spx.cluster.similarity.gauss_sim_from_sqrd_dist_mat(sqrt_dist_mat, scale);
% Run clustering on the dataset
clusterer = spx.cluster.spectral.Clustering(sim_mat);
clusterer.NumClusters = num_clusters;
%cluster_labels = clusterer.cluster_symmetric();
cluster_labels = clusterer.cluster_random_walk();
figure;
colors = [1,0,0;0,1,0;0,0,1;1,1,0;1,0,1;0,1,1;0,0,0];
hold on;
axis equal;
for c=1:num_clusters
    % Identify points in this cluster
    points = raw_data(cluster_labels == c, :);
    X = points(:, 1);
    Y = points(:, 2);
    plot(X, Y, '.','Color',colors(c,:), 'MarkerSize',16);
end
export_fig images/demo_sc_1_clustered.png -r120 -nocrop;

author: indigits

repo: sparse-plex


text:
function [Fn,Vn]=patchExtend(varargin)

% function [Fn,Vn]=patchExtend(F,V,Eb,optionStruct)
% ------------------------------------------------------------------------
% This function extends the input surface geometry allong the edges Eb
% using the local (or averaged) surface direction at these edges.
%
% The input consists of: the input surface (patch) data, defined by the
% faces F and vertices V, the boundary edges Eb to extend the surface
% over, and an optionStruct structure containing all remaining and optional
% inputs.
%
% The input structure may contain the following fields:
%
% optionStruct.extendDistance -> %The distance to extend the surface by
%
% optionStruct.numSteps -> The number of steps to use in the extension
% direction, default is the distance divided by the average edge length.
%
% optionStruct.plotOn -> option to turn on/off plotting, default is 0
%
% optionStruct.meshType -> the mesh type ('tri' or 'quad'), default is
% 'tri'.
%
% optionStruct.globalDirection -> the global (should be seen as mean)
% direction for extension, default is empty ([]) and is instead based on
% the input mesh directions
%
% optionStruct.extendMethod=1; %Method to use for extruding (see below),
% default is 1.
%
% Five extend methods have been implemented (default is 1):
%
% extendMethod=1 -> Equal offset allong local direction
% extendMethod=2 -> Equal offset wrst mean direction allong local direction
% extendMethod=3 -> Varying offset allong local ending planar wrt mean direction
% extendMethod=4 -> Equal offset allong mean direction
% extendMethod=5 -> Equal offset allong mean direction ending planar wrt mean direction
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% Created: 2020/04/01
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case 2
        F=varargin{1};
        V=varargin{2};
        Eb=[];
        optionStruct=[];
    case 3
        F=varargin{1};
        V=varargin{2};
        Eb=varargin{3};
        optionStruct=[];
    case 4
        F=varargin{1};
        V=varargin{2};
        Eb=varargin{3};
        optionStruct=varargin{4};
end

%Check for empty boundary edges
if isempty(Eb)
    Eb=patchBoundary(F,V); %use all boundary edges
end

%Check boundary point spacing
pointSpacing=mean(patchEdgeLengths(Eb,V));

%Create default option structure
defaultOptionStruct.numSteps=[];
defaultOptionStruct.plotOn=0;
defaultOptionStruct.meshType=[];
defaultOptionStruct.globalDirection=[];
defaultOptionStruct.extendDistance=[]; %Extend distance
defaultOptionStruct.extendMethod=1; %Method to use for extruding

%Parse option structure
if isempty(optionStruct)
    optionStruct=defaultOptionStruct;
else
    %Check optionStruct against default
    [optionStruct]=structComplete(optionStruct,defaultOptionStruct,1); %Complement provided with default if missing or empty
end

%Access input parameters
numSteps=optionStruct.numSteps;
plotOn=optionStruct.plotOn;
meshType=optionStruct.meshType;
ns_mean=optionStruct.globalDirection;
extendDistance=optionStruct.extendDistance;
extendMethod=optionStruct.extendMethod;

if isempty(extendDistance)
    extendDistance=pointSpacing;
end

if isempty(numSteps)
    numSteps=ceil(extendDistance./pointSpacing);
end
if numSteps<2
    numSteps=2;
end

if isempty(meshType)
    switch size(F,2)
        case 3
            meshType='tri';
        case 4
            meshType='quad';
        otherwise
            meshType='quad';
    end
end

%%

indEdgePoints=unique(Eb(:)); %Vertex indices of all edge points

[~,~,N]=patchNormal(F,V); %Get surface vertex normal vectors
n=N((indEdgePoints),:); %Surface normal directions at edge start points
[Eb,V1]=patchCleanUnused(Eb,V);
numEdgePoints=size(V1,1);

[~,~,ne]=edgeVec(Eb,V1); %Get edge vectors
ne=vecnormalize(ne); %Edge vectors

ns=vecnormalize(cross(ne,n)); %Allong surface extrude directions orthogonal to both
if isempty(ns_mean)
    ns_mean=vecnormalize(mean(ns,1)); %Mean extrude direction
end

switch extendMethod
    case 1 %Equal offset allong local direction
        V2=V1+extendDistance*ns;
    case 2 %Equal offset wrst mean direction allong local direction
        f=extendDistance./dot(ns,ns_mean(ones(size(ns,1),1),:),2);
        V2=V1+f.*ns;
    case 3 %Varying offset allong local ending planar wrt mean direction
        d=dot(ns,ns_mean(ones(size(ns,1),1),:),2);
        f=extendDistance./d;
        dd=dot(V1-mean(V1,1),ns_mean(ones(numEdgePoints,1),:),2);
        [~,indMax]=max(dd);
        dp=dot(V1-V1(indMax,:),ns_mean(ones(numEdgePoints,1),:),2);
        fp=dp./d;
        V2=V1-fp.*ns+f.*ns;
    case 4 %Equal offset allong mean direction
        V2=V1+extendDistance*ns_mean(ones(numEdgePoints,1),:);
    case 5 %Equal offset allong mean direction ending planar wrt mean direction
        dd=dot(V1-mean(V1,1),ns_mean(ones(numEdgePoints,1),:),2);
        [~,indMax]=max(dd);
        fp=dot(V1-V1(indMax,:),ns_mean(ones(numEdgePoints,1),:),2);
        V2=V1-fp.*ns_mean(ones(numEdgePoints,1),:)+extendDistance*ns_mean(ones(numEdgePoints,1),:);
    otherwise
        error('Invalid extrude option provided, use 1-5');
end

%Create quadrilateral faces
X=linspacen(V1(:,1),V2(:,1),numSteps);
Y=linspacen(V1(:,2),V2(:,2),numSteps);
Z=linspacen(V1(:,3),V2(:,3),numSteps);

Vn=[X(:) Y(:) Z(:)]; %Vertex set
Fn=repmat([fliplr(Eb) Eb],numSteps-1,1);
q=repmat(0:(numSteps-2),size(Eb,1),1);
q=q(:);
Q=[q q q+1 q+1]*size(V1,1);
Fn=Fn+Q;

switch meshType
    case 'tri'
        [Fn,Vn]=quad2tri(Fn,Vn,'a');
    case 'quad'

end

%%
if plotOn==1
    cFigure;
    subplot(1,2,1); hold on;
    title(['Extend method: ',num2str(extendMethod)]);
    hp1(1)=gpatch(F,V,'kw','none',0.5);
    hp1(2)=gpatch(Fn,Vn,'bw','b',0.5,1);
    hp1(3)=quiverVec(V1,n,pointSpacing,'r');
    hp1(4)=quiverVec(V1,ne,pointSpacing,'g');
    hp1(5)=quiverVec(V1,ns,pointSpacing,'b');
    hp1(6)=quiverVec(mean(V1,1),ns_mean,extendDistance+5*pointSpacing,'k');
    legend(hp1,{'Input surface','Extended surface','Vertex surface normal',...
               'vertex edge vector','vertex extrude vector','Global or mean direction'},...
               'Location','SouthOutSide');
    axisGeom; camlight headlight;

    subplot(1,2,2); hold on;
    hp2(1)=gpatch(F,V,'kw','k',1);
    hp2(2)=gpatch(Fn,Vn,'bw','k',1);
    legend(hp2,{'Input surface','Extended surface'},'Location','SouthOutSide');
    axisGeom; camlight headlight;

    gdrawnow;
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% SP_H1_ERROR: Evaluate the error in H^1 norm.
%
%   [errh1, errl2, errh1s] = sp_h1_error (space, msh, u, uex, graduex)
%
% INPUT:
%
%    space:   object defining the space of discrete functions (see sp_multipatch)
%    msh:     object defining the domain partition and the quadrature rule (see msh_multipatch)
%    u:       vector of dof weights
%    uex:     function handle to evaluate the exact solution
%    graduex: function handle to evaluate the gradient of the exact solution
%
% OUTPUT:
%
%     errh1:  error in H^1 norm
%     errl2:  error in L^2 norm
%     errh1s: error in H^1 seminorm
%
% Copyright (C) 2015 Rafael Vazquez
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

function [errh1, errl2, errh1s] = sp_h1_error (space, msh, u, uex, graduex)

  if (space.npatch ~= msh.npatch)
    error ('The number of patches does not coincide')
  end

  if (numel(u) ~= space.ndof)
    error ('Wrong size of the vector of degrees of freedom')
  end

  for iptc = 1:msh.npatch
    if (isempty (space.dofs_ornt))
      u_ptc = u(space.gnum{iptc});
    else
      u_ptc = u(space.gnum{iptc}) .* space.dofs_ornt{iptc}.';
    end
    [error_h1(iptc), error_l2(iptc), error_h1s(iptc)] = ...
        sp_h1_error (space.sp_patch{iptc}, msh.msh_patch{iptc}, u_ptc, uex, graduex);
  end
  errl2 = sqrt (sum (error_l2 .* error_l2));
  errh1 = sqrt (sum (error_h1 .* error_h1));
  errh1s = sqrt (sum (error_h1s .* error_h1s));
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function r4_uniform_01_test ( )

%*****************************************************************************80
%
%% R4_UNIFORM_01_TEST tests R4_UNIFORM_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 July 2006
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R4_UNIFORM_01_TEST\n' );
  fprintf ( 1, '  R4_UNIFORM_01 computes pseudorandom values\n' );
  fprintf ( 1, '  in the interval [0,1].\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The initial seed is %d\n', seed );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ x, seed ] = r4_uniform_01 ( seed );
    fprintf ( 1, '  %6d  %14f\n', i, x );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function val = v_dtiInit
% Validate dti raw preprocessing of GE data
%
%    val = v_dtiInit()
%
% This script takes a very long time to run and should not be used as part
% of the usual unit testing.  LM Perry is thinking about shortening this.
% And we should be writing other test_dtiXXXX scripts for unit testing.
%
% This function checks that the average FA and B0 values resulting from
% dtiRawPreprocess are consistent with the expected value on different
% platforms (LINUX, WINDOWS, etc.).
%
% Requires vistadata and vistasoft on your path. The directory
% fullfile(mrvDataRootPath,'validate','dwi') contains all the necessary
% data files.
%
% This function does the following:
%
% (1) run dtiInit on the files in the raw folder. This will align to the T1
%     and create a B0, BVECS and FA values (it creates more stuff but we
%     only focus on these for the moment).
%
% (2) Show a montage of Alignment i.e., the T1 and DTI overalyed. Check for
%     LR flips and correct alignment.
%
% (3) Compute FA, Mean diffusivity, radial diffusivity nd Axial diffusivity
%     across the brain and check the value obtained on different
%     platforms. THis will be done using: [fa,md,rd,ad] = dtiComputeFA(eigVal)
%
% (5) Load the stored FA, MD, RD, AD for the whole brain in:
%     'dwi/storedValues.mat'
%
% (6) Recompute them from the data
%
% (7) Compute the difference between the stored and the recomputed ones.
%
% (8) Remove the data that was created by this code if there is not a
%     significant difference in the computed values.
%
% Example:
%  This function was meant to be called by mrvValidate, which will
%  compute the difference between the stored values and the re-computed
%  values.
%
%   mrvValidate([],[], 'test_dtiInit');
%
% See also: mrvValidateAll.m
%
% Adapted from v_dtiRawPreprocessGE
%
% (C) Stanford VISTA, 2011
%
%
%%
% This function takes a while to run. We want validation functions to be
% fast. Here we warn the user that this test will take a while and confirm
% that they want to continue.

prompt = 'You are about to process DWI data - this will take ~1hr to run. Do you want to continue?';
resp = questdlg(prompt,'Confirm','Yes','Cancel','Cancel');
if(strcmp(resp,'Cancel'))
    disp('canceled.');
    assertEqual(0,1,'NOTICE: test_dtiInit.m was canceled by the user.');
    return;
end

%% Get the data pathdata path
dataDir = fullfile(mrvDataRootPath,'validate','dwi');
rawDwi  = fullfile(mrvDataRootPath,'validate','dwi','raw','dti.nii.gz');
t1      = fullfile(mrvDataRootPath,'validate','dwi','t1.nii.gz');

%% Retain original directory, change to data directory
curDir = pwd;
cd(dataDir);

%% Run Preprocess
% setting clobber to 'always' = 1, so that output files will be silently
% replaced.
params  = dtiInitParams('phaseEncodeDir',2,'clobber',1);
dt6File = dtiInit(rawDwi,t1,params);

%% Show alignment to t1
% This is automatically computed by dtiRawPreprocess.m
t1pdd = fullfile(mrvDirup(dt6File),'t1pdd.png');
imshow(imread(t1pdd),'InitialMagnification','fit');

%% Compute mean FA, RD, MD, AD values and check them with the stored one.

% Load the stored values
ref = load(fullfile(dataDir,'test_dtiInit.mat'));

% Load the dt6 file.
dt = dtiLoadDt6(dt6File);

% Extract the eigen values.
eigVal = dt.dt6;

% Compute the fractional, mean, radial and axial diffusivity
[fa,md,rd,ad] = dtiComputeFA(eigVal);

% Compute the mean across the whole brain
val.fa = nanmean(fa(:));
val.md = nanmean(md(:));
val.rd = nanmean(rd(:));
val.ad = nanmean(ad(:));

%% Test the results

% Check the mean FA
assertEqual(ref.val.fa, val.fa);
% Check the mean MD
assertEqual(ref.val.md, val.md);
% Check the mean RD
assertEqual(ref.val.rd, val.rd);
% Check the mean AD
assertEqual(ref.val.ad, val.ad);

%% Go back to the original directory, done!
%
cd(curDir)

return

%% This is to save a new mat file with the values.
% save('test_dtiInit.mat','val');

%% old code
% mrvValidate will use the output of this function to compare the
% recomputed values to the stored values.
% fr example doing something like this:
%
% check the computed with the stored values
% val.faErr = diff([mean_fa,meanVals.fa]);
% val.mdErr = diff([mean_md,meanVals.md]);
% val.rdErr = diff([mean_rd,meanVals.rd]);
% val.adErr = diff([mean_ad,meanVals.ad]);
%
% show results on matlab output
% errFields = fields(val);
% meanFields = fields(meanVals);
% for i = 1:length(fields(val))
%  fprintf('[%s] Error in ''%s'': %2.8f\n',mfilename, meanFields{i}, val.(errFields{i}));
% end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function x = vec( x )

% VEC   CVX implementation of vec

s = x.size_;
n = prod(s);
if s(1) ~= n,
	x = cvx( [ n, 1 ], x.basis_ );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
%kLine 'Line Object Menuform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros Line.pane file
%
% Parameters:
%
% Example: kLine( {})
%
% Khoros helpfile follows below:
% .begin tagged
%
% .item "The Line Visual Object"
% A line visual object supports the display of a line annotation.  The endpoints
% of the line are specified in world coordinates.  The  world coordinates are
% specified by the parent.  By default, these world coordinates are from 0 to 1
% unless otherwise specified.  So, for example, if the line is created within
% a manager object, its world coordinates are from 0 to 1.  However, if the line
% is created over an image, the world coordinates will be dictated by the size of
% the image; if it created over a plot, the world coordinates will be dictated
% by the plot.
%
% .item "Line Width"
% This list selection allows you to set the line width used by
% the line object. Line widths range from "Extra Fine" to "Extra Wide".
%
% .item "Line Type"
% This list selection lets you set the line type used on
% the line object. There are seven line types, including "Solid",
% "Dotted", "Dot Dash", "Short Dash", "Long Dash", "Odd Dash", and
% "Grid Dotted".
%
% .item "Line Color"
% This stringlist selection lets you set the foreground color
% of the line object, using the color name.
%
% .item "X Begin"
% This double value specifies the x location of the beginning of
% the line segment in world coordinates.
%
% .item "X End"
% This double value specifies the x location of the end of
% the line segment in world coordinates.
%
% .item "Y Begin"
% This double value specifies the y location of the beginning of
% the line segment in world coordinates.
%
% .item "Y End"
% This double value specifies the y location of the end of
% the line segment in world coordinates.
%
% .end tagged

function varargout = kLine(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kLine(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={};
maxval={};
minval={};
istoggle=[];
was_set=istoggle * 0;
paramtype={};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; 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(0);
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
callKhoros([w 'Line" '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function s = yj_to_s_common ( y, j )

%*****************************************************************************80
%
%% YJ_TO_S_COMMON writes a Common YJ date into a string.
%
%  Format:
%
%    CE YYYY/JJJ
%    BCE YYYY/JJJ
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, J, the YJ date.
%
%    Output, string S, the representation of the date.
%
  if ( y < 0 )
    s = sprintf ( 'BCE %d/%d', -y, d );
  else
    s = sprintf ( 'CE %d/%d', y, d );
  end

  return
end


text:
function [ r, seed ] = r4_uniform_01 ( seed )

%*****************************************************************************80
%
%% R4_UNIFORM_01 returns a unit pseudorandom R4.
%
%  Discussion:
%
%    This routine implements the recursion
%
%      seed = 16807 * seed mod ( 2**31 - 1 )
%      r = seed / ( 2**31 - 1 )
%
%    The integer arithmetic never requires more than 32 bits,
%    including a sign bit.
%
%    If the initial seed is 12345, then the first three computations are
%
%      Input     Output      R4_UNIFORM_01
%      SEED      SEED
%
%         12345   207482415  0.096616
%     207482415  1790989824  0.833995
%    1790989824  2035175616  0.947702
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley Interscience, page 95, 1998.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer SEED, the integer "seed" used to generate
%    the output random number.  SEED should not be 0.
%
%    Output, real R, a random value between 0 and 1.
%
%    Output, integer SEED, the updated seed.  This would
%    normally be used as the input seed on the next call.
%
  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R4_UNIFORM_01 - Fatal error!' );
  end

  seed = floor ( seed );

  seed = mod ( seed, 2147483647 );

  if ( seed < 0 )
    seed = seed + 2147483647;
  end

  k = floor ( seed / 127773 );

  seed = 16807 * ( seed - k * 127773 ) - k * 2836;

  if ( seed < 0 )
    seed = seed + 2147483647;
  end

  r = seed * 4.656612875E-10;

  return
end


text:
function call_price=american_call_futures_currcy_bin(S, K, r, r_f, sigma, time, no_steps)

%--------------------------------------------------------------------------
%
% DESCRIPTION:
%
% Pricing a futures currency option using a binomial approximation
%
%
% Reference:
%
% John Hull, "Options, Futures and other Derivative Securities",
% Prentice-Hall, second edition, 1993.
%
%--------------------------------------------------------------------------
%
% INPUTS:
%
%  S:           spot price
%  K:           exercice price
%  r:           domestic interest rate
%  r_f:         foreign interest rate
%  sigma:       volatility
%  time:        time to maturity
%  no_steps:    number of steps in binomial tree
%
%--------------------------------------------------------------------------
%
% OUTPUT:
%
% call_price: price of a call option
%
%--------------------------------------------------------------------------
%
% Author:  Paolo Z., February 2012
%
%--------------------------------------------------------------------------

exchange_rates=zeros(no_steps+1,1);

call_values=zeros(no_steps+1,1);

t_delta= time/no_steps;
Rinv = exp(-r*(t_delta));
u = exp(sigma*sqrt(t_delta));
d = 1.0/u;
uu= u*u;
pUp = (exp((r-r_f)*t_delta)-d)/(u-d);
pDown = 1.0 - pUp;
exchange_rates(1) = S*(d^no_steps);

for ( i=2:(no_steps+1) )
    exchange_rates(i) = uu*exchange_rates(i-1);
end

for ( i=1:(no_steps+1) )
    call_values(i) = max(0.0, (exchange_rates(i)-K));
end

for ( step=no_steps:-1:1 )
    for ( i=1:1:(step) )
        exchange_rates(i) = d*exchange_rates(i+1);
        call_values(i) = (pDown*call_values(i)+pUp*call_values(i+1))*Rinv;
        call_values(i) = max(call_values(i), exchange_rates(i)-K);
    end
end

call_price = call_values(1);

author: Sable

repo: mcbench-benchmarks


text:
% Program for Bus Power Injections, Line & Power flows (p.u)...

%function [Pi Qi Pg Qg Pl Ql] = loadflow(nb,V,del,BMva)
function [Load_Flow_M Line_Flow_M] = loadflow(nb,V,del,BMva,Line_Data,Bus_Data)

Y = ybusppg(nb,Line_Data);                % Calling Ybus program..
%lined = linedatas(nb);          % Get linedats.
lined = Line_Data;
busd = Bus_Data;
%busd = busdatas(nb);            % Get busdatas..

Vm = pol2rect(V,del);           % Converting polar to rectangular..
Del = 180/pi*del;               % Bus Voltage Angles in Degree...
fb = lined(:,1);                % From bus number...
tb = lined(:,2);                % To bus number...
nl = length(fb);                % No. of Branches..
Pl = busd(:,7);                 % PLi..
Ql = busd(:,8);                 % QLi..

Iij = zeros(nb,nb);
Sij = zeros(nb,nb);
Si = zeros(nb,1);

% Bus Current Injections..
 I = Y*Vm;
 Im = abs(I);
 Ia = angle(I);

%Line Current Flows..
for m = 1:nl
    p = fb(m); q = tb(m);
    Iij(p,q) = -(Vm(p) - Vm(q))*Y(p,q); % Y(m,n) = -y(m,n)..
    Iij(q,p) = -Iij(p,q);
end
Iij = sparse(Iij);
Iijm = abs(Iij);
Iija = angle(Iij);

% Line Power Flows..
for m = 1:nb
    for n = 1:nb
        if m ~= n
            Sij(m,n) = Vm(m)*conj(Iij(m,n))*BMva;
        end
    end
end
Sij = sparse(Sij);
Pij = real(Sij);
Qij = imag(Sij);

% Line Losses..
Lij = zeros(nl,1);
for m = 1:nl
    p = fb(m); q = tb(m);
    Lij(m) = Sij(p,q) + Sij(q,p);
end
Lpij = real(Lij);
Lqij = imag(Lij);

% Bus Power Injections..
for i = 1:nb
    for k = 1:nb
        Si(i) = Si(i) + conj(Vm(i))* Vm(k)*Y(i,k)*BMva;
    end
end
Pi = real(Si);
Qi = -imag(Si);
Pg = Pi+Pl;
Qg = Qi+Ql;

Load_Flow_M = zeros(nb+1,9);

for m = 1:nb
    Load_Flow_M(m,1) = m;
    Load_Flow_M(m,2) = V(m);
    Load_Flow_M(m,3) = Del(m);
    Load_Flow_M(m,4) = Pi(m);
    Load_Flow_M(m,5) = Qi(m);
    Load_Flow_M(m,6) = Pg(m);
    Load_Flow_M(m,7) = Qg(m);
    Load_Flow_M(m,8) = Pl(m);
    Load_Flow_M(m,9) = Ql(m);
end
m = m + 1;
    Load_Flow_M(m,1) = NaN;
    Load_Flow_M(m,2) = NaN;
    Load_Flow_M(m,3) = NaN;
    Load_Flow_M(m,4) = sum(Pi);
    Load_Flow_M(m,5) = sum(Qi);
    Load_Flow_M(m,6) = sum(Pi+Pl);
    Load_Flow_M(m,7) = sum(Qi+Ql);
    Load_Flow_M(m,8) = sum(Pl);
    Load_Flow_M(m,9) = sum(Ql);

Line_Flow_M = zeros(nl+1,10);

for m = 1:nl
    p = fb(m); q = tb(m);
    Line_Flow_M(m,1) = full(p);
    Line_Flow_M(m,2) = full(q);
    Line_Flow_M(m,3) = full(Pij(p,q));
    Line_Flow_M(m,4) = full(Qij(p,q));
    Line_Flow_M(m,5) = full(q);
    Line_Flow_M(m,6) = full(p);
    Line_Flow_M(m,7) = full(Pij(q,p));
    Line_Flow_M(m,8) = full(Qij(q,p));
    Line_Flow_M(m,9) = Lpij(m);
    Line_Flow_M(m,10) = Lqij(m);
end
m = m+1;
    Line_Flow_M(m,1) = NaN;
    Line_Flow_M(m,2) = NaN;
    Line_Flow_M(m,3) = NaN;
    Line_Flow_M(m,4) = NaN;
    Line_Flow_M(m,5) = NaN;
    Line_Flow_M(m,6) = NaN;
    Line_Flow_M(m,7) = NaN;
    Line_Flow_M(m,8) = NaN;
    Line_Flow_M(m,9) = sum(Lpij);
    Line_Flow_M(m,10) = sum(Lqij);


text:
function varargout = coeffs2( f, m, n )
%COEFFS2   Double Fourier coefficients of a SPHEREFUN.
%
%   X = COEFFS2( F ) returns the 2D Fourier modes of the spherefun, viewed
%   as a doubly periodic function.
%
%   [C, D, R] = COEFFS2( F ) returns a low rank approximation to the 2D
%   Fourier modes.
%
%   X = COEFFS2(F, M, N) returns bivariate coefficients with N Fourier
%   modes in latitude and M Fourier modes in longitude.

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

% empty check
if isempty(f)
    varargout = {};
    return
end

% Calculate the CDR decomposition:
[C, D, R] = cdr(f);

if ( nargin == 1 )
    % Find the  coefficients of each slice:
    U = C.coeffs;
    R = R.coeffs;
else
    if ( nargin == 2 )
        n = m;
    end
    % Find the coefficients of each slice:
    U = trigtech.alias(C.coeffs, n);
    R = trigtech.alias(R.coeffs, m);
end

% Prepare the output. Keep in low rank form if nargin > 1.
if ( nargout <= 1 )
    varargout = { U*D*R.' };
elseif ( nargout <= 3 )
    varargout = { U, D, R };
else
    error('SPHEREFUN:COEFFS:NARGOUT',...
            'Too many output arguments')
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function a = ge( x, y )

%Disciplined convex programming information for GE (>=):
%   The right-hand side of a less-than constraint must be convex. The
%   left-hand side must be concave. Of course, real constant and affine
%   expressions are both convex and concave and can be used on either
%   side as well.
%
%Disciplined geometric programming information for GE (>=):å
%   The right-hand side of a less-than constraint must be log-convex---
%   including positive constants, monomials, posynomials, generalized
%   posynomials, and products thereof. The left-hand side must be log-
%   concave---including positive constants, monomials, reciprocals of
%   log-convex expressions, and products thereof.
%
%Note that CVX does not distinguish between strict greater-than (>) and
%greater-than-or-equal (<=) constraints; they are treated identically.
%Feasible interior-point solvers tend to return points which satisfy
%strict inequality, but not all solvers do.

evalin( 'caller', 'cvx_verify' );
b = cvx_pushcnstr( x, y, '>=' );
if nargout, a = b; end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
%install;

load('..\mat\LPnetlib\lp_agg.mat')

A = Problem.A;
p = colamd(A');
A = A(p,:);
w = (0.5+rand(size(A,2),1)).^5;
R = chol((A*(w.*A')));
x = rand(size(A,1),1);
z = R\(R'\x);
d = full(diag(R));
w = diag(sparse(w));

o = AdaptiveChol(A, 1e-8);
acc = o.factorize(w)
z2 = o.solve(ddouble(x), ddouble(w));
z3 = o.solve(ddouble(x), ddouble(w), 2);
z4 = o.solve(ddouble(x), ddouble(w), 3);
z5 = o.solve(ddouble(x), ddouble(w), 4);
d2 = o.diagonal();

norm(z2-z)
norm(z2-z3)
norm(z3-z4)
norm(z4-z5)

%norm(d2-d)
ls = o.leverageScore(100);
max(ls)

author: opencobra

repo: cobratoolbox


text:
function val=BesselI1RatioApprox(kappa,algorithm)
%%BESSELI1RATIOAPPROX Approximate the ratio of modified Bessel functions of
%              the first kind of the form x=I_{1}(kappa)/I_{0}(kappa).
%              For a more exact solution, use the function BesseliRatio,
%              which can also handle different subscripts (orders) for the
%              modified Bessel functions of the first kind.
%
%INPUTS: kappa The real argument of the modified Bessel function of the
%              first kind kappa>=0.
%    algorithm A parameter selecting the approximation that is used.
%              Possible values are:
%              0 (The default if omitted or an empty matrix is passed) Use
%                the approximation used in Equations 8 and 10 of [1].
%                Equation 8 in [1] is the first two terms of an expansion
%                given on page 290 of [2].
%              1 Use the approximations in Equations 5 and 7 of [3].
%                Equation 5 of [3] is related to a Taylor series expansion
%                given on page 289 of [2].
%
%OUTPUTS: x An approximation to the ratio I_{1}(kappa)/I_{0}(kappa).
%
%REFERENCES:
%[1] G. Stienne, S. Reboul, M. Azmani, J. B. Choquel, and M. Benjelloun, A
%    multi-temporal multi-sensor circular fusion filter," Information
%    Fusion, vol. 18, pp. 86-100, Jul. 2014.
%[2] S. R. Jammalamadaka and A. SenGupta, Topics in Circular Statistics.
%    Singapore: World Scientific, 2001.
%[3] G. Stienne, S. Reboul, J. B. Choquel, and M. Benjelloun, "Circular
%    data processing tools applied to a phase open loop architecture for
%    multi-channel signals tracking," in Position Location and Navigation
%    Symposium, Myrtle Beach, SC, 23-26 Apr. 2012, pp. 633-642.
%
%April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<2||isempty(algorithm))
        algorithm=0;
    end

    switch(algorithm)
        case 0%The method of [1].
            if(kappa>=0.6)
                val=1-1/(2*kappa);
            else
                val=exp(-1/(2*kappa));
            end
        case 1
            if(kappa>=0.6)
                val=(1-3/(8*kappa)-15/(128*kappa^2))/(1+1/(8*kappa)+9/(128*kappa^2));
            else
                 val=kappa/2;
            end
        otherwise
            error('Unknown Algorithm Specified.')
    end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function SO3VF = rotate(SO3VF,rot,varargin)
% rotate a SO3 vector field by a rotation
%
% Syntax
%   SO3VF = rotate(SO3VF,rot)
%
% Input
%  SO3VF - @SO3VectorField
%  rot  - @rotation
%
% Output
%  SO3VF - @SO3VectorFieldHandle
%

SO3VF = SO3VectorFieldHandle(@(q) SO3VF.eval(inv(rot)*q),SO3VF.CS,SO3VF.SS);

author: mtex-toolbox

repo: mtex

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


text:
function xBnd = SmoothBnd(x,alpha,Bnd)
%
% xBnd = SmoothBnd(x,alpha,Bnd)
%
% This function is used to smoothly bound the input x;
%
% INPUTS:
%   x = a vector of real numbers to apply bounds
%   alpha = a positive smoothing parameter. Small values correspond to
%           little smoothing
%   Bnd = the desired bounds of the function, expressed as a 2-element row
%         vector. If ommitted it will default to [0,1];
%
% Written by Matthew Kelly
% October 2013
% Cornell University
%

if nargin==2
    Low = 0;
    Upp = 1;
else
    Low = Bnd(1);
    Upp = Bnd(2);
end

infTest1 = exp(max(x-Low)/alpha);
infTest2 = exp(max(-(x-Upp))/alpha);

if isinf(infTest1) || isinf(infTest2)  %Then there is a sharp transition
    xBnd = x;
    xBnd(x<Low) = Low;
    xBnd(x>Upp) = Upp;
else
    %Apply bounding to each part of the input:
    xLow = Low + alpha*log(exp((x-Low)/alpha)+1);
    xUpp = Upp - alpha*log(exp(-(x-Upp)/alpha)+1);

    %Combine:
    xBnd = xLow + xUpp - x;
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function h = image_histogram( im, mask, make_plot )
% IMAGE_HISTOGRAM  Calculate the histogram of the given image
%
% h = image_histogram( im, mask )
%
% calculate the histogram of the given image, optionally restricted on the mask
% h(1) is the number of pixels of value 0, etc.

% F. Nedelec, Dec. 2007
%
% Ramon Casero <rcasero@gmail.com>: Minor edits

%%compatibility with tiffread:
if ( isfield(im,'data') )
    im = double( im.data );
end

if nargin < 3
    make_plot = 0;
end

%%
if nargin < 2 || isempty(mask)

    max_val = ceil(max(reshape( im, numel(im), 1 )));

    h = histc( reshape( im, numel(im), 1 ), 0:max_val);
    %h = sum( histc( im, 0:max_val ), 2);

else

    if any( size(im) ~= size(mask) )
        error('Image and mask must be have the same size');
    end
    if min(reshape(mask, numel(mask), 1 )) < 0
        error('Values of the mask should be non-negative');
    end
    if max(reshape(mask, numel(mask), 1 )) > 1
        error('Values of the mask should be lower or equal to 1');
    end

    imm = im .* mask - ( 1-mask );
    val = reshape( imm, numel(imm), 1 );
    max_val = ceil(max(val));

    h = histc(val, 0:max_val);

end

%% Make a figure to display the histogram
if make_plot

    x = (0:size(h)-1)';
    figure('Name',inputname(1), 'Position', [100 150 800 300]);
    axes('Position', [0.05 0.1 0.9 0.8] );
    xlim([0 size(h,1)]);
    plot(x, h, 'g.' );
    title('Histogram of pixel values');

end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function volROI = flat2volROI(flatROI,flatView,volView)
%
% volROI = flat2volROI(flatROI,flatView,volView)
%
% Creates a volume ROI from a flat ROI by looking up the
% corresponding coords.  The tricky thing about this function
% is that FLAT -> VOL is a one to many mapping, and we need to
% find all the corresponding volume coords.
%
% flatROI and volROI are ROI structures, like those found in
% view.ROIs
%
% volView must be the VOLUME structure.
% flatView must be the FLAT structure.
%
% djh, 8/98.
%
% djh, 8/4/99.  Need to round the coords because we no longer do it in
%               loadGLocs.  Also, remove duplicate flat coords.
% djh, 2/2001. Use intersectCols instead of coords2Indices

global mrSESSION

% check that the ROI fields are properly set
flatROI = roiCheck(flatROI);

% This should be much fastter but its busted because intersect removes duplicates and
% we need to keep the duplicates here because flat->gray is a one-to-many mapping.
%
% coords = [];
% for h = 1:2
%     ROIcoords = flatROI.coords([1:2],find(flatROI.coords(3,:) == h));
%     if ~isempty(ROIcoords)
%         [tmp,ROIIndices,flatIndices] = intersectCols(ROIcoords,round(flatView.coords{h}));
%         coords = [coords, flatView.grayCoords{h}(:,flatIndices)];
%     end
% end

coords = [];
flatImSize = flatView.ui.imSize;
for h = 1:2
  ROIcoords = flatROI.coords([1:2],find(flatROI.coords(3,:) == h));
  if ~isempty(ROIcoords)
    ROIIndices = coords2Indices(ROIcoords,flatImSize);
    flatIndices = coords2Indices(round(flatView.coords{h}),flatImSize);
    bothIndices = intersect(ROIIndices,flatIndices);
    for id = 1:length(bothIndices)
      indices = find(flatIndices == bothIndices(id));
      coords = [coords, flatView.grayCoords{h}(:,indices)];
    end
  end
end

% Remove duplicates (there shouldn't be any, but just to make sure).
coords = intersectCols(coords,coords);

% Set the fields
volROI = flatROI;
volROI.coords = coords;
volROI.viewType = volView.viewType;

return


text:
function [b,a]=v_potsband(fs)
%V_POTSBAND Design filter for 300-3400 telephone bandwidth [B,A]=(FS)
%
%Input: FS=sample frequency in Hz
%
%Output: B/A is a discrete time bandpass filter with a passband gain of 1
%
%The filter meets the specifications of G.151 for any sample frequency
%and has a gain of -3dB at the passband edges.

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

szp=[0.19892796195357i; -0.48623571568937+0.86535995266875i];
szp=[[0; -0.97247143137874] szp conj(szp)];
% s-plane zeros and poles of high pass 3'rd order chebychev2 filter with -3dB at w=1
zl=2./(1-szp*tan(300*pi/fs))-1;
al=real(poly(zl(2,:)));
bl=real(poly(zl(1,:)));
sw=[1;-1;1;-1];
bl=bl*(al*sw)/(bl*sw);
zh=2./(szp/tan(3400*pi/fs)-1)+1;
ah=real(poly(zh(2,:)));
bh=real(poly(zh(1,:)));
bh=bh*sum(ah)/sum(bh);
b=conv(bh,bl);
a=conv(ah,al);

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function writeParaviewWave(waves, t, numPointsX, numPointsY, domainSize, model, simdate, mooring, paraviewPath,TimeBodyParav,g)
% Method to write ``vtp`` Paraview visualization files for the waveClass.
% Executed by paraviewVisualization.m when simu.paraview.option=1 in the
% wecSimInputFile.m
%
% Parameters
% ------------
%   waves : waveClass
%       Instance of the waveClass that is being written to Paraview files.
%   t : float vector
%       Wave time vector
%   numPointsX : int
%       Number of points for visualization in the x direction
%   numPointsY : int
%       Number of points for visualization in the y direction
%   domainSize : int
%       Length of the wave field in meters
%   model : string
%       The simMechanics ``.slx`` filename
%   simdate : string
%       Date and time of the simulation
%   mooring : int
%       MoorDyn flag
%   paraviewPath : directory
%       Directory the Paraview files were saved
%   TimeBodyParav : float vector
%       Paraview time vector
%   g : float
%       Gravitational acceleration constant
%

% ground plane
filename = [paraviewPath, filesep 'ground.txt'];
fid = fopen(filename, 'w');
fprintf(fid,[num2str(domainSize) '\n']);
fprintf(fid,[num2str(waves.waterDepth) '\n']);
fprintf(fid,[num2str(mooring) '\n']);
fclose(fid);
% wave
x = linspace(-domainSize, domainSize, numPointsX);
y = linspace(-domainSize, domainSize, numPointsY);
[X,Y] = meshgrid(x,y);
lx = length(x);
ly = length(y);
numVertex = lx * ly;
numFace = (lx-1) * (ly-1);
for it = 1:length(t)
    % open file
    filename = [paraviewPath, filesep 'waves' filesep 'waves_' num2str(it) '.vtp'];
    fid = fopen(filename, 'w');
    % calculate wave elevation
    Z = waveElevationGrid (waves, t(it), X, Y, it, g);                % write header
    fprintf(fid, '<?xml version="1.0"?>\n');
    fprintf(fid, ['<!-- WEC-Sim Visualization using ParaView -->\n']);
    fprintf(fid, ['<!--   model: ' model ' - ran on ' simdate ' -->\n']);
    fprintf(fid, ['<!--   wave:  ' waves.type ' -->\n']);
    fprintf(fid, ['<!--   time:  ' num2str(TimeBodyParav(it)) ' -->\n']);
    fprintf(fid, '<VTKFile type="PolyData" version="0.1">\n');
    fprintf(fid, '  <PolyData>\n');
    % write wave info
    fprintf(fid,['    <Piece NumberOfPoints="' num2str(numVertex) '" NumberOfPolys="' num2str(numFace) '">\n']);
    % write points
    fprintf(fid,'      <Points>\n');
    fprintf(fid,'        <DataArray type="Float32" NumberOfComponents="3" format="ascii">\n');
    for jj = 1:length(y)
        for ii = 1:length(x)
            pt = [X(jj,ii), Y(jj,ii), Z(jj,ii)];
            fprintf(fid, '          %5.5f %5.5f %5.5f\n', pt);
        end; clear ii
        clear pt
    end; clear jj
    fprintf(fid,'        </DataArray>\n');
    fprintf(fid,'      </Points>\n');
    % write squares connectivity
    fprintf(fid,'      <Polys>\n');
    fprintf(fid,'        <DataArray type="Int32" Name="connectivity" format="ascii">\n');
    for jj = 1:ly-1
        for ii = 1:lx-1
            p1 = (jj-1)*lx + (ii-1);
            p2 = p1+1;
            p3 = p2 + lx;
            p4 = p1 + lx;
            fprintf(fid, '          %i %i %i %i\n', [p1,p2,p3,p4]);
        end; clear ii
    end; clear jj
    fprintf(fid,'        </DataArray>\n');
    fprintf(fid,'        <DataArray type="Int32" Name="offsets" format="ascii">\n');
    fprintf(fid, '         ');
    for ii = 1:numFace
        n = ii * 4;
        fprintf(fid, ' %i', n);
    end; clear ii n
    fprintf(fid, '\n');
    fprintf(fid,'        </DataArray>\n');
    fprintf(fid, '      </Polys>\n');
    % end file
    fprintf(fid, '    </Piece>\n');
    fprintf(fid, '  </PolyData>\n');
    fprintf(fid, '</VTKFile>');
    % close file
    fclose(fid);
end; clear it
clear  numPoints numVertex numFace x y lx ly X Y Z fid filename p1 p2 p3 p4
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function p11_default ( dim_num )

%*****************************************************************************80
%
%% P11_DEFAULT sets default values for problem 11.
%
%  Discussion:
%
%    If a problem uses vector parameters, then the spatial dimension
%    DIM_NUM is needed as input, so that the vector parameters can be
%    properly dimensioned.
%
%    Every problem keeps a count of the number of function calls.  Calling
%    this routine causes that count to be reset to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
  i4 = 0;
  p11_i4 ( 'D', '*', i4 );

  return
end


text:
%kTimer 'Timer Object Menuform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros Timer.pane file
%
% Parameters:
%
% Example: kTimer( {})
%
% Khoros helpfile follows below:
% .begin tagged
%
% .item "The Timer Visual Object"
% The timer object is provided for use as a stopwatch.  It displays the current
% time (including seconds) according to the system clock of the computer.  It can
% be started and stopped as desired, and the update time can be specified
% to fractions of a second.  The timer object is a specialized annotation
% subclassed from the string visual object.  The (x, y) location of the lower
% left corner of the timer timer may be specified in device coordinates or
% world coordinates by the application.  The menuform, however, displays the
% world coordinate location of the timer, regardless of whether the application
% used world or device coordinates to specify that location.  Attributes of the
% timer object that can be interactively specified in addition to its world
% coordinate placement include its start time, and update time, in addition to
% its justification, style, font, and colors.
%
% .item "Update Time"
% This double selection specifies, in fractions of a second, how often
% the timer is updated.  For example, an update time of 5 would cause the timer
% to update every 5 seconds, while an update time of 0.25 would cause the timer
% to update 4 times every second.
%
% .item "Start Time"
% This double selection gives the start time of the timer.  By default it
% is time = 0.0, but it can be set to any number greater than zero.
%
% .item "Justification"
% This attribute specifies the justification of the string displaying the
% time "with respect to the width and height of the timer object itself,
% not its parent".
% Thus, if the timer object is just large enough to accomodate the time,
% justification will have little (if any) effect, as the time string has no
% room to "move" within its boundaries.  However, if the width and height are
% set to "larger" than what is necessary to display the time, the
% justification will take effect.  Justification may be set to "Center",
% "Top", "Bottom", "Left", "Right", "TopRight", "TopLeft", "BottomRight",
% or "BottomLeft".
%
% .item "Fill Background"
% When set to true, this attribute indicates that the background of the timer
% object should appear in the fill color (background color).
%
% .item "String Style"
% The date can be displayed in one of four styles: plain, emphasized, embossed
% in, or embossed out.  When the "plain" style is used, the date is drawn once,
% in the foreground color. When the date is "emphasized", it is double-drawn
% using both the foreground and highlight colors;  this causes it to have a bold,
% or emphasized, effect. When the "embossed in" style is used, the date will
% appear embossed "into" the background; both the foreground and highlight
% colors are used in embossing.  The "embossed out" style also causes the date
% to appear embossed, but the date will seem to come "out of" rather than
% "into" the background.
%
% .item "Font"
% Some of the fonts that are most commonly used with the date object are
% available in the stringlist menu;  or, you may provide your own in the text
% parameter box.
% .so  $DESIGN/manual/shared/xvobjects/general/FontName.sec
%
% .item "Foreground Color"
% This stringlist selection lets you set the foreground color of the timer object.
% The foreground color is always used to draw the time string,
% regardless of the "String Style".
%
% .item "Background Color"
% This stringlist selection lets you set the background color (fill color) of
% the timer object.  It is only used when "Fill Background" is set to true.
%
% .item "Highlight Color"
% This stringlist selection lets you set the highlight color; the highlight
% color is an additional color used in drawing the time string when the
% "String Style" is set to "Emphasized", "Embossed In", or "Embossed Out".
%
% .item "X Placement, Y Placement"
% This pair of double values specifies the (x, y) location of the lower left
% hand corner of the timer visual object in world coordinates.
%
% .end tagged

function varargout = kTimer(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kTimer(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={};
maxval={};
minval={};
istoggle=[];
was_set=istoggle * 0;
paramtype={};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; 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(0);
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
callKhoros([w 'Timer" '],Inputs,narglist);


text:
function varargout = edgesEvalDL( model, varargin )
% Run and evaluate structured edge detector with deep contour features on BSDS500.
% The code is modified from Structured Edge Detection Toolbox

% get default parameters
dfs={'dataType','test', 'name','', 'opts',{}, 'show',0, ...
  'pDistr',{{'type','parfor'}}, 'cleanup',0, 'thrs',99, 'maxDist',.0075};
p=getPrmDflt(varargin,dfs,1);

% load model and update model.opts acoording to opts
if( ischar(model) ), model=load(model); model=model.model; end
for i=1:length(p.opts)/2, model.opts.(p.opts{i*2-1})=p.opts{i*2}; end
rgbd=model.opts.rgbd; model.opts.nms=1;

% get list of relevant directories (image, depth, gt, results)
name = [model.opts.modelFnm,p.name];
imgDir = fullfile(model.opts.bsdsDir,'images',p.dataType);
depDir = fullfile(model.opts.bsdsDir,'depth',p.dataType);
resDir = fullfile(model.opts.modelDir,p.dataType,name);

gtDir  = fullfile(model.opts.bsdsDir,'groundTruth',p.dataType);
evaDir = resDir;

assert(exist(imgDir,'dir')==7); assert(exist(gtDir,'dir')==7);
% run edgesDetect() on every image in imgDir and store in resDir
if( ~exist(fullfile([evaDir '-eval'],'eval_bdry.txt'),'file') )
  if(~exist(resDir,'dir')), mkdir(resDir); end
  ids=dir(imgDir); ids=ids([ids.bytes]>0); ids={ids.name}; n=length(ids);
  ext=ids{1}(end-2:end); for i=1:n, ids{i}=ids{i}(1:end-4); end
  res=cell(1,n); for i=1:n, res{i}=fullfile(resDir,[ids{i} '.png']); end
  do=false(1,n); for i=1:n, do(i)=~exist(res{i},'file'); end
  ids=ids(do); res=res(do); m=length(ids);
  for i=1:m, id=ids{i};
    I = imread(fullfile(imgDir,[id '.' ext])); D=[];
    if(rgbd), D=single(imread(fullfile(depDir,[id '.png'])))/1e4; end
    if(rgbd==1), I=D; elseif(rgbd==2), I=cat(3,single(I)/255,D); end
    model.imgId = id;
    tic;
    E=edgesDetectDL(I,model); imwrite(uint8(E*255),res{i});
    toc;
  end
end

% perform actual evaluation using edgesEvalDir
varargout=cell(1,max(1,nargout));
[varargout{:}] = edgesEvalDir('resDir',resDir,'gtDir',gtDir,...
  'pDistr',p.pDistr,'cleanup',p.cleanup,'thrs',p.thrs,'maxDist',p.maxDist);
if( p.show ), figure(p.show); edgesEvalPlot(evaDir,name); end

end

repo: DeepContour


text:
function deviance = calc_deviance(trueStructNum,testStructNum,devMargin,planC)
% function deviance = calc_deviance(trueStructNum,testStructNum,devMargin,planC)
%
% Calculates the deviance of testStructNum from trueStructNum.
%
% Deviance = (volume of the true structure missed by the
% test structure + excess volume of the test structure over the true
% structure) / volume of the true structure.
%
% The true structure is 3-d contracted by the devMargin while computing the
% volume of the true structure missed by the test structure.
%
% The true structure is 3-d expanded by the devMargin while computing the
% excess volume of the test structure over the true structure.
%
% Example call:
% trueStructNum = 1; % structure index in planC
% testStructNum = 2; % structure index in planC
% devMargin = 0.2; % cm
% deviance = calc_deviance(trueStructNum,testStructNum,devMargin,planC);
%
% APA, 9/18/2017

if ~exist('planC','var')
    global planC
end
indexS = planC{end};

numStructs = length(planC{indexS.structures});

% Expand the test structure by amount equal to the devMargin
planC = createExpandedStructure(trueStructNum, devMargin, planC);
expandedStructNum = numStructs + 1;

% Contract the test structure by amount equal to the devMargin
planC = createExpandedStructure(trueStructNum, -devMargin, planC);
contracttructNum = numStructs + 2;

% Create a structure that's equal to the excess test volume
planC = createDifferenceStructure(testStructNum,expandedStructNum, planC);
excessTestStructNum = numStructs + 3;

% Create a structure that's equal to the excess true volume
planC = createDifferenceStructure(contracttructNum,testStructNum, planC);
excessTrueStructNum = numStructs + 4;

% Calculate the excess test and true volumes in cc
excessTestVol = getStructureVol(excessTestStructNum,planC);
excessTrueVol = getStructureVol(excessTrueStructNum,planC);

% delete the intermediate structures
for strToDelete = numStructs+4:-1:numStructs+1
    planC = deleteStructure(planC,strToDelete);
end

% Add them up
deviance = excessTestVol + excessTrueVol;

% Normalize by the true volume
trueVol = getStructureVol(trueStructNum,planC);

deviance = deviance / trueVol;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function Q = M2Q(M)
% Convert from rotation matrix to quaternion form
% See: http://skal.planet-d.net/demo/matrixfaq.htm
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: M2Q.m 7147 2017-08-03 14:07:01Z spm $

d = diag(M(1:3,1:3));
t = sum(d) + 1;
if t>0.5
    s = sqrt(t)*2;
    Q = [(M(3,2)-M(2,3))/s (M(1,3)-M(3,1))/s (M(2,1)-M(1,2))/s 0.25*s]';
else
    t = find(d==max(d));
    t = t(1);
    switch(t)
    case 1
        s = 2*sqrt(1 + M(1,1) - M(2,2) - M(3,3));
        Q = [0.25*s (M(1,2)+M(2,1))/s (M(3,1)+M(1,3))/s (M(3,2)-M(2,3))/s]';
    case 2
        s = 2*sqrt(1 + M(2,2) - M(1,1) - M(3,3));
        Q = [(M(1,2)+M(2,1))/s 0.25*s (M(2,3)+M(3,2))/s (M(1,3)-M(3,1))/s ]';
    case 3
        s = 2*sqrt(1 + M(3,3) - M(1,1) - M(2,2));
        Q = [(M(3,1)+M(1,3))/s (M(2,3)+M(3,2))/s 0.25*s (M(2,1)-M(1,2))/s]';
    end
end
if Q(4)<0, Q = -Q; end % w must be +ve

author: spm


text:
function varargout = gradientFD(f, dn, dim, deriv_order, accuracy_order)
%GRADIENTFD Calculate the gradient using a finite-difference method.
%
% DESCRIPTION:
%       gradientFD calculates the gradient of an n-dimensional input matrix
%       using the finite-difference method. For one-dimensional inputs, the
%       gradient is always computed along the non-singleton dimension. For
%       higher dimensional inputs, the gradient for singleton dimensions is
%       returned as 0. For elements in the centre of the grid, the gradient
%       is computed using centered finite-differences. For elements on the
%       edge of the grid, the gradient is computed using forward or
%       backward finite-differences. The order of accuracy of the
%       finite-difference approximation is controlled by accuracy_order
%       (default = 2). The calculations are done using sparse
%       multiplication, so the input matrix is always cast to double
%       precision.
%
% USAGE:
%       fx = gradientFD(f, dx)
%       fx = gradientFD(f, dx, [], deriv_order)
%       fx = gradientFD(f, dx, [], deriv_order, accuracy_order)
%       fn = gradientFD(f, dn, dim)
%       fn = gradientFD(f, dn, dim, deriv_order, accuracy_order)
%       [fx, fy] = gradientFD(f, dn)
%       [fx, fy] = gradientFD(f, dn, [], deriv_order, accuracy_order)
%       [fx, fy, fz, ...] = gradientFD(f, dn)
%       [fx, fy, fz, ...] = gradientFD(f, dn, [], deriv_order, accuracy_order)
%
% INPUTS:
%       f           - matrix or vector to find the gradient of
%       dn          - array of values for the grid point spacing in each
%                     dimension. If a value for dim is given, dn is the
%                     spacing in dimension dim.
%
% OPTIONAL INPUTS:
%       dim             - optional input to specify a single dimension over
%                         which to compute the gradient for n-dimension
%                         input functions
%       deriv_order     - order of the derivative to compute, e.g., use 1
%                         to compute df/dx, 2 to compute df^2/dx^2, etc.
%                         (default = 1)
%       accuracy_order  - order of accuracy for the finite difference
%                         coefficients. Because centered differences are
%                         used, this must be set to an integer multiple of
%                         2 (default = 2)
%
% OUTPUTS:
%       fx, fy, ... - gradient in the each dimension, where x corresponds
%                     to dim = 1, y corresponds to dim = 2 etc
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 24th August 2012
%       last update - 3rd September 2012
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also getFDMatrix, gradient, gradientSpect

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

% % display warning if input is not double precision
% if ~isa(f, 'double')
%     disp('gradientFD: converting input to double precision');
% end

% force input to be in double precision (sparse operations are only
% supported in double or logical formats)
f = double(f);

% get size of the input function
sz = size(f);

% check dimension of input
f_dims = numDim(f);

% check for user input for accuracy_order
if nargin < 5 || isempty(accuracy_order)
    accuracy_order = 2;
elseif rem(accuracy_order, 2)
    error('Input for accuracy_order must be an integer multiple of 2');
end

% check for user input for deriv_order
if nargin < 4 || isempty(deriv_order)
    deriv_order = 1;
elseif ~((deriv_order > 0) && (deriv_order == round(deriv_order)))
    error('Input for deriv_order must be an integer > 0');
end

% check for user input for dim
if (nargin < 3) || isempty(dim)

    % check if input is 1D
    if f_dims == 1
        % only compute gradient over required dimension
        if sz(1) == 1
            dim_array = 2;
        else
            dim_array = 1;
        end
    else
        % otherwise compute gradient over all dimensions
        dim = 0;
        dim_array = 1:f_dims;

        % check for the correct number of dn values
        if length(dn) ~= length(sz)
            error([num2str(length(sz)) ' values for dn must be specified for a ' num2str(length(sz)) '-dimensional input matrix']);
        end
    end
else
    dim_array = dim;
end

% only allow 1, 2, or 3D inputs (only these cases are implemented)
if f_dims > 3
    error('Input for f must have 1, 2, or 3 dimensions');
end

% if input is 1D, only allow it to be a row or column vector
if (f_dims == 1) && (length(sz) > 2)
    error('1D inputs for f must be row or column vectors');
end

% if input is 2D, only allow it to be a regular matrix
if (f_dims == 2) && (length(sz) > 2)
    error('2D inputs for f must be of size m x n');
end

% if input is 2D or 3D, check dim input isn't bigger than the matrix size
if (f_dims > 1) && (dim > f_dims)
    error('Input for dim cannot be greater than the number of dimensions of f');
end

% set output argument index
argout_index = 1;

% loop through required dimensions
for dim_index = dim_array

    % set dimension
    dim = dim_index;

    % check if a single dn value is given, if not, extract the required
    % value
    if numel(dn) ~= 1
        dn_val = dn(dim);
    else
        dn_val = dn;
    end

    % get the grid size along the specified dimension, or the longest
    % dimension if 1D
    if nargout == 1 && (max(sz) == prod(sz))
        [Nx, dim] = max(sz);
    else
        Nx = sz(dim);
    end

    % get the finite-difference matrix
    FDM = getFDMatrix(Nx, dn_val, deriv_order, accuracy_order);

    % compute derivatives of 1D or 2D matrix
    if ndims(f) < 3
        switch dim

            % derivative along dimension 1 (columns)
            case 1
                varargout{argout_index} = FDM * f;

            % derivative along dimension 2 (rows)
            case 2
                varargout{argout_index} = f * FDM.';

        end

    % compute derivatives of 3D matrix
    else
        switch dim

            % derivative along dimension 1
            case 1

                % preallocate output matrix
                varargout{argout_index} = zeros(size(f));

                % loop through z-plane
                for dim3_index = 1:sz(3)
                    varargout{argout_index}(:, :, dim3_index) = FDM * f(:, :, dim3_index);
                end

            % derivative along dimension 2
            case 2

                % preallocate output matrix
                varargout{argout_index} = zeros(size(f));

                % rotate FDM
                FDM = FDM.';

                % loop through z-plane
                for dim3_index = 1:sz(3)
                    varargout{argout_index}(:, :, dim3_index) = f(:, :, dim3_index) * FDM;
                end

            % derivative along dimension 3
            case 3

                % preallocate output matrix
                varargout{argout_index} = zeros(size(f));

                % rotate FDM
                FDM = FDM.';

                % loop through x-plane
                for dim1_index = 1:sz(1)
                    varargout{argout_index}(dim1_index, :, :) = squeeze(f(dim1_index, :, :)) * FDM;
                end

        end

    end

    % increment output argument index
    argout_index = argout_index + 1;

end

repo: Aero-matlab


text:
function patches = extract_patches(image, keypoints, patch_radius)
% EXTRACT_PATCHES - Extract local patches for all keypoints in an image.
%   image:
%       Gray-scale single image containing the keypoints.
%   keypoints:
%       The keypoints with each row as x, y, scale, orientation.
%   patch_radius:
%       The desired patch radius in pixels.
%
%   patches:
%       The local patches around each keypoint with normalized intensities.
%
% Copyright 2017: Johannes L. Schoenberger <jsch at inf.ethz.ch>

if size(keypoints, 1) == 0
    patches = zeros(0, 2 * patch_radius + 1, 2 * patch_radius + 1);
    return
end

[~, patches, ~] = vl_covdet(image, ...
                            'frames', vl_frame2oell(keypoints'), ...
                            'descriptor', 'patch', ...
                            'patchresolution', patch_radius);
patches = reshape(patches, [2 * patch_radius + 1, 2 * patch_radius + 1, ...
                            size(keypoints, 1)]);
patches = permute(patches, [3, 1, 2]);

end

author: ahojnnes

repo: local-feature-evaluation


text:
function disp(tt,name,varargin)
%Command window display of a QTT_TUCKER
%
disp_mode_sizes=false;
for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'modes'
            disp_mode_sizes=varargin{i+1};
        otherwise
            error('Unrecognized option: %s\n',varargin{i});
    end
end

if ~exist('name','var')
    name = 'ans';
end
fprintf('%s is a QTT-tucker tensor: \n',name);
fprintf('Tucker ranks: \n');
sz=size(tt.core);
rk=rank(tt.core);
for i=1:tt.dphys-1
   fprintf('%d-',sz(i));
end
fprintf('%d\n',sz(tt.dphys));
fprintf('TT-ranks for the core: \n');
for i=1:tt.dphys
   fprintf('%d-',rk(i));
end
fprintf('%d\n',rk(tt.dphys+1));
tk=tt.tuck;
fprintf('TT-ranks/mode sizes for the factors: \n');
for i=1:tt.dphys
   di=ndims(tk{i});
   ri=rank(tk{i});
   ni=size(tk{i});
   fprintf('Factor %d: ',i);
   for j=1:di
      fprintf('%d-',ri(j));
   end
   fprintf('%d\n',ri(di+1));
   if ( disp_mode_sizes )
     fprintf('Modes   : ');
     for j=1:di-1
       fprintf('%d-',ni(j));
     end
     fprintf('%d\n',ni(di));
   end
end
% %fprintf('r(1)=%d \n', r(1));
% for i=1:d
%    fprintf('r(%d)=%d \t n(%d)=%d \n',i,r(i),i,n(i));
% end
% fprintf('r(%d)=%d \n',d+1,r(d+1));

author: oseledets

repo: TT-Toolbox


text:
%% Settings

robotName='iCubGenova04'; %% Name of the robot

testFrames = {'l_hand', 'r_hand', 'head'}; %% The frames to display

meshFilePrefix = [getenv('ROBOTOLOGY_SUPERBUILD_INSTALL_PREFIX') '/share']; %% Path to the model meshes

modelPath = [getenv('ROBOTOLOGY_SUPERBUILD_INSTALL_PREFIX') '/share/iCub/robots/' robotName '/'];  %% Path to the robot model

fileName='model.urdf'; %% Name of the urdf file

jointOrder={     %% The list and the order of the joints
    'r_hip_pitch';
    'r_hip_roll';
    'r_hip_yaw';
    'r_knee';
    'r_ankle_pitch';
    'r_ankle_roll';
    'l_hip_pitch';
    'l_hip_roll';
    'l_hip_yaw';
    'l_knee';
    'l_ankle_pitch';
    'l_ankle_roll';
    'torso_pitch';
    'torso_roll';
    'torso_yaw';
    'r_shoulder_pitch';
    'r_shoulder_roll';
    'r_shoulder_yaw';
    'r_elbow';
    'r_wrist_prosup';
    'r_wrist_pitch';
    'r_wrist_yaw';
    'l_shoulder_pitch';
    'l_shoulder_roll';
    'l_shoulder_yaw';
    'l_elbow';
    'l_wrist_prosup';
    'l_wrist_pitch';
    'l_wrist_yaw';
    'neck_pitch';
    'neck_roll';
    'neck_yaw'
    };

groundFrame = 'l_sole'; %% The frame where to place the ground

rootLink = 'root_link'; %% The base link

joints_positions=zeros(length(jointOrder),1); %% The joints'position

printModel = true; %% Define if the model has to be printed on the command window

%% No need to edit from here

% Main variable of iDyntreeWrappers used for many things including updating
% robot position and getting world to frame transforms
KinDynModel = iDynTreeWrappers.loadReducedModel(jointOrder, rootLink, modelPath,fileName,false);

if printModel
    KinDynModel.kinDynComp.model().toString() %%Print model
end

% Set initial position of the robot
iDynTreeWrappers.setRobotState(KinDynModel,eye(4),joints_positions,zeros(6,1),zeros(size(joints_positions)),[0,0,-9.81]);

groundFrameTransform = iDynTreeWrappers.getWorldTransform(KinDynModel, groundFrame);

world_H_base = groundFrameTransform^-1;

iDynTreeWrappers.setRobotState(KinDynModel,world_H_base,joints_positions,zeros(6,1),zeros(size(joints_positions)),[0,0,-9.81]);

% Prepare figure, handles and variables required for the update, some extra
% options are commented.
[visualizer,objects]=iDynTreeWrappers.prepareVisualization(KinDynModel,meshFilePrefix,...
    'color',[1,1,1],'transparency',1, 'name', ['Plot frame ', robotName], 'reuseFigure', 'name');
xlim([-0.5, 0.5])
ylim([-0.5, 0.5])
zlim([-0.05, 1.4])

for t = 1 : length(testFrames)
    frameTransform = iDynTreeWrappers.getWorldTransform(KinDynModel, testFrames{t});
    testFramesObjects{t} = iDynTreeWrappers.plotFrame(frameTransform, 0.2, 5);
end

author: robotology

repo: whole-body-controllers


text:
% edu.stanford.covert.io.editMatFileHeader and
% edu.stanford.covert.io.readMatFileHeader test cases
%
% Author: Jonathan Karr
% Author: Jared Jacobs
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 9/14/2010
classdef editMatFileHeaderTest < TestCase
    methods
        function this = editMatFileHeaderTest(name)
            this = this@TestCase(name);
        end

        function testEditMatFileTimestamp(~)
            %generate mat file
            matFile = [tempname '.mat'];
            save(matFile, 'matFile', 'matFile');
            timestamp = edu.stanford.covert.io.readMatFileHeader(matFile);
            assertFalse(733899 == timestamp);

            %change timestamp
            edu.stanford.covert.io.editMatFileHeader(matFile, 733899, 'PCWIN', '5.0');
            [timestamp, platform, version] = ...
                edu.stanford.covert.io.readMatFileHeader(matFile);
            assertEqual(733899, timestamp);
            assertEqual('PCWIN', platform);
            assertEqual('5.0', version);

            %cleanup
            delete(matFile);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell