text:
% Construction Data_MI
% MI: 50 trials
% 4-class (Cylindrical/Spherical/LeftTwist/RightTwist)
%% Initializing
clc; close all; clear all;
%% file
dd='C:\Users\Doyeunlee\Desktop\Journal_dylee\1_ConvertedData\';
fsldx=2;
fs={'100','250','1000'};
ref='FCz';
channelMatrix={'F3','F1','Fz','F2','F4';
    'FC3','FC1','FCz','FC2','FC4';
    'C3','C1', 'Cz', 'C2', 'C4';
    'CP3','CP1','CPz','CP2','CP4';
    'P3','P1','Pz','P2','P4'};

trainRatio=0.8;
modldx=1;
modality={'MI','realMove'};

filelist={'sub01','sub02','sub03','sub04','sub05'};

%%
for sub = 1:length(filelist)
    ival = [0 3000];

    % Data load - MI
    [cntGrasp,mrkGrasp,mntGrasp]=eegfile_loadMatlab([dd 'realMove' '\' fs{fsldx} '\' filelist{sub} '_multigrasp_' 'realMove']);
    [cntTwist,mrkTwist,mntTwist]=eegfile_loadMatlab([dd 'realMove' '\' fs{fsldx} '\' filelist{sub} '_twist_' 'realMove']);

    % theta, mu, and beta band
    filterBank = {'[4 8]','[8 13]','[13 40]'};

    for filt = 1:length(filterBank)
        filtBank = {filt};

        % Butterworth
        cntGrasp=proc_filtButter(cntGrasp,5,[4 40]);
        cntTwist=proc_filtButter(cntTwist,5,[4 40]);

        % CAR
        cntGrasp = proc_commonAverageReference(cntGrasp);
        cntTwist = proc_commonAverageReference(cntTwist);

        epoGrasp=cntToEpo(cntGrasp,mrkGrasp,ival);
        epoTwist=cntToEpo(cntTwist,mrkTwist,ival);

        epoGrasp=proc_selectChannels(epoGrasp,{'F3','F1','Fz','F2','F4',...
            'FC3','FC1','FCz','FC2','FC4',...
            'C3','C1', 'Cz', 'C2', 'C4', ...
            'CP3','CP1','CPz','CP2','CP4',...
            'P3','P1','Pz','P2','P4'});
        epoTwist=proc_selectChannels(epoTwist,{'F3','F1','Fz','F2','F4',...
            'FC3','FC1','FCz','FC2','FC4',...
            'C3','C1', 'Cz', 'C2', 'C4', ...
            'CP3','CP1','CPz','CP2','CP4',...
            'P3','P1','Pz','P2','P4'});

        % Retrieve trial amounts
        %% 4-class: Forward, Cylindrical, LeftTwist, Rest
        epoCylindrical=proc_selectClasses(epoGrasp,{'Cylindrical'});
        epoSpherical=proc_selectClasses(epoGrasp,{'Spherical'});
        epoLeftTwist=proc_selectClasses(epoTwist,{'LeftTwist'});
        epoRightTwist=proc_selectClasses(epoTwist,{'RightTwist'});

        % append
        epo = proc_appendEpochs(epoCylindrical,epoSpherical);
        epo = proc_appendEpochs(epo,epoLeftTwist);
        epo = proc_appendEpochs(epo,epoRightTwist);

        % Data shuffling
        epoCylindrical.x=datasample(epoCylindrical.x,size(epoCylindrical.x,3),3,'Replace',false);
        epoSpherical.x=datasample(epoSpherical.x,size(epoSpherical.x,3),3,'Replace',false);
        epoLeftTwist.x=datasample(epoLeftTwist.x,size(epoLeftTwist.x,3),3,'Replace',false);
        epoRightTwist.x=datasample(epoRightTwist.x,size(epoRightTwist.x,3),3,'Replace',false);

        % CSP
        [csp_fv,csp_w,csp_eig]=proc_multicsp(epo,3);
        proc=struct('memo','csp_w');

        proc.train=['[fv,csp_w]= proc_multicsp(fv, 3);' ...
            'fv=proc_variance(fv);' ...
            'fv=proc_logarithm(fv);'];

        proc.apply=['fv=proc_linearDerivation(fv, csp_w);','fv=proc_variance(fv);','fv=proc_logarithm(fv);'];

        % Setting training Data
        minTrial=min([size(epoCylindrical.x,3), ...
            size(epoSpherical.x,3), ...
            size(epoLeftTwist.x,3), size(epoRightTwist.x,3)]);
        trainNum=minTrial*trainRatio;

        % training
        epoCylindricalTrain=proc_selectEpochs(epoCylindrical,1:trainNum);
        epoSphericalTrain=proc_selectEpochs(epoSpherical,1:trainNum);
        epoLeftTwistTrain=proc_selectEpochs(epoLeftTwist,1:trainNum);
        epoRightTwistTrain=proc_selectEpochs(epoRightTwist,1:trainNum);

        % test
        epoCylindricalTest=proc_selectEpochs(epoCylindrical,trainNum+1:size(epoCylindrical.x,3));
        epoSphericalTest=proc_selectEpochs(epoSpherical,trainNum+1:size(epoSpherical.x,3));
        epoLeftTwistTest=proc_selectEpochs(epoLeftTwist,trainNum+1:size(epoLeftTwist.x,3));
        epoRightTwistTest=proc_selectEpochs(epoRightTwist,trainNum+1:size(epoRightTwist.x,3));

        % Concat
        epoTrain=proc_appendEpochs(epoCylindricalTrain,epoSphericalTrain);
        epoTrain=proc_appendEpochs(epoTrain,epoLeftTwistTrain);
        epoTrain=proc_appendEpochs(epoTrain,epoRightTwistTrain);
        Train=epoTrain;

        epoTest=proc_appendEpochs(epoCylindricalTest,epoSphericalTest);
        epoTest=proc_appendEpochs(epoTest,epoLeftTwistTest);
        epoTest=proc_appendEpochs(epoTest,epoRightTwistTest);
        Test=epoTest;

        %%
        trainX=epoTrain.x;
        testX=epoTest.x;
        trainY=epoTrain.y;testY=epoTest.y;
        trainY=vec2ind(trainY)-1;
        testY=vec2ind(testY)-1;
        str= char('_train');
        str2=char('_test');
        str3=char('.mat');
%         str4=char();
        save(['C:\Users\Doyeunlee\Desktop\AAAI_20190717\Dataconstruction_filter\realMove_train/',num2str(sub),filterBank{filt},str3],'trainX','trainY');
        save(['C:\Users\Doyeunlee\Desktop\AAAI_20190717\Dataconstruction_filter\realMove_test/',num2str(sub),filterBank{filt},str3],'testX','testY');
    end % filterBank
end % subject

disp('Done');

author: PatternRecognition

repo: OpenBMI


text:
function []=Main(result,speed,options,problem,sgd)
global mones  mzeros convert  in_size gate_size numMmcell globalData globalMask usegpu...
    out_size share_size share_size2  in ingate cellstate cells outgate...
    node_outgateInit cellinInit node_cellbiasInit delta_outInit cellstatusInit signum

 signum = 3.5;
if problem.batchsize>problem.numsamples
    problem.batchsize=problem.numsamples;
end

[mones,mzeros,convert,usegpu]=gputype(speed.usegpu);
[data,mask,test,masktest ]=feval(['gen' problem.name],problem);
problem.numsamples=size(data,1);

[in_size,gate_size,out_size,share_size,share_size2,numMmcell,W,in,ingate,cellstate,cells,outgate,...
    node_outgateInit,cellinInit,node_cellbiasInit,delta_outInit,cellstatusInit ]=netInit(problem);
if problem.continualPredict
    Bpfunc=@batch_equal_nomask_lstm;
else
    Bpfunc=@batch_cell_lstm;
end

if speed.gradientchecking

    for i=1:fix(problem.numsamples/problem.batchsize)

        bashdata{i} =  data(  (i-1)*problem.batchsize + 1:i*problem.batchsize , : , :   );
        maskdata{i}= mask(  (i-1)*problem.batchsize + 1:i*problem.batchsize , : , :   );
    end

    globalData=bashdata{1};
    globalMask=maskdata{1};
    [ err,dw,inLL,right ]=Bpfunc(W);

    [numgrad]=computeNumericalGradient(Bpfunc,W);
    disp( [numgrad dw])
    diff = norm(numgrad-dw)/norm(numgrad+dw);
    disp(diff);
    fprintf('Norm of the difference between numerical and analytical gradient (should be < 1e-9)\n\n');
    return
end

if speed.usecluster
    serveroptions.slavecount=speed.numcluster;
    serveroptions.slavedir='/usr/local/anew-lstm';
    server=Server(serveroptions);

    markers = floor(linspace(1, size(data,2), server.slaveCount+1));
    markers(end) = size(data,2);
    markers(1) = 0;
    trainStarts = cell(1,server.slaveCount);
    trainEnds =cell(1,server.slaveCount);
    for i=1:server.slaveCount
        trainStarts{i} = markers(i)+1;
        trainEnds{i} = markers(i+1);
        fprintf('i = %d, trainStart = %d, trainEnd = %d\n', i, trainStarts{i}, trainEnds{i});
    end

    options.Method = 'lbfgs';

    for i=1:fix(problem.numsamples/problem.batchsize/server.slaveCount )
        bashdata{i} =  data(  (i-1)*problem.batchsize*server.slaveCount  + 1:...
            i*problem.batchsize*server.slaveCount  , : , :   );
        maskdata{i}= mask(  (i-1)*problem.batchsize*server.slaveCount + 1:...
            i*problem.batchsize*server.slaveCount  , : , :   );
    end
    clear data

    record=1;
    timer=tic;
    for epoch =1:1
        for i=1:fix(problem.numsamples/problem.batchsize/server.slaveCount )
            [Winit]=  server.rpc('clientLoadDataMinibatchNomask_ref',...
                trainStarts, trainEnds, problem.batchsize  ,...
                problem,speed,bashdata{i},maskdata{i},{1:server.slaveCount});
            if i==1&&epoch==1
                W=Winit{1};
            end
            for repeat = 1:32
                [W, cost] = minFunc( @server_batch_cell_lstm,W,options,server);
                if  toc(timer)>1
                    timearray{1}( record )= toc(timer)
                    [errorarray{1}(record),~,~,rightarray{1}(record) ]=testmodel(W,test,masktest ,problem);
                    timer=tic;
                    record=record+1;
                end
            end
        end
    end
    figure
    for method = 1:1
        plot(cumsum(timearray{method}), rightarray{method} , result.linspec{method}); hold on
    end
    legend('server  ' );
    xlabel('time(seconds)');
    ylabel('test obj');

    saveas(gcf,[ pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest) problem.name '.fig']);
    saveas(gcf,[ pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest)  problem.name '.eps'],'epsc');
    eps2pdf([pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest)  problem.name '.eps']);

    save('server.mat','timearray','errorarray');

elseif ~speed.usecluster

    for i=1:fix(problem.numsamples/problem.batchsize)

        bashdata{i} =  data(  (i-1)*problem.batchsize + 1:i*problem.batchsize , : , :   );
        maskdata{i}= mask(  (i-1)*problem.batchsize + 1:i*problem.batchsize , : , :   );
    end

    Winit=W;

   timer=tic;
   record=1;
    for i=1:fix(problem.numsamples/problem.batchsize)
        for inrepeat=1:20

            globalData =  bashdata{i};
            globalMask= maskdata{i};

            options.Method = 'lbfgs';

            [W, cost] = minFunc( Bpfunc,W,options);
            if toc(timer)>2

                timearray{1}( record  )= toc(timer);
                [errorarray{1}(record),~,~,rightarray{1}(record) ]=testmodel(W,test,masktest ,problem);
                disp([ 'lbfgs error '  num2str(errorarray{1}(record) )]  )
                record=record+1;
                timer=tic;
            end
        end
    end

        W=Winit;
        momentum=sgd.momentum;
        alpha=sgd.alpha;
        oldGradient=0;

        record =1;
        timer=tic;
        for repeat =1 :5
            for inner =1:50

                for  i=1:fix(problem.numsamples/problem.batchsize)

                    globalData =  bashdata{i};
                    globalMask= maskdata{i};

                    [~,dw,~,~]=Bpfunc(W);
                    oldGradient = alpha*dw + momentum* oldGradient;
                    W= W - alpha* oldGradient;
                end
                timearray{2}(record  )= toc(timer);
                [errorarray{2}(record),~,~, rightarray{2}(record) ]=testmodel(W,test,masktest ,problem);
                record=record+1;
                timer=tic;
            end
            disp([ 'sgd error '  num2str(errorarray{2}(record-1) )]  )
        end

    figure
    for method = 1:2
        plot(cumsum(timearray{method}), rightarray{method} , result.linespec{method}); hold on
    end
    legend('lbfgs','sgd');
    xlabel('time(seconds)');
    ylabel('test obj');

    figure
    for method = 1:2
        plot(cumsum(timearray{method}),errorarray{method} , result.linespec{method}); hold on
    end
    legend('lbfgs','sgd');
    xlabel('time(seconds)');
    ylabel('test e');

    if result.savepic
        saveas(gcf,[ pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest) problem.name '.fig']);
        saveas(gcf,[ pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest)  problem.name '.eps'],'epsc');
     %   eps2pdf([pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest)  problem.name '.eps']);
    end

