text:
%%*******************************************************************
%% Prod2: compute the block diagonal matrix A*B
%%
%%    C = Prod2(blk,A,B,options);
%%
%% INPUT:  blk = a cell array describing the block structure of A and B
%%         A,B = square matrices or column vectors.
%%
%%   options = 0  if no special structure
%%             1  if C is symmetric
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

function  C = Prod2(blk,A,B,options)

global spdensity

if (nargin == 3); options = 0; end;
iscellA = iscell(A); iscellB = iscell(B);
%%
if (~iscellA && ~iscellB)
    if (size(blk,1) > 1);
        error('Prod2: blk and A,B are not compatible');
    end;
    if strcmp(blk{1},'s')
        numblk = length(blk{2});
        isspA = issparse(A); isspB = issparse(B);
        if (numblk > 1)
            if ~isspA; A=sparse(A); isspA=1; end
            if ~isspB; B=sparse(B); isspB=1; end
        end
        %%use_matlab = (options==0 && ~isspA && ~isspB) || (isspA && isspB);
        use_matlab = (~isspA && ~isspB) || (isspA && isspB);
        if (use_matlab)
            C = A*B;
            if (options==1); C = 0.5*(C+C'); end;
        else
            C = mexProd2(blk,A,B,options);
        end
        checksparse = (numblk==1) && (isspA || isspB);
        if (checksparse)
            n2 = sum(blk{2}.*blk{2});
            if (mexnnz(C) <= spdensity*n2);
                if ~issparse(C); C = sparse(C); end;
            else
                if issparse(C); C = full(C); end;
            end
        end
    elseif (strcmp(blk{1},'q') || strcmp(blk{1},'l') || strcmp(blk{1},'u'))
        C = A.*B;
    end
else
    error('Prod2: A,B must be matrices');
end
%%*******************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% this script is to add metabolite identifiers (#174) that will be used for Met Atlas web-portal.

% load data
metIDs = importTsvFile('metabolites.tsv');

% get compartment id for non-HMR mets that will be used later
compID = cellfun(@(x) regexprep(x, '^.+_(\w)$', '$1'), metIDs.mets, 'UniformOutput', false);
matche2NonHMR = regexp(metIDs.mets, '_\w$', 'match');
emptymask = cellfun('isempty', matche2NonHMR);
compID(emptymask) = {''};    % empty HMR met ids

% generate MA met ids from HMR ids only with "MA-" prefix
metMAID = cellfun(@(x) regexprep(x, '^m(\d\d\d\d\d\w)$', 'MA-M$1'), metIDs.mets, 'UniformOutput', false);

% check consistency between HMR ids and empty elemnets of non-HMR ids
matche2HMR = regexp(metMAID, '^MA-M\d\d\d\d\d\w$', 'match');
ind2HMR = ~cellfun('isempty', matche2HMR);
check = find(ind2HMR ~= emptymask);
metIDs.mets(check)
%    {'temp001c'}
%    {'temp001s'}
% these two mets are from HMR but in non-standard format, will be treated as
% non-HMR ids
compID(check) = {'c';'s'};

% kepp HMR ids in metMAID and clean others
ind_MAID = startsWith(metMAID, 'MA-M');
metMAID(~ind_MAID) = {''};
metMAIDNoComp = cellfun(@(x) regexprep(x, '^(.+)\w$', '$1'), metMAID, 'UniformOutput', false);

% get index of nonHMR ids from compID and remove their comppart id
nonHMRidInd = ~cellfun('isempty', compID);
nonHMRidNoComp = cellfun(@(x) regexprep(x, '^(.+)_\w$', '$1'), metIDs.mets(nonHMRidInd), 'UniformOutput', false);

% get unique list of nonHMR ids
[uniqueID, ia, ic] = unique(nonHMRidNoComp);
%isequal(nonHMRidNoComp, uniqueID(ic))
%isequal(uniqueID, nonHMRidNoComp(ia))

% prepare standard format MA ids for nonHMR mets
idsToAdd = cell(10000,1);
for i=1:10000
    idsToAdd{i} = strcat('MA-M',sprintf('%05d', i));
end
idsToAdd = setdiff(idsToAdd, metMAIDNoComp);     % exclude existing ones
idsToAdd = idsToAdd(1:length(uniqueID));

% get ids to fill in
idsToFill = idsToAdd(ic);

% append compartment id to idsToFill
idsToFill = strcat(idsToFill,compID(nonHMRidInd));

% complete MA met ids
metMAID(nonHMRidInd) = idsToFill;

% format check
all(cell2mat(regexp(metMAID, '^MA-M\d\d\d\d\d\w$')))
% unique check
isequal(length(metMAID), length(sort(metMAID)))
% empty check
B = cellfun('isempty', metMAID);
all(B(:) == 0)

% add new field and output
metIDs.metMAID = metMAID;
exportTsvFile(metIDs,'metabolites.tsv');

author: SysBioChalmers

repo: Human-GEM


text:
function [ vX, mP ] = ApplyKalmanFilterIteration( vX, mP, vZ, hF, hH, hMf, hMh, mQ, mR )
% ----------------------------------------------------------------------------------------------- %
% [ vX, mP ] = ApplyKalmanFilterIteration( vX, mP, vZ, hF, hH, mQ, mR, mF, mH )
%   Applies iteration of the Kalman Filter (Prediction + Update). Supports
%   both Linear and Non Linear Mode (Extended Kalman Filter).
% Input:
%   - vX            -   Input State Vector.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - mP            -   Error Covariance.
%                       Structure: Matrix.
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - vZ            -   Measurement Vector.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - hF            -   Model Function.
%                       Structure: Function Handler.
%                       Type: Function Handler.
%                       Range: NA.
%   - hH            -   Measurement Function.
%                       Structure: Function Handler.
%                       Type: Function Handler.
%                       Range: NA.
%   - mQ            -   Process Noise Covariance.
%                       Structure: Matrix.
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - mR            -   Measurement Noise Covariance.
%                       Structure: Matrix.
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - hMf           -   Model Matrix Function.
%                       Generates the matrix `mF` from the state vector vX.
%                       Structure: Function Handler.
%                       Type: Function Handler.
%                       Range: NA.
%   - hMh           -   Measurement Matrix Function.
%                       Generates the matrix `mH` from the state vector vX.
%                       The state vector is hF(vX).
%                       Structure: Function Handler.
%                       Type: Function Handler.
%                       Range: NA.
% Output:
%   - vX            -   Updated State Vector.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - mP            -   Error Covariance.
%                       Structure: Matrix.
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  Extended Kalman Filter (Wikipedia) - https://en.wikipedia.org/wiki/Extended_Kalman_filter.
%   2.  Kalman Filter (Wikipedia) - https://en.wikipedia.org/wiki/Kalman_filter.
% Remarks:
%   1.  If the Complex Mode is selected the function must return complex
%       values in order to work. For instance, if the input function is
%       'norm(vX)' use 'sum(vX .^ 2)' and if the input function is
%       'sum(abs(vX))' use 'sum(sqrt(vX .^ 2))'.
%   2.  This implementation use the Joseph Form of the Covariance Update.
%       While it requires more computational effort it is more stable and
%       guaranteed to generate PSD matrix.
% TODO:
%   1.  U.
% Release Notes:
%   -   1.1.000     25/07/2021  Royi Avital
%       *   Added support for a fnction handle to calculate the matrices mF
%           and mH.
%       *   Add `arguments` block.
%       *   Added a step to ensure symmetry of mP.
%   -   1.0.000     22/08/2018  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

arguments
    vX (:, 1) {mustBeNumeric, mustBeReal}
    mP (:, :) {mustBeNumeric, mustBeReal, mustBePdMatrix}
    vZ (:, 1) {mustBeNumeric, mustBeReal}
    hF (1, 1) {mustBeFunctionHandler}
    hH (1, 1) {mustBeFunctionHandler}
    hMf (1, 1) {mustBeFunctionHandler}
    hMh (1, 1) {mustBeFunctionHandler}
    mQ (:, :) {mustBeNumeric, mustBeReal, mustBePdMatrix}
    mR (:, :) {mustBeNumeric, mustBeReal, mustBePdMatrix}
end

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

mI = eye(size(vX, 1));

% Prediction Step
vX = hF(vX);
mF = hMf(vX);
mP = (mF * mP * mF.') + mQ;

mH = hMh(vX);

% Update Step
vY = vZ - hH(vX);
mS = (mH * mP * mH.') + mR;
mK = (mP * mH.') / mS;
vX = vX + (mK * vY);
% Joseph Form of the Covariance Update (Numerically more stable)
mT = mI - (mK * mH);
mP = (mT * mP * mT.') + (mK * mR * mK.'); %<! Joseph Form
mP = (mP.' + mP) / 2; %<! Ensure symmetry

end

function [ ] = mustBePdMatrix( mX )
% https://www.mathworks.com/matlabcentral/answers/107552
if(issymmetric(mX) && (all(eig(mX) > 0)))

else
    errorId     = 'mustBePdMatrix:notPositiveDefiniteMatrix';
    errorMsg    = 'The input must be a Positive Definite Matrix';
    throwAsCaller(MException(errorId, errorMsg));
end

end

function [ ] = mustBeFunctionHandler( hF )
% https://www.mathworks.com/matlabcentral/answers/107552
if(~isa(hF, 'function_handle'))
    errorId     = 'mustBeFunctionHandler:notFunctionHandler';
    errorMsg    = 'The input must be a Function Handler';
    throwAsCaller(MException(errorId, errorMsg));
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
% gsn - training an GSN (stochastic backprop)
% Copyright (C) 2013 KyungHyun Cho
%
% 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function [G] = gsn(G, patches, n_walkback, valid_patches, valid_portion);

if nargin < 3
    n_walkback = G.structure.n_layers * 2;
end

valid_err = 0;

if nargin < 4
    early_stop = 0;
    valid_patches = [];
    valid_portion = 0;
else
    early_stop = 1;
    valid_err = -Inf;
    valid_best_err = -Inf;
    valid_violate_cnt = 0;
end

actual_lrate = G.learning.lrate;

n_samples = size(patches, 1);

layers = G.structure.layers;
n_layers = length(layers);

if layers(1) ~= size(patches, 2)
    error('Data is not properly aligned');
end

minibatch_sz = G.learning.minibatch_sz;
n_minibatches = ceil(n_samples / minibatch_sz);

n_epochs = G.iteration.n_epochs;

momentum = G.learning.momentum;
weight_decay = G.learning.weight_decay;

biases_grad = cell(n_layers, 1);
W_grad = cell(n_layers, 1);
biases_grad_old = cell(n_layers, 1);
W_grad_old = cell(n_layers, 1);
for l = 1:n_layers
    biases_grad{l} = zeros(size(G.biases{l}))';
    if l < n_layers
        W_grad{l} = zeros(size(G.W{l}));
    end
    biases_grad_old{l} = zeros(size(G.biases{l}))';
    if l < n_layers
        W_grad_old{l} = zeros(size(G.W{l}));
    end
end

min_recon_error = Inf;
min_recon_error_update_idx = 0;
stopping = 0;

do_normalize = G.do_normalize;
do_normalize_std = G.do_normalize_std;

if G.data.binary == 0
    if do_normalize == 1
        % make it zero-mean
        patches_mean = mean(patches, 1);
        patches = bsxfun(@minus, patches, patches_mean);
    end

    if do_normalize_std ==1
        % make it unit-variance
        patches_std = std(patches, [], 1);
        patches = bsxfun(@rdivide, patches, patches_std);
    end
end

anneal_counter = 0;
actual_lrate0 = actual_lrate;

rerr_ma = 0;

% try
%     use_gpu = gpuDeviceCount;
% catch errgpu
%     use_gpu = false;
%     disp(['Could not use CUDA. Error: ' errgpu.identifier])
% end
use_gpu = 0;

% figure;

for step=1:n_epochs
    if G.verbose
        fprintf(2, 'Epoch %d/%d: ', step, n_epochs)
    end
    if use_gpu
        % push
        G = push_to_gpu (G);
    end

    for mb=1:n_minibatches
        G.iteration.n_updates = G.iteration.n_updates + 1;

        mb_start = (mb - 1) * minibatch_sz + 1;
        mb_end = min(mb * minibatch_sz, n_samples);

        % p_0
        v0 = patches(mb_start:mb_end, :);
        mb_sz = size(v0,1);

        if use_gpu > 0
            v0 = gpuArray(single(v0));
        end

        % add error
        v0_clean = v0;

        % forward pass
        h0 = cell(n_layers, n_walkback * 2 + 1);
        v1clean = cell(1, n_walkback * 2 + 1);
        occupied = zeros(n_layers, n_walkback * 2 + 1);

        h0{1,1} = v0;
        v1clean{1,1} = v0;
        occupied(1,1) = 1;

        for wi = 1:(n_walkback * 2 + 1)
            for l = 1:min(wi,n_layers)
                if sum(occupied(1, 2:end)) == n_walkback
                    break;
                end

                if l == 2
                    if G.data.binary == 0 && G.noise.level > 0
                        if use_gpu
                            h0{l-1,wi-1} = h0{l-1,wi-1} + G.noise.level * parallel.gpu.GPUArray.randn(size(h0{l-1,wi-1}));
                        else
                            h0{l-1,wi-1} = h0{l-1,wi-1} + G.noise.level * randn(size(h0{l-1,wi-1}));
                        end
                    end

                    if G.noise.drop > 0
                        mask = binornd(1, 1-G.noise.drop, size(h0{l-1,wi-1}));
                        if G.data.binary
                            sandp = binornd(1, 0.5, size(h0{l-1,wi-1}));
                        else
                            sandp = zeros(size(h0{l-1,wi-1}));
                        end
                        h0{l-1,wi-1} = h0{l-1,wi-1} .* mask + (1 - mask) .* sandp;
                        clear mask;
                    end
                end

                if l == 1 && wi == 1
                    continue;
                end

                prev_exist = 0;

                if l > 1 && occupied(l-1, wi-1) == 1
                    prev_exist = prev_exist + 1;
                end
                if l < n_layers && occupied(l+1, wi-1) == 1
                    prev_exist = prev_exist + 1;
                end

                if prev_exist == 0
                    continue;
                end

                occupied(l, wi) = 1;
                h0{l, wi} = zeros(mb_sz, G.structure.layers(l));
                if l > 1 && occupied(l-1, wi-1) == 1
                    h0{l, wi} = h0{l, wi} + h0{l-1, wi-1} * G.W{l-1};
                end
                if l < n_layers && occupied(l+1, wi-1) == 1
                    h0{l, wi} = h0{l, wi} + h0{l+1, wi-1} * G.W{l}';
                end
                h0{l, wi} = bsxfun(@plus, h0{l, wi}, G.biases{l}');

                % pre-sigmoid noise
                if l > 1
                    if G.hidden.add_noise(l)
                        if use_gpu
                            h0{l, wi} = h0{l, wi} + G.hidden.noise_level * parallel.gpu.GPUArray.randn(mb_sz, G.structure.layers(l));
                        else
                            h0{l, wi} = h0{l, wi} + G.hidden.noise_level * randn(mb_sz, G.structure.layers(l));
                        end
                    end
                end

                if l == 1
                    if G.data.binary
                        h0{l, wi} = sigmoid(h0{l, wi});
                    end
                    v1clean{l, wi} = h0{l, wi};
                else
                    h0{l, wi} = sigmoid(h0{l, wi}, G.hidden.use_tanh);
                end

                % FIXME: unable to train a GSN with post-sigmoid noise.
%                % post-sigmoid noise
%                if l > 1
%                    if G.hidden.add_noise(l)
%                        if use_gpu
%                            h0{l, wi} = h0{l, wi} + G.hidden.noise_level * parallel.gpu.GPUArray.randn(mb_sz, G.structure.layers(l));
%                        else
%                            h0{l, wi} = h0{l, wi} + G.hidden.noise_level * randn(mb_sz, G.structure.layers(l));
%                        end
%                    end
%                end
            end

            if sum(occupied(1, 2:end)) == n_walkback
                break;
            end
        end

        % reset gradients
        for l = 1:n_layers
            biases_grad{l} = 0 * biases_grad{l};
            if l < n_layers
                W_grad{l} = 0 * W_grad{l};
            end
        end

        % error backprop
        delta = cell(n_layers, n_walkback * 2 + 1);

        rerr = 0;
        for wi = 2:(n_walkback * 2 + 1)
            if occupied(1, wi)
                rerr = rerr + mean(sum((v1clean{1, wi} - v0_clean).^2,2));
            end
        end

        if use_gpu > 0
            rerr = gather(rerr);
        end

        rerr_ma = rerr * 0.1 + rerr_ma * 0.9;

        G.signals.recon_errors = [G.signals.recon_errors rerr];

        for wi = (n_walkback * 2 + 1):-1:1
            for l = 1:min(wi,n_layers)
                if occupied(l, wi) == 0
                    continue;
                end

                delta{l, wi} = zeros(mb_sz, G.structure.layers(l));

                if wi < (n_walkback * 2 + 1)
                    if l > 1 && occupied(l-1,wi+1) == 1
                        delta{l, wi} = delta{l, wi} + delta{l-1,wi+1} * G.W{l-1};
                    end
                    if l < n_layers && occupied(l+1,wi+1) == 1
                        delta{l, wi} = delta{l, wi} + delta{l+1,wi+1} * G.W{l}';
                    end
                end

                if l == 1
                    if G.data.binary && wi < (n_walkback * 2 + 1)
                        delta{l, wi} = delta{l, wi} .* dsigmoid(h0{l, wi});
                    end

                    delta{l, wi} = delta{l, wi} + (v1clean{l, wi} - v0_clean);
                    %delta{l, wi} = delta{l, wi} + (h0{l, wi} - v0_clean);
                else
                    delta{l, wi} = delta{l, wi} .* dsigmoid(h0{l, wi}, G.hidden.use_tanh);
                end

                if wi > 1
                    biases_grad{l} = biases_grad{l} + mean(delta{l, wi}, 1);
                    if l > 1 && occupied(l-1,wi-1) == 1
                        W_grad{l-1} = W_grad{l-1} + (h0{l-1, wi-1}' * delta{l, wi})/size(v0,1);
                    end
                    if l < n_layers && occupied(l+1,wi-1) == 1
                        W_grad{l} = W_grad{l} + (delta{l, wi}' * h0{l+1, wi-1})/size(v0,1);
                    end
                end
            end
        end

        % learning rate
        if G.adadelta.use
            % update
            for l = 1:n_layers
                biases_grad_old{l} = (1 - momentum) * biases_grad{l} + momentum * biases_grad_old{l};
                if l < n_layers
                    W_grad_old{l} = (1 - momentum) * W_grad{l} + momentum * W_grad_old{l};
                end
            end

            if G.iteration.n_updates == 1
                adamom = 0;
            else
                adamom = G.adadelta.momentum;
            end

            for l = 1:n_layers
                if l < n_layers
                    G.adadelta.gW{l} = adamom * G.adadelta.gW{l} + (1 - adamom) * W_grad_old{l}.^2;
                end

                G.adadelta.gbiases{l} = adamom * G.adadelta.gbiases{l} + (1 - adamom) * biases_grad_old{l}.^2';
            end

            for l = 1:n_layers
                dbias = -(biases_grad_old{l}' + ...
                    weight_decay * G.biases{l}) .* (sqrt(G.adadelta.biases{l} + G.adadelta.epsilon) ./ ...
                    sqrt(G.adadelta.gbiases{l} + G.adadelta.epsilon));
                G.biases{l} = G.biases{l} + dbias;

                G.adadelta.biases{l} = adamom * G.adadelta.biases{l} + (1 - adamom) * dbias.^2;
                clear dbias;

                if l < n_layers
                    dW = -(W_grad_old{l} + ...
                        weight_decay * G.W{l}) .* (sqrt(G.adadelta.W{l} + G.adadelta.epsilon) ./ ...
                        sqrt(G.adadelta.gW{l} + G.adadelta.epsilon));
                    G.W{l} = G.W{l} + dW;

                    G.adadelta.W{l} = adamom * G.adadelta.W{l} + (1 - adamom) * dW.^2;

                    clear dW;
                end

            end
        else
            if G.learning.lrate_anneal > 0 && (step >= G.learning.lrate_anneal * n_epochs)
                anneal_counter = anneal_counter + 1;
                actual_lrate = actual_lrate0 / anneal_counter;
            else
                if G.learning.lrate0 > 0
                    actual_lrate = G.learning.lrate / (1 + G.iteration.n_updates / G.learning.lrate0);
                else
                    actual_lrate = G.learning.lrate;
                end
                actual_lrate0 = actual_lrate;
            end

            G.signals.lrates = [G.signals.lrates actual_lrate];

            % update
            for l = 1:n_layers
                biases_grad_old{l} = (1 - momentum) * biases_grad{l} + momentum * biases_grad_old{l};
                if l < n_layers
                    W_grad_old{l} = (1 - momentum) * W_grad{l} + momentum * W_grad_old{l};
                end
            end

            for l = 1:n_layers
                G.biases{l} = G.biases{l} - actual_lrate * (biases_grad_old{l}' + weight_decay * G.biases{l});
                if l < n_layers
                    G.W{l} = G.W{l} - actual_lrate * (W_grad_old{l} + weight_decay * G.W{l});
                end
            end
        end

        if G.verbose == 1
            fprintf(2, '.%f.', rerr);
        end

        if use_gpu > 0
            clear v0 h0 h0d h0e v0_clean vr hr deltae deltad v1clean
        end

        if early_stop
            n_valid = size(valid_patches, 1);
            rndidx = randperm(n_valid);
            v0valid = gpuArray(single(valid_patches(rndidx(1:round(n_valid * valid_portion)),:)));
            valid_sz = size(v0valid, 1);

            valid0 = cell(n_layers, 1);
            valid0{1} = v0valid;
            for l = 2:n_layers
                valid0{l} = zeros(valid_sz, G.structure.layers(l));
            end

            valid0 = gsn_sample(G, valid0, 1);

            rerr = mean(sum((v0valid - valid0{1}).^2, 1));

            if use_gpu > 0
                rerr = gather(rerr);
            end

            G.signals.valid_errors = [G.signals.valid_errors rerr];

            if valid_err == -Inf
                valid_err = rerr;
                valid_best_err = rerr;
                valid_violate_cnt = 0;

                M_best = G;
                M_best = pull_from_gpu (M_best);
            else
                prev_err = valid_err;
                valid_err = 0.99 * valid_err + 0.01 * rerr;

                if step > G.valid_min_epochs
                    if (1.1 * valid_best_err) < valid_err
                        fprintf(2, 'Early-stop! %f, %f\n', valid_err, valid_best_err);
                        stopping = 1;
                        break;
                    end

                    if valid_best_err <= valid_err
                        valid_violate_cnt = valid_violate_cnt + 1;
                        if valid_violate_cnt > (n_minibatches * G.valid_min_epochs)
                            fprintf(2, 'Unable to improve! %f, %f\n', valid_err, valid_best_err);
                            stopping = 1;
                            break;
                        end
                    else
                        valid_violate_cnt = 0;
                    end

                end

                if valid_err < valid_best_err
                    valid_best_err = valid_err;

                    G_best = G;
                    G_best = pull_from_gpu (G_best);
                end
            end
        else
            if G.stop.criterion > 0
                if G.stop.criterion == 1
                    if min_recon_error > G.signals.recon_errors(end)
                        min_recon_error = G.signals.recon_errors(end);
                        min_recon_error_update_idx = G.iteration.n_updates;
                    else
                        if G.iteration.n_updates > min_recon_error_update_idx + G.stop.recon_error.tolerate_count
                            fprintf(2, '\nStopping criterion reached (recon error) %f > %f\n', ...
                                G.signals.recon_errors(end), min_recon_error);
                            stopping = 1;
                            break;
                        end
                    end
                else
                    error ('Unknown stopping criterion %d', G.stop.criterion);
                end
            end
        end

        if length(G.hook.per_update) > 1
            err = G.hook.per_update{1}(G, G.hook.per_update{2});

            if err == -1
                stopping = 1;
                break;
            end
        end

    end

    if use_gpu > 0
        % pull
        G = pull_from_gpu (G);
    end

    if length(G.hook.per_epoch) > 1
        err = G.hook.per_epoch{1}(G, G.hook.per_epoch{2});

        if err == -1
            stopping = 1;
        end
    end

    if stopping == 1
        break;
    end

    if G.verbose == 1
        fprintf(2, '\n');
    end

    fprintf(2, 'Epoch %d/%d - recon_error: %f(%f) valid_error: %f\n', step, n_epochs, rerr_ma, rerr, valid_err);
end

if use_gpu > 0
    % pull
    G = pull_from_gpu (G);
end

if early_stop
    G = G_best;
end

end

function [G] = push_to_gpu (G)
    n_layers = length(G.structure.layers);

    % push
    for l = 1:n_layers
        if l < n_layers
            G.W{l} = gpuArray(single(G.W{l}));
        end
        G.biases{l} = gpuArray(single(G.biases{l}));
    end

    if G.adadelta.use
        for l = 1:n_layers
            if l < n_layers
                G.adadelta.gW{l} = gpuArray(single(G.adadelta.gW{l}));
                G.adadelta.W{l} = gpuArray(single(G.adadelta.W{l}));
            end
            G.adadelta.gbiases{l} = gpuArray(single(G.adadelta.gbiases{l}));
            G.adadelta.biases{l} = gpuArray(single(G.adadelta.biases{l}));
        end
    end
end

function [G] = pull_from_gpu (G)
    n_layers = length(G.structure.layers);

    for l = 1:n_layers
        if l < n_layers
            G.W{l} = gather(G.W{l});
        end
        G.biases{l} = gather(G.biases{l});
    end

    if G.adadelta.use
        for l = 1:n_layers
            if l < n_layers
                G.adadelta.W{l} = gather(G.adadelta.W{l});
                G.adadelta.gW{l} = gather(G.adadelta.gW{l});
            end
            G.adadelta.biases{l} = gather(G.adadelta.biases{l});
            G.adadelta.gbiases{l} = gather(G.adadelta.gbiases{l});
        end
    end

end

author: kyunghyuncho

repo: deepmat


text:
function ori = accumarray(subs,ori,varargin)
% accumarray for orientation
%
% Syntax
%   ori = accumarray(subs,ori)
%
% Input
%  subs -
%  ori - @orienrtation
%
% Output
%  ori - @orientation

% find a reference orientation for each class
ref = accumarray(subs,1:length(ori),[],@(x) x(1));
ori_ref = ori.subSet(ref(subs));

ori = project2FundamentalRegion(ori,ori_ref);

a = accumarray(subs,ori.a,varargin{:});
b = accumarray(subs,ori.b,varargin{:});
c = accumarray(subs,ori.c,varargin{:});
d = accumarray(subs,ori.d,varargin{:});

% normalize
s = sqrt(a.^2 + b.^2 + c.^2 + d.^2);
ori.a = a ./ s; ori.b = b ./ s; ori.c = c ./ s; ori.d = d ./ s;

ori.i = ori_ref.i(ref);

author: mtex-toolbox

repo: mtex

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


text:
classdef BaseApp < matlab.apps.AppBase & matlab.mixin.SetGetExactNames & ...
        wt.mixin.ErrorHandling
    % Base class for Widgets Toolbox apps

    % Copyright 2020-2021 The MathWorks, Inc.

    %% Properties
    properties (AbortSet)

        % Name of the app
        Name (1,1) string = "My App"

    end %properties

    properties (SetAccess = protected)

        % Model class for App preferences
        %(may subclass wt.model.Preferences to add more prefs)
        Preferences (1,1) wt.model.Preferences

        % Name of group to store preferences (defaults to class name)
        PreferenceGroup (1,1) string

    end %properties

    properties (AbortSet, Dependent)

        % Position of the app window
        Position

        % Visibility of the app window
        Visible

        % State of the app window
        WindowState

    end %properties

    %% Internal properties
    properties (Hidden, Transient, NonCopyable, SetAccess = immutable)

        % Figure window of the app
        Figure matlab.ui.Figure

        % Primary grid to place contents
        Grid matlab.ui.container.GridLayout

    end %properties

    properties (Transient, NonCopyable, ...
            Access = {?wt.apps.BaseApp, ?wt.test.BaseWidgetTest} )

        % Last used folder (for file operations)
        LastFolder (1,1) string = pwd

        % Is setup complete?
        SetupComplete (1,1) logical = false;

    end %properties

    %% Abstract methods (subclass must implement these)
    methods (Abstract, Access = protected)

        setup(app)

        update(app)

    end %methods

    %% Debugging Methods
    methods

        function forceUpdate(app)
            % Forces update to run (For debugging only!)

            app.update();

        end %function

    end %methods

    %% Constructor / destructor
    methods (Access = public)

        function app = BaseApp(varargin)
            % Constructor

            % Create the figure and hide until components are created
            app.Figure = uifigure( ...
                'AutoResizeChildren','off',...
                'Units','pixels', ...
                'DeleteFcn',@(h,e)delete(app), ...
                'CloseRequestFcn',@(h,e)close(app), ...
                'Visible','off');

            % Create MainGridLayout
            app.Grid = uigridlayout(app.Figure,[1 1]);
            app.Grid.Padding = [0 0 0 0];

            % Check for preference input and assign it first, in case
            % Preferences was subclassed
            % [splitArgs,varargin] = app.splitArgs('Preferences', varargin{:});
            % if ~isempty(splitArgs)
            %     app.Preferences = splitArgs{2};
            % end

            % Retrieve preferences
            app.loadPreferences();

            % Load last figure position
            % Note app.Figure.Position is inner position, where
            % app.Position is outer position. Inner position is the same
            % regardless of any menubar or toolbar, while outer position
            % changes if either is added. Use inner position for this
            % purpose so it does not depend on if/when any menubar or
            % toolbar is added to the figure.
            app.Figure.Position = app.getPreference('Position',[100 100 1000 700]);

            % Set up components
            app.setup_internal();
            app.setup();

            % Set any P-V pairs
            if ~isempty(varargin)
                set(app, varargin{:});
            end

            % Register the app with App Designer
            registerApp(app, app.Figure)

            % Ensure it's on screen
            app.moveOnScreen();

            % Mark the setup complete
            app.SetupComplete = true;

            % Update the app
            app.update();

            % Update the title
            app.updateTitle();

            % Force drawing to finish
            drawnow

            % Now, make it visible
            app.Figure.Visible = 'on';

        end %function

        function delete(app)
            % Destructor

            % Store last position in preferences
            if isscalar(app.Figure) && isvalid(app.Figure)
                app.setPreference('Position',app.Figure.Position)
            end

            % Save preferences
            app.savePreferences();

            % Now, delete the figure
            delete(app.Figure)

        end %function

    end %methods

    %% Public Methods
    methods

        function close(app)
            % Triggered on figure closed

            app.delete();

        end %function

        function selection = promptYesNoCancel(app, message, title, default, icon)
            % Prompt the user with a yes/no/cancel selection

            % Define arguments
            arguments
                app (1,1) wt.apps.BaseApp
                message (1,1) string = "Are you sure?"
                title (1,1) string = ""
                default (1,1) string = "Cancel"
                icon (1,1) string = "question"
            end

            % Launch the prompt
            selection = uiconfirm(app.Figure, message, title,...
                "Options",["Yes","No","Cancel"],...
                "DefaultOption",default,...
                "CancelOption","Cancel",...
                "Icon",icon);

        end %function

        function filePath = promptToSaveAs(app, filePath, filter, title)
            % Prompt the user to save a file

            % Define arguments
            arguments
                app (1,1) wt.apps.BaseApp
                filePath (1,1) string = pwd
                filter = ["*.mat","MATLAB MAT File"];
                title (1,1) string = "Save as"
            end

            % Prompt for the file
            [fileName,pathName] = uiputfile(filter, title, filePath);

            % Did the user cancel?
            if isequal(fileName,0)
                filePath = string.empty(0);
            else
                filePath = fullfile(pathName,fileName);
                app.LastFolder = pathName;
            end %if isequal(fileName,0)

        end %function

        function filePath = promptToLoad(app, filter, title)
            % Prompt the user to load a file

            % Define arguments
            arguments
                app (1,1) wt.apps.BaseApp
                filter = ["*.mat","MATLAB MAT File"];
                title (1,1) string = "Open"
            end

            % Prompt for the file
            [fileName,pathName] = uigetfile(filter, title, app.LastFolder);

            % Did the user cancel?
            if isequal(fileName,0)
                filePath = string.empty(0);
            else
                filePath = fullfile(pathName,fileName);
                app.LastFolder = pathName;
            end %if isequal(fileName,0)

        end %function

    end %methods

    %% Sealed Public methods
    methods (Sealed)

        function value = getPreference(app,propName,defaultValue)
            % Get an app preference from the Preferences object

            if isprop(app.Preferences, propName)
                value = app.Preferences.(propName);
            elseif nargin>2
                value = defaultValue;
            else
                value = [];
            end

        end %function

        function setPreference(app,propName,value)
            % Set an app preference in the Preferences object

            if ~isprop(app.Preferences,propName)
                addprop(app.Preferences,propName);
            end
            app.Preferences.(propName) = value;

        end %function

        function moveOnScreen(app)
            % Ensure the figure is placed on screen

            if strcmp(app.Figure.Units,'pixels')

                % Get the corners of each screen
                g = groot;
                screenPos = g.MonitorPositions;
                screenCornerA = screenPos(:,1:2);
                screenCornerB = screenPos(:,1:2) + screenPos(:,3:4) - 1;

                % Buffer for title bar
                titleBarHeight = 30;

                % Get the corners of the figure (bottom left and top right)
                figPos = app.Figure.OuterPosition;
                figCornerA = figPos(1:2);
                figCornerB = figPos(1:2) + figPos(:,3:4) - 1;

                % Are the corners on any screen?
                aIsOnScreen = all( figCornerA >= screenCornerA & ...
                    figCornerA <= screenCornerB, 2 );
                bIsOnScreen = all( figCornerB >= screenCornerA & ...
                    figCornerB <= screenCornerB, 2);

                % Are corners on a screen?

                % Are both corners fully on any screen?
                if any(aIsOnScreen) && any(bIsOnScreen)
                    % Yes - do nothing

                elseif any(bIsOnScreen)
                    % No - only upper right corner is on a screen

                    % Calculate the adjustment needed, and make it
                    figAdjust = max(figCornerA, screenCornerA(bIsOnScreen,:)) ...
                        - figCornerA;
                    figPos(1:2) = figPos(1:2) + figAdjust;

                    % Ensure the upper right corner still fits
                    figPos(3:4) = min(figPos(3:4), ...
                        screenCornerB(bIsOnScreen,:) - figPos(1:2) - [0 titleBarHeight] + 1);

                    % Move the figure
                    app.Figure.Position = figPos;

                elseif any(aIsOnScreen)
                    % No - only lower left corner is on a screen

                    % Calculate the adjustment needed, and make it
                    figAdjust = min(figCornerB, screenCornerB(aIsOnScreen,:)) ...
                        - figCornerB;
                    figPos(1:2) = max( screenCornerA(aIsOnScreen,:),...
                        figPos(1:2) + figAdjust );

                    % Ensure the upper right corner still fits
                    figPos(3:4) = min(figPos(3:4), ...
                        screenCornerB(aIsOnScreen,:) - figPos(1:2) - [0 titleBarHeight] + 1);

                    % Move the figure
                    app.Figure.Position = figPos;

                else
                    % No - Not on any screen

                    % This is slower, but uncommon anyway
                    movegui(app.Figure,'onscreen');

                end %if any( all(aIsOnScreen,2) & all(bIsOnScreen,2) )

            else

                % This is slower, but uncommon anyway
                movegui(app.Figure,'onscreen');

            end %if strcmp(app.Figure.Units,'pixels')

        end %function

    end %methods

    %% Protected Methods
    methods (Access = protected)

        function setup_internal(~)
            % Preform internal pre-setup necessary

            % This is used for session managed apps

        end %function

        function loadPreferences(app)
            % Load stored preferences

            app.Preferences.load(app.PreferenceGroup);

        end %function

        function savePreferences(app)
            % Save preferences

            app.Preferences.save(app.PreferenceGroup);

        end %function

        function updateTitle(app)
            % Update the figure title

            app.Figure.Name = app.Name;

        end %function

    end %methods

    %% Accessors
    methods

        function set.Name(app,value)
            app.Name = value;
            app.updateTitle();
        end

        function value = get.PreferenceGroup(app)
            value = app.PreferenceGroup;
            if isempty(value)
                value = class(app);
            end
            value = matlab.lang.makeValidName(value);
        end

        function value = get.Position(app)
            value = app.Figure.Position;
        end

        function set.Position(app,value)
            app.Figure.Position = value;
        end

        function value = get.Visible(app)
            value = app.Figure.Visible;
        end

        function set.Visible(app,value)
            app.Figure.Visible = value;
        end

        function value = get.WindowState(app)
            value = app.Figure.WindowState;
        end

        function set.WindowState(app,value)
            app.Figure.WindowState = value;
        end

    end %methods

end % classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ r, seed ] = r8col_uniform_abvec ( m, n, a, b, seed )

%*****************************************************************************80
%
%% R8COL_UNIFORM_ABVEC fills an R8COL with scaled pseudorandom numbers.
%
%  Discussion:
%
%    An R8COL is an array of R8 values, regarded as a set of column vectors.
%
%    The user specifies a minimum and maximum value for each row.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in
%    the array.
%
%    Input, real A(M), B(M), the lower and upper limits.
%
%    Input/output, integer SEED, the "seed" value, which
%    should NOT be 0.  On output, SEED has been updated.
%
%    Output, real R(M,N), the array of pseudorandom values.
%
  i4_huge = 2147483647;

  for j = 1 : n

    for i = 1 : m

      k = floor ( seed / 127773 );

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

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

      r(i,j) = a(i) + ( b(i) - a(i) ) * seed * 4.656612875E-10;

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   MoveAbsJ:
%   Moves robot to a particular joint coordinates
%   Example: MoveAbsJ(robot, joint_coord, 'v1000', 'z100', gripper, wobj0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

global configuration robot

fprintf('\nCall to MoveAbsJ');

%obtain current joint coordinates
q_current=robot.q;

%obtain target joint coordinates
q_final= joint_coord(1:robot.DOF)';

%test joint limits for the final values
test_joints(robot, q_final);

speed = obtain_joint_speed(robot, speeddata);
%when fine is specified, the radius is zero, otherwise,
%the radius represents a sphere around the target where the
%movement is changed to the next target
radius = obtain_zone_data(zonedata);

%find time to perform movement
[actual_speed, tmax]=synchronize(q_current, q_final, speed, robot.accelmax);

%in this case, if the total time for the movement is lower than the
%configuration.delta_time time then 10 points are interpolated
if tmax <= configuration.delta_time
    t = [0:tmax/10:tmax]';
else
    %local time for the planning, normal case
    t = [0:configuration.delta_time:tmax]';
end

if radius==0
    %compute a smooth trajectory in q and qd
    %in this case, the final speed is zero
    [q, qd, qdd] = compute_joint_trajectory_indep(robot, q_current, q_final, robot.qd, zeros(1,robot.DOF), t);
else
    %in this case, the final speed is selected as the speed selected by the
    %variable speeddata
    [q, qd, qdd] = compute_joint_trajectory_indep(robot, q_current, q_final, robot.qd, actual_speed, t);
end

T=directkinematic(robot, q_final);
%finds the first joint coordinates that are inside a radius r of the target
%point
index = find_first_in_zone_data(robot, q, T, radius);

%the robot performs the movement until the index found. The coordinates, joint speed and acceleratin
%are stored and used in the planning of the next point
robot.q=q(:,index);
robot.qd=qd(:,index);
robot.qdd=qdd(:,index);

%store all the trajectory for plotting
%the joint trajectories, speeds and acceleration of susequent movements are
%store here
robot.q_vector=[q(:, 1:index)];
robot.qd_vector=[qd(:, 1:index)];
robot.qdd_vector=[qdd(:, 1:index)];

%a global time for the planning is computed.
%in this way, the total trajectory of different movements can be plotted
if length(robot.time)==0
    tend = 0;
else
    tend = robot.time(end);
end
t = t + tend;
%store total time
%robot.time=[robot.time t(1:index)'];
robot.time=[t(1:index)'];

%Test whether there are joints outside mechanical limits
test_joint_limits(robot);

%Plot position, velocity and acceleration
%plot_joint_data(robot);
%Now, animate the robot in 3D
animate(robot, q);

%Plot the trajectory
pp = zeros(3,size(q,2));
for i=1:size(q,2),
    T=directkinematic(robot, q(:,i));
    pp(1:3,i)=T(1:3,4);
end
plot3(pp(1,:),pp(2,:),pp(3,:), 'k', 'LineWidth', 2)

repo: ARTE


text:
function [ALLEEG_out cfg] = pop_pre_prepData(ALLEEG,typeproc,varargin)
%
% Preprocess EEG dataset(s) for connectivity analysis. See [1] for
% mathematical details on preprocessing steps.
%
%
% Input:
%
%   ALLEEG:         Array of EEGLAB datasets to preprocess.
%   typeproc:       Reserved for future use. Use 0
%
% Optional:
%
%   <'Name',value> pairs as defined in pre_prepData()
%
% Output:
%
%   ALLEEG:         Prepocessed EEG structure(s)
%   cfg:            Argument specification structure.
%
%
% See Also: pre_prepData()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Section 6.5.1
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% Author: Tim Mullen 2009, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu
%
% Revised Jan 2010.

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

if nargin<2
    typeproc = 0;
end

% set default output
ALLEEG_out = ALLEEG;
cfg = [];

% generate splash screen
% initialize SIFT, etc
StartSIFT(~strcmpi(typeproc,'nogui'),true);

fcnName     = strrep(mfilename,'pop_','');
fcnHandle   = str2func(fcnName);

% check if we've applied SIFT to this dataset before
res = hlp_checkeegset(ALLEEG,{'cat'});
if isempty(res) && isfield(ALLEEG(1).CAT.configs,fcnName)
    % get default configuration (from prior use) and merge with varargin
    varargin = [hlp_struct2varargin(ALLEEG(1).CAT.configs.(fcnName)) varargin];
end

if strcmpi(typeproc,'nogui')
    % get the default config from function and overload supplied args
    cfg = arg_tovals(arg_report('rich',fcnHandle,[{'EEG',ALLEEG(1)},varargin]),false);
else
    % render the GUI
    [PGh figh] = feval(['gui_' fcnName],ALLEEG(1),varargin{:});

    if isempty(PGh)
        % user chose to cancel
        return;
    end

    % get the specification of the PropertyGrid
    ps = PGh.GetPropertySpecification;
    cfg = arg_tovals(ps,false);
end

drawnow;

if strcmpi(typeproc,'cfg_only')
    return;
end

% initialize progress bar
if cfg.verb==2 && length(ALLEEG)>1
    waitbarTitle = 'Preprocessing datasets';

    multiWaitbar(waitbarTitle,'Reset');
    multiWaitbar(waitbarTitle,'ResetCancel',true);
    multiWaitbar(waitbarTitle,...
                 'Color', [0.8 0.0 0.1],  ...
                 'CanCancel','on',        ...
                 'CancelFcn',@(a,b)disp('[Cancel requested. Please wait...]'));
end

% re-initialize output
clear ALLEEG_out;

% preprocess datasets
for cnd=1:length(ALLEEG)

    % execute the low-level function
    [ALLEEG_out(cnd)] = feval(fcnHandle,'EEG',ALLEEG(cnd),cfg);

    if ~isempty(cfg)
        % store the configuration structure
        ALLEEG_out(cnd).CAT.configs.(fcnName) = cfg;
    end

    if cfg.verb==2 && length(ALLEEG)>1
        % update waitbar
        drawnow;
        cancel = multiWaitbar(waitbarTitle,cnd/length(ALLEEG));
        if cancel && hlp_confirmWaitbarCancel(waitbarTitle)
            % restore original dataset
            ALLEEG_out = ALLEEG;
            break;
        end
    end

end

% cleanup progress bar
if cfg.verb==2 && length(ALLEEG)>1
    multiWaitbar(waitbarTitle,'Close');
end

repo: REST


text:
function [counter, listofimages] = HandleImgFolder(pth)
%
% read all images - and find problematic ones
%
% Usage:
% counter = HandleImgFolder( pth )
%
% Inputs:
% pth - full path on local disk where LM (or SUN) images are.
% this should be the same path as "destImages" used in SUNinstall function
%
% Outputs:
% counter - number of problematic images encountered and deleted
%
% Note: problematic images will be deleted, their names (full path) will be
% printed to stdout.
%
% written by: Shai Bagon, Oct 2012.
%
% modified by Torralba to return a list of problematic images
%

counter = 0;
listofimages  = {};

fls = dir(fullfile(pth,'*'));

imgfrmt = imformats();

for fi=3:numel(fls) % skip first '.' and '..'
    fnm = fullfile(pth, fls(fi).name);
    % process sub folder
    if fls(fi).isdir
        disp(fls(fi))
        [c,l] = HandleImgFolder( fnm );
        counter = counter + c;
        listofimages = [listofimages l];
        continue;
    end

    % process image file
    try
        fbit = arrayfun( @(x) x.isa(fnm), imgfrmt);
        if any( fbit )
            % switch off warnings
            lastwarn('');
            imgfrmt(fbit).read( fnm );
            [wm wi] = lastwarn();
            if ~isempty(wi)
                lastwarn(''); % clear warning
                error('LMVerifyDownload:HandleFolder:warn',...
                    'Got warning %s\nthrowing a custom error', wm);
            end
        end
    catch em
        fprintf(1, '\tProblematic image %s \n', fnm);
        %delete(fnm);
        counter = counter + 1;
        listofimages{counter} = fnm;
    end
end

author: CSAILVision

repo: LabelMeToolbox


text:
function[s,t,E] = sta(data_spk,data_lfp,smp,plt,w,T,D,err)
% Spike Triggered Average
%     Usage: [s,t,E] = sta(data_spk,data_lfp,smp,plt,w,T,D,err)
%
% Inputs
%
% Note that all times have to be consistent. If data_spk
% is in seconds, so must be sig and t. If data_spk is in
% samples, so must sig and t. The default is seconds.
%
% data_spk    - strucuture array of spike times data
%               or NaN padded matrix
% data_lfp    - array of lfp data(samples x trials)
%
% Optional...
% plt 'n'|'r' etc
% width kernel smoothing in s
% T = [-0.1 0.1] - extract this range about each spk
% D = plot spike triggered average out to [D1 D2]
% err = calcluate error bars (bootstrap)
%
% Outputs:
%
% s  spike triggered average
% t  times
% E  bootstrap standard err

if nargin < 3;error('Require spike, lfp and lfp times');end
if isstruct(data_spk)
   [data_spk]=padNaN(data_spk); % create a zero padded data matrix from input structural array
   sz=size(data_spk);
   if sz(1)>sz(2); data_spk=data_spk'; end;% transpose data to get in form compatible with Murray's routine
else
   sz=size(data_spk);
   if sz(1)>sz(2); data_spk=data_spk'; end;% transpose data to get in form compatible with Murray's routine
end
sz=size(data_lfp);
if sz(1)>sz(2); data_lfp=data_lfp'; end;% transpose data to get in form compatible with Murray's routine
verbose = 1;
t = smp;
if nargin < 4; plt = 'r'; end
if nargin < 5; w = 0; end
if nargin < 6; T = [min(t) max(t)]; end
if nargin < 7; D = 0.25*[-1 1]; end
if nargin < 8; err = 1;end

if isempty(plt); plt = 'r'; end
if isempty(w); w = 0; end
if isempty(T); T = [min(t) max(t)]; end
if isempty(D); D = 0.25*[-1 1]; end
if isempty(err); err = 1;end

if w > (T(2)-T(1))/2
  disp('Smoothing > data segment : should be in seconds : turn off smoothing')
  w = 0;
end

sz = size(data_spk);
NT = sz(1);
mlfp = 0;
slfp = 0;
Nspk = 0;
smp = t(2)-t(1);
if D(1) <= 0 && D(2) >= 0
  t1 = [D(1):smp:(-smp+eps) 0:smp:D(2)+eps];
else
  t1 = (round(D(1)/smp)*smp):smp:D(2);
end

% count up the spikes...

if err
  for n=1:NT
    indx = find(t>T(1)&t<T(2));
%     lfp = data_lfp(n,indx);
    spk = data_spk(n,data_spk(n,:)>T(1) && data_spk(n,:)<T(2) && data_spk(n,:)~=0);
    tt = t(indx);
    if ~isempty(spk) > 0
      ND = length(spk);
      for s=1:ND
        spktime = spk(s);
        t0 = tt-spktime + eps;
        if min(t0)>(D(1)-smp) || max(t0)<(D(2)+smp); break;end
        Nspk = Nspk + 1;
      end
    end
  end
  Err = zeros(Nspk,length(t1));
  Nspk = 0;
end

for n=1:NT
  indx = find(t>T(1)&t<T(2));
  lfp = data_lfp(n,indx);
  spk = data_spk(n,data_spk(n,:)>T(1) && data_spk(n,:)<T(2) && data_spk(n,:)~=0);
  tt = t(indx);
  if ~isempty(spk) > 0
    ND = length(spk);
    for s=1:ND
      spktime = spk(s);
      t0 = tt-spktime + eps;
      if min(t0) < (D(1)-smp) && max(t0) > (D(2)+smp);
        indx = find(t0<D(1));
        indx = indx(length(indx));
        offset = (t0(indx)-D(1))/smp;
        indx = indx:(indx+length(t1)-1);
        lfp_t1 = lfp(indx) + (lfp(indx+1)-lfp(indx))*offset;
        Nspk = Nspk + 1;
        mlfp = mlfp + lfp_t1;
        slfp = slfp + lfp_t1.^2;
        Err(Nspk,:) = lfp_t1;
      end
    end
  end
end
if Nspk == 0
  if verbose;disp('No spikes in interval');end
  t = t1;
  s = zeros(length(t),1);
  E = zeros(length(t),1);
  return
end
mlfp = mlfp/Nspk;
slfp = slfp/Nspk;
stdlfp = sqrt((slfp - mlfp.^2)/Nspk);

% local smoother...

N = fix(w/smp);
if N > 5
  mlfp = locsmooth(mlfp,N,fix(N/2));
end

% bootstrap errorbars...

if err == 1;
  Nboot = 20;
  blfp = 0;
  slfp = 0;
  for n = 1:Nboot
    indx = floor(Nspk*rand(1,Nspk)) + 1;
    lfptmp = mean(Err(indx,:));
    if N > 5
      lfptmp = locsmooth(lfptmp,N,fix(N/2));
    end
    blfp = blfp + lfptmp;
    slfp = slfp + lfptmp.^2;
  end
  stdlfp = sqrt((slfp/Nboot - blfp.^2/Nboot^2));
end

s = mlfp-mean(mlfp);
E = stdlfp;
t = t1;

%cols = 'krbgycm';
if plt == 'n';return;end
plot(t1,s,plt)
xax = get(gca,'xlim');
%yax = get(gca,'ylim');
if err == 1
  me = real(2*mean(stdlfp));
  line(xax,me*[1 1],'color','b')
  line(xax,-me*[1 1],'color','b')
  line(xax,0*[1 1],'color','k')
%  errorbar(0.1*xax(2)+0.9*xax(1),0.1*yax(2)+0.9*yax(1), ...
%         mean(stdlfp),'k')
%plot(0.1*xax(2)+0.9*xax(1),0.1*yax(2)+0.9*yax(1),'k.')
end

title(['spike triggered average : ' ...
      num2str(Nspk) ' used : '     ...
      ' Errorbars are two std err']);
%line(get(gca,'xlim'),mean(mlfp)*[1 1],'color','k')
line([0 0],get(gca,'ylim'),'color','k')
hold off

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function  sL = subBoundaryLength(grains,varargin)
% subgrain boundary length per grain
%
% Input
%  grains - @grain2d
%
% Output
%  l - number of boundary segments
%
% Syntax
%   l = grains.subBoundaryLength
%

grainIds = grains.innerBoundary.grainId;

sL = grains.innerBoundary.segLength;

% if varargin is a logical use for selection of relevant innerBoundary
if nargin>1 && islogical(varargin{1})
    grainIds = grainIds(varargin{1},:);
    sL = sL(varargin{1});
end

ind = (diff(grainIds,1,2)~=0);
grainIds(ind,:) = [];
sL(ind,:) = [];

sL = accumarray(grainIds(:,1),sL,[max(grains.id) 1]);

% convert from id to ind
sL = sL(grains.id);


text:
function exportImg(file_name, data_map, c_map)
% Save on a file a matrix provided in "data"
%
% SYNTAX:
%   exportImg(file_name, data_map, <col_map>)

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

    narginchk(2,3);

    if islogical(data_map)
        imwrite(data_map, file_name);
    else
        if nargin < 3
            c_map = jet(256);
        end

        % Normalize data
        im_alpha = isnan(data_map);
        data_map = data_map - min(data_map(:));
        data_map = max(2, ceil(data_map ./ max(data_map(:)) .* (size(c_map, 1) -1)) + 1);
        data_map(im_alpha) = 1;
        c_map = [[0, 0, 0]; c_map];

        im = reshape([c_map(data_map(:), 1); c_map(data_map(:), 2); c_map(data_map(:), 3)], size(data_map, 1), size(data_map, 2), 3);
        imwrite(im, file_name, 'Alpha', 1-im_alpha);
    end
end

author: goGPS-Project


text:
classdef comparison

methods(Static)

        function result = matching_atoms_ratio(original, new, options)
            % Returns distance betweeen two dictionaries
            % in terms of how many atoms from the original
            % dictionary could be
            % found in the new dictionary
            % Assumes that all atoms in dictionary are normalized.
            % Works well if the dictionaries have low coherence.
            distance_threshold=0.01;
            if nargin > 2
                if isfield(options, 'distance_threshold')
                    distance_threshold = options.distance_threshold;
                end
            end
            % Number of atoms which could be
            % identified properly
            atomsFound=0;
            % We compute inner product of each
            % atom in one dictionary with all
            % atoms in other dictionary
            innerProducts =abs(original'*new);
            % Find the number of atoms in original dictionary
            numAtoms = size(original,2);
            for i=1:1:numAtoms
                %  Trying to find a match for the i-th atom in original dictionary.
                % We find its highest inner product
                % with atoms in new dictionary
                max_similarity = max(innerProducts(i,:));
                % The distance is 1  - max_similarity
                min_distance =1- max_similarity;
                % If distance is less than the threshold
                % then the atom has been found.
                match_found = (min_distance<distance_threshold);
                atomsFound=atomsFound + match_found;
            end;
            % Find the ratio of found atoms with total atoms
            result= atomsFound/numAtoms;
        end

end

end

author: indigits

repo: sparse-plex


text:
function [Y, Phi] = modelOut(model, X, varargin)

% MODELOUT Give the output of a model for given X.
% FORMAT
% DESC gives the output of the model for a given input X. For
% latent variable models it gives a position in data space given a
% position in latent space.
% ARG model : structure specifying the model.
% ARG X : input location(s) for which output is to be computed.
% RETURN Y : output location(s) corresponding to given input
% locations.
%
% FORMAT
% DESC gives the output of the model for a given input X. For
% latent variable models it gives a position in data space given a
% position in latent space.
% ARG model : structure specifying the model.
% ARG X : input location(s) for which output is to be computed.
% RETURN Phi : output basis function(s) corresponding to given input
% RETURN Y : output location(s) corresponding to given input
% locations.
%
% SEEALSO : modelCreate
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006
%
% MODIFICATIONS : Cark Henrik Ek, 2008

% MLTOOLS

fhandle = str2func([model.type 'Out']);
if nargout > 1
  [Y, Phi] = fhandle(model, X, varargin{:});
else
  Y = fhandle(model, X, varargin{:});
end
if(isfield(model,'indexOut')&&~isempty(model.indexOut))
  Y(:,setdiff(1:1:size(Y,2),model.indexOut)) = NaN;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function e = R2e(R)

% R2E  Rotation matrix to Euler angles conversion.
%   R2E(R) returns the euler angles vector [roll;pitch;yaw] corresponding
%   to the orientation of the rotation matrix R. The result is such that
%   E2R(R2E(R)) = R and R2E(E2R(E)) = E.
%
%   See also E2R, FRAME.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

s = whos('R');

if (strcmp(s.class,'sym'))
    roll  = atan(R(3,2)/R(3,3));
    pitch = asin(-R(3,1));
    yaw   = atan(R(2,1)/R(1,1));
else
    roll  = atan2(R(3,2),R(3,3));
    pitch = atan2(-R(3,1), sqrt(R(1,1)^2+R(2,1)^2));
    yaw   = atan2(R(2,1),R(1,1));
end

e = [roll;pitch;yaw];

% ========== 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:
classdef prtClassMaryLogDiscSparseKernel < prtClassMaryLogDisc

    properties (SetAccess=private)
    end

    properties (SetAccess = protected)
        % Lambda
        lambda = 10;  %2; for dataBiModal;
        k
    end

    methods

        function self = prtClassMaryLogDiscSparseKernel(varargin)

            self = prtUtilAssignStringValuePairs(self,varargin{:});
            self.wChangeTolerance = 1;
        end

    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,dataSet)
            %self = trainAction(self,dataSet)

            x = dataSet.getObservations;
            %x = cat(2,ones(size(x,1),1),x); %DC component
            k = prtKernelRbf;
            k = k.train(dataSet);
            self.k = k;

            x = getObservations(k.run(dataSet));
            y = dataSet.getTargetsAsBinaryMatrix;

            nClasses = dataSet.nClasses;
            d = size(x,2);

            %random initialization; last set of weights is set to 0
            numWeights = d*(nClasses-1);
            weightMatrix = randn(nClasses-1,d);
            weightMatrix = cat(1,weightMatrix,zeros(1,size(weightMatrix,2)));
            weightMatrixOld = weightMatrix;

            %Can calculate B matrix outside loop; makes life fast
            xx = x'*x;
            B = kron(-1/2*(eye(nClasses-1)-ones(nClasses-1)/dataSet.nClasses),xx);

            self.converged = false;

            for j = 1:self.maxIter
                for k = 1:numWeights
                    psi = (weightMatrix*x')';
                    py = bsxfun(@rdivide,exp(psi),sum(exp(psi),2));

                    %label error
                    yError = y-py;

                    if j == 1
                        kronIndex1 = repmat(1:size(yError,2),size(x,2),1);
                        kronIndex2 = repmat(1:size(x,2),1,size(y,2));
                        xRepmat = x(:,kronIndex2);
                    end
                    %Can we speed this up?  KRON is needlessly slow
                    %                     g = 0;
                    %                     for i = 1:size(yError,1)
                    %                         g = g + kron(yError(i,:),x(i,:));
                    %                     end
                    %                     g = g(:);
                    %                     tic;
                    g = sum(yError(:,kronIndex1(:)).*xRepmat);
                    g = g(:);

                    wVec = weightMatrix(1:end-1,:)';
                    wVec = wVec(:);
                    wVec(k) = prtUtilSoft(wVec(k)-g(k)/B(k,k),-self.lambda./B(k,k));
                    weightMatrix(1:end-1,:) = reshape(wVec,size(weightMatrix(1:end-1,:),2),size(weightMatrix(1:end-1,:),1))';
                end
                stem(weightMatrix(1:end-1,:)');
                title(j);
                drawnow;
                if norm(weightMatrix(:)-weightMatrixOld(:)) < self.wChangeTolerance
                    self.converged = true;
                    break;
                end
                weightMatrixOld = weightMatrix;

            end
            self.wMat = weightMatrix;
        end

        function ClassifierResults = runAction(self,dataSet)
            %ClassifierResults = runAction(self,DataSet)

            x = getObservations(self.k.run(dataSet));

            psi = (self.wMat*x')';
            y = bsxfun(@rdivide,exp(psi),sum(exp(psi),2));

            ClassifierResults = dataSet.setObservations(y);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [D]=d2a1(A,B)
%D2A1 Pairwise squared L2 distance along 1st axis.
% D = D2A1(A,B) Computes squared L2 (Euclidean) distance
% between all pairs of d-dimensional points in A and B.
%
% Inputs:
% A - d-by-m matrix of m d-dimensional points;
% B - d-by-n matrix of n d-dimensional points.
%
% Outputs:
% D - m-by-n matrix of pairwise distances.
%
% See also D2A2.

if nargin < 2
    D = full(A'*A);
    d = diag(D);
    D = bsxfun(@minus, d, 2*D);
    D = bsxfun(@plus, d', D);
    D = max(D, 0);
else
    D = full(A'*B);
    D = bsxfun(@minus, sum(B.^2,1), 2*D);
    D = bsxfun(@plus, sum(A.^2,1)', D);
    D = max(D, 0);
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function seq=seqa(a,b,c);
% PURPOSE: produce a sequence of values
% -----------------------------------------------------
% USAGE: y = seqa(a,b,c)
%  where    a = initial value in sequence
%           b = increment
%           c = number of values in the sequence
% -----------------------------------------------------
% RETURNS: a sequence, (a:b:(a+b*(c-1)))' in MATLAB notation
% -----------------------------------------------------
% NOTE: a Gauss compatability function
% -----------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jpl@jpl.econ.utoledo.edu

% seqa Gauss eqivalent of seqa(a,b,c)
seq=(a:b:(a+b*(c-1)))';
return;

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function OutputFiles = db_set_noisecov(iSrcStudy, Target, isDataCov, ReplaceFile)
% DB_SET_NOISECOV: Apply a noise covariance node to other studies
%
% USAGE:  OutputFiles = db_set_noisecov(iSrcStudy, iDestStudies, isDataCov=0, ReplaceFile=[ask])  : Apply to the target studies
%         OutputFiles = db_set_noisecov(iSrcStudy, 'AllConditions')                               : Apply to all the conditons in the same subject
%         OutputFiles = db_set_noisecov(iSrcStudy, 'AllSubjects')                                 : Apply to all the conditons in all the subjects

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009-2016

% Parse inputs
if (nargin < 4) || isempty(ReplaceFile)
    ReplaceFile = [];
end
if (nargin < 3) || isempty(isDataCov)
    isDataCov = 0;
end
OutputFiles = {};

% ===== GET SOURCE STUDY =====
% Get source study
sSrcStudy = bst_get('Study', iSrcStudy);
% Check for noise covariance matrix
if isempty(sSrcStudy.Channel)
    error('This operation requires that you have a channel file defined in the source and destination studies.');
end
if isempty(sSrcStudy.NoiseCov) || (~isDataCov && isempty(sSrcStudy.NoiseCov(1).FileName)) || (isDataCov && ((length(sSrcStudy.NoiseCov) < 2) || isempty(sSrcStudy.NoiseCov(2).FileName)))
    error('No noise covariance available in the source study.');
end

% ===== GET TARGET STUDIES =====
if isnumeric(Target)
    % Destination studies are passed in argument
    iDestStudies = Target;
elseif strcmpi(Target, 'AllConditions')
    % Get all the studies for this subject
    [sDestStudies, iDestStudies] = bst_get('StudyWithSubject', sSrcStudy.BrainStormSubject);
elseif strcmpi(Target, 'AllSubjects')
    % Get the whole database
    ProtocolSubjects = bst_get('ProtocolSubjects');
    % Get list of subjects (sorted alphabetically => same order as in the tree)
    [uniqueSubjects, iUniqueSubjects] = sort({ProtocolSubjects.Subject.Name});
    % Process each subject
    iDestStudies = [];
    for iSubj = 1:length(uniqueSubjects)
        % Get subject filename
        iSubject = iUniqueSubjects(iSubj);
        SubjectFile = ProtocolSubjects.Subject(iSubject).FileName;
        % Get all the studies for this subject
        [sStudies, iStudies] = bst_get('StudyWithSubject', SubjectFile, 'intra_subject', 'default_study');
        iDestStudies = [iDestStudies, iStudies];
    end
else
    return;
end
% Get the channels studies for those studies
[sDestChannels, iDestChanStudies] = bst_get('ChannelForStudy', iDestStudies);
% Unique list
iDestChanStudies = unique(iDestChanStudies);
% Remove the one the NoiseCov file comes from
iDestChanStudies = setdiff(iDestChanStudies, iSrcStudy);
% If there are no targets
if isempty(iDestChanStudies)
    error('No destination studies.');
end

% ===== READ SOURCE DATA =====
% Load noisecov file
if isDataCov
    SrcNoiseCovMat = load(file_fullpath(sSrcStudy.NoiseCov(2).FileName));
else
    SrcNoiseCovMat = load(file_fullpath(sSrcStudy.NoiseCov(1).FileName));
end
% Load channel file
SrcChannelMat = in_bst_channel(sSrcStudy.Channel.FileName);

% ===== APPLY NOISECOV =====
% Process each target study
for i = 1:length(iDestChanStudies)
    % Get destination study
    iDestStudy = iDestChanStudies(i);
    sDestStudy = bst_get('Study', iDestStudy);
    % If no channel defined: ignore destination channel
    if isempty(sDestStudy.Channel) || isempty(sDestStudy.Channel.FileName)
        disp(['BST> Warning: No channel file in ' bst_fileparts(sDestStudy.FileName) ': Ignoring condition...']);
        continue;
    end
    % Load destination channel file
    DestChannelMat = in_bst_channel(sDestStudy.Channel.FileName);
    % Intialize destination noise covariance matrix
    nbDestChan = length(DestChannelMat.Channel);
    DestNoiseCovMat = SrcNoiseCovMat;
    DestNoiseCovMat.NoiseCov = zeros(nbDestChan);

    % Find the indices of src channels in dest channels
    iSrcChan = [];
    iDestChan = [];
    for iChan = 1:nbDestChan
        iTmp = find(strcmpi({SrcChannelMat.Channel.Name}, DestChannelMat.Channel(iChan).Name));
        if ~isempty(iTmp) && ~isempty(DestChannelMat.Channel(iChan).Name)
            iDestChan(end+1) = iChan;
            iSrcChan(end+1) = iTmp(1);
        end
    end

    % If some channels were the same in both files
    if ~isempty(iDestChan)
        DestNoiseCovMat.NoiseCov(iDestChan,iDestChan) = SrcNoiseCovMat.NoiseCov(iSrcChan,iSrcChan);
        % Other fields
        if isfield(SrcNoiseCovMat, 'FourthMoment') && ~isempty(SrcNoiseCovMat.FourthMoment)
            DestNoiseCovMat.FourthMoment = zeros(nbDestChan);
            DestNoiseCovMat.FourthMoment(iDestChan,iDestChan) = SrcNoiseCovMat.FourthMoment(iSrcChan,iSrcChan);
        end
        if isfield(SrcNoiseCovMat, 'nSamples') && ~isempty(SrcNoiseCovMat.nSamples)
            DestNoiseCovMat.nSamples = zeros(nbDestChan);
            DestNoiseCovMat.nSamples(iDestChan,iDestChan) = SrcNoiseCovMat.nSamples(iSrcChan,iSrcChan);
        end
        % Copy this structure in all the target studies
        NewFile = import_noisecov(iDestStudy, DestNoiseCovMat, ReplaceFile, isDataCov);
        if ~isempty(NewFile) && ischar(NewFile)
            OutputFiles{end+1} = NewFile;
        end
    end
end

author: brainstorm-tools


text:
%compute tail speed in the direction perpendicular to tail-head direction
function [data,units]=compute_veltailperth(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
veltailperth=cell(1,numlarvae);
%tailheadangperp=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    %tailheadangperp{1,i}=trx(larva).tailheadang-pi/2;
    veltailperth{1,i}=trx(larva).velmagtail.*(cos(trx(larva).velangtail).*cos(trx(larva).tailheadang(1,1:end-1)+pi/2)+sin(trx(larva).velangtail).*sin(trx(larva).tailheadang(1,1:end-1)+pi/2));
end

units=parseunits('mm/s');
data=veltailperth;

author: kristinbranson

repo: JAABA


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Demonstration of Basic Utilities of JSONlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rngstate = rand ('state');
randseed=hex2dec('623F9A9E');
clear data2json json2data

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a simple scalar value \n')
fprintf(1,'%%=================================================\n\n')

data2json=pi
savebj('',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty array \n')
fprintf(1,'%%=================================================\n\n')

data2json=[]
savebj('empty',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an ampty string \n')
fprintf(1,'%%=================================================\n\n')

data2json=''
savebj('emptystr',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a simple row vector \n')
fprintf(1,'%%=================================================\n\n')

data2json=1:3
savebj('',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a simple column vector \n')
fprintf(1,'%%=================================================\n\n')

data2json=(1:3)'
savebj('',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a string array \n')
fprintf(1,'%%=================================================\n\n')

data2json=['AC';'EG']
savebj('',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a string with escape symbols \n')
fprintf(1,'%%=================================================\n\n')

data2json=sprintf('AB\tCD\none"two')
savebj('str',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a mix-typed cell \n')
fprintf(1,'%%=================================================\n\n')

data2json={'a',true,[2;3]}
savebj('',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 3-D array in nested array form\n')
fprintf(1,'%%=================================================\n\n')

data2json=reshape(1:(2*4*6),[2,4,6]);
savebj('',data2json,'NestArray',1)
json2data=loadbj(ans)
% if(any(json2data(:)~=data2json(:)) || any(size(json2data)~=size(data2json)))
%     warning('conversion does not preserve original data');
% end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 3-D array in annotated array form\n')
fprintf(1,'%%=================================================\n\n')

data2json=reshape(1:(2*4*6),[2,4,6]);
savebj('',data2json,'NestArray',0)
json2data=loadbj(ans)
if(any(json2data(:)~=data2json(:)) || any(size(json2data)~=size(data2json)))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 4-D array in annotated array form\n')
fprintf(1,'%%=================================================\n\n')

data2json=reshape(1:(2*4*3*2),[2,4,3,2]);
savebj('',data2json,'NestArray',0)  % nestarray for 4-D or above is not working
json2data=loadbj(ans)
if(any(json2data(:)~=data2json(:)) || any(size(json2data)~=size(data2json)))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 3-D array in nested array form (JSONLab 1.9)\n')
fprintf(1,'%%=================================================\n\n')

data2json=reshape(1:(2*4*6),[2,4,6]);
savebj('',data2json,'NestArray',1,'FormatVersion',1.8)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 3-D array in annotated array form (JSONLab 1.9 or earlier)\n')
fprintf(1,'%%=================================================\n\n')

data2json=reshape(1:(2*4*6),[2,4,6]);
savebj('',data2json,'NestArray',0,'FormatVersion',1.8)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a complex number\n')
fprintf(1,'%%=================================================\n\n')

data2json=1+2i
savebj('',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a complex matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=magic(6);
data2json=data2json(:,1:3)+data2json(:,4:6)*1i
savebj('',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  MATLAB special constants\n')
fprintf(1,'%%=================================================\n\n')

data2json=[NaN Inf -Inf]
savebj('specials',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a real sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=sprand(10,10,0.1)
savebj('sparse',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a complex sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=data2json-data2json*1i
savebj('complex_sparse',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an all-zero sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse(2,3);
savebj('all_zero_sparse',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse([]);
savebj('empty_sparse',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty 0-by-0 real matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=[];
savebj('empty_0by0_real',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty 0-by-3 real matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=zeros(0,3);
savebj('empty_0by3_real',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse real column vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse([0,3,0,1,4]');
savebj('sparse_column_vector',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse complex column vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=data2json-1i*data2json;
savebj('complex_sparse_column_vector',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse real row vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse([0,3,0,1,4]);
savebj('sparse_row_vector',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse complex row vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=data2json-1i*data2json;
savebj('complex_sparse_row_vector',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a structure\n')
fprintf(1,'%%=================================================\n\n')

data2json=struct('name','Think Different','year',1997,'magic',magic(3),...
                 'misfits',[Inf,NaN],'embedded',struct('left',true,'right',false))
savebj('astruct',data2json,struct('ParseLogical',1))
json2data=loadbj(ans)
class(json2data.astruct.embedded.left)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a structure array\n')
fprintf(1,'%%=================================================\n\n')

data2json=struct('name','Nexus Prime','rank',9);
data2json(2)=struct('name','Sentinel Prime','rank',9);
data2json(3)=struct('name','Optimus Prime','rank',9);
savebj('Supreme Commander',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a cell array\n')
fprintf(1,'%%=================================================\n\n')

data2json=cell(3,1);
data2json{1}=struct('buzz',1.1,'rex',1.2,'bo',1.3,'hamm',2.0,'slink',2.1,'potato',2.2,...
              'woody',3.0,'sarge',3.1,'etch',4.0,'lenny',5.0,'squeeze',6.0,'wheezy',7.0);
data2json{2}=struct('Ubuntu',['Kubuntu';'Xubuntu';'Lubuntu']);
data2json{3}=[10.04,10.10,11.04,11.10]
savebj('debian',data2json,struct('FloatFormat','%.2f'))
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  invalid field-name handling\n')
fprintf(1,'%%=================================================\n\n')

json2data=loadbj(savebj('',loadjson('{"ValidName":1, "_InvalidName":2, ":Field:":3, "项目":"绝密"}')))

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a function handle\n')
fprintf(1,'%%=================================================\n\n')

data2json=@(x) x+1
savebj('handle',data2json)
json2data=loadbj(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 2D cell array\n')
fprintf(1,'%%=================================================\n\n')

data2json={{1,{2,3}},{4,5},{6};{7},{8,9},{10}};
savebj('data2json',data2json)
json2data=loadbj(ans)  % only savebj works for cell arrays, loadbj has issues

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 2D struct array\n')
fprintf(1,'%%=================================================\n\n')

data2json=repmat(struct('idx',0,'data','structs'),[2,3])
for i=1:6
    data2json(i).idx=i;
end
savebj('data2json',data2json)
json2data=loadbj(ans)

if(exist('datetime'))
    fprintf(1,'\n%%=================================================\n')
    fprintf(1,'%%  datetime object \n')
    fprintf(1,'%%=================================================\n\n')

    data2json=datetime({'8 April 2015','9 May 2015'}, 'InputFormat','d MMMM yyyy')
    savebj('',data2json)
    json2data=loadbj(ans)
end

if(exist('containers.Map'))
    fprintf(1,'\n%%=================================================\n')
    fprintf(1,'%%  a container.Maps object \n')
    fprintf(1,'%%=================================================\n\n')

    data2json=containers.Map({'Andy','William','Om'},{21,21,22})
    savebj('',data2json)
    json2data=loadbj(ans)
end

if(exist('istable'))
    fprintf(1,'\n%%=================================================\n')
    fprintf(1,'%%  a table object \n')
    fprintf(1,'%%=================================================\n\n')

    Names={'Andy','William','Om'}';
    Age=[21,21,22]';
    data2json=table(Names,Age)
    savebj('table',table(Names,Age))
    json2data=loadbj(ans)
end

if(exist('bandwidth'))
    fprintf(1,'\n%%=================================================\n')
    fprintf(1,'%%  use _ArrayShape_ \n')
    fprintf(1,'%%=================================================\n\n')

    lband=2;
    uband=3;
    data2json=spdiags(true(8,lband+uband+1),-uband:lband,5,8);
    data2json=full(double(data2json));
    data2json(data2json~=0)=find(data2json)

    savebj('',data2json,'usearrayshape',1)
    json2data=loadbj(ans,'fullarrayshape',1)

    savebj('',tril(data2json),'usearrayshape',1)
    json2data=loadbj(ans,'fullarrayshape',1)

    savebj('',triu(data2json+1i*data2json),'usearrayshape',1)
    json2data=loadbj(ans,'fullarrayshape',1)

    savebj('',tril(triu(int8(data2json))),'usearrayshape',1)
    json2data=loadbj(ans,'fullarrayshape',1)

    savebj('',data2json(:,1:5)+data2json(:,1:5)','usearrayshape',1)
    json2data=loadbj(ans,'fullarrayshape',1)
end

try
    val=zlibencode('test');
    fprintf(1,'\n%%=================================================\n')
    fprintf(1,'%%  a 2-D array in compressed array format\n')
    fprintf(1,'%%=================================================\n\n')

    data2json=eye(10);
    data2json(20,1)=1;
    savebj('',data2json,'Compression','zlib','CompressArraySize',0)  % nestarray for 4-D or above is not working
    json2data=loadbj(ans)
    if(any(json2data(:)~=data2json(:)) || any(size(json2data)~=size(data2json)))
        warning('conversion does not preserve original data');
    end
catch
end

rand ('state',rngstate);

author: fangq

repo: jsonlab

save_path: github-repos/MATLAB/fangq-jsonlab


text:
function M = rmCompareModelsGUI_sortVoxels(criterion, modelNum, M);
% Sort voxels in rmCompareModelsGUI.
%
% M = rmCompareModelsGUI_sortVoxels([criterion], [modelNum=dialog], [M=get from cur figure]);
%
% INPUTS:
%	criterion: can be one of 'ecc', 'pol', 'sigma', 'varexp', or 'beta'.
%	Will sort the voxels in ascending order for the specified parameter and
%	model #.
%	[default: get criterion from dialog]
%
%	modelNum: index of model in the comparison GUI from which to take the
%	voxel data to use for sorting.
%	[default: get model number from dialog]
%
%	M: comparison GUI struct, created by rmCompareModelsGUI_getData.
%	[default: get from cur figure, assumes you've already opened a
%	comparison GUI.]
%
%
%
% ras, 09/2008.
if notDefined('M'), M = get(gcf, 'UserData'); end
if notDefined('modelNum'),
	modelNum = listdlg('PromptString', 'Sort by Which Model?', ...
						'ListString', M.dtList, 'InitialValue', 1);
end

if notDefined('criterion'),
	cList = {'Variance Explained' 'Eccentricity' 'Polar Angle' ...
			 'x0' 'y0' 'pRF Size (sigma)' 'Scale Factor (beta)'};
	fieldList = {'varexp' 'ecc' 'pol' 'x0' 'y0' 'sigma' 'beta'};
	cNum = listdlg('PromptString', 'Sort by Which Criterion?', ...
						'ListString', cList, 'InitialValue', 1);
	criterion = fieldList{cNum};
end

% lookup for alternate names for criteria
switch lower(criterion)
	case 'eccentricity', criterion = 'ecc';
	case 'polar angle', criterion = 'pol';
	case 'prf size (sigma)', criterion = 'sigma';
	case 'scale factor (beta)', criterion = 'beta';
end

% get the values for the criterion
vals = M.(criterion){modelNum};

% special case: for beta, we want the first beta value (for the pRF, not
% the trend terms) in the model
if isequal(criterion, 'beta')
	vals = vals(:,1);
end

% sort the values, getting I -- the new voxel order
[sortedVals I] = sort(vals);

% reorder the ROI coordinates and data fields to match this:
M.roi.coords = M.roi.coords(:,I);
for f = {'x0' 'y0' 'sigma' 'pol' 'ecc' 'varexp'}
	for m = 1:M.nModels
		M.(f{1}){m} = M.(f{1}){m}(I);
	end
end

for m = 1:M.nModels
	M.beta{m} = M.beta{m}(I,:);
end

% reoder time series
for m = 1:M.nModels
	M.tSeries{m} = M.tSeries{m}(:,I);
end

% udpate the GUI
set(gcf, 'UserData', M);
% rmCompareModelsGUI_update;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [odf,alpha] = calcFEMODF(pf,varargin)
% PDF to ODF inversion
%
% *calcFEMODF* is one of the main function of the MTEX toolbox.
% It estimates an ODF from given Polefigure intensities by
% <PoleFigure2ODF.html fitting an ODF that consists of a large number of unimodal ODFs to the data>.
% It does so by minimizing a least squares functional. The command
% *calcODF* supports <PoleFigure2ODFGhostCorrection.html automatic ghost correction> and
% <PoleFigureDubna.html the zero range method>.
% The function *calcFEMODF* has several options to control convergence,
% resolution, smoothing, etc. See below for a complete description.
%
%
% Input
%  pf - @PoleFigure
%
% Options
%  resolution     - localization grid for the ansatz fucntions (default = 3/2 resolution(pf))
%  iterMax        - maximum number of iterations (default = 11)
%  regularisation - weighting coefficient lambda (default = 0)
%
% Flags
%  zeroRange         - apply zero range method (default = )
%  noGhostCorrection - omit ghost correction
%
% Output
%  odf    - reconstructed @FEMODF
%  alpha  - scaling factors, calculated during reconstruction
%
% See also
% PoleFigure2odf ODF_demo PoleFigureSimulation_demo
% PoleFigure.load ImportPoleFigureData examples_index

tic
vdisp('------ MTEX -- PDF to ODF inversion ------------------',varargin{:})

% ------------------- get input--------------------------------------------

CS = pf.CS; SS = pf.SS;

% generate FEM discretization of orientation space
res = get_option(varargin,'resolution',pf.resolution);
%ori = equispacedSO3Grid(CS,SS,'resolution',min(res,10*degree));
%DSO3 = DelaunaySO3(ori);
DSO3 = varargin{1};

% zero range method - TODO
%if check_option(varargin,'zero_range'), S3G = zero_range(pf,S3G,varargin{:});end

vdisp('Setting up matrices',varargin{:});
% compute matrices
for ipf = 1:length(pf)
  h = pf(ipf).h;
  vdisp(char(h),varargin{:});
  M{ipf} = sparse(length(pf(ipf).r),length(DSO3));
  for ih = 1:length(h)
    M{ipf} = M{ipf} + pf(ipf).c(ih) .* DSO3.pdfMatrix(h(ih),pf(ipf).r,varargin{:});
  end

  b{ipf} = pf(ipf).intensities(:);
end

M = vertcat(M{:});
b = vertcat(b{:});

vdisp(['starting solver'],varargin{:});
vdisp(['matrix is ' sizestr(M) ],varargin{:});
% solve the linear system of equation
c = lsqnonneg(M,b);
%c = M \ b;

% set up FEMODF
odf = femODF(DSO3,'weights',c);

end


text:
function yi = LagrangeInter(x,y,xi)

% SYNTAX:
%   yi = LagrangeInter(x,y,xi);
%
% INPUT:
%   x,y = row-vectors of (n+1) data values (x,y)
%   xi  = a row-vector of x-values, where interpolation is to be found (could be a single value)
%
% OUTPUT:
%   yi = a row-vector of interpolated y-values
%
% DESCRIPTION:
%   Lagrange interpolation algorithm.

%----------------------------------------------------------------------------------------------
%
% Author: Dmitry Pelinovsky
% Available online at: http://dmpeli.mcmaster.ca/Matlab/Math4Q3/Lecture2-1/LagrangeInter.m
%
%----------------------------------------------------------------------------------------------

n = length(x) - 1; % the degree of interpolation polynomial
ni = length(xi); % the number of x-values, where interpolation is to be found

L = ones(n+1,ni); % the matrix for Lagrange interpolating polynomials L_(n,k)(x)
                    % has (n+1) rows for each polynomial at k = 0,1,...,n
                    % has ni column for each x-value of xi

% Note: the algorithm uses the MATLAB capacities for matrix handling!
% The two nested loops below are designed to compute in parallel
% the values of Lagrange interpolating polynomials at each x-value of xi !

for k = 0 : n % start the outer loop through the data values for x

    for kk = 0 : (k-1) % start the inner loop through the data values for x (if k = 0 this loop is not executed)
        L(kk+1,:) = L(kk+1,:).*(xi - x(k+1))/(x(kk+1)-x(k+1)); % see the Lagrange interpolating polynomials
    end % end of the inner loop

    for kk = k+1 : n % start the inner loop through the data values (if k = n this loop is not executed)
        L(kk+1,:) = L(kk+1,:).*(xi - x(k+1))/(x(kk+1)-x(k+1));
    end % end of the inner loop

end % the end of the outer loop

% Now multiply the values for Lagrange interpolating polynomials by the data values for y
yi = y * L; % use matrix multiplication of row-vector y and the matrix L, the output is the row-vector


text:
function value = r8_log ( x )

%*****************************************************************************80
%
%% R8_LOG evaluates the logarithm of an R8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the logarithm of X.
%
  persistent alncen
  persistent alncs
  persistent center
  persistent nterms

  aln2 = 0.06814718055994530941723212145818;

  if ( isempty ( nterms ) )

    alncs = [ ...
      +0.13347199877973881561689386047187E+01, ...
      +0.69375628328411286281372438354225E-03, ...
      +0.42934039020450834506559210803662E-06, ...
      +0.28933847795432594580466440387587E-09, ...
      +0.20512517530340580901741813447726E-12, ...
      +0.15039717055497386574615153319999E-15, ...
      +0.11294540695636464284521613333333E-18, ...
      +0.86355788671171868881946666666666E-22, ...
      +0.66952990534350370613333333333333E-25, ...
      +0.52491557448151466666666666666666E-28, ...
      +0.41530540680362666666666666666666E-31 ]';

    center = [ 1.0, 1.25, 1.50, 1.75 ];

    alncen = [ ...
      +0.0, ...
      +0.22314355131420975576629509030983, ...
      +0.40546510810816438197801311546434, ...
      +0.55961578793542268627088850052682, ...
      +0.69314718055994530941723212145817 ]';

    nterms = r8_inits ( alncs, 11, 28.9 * r8_mach ( 3 ) );

  end

  if ( x <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_LOG - Fatal error!\n' );
    fprintf ( 1, '  X <= 0.0\n' );
    error ( 'R8_LOG - Fatal error!' )
  end

  [ y, n ] = r8_upak ( x );

  xn = n - 1;
  y = 2.0 * y;
  ntrval = r8_aint ( 4.0 * y - 2.5 );

  if ( ntrval == 5 )
    t = ( ( y - 1.0 ) - 1.0 ) / ( y + 2.0 );
  elseif ( ntrval < 5 )
    t = ( y - center(ntrval) ) / ( y + center(ntrval) );
  end

  t2 = t * t;
  value = 0.625 * xn + ( aln2 * xn + alncen(ntrval) ...
    + 2.0 * t + t * t2 ...
    * r8_csevl ( 578.0 * t2 - 1.0, alncs, nterms ) );

  return
end


text:
function a = moler3_llt ( alpha, n )

%*****************************************************************************80
%
%% MOLER1_LLT returns the Cholesky factor of the MOLER1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for j = 1 : n
    a(j,j) = 1.0;
    for i = j + 1 : n
      a(i,j) = alpha;
    end
  end

  return
end


text:
function K = k11(x, xp, hyp, i)

K = k11.ku1u1(x, xp, hyp, i);

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM