text:
function [ddCs, rhsCs] = electrodeConcentration(dCs,cs_barrato,T,jflux,param)
% electrodeConcentration describes the ODE of the concentration of lithium ions within the
% electrodes.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

if(param.SolidPhaseDiffusion==1 || param.SolidPhaseDiffusion==2)
    % Cathode
    rhsCs_p  =((-3/param.Rp_p)*jflux(1:param.Np));
    ddCs_p   = dCs(1:param.Np) - rhsCs_p;

    % Anode
    rhsCs_n  = ((-3/param.Rp_n)*jflux(param.Np+1:end));
    ddCs_n   = dCs(param.Np+1:end) - rhsCs_n;
else
    switch param.SolidPhaseDiffusionNumericalScheme
        % Use the FDM method for the solid phase diffusion
        case 1
            [rhsCs_p, rhsCs_n, ddCs_p, ddCs_n] = FDM9orderElectrodeDiffusion(T, cs_barrato, jflux, dCs, param);
        % Use the spectral method for the discretization of the solid phase
        % diffusion
        case 2
            [rhsCs_p, rhsCs_n, ddCs_p, ddCs_n] = spectralMethodElectrodeDiffusion(T, cs_barrato, jflux, dCs, param);
    end
end
rhsCs   = [rhsCs_p;rhsCs_n];
ddCs    = [ddCs_p;ddCs_n];
end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function pdf = english_word_length_pdf ( x )

%*****************************************************************************80
%
%% ENGLISH_WORD_LENGTH_PDF evaluates the English Word Length PDF.
%
%  Discussion:
%
%    PDF(A,B;X) = B(X) if 1 <= X <= A
%                = 0    otherwise
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 July 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Henry Kucera, Winthrop Francis,
%    Computational Analysis of Present-Day American English,
%    Brown University Press, 1967.
%
%  Parameters:
%
%    Input, integer X, the word length whose probability is desired.
%
%    Output, real PDF, the value of the PDF.
%
  word_length_max = 27;

  pdf_vec = [ ...
    0.03160, ...
    0.16975, ...
    0.21192, ...
    0.15678, ...
    0.10852, ...
    0.08524, ...
    0.07724, ...
    0.05623, ...
    0.04032, ...
    0.02766, ...
    0.01582, ...
    0.00917, ...
    0.00483, ...
    0.00262, ...
    0.00099, ...
    0.00050, ...
    0.00027, ...
    0.00022, ...
    0.00011, ...
    0.00006, ...
    0.00005, ...
    0.00002, ...
    0.00001, ...
    0.00001, ...
    0.00001, ...
    0.00001, ...
    0.00001 ];
  pdf_sum = 0.99997;

  if ( 1 <= x & x <= word_length_max )
    pdf = pdf_vec(x) / pdf_sum;
  else
    pdf = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [pvec, pstruct] = tapas_bayes_optimal_whichworld_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% This is a dummy function

pvec    = [];
pstruct = struct;

return;

author: translationalneuromodeling

repo: tapas


text:
function [ im_data ] = prepare_img_deng( im_data, mean_data, im_size_h,im_size_w )
%PREPARE_IMG Summary of this function goes here
%   Detailed explanation goes here
% Convert an image returned by Matlab's imread to im_data in caffe's data
% format: W x H x C with BGR channels
im_data = im_data(:, :, [3, 2, 1]);  % permute channels from RGB to BGR
im_data = permute(im_data, [2, 1, 3]);  % flip width and height
im_data = single(im_data);  % convert from uint8 to single
im_data = imresize(im_data, [im_size_w im_size_h ], 'bilinear');  % resize im_data
im_data = im_data - mean_data;  % subtract mean_data (already in W x H x C, BGR)
end

repo: Learning-via-Translation


text:
function save_vt_history_data(fileName, data)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     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/>.

    fileID = fopen(fileName,'w');
    len = length(data(:,1));
    for ind = 1 : len
        fprintf(fileID,'%d %d \r\n',ind,data(ind, :));
    end
    fclose(fileID);
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function [isSame, nDiff, commonFields] = isSameCobraModel(model1, model2, printLevel)
% Checks if two COBRA models are the same
%
% USAGE:
%
%    [isSame, nDiff, commonFields] = isSameCobraModel(model1, model2, printLevel)
%
% INPUTS:
%    model1:          COBRA model structure 1
%    model2:          COBRA model structure 2
%
% OPTIONAL INPUTS:
%    printLevel:      Whether to provide additional output (0 - no output (default), 1 - output).
%
% OUTPUTS:
%    isSame:          True if all common fields are identical, else false
%    nDiff:           Number of differences between the two models for each field
%    commonFields:    List of common fields
%
% .. Authors:
%     - Markus Herrgard 9/14/07
%     - CI integration: Laurent Heirendt

%TODO this function needs updating to use structeq.m

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

isSame = true;

fields1 = fieldnames(model1);
fields2 = fieldnames(model2);
onlyIn1 = setdiff(fields1, fields2);
onlyIn2 = setdiff(fields2, fields1);
if printLevel > 0
    if ~isempty(onlyIn1)
        fprintf('The following fields are only present in the first input model:\n');
        for i = 1:numel(onlyIn1)
            fprintf('%s\n',onlyIn1{i});
        end
    end
    if ~isempty(onlyIn2)
        fprintf('The following fields are only present in the second input model:\n');
        for i = 1:numel(onlyIn2)
            fprintf('%s\n',onlyIn2{i});
        end
    end
end

commonFields = intersect(fields1, fields2);
commonFields = commonFields(~strcmpi('description', commonFields));

if (~isempty(onlyIn1) || ~isempty(onlyIn2))
    isSame = false;
end

% initialize variables
nFields = length(commonFields);
nDiff = zeros(nFields,1);

% loop through all fields
for i = 1:nFields
    fieldName = commonFields{i};
    value1 = getfield(model1, fieldName);
    value2 = getfield(model2, fieldName);

    if 0 %debugging code
            if strcmp(fieldName,'rxnConfidenceScores')
                pause(0.1);
            end
    end
    % replace all whitespaces
    if iscellstr(value1)
        value1 = regexprep(value1, '[^\w'']', '');
        value2 = regexprep(value2, '[^\w'']', '');
    end

    if isnumeric(value1)
        nDiff(i) = sum(sum(~((value1 == value2) | (isnan(value1) & isnan(value2))) ));
    elseif iscellstr(value1)
        if 0 %debugging code
        for i=1:length(value1)
            if class(value1{i})~=class(value2{i})
                pause(0.1)
            end
            if length(value1{i})~=length(value2{i})
                pause(0.1)
            end
        end
        end
        nDiff(i) = sum(~strcmp(value1, value2));
    elseif ischar(value1)
        nDiff(i) = ~strcmp(value1, value2);
    end

    if (nDiff(i) > 0)
        isSame = false;
    end
    if printLevel > 0
        if nDiff(i) > 0
            fprintf('Field %s differs in %d positions between the models\n',fieldName,nDiff(i));
        end
    end

end

author: opencobra

repo: cobratoolbox


text:
function vwap = getVWAP(price, volume, dates)
%GETVWAP: calculate the Volume Weighted Average Price at the end of each
%day, given intra daily data of the closing price and volume.
%
%   VWAP = GETVWAP(price, volume, dates) returns the Volume Weighted
%   Average Price (VWAP) at the end of the day. The input consists of the
%   intra-daily price, volume and dates (in formatted form).
%
%   vwap:           is a vector of the VWAP prices at the end of each
%                   unique day, conditional on the dates.
%
%  $Date: 04/10/2012$
%
% -------------------------------------------------------------------------

% Not all securities publish their volume shares, i.e. sometimes the volume
% vector is empty.
if sum(volume) == 0, error('getVWAP:InvalidInput','No historical intra-daily data of volume'); end
if size(price,2) > 1 || size(volume,2) > 1 || size(dates,2) > 1, error('getVWAP:InvalidInput','Price, volume and dates should be a row vector'); end

% FIND THE UNIQUE DAYS. TWO OPTIONS ARE POSSIBLE:
% -- 1 --   make use of the included GETUNIQUEDAYELEMENTS code (general
%           framework, adaptable for multi purposes, external)

uniqueDays = getUniqueDayElements(dates);
k = size(unique(day(dates)),1);
vwap = zeros(1, k );

% -- 2 --   get rid of the [EXTERNAL] GETUNIQUEDAYELEMENTS code and use the
%           following [INTERNAL] method:

    % k = size(unique(day(dates)),1);
    % vwap = zeros(1, k );
    % uniqueDays = zeros(1, k ); iter = 1; uniqueDays(iter) = day(dates(1));
    %
    % % find the unique days
    % for i = 2:size(dates,1);
    %     if uniqueDays(iter) ~= day(dates(i));
    %         iter = iter + 1;
    %         uniqueDays(iter) = day(dates(i));
    %     end
    % end

% calculate vwap
for i = 1:size(uniqueDays,2)
    dayi = day(dates)==uniqueDays(i);
    vwap(i) = sum( volume(dayi) ...
        .*price( dayi ) ) / sum( volume( dayi ) );
end

author: Sable

repo: mcbench-benchmarks


text:
function [distAll, scoresAll, pointsAll] = assignGTinferenceMulti(keypointsAll,annolist,pidxsAll,parts,thresh_gt,marg_scores,nTrials)

if (nargin < 5)
    thresh_gt = 0.5;
end

scoresAll = cell(length(pidxsAll),1);
pointsAll = cell(length(pidxsAll),1);
distAll = cell(length(pidxsAll),1);

for i = 1:length(pidxsAll)
    scoresAll{i} = cell(length(annolist),1);
    pointsAll{i} = cell(length(annolist),1);
    distAll{i} = cell(length(annolist),1);
end

for imgidx = 1:length(annolist)

%     figure(200); clf; imagesc(imread(annolist(imgidx).image.name));
%     hold on; axis equal;

    det_c = {'g','b'};
    for trialidx = 1:nTrials/2
        if (isempty(keypointsAll(imgidx).det{trialidx}))
            continue;
        end
        dist = cell(2,1);
        for detidx = 1:2

            dist{detidx} = inf(length(pidxsAll),length(annolist(imgidx).annorect));
            for i = 1:length(pidxsAll)
                pidx = pidxsAll(i);
                % part is a joint
                assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
                jidx = parts(pidx+1).pos(1);

                det = keypointsAll(imgidx).det{trialidx}{jidx+1};
                assert(~isempty(det));
%                 [val,id] = max(det(:,4+detidx));
                id = detidx;
                val = det(id,4+detidx);
                pp = det(id,1:2);

                if (isnan(val))
                    val = -inf;
                end

                if (isnan(pp(1)))
                    assert(isnan(pp(2)));
                    pp = [0 0];
                end

                pointsAll{i}{imgidx} = [pointsAll{i}{imgidx}; pp];
                if (marg_scores)
                    scoresAll{i}{imgidx} = [scoresAll{i}{imgidx}; val];
                else
                    scoresAll{i}{imgidx} = [scoresAll{i}{imgidx}; det(id,3)];
                end

                for ridx = 1:length(annolist(imgidx).annorect)
                    rect = annolist(imgidx).annorect(ridx);
                    refDist = util_get_head_size(rect);
                    if (isfield(rect, 'annopoints') && isfield(rect.annopoints, 'point'))
                        p = util_get_annopoint_by_id(rect.annopoints.point, jidx);
                        if (~isempty(p))
                            dist{detidx}(i,ridx) = norm([p.x p.y] - pp)/refDist;
                            %                         plot(p.x,p.y,'ro','MarkerFaceColor','y','MarkerEdgeColor','k','MarkerSize',10);
                        else
                            dist{detidx}(i,ridx) = nan;
                        end
                    end
                end
                %             plot(pp(1),pp(2),'ro','MarkerFaceColor',det_c{detidx},'MarkerEdgeColor','k','MarkerSize',10);
            end
        end

        for detidx = 1:2
            pckSum = zeros(1,size(dist{detidx},2));
            idxsNotNaN = ~isnan(dist{detidx});
            for i=1:size(dist{detidx},2)
                pckSum(i) = sum(dist{detidx}(idxsNotNaN(:,i),i) <= thresh_gt);
            end
            [val,idx] = max(pckSum);
            idx2 = setdiff(1:size(dist{detidx},2),idx);
            dist{detidx}(idxsNotNaN(:,idx2),idx2) = inf;
            for i = 1:length(pidxsAll)
                distAll{i}{imgidx} = [distAll{i}{imgidx}; dist{detidx}(i,:)];
            end
        end
    end