end

end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function colorDetectHSV(fileName, hsvVal, tol)

%
% function colorDetectHSV(fileName, hsvVal, tol)
%
% This function is used for detecting a specified hsv value in images.
%
% ARGUMENTS:
% fileName: the name of the jpg file to be loaded
% hsvVal: 3x1 array containing the HSV values to be detected
% tol: 1x1 or 2x1 or 3x1 array containing the tolerance (i.e. the maximum
% distance - in each hsv coefficient - of each pixel from hsvVal).
%
% Example:
% colorDetectHSV('train/face07.jpg', median(HSV), [0.05 0.05 0.1]);
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Theodoros Giannakopoulos - January 2008
% www.di.uoa.gr/~tyiannak
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

RGB = imread(fileName);

HSV = rgb2hsv(RGB);

% find the difference between required and real H value:
diffH = abs(HSV(:,:,1) - hsvVal(1));

[M,N,t] = size(RGB);
I1 = zeros(M,N); I2 = zeros(M,N); I3 = zeros(M,N);

T1 = tol(1);

I1( find(diffH < T1) ) = 1;

if (length(tol)>1)
    % find the difference between required and real S value:
    diffS = abs(HSV(:,:,2) - hsvVal(2));
    T2 = tol(2);
    I2( find(diffS < T2) ) = 1;
    if (length(tol)>2)
        % find the difference between required and real V value:
        difV = HSV(:,:,3) - hsvVal(3);
        T3 = tol(3);
        I3( find(diffS < T3) ) = 1;
        I = I1.*I2.*I3;
    else
        I = I1.*I2;
    end
else
    I = I1;
end

subplot(1,2,1),imshow(RGB); title('Original Image');
subplot(1,2,2),imshow(I,[]); title('Detected Areas');

author: Sable

repo: mcbench-benchmarks


text:
function x = solvedbi_sm(ah, rho, b, c)

% solvedbi_sm -- Solve a diagonal block linear system with a scaled identity
%                term using the Sherman-Morrison equation
%
%         The solution is obtained by independently solving a set of linear
%         systems of the form (see wohlberg-2016-efficient)
%
%                  (a a^H + rho I) x = b
%
%         In this equation inner products and matrix products are taken along
%         the 3rd dimension of the corresponding multi-dimensional arrays; the
%         solutions are independent over the 1st and 2nd (and 4th, if
%         non-singleton) dimensions.
%
% Usage:
%       x = solvedbi_sm(ah, rho, b, c);
%
% Input:
%       ah          Multi-dimensional array containing a^H
%       rho         Scalar rho
%       b           Multi-dimensional array containing b
%       c           Multi-dimensional array containing pre-computed quantities
%                   a^H / (a^H a + rho)
%
% Output:
%       x           Multi-dimensional array containing linear system solution
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2014-12-18
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

a = conj(ah);

if nargin < 4 || isempty(c),
  c = bsxfun(@rdivide, ah, sum(ah.*a, 3) + rho);
end

cb = sum(bsxfun(@times, c, b), 3);
clear ah c;
cba = bsxfun(@times, cb, a);
clear a cb;
x = (b - cba) / rho;

return

author: thfylsty


text:
%% flandmark_learn.m
% TODO - add description
%
% 14-07-10 Michal Uricar
% 17-05-11 Michal Uricar
% 02-08-11 Michal Uricar, flandmark2 - detection of corners of eyes and mouth
% 21-03-12 Michal Uricar, LFW annotation in one file
% 13-09-12 Michal Uricar, refactoring

clc;
close all; clearvars;

%% User defined options

% setup paths to face database (and benchmarks)
ROOT_FACE_FOLDER = '../LFW/';                           % path to folder with database images
DB_ANN = '../data/LFW_annotation.txt';                  % annotation file for the database
%TRN_EXAMPLES = '../benchmarks/LFW_TRN_annotation.txt';  % benchmark db division, leave empty to create different division
%VAL_EXAMPLES = '../benchmarks/LFW_VAL_annotation.txt';  % benchmark db division, leave empty to create different division
%TST_EXAMPLES = '../benchmarks/LFW_TST_annotation.txt';  % benchmark db division, leave empty to create different division
TRN_EXAMPLES = [];
VAL_EXAMPLES = [];
TST_EXAMPLES = [];

% options configuration
% order of components in LFW database
options.compnames = { 'center', 'canthus-rl', 'canthus-lr', 'mouth-corner-r', 'mouth-corner-l', ...
                      'canthus-rr', 'canthus-ll', 'nose' };
options.bw = [40; 40];
options.components = [20 10 10 10 10 10 10 10;
                      20 10 10 10 10 10 10 10];
options.comselect = [10, 4, 5, 8, 9, 3, 6, 10];
% number of components
options.M = size(options.components, 2);
% base window margin (percents)
options.bw_margin = [20; 20];
% path to LFW database (where images and anotations are stored) must contain subdir /mgt_v2 with xml annotations
options.image_path = ROOT_FACE_FOLDER;

%% Timestamp

fprintf(1,'Started on %s\n\n', datestr(now));

%% Add paths

	addpath('./code/');

%% Check project dir structure

    fprintf('Checking project dir structure...\n');

    DIRS = {'./code/img/', './code/results/', './code/MAT'};

    for i = 1 : numel(DIRS)
        fprintf('%s ... ', DIRS{i});
        if (~exist(DIRS{i}, 'dir'))
            fprintf('MISSING\n');
            fprintf('Creating dir %s...\n', DIRS{i});
            mkdir(DIRS{i});
        else
            fprintf('OK\n');
        end;
    end;

    DIRS = {'./bmrm', './libqp', './code/Functions', ...
            ROOT_FACE_FOLDER};

    for i = 1 : numel(DIRS)
        fprintf('%s ... ', DIRS{i});
        if (~exist(DIRS{i}, 'dir'))
            fprintf('MISSING\n');
            error('??? Corrupted project structure.');
        else
            fprintf('OK\n');
        end;
    end;

%% Save above specified user data

    options.image_path = ['../' ROOT_FACE_FOLDER];

    save('./code/MAT/paths.mat', 'ROOT_FACE_FOLDER', 'TRN_EXAMPLES', 'VAL_EXAMPLES', 'TST_EXAMPLES', 'DB_ANN');
    save('./code/MAT/options.mat', 'options');

%% Compile mex files if neccessary

    fprintf('Checking for neccessary mex-files...\n');

    if (~exist(['../matlab_toolbox/mex/lbp2sparse.' mexext], 'file') || ...
        ~exist(['../matlab_toolbox/mex/argmax_mex.' mexext], 'file') || ...
        ~exist(['../matlab_toolbox/mex/getNormalizedFrame.' mexext], 'file') || ...
        ~exist(['../matlab_toolbox/mex/flandmark_detector.' mexext], 'file'))
        error(['??? Mex-files not found! Please compile the flandmark with BUILD_MATLAB_BINDINGS turned on'...
               ' or run /matlab_toolbox/flandmark_compilemex.m to create them.\n']);
    end

    fprintf('Done.\n');

%% Collect statistics about components positions

    % check if statistics already exist
    if (exist('./code/results/statistics_orig.mat', 'file'))
        fprintf('Statistics found. Skipping...\n');
    else
        fprintf('Collecting statistics about components...\n');

        run './code/runstats.m';
    end;

%% Division of database into 3 parts (60% traning, 20% validating, 20% testing)

    % check if division already done.
    if (exist('./code/MAT/TRN.mat', 'file') && ...
            exist('./code/MAT/VAL.mat', 'file') && ...
            exist('./code/MAT/TST.mat', 'file'))
        fprintf('Division already done. Skipping...\n');
    else
        fprintf('Dividing database into 3 parts...\n');

        run './code/dbDivision.m';
    end;

%% Create data for BMRM solver (TRN, VAL, TST)

    %% Functions selection

    fprintf('Prepare functions configurations\n');
    f_argmax = @argmax2;
    f_risk = @risk2;
%     f_risk = @risk2_unnormalized;
    f_buildPsi = @buildPsi2;
    f_createMapTable = @createMapTable2;
    displacement = true;
    save('./code/MAT/functions_gdisp.mat', 'f_argmax', 'f_risk', ...
        'f_buildPsi', 'f_createMapTable', 'displacement');

    %% Configuration - nImages, height_of_pyramid
    fprintf('Looking for config file...\n');
    if (~exist('./code/MAT/config.mat', 'file'))
        fprintf('Not found. Creating...\nconfig = \n');
        config.nImages = 'all';
        config.heightOfPyramid = 4;
        disp(config);
        save('./code/MAT/config.mat', 'config');
    else
        fprintf('File %s found.\n', './code/MAT/config.mat');
    end;

    fprintf('Prepare data for LBP\n');
    run './code/prepDataForLBP.m';

%% Call BMRM solver

    % prepare lamdas
    fprintf('Looking for file with lambda values...\n');
    if (~exist('./code/MAT/lambdas.mat', 'file'))
        lambdas = [1e-3, 1e-2, 1e-1, 1];
        fprintf('Not found. Creating...\n');
        save('./code/MAT/lambdas.mat', 'lambdas');
        disp(lambdas);
    else
        fprintf('Lambda values found:');
        load('./code/MAT/lambdas.mat');
        disp(lambdas);
    end;

    % prepare normalization coefficients for each image in database
    if (~exist('./code/MAT/kappas.mat', 'file'))
	run('./code/computeKappas.m');
    end;

    fprintf('Calling BMRM solver...\n');
    run './code/test_bmrm.m';

%% Validation

    if (exist('./code/results/W.mat', 'file'))
        fprintf('Best vector W already found. Skipping...\n');
    else
        fprintf('Running validaton...\n');

        run './code/testValidation.m';
    end;

%% Test

    if (exist('./code/results/errors.mat', 'file'))
        fprintf('Test already passed. Skipping...\n');
    else
        fprintf('Running test...\n');

        run './code/test1.m';
    end;

    fprintf('Show histograms.\n');
    run './code/test2.m';

%% Generate model.mat (for use with flandmark_load_model - binary dat file creation)

    run './code/detectorModelBuild.m';

%% Timestamp

 fprintf(1,'Finished on %s.\n', datestr(now));

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function vl_demo_kdtree
% VL_DEMO_KDTREE
%   Demonstrates the construction of a kd-tree for uniformly
%   distributed 2-D data and data distributed on a circle. Both the
%   median and mean threshold selection methods are illustrated.

  randn('state',0) ;
  rand('state',0) ;

  % uniform 2D points
  Q = single(rand(2,1)) ;
  X = single(rand(2, 100)) ;

  % kd-tree with median and mean thresholds
  figure(1) ; clf ; do('thresholdmethod', 'mean') ;
  figure(2) ; clf ; do('thresholdmethod', 'median') ;

  % 2D points on a circle
  thr = 2*pi*rand(1,100) ;
  X = single(.4*[cos(thr); sin(thr)] + .5) ;
  Q = single(Q) ;

  % kd-tree with median and mean thresholds
  figure(3) ; clf ; do('thresholdmethod', 'mean') ;
  figure(4) ; clf ; do('thresholdmethod', 'median') ;

  vl_demo_print(1, 'kdtree_uniform_mean') ;
  vl_demo_print(2, 'kdtree_uniform_median') ;
  vl_demo_print(3, 'kdtree_circle_mean') ;
  vl_demo_print(4, 'kdtree_circle_median') ;

  function do(varargin)
    % build the forest
    kdforest = vl_kdtreebuild (X, 'verbose', varargin{:}) ;

    % query the kd-tree
    [i, d] = vl_kdtreequery (kdforest, X, Q, 'numneighbors', 10, 'verbose') ;
    vl_plotframe(X, 'ro') ;
    hold on ;
    xl = [0, +1] ;
    yl = [0, +1] ;
    axis equal ;
    xlim(xl) ;
    ylim(yl) ;
    vl_demo_kdtree_plot(kdforest.trees(1), 1, xl, yl) ;

    % plot query point and returns
    vl_plotframe(Q,'b*') ;
    for k=1:length(i)
      vl_plotframe([Q ; sqrt(d(k))],'b:','linewidth',1) ;
      vl_plotframe(X(:, i(k)), 'bx') ;
    end
    set(gca,'xtick',[],'ytick',[]) ;
  end

end

author: Cloud-CV

repo: object-proposals


text:
function tapas_sem_prosa_llh_check_input(t, a, u, theta)
%% Checks the input to the likelihood function
%
% aponteeduardo@gmail.com
% copyright (C) 2015
%

DIM_THETA = tapas_sem_prosa_ndims();

ns = size(t, 1);
np = size(theta, 2);

check_matrix(t, [ns, 1], 't');
check_matrix(a, [ns, 1], 'a');
check_matrix(u, [ns, 1], 'u');

if np > 1
    assert(np == ns, sprintf('tapas:sem:llh:input:size'), ...
        'Number fo parameters should be %d, instead %d', ns, np);
end

check_matrix(theta, [2 * DIM_THETA, np], 'theta');
%check_delays(t, a, u, theta);

end

function check_matrix(mat, ns, field)
%% Checks a matrix in theta.
%
% Input
%   mat -- Matrix
%   ns -- Expected size of the matrix
%   field -- String with the name of the parameter

assert(isnumeric(mat), ...
    sprintf('tapas:sem:llh:input:%s:not_numeric', field), ...
    '%s should be numeric', field);
assert(isreal(mat), ...
    sprintf('tapas:sem:llh:input:%s:not_real', field), ...
    '%s should be real', field);
assert(~issparse(mat), ...
    sprintf('tapas:sem:llh:input:%s:sparse', field), ...
    '%s should not be sparse', field);
assert(ndims(mat) == 2, ...
    sprintf('tapas:sem:llh:input:%s:ndim', field), ...
    '%s should have 2 dimensions', field)
assert(all(size(mat) == ns), ...
    sprintf('tapas:sem:llh:input:%s:dsize', field), ...
    '%s should have dimensions [%d,%d], instead [%d,%d].', ...
    field, ns(1), ns(2), size(mat, 1), size(mat, 2));
%assert(all(mat(:) >= 0 ), ...
%    sprintf('tapas:sem:llh:input:%s:range', field), ...
%    '%s should be positive', field)
end

function check_delays(t, a, u, theta)
% Check if the delays make a trial implaussible

[PROSACCADE, ANTISACCADE] = tapas_sem_action_codes();

assert(all(t > theta(13)), 'tapas:sem:llh:input:main_delay', ...
    'Main delay %0.5f is higher than at least one saccade', theta(13));

assert(all(t(a == ANTISACCADE & u == PROSACCADE) > theta(13) + theta(14)),...
    'tapas:sem:llh:input:pro_delay', ...
    ['Delay of antisaccade in a prosaccade trial %0.5f is higher than at '...
    'least one saccade'], theta(13) + theta(14));

assert(all(t(a == ANTISACCADE & u == PROSACCADE) > theta(13) + theta(15)),...
    'tapas:sem:llh:input:anti_delay', ...
    ['Delay of antisaccade in an antisaccade trial %0.5f is higher than at ' ...
   'least one saccade'], theta(13) + theta(15));

end

author: translationalneuromodeling

repo: tapas


text:
function [Signal] = CliffFrontRightSignalStrengthRoomba(serPort);
%[Signal] = CliffFrontRightSignalStrengthRoomba(serPort)
%Displays the strength of the front right cliff sensor's signal.
%Ranges between 0-100 percent signal

% By; Joel Esposito, US Naval Academy, 2011

%Initialize preliminary return values
Signal = nan;

try

%Flush Buffer
N = serPort.BytesAvailable();
while(N~=0)
fread(serPort,N);
N = serPort.BytesAvailable();
end

warning off
global td
fwrite(serPort, [142]);  fwrite(serPort,30);

Strength =  fread(serPort, 1, 'uint16');
Signal=(Strength/4095)*100;

pause(td)
catch
    disp('WARNING:  function did not terminate correctly.  Output may be unreliable.')
end

author: petercorke

repo: robotics-toolbox-matlab


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

% DTFT properties

%linearity

syms w
n= 0:10;
x1=0.8.^n;
x2=0.7.^n;
a=3; b=4;
L=a*x1+b*x2;
Left=sum(L.*exp(-j*w*n));
ezplot(abs(Left),[-pi pi])
title('| DTFT(ax_1[n]+bx_2[n]) |')
figure
w1=-pi:.01:pi;
LL=subs(Left,w,w1)
plot(w1,angle(LL));
xlim([-pi pi])
title('Phase of DTFT')

figure
X1= sum(x1.*exp(-j*w*n));
X2= sum(x2.*exp(-j*w*n));
Right=a*X1+b*X2 ;
ezplot(abs(Right),[-pi pi])
title('Magnitude of right side')
figure
w1=-pi:.01:pi;
RR=subs(Right,w,w1)
plot(w1,angle(RR));
xlim([-pi pi])
title('Phase of right part')


text:
function status = test_interpolation_methods(modus)
%TEST_INTERPOLATION_METHODS demonstrates the difference between interpolation
%methods
%
%   Usage: status = test_interpolation_methods(modus)
%
%   Input parameters:
%       modus   - 0: numerical (not available)
%                 1: visual
%
%   Output parameters:
%       status  - true or false

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

status = false;

%% ===== Checking of input parameters ====================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);
if ~modus
    warning('%s: numerical modus not available.',upper(mfilename));
    status = true;
    return;
end

%% ===== Settings ========================================================
conf = struct;
conf.ir.useinterpolation = true;

%% ===== Main ============================================================
% Check if HRTF data set with low frequency correction is available,
% download otherwise
basepath = get_sfs_path();
hrtf = 'KEMAR_HRTFs_lfcorr.sofa';
hrtf_file = [basepath '/data/HRTFs/' hrtf];
if ~exist(hrtf_file,'file')
    disp('Download');
    url = ['https://github.com/spatialaudio/lf-corrected-kemar-hrtfs/', ...
           'blob/master/' hrtf '?raw=true'];
    download_file(url,hrtf_file);
end
% Load HRTF data set
hrtf = SOFAload(hrtf_file);

%% ===== Interpolate shifted Dirac impulses ==============================
% interpolation points
x0 = [1 3; 0 0; 0 0]/180*pi;
% weights (for target point: xs = [2; 0; 0]/180*pi )
weights = [.5; .5];
N = 50; %length impulse responses

% 1. Interpolate between Dirac impulses with one sample in between
h1 = zeros(1,N);
h1(3) = 1;
h2 = zeros(1,N);
h2(5) = 1;

irs1(1,1,:) = h1;
irs1(2,1,:) = h2;
irs1(:,2,:) = irs1(:,1,:); %redundant second channel

conf.ir.interpolationmethod = 'simple';
h_int1_simple = interpolate_ir(irs1,weights,conf);
conf.ir.interpolationmethod = 'freqdomain';
h_int1_fd = interpolate_ir(irs1,weights,conf);
conf.ir.interpolationmethod = 'timedomain';
h_int1_td = interpolate_ir(irs1,weights,conf);

% 2. Interpolate between neighbouring Dirac impulses at impulse reponse start
h1 = zeros(1,N);
h1(3) = 1;
h3 = zeros(1,N);
h3(4) = 1;

irs2(1,1,:) = h1;
irs2(2,1,:) = h3;
irs2(:,2,:) = irs2(:,1,:); %redundant second channel

conf.ir.interpolationmethod = 'simple';
h_int2_simple = interpolate_ir(irs2,weights,conf);
conf.ir.interpolationmethod = 'freqdomain';
h_int2_fd = interpolate_ir(irs2,weights,conf);
conf.ir.interpolationmethod = 'timedomain';
h_int2_td = interpolate_ir(irs2,weights,conf);

% 3. Interpolate between neighbouring Dirac impulses in middle of impulse response
h4 = zeros(1,N);
h4(25) = 1;
h5 = zeros(1,N);
h5(26) = 1;

irs3(1,1,:) = h4;
irs3(2,1,:) = h5;
irs3(:,2,:) = irs3(:,1,:);

conf.ir.interpolationmethod = 'simple';
h_int3_simple = interpolate_ir(irs3,weights,conf);
conf.ir.interpolationmethod = 'freqdomain';
h_int3_fd = interpolate_ir(irs3,weights,conf);
conf.ir.interpolationmethod = 'timedomain';
h_int3_td = interpolate_ir(irs3,weights,conf);

% Plots
% impulse responses
ir_list = {
    irs1(1,1,:)         , irs2(1,1,:)          , irs3(1,1,:)
    irs1(2,1,:)         , irs2(2,1,:)          , irs3(2,1,:)
    h_int1_simple(1,1,:), h_int2_simple(1,1,:) , h_int3_simple(1,1,:)
    h_int1_fd(1,1,:)    , h_int2_fd(1,1,:)     , h_int3_fd(1,1,:)
    h_int1_td(1,1,:)    , h_int2_td(1,1,:)     , h_int3_td(1,1,:)
    };  % list of IRs

col = {'k' 'b' 'r' 'm' 'c'};

titlelabels = {
    'Interpolation between Dirac impulses with one sample in between'
    'Interpolation of neighbouring Dirac impulses at impulse response start'
    'Interpolation of neighbouring Dirac impulses in middle of impulse response'
    };

t = 0:N-1;  % time-axis in samples

for cdx=1:size(ir_list,2)
    figure
    for rdx=1:size(ir_list,1)
        plot(t,squeeze(ir_list{rdx, cdx}), col{rdx})
        hold on
    end
    hold off
    grid
    xlabel('samples')
    ylabel('amplitude')
    legend('h_1','h_2','simple interp','freqdomain interp','timedomain interp')
    title(titlelabels{cdx})
end