end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function weight = olmar1_kernel(data, weight_o, epsilon, W)
% This program output the final portfolio the OLMAR-1 algorithm
% OLMAR-1 has only one experts, thus, we go expert directly.
% If we want buy and hold, combine the experts here.
%
% function [weight] = olmar1_expert(data, weight_o, epsilon, W);
%
% weight: final portfolio, used for next rebalance
%
% data: market sequence vectors
% weight_o: last portfolio
% epsilon: mean reversion threshold
% W: window size for calculating moving average
%
% Example: [weight] = olmar1_expert(data, weight_o, epsilon, W);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

weight = olmar1_expert(data, weight_o, epsilon, W);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function test_suite = test_partitions
% tests for partitioning functions
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_nfold_partitioner()
    ds=cosmo_synthetic_dataset('nchunks',5,'ntargets',4);

    p=cosmo_nfold_partitioner(ds);
    assertEqual(p, cosmo_nfold_partitioner(ds.sa.chunks));

    fns={'train_indices';'test_indices'};
    assertEqual(fns, fieldnames(p));
    for k=1:5
        test_indices=(k-1)*4+(1:4);
        train_indices=setdiff(1:20,test_indices);
        for j=1:2
            fn=fns{j};
            if j==1
                v=train_indices;
            else
                v=test_indices;
            end
            w=p.(fn);
            assertEqual(w{k}, v');
        end
    end

function test_nchoosek_partitioner()
    ds=cosmo_synthetic_dataset('nchunks',5,'ntargets',4);

    p=cosmo_nfold_partitioner(ds);
    q=cosmo_nchoosek_partitioner(ds,1);
    assertEqual(p,q);

    q=cosmo_nchoosek_partitioner(ds,.2);
    assertEqual(p,q);

    p=cosmo_nchoosek_partitioner(ds,3);
    q=cosmo_nchoosek_partitioner(ds,.6);

    assertEqual(p,q);
    assertFalse(isequal(p, cosmo_nchoosek_partitioner(ds,.4)));

    q2=cosmo_nchoosek_partitioner(ds.sa.chunks,3);
    assertEqual(p,q2);

    fns={'train_indices';'test_indices'};
    for j=1:2
        fn=fns{j};
        counts=zeros(20,1);

        v=p.(fn);
        assertEqual(size(v),[1 10]);

        for k=1:numel(v)
            w=v{k};
            counts(w)=counts(w)+1;
        end
        assertEqual(counts,ones(20,1)*j*2+2);
    end

function test_nchoosek_partitioner_half()
    offsets=[0 floor(rand()*10+20)];
    for offset=offsets
        for nchunks=2:4:10

            ds=cosmo_synthetic_dataset('nchunks',nchunks,'ntargets',3);
            ds.sa.chunks=ds.sa.chunks+offset;

            p=cosmo_nchoosek_partitioner(ds,'half');
            combis=nchoosek(1:nchunks,nchunks/2);

            n=size(combis,1);
            assertEqual(numel(p.train_indices),n/2);
            assertEqual(numel(p.test_indices),n/2);

            for k=1:n/2
                tr_idx=find(cosmo_match(ds.sa.chunks-offset,...
                                                        combis(k,:)));
                te_idx=find(cosmo_match(ds.sa.chunks-offset,...
                                                        combis(n-k+1,:)));

                assertEqual(p.train_indices{k},tr_idx);
                assertEqual(p.test_indices{k},te_idx);
            end
        end
    end

function test_nchoosek_partitioner_grouping()
    for nchunks=[2 5]
        ds=cosmo_synthetic_dataset('nchunks',nchunks,'ntargets',6);
        ds.sa.modality=mod(ds.sa.targets,2)+1;
        ds.sa.targets=floor(ds.sa.targets/2);

        for n_test=1:(nchunks-1)
            for moda_idx=1:4

                if moda_idx==3
                    modas={1 2};
                    moda_arg={1 2};
                elseif moda_idx==4
                    modas={1,2};
                    moda_arg=[];
                else
                    modas={moda_idx};
                    moda_arg=moda_idx;
                end

                n_moda=numel(modas);

                p=cosmo_nchoosek_partitioner(ds,n_test,'modality',...
                                                            moda_arg);
                combis=nchoosek(1:nchunks,n_test);
                n_combi=size(combis,1);

                n_folds=numel(p.train_indices);
                assertEqual(numel(p.test_indices),n_folds);
                assertEqual(n_folds,n_combi*n_moda);

                % each fold must be present exactly once
                visited_count=zeros(1,n_folds);
                for m=1:n_moda
                    for j=1:n_combi
                        tr_msk=~cosmo_match(ds.sa.chunks,combis(j,:)) & ...
                                    ~cosmo_match(ds.sa.modality,modas{m});
                        te_msk=cosmo_match(ds.sa.chunks,combis(j,:)) & ...
                                    cosmo_match(ds.sa.modality,modas{m});
                        tr_idx=find_fold(p.train_indices,tr_msk);
                        te_idx=find_fold(p.test_indices,te_msk);
                        assertEqual(tr_idx,te_idx);
                        visited_count(tr_idx)=visited_count(tr_idx)+1;
                    end
                end

                assertEqual(visited_count,ones(1,n_folds));

                % also possible with indices
                p2=cosmo_nchoosek_partitioner(ds.sa.chunks,n_test,...
                                            ds.sa.modality,moda_arg);
                assertEqual(p,p2);

            end
        end
    end

function pos=find_fold(folds, msk)
    idxs=find(msk);
    n=numel(folds);

    pos=[];
    for k=1:n
        if isequal(sort(folds{k}(:)),sort(idxs(:)))

            % no duplicates
            assert(isempty(pos));
            pos=k;
        end
    end
    assert(~isempty(pos));

function assert_disjoint(vs,i,j)
    common=intersect(vs(i),vs(j));
    if ~isempty(common)
        assertFalse(true,sprintf('element in common: %d', common(1)));

    end
function test_nchoosek_partitioner_exceptions()
    aet=@(varargin)assertExceptionThrown(@()...
                        cosmo_nchoosek_partitioner(varargin{:}),'');
    for nchunks=2:3
        ds=cosmo_synthetic_dataset('nchunks',nchunks,'ntargets',4);

        aet(ds,-1);
        aet(ds,0);
        aet(ds,1.01);
        aet(ds,[1 1]);
        aet(ds,.99);
        aet(struct,1);
        aet(ds,'foo');
        aet(ds,.5,'foo');
        aet(ds,struct);
        aet(ds,1,1,1);
        aet(ds.sa.chunks,1,'chunks',1);
        aet(ds.sa.chunks,1,'chunks',1,'chunks');

        ds.sa.modality=3; % size mismatch
        aet(ds,1,'modality',1);

    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function increaseHumanGEMVersion(bumpType)
% increaseHumanGEMVersion
%   Increase version for the humanGEM respositories
%
% Input:
%   bumpType      the value has to be one selected among 'major', 'minor' or 'patch'
%
% Usage: increaseHumanGEMVersion(bumpType)
%

%Check if in main:
currentBranch = git('rev-parse --abbrev-ref HEAD');
if ~strcmp(currentBranch,'main')
    error('ERROR: not in main')
end

%Get model path
[ST, I]=dbstack('-completenames');
modelPath=fileparts(fileparts(fileparts(ST(I).file)));

%Bump version number:
versionFile=fullfile(modelPath,'version.txt');
fid = fopen(versionFile,'r');
    oldVersion = fscanf(fid, '%s');
fclose(fid);
oldVersion = str2double(strsplit(oldVersion,'.'));
newVersion = oldVersion;
switch bumpType
    case 'major'
        newVersion(1) = newVersion(1) + 1;
        newVersion(2) = 0;
        newVersion(3) = 0;
    case 'minor'
        newVersion(2) = newVersion(2) + 1;
        newVersion(3) = 0;
    case 'patch'
        newVersion(3) = newVersion(3) + 1;
    otherwise
        error('ERROR: invalid input. Use either "major", "minor" or "patch"')
end
newVersion = num2str(newVersion,'%d.%d.%d');

%Check if history has been updated:
%fid     = fopen('../../history.md','r');
%history = fscanf(fid,'%s');
%fclose(fid);
%if ~contains(history,['human' newVersion ':'])
%    error('ERROR: update history.md first')
%end
%To be included

%Load model:
ymlFile=fullfile(modelPath,'model','Human-GEM.yml');
ihuman = importYaml(ymlFile);

%Include tag and save model:
ihuman.version = newVersion;

%Export model to multiple formats
exportHumanGEM(ihuman,'Human-GEM',modelPath,{'mat', 'txt', 'xml', 'yml', 'xlsx'});

%Update version file:
fid = fopen(versionFile,'wt');
fprintf(fid,newVersion);
fclose(fid);

%Update readme file:
readmeFile=fullfile(modelPath,'README.md');
content = fileread(readmeFile);
content = strrep(content,'{{DATE}}',datestr(now,29));
content = strrep(content,'{{nRXN}}',num2str(length(ihuman.rxns)));
content = strrep(content,'{{nMET}}',num2str(length(ihuman.mets)));
content = strrep(content,'{{nGENE}}',num2str(length(ihuman.genes)));
fid = fopen(readmeFile,'wt');
fwrite(fid,content);
fclose(fid);

end

author: SysBioChalmers

repo: Human-GEM


text:
% eeg_lat2point() - convert latencies in time units relative to the
%                   time locking event of an eeglab() data epoch to
%                   latencies in data points (assuming concatenated epochs).
% Usage:
%       >> [newlat] = eeg_lat2point( lat_array, epoch_array,...
%                                 srate, timelimits, timeunit);
%       >> [newlat] = eeg_lat2point( lat_array, epoch_array,...
%                                 srate, timelimits, '','outrange',1);
% Inputs:
%   lat_array   - latency array in 'timeunit' units (see below)
%   epoch_array - epoch number for each latency
%   srate       - data sampling rate in Hz
%   timelimits  - [min max] epoch timelimits in 'timeunit' units (see below)
%   timeunit    - time unit relative to seconds. Default is 1 = seconds.
%
% Optional inputs:
%   outrange    - [1/0] Replace the points out of the range with the value of
%                 the maximun point in the valid range or raise an error.
%                 Default [1] : Replace point.
%
% Outputs:
%   newlat      - converted latency values in points assuming concatenated
%                 data epochs (see eeglab() event structure)
%   flag        - 1 if any point out of range was replaced.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2 Mai 2002
%
% See also: eeg_point2lat(), eeglab()

% Copyright (C) 2 Mai 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [newlat,flag] = eeg_lat2point( lat_array, epoch_array, srate, timewin, timeunit, varargin);
% -------------------------------------------------------------------------
try
    options = varargin;
    if ~isempty( varargin ),
        for i = 1:2:numel(options)
            g.(options{i}) = options{i+1};
        end
    else g= []; end;
catch
    error('std_checkdatasession() error: calling convention {''key'', value, ... } error');
end;
 try, g.outrange; catch, g.outrange = 1; end; %

 flag = 0;
% -------------------------------------------------------------------------

if nargin <4
    help eeg_lat2point;
    return;
end;
if nargin <5 | isempty(timeunit)
	timeunit = 1;
end;

if length(lat_array) ~= length(epoch_array)
	if length(epoch_array)~= 1
		disp('eeg_lat2point: latency and epochs must have the same length'); return;
	else
		epoch_array = ones(1,length(lat_array))*epoch_array;
	end;
end;
if length(timewin) ~= 2
    disp('eeg_lat2point: timelimits must have length 2'); return;
end;
if iscell(epoch_array)
	epoch_array = [ epoch_array{:} ];
end;
if iscell(lat_array)
	lat_array = [ lat_array{:} ];
end

timewin = timewin*timeunit;
pnts = (timewin(2)-timewin(1))*srate+1;
newlat  = (lat_array*timeunit-timewin(1))*srate+1 + (epoch_array-1)*pnts;

% Detecting points out of range (RMC)
% Note: This is neccesary since the double precision  multiplication could lead to the
% shifting in one sample out of the valid range

if  and(~isempty(newlat),~isempty(epoch_array)) && max(newlat(:)) > max((epoch_array)*pnts)
    if g.outrange == 1
        IndxOut = find(newlat(:) > max((epoch_array)*pnts));
        newlat(IndxOut) = max((epoch_array)*pnts);
        flag = 1;
        warning('eeg_lat2point(): Points out of range detected. Points replaced with maximum value');
    elseif g.outrange == 0
        error('Error in eeg_lat2point(): Points out of range detected');
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function aamax = ivec_amax ( n, a )

%*****************************************************************************80
%
%% IVEC_AMAX returns the largest magnitude in an I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer A(N), the vector to be searched.
%
%    Output, integer AAMAX, the value of the entry of largest magnitude.
%
  aamax = max ( abs ( a(1:n) ) );

  return
end


text:
%% Camera Calibration using ChArUco Boards Demo
% Calibration using a ChArUco board.
%
% To capture a frame for calibration, press 'c',
% If input comes from video, press any key for next frame
% To finish capturing, press 'ESC' key and calibration starts.
%
% Sources:
%
% * <https://docs.opencv.org/3.1.0/da/d13/tutorial_aruco_calibration.html>
% * <https://github.com/opencv/opencv_contrib/blob/3.1.0/modules/aruco/samples/calibrate_camera_charuco.cpp>
%

%% Parameters

% options
vidFile = '';                   % Use video file instead of camera as input
squaresX = 5;                   % Number of squares in X direction
squaresY = 7;                   % Number of squares in Y direction
squareLength = 60;              % Square side length (in pixels)
markerLength = 30;              % Marker side length (in pixels)
dictionaryId = '6x6_250';       % Dictionary id
refindStrategy = false;         % Apply refined strategy
showChessboardCorners = false;  % Show detected chessboard corners after calibration
calibrationFlags = {
    'UseIntrinsicGuess',false, ...
    'FixAspectRatio',false, ...     % Fix aspect ratio (fx/fy)
    'ZeroTangentDist',false, ...    % Assume zero tangential distortion
    'FixPrincipalPoint',false       % Fix the principal point at the center
};
aspectRatio = 1;                    % Fix aspect ratio (fx/fy) to this value

% FixAspectRatio in camera parameters
if calibrationFlags{4}
    camMatrix = eye(3);
    camMatrix(1,1) = aspectRatio;
    calibrationFlags = [calibrationFlags, 'CameraMatrix',camMatrix];
end

% marker detector parameters
detectorParams = struct();
if false
    %detectorParams.nMarkers = 1024;
    detectorParams.adaptiveThreshWinSizeMin = 3;
    detectorParams.adaptiveThreshWinSizeMax = 23;
    detectorParams.adaptiveThreshWinSizeStep = 10;
    detectorParams.adaptiveThreshConstant = 7;
    detectorParams.minMarkerPerimeterRate = 0.03;
    detectorParams.maxMarkerPerimeterRate = 4.0;
    detectorParams.polygonalApproxAccuracyRate = 0.05;
    detectorParams.minCornerDistanceRate = 0.05;
    detectorParams.minDistanceToBorder = 3;
    detectorParams.minMarkerDistanceRate = 0.05;
    detectorParams.cornerRefinementMethod = 'None';
    detectorParams.cornerRefinementWinSize = 5;
    detectorParams.cornerRefinementMaxIterations = 30;
    detectorParams.cornerRefinementMinAccuracy = 0.1;
    detectorParams.markerBorderBits = 1;
    detectorParams.perspectiveRemovePixelPerCell = 8;
    detectorParams.perspectiveRemoveIgnoredMarginPerCell = 0.13;
    detectorParams.maxErroneousBitsInBorderRate = 0.04;
    detectorParams.minOtsuStdDev = 5.0;
    detectorParams.errorCorrectionRate = 0.6;
end

% create charuco board
dictionary = {'Predefined', dictionaryId};
board = {squaresX, squaresY, squareLength, markerLength, dictionary};

%% Input source
if ~isempty(vidFile) && exist(vidFile, 'file') == 2
    vid = cv.VideoCapture(vidFile);
    waitTime = 1;     % 1 sec
else
    vid = cv.VideoCapture(0);
    waitTime = 0.01;  % 10 msec
end
if ~vid.isOpened(), error('failed to initialize VideoCapture'); end

%% Collect

% collect data from each frame
allCorners = {};
allIds = {};
allImgs = {};
imgSize = [];
hImg = []; hFig = [];
while true
    % grab frame
    img = vid.read();
    if isempty(img), break; end

    % detect markers
    [corners, ids, rejected] = cv.detectMarkers(img, dictionary, ...
        'DetectorParameters',detectorParams);

    % refined strategy to detect more markers
    if refindStrategy
        [corners, ids, rejected] = cv.refineDetectedMarkers(img, ...
            ['CharucoBoard',board], corners, ids, rejected);
    end

    % interpolate charuco corners
    if ~isempty(ids)
        [charucoCorners, charucoIds] = cv.interpolateCornersCharuco(...
            corners, ids, img, board);
    end

    % draw results
    out = img;
    if ~isempty(ids)
        out = cv.drawDetectedMarkers(out, corners);  % 'IDs',ids
        if ~isempty(charucoCorners)
            out = cv.drawDetectedCornersCharuco(out, charucoCorners, ...
                'IDs',charucoIds);
        end
    end
    out = cv.putText(out, ['Press "c" to add current frame. ', ...
        '"ESC" to finish and calibrate'], [10 20], ...
        'FontScale',0.5, 'Color',[255 0 0], 'Thickness',2);

    if isempty(hImg)
        hImg = imshow(out);
        hFig = ancestor(hImg, 'figure');
        set(hFig, 'KeyPressFcn',@(o,e) setappdata(o, 'key',e.Key));
        setappdata(hFig, 'key','');
    elseif ishghandle(hImg)
        set(hImg, 'CData',out);
    else
        break;
    end
    drawnow; pause(waitTime);

    % collect frame
    switch getappdata(hFig, 'key')
        case {'space', 'return', 'c'}
            if ~isempty(ids)
                fprintf('Frame captured at %s\n', datestr(now()));
                allCorners{end+1} = corners;
                allIds{end+1} = ids;
                allImgs{end+1} = img;
                imgSize = size(img);
            else
                disp('frame skipped, no corners detected!');
            end
        case {'escape', 'q'}
            disp('Finished collecting frames.');
            break;
        case 'p'
            pause(5);
    end
    setappdata(hFig, 'key','');
end
vid.release();

%% Calibration
disp('Calibrating...')

% calibrate camera using aruco markers
if isempty([allIds{:}]), error('Not enough captures for calibration'); end
[camMatrix, distCoeffs, arucoRepErr] = ...
    cv.calibrateCameraAruco([allCorners{:}], [allIds{:}], ...
        cellfun(@numel, allCorners), ['CharucoBoard',board], ...
        imgSize([2 1]), calibrationFlags{:});

% prepare data for charuco calibration
allCharucoCorners = cell(size(allCorners));
allCharucoIds = cell(size(allCorners));
for i=1:numel(allCorners)
    % interpolate using camera parameters
    [allCharucoCorners{i}, allCharucoIds{i}] = cv.interpolateCornersCharuco(...
        allCorners{i}, allIds{i}, allImgs{i}, board, ...
        'CameraMatrix',camMatrix, 'DistCoeffs',distCoeffs);
    if isempty(allCharucoIds{i})
        warning('interpolateCornersCharuco found no corners');
    end
end

% drop images where it failed to interpolate charuco corners
idx = cellfun(@isempty, allCharucoIds);
allCharucoIds(idx) = [];
allCharucoCorners(idx) = [];
allImgs(idx) = [];

% calibrate camera using charuco
if numel(allCharucoCorners) < 4, error('Not enough corners for calibration'); end
[camMatrix, distCoeffs, repError, rvecs, tvecs] = ...
    cv.calibrateCameraCharuco(allCharucoCorners, allCharucoIds, board, ...
        imgSize([2 1]), calibrationFlags{:});

% calibration results
fprintf('Calibration Time: %s\n', datestr(now()));
fprintf('Image Width: %d, Image Height: %d\n', imgSize(2), imgSize(1));
disp('Flags:'); cellfun(@disp,calibrationFlags);
if calibrationFlags{4}, fprintf('Aspect Ratio: %f\n', aspectRatio); end
disp('Camera Matrix:'); disp(camMatrix)
disp('Distortion Coefficients:'); disp(distCoeffs)
fprintf('Reprojection Error: Charuco = %f, Aruco = %f\n', repError, arucoRepErr);
save camera_parameters.mat -mat camMatrix distCoeffs

% show interpolated charuco corners for debugging
if showChessboardCorners
    axisLength = 0.5 * min(squaresX, squaresY) * squareLength;
    outImgs = {};
    for i=1:numel(allImgs)
        if ~isempty(allIds{i}) && ~isempty(allCharucoCorners)
            img = cv.drawDetectedCornersCharuco(...
                allImgs{i}, allCharucoCorners{i}, 'IDs',allCharucoIds{i});
            img = cv.drawAxis(img, camMatrix, distCoeffs, ...
                rvecs{i}, tvecs{i}, axisLength);
            outImgs{end+1} = img;
        end
    end
    if ~mexopencv.isOctave() && mexopencv.require('images')
        %HACK: IMPLAY not implemented in Octave
        implay(cat(4, outImgs{:}), 1);
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [results, init_target_sz] = tracker(params)

%% Initialisation
% Get sequence information
[seq, im] = get_sequence_info(params.seq);
params = rmfield(params, 'seq');
if isempty(im)
    seq.rect_position = [];
    [~, results] = get_sequence_results(seq);
    return;
end

% Initialise position
pos = seq.init_pos(:)';
target_sz = seq.init_sz(:)';
params.init_sz = target_sz;

% Feature settings and default parameters
features = params.t_features;
params = init_default_params(params);
params.use_mexResize = true;

% Global feature parameters
if isfield(params, 't_global')
    global_fparams = params.t_global;
else
    global_fparams = [];
end
global_fparams.use_gpu = params.use_gpu;
global_fparams.gpu_id = params.gpu_id;
global_fparams.use_mexResize = true;

% Data types
if params.use_gpu
    params.data_type = zeros(1, 'single', 'gpuArray');
else
    params.data_type = zeros(1, 'single');
end
params.data_type_complex = complex(params.data_type);
global_fparams.data_type = params.data_type;

init_target_sz = target_sz;

% Check if color image
if size(im,3) == 3
    if all(all(im(:,:,1) == im(:,:,2)))
        is_color_image = false;
    else
        is_color_image = true;
    end
else
    is_color_image = false;
end

if size(im,3) > 1 && is_color_image == false
    im = im(:,:,1);
end

% Calculate search area and initial scale factor
search_area = prod(init_target_sz * params.search_area_scale);
if search_area > params.max_image_sample_size
    currentScaleFactor = sqrt(search_area / params.max_image_sample_size);
elseif search_area < params.min_image_sample_size
    currentScaleFactor = sqrt(search_area / params.min_image_sample_size);
else
    currentScaleFactor = 1.0;
end

% target size at the initial scale
base_target_sz = target_sz / currentScaleFactor;

% window size, taking padding into account
switch params.search_area_shape
    case 'proportional'
        img_sample_sz = floor(base_target_sz * params.search_area_scale);     % proportional area, same aspect ratio as the target
    case 'square'
        img_sample_sz = repmat(sqrt(prod(base_target_sz * params.search_area_scale)), 1, 2); % square area, ignores the target aspect ratio
    case 'fix_padding'
        img_sample_sz = base_target_sz + sqrt(prod(base_target_sz * params.search_area_scale) + (base_target_sz(1) - base_target_sz(2))/4) - sum(base_target_sz)/2; % const padding
    case 'custom'
        img_sample_sz = [base_target_sz(1)*2 base_target_sz(2)*2];
end

[features, global_fparams, feature_info] = init_features(features, global_fparams, is_color_image, img_sample_sz, 'exact');

% Set feature info
img_support_sz = feature_info.img_support_sz;
feature_sz = unique(feature_info.data_sz, 'rows', 'stable');
feature_cell_sz = unique(feature_info.min_cell_size, 'rows', 'stable');
num_feature_blocks = size(feature_sz, 1);
feature_extract_info = get_feature_extract_info(features);

% Size of the extracted feature maps
feature_sz_cell = mat2cell(feature_sz, ones(1,num_feature_blocks), 2);
filter_sz = feature_sz;
filter_sz_cell = permute(mat2cell(filter_sz, ones(1,num_feature_blocks), 2), [2 3 1]);

% The size of the label function DFT. Equal to the maximum filter size
[output_sz, k1] = max(filter_sz, [], 1);
k1 = k1(1);

% Get the remaining block indices
block_inds = 1:num_feature_blocks;
block_inds(k1) = [];

% Construct the Gaussian label function
yf = cell(numel(num_feature_blocks), 1);
for i = 1:num_feature_blocks
    sz = filter_sz_cell{i};
    output_sigma = sqrt(prod(floor(base_target_sz/feature_cell_sz(i)))) * params.output_sigma_factor;
    rg           = circshift(-floor((sz(1)-1)/2):ceil((sz(1)-1)/2), [0 -floor((sz(1)-1)/2)]);
    cg           = circshift(-floor((sz(2)-1)/2):ceil((sz(2)-1)/2), [0 -floor((sz(2)-1)/2)]);
    [rs, cs]     = ndgrid(rg,cg);
    y            = exp(-0.5 * (((rs.^2 + cs.^2) / output_sigma^2)));
    yf{i}           = fft2(y);
end

% Compute the cosine windows
cos_window = cellfun(@(sz) hann(sz(1))*hann(sz(2))', feature_sz_cell, 'uniformoutput', false);

% Define spatial regularization windows
reg_window = cell(num_feature_blocks, 1);
for i = 1:num_feature_blocks
    reg_scale = floor(1.3*base_target_sz/params.feature_downsample_ratio(i));
    use_sz = filter_sz_cell{i};
    reg_window{i} = ones(use_sz) * 1e5;
    range = zeros(numel(reg_scale), 2);

    % determine the target center and range in the regularization windows
    for j = 1:numel(reg_scale)
        range(j,:) = [0, reg_scale(j) - 1] - floor(reg_scale(j) / 2);
    end
    center = floor((use_sz + 1)/ 2) + mod(use_sz + 1,2);
    range_h = (center(1)+ range(1,1)) : (center(1) + range(1,2));
    range_w = (center(2)+ range(2,1)) : (center(2) + range(2,2));

    reg_window{i}(range_h, range_w) = 1e-3;
end
params.reg_window = reg_window;

% Pre-computes the grid that is used for socre optimization
ky = circshift(-floor((filter_sz_cell{1}(1) - 1)/2) : ceil((filter_sz_cell{1}(1) - 1)/2), [1, -floor((filter_sz_cell{1}(1) - 1)/2)]);
kx = circshift(-floor((filter_sz_cell{1}(2) - 1)/2) : ceil((filter_sz_cell{1}(2) - 1)/2), [1, -floor((filter_sz_cell{1}(2) - 1)/2)])';
newton_iterations = params.newton_iterations;

% Use the translation filter to estimate the scale
nScales = params.number_of_scales;
scale_step = params.scale_step;
scale_exp = (-floor((nScales-1)/2):ceil((nScales-1)/2));
scaleFactors = scale_step .^ scale_exp;

if nScales > 0
    min_scale_factor = scale_step ^ ceil(log(max(5 ./ img_support_sz)) / log(scale_step));
    max_scale_factor = scale_step ^ floor(log(min([size(im,1) size(im,2)] ./ base_target_sz)) / log(scale_step));
end

seq.time = 0;

% Define the learning variables
theta_pre_f = cell(num_feature_blocks, 1);
ctheta_f = cell(num_feature_blocks, 1);
L = cell(num_feature_blocks, 1);
scores_fs_feat = cell(1,1,num_feature_blocks);
while true
    % Read image
    if seq.frame > 0
        [seq, im] = get_sequence_frame(seq);
        if isempty(im)
            break;
        end
        if size(im,3) > 1 && is_color_image == false
            im = im(:,:,1);
        end
    else
        seq.frame = 1;
    end

    tic();

    %% Target localization stage

    if seq.frame > 1
        old_pos = inf(size(pos));
        iter = 1;

        %translation search
        while iter <= params.refinement_iterations && any(old_pos ~= pos)
            % Extract features at multiple resolutions
            sample_pos = round(pos);
            det_sample_pos = sample_pos;
            sample_scale = currentScaleFactor*scaleFactors;
            xt = extract_features(im, sample_pos, sample_scale, features, global_fparams, feature_extract_info);

            % Do windowing of features
            xtw = cellfun(@(feat_map, cos_window) bsxfun(@times, feat_map, sqrt(cos_window)), xt, cos_window, 'uniformoutput', false);

            % Compute the fourier series
            xtf = cellfun(@fft2, xtw, 'uniformoutput', false);

            scores_fs_feat{k1} = gather(sum(bsxfun(@times, conj(ctheta_f{k1}), xtf{k1}), 3));
            scores_fs_sum = scores_fs_feat{k1};
            for k = block_inds
                scores_fs_feat{k} = gather(sum(bsxfun(@times, conj(ctheta_f{k}), xtf{k}), 3));
                scores_fs_feat{k} = resizeDFT2(scores_fs_feat{k}, output_sz);
                scores_fs_sum = scores_fs_sum +  scores_fs_feat{k};
            end
            scores_fs = permute(gather(scores_fs_sum), [1 2 4 3]);

            responsef_padded = resizeDFT2(scores_fs, output_sz);
            response = ifft2(responsef_padded, 'symmetric');
            [disp_row, disp_col, sind] = resp_newton(response, responsef_padded, newton_iterations, ky, kx, output_sz);

            % Compute the translation vector in pixel-coordinates and round
            % to the closest integer pixel.
            translation_vec = [disp_row, disp_col] .* (img_support_sz./output_sz) * currentScaleFactor * scaleFactors(sind);
            scale_change_factor = scaleFactors(sind);

            % update position
            old_pos = pos;
            if sum(isnan(translation_vec))
                pos = sample_pos;
            else
                pos = sample_pos + translation_vec;
            end

            if params.clamp_position
                pos = max([1 1], min([size(im,1) size(im,2)], pos));
            end

            % Update the scale
            currentScaleFactor = currentScaleFactor * scale_change_factor;

            % Adjust to make sure we are not to large or to small
            if currentScaleFactor < min_scale_factor
                currentScaleFactor = min_scale_factor;
            elseif currentScaleFactor > max_scale_factor
                currentScaleFactor = max_scale_factor;
            end

            iter = iter + 1;
        end
    end

    %% Model learning and update stage

    % extract image region for training sample
    sample_pos = round(pos);
    xl = extract_features(im, sample_pos, currentScaleFactor, features, global_fparams, feature_extract_info);

    % do windowing of features
    xlw = cellfun(@(feat_map, cos_window) bsxfun(@times, feat_map, cos_window), xl, cos_window, 'uniformoutput', false);

    % compute the fourier series
    xlf = cellfun(@fft2, xlw, 'uniformoutput', false);

    [ctheta_f, theta_pre_f, L] = train_filter(ctheta_f, xlf, yf, theta_pre_f, params, output_sz, seq);

    % Update target size
    target_sz = base_target_sz * currentScaleFactor;

    % save position and time
    tracking_result.center_pos = double(pos);
    tracking_result.target_size = double(target_sz);
    seq = report_tracking_result(seq, tracking_result);

    seq.time = seq.time + toc();

    %% Visualization

    if params.visualization
        rect_position_vis = [pos([2,1]) - (target_sz([2,1]) - 1)/2, target_sz([2,1])];
        im_to_show = double(im)/255;
        if size(im_to_show,3) == 1
            im_to_show = repmat(im_to_show, [1 1 3]);
        end
        if seq.frame == 1
            fig_handle = figure('Name', 'Tracking');
            imagesc(im_to_show);
            hold on;
            rectangle('Position',rect_position_vis, 'EdgeColor','g', 'LineWidth',2);
            text(10, 10, int2str(seq.frame), 'color', [0 1 1]);
            hold off;
            axis off;axis image;set(gca, 'Units', 'normalized', 'Position', [0 0 1 1])

        else
            % Do visualization of the sampled confidence scores overlayed
%             resp_sz = round(img_support_sz*currentScaleFactor*scaleFactors(sind));
%             xs = floor(det_sample_pos(2)) + (1:resp_sz(2)) - floor(resp_sz(2)/2);
%             ys = floor(det_sample_pos(1)) + (1:resp_sz(1)) - floor(resp_sz(1)/2);
%             spatial_feature_display = circshift(imresize(L{1},[numel(xs),numel(ys)]);
            figure(fig_handle);
            imagesc(im_to_show);
%             hold on;
%             resp_handle = imagesc(xs, ys, spatial_feature_display); colormap hsv;
%             alpha(resp_handle, 0.5);
            rectangle('Position',rect_position_vis, 'EdgeColor','g', 'LineWidth',2);
            text(10, 10, int2str(seq.frame), 'color', [0 1 1]);
            hold off;

        end

        drawnow

    end

end

[seq, results] = get_sequence_results(seq);

author: XU-TIANYANG

repo: LADCF

save_path: github-repos/MATLAB/XU-TIANYANG-LADCF


text:
classdef Quad < handle
    %QUAD Summary of this class goes here
    %   Detailed explanation goes here

    properties
    V00,V01,V10,V11;
    end

    methods
        function obj =Quad(V00,V01,V10,V11)
            obj.V00 = V00;
            obj.V01 = V01;
            obj.V10 = V10;
            obj.V11 = V11;
        end
        function a = isPointIn(obj,pt)
           in1 = isPointInTriangular(pt,obj.V00,obj.V01,obj.V11);
           in2 = isPointInTriangular(pt,obj.V00,obj.V10,obj.V11);
           if in1==1 || in2==1
               a = 1;
           else
               a = 0;
           end
        end

        function a = isPointsIn(obj,pts)
            in1 = isPointsInTriangular(pts(:,1),pts(:,2),obj.V00,obj.V01,obj.V11);
            in2 = isPointsInTriangular(pts(:,1),pts(:,2),obj.V00,obj.V10,obj.V11);

            a = double(in1|in2);

        end

        function [coefficients,bool] = getBilinearCoordinates(obj,pt)

            a_x = obj.V00.x - obj.V01.x - obj.V10.x + obj.V11.x;
            b_x = -obj.V00.x + obj.V01.x;
            c_x = -obj.V00.x + obj.V10.x;
            d_x = obj.V00.x - pt.x;

            a_y = obj.V00.y - obj.V01.y - obj.V10.y + obj.V11.y;
            b_y = -obj.V00.y + obj.V01.y;
            c_y = -obj.V00.y + obj.V10.y;
            d_y = obj.V00.y - pt.y;

            bigA = -a_y*b_x + b_y*a_x;
            bigB = -a_y*d_x - c_y*b_x + d_y*a_x +b_y*c_x;
            bigC = -c_y*d_x + d_y*c_x;

            tmp1 = -1;
            tmp2 = -1;
            tmp3 = -1;
            tmp4 = -1;

            if bigB*bigB - 4*bigA*bigC >= 0.0
                if abs(bigA) >= 0.000001
                    tmp1 = ( -bigB + sqrt(bigB*bigB - 4*bigA*bigC) ) / ( 2*bigA );
                    tmp2 = ( -bigB - sqrt(bigB*bigB - 4*bigA*bigC) ) / ( 2*bigA );
                else
                    tmp1 = -bigC/bigB;
                end

                if ( tmp1 >= -0.999999 && tmp1 <= 1.000001)

                    tmp3 = -(b_y*tmp1 + d_y) / (a_y*tmp1 + c_y);
                    tmp4 = -(b_x*tmp1 + d_x) / (a_x*tmp1 + c_x);
                    if tmp3 >= -0.999999 && tmp3 <= 1.000001
                        k1 = tmp1;
                        k2 = tmp3;
                    else if tmp4 >= -0.999999 && tmp4 <= 1.000001
                            k1 = tmp1;
                            k2 = tmp4;
                        end
                    end
                end
                if ( tmp2 >= -0.999999 && tmp2 <= 1.000001)

                    if tmp3 >= -0.999999 && tmp3 <= 1.000001
                        k1 = tmp2;
                        k2 = tmp3;
                    else if tmp4 >= -0.999999 && tmp4 <= 1.000001
                            k1 = tmp2;
                            k2 = tmp4;
                        end
                    end
                end
            end

            if k1>=-0.999999 && k1<=1.000001 && k2>=-0.999999 && k2<=1.000001

                coe1 = (1.0-k1)*(1.0-k2);
                coe2 = k1*(1.0-k2);
                coe3 = (1.0-k1)*k2;
                coe4 = k1*k2;

                coefficients(1) = coe1;
                coefficients(2) = coe2;
                coefficients(3) = coe3;
                coefficients(4) = coe4;

                bool = 1;
            else
                bool = 0;
            end

        end

        function minx = getMinX(obj)
            minx = min(obj.V00.x,obj.V01.x);
            minx = min(minx, obj.V10.x);
            minx = min(minx,obj.V11.x);
        end
        function maxx = getMaxX(obj)
            maxx = max(obj.V00.x,obj.V01.x);
            maxx = max(maxx,obj.V10.x);
            maxx = max(maxx,obj.V11.x);
        end
        function miny = getMinY(obj)
            miny = min(obj.V00.y,obj.V01.y);
            miny = min(miny,obj.V10.y);
            miny = min(miny,obj.V11.y);
        end
        function maxy = getMaxY(obj)
           maxy = max(obj.V00.y,obj.V01.y);
           maxy = max(maxy,obj.V10.y);
           maxy = max(maxy,obj.V11.y);
        end

    end
end

function a = isPointInTriangular(pt,V0,V1,V2)
    lambda1 = ((V1.y-V2.y)*(pt.x-V2.x) + (V2.x-V1.x)*(pt.y-V2.y)) / ((V1.y-V2.y)*(V0.x-V2.x) + (V2.x-V1.x)*(V0.y-V2.y));
	lambda2 = ((V2.y-V0.y)*(pt.x-V2.x) + (V0.x-V2.x)*(pt.y-V2.y)) / ((V2.y-V0.y)*(V1.x-V2.x) + (V0.x-V2.x)*(V1.y-V2.y));
    lambda3 = 1-lambda1-lambda2;
    if lambda1 >= 0.0 && lambda1 <= 1.0 && lambda2 >= 0.0 && lambda2 <= 1.0 && lambda3 >= 0.0 && lambda3 <= 1.0
        a = 1; %true
    else
        a = 0; %false
    end
end

function a = isPointsInTriangular(ptx,pty,V0,V1,V2)
    y12 = V1.y - V2.y;
    x21 = V2.x - V1.x;
    x02 = V0.x - V2.x;
    x21 = V2.x - V1.x;
    y02 = V0.y - V2.y;
    y20 = V2.y - V0.y;
    x12 = V1.x - V2.x;
    s1 = y12*x02 + x21*y02;
    s2 = y20*x12 + x02*y12;

    lambda1 = (y12.*(ptx-V2.x) + x21.*(pty-V2.y)) ./ s1;
    lambda2 = (y20.*(ptx-V2.x) + x02.*(pty-V2.y)) ./ s2;
    lambda3 = 1-lambda1-lambda2;

    condition1 = lambda1>=0.0;
    condition2 = lambda1<=1.0;
    condition3 = lambda2>=0.0;
    condition4 = lambda2<=1.0;
    condition5 = lambda3>=0.0;
    condition6 = lambda3<=1.0;

    a = condition1 & condition2 & condition3 & condition4 & condition5 & condition6;

end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
filename='ImprovedBridge_hexahedra';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance','perimeter'};
weights = [1 0.1];
constraint = {'volume'};
optimizer = 'HAMILTON-JACOBI'; incrementFactor = 1;
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';
showBC = true;

HJiter0 = 1;
e2 = 1;
N_holes = [12 5 5];
R_holes = 0.2;
phase_holes = [0 0 0];
nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 0.7;
optimality_initial = 5e-2;
constr_initial = 5e-2;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [Xnorm] = wL1norm(X,lambda)
% ||X||_wL1 = sum_i \lambda_ij|X_ij|
Xnorm = sum(sum(lambda.*abs(X),2));

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function [motionC, measC] = factorCoord(Lmk,Frm,Fac)

nframes  = sum([Frm.used]);
nlmks    = sum([Lmk.used]);

% Motion
motionC = zeros(nframes,3);

% Measurements
measC = zeros(nframes+nlmks,3);

for fac = find([Fac.used])

    frames = Fac(fac).frames;
    lmk = Fac(fac).lmk;

    switch Fac(fac).type
        case'motion'
            motionC(frames(1),1:3) = Frm(frames(1)).state.x(1:3)';
            motionC(frames(2),1:3) = Frm(frames(2)).state.x(1:3)';
        case 'absolute'
        case 'measurement'
            measC(frames,1:3) = Frm(frames).state.x(1:3)';
            measC(lmk+nframes,1:3) = Lmk(lmk).state.x(1:3)';
    end

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function spiral_test03 ( )

%*****************************************************************************80
%
%% SPIRAL_TEST03 generates a field on a regular grid and plots it.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPIRAL_TEST03:\n' );
  fprintf ( 1, '  Generate a spiral velocity field on a regular grid.\n' );
  fprintf ( 1, '  Store in GNUPLOT data and command files.\n' );

  x_lo = 0.0;
  x_hi = 1.0;
  x_num = 21;

  y_lo = 0.0;
  y_hi = 1.0;
  y_num = 21;

  [ x, y ] = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi );

  n = x_num * y_num;
  c = 1.0;

  [ u, v ] = uv_spiral ( n, x, y, c );

  header = 'spiral';
  s = 0.05;
  spiral_gnuplot ( header, n, x, y, u, v, s );

  return
end

repo: jburkardt-m


text:
function fake_output = mmx(varargin)
%MMX - Multithreaded matrix operations on N-D matrices
%    MMX treats an N-D matrix of double precision values as a set of pages
%    of 2D matrices, and performs various matrix operations on those pages.
%    MMX uses multithreading over the higher dimensions to achieve good
%    performance. Full singleton expansion is available for most operations.
%
%    C = MMX('mult', A, B) is equivalent to the matlab loop
%    for i=1:N,
%        C(:,:,i) = A(:,:,i) * B(:,:,i);
%    end
%    Singleton expansion is enabled on all dimensions so for example if
%    A = randn(5,4,3,10,1);
%    B = randn(4,6,3,1 ,6);
%    C = zeros(5,6,3,10,6);
%    then C = mmx('mult',A,B) equivalent to
%    for i = 1:3
%       for j = 1:10
%          for k = 1:6
%             C(:,:,i,j,k) = A(:,:,i,j,1) * B(:,:,i,1,k);
%          end
%       end
%    end
%
%    C = MMX('mult', A, B, mod) and where mod is a modifier string, will
%    transpose one or both of A and B. Possible values for mod are
%    'tn', 'nt' and  'tt' where 't' stands for 'transposed' and 'n' for
%    'not-transposed'. For example
%    >> size(mmx('mult',randn(4,2),randn(4,2),'tn'))
%    ans =   2     2
%
%    C = MMX('square', A, [])     will perform C = A*A'
%    C = MMX('square', A, [],'t') will perform C = A'*A
%
%    C = MMX('square', A, B)       will perform C = 0.5*(A*B'+B*A')
%    C = MMX('square', A, B, 't')  will perform C = 0.5*(A'*B+B'*A)
%
%    C = MMX('chol',   A, []) will perform C = chol(A)
%
%    C = MMX('backslash', A, B) will perform C = A\B
%    Unlike other commands, 'backslash' does not support singleton
%    expansion. If A is square, mmx will use LU factorization, otherwise it
%    will use QR factorization. In the underdetermined case, (i.e. when
%    size(A,1) < size(A,2)), mmx will give the least-norm solution which
%    is equivalent to C = pinv(A)*B, unlike matlab's mldivide.
%
%    C = MMX('backslash', A, B, 'U') or MMX('backslash', A, B, 'L') will
%    perform C = A\B assuming that A is upper or lower triangular,
%    respectively.
%
%    C = MMX('backslash', A, B, 'P') will perform C = A\B assuming that A
%    is symmetric-positive-definite.
%
%    MMX(n) does thread control: mmx will automatically start a number of
%    threads equal to the number of available processors, however the
%    number can be set manually to n using the command mmx(n). mmx(0) will
%    clear the threads from memory.
%
%    IMPORTANT NOTE: The functions which assume special types of square
%    matrices as input ('chol' and 'backslash' for 'U','L' or 'P'
%    modifiers) do not check that the inputs are indeed what you say they
%    are, and produce no error if they are not. Caveat computator.
%
%    COMPILATION: To compile run 'build_mmx'. Type 'help build_mmx' to read
%    about compilation issues and options

error(sprintf('MEX file not found.\nTry ''build_mmx''.\nType ''help mmx'' for details.'));

author: NicolasBoumal

repo: manopt


text:
%% Birefrigence
%
% Birefringence is the optical property of a material having a refractive
% index that depends on the polarization and propagation direction of
% light. It is one of the oldest methods to determine orientations of
% crystals in thin sections of rocks.

%% Import Olivine Data
% In order to illustarte the effect of birefringence lets consider a
% olivine data set.

mtexdata olivine

% reconstruct grains
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'));
ebsd(grains(grains.grainSize < 5)) = [];
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'));

% some data denoising
grains = smooth(grains,5);

F = splineFilter;
ebsd = smooth(ebsd,F,'fill',grains);

%%

% plot the olivine phase
plot(ebsd('olivine'),ebsd('olivine').orientations);
hold on
plot(grains.boundary,'lineWidth',2)
hold off

gg = grains(grains.grainSize > 100);
gg = gg('o')
cS = crystalShape.olivine;
hold on
plot(gg,0.8*cS,'FaceColor','none')
hold off

%% The refractive index tensor
%
% The refractive index of a material describes the dependence of the speed
% of light with respect to the propagation direction and the polarization
% direction. In a linear world this relation ship is modeled by a symmetric
% rank 2 tensor - the so called refractive index tensor, which is usually
% given by it principle values: n_alpha, n_beta and n_gamma. In
% orthorhombic minerals such as olivine the principal values are parallel
% to the crystallographic axes. Care has to be applied when associating the
% principle values with the correct axes.

%%
% For Forsterite the priniple refractive values are
n_alpha = 1.635; n_beta = 1.651; n_gamma = 1.670;

%%
% with the largest refractive index n_gamma beeing alligned with the
% a-axis, the intermediate index n_beta with the c-axis and the smallest
% refractive index n_alpha with the b-axis. Hence, the refractive index
% tensor for Forsterite takes the form

cs = ebsd('olivine').CS;
rI_Fo = refractiveIndexTensor(diag([ n_gamma  n_alpha  n_beta]),cs)

%%
% For Fayalite the priniple refractive values

n_alpha = 1.82; n_beta = 1.869; n_gamma = 1.879;

%%
% are aligned to the crystallograhic axes in an analogous way. Which leads
% to the refractive index tensor

rI_Fa = refractiveIndexTensor(diag([ n_gamma  n_alpha  n_beta]),cs)

%%
% The refractive index of composite materials like Olivine can now be
% modelled as the weighted sum of the of the refractive index tensors of
% Forsterite and Fayalite. Lets assume that the relative Forsterite content
% (atomic percentage) is sgiven my

XFo = 0.86; % 86 percent Forsterite

%%
% Then is refractive index tensor becomes

rI = XFo*rI_Fo + (1-XFo) * rI_Fa

%% Birefringence
% The birefringence describes the difference |n| in diffraction index
% between the fastest polarization direction |pMax| and the slowest
% polarization direction |pMin| for a given propagation direction |vprop|.

% lets define a propagation direction
vprop = Miller(1,1,1,cs);

% and compute the birefringence
[dn,pMin,pMax] = rI.birefringence(vprop)

%%
% If the polarization direction is ommited the results are spherical
% functions which can be easily visualized.

% compute the birefringence as a spherical function
[dn,pMin,pMax] = rI.birefringence

% plot it
plot3d(dn,'complete')
mtexColorbar

% and on top of it the polarization directions
hold on
quiver3(pMin,'color','white')
quiver3(pMax)
hold off

%% The Optical Axis
% The optial axes are all directions where the birefringence is zero

% compute the optical axes
vOptical = rI.opticalAxis

% and check the birefringence is zero
rI.birefringence(rI.opticalAxis)

% annotate them to the birefringence plot
hold on
arrow3d(vOptical,'antipodal','facecolor','red')
hold off

%% Spectral Transmission
% If white light with a certain polarization is transmited though a crystal
% with isotropic refrative index the light changes wavelength and hence
% appears collored. The resulting color depending on the propagation
% direction, the polarization direction and the thickness can be computed
% by

vprop = Miller(1,1,1,cs);
thickness = 30000;
p =  Miller(-1,1,0,cs);
rgb = rI.spectralTransmission(vprop,thickness,'polarizationDirection',p)

%%
% Effectively, the rgb value depend only on the angle tau between the
% polariztzion direction and the slowest polarization direction |pMin|.
% Instead of the polarization direction this angle may be specified
% directly

rgb = rI.spectralTransmission(vprop,thickness,'tau',30*degree)

%%
% If the angle tau is fixed and the propagation direction is ommited as
% input MTEX returns the rgb values as a spherical function. Lets plot
% these functions for different values of tau.

newMtexFigure('layout',[1,3]);

mtexTitle('$\tau = 15^{\circ}$')
plot(rI.spectralTransmission(thickness,'tau',15*degree),'rgb')

nextAxis
mtexTitle('$\tau = 30^{\circ}$')
plot(rI.spectralTransmission(thickness,'tau',30*degree),'rgb')

nextAxis
mtexTitle('$\tau = 45^{\circ}$')
plot(rI.spectralTransmission(thickness,'tau',45*degree),'rgb')

drawNow(gcm,'figSize','normal')

%%
% Usually, the polarization direction is chosen at angle phi = 90 degree of
% the analyzer. The following plots demonstrate how to change this angle

newMtexFigure('layout',[1,3]);

mtexTitle('$\tau = 15^{\circ}$')
plot(rI.spectralTransmission(thickness,'tau',45*degree,'phi',30*degree),'rgb')

nextAxis
mtexTitle('$\tau = 30^{\circ}$')
plot(rI.spectralTransmission(thickness,'tau',45*degree,'phi',60*degree),'rgb')

nextAxis
mtexTitle('$\tau = 45^{\circ}$')
plot(rI.spectralTransmission(thickness,'tau',45*degree,'phi',90*degree),'rgb')

drawNow(gcm,'figSize','normal')

%% Spectral Transmission at Thin Sections
% All the above computations have been performed in crystal coordinates.
% However, in practical applications the direction of the polarizer as well
% as the propagation direction are given in terms of specimen coordinates.

% the propagation direction
vprop = vector3d.Z;

% the direction of the polarizer
polarizer = vector3d.X;

% the thickness of the thin section
thickness = 22800;

%%
% As usal we have two options: Either we transform the refractive index
% tensor into specimen coordinates or we transform the polarization
% direction and the propagation directions into crystal coordinates.
% Lets start with the first option:

% extract the olivine orientations
ori = ebsd('olivine').orientations;

% transform the tensor into a list of tensors with respect to specimen
% coordinates
rISpecimen = ori * rI;

% compute RGB values
rgb = rISpecimen.spectralTransmission(vprop,thickness,'polarizationDirection',polarizer);

% colorize the EBSD maps according to spectral transmission
plot(ebsd('olivine'),rgb)

%%
% and compare it with option two:

% transfom the propation direction and the polarizer direction into a list
% of directions with respect to crystal coordinates
vprop_crystal = ori \ vprop;
polarizer_crystal = ori \ polarizer;

% compute RGB values
rgb = rI.spectralTransmission(vprop_crystal,thickness,'polarizationDirection',polarizer_crystal);

% colorize the EBSD maps according to spectral transmission
plot(ebsd('olivine'),rgb)

%% Spectral Transmission as a color key
% The above computations can be automized by defining a spectral
% transmission color key.

% define the colorKey
colorKey  = spectralTransmissionColorKey(rI,thickness);

% the following are the defaults and can be ommited
colorKey.propagationDirection = vector3d.Z;
colorKey.polarizer = vector3d.X;
colorKey.phi = 90 * degree;

% compute the spectral transmission color of the olivine orientations
rgb = colorKey.orientation2color(ori);

plot(ebsd('olivine'), rgb)

%%
% As usual we me visualize the color key as a colorization of the
% orientation space, e.g., by plotting it in sigma sections:

plot(colorKey,'sigma')

%% Circular Polarizer
% In order to simulate we a circular polarizer we simply set the polarizer
% direction to empty, i.e.

colorKey.polarizer = [];

% compute the spectral transmission color of the olivine orientations
rgb = colorKey.orientation2color(ori);

plot(ebsd('olivine'), rgb)

%% Illustrating the effect of rotating polarizer and analyser simultanously

colorKey.polarizer = vector3d.X;
figure
plotHandle = plot(ebsd('olivine'),colorKey.orientation2color(ori),'micronbar','off');
hold on
plot(grains.boundary,'lineWidth',2)
hold off
textHandle = text(750,50,[num2str(0,'%10.1f') '\circ'],'fontSize',15,...
  'color','w','backGroundColor', 'k');

% define the step size in degree
stepSize = 2.5;

for omega = 0:stepSize:90-stepSize

  % update polarsation direction
  colorKey.polarizer = rotate(vector3d.X, omega * degree);

  % update rgb values
  plotHandle.FaceVertexCData = colorKey.orientation2color(ori);

  % update text
  textHandle.String = [num2str(omega,'%10.1f') '\circ'];

  drawnow

end

author: mtex-toolbox

repo: mtex

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


text:
function value = p10_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P10_F evaluates the integrand for problem 10.
%
%  Dimension:
%
%    DIM_NUM arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integrand:
%
%    sum ( abs ( x(1:dim_num) - 0.5 ) )
%
%  Exact Integral:
%
%    DIM_NUM / 4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%    Thomas Patterson,
%    [Integral #4],
%    On the Construction of a Practical Ermakov-Zolotukhin
%    Multiple Integrator,
%    in Numerical Integration: Recent Developments, Software
%    and Applications,
%    edited by Patrick Keast, Graeme Fairweather,
%    D. Reidel, 1987, pages 269-290,
%    LC: QA299.3.N38.
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  value(1:point_num) = 0.0;

  for point = 1 : point_num
    value(point) = sum ( abs ( x(1:dim_num,point) - 0.5 ) );
  end

  p10_i4 ( 'I', '#', point_num );

  return
end

author: johannesgerer


text:
function out = fatrix2_imask_array(st)
if isempty(st.imask)
	out = true([st.idim 1]);
else
	out = st.imask;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function g = sinh(f)
%SINH   Hyperbolic sine of a BALLFUN.
%   SINH(F) computes the hyperbolic sine of the BALLFUN F.
%
% See also SIN, COS.

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

g = compose( f, @sinh );
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function quad_error = monomial_quadrature ( dim_num, expon, point_num, ...
  weight, x, rule )

%*****************************************************************************80
%
%% MONOMIAL_QUADRATURE applies a quadrature rule to a monomial.
%
%  Discussion:
%
%    This routine assumes that the integral being approximated is that of
%    a multidimensional monomial, integrated over the [-1,+1] hypercube,
%    with a Legendre weight (that is, w(x) = 1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer EXPON(DIM_NUM), the exponents.
%
%    Input, integer POINT_NUM, the number of points in the rule.
%
%    Input, real WEIGHT(POINT_NUM), the quadrature weights.
%
%    Input, real X(DIM_NUM,POINT_NUM), the quadrature points.
%
%    Input, integer RULE, the index of the rule.
%    1, "CC", Clenshaw Curtis Closed Fully Nested rule.
%    2, "F1", Fejer 1 Open Fully Nested rule.
%    3, "F2", Fejer 2 Open Fully Nested rule.
%    4, "GP", Gauss Patterson Open Fully Nested rule.
%    5, "GL", Gauss Legendre Open Weakly Nested rule.
%    6, "GH", Gauss Hermite Open Weakly Nested rule.
%    7, "LG", Gauss Laguerre Open Non Nested rule.
%
%    Output, real QUAD_ERROR, the quadrature error.
%

%
%  Get the exact value of the integral of the monomial.
%
  if ( 1 <= rule & rule <= 5 )
    exact = monomial_integral_legendre ( dim_num, expon );
  elseif ( rule == 6 )
    exact = monomial_integral_hermite ( dim_num, expon );
  elseif ( rule == 7 )
    exact = monomial_integral_laguerre ( dim_num, expon );
  end
%
%  Evaluate the monomial at the quadrature points.
%
  value = monomial_value ( dim_num, point_num, x, expon );
%
%  Compute the quadrature sum.
%
  quad = weight * value';
%
%  Absolute error if EXACT = 0, relative error otherwise:
%
  if ( exact == 0.0 )
    quad_error = abs ( quad - exact );
  else
    quad_error = abs ( quad - exact ) / abs ( exact );
  end

  return
end


text:
function p = prior_sqinvgamma(varargin)
%PRIOR_SQINVGAMMA  Gamma prior structure for square inverse of the parameter
%
%  Description
%    P = PRIOR_SQINVGAMMA('PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    creates Gamma prior structure for square inverse of the
%    parameter in which the named parameters have the specified
%    values. Any unspecified parameters are set to default values.
%
%    P = PRIOR_SQINVGAMMA(P, 'PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    modify a prior structure with the named parameters altered
%    with the specified values.
%
%    Parametrisation is done by Bayesian Data Analysis,
%    second edition, Gelman et.al. 2004.
%
%    Parameters for Gamma prior [default]
%      sh       - shape [4]
%      is       - inverse scale [1]
%      sh_prior - prior for sh [prior_fixed]
%      is_prior - prior for is [prior_fixed]
%
%  See also
%    PRIOR_*

% Copyright (c) 2000-2001,2010,2012 Aki Vehtari
% Copyright (c) 2010 Jaakko Riihim�ki

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'PRIOR_SQINVGAMMA';
  ip.addOptional('p', [], @isstruct);
  ip.addParamValue('sh',4, @(x) isscalar(x) && x>0);
  ip.addParamValue('sh_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('is',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('is_prior',[], @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'SqInv-Gamma';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'SqInv-Gamma')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('sh',ip.UsingDefaults)
    p.sh = ip.Results.sh;
  end
  if init || ~ismember('is',ip.UsingDefaults)
    p.is = ip.Results.is;
  end
  % Initialize prior structure
  if init
    p.p=[];
  end
  if init || ~ismember('sh_prior',ip.UsingDefaults)
    p.p.sh=ip.Results.sh_prior;
  end
  if init || ~ismember('is_prior',ip.UsingDefaults)
    p.p.is=ip.Results.is_prior;
  end

  if init
    % set functions
    p.fh.pak = @prior_sqinvgamma_pak;
    p.fh.unpak = @prior_sqinvgamma_unpak;
    p.fh.lp = @prior_sqinvgamma_lp;
    p.fh.lpg = @prior_sqinvgamma_lpg;
    p.fh.recappend = @prior_sqinvgamma_recappend;
  end

end

function [w, s] = prior_sqinvgamma_pak(p)

  w=[];
  s={};
  if ~isempty(p.p.sh)
    w = log(p.sh);
    s=[s; 'log(SqInv-Gamma.sh)'];
  end
  if ~isempty(p.p.is)
    w = [w log(p.is)];
    s=[s; 'log(SqInv-Gamma.is)'];
  end
end

function [p, w] = prior_sqinvgamma_unpak(p, w)

  if ~isempty(p.p.sh)
    i1=1;
    p.sh = exp(w(i1));
    w = w(i1+1:end);
  end
  if ~isempty(p.p.is)
    i1=1;
    p.is = exp(w(i1));
    w = w(i1+1:end);
  end
end

function lp = prior_sqinvgamma_lp(x, p)

  lJ = -log(x)*3 + log(2);  % log(-2/x^3) log(|J|) of transformation
  xt = x.^-2;               % transformation
  lp = sum(-p.is.*xt + (p.sh-1).*log(xt) +p.sh.*log(p.is)  -gammaln(p.sh) +lJ);

  if ~isempty(p.p.sh)
    lp = lp + p.p.sh.fh.lp(p.sh, p.p.sh) + log(p.sh);
  end
  if ~isempty(p.p.is)
    lp = lp + p.p.is.fh.lp(p.is, p.p.is) + log(p.is);
  end
end

function lpg = prior_sqinvgamma_lpg(x, p)

  lJg = -3./x;              % gradient of log(|J|) of transformation
  xt  = x.^-2;              % transformation
  xtg = -2/x.^3;            % derivative of transformation
  lpg = xtg.*((p.sh-1)./xt - p.is) + lJg;

  if ~isempty(p.p.sh)
    lpgsh = (sum(-digamma1(p.sh) + log(p.is) + log(x)) + p.p.sh.fh.lpg(p.sh, p.p.sh)).*p.sh + 1;
    lpg = [lpg lpgsh];
  end
  if ~isempty(p.p.is)
    lpgis = (sum(p.sh./p.is+x) + p.p.is.fh.lpg(p.is, p.p.is)).*p.is + 1;
    lpg = [lpg lpgis];
  end

end

function rec = prior_sqinvgamma_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
  if ~isempty(p.p.sh)
    rec.sh(ri,:) = p.sh;
  end
  if ~isempty(p.p.is)
    rec.is(ri,:) = p.is;
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function s=compressShape(w, forceCompression)
%compressShape Compress a gradient or pulse shape.
%   s=compressShape(w) Compress the waveform using a run-length compression
%   scheme on the derivative. This strategy encodes constant and linear
%   waveforms with very few samples. A structure is returned with the
%   fields:
%     num_samples - the number of samples in the uncompressed waveform
%     data - containing the compressed waveform
%
%   See also decompressShape

if nargin<2
    forceCompression=false;
end

if any(~isfinite(w))
    error('compressShape() received infinite samples');
end

if ~forceCompression && length(w) <= 4 % avoid compressing very short shapes
    s.num_samples=length(w);
    s.data = w(:)';
    return;
end

% %MZ: old code with implicit quantization
% data = [w(1); diff(w(:))];
% maskChanges = [true; abs(diff(data))>1e-8];   % TRUE if values change
% vals = data(maskChanges);                     % Elements without repetitions

% MZ: explicit quantization with error correction
quant_fac=1e-7; % single precision floating point has ~7.25 decimal places
ws=w./quant_fac;
datq=round([ws(1); diff(ws(:))]);
qerr=ws(:)-cumsum(datq);
qcor=[0; diff(round(qerr))];
datd=datq+qcor;
maskChanges=[true; diff(datd)~=0];
vals=datd(maskChanges).*quant_fac;            % Elements without repetitions

k = find([maskChanges', true]);               % Indices of changes
n = diff(k)';                                 % Number of repetitions

% Encode in Pulseq format
nExtra=n-2;
vals2=vals;
vals2(nExtra<0)=nan;
nExtra(nExtra<0)=nan;
v=[vals vals2 nExtra]';
v=v(isfinite(v));
v(abs(v)<1e-10)=0;
s.num_samples = length(w);
% decide whether compression makes sense, otherwise store the original
if forceCompression || s.num_samples > length(v)
    s.data = v';
else
    s.data = w(:)';
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function Offspring = BinaryMutation(Offspring)
% Unbalanced binary mutation

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

    [N,D] = size(Offspring);
    MOne  = mean(Offspring,2);
    r     = min(min(1/D,2*MOne),2*(1-MOne));
    rate1 = repmat(r./2./MOne,1,D);
    rate0 = repmat(r./2./(1-MOne),1,D);
    rate  = zeros(N,D);
    rate(Offspring)  = rate1(Offspring);
    rate(~Offspring) = rate0(~Offspring);
    exchange            = rand(N,D) < rate;
    Offspring(exchange) = ~Offspring(exchange);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function x=logistic(n,level,a,x0)
%Syntax: x=logistic(n,level,a,x0)
%________________________________
%
% Simulation of the Quadratic map.
%    x'=ax(1-x)
%
% x is the simulated time series.
% n is the number of the simulated points.
% level is the noise standard deviation divided by the standard deviation of the
%   noise-free time series. We assume Gaussian noise with zero mean.
% a is the parameter.
% x0 is the initial value for x.
%
%
% Reference:
%
% May R M (1976): Simple mathematical modelswith very complicated dynamics.
% Nature 261: 459-467
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 17 Nov 2001

if nargin<1 | isempty(n)==1
    n=500;
else
    % n must be scalar
    if sum(size(n))>2
        error('n must be scalar.');
    end
    % n must be positive
    if n<0
        error('n must be positive.');
    end
    % n must be an integer
    if round(n)-n~=0
        error('n must be an integer.');
    end
end

if nargin<2 | isempty(level)==1
    level=0;
else
    % level must be a scalar
    if sum(size(level))>2
        error('level must be scalar.');
    end
    % level must be positive
    if level<0
        error('level must be positive.');
    end
end

if nargin<3 | isempty(a)==1
    a=4;
else
    % a must be scalar
    if sum(size(a))>2
        error('a must be scalar.');
    end
end

if nargin<4 | isempty(x0)==1
    x0=0.1;
else
    % x0 must be scalar
    if sum(size(x0))>2
        error('x0 must be scalar.');
    end
end

% Initialize
x(1,1)=a*x0*(1-x0);

% Simulate
for i=2:n
    x(i,1)=a*x(i-1,1)*(1-x(i-1,1));
end

% Add normal white noise
x=x+randn(n,1)*level*std(x);

author: Sable

repo: mcbench-benchmarks


text:
% imbibition capillary pressure curve
% Written by Ali A. Eftekhari
% Note that pc_max_o is specified as pc_min in the input json files for
% slightly more consistency witht the polynomial Pc
function res=pc_imb(sw, pce_w, pce_o, swc, sor, labda_w, labda_o, pc_max_w, pc_max_o)
  pc1=pc_drain(sw, pce_w, swc, labda_w, pc_max_w);
  pc2=pc_drain(1-sw, pce_o, sor, labda_o, pc_max_o);
  res=pc1-pc2;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
clear all; close all; clc;

t=linspace(-pi,pi,10000);
delay=-pi/2;
amp1=0.8;
amp2=1;
slowness=150;

for ii=1:1:length(t)

    A=amp1*sin(amp1*2*pi*1*(ii/slowness));
    B=amp2*sin(amp2*2*pi*1*(ii/slowness)+delay);

    sin_graph=amp1*sin(amp1*2*pi*1*(t+ii/slowness));
    cos_graph=amp2*sin(amp2*2*pi*1*(t+ii/slowness)+delay);
    subplot(311);
    plot(t,sin_graph);hold on;
    plot(0,A,'r*','Markersize',10);hold off;title('Signal A : X-Axis');
    xlim([-pi pi]);ylim([-1 1]);
    set(gcf,'Position',[250,100, 800, 800]);
    grid on;
    subplot(312);
    T=num2str(delay);
    Title_for_B=strcat('Signal B with Delay = ',T,' : Y-Axis');
    plot(t,cos_graph);hold on;title(Title_for_B);
    xlim([-pi pi]);ylim([-1 1]);
    plot(0,B,'r*','MarkerSize',10);hold off;
    grid on;
    subplot(338);
    T2=num2str(amp1/amp2);
    Title_for_C=strcat('Ratio A/B = ',T2);
    plot(sin_graph,cos_graph);title(Title_for_C);
    hold on;
    plot(A,B,'r*','MarkerSize',10);hold off;xlim([-1 1]); ylim([-1 1]);
    grid on;
%     drawnow;
    pause(0.00000001)

end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function test_bug2560

% MEM 2gb
% WALLTIME 00:15:00
% DEPENDENCY ft_sourceplot ft_plot_slice

% this function is meant to test the correctness of the 'ortho' method in
% sourceplot. It relies on ft_plot_ortho and ft_plot_slice. there may be a
% slight problem in dealing with the functional data/mask

% investigate the issue by creating some data
[x,y,z] = ndgrid(1:7,1:8,1:9);
inside  = zeros(7,8,9);
inside(2:6,2:7,2:8) = 1;
outside = find(inside==0);
inside  = find(inside==1);

source = [];
source.pos = [x(:) y(:) z(:)]; clear x y z
source.dim = [7 8 9];
source.inside  = inside;
source.outside = outside;
source.avg.pow = zeros(7,8,9);
source.avg.pow(inside) = 1;

cfg = [];
cfg.funparameter = 'avg.pow';
cfg.method       = 'ortho';
cfg.funcolormap  = 'jet';
ft_sourceplot(cfg, source);

% the code above confirms an the hunch I had: the inside mask is shifted
% relative to the functional data


text:
function writeAllLabeledImages(imdir, imsegs, pg, outdir)
% writeAllLabeledImages(imdir, imsegs, pg, outdir)
%
% Write label images to specified outdir
%

DO_GT = 0;
if isempty(pg)
    disp('Writing ground truth')
    DO_GT = 1;
end

for f = 1:numel(imsegs)

    if DO_GT
        nclasses =  numel(imsegs(f).label_names);
        pg{f} = zeros(imsegs(f).nseg, nclasses);
        lab = imsegs(f).labels;
        ind = find(lab>0);
        pg{f}(ind + imsegs(f).nseg*(lab(ind)-1)) = 1;
    end

    im = im2double(imread([imdir '/' imsegs(f).imname]));
    disp([num2str(f) ': ' imsegs(f).imname])
    [vc, hc] = splitpg(pg{f});

    lim = APPgetLabeledImage2(im, imsegs(f), vc, hc);

    %system(['cp ' imdir '/' imsegs(f).imname ' ' outdir '/' imsegs(f).imname]);
    imwrite(lim, [outdir '/' strtok(imsegs(f).imname, '.') '.l.jpg'], 'Quality', 90);
end

author: zouchuhang

repo: LayoutNet


text:
function subset_check ( n, t )

%*****************************************************************************80
%
%% SUBSET_CHECK checks a subset.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of elements in the master set.
%    N must be positive.
%
%    Input, integer T(N), the subset.  If T(I) = 0, item I is
%    not in the subset; if T(I) = 1, item I is in the subset.
%
  if ( n < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SUBSET_COLEX_RANK - Fatal error!\n' );
    fprintf ( 1, '  N = %d < 1.\n', n );
    error ( 'SUBSET_COLEX_RANK - Fatal error!' );
  end

  for i = 1 : n

    if ( t(i) ~= 0 && t(i) ~= 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SUBSET_COLEX_RANK - Fatal error!\n' );
      fprintf ( 1, '  T(%d) = %d, but must be 0 or 1.\n', i, t(i) );
      error ( 'SUBSET_COLEX_RANK - Fatal error!' );
    end

  end

  return
end


text:
function b = isappr(x,y,varargin)
% check double == double

b = isnull(abs(x-y),varargin{:});


text:
function [llh] = tapas_sem_multiv_llh(data, theta, ptheta)
%% Compute likelihood of the eye movement model in a hierarchical format.
%
%   Input
%
%   data        -- Stucture array of dimension Nx1 with fields y and u
%   theta       -- Structure array of dimension 1 X 1
%   ptheta      -- Priors
%
%   Ouput
%
%   llh         -- Matrix of dimensions NxM where N is the number fo subjects
%                   and M the number of chains.
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%

ns = size(theta.y, 1);
nc = size(theta.y, 2);

atheta = cell(ns, nc);

for i = 1:ns
    for j = 1:nc
        atheta{i, j} = ptheta.model.p0 + ptheta.model.jm * theta.y{i, j};
    end
end

llh =  ptheta.model.llh(data, atheta);

end

author: translationalneuromodeling

repo: tapas


text:
function enter_ax_an_angs

global phi theta psi rad deg x y z gamma delta alpha

naner=false;

gamman=str2num(get(gamma,'String'));
if length(gamman)==0
    naner=true;
end

deltan=str2num(get(delta,'String'));
if length(deltan)==0
    naner=true;
end

alphan=str2num(get(alpha,'String'));
if length(alphan)==0
    naner=true;
end

if naner
    nan_error;
else

    if get(deg,'Value')
        gamman=pi*gamman/180;
        deltan=pi*deltan/180;
        alphan=pi*alphan/180;
    end

    an=ax_an_bounding(gamman,deltan,alphan);

    if an(1)
        gamman=an(2);
        deltan=an(3);
        alphan=an(4);
        if get(deg,'Value')
            set(gamma,'String',num2str(180*gamman/pi));
            set(delta,'String',num2str(180*deltan/pi));
            set(alpha,'String',num2str(180*alphan/pi));
        else
            set(gamma,'String',num2str(gamman));
            set(delta,'String',num2str(deltan));
            set(alpha,'String',num2str(alphan));
        end
    end

    an=axan2euler(gamman,deltan,alphan); % convert to Euler angles
    an1=an{1};
    phin=an1{1};
    thetan=an1{2};
    psin=an1{3};

    if get(deg,'Value')
        set(phi,'String',num2str(180*phin/pi));
        set(theta,'String',num2str(180*thetan/pi));
        set(psi,'String',num2str(180*psin/pi));
    else
        set(phi,'String',num2str(phin));
        set(theta,'String',num2str(thetan));
        set(psi,'String',num2str(psin));
    end

    vn=an{2};
    set(x,'string',num2str(vn(1)));
    set(y,'string',num2str(vn(2)));
    set(z,'string',num2str(vn(3)));

    %enter_eul_angs_1;
    enter_eul_angs(true,gamman,deltan,alphan); % with not set v

end


text:
function infoList = removeDuplicateNames(infoList)
% Accepts a cell array of strings, with info seperated
% by a '|' within each cell, and removes the duplicate names in each cell.
% Also orders the names by size from smallest to largest. Best used after
% `fixNames`.
% Called by `cleanTmodel`, `verifyModel`, `addSEEDInfo`.
%
% USAGE:
%
%    infoList = removeDuplicateNames(infoList)
%
% INPUTS:
%    infoList:   Cell array of strings.
%
% OUTPUTS:
%    infoList:   After removal
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

for iInfo = 1:length(infoList) % Get rid of those duplicates!
    info = infoList{iInfo} ;
    newInfo = '' ;
    % Break information into parts.
    pipePos = [0 strfind(info, '|') length(info) + 1] ;
    newInfoCell = {[]} ;
    newInfoLength = 0 ;
    for iP = 1:length(pipePos) - 1
        nowInfo = info(pipePos(iP) + 1:pipePos(iP + 1) - 1) ;
        % Check to see if the information exists.
        if isempty(find(strcmp(nowInfo, newInfoCell), 1))
            % If not, record the data and it's length
            if isempty(newInfoCell{1})
                newInfoCell{1} = nowInfo ;
                newInfoLength(1) = length(nowInfo) ;
            else
                newInfoCell{length(newInfoCell) + 1, 1} = nowInfo ;
                newInfoLength(length(newInfoCell), 1) = length(nowInfo) ;
            end
        end
    end

    % Sort the data based on length.
    [null, lengthOrder] = sort(newInfoLength) ;
    newInfoCell = newInfoCell(lengthOrder) ;

    % Assemble cell array into single string.
    newInfo = newInfoCell{1} ;
    for iInfoPart = 2:length(newInfoCell)
        newInfo = [newInfo '|' newInfoCell{iInfoPart}] ;
    end

    % Replace the information the updated string.
    infoList{iInfo} = newInfo ;
end

author: opencobra

repo: cobratoolbox


text:
function [blk, At, b] = sdpt3_create_hpsd_blk(n)
%SDPT3_CREATE_HPSD_BLK Creates a Hermitian PSD matrix variable for the
%SDPT3 solver.
%Note that
%   W \succeq 0 \iff U = [ Re(W)  Im(W) ] \succeq 0
%                      [ -Im(W) Re(W) ]
%We declare a 2nx2n symmetric matrix variable U and enforce additional
%constraints according to the structure of U:
%   1. U(i+n,j+n) = U(i,j) \forall 1 <= i <= j <= n
%   2. U(i,j+n) = -U(j,i+n) \forall 1 <= i <= j <= n
%Syntax:
%   [blk, At, b] = SDPT3_CREATE_HPSD_BLK(n);
%Inputs:
%   n - Dimension of the matrix.
%Outputs:
%   blk - Block definition: {'s', 2*n}.
%   At - Transposed constraint matrix (as a submatrix of the full
%        constraint matrix).
%   b - Constraint value vector (as a subvector of the full constraint
%       value vector).
blk = {'s', 2*n};
n_c = n * (n + 1);
n_var = n * (2 * n + 1);
n_sub = n * (n + 1) / 2;
n_nz = 2 * n_c - n;
% triplets for the sparse constraint matrix
idx_row = zeros(n_nz, 1);
idx_col = zeros(n_nz, 1);
vals = zeros(n_nz, 1);
i_c = 0;
i_nz = 0;
% constraints are expressed by
%   A svec(W) = b
%   W(i,j) --> svec(W)(i + j(j-1)/2)
% Re(W) = Re(W)
for jj = 1:n
    for ii = 1:jj
        i_c = i_c + 1;
        i_nz = i_nz + 1;
        idx_row(i_nz) = i_c;
        idx_col(i_nz) = ii + jj*(jj-1)/2;
        vals(i_nz) = 1;
        i_nz = i_nz + 1;
        idx_row(i_nz) = i_c;
        idx_col(i_nz) = ii + n + (jj + n)*(jj + n - 1)/2;
        vals(i_nz) = -1;
    end
end
% Im(W) = -Im(W)^T
% diagonals should be zeros
for ii = 1:n
    i_c = i_c + 1;
    i_nz = i_nz + 1;
    idx_row(i_nz) = i_c;
    idx_col(i_nz) = ii + (ii + 2*n)*(ii - 1)/2 + n_sub;
    vals(i_nz) = 1;
end
% skew symmetry
for jj = 2:n
    for ii = 1:jj-1
        i_c = i_c + 1;
        i_nz = i_nz + 1;
        idx_row(i_nz) = i_c;
        idx_col(i_nz) = ii + (jj + 2*n)*(jj - 1)/2 + n_sub;
        vals(i_nz) = 1;
        i_nz = i_nz + 1;
        idx_row(i_nz) = i_c;
        idx_col(i_nz) = jj + (ii + 2*n)*(ii - 1)/2 + n_sub;
        vals(i_nz) = 1;
    end
end
At = sparse(idx_row, idx_col, vals, n_c, n_var)';
b = sparse(n_c, 1);
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
%
% write_matrix(FID,data)
%
% This function write a matrix to the specified file. Both two dimensional
% matrix and one dimensional vector are accepted.
%
% Author: Xiao Xiong
% Created: 6 Feb, 2007
% Last modified: 6 Feb, 2007

function write_matrix_file(filename, data)

FID = fopen(filename, 'w');
write_matrix(FID, data);
fclose(FID);

author: singaxiong

repo: SignalGraph


text:
function [alpha,ro] = fit_corr(corr)

kmax = length(corr);

[x eval] = fminsearch(@fit_model, [0, 1/2], optimset('TolX',1e-8),kmax,corr);

alpha = x(1);
ro = x(2);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function obj = cfg_inv_out(varargin)

% function obj = cfg_inv_out(varargin)
% Auxiliary class to mark invalid (i.e. not yet available) outputs of
% cfg_exbranch'es. An object of this type will be assigned automatically
% to a cfg_exbranch'es .jout field. resolve_deps will not resolve outputs
% that consist of a cfg_inv_out object.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_inv_out.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

% Dummy field - nothing should be stored here
obj.invalid = true;
obj = class(obj, 'cfg_inv_out');

author: spm


text:
function pos=read_pascal(cls,year,keypoints,segmentations)

% Read pascal to generate filtered instances of classes

globals;
dataset_fg = 'trainval';
conf.year=year;
conf.dev_kit=[BASE_DIR 'data/VOC2011/VOCdevkit'];
VOCopts    = get_voc_opts(conf);

% Setup keypoints and segmentations voc_id
kps_voc_id = cellfun(@(x,y) [x, '_', y], keypoints.voc_image_id, arrayfun(@num2str, keypoints.voc_rec_id, 'Uniform', false), 'UniformOutput', false);
segs_voc_id = cellfun(@(x,y) [x, '_', y], segmentations.voc_image_id, arrayfun(@num2str, segmentations.voc_rec_id, 'Uniform', false), 'UniformOutput', false);

  % Positive examples from the foreground dataset
  ids      = textread(sprintf(VOCopts.imgsetpath, dataset_fg), '%s');
  pos      = [];
  numpos   = 0;
  for i = 1:length(ids);
    if(mod(i,500)==0)
        fprintf('Parsed ids %d/%d\n',i,length(ids));
    end
    % Parse record and exclude difficult examples
    rec           = PASreadrecord(sprintf(VOCopts.annopath, ids{i}));
    clsinds       = strmatch(cls, {rec.objects(:).class}, 'exact');
    count         = length(clsinds(:));
    % Skip if there are no objects in this image
    if count == 0
      continue;
    end

    diff          = [rec.objects(clsinds).difficult];
    trunc         = [rec.objects(clsinds).truncated];
    if(~isempty(clsinds) && isfield(rec.objects(clsinds(1)),'occluded'))
        occl = [rec.objects(clsinds).occluded];
    else
        occl = true(size(clsinds));
    end

    clsinds2del = diff | trunc | occl; % Change this for different datasets
    %clsinds2del = diff;
    clsinds(clsinds2del)=[];

    count         = length(clsinds(:));
    % Skip if there are no objects in this image
    if count == 0
      continue;
    end

    % Create one entry per bounding box in the pos array
    for j = clsinds(:)'
      numpos = numpos + 1;
      j_voc_id = [rec.filename(1:end-4) '_' num2str(j)];
      ki = find(ismember(kps_voc_id,j_voc_id));
      si = find(ismember(segs_voc_id,j_voc_id));
      bbox   = rec.objects(j).bbox;
      pos(numpos).imsize = [rec.size.height rec.size.width];
      pos(numpos).voc_image_id = rec.filename(1:end-4);
      pos(numpos).voc_rec_id = j;
      %pos(numpos).im      = [VOCopts.datadir rec.imgname];
      pos(numpos).pascal_bbox   = [bbox(1:2) bbox(3)-bbox(1)+1 bbox(4)-bbox(2)+1]; %(x,y,w,h)
      pos(numpos).view    = rec.objects(j).view;
      pos(numpos).kps     = squeeze(keypoints.coords(ki,:,:));
      pos(numpos).part_names  = keypoints.labels;
      pos(numpos).bbox        = keypoints.bbox(ki,:);
      pos(numpos).poly_x      = segmentations.poly_x{si};
      pos(numpos).poly_y      = segmentations.poly_y{si};
      pos(numpos).class       = cls;
    end
  end
  pascal_data= pos;
  fname = [cachedir 'pascal_filt_data/pascal_' cls '_filt.mat'];
  save(fname,'pascal_data');
end

function VOCopts = get_voc_opts(conf)
% cache VOCopts from VOCinit
persistent voc_opts;

key = conf.year;
if isempty(voc_opts) || ~voc_opts.isKey(key)
  if isempty(voc_opts)
    voc_opts = containers.Map();
  end
  tmp = pwd;
  cd(conf.dev_kit);
  addpath([cd '/VOCcode']);
  VOCinit;
  cd(tmp);
  voc_opts(key) = VOCopts;
end
VOCopts = voc_opts(key);
end

repo: CategoryShapes


text:
function line_num = sphere_cubed_line_num ( n )

%*****************************************************************************80
%
%% SPHERE_CUBED_LINE_NUM counts lines on a cubed sphere grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of sections into which each face of
%    the cube is to be divided.
%
%    Output, integer LINE_NUM, the number of lines.
%
  line_num = 0;
%
%  If N = 1, the corners form 12 lines.
%
  if ( n == 1 )
    line_num = 12;
    return
%
%  If 1 < N, each of 8 corners connects to three neighboring edges.
%
  else
    line_num = line_num + 8 * 3;
  end
%
%  If 2 < N, then each of the 12 edges includes lines.
%
  if ( 2 < n )
    line_num = line_num + 12 * ( n - 2 );
  end
%
%  Lines that belong to one of the six faces.
%
  if ( 1 < n )
    line_num = line_num + 6 * 2 * n * ( n - 1 );
  end

  return
end

repo: jburkardt-m