%% ===== Interpolate HRIRs ===============================================
% 1. Interpolate between close HRIRs
idx0 = 181; %index for 0° azimuth
idx1 = 182; %index for 1° azimuth
idx2 = 183; %index for 2° azimuth
x0_close = [hrtf.SourcePosition(idx0,:).' hrtf.SourcePosition(idx2,:).']/180*pi;
% weights (for target point: xs_close = hrtf.SourcePosition(idx1,:).'/180*pi )
weights_close = [.5; .5];
hrir_close = [hrtf.Data.IR(idx0,:,:); hrtf.Data.IR(idx2,:,:)];
hrir_close_ref = hrtf.Data.IR(idx1,:,:);

conf.ir.interpolationmethod = 'simple';
hrir_close_simple = interpolate_ir(hrir_close,weights_close,conf);
conf.ir.interpolationmethod = 'freqdomain';
hrir_close_fd = interpolate_ir(hrir_close,weights_close,conf);
conf.ir.interpolationmethod = 'timedomain';
hrir_close_td = interpolate_ir(hrir_close,weights_close,conf);

% 2. Interpolate between distant HRIRs
idx0 = 181; %index for 0° azimuth
idx30 = 211; %index for 30° azimuth
idx60 = 241; %index for 60° azimuth
x0_dist = [hrtf.SourcePosition(idx0,:).' hrtf.SourcePosition(idx60,:).']/180*pi;
% weights (for target point: xs_dist = hrtf.SourcePosition(idx30,:).'/180*pi )
weights_dist = [.5; .5];
hrir_dist = [hrtf.Data.IR(idx0,:,:); hrtf.Data.IR(idx60,:,:)];
hrir_dist_ref = hrtf.Data.IR(idx30,:,:);

conf.ir.interpolationmethod = 'simple';
hrir_dist_simple = interpolate_ir(hrir_dist,weights_dist,conf);
conf.ir.interpolationmethod = 'freqdomain';
hrir_dist_fd = interpolate_ir(hrir_dist,weights_dist,conf);
conf.ir.interpolationmethod = 'timedomain';
hrir_dist_td = interpolate_ir(hrir_dist,weights_dist,conf);

% Plots
hrir_list = {
    hrir_close(1,1,:)       , hrir_dist(1,1,:)
    hrir_close(2,1,:)       , hrir_dist(2,1,:)
    hrir_close_ref(1,1,:)   , hrir_dist_ref(1,1,:)
    hrir_close_simple(1,1,:), hrir_dist_simple(1,1,:)
    hrir_close_fd(1,1,:)    , hrir_dist_fd(1,1,:)
    hrir_close_td(1,1,:)    , hrir_dist_td(1,1,:)
    };  % list of HRIRs

col = {'k' 'b' 'g' 'r' 'm' 'c'};
titlelabels = {'close', 'distant'};
legendpos = {'northeast', 'southwest'};

t = 0:hrtf.API.N-1;  % time-axis in samples
f = (0:hrtf.API.N-1)/hrtf.API.N*hrtf.Data.SamplingRate;  % frequency axis

for pdx=1:2  % impulse responses or magnitude spectrum
    for cdx=1:size(hrir_list,2)
        figure
        for rdx=1:size(hrir_list,1)
            if pdx == 1
                % impulse responses
                plot(t,squeeze(hrir_list{rdx, cdx}),col{rdx})
            else
                % magnitude spectra
                mag_spectrum = db(abs(fft(squeeze(hrir_list{rdx, cdx}))));
                semilogx(f(2:end), mag_spectrum(2:end), col{rdx});
            end
            hold on
        end
        hold off

        if pdx == 1
            % impulse responses
            xlabel('samples')
            ylabel('amplitude')
            axis([0 160 -0.6 0.6])
        else
            % magnitude spectra
            xlabel('frequency in Hz')
            ylabel('amplitude in dB')
            axis([f(2) hrtf.Data.SamplingRate/2 -60 20])
        end
        legend('hrir_1','hrir_2','hrir_{ref}','simple interp', ...
            'freqdomain interp', 'timedomain interp','Location', legendpos{pdx})
        title(['Interpolation of ' titlelabels{cdx} ' HRIRs'])
    end
end

status = true;

author: sfstoolbox

repo: sfs-matlab


text:
function F = rdivide(F, G)
%./   Pointwise right divide for CHEBFUN3V objects.
%
% See also CHEBFUN3V/LDIVIDE.

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

% Empty check:
if ( isempty(F) || isempty(G) )
   F = chebfun3v;
   return
end

% Componentwise divide.
if ( isa(F, 'chebfun3v') && isa(G, 'chebfun3v') )
    nF = F.nComponents;
    nG = F.nComponents;
    if ( nF ~= nG )
        error('CHEBFUN:CHEBFUN3V:rdivide:dim','Dimension mismatch.');
    end
    for jj = 1:nF
        F.components{jj} = rdivide(F.components{jj}, G.components{jj});
    end

elseif ( isa(F, 'chebfun3v') && ( isa(G, 'chebfun3') || isa(G, 'double') ) )
    nF = F.nComponents;
    for jj = 1:nF
        F.components{jj} = rdivide(F.components{jj}, G);
    end

elseif  ( isa(G, 'chebfun3v') && ( isa(F, 'chebfun3') || isa(F, 'double') ) )
    H = G;
    nG = G.nComponents;
    for jj = 1:nG
        H.components{jj} = rdivide(F, H.components{jj});
    end
    F = H;

else
   error('CHEBFUN:CHEBFUN3V:rdivide:inputs','Unrecognized input arguments.')
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [s, cfg] = ft_statfun_mean(cfg, dat, design)

% FT_STATFUN_MEAN demonstrates how to compute the mean over all conditions in the
% data. Since this does NOT depend on the experimental design, it cannot be used for
% testing for differences between conditions.
%
% This function serves as an example for a statfun. You can use such a function with
% the statistical framework in FieldTrip using FT_TIMELOCKSTATISTICS,
% FT_FREQSTATISTICS or FT_SOURCESTATISTICS to perform a statistical test, without
% having to worry about the representation of the data.
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS, and see FT_STATFUN_DIFF for a similar example

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

if numel(unique(design(cfg.ivar,:)))>1
  ft_warning('inappropriate design, it should only contain one condition');
end

% return the means as the statistic of interest
s.stat = mean(dat,2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function annolist_orig = unflatten_annolist(annolist_flat,annolist_orig)

n = 0;
for imgidx = 1:length(annolist_orig)
    rect_orig = annolist_orig(imgidx).annorect;
    for ridx = 1:length(rect_orig)
        if (isfield(rect_orig(ridx),'objpos') && ~isempty(rect_orig(ridx).objpos))
            n = n + 1;
            rect_orig(ridx) = annolist_flat(n).annorect;
        end
    end
    annolist_orig(imgidx).annorect = rect_orig;
end

end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function out = webcam
%WEBCAM M-Code for creating a video input object.
%
%   This is the machine generated representation of a video input object.
%   This M-file, WEBCAM.M, was generated from the OBJ2MFILE function.
%   A MAT-file is created if the object's UserData property is not
%   empty or if any of the callback properties are set to a cell array
%   or to a function handle. The MAT-file will have the same name as the
%   M-file but with a .MAT extension. To recreate this video input object,
%   type the name of the M-file, webcam, at the MATLAB command prompt.
%
%   The M-file, WEBCAM.M and its associated MAT-file, WEBCAM.MAT (if
%   it exists) must be on your MATLAB path.
%
%   Example:
%       vidobj = webcam;
%
%   See also VIDEOINPUT, IMAQDEVICE/PROPINFO, IMAQHELP, PATH.
%

% Copyright 2011 The MathWorks, Inc.

% Check if we can check out a license for the Image Acquisition Toolbox.
canCheckoutLicense = license('checkout', 'Image_Acquisition_Toolbox');

% Check if the Image Acquisition Toolbox is installed.
isToolboxInstalled = exist('videoinput', 'file');

if ~(canCheckoutLicense && isToolboxInstalled)
    % Toolbox could not be checked out or toolbox is not installed.
    errID = 'imaq:obj2mfile:invalidToolbox';
    error(errID, 'Image Acquisition Toolbox is not installed.');
end

% Device Properties.
adaptorName = 'winvideo';
deviceID = 1;
vidFormat = 'RGB24_320x240';
tag = '';

% Search for existing video input objects.
existingObjs1 = imaqfind('DeviceID', deviceID, 'VideoFormat', vidFormat, 'Tag', tag);

if isempty(existingObjs1)
    % If there are no existing video input objects, construct the object.
    vidObj1 = videoinput(adaptorName, deviceID, vidFormat);
else
    % There are existing video input objects in memory that have the same
    % DeviceID, VideoFormat, and Tag property values as the object we are
    % recreating. If any of those objects contains the same AdaptorName
    % value as the object being recreated, then we will reuse the object.
    % If more than one existing video input object contains that
    % AdaptorName value, then the first object found will be reused. If
    % there are no existing objects with the AdaptorName value, then the
    % video input object will be created.

    % Query through each existing object and check that their adaptor name
    % matches the adaptor name of the object being recreated.
    for i = 1:length(existingObjs1)
        % Get the object's device information.
        objhwinfo = imaqhwinfo(existingObjs1{i});
        % Compare the object's AdaptorName value with the AdaptorName value
        % being recreated.
        if strcmp(objhwinfo.AdaptorName, adaptorName)
            % The existing object has the same AdaptorName value as the
            % object being recreated. So reuse the object.
            vidObj1 = existingObjs1{i};
            % There is no need to check the rest of existing objects.
            % Break out of FOR loop.
            break;
        elseif(i == length(existingObjs1))
            % We have queried through all existing objects and no
            % AdaptorName values matches the AdaptorName value of the
            % object being recreated. So the object must be created.
            vidObj1 = videoinput(adaptorName, deviceID, vidFormat);
        end %if
    end %for
end %if

% Configure vidObj1 properties.
set(vidObj1, 'ErrorFcn', @imaqcallback);
set(vidObj1, 'FramesAcquiredFcnCount', 100);
set(vidObj1, 'FramesAcquiredFcn', @(h, e) flushdata(h));
set(vidObj1, 'TriggerRepeat', Inf);

out = vidObj1 ;


text:
function [A, iter, alpha] = nesterov_mnls(X, B, A_star, alpha, max_iter, func_type)
% Accelerated proximal gradient algorithm for nonnegative constraint
%   min_A 1/2 | X - A*B' |^2_F,
%   s.t. A >= 0.
%

    tol_1 = 1e-2;
    tol_2 = 1e-2;

    R = size(B,2);
    Z = B' * B;

    if strcmp(func_type, 'smooth') || strcmp(func_type, 'stochastic') ||  strcmp(func_type, 'strong_alpha_beta')
        s = svd(Z);
        L = max(s);
        mu = min(s);
    else
        L = norm(Z, 'fro');
    end

    if strcmp(func_type, 'smooth') || strcmp(func_type, 'stochastic')
        W = - X * B;
        q = mu/L;
    elseif strcmp(func_type, 'strong_alpha_beta')
        lambda = g_lambda(L, mu);
        Z = Z + lambda * eye(R);
        W = - X * B - lambda * A_star;
        L = L + lambda;
        q = (mu + lambda)/L;
    else
        W = - X * B;
    end

    Y = A_star;
    A_prev = Y;

    %alpha = 1;
    iter = 0;
    while iter < max_iter
        grad = W + Y * Z;

        % terminate the loop
        cond_1 = max(max( abs(grad .* Y)));
        cond_2 = min((min (grad)));
        if cond_1 < tol_1 && cond_2 > - tol_2 && iter > 0
            break
        else
            A = Y - 1/L * grad;
            A = max(A, 0);
            if strcmp(func_type, 'smooth') || strcmp(func_type, 'strong_alpha_beta')
                new_alpha = update_alpha(alpha, q);
                beta = alpha * (1 - alpha) / (alpha^2 + new_alpha);
                alpha = new_alpha;
            elseif strcmp(func_type, 'stochastic')
                beta = (1-sqrt(mu * 1/L))/(1+sqrt(mu * 1/L));
            elseif strcmp(func_type, 'basic')
                new_alpha = (1+sqrt(4 * alpha^2 + 1))/2;
                beta = (alpha-1)/new_alpha;
                alpha = new_alpha;
            else
            end
            Y = A + beta * (A - A_prev);

            A_prev = A;
            iter = iter + 1;
        end
    end

    A = Y;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function plotData (data_in,box,colordata, sizedata)
% PLOTDATA plots
% plotData (data_in,box,colordata, sizedata)
% data_in - N-by-2 matrix, each row is x-y coordinates of datapoint
% box - [xlim1 xlim2 ylim1 ylim2] - boundry of data
% colordata - color of ploted datapoints
% sizedata - size of ploted datapoints

if nargin<3
    colordata = 'b';
end

if nargin<4
    sizedata = 1;
end

if ~isempty(data_in)
    scatter(data_in(:,1), data_in(:,2),sizedata,colordata);

    if exist('box', 'var')
        if ~isempty(box)
            xlim([box(1), box(2)]); ylim([box(3) box(4)]);
        end
    end
end
set (gca, 'DataAspectRatio',[1 1 1]);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


function status = testPatches(k)
% TESTPATCHES Test suite for patches
%
% See also: ACID, matlab2tikz_acidtest

testfunction_handles = {
    @patch01;
    @patch02;
    @patch03;
    @patch04;
    @patch05;
    @patch06;
    @patch07;
    @patch08;
    };

numFunctions = length( testfunction_handles );

if nargin < 1 || isempty(k) || k <= 0
    status = testfunction_handles;
    return;  % This is used for querying numFunctions.

elseif (k<=numFunctions)
    status = testfunction_handles{k}();
    status.function = func2str(testfunction_handles{k});

else
    error('patchTests:outOfBounds', ...
        'Out of bounds (number of testfunctions=%d)', numFunctions);
end

end

% =========================================================================
function p = patch00()
% DO NOT INCLUDE IN ACID LIST
% Base patch plot for following tests
xdata = [2 2 0 2 5; 2 8 2 4 5; 8 8 2 4 8];
ydata = [4 4 4 2 0; 8 4 6 2 2; 4 0 4 0 0];
zdata = ones(3,5)*2;
p     = patch(xdata,ydata,zdata);
end
% =========================================================================
function stat = patch01()
stat.description = 'Set face color red';

p = patch00();
set(p,'FaceColor','r')
end
% =========================================================================
function stat = patch02()
stat.description = 'Flat face colors scaled in clim [0,40]';

p = patch00();
set(gca,'CLim',[0 40])
cdata = [15 30 25 2 60];
set(p,'FaceColor','flat','CData',cdata,'CDataMapping','scaled')
end
% =========================================================================
function stat = patch03()
stat.description = 'Flat face colors direct in clim [0,40]';

p = patch00();
set(gca,'CLim',[0 40])
cdata = [15 30 25 2 60];
set(p,'FaceColor','flat','CData',cdata,'CDataMapping','direct')
end
% =========================================================================
function stat = patch04()
stat.description = 'Flat face colors with 3D (truecolor) CData';

p = patch00();
cdata(:,:,1) = [0 0 1 0 0.8];
cdata(:,:,2) = [0 0 0 0 0.8];
cdata(:,:,3) = [1 1 1 0 0.8];
set(p,'FaceColor','flat','CData',cdata)
end
% =========================================================================
function stat = patch05()
stat.description = 'Flat face color, scaled edge colors in clim [0,40]';

p = patch00();
set(gca,'CLim',[0 40])
cdata = [15 30 25 2 60; 12 23 40 13 26; 24 8 1 65 42];
set(p,'FaceColor','flat','CData',cdata,'EdgeColor','flat','LineWidth',5,'CDataMapping','scaled')
end
% =========================================================================
function stat = patch06()
stat.description = 'Flat face color, direct edge colors in clim [0,40]';

p = patch00();
set(gca,'CLim',[0 40])
cdata = [15 30 25 2 60; 12 23 40 13 26; 24 8 1 65 42];
set(p,'FaceColor','flat','CData',cdata,'EdgeColor','flat','LineWidth',5,'CDataMapping','direct')
end
% =========================================================================
function stat = patch07()
stat.description = 'Flat face color with 3D CData and interp edge colors';

p = patch00();
cdata(:,:,1) = [0 0 1 0 0.8;
                0 0 1 0.2 0.6;
                0 1 0 0.4 1];
cdata(:,:,2) = [0 0 0 0 0.8;
                1 1 1 0.2 0.6;
                1 0 0 0.4 0];
cdata(:,:,3) = [1 1 1 0 0.8;
                0 1 0 0.2 0.6;
                1 0 1 0.4 0];
set(p,'FaceColor','flat','CData',cdata,'EdgeColor','interp','LineWidth',5)
end
% =========================================================================
function stat = patch08()
stat.description = 'Interp face colors, flat edges, scaled CData in clims [0,40]';

p = patch00();
set(gca,'CLim',[0 40])
cdata = [15 30 25 2 60; 12 23 40 13 26; 24 8 1 65 42];
set(p,'FaceColor','interp','CData',cdata,'EdgeColor','flat','LineWidth',5,'CDataMapping','scaled')
end
% =========================================================================


text:
clear;
addpath(genpath('../'));

output_dir = './CLM_res/';
if(~exist(output_dir, 'dir'))
    mkdir(output_dir)
end
%% select database and load bb initializations
db_root = 'E:\datasets\300VW\300VW_Dataset_2015_12_14/';
bb_root = './300VW_dets_mtcnn/';
extra_dir = 'E:\datasets\300VW\300VW_Dataset_2015_12_14\extra';
[ vid_locs, bboxes, gts_all, invalid_frames ] = CollectTestData(db_root, bb_root, extra_dir);

%% loading the patch experts and the PDM
[ patches, pdm, clmParams ] = Load_CLM_wild();

% Defining the view orientations that will be evaluated
%views = [0,0,0];
views = [0,0,0; 0,-30,0; 0,30,0; 0,-55,0; 0,55,0; 0,0,30; 0,0,-30];
views = views * pi/180;

num_points = numel(pdm.M) / 3;

verbose = true;

%% Select video
for i=1:numel(vid_locs)

    vid = VideoReader(vid_locs{i});

    bounding_boxes = bboxes{i};
    preds = [];
    n_frames = size(bounding_boxes,1);
    for f=1:n_frames
        input_image = readFrame(vid);

        %% Initialize from detected bounding boxes every 30 frames
        if (mod(f-1, 30) == 0)
            ind = min(f, size(bounding_boxes,1));
            bb = bounding_boxes(ind, :);
            % If no face detected use the closest detected BB
            if(bb(3) == 0)
               ind_next = ind + find(bounding_boxes(ind+1:end,3), 1);
               if(isempty(bounding_boxes(ind_next)) || bounding_boxes(ind_next,3)==0)
                   ind_next = find(bounding_boxes(1:ind,3), 1, 'last');
               end
               bb = bounding_boxes(ind_next, :);
            end
            bb(3) = bb(1) + bb(3);
            bb(4) = bb(2) + bb(4);
            reset = true;
        else
            reset = false;
        end

        % have a multi-view version for initialization, otherwise use
        % previous shape
        if(reset)
            clmParams.window_size = [25,25; 23,23; 21,21; 21,21];
            clmParams.numPatchIters = 4;
            clmParams.startScale = 1;

            [shape,g_param,l_param,lhood,lmark_lhood,view_used] = Fitting_from_bb_multi_hyp(input_image, [], bb, pdm, patches, clmParams, views);

        else
            clmParams.window_size = [23,23; 21,21; 19,19; 17,17];
            clmParams.numPatchIters = 3;
            clmParams.startScale = 2;
            [shape,g_param,l_param,lhood,lmark_lhood,view_used] = Fitting_from_bb(input_image, [], bb, pdm, patches, clmParams, 'gparam', g_param, 'lparam', l_param);
        end

        preds = cat(3, preds, shape);

        %% plot the result
        imshow(input_image);
        hold on;
        plot(shape(:,1), shape(:,2), '.r');
%         rectangle('Position', [bb(2), bb(1), bb(4), bb(3)]);
        hold off;
        drawnow expose

    end

    %% Grab the ground truth
    gt_landmarks = gts_all{i};

    % Remove unreliable frames
    if(~isempty(invalid_frames{i}))
        preds(:,:,int32(invalid_frames{i}))=[];
    end

    if(size(gt_landmarks,3) ~= size(preds,3))
        fprintf('something went wrong with vid %d\n', i);
    end

    [vid_name,~,~] = fileparts(vid_locs{i});
    [~,vid_name,~] = fileparts(vid_name);
    save([output_dir, '/', vid_name], 'preds', 'gt_landmarks');
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [date_time] = TimestampToString(t)

% TimestampToString - convert a UNIX timestamp to a human-readable string
%
% [date_time] = TimestampToString(t)
%
% INPUTS:
%   t: UNIX timestamp
%
% OUTPUTS:
%   date_time: string representing t in the format yyyy/mm/dd HH:MM:SS.FFF

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2016 University of Oxford
% Authors:
%  Geoff Pascoe (gmp@robots.ox.ac.uk)
%  Will Maddern (wm@robots.ox.ac.uk)
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  microseconds_per_day = 864e8; % 60 * 60 * 24 * 1000000
  unix_epoch = datenum('1970', 'yyyy');
  date_time = datestr(t / microseconds_per_day + unix_epoch, ...
    'yyyy/mm/dd HH:MM:SS.FFF');

end

author: ori-mrg

repo: robotcar-dataset-sdk


text:
function [ value, jsr ] = shr3 ( jsr )

%*****************************************************************************80
%
%% SHR3 evaluates the SHR3 generator for integers.
%
%  Discussion:
%
%    The input and output values are all unsigned 32 bit integers.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2008
%
%  Author:
%
%    Original C version by George Marsaglia, Wai Wan Tsang
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    George Marsaglia, Wai Wan Tsang,
%    The Ziggurat Method for Generating Random Variables,
%    Journal of Statistical Software,
%    Volume 5, Number 8, October 2000, seven pages.
%
%  Parameters:
%
%    Input, uint32 JSR, the seed.
%
%    Output, uint32 SHR3, the value of the SHR3 generator.
%
%    Output, uint32 JSR, the updated seed.
%
  jsr_input = uint32 ( jsr );

  jsr = uint32 ( bitxor ( jsr, uint32 ( bitshift ( jsr, + 13 ) ) ) );
  jsr = uint32 ( bitxor ( jsr, uint32 ( bitshift ( jsr, - 17 ) ) ) );
  jsr = uint32 ( bitxor ( jsr, uint32 ( bitshift ( jsr, +  5 ) ) ) );
%
%  We simply want to add JSR_INPUT and JSR and ignore the overflow digits.
%  But if we add two UINT32's and overflow, MATLAB sets the result to the
%  maximum 32 bit integer!
%
%  To avoid this, until I get smarter, we convert both UINT32's to doubles,
%  add them, subtract the maximum 32 bit integer + 1 (the overflow digit),
%  and convert back.  How ghastly.
%
  goof = double ( jsr_input ) + double ( jsr );
  if ( 4294967296 <= goof )
    goof = goof - 4294967296;
  end
  value = uint32 ( goof );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear;

addpath('../PDM_helpers/');

%% Create the PDM using PCA, from the recovered 3D data
clear
load('Torr_wild');

% need to still perform procrustes though

x = P3(1:end/3,:);
y = P3(end/3+1:2*end/3,:);

% To make sure that PDM faces the right way (positive Z towards the screen)
z = -P3(2*end/3+1:end,:);

inner_ids = 18:68;

%%
[ normX, normY, normZ, meanShape, Transform ] = ProcrustesAnalysis3D(x,y,z, true);

normX = normX(:, inner_ids);
normY = normY(:, inner_ids);
normZ = normZ(:, inner_ids);
meanShape = meanShape(inner_ids,:);
% make it zero centered now?
offsets = mean(meanShape);

normX = normX - offsets(1);
normY = normY - offsets(1);
normZ = normZ - offsets(1);

observations = [normX normY normZ];

[princComp, score, eigenVals] = princomp(observations,'econ');
% Keeping most variation
totalSum = sum(eigenVals);
count = numel(eigenVals);
for i=1:numel(eigenVals)
   if ((sum(eigenVals(1:i)) / totalSum) >= 0.999)
      count = i;
      break;
   end
end

load('pdm_68_multi_pie.mat', 'M');
M_full = M;

V = princComp(:,1:count);
E = eigenVals(1:count);
M = meanShape(:);

% Now normalise it to have actual world scale
lPupil = [(M(20) + M(23))/2; (M(20 + 51) + M(23 + 51))/2];
rPupil = [(M(26) + M(29))/2; (M(26 + 51) + M(29 + 51))/2];

dist = norm(lPupil - rPupil,2);

% average human interocular distance is 65mm
scaling = 65 / dist;

% normalise the mean values as well
M(1:end/3) = M(1:end/3) - mean(M(1:end/3));
M(end/3+1:2*end/3) = M(end/3+1:2*end/3) - mean(M(end/3+1:2*end/3));
M(2*end/3+1:end) = M(2*end/3+1:end) - mean(M(2*end/3+1:end));

M = M * scaling;
E = E * scaling .^ 2;

% orthonormalise V
scalingV = sqrt(sum(V.^2));
V = V ./ repmat(scalingV, numel(M),1);

E = E .* (scalingV') .^ 2;

%% Also align the model to a Multi-PIE one for accurate head orientation
% also align the mean shape model (an aligned 68 point PDM from Multi-PIE dataset)
M_inner = M;

load('pdm_68_multi_pie.mat', 'M');
M_align = M;

n_points = numel(M)/3;
M_inner_full = reshape(M_align, n_points,3);
M_inner_full = M_inner_full(inner_ids, :);

n_points_mouth = numel(inner_ids);

M_inner = reshape(M_inner, n_points_mouth, 3);

% work out the translation and rotation needed
[ R, T ] = AlignShapesKabsch(M_inner, M_inner_full);

% Transform the wild one to be in the same reference as the Multi-PIE one
M_aligned = (R * M_inner')';

M_aligned = M_aligned(:);
V_aligned = V;

% need to align the principal components as well
for i=1:size(V,2)

    V_aligned_pie_curr = (R * reshape(V(:,i), n_points_mouth, 3)')';
    V_aligned(:,i) = V_aligned_pie_curr(:);

end

V = V_aligned;
M = M_aligned;

save('pdm_51_inner.mat', 'E', 'M', 'V');

writePDM(V, E, M, 'pdm_51_inner.txt');

% also save this to model location
if(exist('../../models/pdm/', 'file'))
    save('../../models/pdm/pdm_51_inner.mat', 'E', 'M', 'V');
end


text:
function [outC,ptListC] = stackHDF5Files(outPath,passedScanDim)
% stackHDF5Files.m
%
% Reads .H5 files with mask slices and returns 3D stacks
%
% AI 9/13/19
%--------------------------------------------------------------------------
%INPUTS:
% outPath       : Path to generated H5 files
%                 Note: Assumes output filenames are of the form: prefix_slice# if
%                 passedScanDim = '2D' and of the form prefix_3D if passedScanDim = '3D'.
%------------------------------------------------------------------------
% RKP 9/18/2019 Updates for compatibility with testing pipeline

dirS = dir(fullfile(outPath,'outputH5','*.h5'));
fileNameC = {dirS.name};
ptListC = unique(strtok(fileNameC,'_'));
outC = cell(length(ptListC),1);

for p = 1:length(ptListC)

    %Get mask filenames for each pt
    matchIdxV = find(strcmp(strtok(fileNameC,'_'),ptListC{p}));

    switch(passedScanDim)
        case '3D'

            fileName = fullfile(outPath,'outputH5',fileNameC{matchIdxV});
            mask3M = h5read(fileName,'/mask');
            mask3M = permute(mask3M,[3 2 1]);

        case '2D'
            %Stack files
            mask3M = [];
            for s = 1: length(matchIdxV)

                slcName = fullfile(outPath,'outputH5',...
                          fileNameC{matchIdxV(s)});
                idx = strfind(slcName,'_slice');
                slcNum = str2double(slcName(idx+7:end-3));
                labelM = h5read(slcName,'/mask').';
                mask3M(:,:,slcNum) = labelM;

            end
    end

    outC{p} = mask3M;

end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [mu, U, w, likelihood] = prtUtilPenalizedIrls(y,cPhi,mu,A)
%% [mu, U, w, likelihood] = prtUtilPenalizedIrls(y,cPhi,mu,A)

nMaxIterations = 100; % Max outer iterations
gThreshold = 1e-6; % Termination criterion
nMaxStepSizeIterations = 20; % Per outer iteration how many steps to take
                             % Steps decrease by a factor of 1/2 each time

% Get initial error just so we know where we started
[dataNegativeLogLikelihood, yHat] = modelError(y,cPhi,mu);
weightNegativeLogLikelihood = (diag(A)'*(mu.^2))/2; % mu'*A*mu
oldTotalNegativeLogLikelihood = dataNegativeLogLikelihood + weightNegativeLogLikelihood;

totalNegativeLogLikelihood = oldTotalNegativeLogLikelihood;

for iteration = 1:nMaxIterations
    % Calculate find direction to go in

    % Error gradient
    e = (y - yHat); % Eq. 10 of Nabney, 1999
    g = cPhi'*e - diag(A).*mu; % In paragraph above Eq. 10 of Nabney, 1999

    % Compute the likelihood-dependent analogue of the noise precision.
    w = yHat.*(1-yHat); % In paragraph above Eq. 10 of Nabney, 1999

    % Compute the Hessian
    cPhiBeta = bsxfun(@times,cPhi,w);
    H = (cPhiBeta'*cPhi + A);

    % Invert Hessian
    [U, pdErr] = chol(H);
    % Make sure its positive definite
    if pdErr
        error('Ill conditioned hessian. Consider modifying your basis')
    end

    % Check convergence
    if all(abs(g)<gThreshold)
        break
    end
    % Use chol decomp trick to invert nicely and calculate direction
    hessianInvDir = U \ (U' \ g); % inv(H)*g

    stepSize = 1;
    for stepIteration = 1:nMaxStepSizeIterations

        muStar = mu + stepSize*hessianInvDir;

        % Compute total error (negative log likelihood
        [dataNegativeLogLikelihood, yHat] = modelError(y,cPhi,muStar);
        weightNegativeLogLikelihood = (diag(A)'*(muStar.^2))/2;

        totalNegativeLogLikelihood = dataNegativeLogLikelihood + weightNegativeLogLikelihood;

        if totalNegativeLogLikelihood < oldTotalNegativeLogLikelihood
            break
        else
            % We made things worse so descrease the step size.
            stepSize = stepSize/2;
        end
    end

    if stepIteration == nMaxStepSizeIterations
        % Went to the max on the inner loop means we are around numerical
        % precision. So exit
        %warning('bad exit : %0.4f',totalNegativeLogLikelihood-oldTotalNegativeLogLikelihood)
        totalNegativeLogLikelihood = oldTotalNegativeLogLikelihood;
        break
    else
        % Keep going
        mu = muStar;
        oldTotalNegativeLogLikelihood = totalNegativeLogLikelihood;
    end
end

likelihood = -totalNegativeLogLikelihood;

function [negativeLogLikelihood, yHat] = modelError(y,cPhi,mu)

yHat = 1 ./ (1+exp(-cPhi*mu));

yHat0 = (yHat==0);
yHat1 = (yHat==1);
negativeLogLikelihood = -(y(~yHat0)'*log(yHat(~yHat0)) + (1-y(~yHat1))'*log(1-yHat(~yHat1)));

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function A = getval(OBJ,PROP)

   % Function to extract a cell array from an
   % object and convert it into an array of values.
   % Basically it applies the {} operator to all
   % values in the array.
   %
   % NOTE: This routine is only valid for properties that returns a single
   % value per waveform.
   %

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   A = [];
   Ao = get(OBJ,PROP);
   for i = 1:length(Ao)
      A = cat(1,A,Ao{i});
   end;
   %A = A';
end

author: geoscience-community-codes

repo: GISMO


text:
% Chapter 6 - Controlling Chaos.
% Program_6b - Controlling Chaos in the Henon Map.
% Copyright Birkhauser 2013. Stephen Lynch.

% Control to a period-1 orbit (Figure 6.6(a)).
clear
alpha=1.2;beta=0.4;
k1=-1.8;k2=1.2;              % The regulator poles.
xstar=0.8358;ystar=xstar;    % The point to be stabilized.
N=400;x=zeros(1,N);y=zeros(1,N);rsqr=zeros(1,N);
x(1)=0.5;y(1)=0.6;           % Initial values.
rsqr(1)=(x(1))^2+(y(1))^2;
for n=1:N
    if n>198                 % Check point is in control region.
    x(n+1)=(-k1*(x(n)-xstar)-k2*(y(n)-ystar)+alpha)+beta*y(n)-(x(n))^2;
    y(n+1)=x(n);
    else
    x(n+1)=alpha+beta*y(n)-(x(n))^2;
    y(n+1)=x(n);
    end
    rsqr(n+1)=(x(n+1))^2+(y(n+1))^2;
end
hold on
axis([0 N 0 6])
plot(1:N,rsqr(1:N))
plot(1:N,rsqr(1:N),'o')
fsize=15;
set(gca,'XTick',0:50:N,'FontSize',fsize)
set(gca,'ytick',[0,6],'FontSize',fsize)
xlabel('n','FontSize',fsize)
ylabel('\it{r^2}','FontSize',fsize)
hold off

% End of Program_6b.

author: Sable

repo: mcbench-benchmarks


text:
% vec - vectorize an array
%
% Copyright(c) 2009-2011 Ryota Tomioka
% This software is distributed under the MIT license. See license.txt
function V=vec(M)
sz=size(M);
V=reshape(M, [prod(sz), 1]);

repo: REST


text:
function pairs = MatchPairs(list1,list2,varargin)

%MatchPairs - Pair nearest values in two lists.
%
%  USAGE
%
%    pairs = MatchPairs(list1,list2,<options>)
%
%    list1          first list
%    list2          second list
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'error'       maximum error (default = 0)
%    =========================================================================
%
%  OUTPUT
%
%    pairs        list of matched pairs (NaN indicates no match)
%
%  NOTE
%
%    This function uses list1 as the reference and is generally not symmetrical.

% ***************** TODO: Change PlotRippleStats to use Match instead of MatchPairs, then remove MatchPairs

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% 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.

% Default values
err = 0;

% Check number of parameters
if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help MatchPairs'' for details).');
end

% Check parameter sizes
if ~isdvector(list1),
	error('Parameter ''list1'' is not a vector (type ''help MatchPairs'' for details).');
end
if ~isdvector(list2),
	error('Parameter ''list2'' is not a vector (type ''help MatchPairs'' for details).');
end
if size(list2,1) == 1,
	list2 = list2';
end

% Parse parameter list
for i = 1:2:length(varargin),
  if ~ischar(varargin{i}),
    error(['Parameter ' num2str(i+2) ' is not a property (type ''help MatchPairs'' for details).']);
  end
  switch(lower(varargin{i})),
    case 'error',
      err = varargin{i+1};
      if ~isdscalar(err,'>0'),
        error('Incorrect value for property ''error'' (type ''help MatchPairs'' for details).');
      end
    otherwise,
      error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help MatchPairs'' for details).']);
  end
end

i = 1;
i2 = 1;
for i1 = 1:length(list1),
	% disp(['(' int2str(i1) ') ' num2str(list1(i1)) ' - (' int2str(i2) ') ' num2str(list2(i2)) ' ?']);
	while list1(i1) > list2(i2) + err,
		% disp([' (' int2str(i1) ') ' num2str(list1(i1)) ' > (' int2str(i2) ') ' num2str(list2(i2)) ' + ' num2str(err)]);
		pairs(i,:) = [NaN list2(i2)];
		i = i+1;
		i2 = i2+1;
	end
	if list1(i1) < list2(i2) - err,
		% disp([' (' int2str(i1) ') ' num2str(list1(i1)) ' < (' int2str(i2) ') ' num2str(list2(i2)) ' - ' num2str(err)]);
		pairs(i,:) = [list1(i1) NaN];
		i = i+1;
	else
		% disp([' (' int2str(i1) ') ' num2str(list1(i1)) ' ~ (' int2str(i2) ') ' num2str(list2(i2)) ' +- ' num2str(err)]);
		pairs(i,:) = [list1(i1) list2(i2)];
		i2 = i2+1;
		i = i+1;
	end
end
n = length(list2)-i2+1;
if n > 0,
	pairs = [pairs ; repmat(NaN,n,1) list2(i2:end)];
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%J.C. Spall, Jan. 2000
%This code implements SPSA with constraints for theta to lie in
%a specified hypercube (i.e., component-wise constraints).  Allows for multiple cases
%for purposes of statistical evaluation based on knowledge of true (noise-free) loss value
%(set cases=1 if user only wants one run).
clear all
global p z sigma;  %declaration of random var. (z) used for normal noise
                 %generation in loss fn. calls given seed above;
                 %also sigma in noise (noise may be dependent on theta)
p=10;
n=1000;						%total no. of loss measurements
cases=1;
alpha =.602;
gamma =.101;
a=.0017;
c=1.9;      				%chosen by standard guidelines
A=50;
lossfinalsq=0;          %variable for cum.(over 'cases')squared loss values
lossfinal=0;            %variable for cum. loss values
theta_lo=-1000*ones(p,1);   %lower bounds on theta
theta_hi=1000*ones(p,1);    %upper bounds on theta
theta_0=6*ones(p,1);
lossfinaleval='loss_test';  %choice of loss function for final perf. evaluations (noise-free)                            % evaluation
loss='loss_test_noise';     %loss function used in algorithm operations (with noise)
rand('seed',31415927)
randn('seed',3111113)
sigma=.1;       %multiplier in loss noise
feval(lossfinaleval,theta_0)
for i=1:cases
  theta=theta_0;
  for k=0:n-1
    ak = a/(k+1+A)^alpha;
    ck = c/(k+1)^gamma;
    delta = 2*round(rand(p,1))-1;
    thetaplus = theta + ck*delta;
    thetaminus = theta - ck*delta;
    yplus=feval(loss,thetaplus);
    yminus=feval(loss,thetaminus);
    ghat = (yplus - yminus)./(2*ck*delta);
    theta=theta-ak*ghat;
    % Two lines below invoke component-wise constraints
    theta=min(theta,theta_hi);
    theta=max(theta,theta_lo);
  end
  lossvalue=feval(lossfinaleval,theta);
  lossfinalsq=lossfinalsq+lossvalue^2;
  lossfinal=lossfinal+lossvalue;
end
% Display results: Mean loss value and standard deviation
%
disp('mean loss value over "cases" runs')
lossfinal/cases
%
if cases > 1
	disp('sample standard deviation of mean loss value')
   sd=((cases/(cases-1))^.5)*(lossfinalsq/cases-(lossfinal/cases)^2)^.5;
   sd=sd/(cases^.5)
else
end


text:
% switch on frame saving
saveframes=1;

% solve the advection equation: du/dt = c*du/dt
% numerical method: RK in time, right differencing
% du_m/dt = c*(u_{m+1}-u_m)/dt

% number of data points
M = 100;

% periodic interval dimensions (i.e. x \in [a,b))
a = 0; b =  2;

% data point spacing
dx = (b-a)/M;

% coordinates of data points ("nodes")
x = linspace(a, b-dx, M);

% start time
StartTime = 0;
EndTime = 8;

% dt chosen small to make space errors dominant
dt = 0.05*dx;

% choose order of RK integrator
s = 4;

% compute number of time steps
Nsteps = ceil((EndTime-StartTime)/dt);

% modify dt to end exactly at EndTime
dt = (EndTime-StartTime)/Nsteps;

% initial condition
u = exactsolution(x);

t = StartTime; % set time variable
fram = 1;
for n=1:Nsteps % begin time stepping

    utilde = u;  % initiate RK loop

    for rk=s:-1:1 % begin RK sub stages

        up1 = [utilde(2:M),utilde(1)];       % u_{m+1}

        % deltaplus
        deltaplus = ...
            (up1-utilde)/(dx);

        % update RK variable
        utilde = u+(dt/rk)*deltaplus;

    end % end RK sub stages

    u = utilde; % finish RK step
    t  = t+dt;   % update time

    if(mod(n,100)==0 |n==Nsteps) % selective plotting

        plot(x, u, 'r-*', 'LineWidth', 1); % plot numerical solution
        hold on;
        xplot = linspace(a, b, 1000); % plot exact solution at lots of points
        xmod = mod(xplot+t,b); % coordinate for exact pulse
        uexact = exactsolution(xmod);
        plot(xplot, uexact, 'k-', 'LineWidth', 1);hold off;
        axis([a b -.5 1.5]);
        legend(sprintf('Numerical solution (M=%d,t=%g)', M,t), 'Exact solution');
        drawnow; pause(0.02);
        if(saveframes == 1)        % save frame to file
            if(fram<10)            fname = sprintf('anim_000%d.ppm', fram);
            elseif(fram<100)   fname = sprintf('anim_00%d.ppm', fram);
            elseif(fram<1000) fname = sprintf('anim_0%d.ppm', fram);
            end
            print('-dppm', fname);
            fram = fram+1;
        end
    end
end
% to make the animated gif: I used ImageMagick's convert
% under cygwin:   convert *ppm deltaplus.gif
xmod = mod(x+t,b);
uexact = exactsolution(xmod);
finalerror = abs(uexact-u); % relative error

repo: Aero-matlab


text:
function model = generateGrRules(model, hyperlinkCommand)
% Generate model.grRules consistent with model.rules and model.genes, with
% an additional function of adding a field *.grRulesLinked in which the
% genes in the rules are hyperlinked with the provided command when
% displayed in the Matlab command window (used by `surfNet.m`).
%
% USAGE:
%    model = generateGrRules(model, hyperlinkCommand)
%
% INPUT:
%    model:            COBRA model
%
% OPTIONAL INPUT:
%    hyperlinkCommand: a string of Matlab command for the genes in the model
%                      e.g., 'fprintf(''%s'')' where %s will be replaced by the gene
%                      (default [], not adding the field for hyperlinked grRules)
%
% OUTPUT:
%    model:            COBRA model updated with *.grRules, and
%                      *.grRulesLinked if `hyperlinkCommand` is a non-empty string

if nargin < 2
    hyperlinkCommand = [];
end
% add *.grRulesLinked or not
addGrRulesLinked = ~isempty(hyperlinkCommand);
% add empty rules if nothing is there
if ~isfield(model, 'rules') && ~isfield(model, 'grRules')
    model.rules = repmat({''}, numel(model.rxns), 1);
end
if ~isfield(model, 'rules') && isfield(model, 'grRules')
    % if *.rules not exist, get from *.grRules
    model = generateRules(model);
end
% if *.genes is not present, just use x(1), ..., x(N) in *.rules as genes
genesFromRules = false;
if ~isfield(model, 'genes') && isfield(model, 'rules')
    genesFromRules = true;
    maxGeneId = 0;
end
% generate model.grRules from model.rules and model.genes
model.grRules = strrep(model.rules, '|', 'or');
model.grRules = strrep(model.grRules, '&', 'and');
model.grRules = strrep(model.grRules, '~', 'not');

if addGrRulesLinked
    model.grRulesLinked = model.grRules;
end
for j = 1:numel(model.grRules)
    % identify the gene ID
    re = regexp(model.grRules{j}, 'x\((\d+)\)', 'tokens');
    if ~isempty(re)
        for k = 1:numel(re)
            id = str2double(re{k}{1});
            if genesFromRules
                % if *.genes not exist, name the gene as x(n) and get the max number of genes
                gene = ['x(', re{k}{1}, ')'];
                maxGeneId = max([maxGeneId, id]);
            else
                % replace x(n) with the corresponding gene in *.genes
                gene = strtrim(regexprep(model.genes{id}, '\s', ''));
                model.grRules{j} = strrep(model.grRules{j}, ['x(', re{k}{1}, ')'], gene);
            end
            if addGrRulesLinked
                % add hyperlinks for the genes in *.grRulesLinked
                linkedGene = printHyperlink(sprintf(hyperlinkCommand, gene), gene, 0, 0);
                model.grRulesLinked{j} = strrep(model.grRulesLinked{j}, ['x(', re{k}{1}, ')'], linkedGene);
            end
        end
    end
end
if genesFromRules
    % assign *.genes if not exist
    if maxGeneId > 0
        model.genes = strcat('x(', cellfun(@num2str, num2cell((1:maxGeneId)'), 'UniformOutput', false), ')');
    else
        model.genes = cell(0, 1);
    end
end
end

author: opencobra

repo: cobratoolbox


text:
function linpack_s_test08 ( )

%*****************************************************************************80
%
%% TEST08 calls SGECO and SGESL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Local Parameters:
%
%    LDA defines the maximum matrix size we will use.
%
  lda = 10;
  n = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST08\n' );
  fprintf ( 1, '  For a general matrix,\n' );
  fprintf ( 1, '  SGECO computes the LU factors and computes\n' );
  fprintf ( 1, '  its reciprocal condition number;\n' );
  fprintf ( 1, '  SGESL solves a factored linear system.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Set the values of the matrix A.
%
  a(1:3,1:3) = [ ...
    1.0, 2.0, 3.0; ...
    4.0, 5.0, 6.0; ...
    7.0, 8.0, 0.0 ];
%
%  Factor the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix.\n' );

  [ a, ipvt, rcond, z ] = sgeco ( a, lda, n );

  fprintf ( 1, '  The reciprocal matrix condition number = %f\n', rcond );

  if ( rcond + 1.0 == 1.0 )
    fprintf ( 1, '  Error!  The matrix is nearly singular!\n' );
    return
  end
%
%  Set a right hand side.
%
  b(1:3) = [ 6.0, 15.0, 15.0 ]';
%
%  Solve the linear system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the linear system.\n' );

  job = 0;
  b = sgesl ( a, lda, n, ipvt, b, job );
%
%  Print the results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution returned by SGESL\n' );
  fprintf ( 1, '  (Should be (1,1,1))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %14f\n', b(i) );
  end
%
%  A second right hand side can be solved without refactoring a.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Call SGESL for a new right hand\n' );
  fprintf ( 1, '  side for the same, factored matrix.\n' );
%
%  Set the right hand side.
%
  b(1:3) = [ 1.0, 4.0, 7.0 ]';
%
%  Solve the system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve a linear system.\n' );

  job = 0;
  b = sgesl ( a, lda, n, ipvt, b, job );
%
%  Print the results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution returned by SGESL\n' );
  fprintf ( 1, '  (should be (1,0,0))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %14f\n', b(i) );
  end
%
%  The transposed problem A'*X = B can be solved by SGESL
%  as well, without any refactoring.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Call SGESL for transposed problem.\n' );
%
%  Set the right hand side.
%
  b(1:3) = [ 6.0, 6.0, -3.0 ]';
%
%  Solve the transposed system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Call SGESL to solve a transposed linear system.\n' );

  job = 1;
  b = sgesl ( a, lda, n, ipvt, b, job );
%
%  Print the results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution returned by SGESL\n' );
  fprintf ( 1, '  (should be (-1,0,1))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %14f\n', b(i) );
  end

  return
end


text:
function [errorCode,ima_stack,metastructs] = pco_sdk_example_read(imacount,segment,exposure_time,triggermode)
% set variables and grab images to internal memory and readout afterwards to a Matlab array
%
%   [ima_stackerrorCode] = pco_sdk_example_read(imacount,segment,exposure_time,triggermode)
%
% * Input parameters :
%    imacount                number of images to grab
%    segment                 segment to use for readout (default=1)
%    exposure_time           camera exposure time (default=10ms)
%    triggermode             camera trigger mode (default=AUTO)
%
% * Output parameters :
%    ima_stack               stack with grabbed images
%    errorCode               ErrorCode returned from pco.camera SDK-functions
%
%grab images into internal camera memory
%and read the images to an image stack
%
%function workflow
%open camera
%set variables
%setup internal memory
%start camera
%grab images
%stop camera
%read images from internal memory to image stack
%close camera
%

glvar=struct('do_libunload',0,'do_close',0,'camera_open',0,'out_ptr',[]);

if(~exist('imacount','var'))
 imacount = 10;
end

if(~exist('exposure_time','var'))
 exposure_time = 10;
end

if(~exist('triggermode','var'))
 triggermode = 0;
end

if(~exist('segment','var'))
 segment = uint16(1);
else
 if((segment<1)||(segment>4))
  segment=1;
 end
end

pco_camera_load_defines();
subfunc=pco_camera_subfunction();
subfunc.fh_lasterr(0);

[errorCode,glvar]=pco_camera_open_close(glvar);
pco_errdisp('pco_camera_setup',errorCode);
disp(['camera_open should be 1 is ',int2str(glvar.camera_open)]);
if(errorCode~=PCO_NOERROR)
 commandwindow;
 return;
end

out_ptr=glvar.out_ptr;

try

subfunc.fh_stop_camera(out_ptr);

pco_camera_resetmem(glvar);

cam_desc=libstruct('PCO_Description');
set(cam_desc,'wSize',cam_desc.structsize);
[errorCode,~,cam_desc] = calllib('PCO_CAM_SDK', 'PCO_GetCameraDescription', out_ptr,cam_desc);
pco_errdisp('PCO_GetCameraDescription',errorCode);

if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_NO_TIMESTAMP)==0)
 subfunc.fh_enable_timestamp(out_ptr,TIMESTAMP_MODE_BINARYANDASCII);
end

%enable MetaData if available
if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_METADATA))
 subfunc.fh_set_metadata_mode(out_ptr,1);
end

subfunc.fh_set_bitalignment(out_ptr,BIT_ALIGNMENT_MSB);
subfunc.fh_set_exposure_times(out_ptr,exposure_time,2,0,2);
subfunc.fh_set_triggermode(out_ptr,triggermode);

errorCode = calllib('PCO_CAM_SDK', 'PCO_ArmCamera', out_ptr);
pco_errdisp('PCO_ArmCamera',errorCode);

%if PCO_ArmCamera does fail no images can be grabbed
if(errorCode~=PCO_NOERROR)
 commandwindow;
 return;
end

subfunc.fh_get_triggermode(out_ptr);
subfunc.fh_set_transferparameter(out_ptr);
subfunc.fh_show_frametime(out_ptr);
ts_bin=subfunc.fh_is_binary_timestamp_enabled(out_ptr);

ima_stack=[];
%get images
[errorCode,imacount]=pco_camera_recmem(imacount,segment,glvar);

if(errorCode==0)
 [errorCode,ima_stack,metadata_stack]=pco_camera_readmem(imacount,1,segment,glvar);
end

if(errorCode==0)
 [~,~,count]=size(ima_stack);

 if(~isempty(metadata_stack))
  libmeta=libstruct('PCO_METADATA_STRUCT');
  [metadatasize,count_m]=size(metadata_stack);
  metastructs(count_m)=get(libmeta);
  for n=1:count_m
   [errorCode_m,metastructs(n)]=subfunc.fh_get_struct_metadata(metadata_stack(:,n),metadatasize);
   if(errorCode_m)
    pco_errdisp('get_struct_metadata',errorCode_m);
    metadata_stack=[];
    break;
   end
  end
  clear libmeta;
 end

 if(count==1)
  m=max(max(ima_stack(10:end-10,10:end-10)));
  disp(['image done maxvalue: ',int2str(m)]);
  if(ts_bin)
   txt=subfunc.fh_print_timestamp_t(ima_stack,1,16);
   disp(['Timestamp data of image: ',txt]);
  end
  if(~isempty(metadata_stack))
   txt=subfunc.fh_print_meta_timestamp(metastructs(1));
   disp(['Timestamp of meta (',num2str(1,'%04d'),'):  ',txt]);
   subfunc.fh_print_meta_struct(metastructs(1));
  end
 else
  disp([int2str(count),' images done']);
  reply = input('Show timestamps? Y/N [Y]: ', 's');
  if((isempty(reply))||(reply(1)== 'Y')||(reply(1)=='y'))
   for n=1:count
    if(ts_bin)
     txt=subfunc.fh_print_timestamp_t(ima_stack(:,:,n),1,16);
     disp(['Timestamp data of image(',num2str(n,'%04d'),'): ',txt]);
    end
    if(~isempty(metadata_stack))
     txt=subfunc.fh_print_meta_timestamp(metastructs(n));
     disp(['Timestamp data of meta (',num2str(n,'%04d'),'): ',txt]);
    end
   end
  end
 end
end

catch ME
 errorCode=subfunc.fh_lasterr();
 txt=blanks(101);
 txt=calllib('PCO_CAM_SDK','PCO_GetErrorTextSDK',pco_uint32err(errorCode),txt,100);

 clearvars -except ME glvar errorCode txt;

 if(glvar.camera_open==1)
  glvar.do_close=1;
  glvar.do_libunload=1;
  pco_camera_open_close(glvar);
 end

 if strfind(ME.identifier,'PCO_ERROR:')
  msg=[ME.identifier,' ',ME.message];
  disp(txt);
  warning('off','backtrace')
  warning(msg)
  for k=1:length(ME.stack)
    disp(['from file ',ME.stack(k).file,' at line ',num2str(ME.stack(k).line)]);
  end
  close();
  clearvars -except errorCode;
  commandwindow;
  return;
 else
  close();
  clearvars -except ME;
  rethrow(ME)
 end
end

pco_camera_resetmem(glvar);

if(glvar.camera_open==1)
 glvar.do_close=1;
 glvar.do_libunload=1;
 pco_camera_open_close(glvar);
end

if(exist('metastructs','var'))
 [~,count]=size(metastructs);
 if(count>3)
  subfunc=pco_camera_subfunction();
  n=3;
  txt=subfunc.fh_print_meta_timestamp(metastructs(n));
  disp(['Timestamp of meta(',num2str(n,'%04d'),'):  ',txt]);
  subfunc.fh_print_meta_struct(metastructs(n));
  clearvars subfunc;
 end
end

clear glvar;
commandwindow;
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function fd1d_heat_implicit_test02 ( )

%*****************************************************************************80
%
%% FD1D_HEAT_IMPLICIT_TEST02 does a problem with known solution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_HEAT_IMPLICIT_TEST02:\n' );
  fprintf ( 1, '  Compute an approximate solution to the time-dependent\n' );
  fprintf ( 1, '  one dimensional heat equation for a problem where we\n' );
  fprintf ( 1, '  know the exact solution.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    dH/dt - K * d2H/dx2 = f(x,t)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Run a simple test case.\n' );
%
%  Heat coefficient.
%
  k = k_test02 ( );
%
%  X_NUM is the number of equally spaced nodes to use between 0 and 1.
%
  x_num = 21;
  x_min = 0.0;
  x_max = 1.0;
  dx = ( x_max - x_min ) / ( x_num - 1 );
  x = linspace ( x_min, x_max, x_num );
%
%  T_NUM is the number of equally spaced time points between 0 and 10.0.
%
  t_num = 26;
  t_min = 0.0;
  t_max = 10.0;
  dt = ( t_max - t_min ) / ( t_num - 1 );
  t = linspace ( t_min, t_max, t_num );
%
%  Get the CFL coefficient.
%
  cfl = fd1d_heat_implicit_cfl ( k, t_num, t_min, t_max, x_num, x_min, x_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of X nodes = %d\n', x_num );
  fprintf ( 1, '  X interval is [%f,%f]\n', x_min, x_max );
  fprintf ( 1, '  X spacing is %f\n', dx );
  fprintf ( 1, '  Number of T values = %d\n', t_num );
  fprintf ( 1, '  T interval is [%f,%f]\n', t_min, t_max );
  fprintf ( 1, '  T spacing is %f\n', dt );
  fprintf ( 1, '  Constant K = %g\n', k );
  fprintf ( 1, '  CFL coefficient = %g\n', cfl );
%
%  Get the system matrix.
%
  a = fd1d_heat_implicit_matrix ( x_num, cfl );

  gmat = zeros ( x_num, t_num );
  hmat = zeros ( x_num, t_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Step            Time       RMS Error\n' );
  fprintf ( 1, '\n' );

  for j = 1 : t_num

    if ( j == 1 )
      h = ic_test02 ( x_num, x, t(j) );
      h = bc_test02 ( x_num, x, t(j), h );
    else
      h = fd1d_heat_implicit ( a, x_num, x, t(j-1), dt, cfl, @rhs_test02, @bc_test02, h );
    end

    g = exact_test02 ( x_num, x, t(j) );
    e = norm ( h(1:x_num) - g(1:x_num) ) / sqrt ( x_num );
    fprintf ( 1, '  %4d  %14.6g  %14.6g\n', j, t(j), e );
    gmat(1:x_num,j) = g(1:x_num);
    hmat(1:x_num,j) = h(1:x_num);

  end
%
%  Plot the data.
%
  figure ( 2 )
  [ tmat, xmat ] = meshgrid ( t, x );
  mesh ( tmat, xmat, hmat );
  title ( 'H(X,T) for TEST02 computed by FD1D\_HEAT\_IMPLICIT' );
  xlabel ( '<-- Time -->' );
  ylabel ( '<-- X -->' );
  zlabel ( '<--  H(X,T) -->' );
%
%  Write the data to files.
%
  r8mat_write ( 'g_test02.txt', x_num, t_num, gmat );
  r8mat_write ( 'h_test02.txt', x_num, t_num, hmat );
  r8vec_write ( 't_test02.txt', t_num, t );
  r8vec_write ( 'x_test02.txt', x_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  G(X,T) written to "g_test02.txt"\n' );
  fprintf ( 1, '  H(X,T) written to "h_test02.txt"\n' );
  fprintf ( 1, '  T values written to "t_test02.txt"\n' );
  fprintf ( 1, '  X values written to "x_test02.txt"\n' );
  return
end


text:
function [D]=distND(V1,V2)

D=zeros(size(V1,1),size(V2,1));

for q=1:1:size(V1,2) %For all dimensions
    A=V1(:,q); %Coordinates of first set
    B=V2(:,q); %Coordinates of second set

    AB=A(:,ones(1,size(B,1)));
    BA=B(:,ones(1,size(A,1)))';

    D=D+(AB-BA).^2;
end
D=sqrt(D);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function styles = get_vis_styles()

    colors = {'r','g','b','c','m','y'};
    markers = {'+','o','s', 'x','.'}; %,'-','s','d'};
    styles = {};
    j = 1;
    for m = 1:length(markers)
        for c = 1:length(colors)
            styles{j} = [colors{c} markers{m}];
            j = j + 1;
        end
    end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [shapeModelOpt] = learnDenseShape(statesDir,class,jobID)

%% Initializations
globals;
params=get_params();

%% Penalty coefficients and max allowed changes
lambda = params.opt.lambda;
maxDelta = params.opt.maxDelta;
stepsize = params.opt.trainStepSize;

%% Loading/Initializing using Visual Hull
shapeFile = fullfile(cachedir,class,sprintf('currentShape%s.mat',jobID));
initShapeFile = fullfile(cachedir,class,sprintf('InitShape%s.mat',jobID));

if(exist(initShapeFile,'file'))
    fprintf('Loading cached Visual Hull initialization from \n%s\n',initShapeFile);
    load(initShapeFile);
else
    [S,V,shapePriorMean,deformationPrior,shapePriorInstance,tri,normals] =...
        initializeShape(statesDir);
    initMesh = struct('vertices',S,'faces',tri);
    fnames = getFileNamesFromDirectory(statesDir,'types',{'.mat'},'mode','path');
    alpha = initializeAlpha(S,V,fnames);
    save(initShapeFile,'S','shapePriorMean','deformationPrior','shapePriorInstance',...
        'fnames','V','alpha','tri','normals','initMesh');
end

%% Iterating to learn shapes
fprintf('\nOptimizing shapes via block co-ordinate descentn\n');
fprintf('Iter\tKPGrad\tOccGrad\tSilGrad\tSIPrior\tShPrior\tDfPrior\tStep\n');
for iter = 1:params.opt.trainIters
    %disp(['Iteration ' int2str(iter)]);
    [gradS,gradV,gradAlpha] = refineShape(S,V,alpha,fnames,lambda,...
        shapePriorMean,deformationPrior,shapePriorInstance,iter);
    if(mod(iter,2))
        %% Updating shape & basis
        step = min([maxDelta(1)/max(max(abs(gradS))),maxDelta(2)...
            /max(max(abs(gradV))),stepsize]);
        fprintf('%3.5f\n',step);
        S = step*gradS + S;
        V = step*gradV + V;

        % Visualization
        if(mod(iter,100)==-1)
            close all;
            subplot(1,2,1);
            plot3(S(:,1),S(:,2),S(:,3),'b.');hold on;axis equal;
            subplot(1,2,2);
            trisurf(tri,S(:,1),S(:,2),S(:,3));axis equal;
            pause();close all;
        end
    else
        %% Normalizing alpha for each image
        mGrad = max(abs(gradAlpha),[],2);
        limGrad = min(mGrad,maxDelta(3)/stepsize);
        scale = limGrad./mGrad;
        gradAlpha = gradAlpha.*repmat(scale,1,size(gradAlpha,2));
        step = min(maxDelta(3)/max(max(abs(gradAlpha))),stepsize);
        fprintf('%3.5f\n',step);
        alpha = step*gradAlpha + alpha;
    end

    %% Renormalizing V
    for k=1:params.opt.numbasis
        den = norm(V(:,k))/sqrt(size(S,1));
        V(:,k)=V(:,k)/den;
        alpha(k,:)=den*alpha(k,:);
    end

end

save(shapeFile,'S','shapePriorMean','deformationPrior','shapePriorInstance',...
    'fnames','V','alpha','tri','normals','initMesh','params');
shapeModelOpt = load(shapeFile);
end

repo: CategoryShapes


text:
function [lb,ub,cand_rows] = find_lp_bounds_interval(F_struc,K)
%find_lp_bounds Internal function to extract upper and lower variable bounds

% special code for the interval data case, to avoid overhead in the default
% function find_lp_bounds

n = size(F_struc,2)-1;
lb = -inf*ones(n,1);
ub = inf*ones(n,1);
cand_rows = [];

if (K.f ~=0)
    A = -F_struc(1:K.f,2:end);
    b = F_struc(1:K.f,1);
    n = size(F_struc,2)-1;
    cand_rows = find(sum(A~=0,2)==1);
    for i = 1:length(cand_rows)
        j = find(A(cand_rows(i),:));
        ub(j)=min(ub(j),b(cand_rows(i))/A(cand_rows(i),j));
        lb(j)=max(lb(j),b(cand_rows(i))/A(cand_rows(i),j));
    end
end

if (K.l ~=0)
    A = -F_struc(K.f+1:K.f+K.l,2:end);
    b = F_struc(K.f+1:K.f+K.l,1);
    n = size(F_struc,2)-1;
    cand_rows = find(sum(A~=0,2)==1);
    for i = 1:length(cand_rows)
        j = find(A(cand_rows(i),:));
        if A(cand_rows(i),j)>0
            ub(j)=min(ub(j),sup(b(cand_rows(i))/A(cand_rows(i),j)));
        else
            lb(j)=max(lb(j),inf_(b(cand_rows(i))/A(cand_rows(i),j)));
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef Scale < nntest
  properties
    x
    a
    b
  end

  properties (TestParameter)
    dim = {1 2 3 4}
  end

  methods (TestClassSetup)
    function data(test,device)
      test.x = test.randn(15,14,3,2) ;
      test.a = test.randn(15,14,3,2) ;
      test.b = test.randn(15,14,3,2) ;
    end
  end

  methods (Test)
    function data_and_parameters(test, dim)
      x = test.x ;
      a = test.a ;
      b = test.b ;

      a = sum(a, dim) ;
      b = sum(b, dim) ;

      scale = dagnn.Scale('size', size(a), 'hasBias', true) ;

      output = scale.forward({x}, {a,b}) ;
      dzdy = test.randn(size(output{1})) ;
      [derInputs, derParams] = scale.backward({x}, {a,b}, {dzdy}) ;

      pick = @(x) x{1} ;
      dzdx = derInputs{1} ;
      dzda = derParams{1} ;
      dzdb = derParams{2} ;

      test.der(@(x) pick(scale.forward({x},{a,b})), x, dzdy, dzdx, 1e-2 * test.range) ;
      test.der(@(a) pick(scale.forward({x},{a,b})), a, dzdy, dzda, 1e-2 * test.range) ;
      test.der(@(b) pick(scale.forward({x},{a,b})), b, dzdy, dzdb, 1e-2 * test.range) ;
    end

    function data_only(test, dim)
      x = test.x ;
      a = test.a ;
      b = test.b ;

      a = sum(a, dim) ;
      b = sum(b, dim) ;

      scale = dagnn.Scale('size', size(a), 'hasBias', true) ;

      output = scale.forward({x,a,b}, {}) ;
      dzdy = test.randn(size(output{1})) ;
      [derInputs, derParams] = scale.backward({x,a,b}, {}, {dzdy}) ;

      pick = @(x) x{1} ;
      dzdx = derInputs{1} ;
      dzda = derInputs{2} ;
      dzdb = derInputs{3} ;

      test.der(@(x) pick(scale.forward({x,a,b},{})), x, dzdy, dzdx, 1e-2 * test.range) ;
      test.der(@(a) pick(scale.forward({x,a,b},{})), a, dzdy, dzda, 1e-2 * test.range) ;
      test.der(@(b) pick(scale.forward({x,a,b},{})), b, dzdy, dzdb, 1e-2 * test.range) ;
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function pixels = getpixels(data, n)

x = data(1);
y = data(2);
m = floor(n/2);

[i0, j0] = toindex(x, y);
r = repmat(-m:m, n, 1);
i = i0 + r;
j = j0 + r';
pixels = [i(:) j(:)];


text:
function [Y, eigvalue] = Eigenmap(W, ReducedDim, bEigs)
%function [Y, eigvalue] = Eigenmap(W, ReducedDim, bEigs)
%
%                   W       -  the affinity matrix.
%           ReducedDim      -  the dimensionality of the reduced subspace.
%                bEigs      -  whether to use eigs to speed up. If not
%                              specified, this function will automatically
%                              decide based on the size of W.
%
%   version 2.1 --November/2011
%   version 2.0 --May/2009
%   version 1.5 --Dec./2005
%   version 1.0 --Aug./2003
%
%   Written by Deng Cai (dengcai AT gmail.com)

MAX_MATRIX_SIZE = 1600; % You can change this number according your machine computational power
EIGVECTOR_RATIO = 0.1; % You can change this number according your machine computational power

[row,col] = size(W);
if row ~= col
    error('W must square matrix!!');
end

nSmp = row;

if ~exist('ReducedDim','var')
    ReducedDim = 10;
end
ReducedDim = min(ReducedDim+1,row);

D_mhalf = full(sum(W,2).^-.5);
D_mhalfMatrix = spdiags(D_mhalf,0,nSmp,nSmp);
W = D_mhalfMatrix*W*D_mhalfMatrix;

W = max(W,W');

dimMatrix = size(W,2);
if ~exist('bEigs','var')
    if (dimMatrix > MAX_MATRIX_SIZE && ReducedDim < dimMatrix*EIGVECTOR_RATIO)
        bEigs = 1;
    else
        bEigs = 0;
    end
end

if bEigs
    option = struct('disp',0);
    [Y, eigvalue] = eigs(W,ReducedDim,'la',option);
    eigvalue = diag(eigvalue);
else
    [Y, eigvalue] = eig(full(W));
    eigvalue = diag(eigvalue);

    [junk, index] = sort(-eigvalue);
    eigvalue = eigvalue(index);
    Y = Y(:,index);
    if ReducedDim < length(eigvalue)
        Y = Y(:, 1:ReducedDim);
        eigvalue = eigvalue(1:ReducedDim);
    end
end

eigIdx = find(abs(eigvalue) < 1e-6);
eigvalue (eigIdx) = [];
Y (:,eigIdx) = [];

nGotDim = length(eigvalue);

idx = 1;
while(abs(eigvalue(idx)-1) < 1e-12)
    idx = idx + 1;
    if idx > nGotDim
        break;
    end
end
idx = idx - 1;

if(idx > 1)  % more than one eigenvector of 1 eigenvalue
    u = zeros(size(Y,1),idx);

    d_m = 1./D_mhalf;
    cc = 1/norm(d_m);
    u(:,1) = cc./D_mhalf;

    bDone = 0;
    for i = 1:idx
        if abs(Y(:,i)' * u(:,1) - 1) < 1e-14
            Y(:,i) = Y(:,1);
            Y(:,1) = u(:,1);
            bDone = 1;
        end
    end

    if ~bDone
        for i = 2:idx
            u(:,i) = Y(:,i);
            for j= 1:i-1
                u(:,i) = u(:,i) - (u(:,j)' * Y(:,i))*u(:,j);
            end
            u(:,i) = u(:,i)/norm(u(:,i));
        end
        Y(:,1:idx) = u;
    end
end

Y = D_mhalfMatrix*Y;

Y(:,1) = [];
eigvalue(1) = [];

author: ZJULearning

repo: MatlabFunc