text:
% DEMHUMANPOSESVARGPLVMMISSING Run the shared variational GPLVM on various kinds of data
% allowing missing inputs.
% DESC Run the shared variational GPLVM on various kinds of data
% allowing missing inputs.
%
% COPYRIGHT: Andreas C. Damianou, Carl Henrik Ek,  2011
% SEEALSO : demHumanPoseSvargplvm1
%
% VARGPLVM

% Human pose data with the whole silhouette
%clear ;close all; experimentNo=404; imageSilhouette=1;initLatent='ppca';
%latentDimPerModel=7;dataSetNames =
%{};toyDataCreate='humanPose';initVardistIters = 380;
%itNo = [500 200 200 200 200 200 200 200 200 200];  indPoints=100;TEMPdemSharedVargplvm

%___

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

% Define constants (in a manner that allows other scripts to parametrize
% this one).
if ~exist('experimentNo') ,  experimentNo = 404;      end
if ~exist('itNo')         ,  itNo = [500 500];              end     % Default: 2000
if ~exist('indPoints')    ,  indPoints = 100;          end     % Default: 49
if ~exist('latentDimPerModel')    ,  latentDimPerModel = 6;          end
% Set to 1 to use dynamics or to 0 to use the standard var-GPLVM
if ~exist('dynUsed')      ,  dynUsed = 0;             end
if ~exist('initVardistIters'), initVardistIters = 180;      end
if ~exist('dynamicKern')   ,     dynamicKern = {'rbf', 'white', 'bias'}; end
if ~exist('mappingKern')   ,  mappingKern = {'rbfard2', 'white', 'bias'}; end

% 0.1 gives around 0.5 init.covars. 1.3 biases towards 0.
if ~exist('vardistCovarsMult'),  vardistCovarsMult=1.3;                  end
% Set to empty value {} to work with toy data

if ~exist('invWidthMultDyn'),    invWidthMultDyn = 100;                     end
if ~exist('invWidthMult'),       invWidthMult = 5;                     end
if ~exist('initX'),     initX ='ppca';   end % That's for the dynamics initialisation
if ~exist('dataType'), dataType = 'humanPose'; end
if ~exist('enableParallelism'), enableParallelism = 1; end
if ~exist('DgtN'), DgtN = false; end
% Create initial X by doing e.g. ppca in the concatenated model.m's or by
% doing ppca in the model.m's separately and concatenate afterwards?
if ~exist('initial_X'), initial_X = 'separately'; end % Other options: 'together'
% Which indices to use for training, rest for test
if ~exist('indTr'), indTr = -1; end

demHumanPosePrepareData

dataSetNames={'silhouette', 'pose'};
% X_init = Xp;
%mappingKern = {'linard2', 'white'};
%mappingKern = {'rbfard2', 'white'};
latentDim = 5; % Anything > 2 and < 10
%xyzankurAnim(Z_test, 3);
numberOfDatasets = length(Yall);

%-- Load datasets
for i=1:numberOfDatasets
    Y = Yall{i};

    dims{i} = size(Y,2);
    N{i} = size(Y,1);
    if indTr == -1
        indTr = 1:N{i};
    end
    %t{i} = linspace(0, 2*pi, size(Yall{i}, 1)+1)'; t{i} = t{i}(1:end-1, 1);

    indTs = setdiff(1:size(Y,1), indTr);
    Ytr{i} = Y(indTr,:);
    Yts{i} = Y(indTs,:);

    d{i} = size(Ytr{i}, 2);
end
% timeStampsTraining = t{1}(indTr,1); %timeStampsTest = t(indTs,1);

t = linspace(0, 2*pi, length(indTr)+1)'; t = t(1:end-1, 1);

% Fix times:
prevSeq = 1;
timeStampsTraining = [];
dt=0.05;
for i=1:length(seq)
    t = ([0:(seq(i)-prevSeq)].*dt)';
    prevSeq = seq(i)+1;
    timeStampsTraining = [timeStampsTraining ;t];
end;

dt=t(2)-t(1);
timeStampsTest = ([0:size(Y_test,1)-1].*dt)';

for i=2:numberOfDatasets
    if N{i} ~= N{i-1}
        error('The number of observations in each dataset must be the same!');
    end
end

model = svargplvmRestorePrunedModel(prunedModel, Ytr);
if isfield(model, 'dynamics') & ~isempty(model.dynamics)
    dynUsed=1;
else
    dynUsed = 0;
end
makePlots=0;

%%
%---------------------------- PREDICTIONS ---------------

% Set to 1 to test on the training data itself, set to 0 to use the test
% dataset.
if ~exist('testOnTraining')
    testOnTraining=0;
end

% 1 is for the HoG image features. 2 is for the pose features.
obsMod = 1; % one of the involved sub-models (possible values: 1 or 2).
infMod = setdiff(1:2, obsMod);

% Find the dimensions that are shared for obsMod and infMod
if ~exist('sharedDims')
    s1 = model.comp{obsMod}.kern.comp{1}.inputScales;
    s2 = model.comp{infMod}.kern.comp{1}.inputScales;
    % Normalise values between 0 and 1
    s1 = s1 / max(s1);
    s2 = s2 / max(s2);

    %  thresh = max(model.comp{obsMod}.kern.comp{1}.inputScales) * 0.001;
    thresh = 0.005;

    retainedScales{obsMod} = find(s1 > thresh);
    %thresh = max(model.comp{infMod}.kern.comp{1}.inputScales) * 0.001;
    retainedScales{infMod} = find(s2  > thresh);
    sharedDims = intersect(retainedScales{obsMod}, retainedScales{infMod});
end

% Find X_* only for the shared dimensions (Xs*):
if ~exist('privateDims')
    privateDims = setdiff(1:model.comp{obsMod}.q, sharedDims);
end

% Number of test points to use
numberTestPoints = 10;
if testOnTraining
    perm = randperm(model.N);
    testInd = perm(1:numberTestPoints);
else
    startInd = 60;
    missingInd = 1:25;
    Y_testOrig = Y_test;
    Y_test(startInd:end, missingInd) = NaN;

    %model.comp{2}.y(startInd:end, missingInd)=NaN; %%%???????????
    disp('# Initializing latent points...');
    indexPresent = setdiff(1:model.comp{1}.d, missingInd);

    Yts{obsMod} = Y_test;
    Yts{infMod} = Z_test;
    testInd = 1:size(Yts{1},1);
    ZpredAll = zeros(size(Z_test));
    indsAll = zeros(size(Z_test,1),1);
end

scrsz = get(0,'ScreenSize');

x_star = zeros(length(testInd), size(model.X,2));
if ~(~testOnTraining & dynUsed) % If we have a test set and dynamics, then we need a different inference procedure
    for i=1:length(testInd)
        curInd = testInd(i);
        fprintf('# Testing indice number %d ', curInd);
        if testOnTraining
            fprintf('taken from the training set\n');
            y_star = model.comp{obsMod}.y(curInd,:);
            x_star(i,:) = model.comp{obsMod}.vardist.means(curInd,:);
            varx_star = model.comp{obsMod}.vardist.covars(curInd,:);
        else
            fprintf('taken from the test set\n');
            y_star = Yts{obsMod}(curInd,:);
            z_star = Yts{infMod}(curInd,:);
            dst = dist2(y_star(indexPresent), Y_test(:,indexPresent));

            [mind, mini(i)] = min(dst);
            %miniAll(i) = mini;
            Init(i,:) = model.vardist.means(mini(i),:);
            vardistx = vardistCreate(model.comp{obsMod}.vardist.means(mini(i),:), model.q, 'gaussian');
            vardistx.covars = model.comp{obsMod}.vardist.covars(mini(i),:);
            model.comp{obsMod}.vardistx = vardistx;
            display=1;
            iters = 250;
            % Find p(X_* | Y_*) which is approximated by q(X_*)
            [x_star(i,:), varx_star, modelUpdated] = vargplvmOptimisePoint(model.comp{obsMod}, vardistx, y_star, display, iters);%%%
        end
    end
    fprintf('# Predicting images from the NN of X_* ');

    for i=1:length(testInd)
        curInd = testInd(i);
        y_star = Yts{obsMod}(curInd,:);
        z_star = Yts{infMod}(curInd,:);

        numberOfNN = 1;
        % Now we selected a datapoint X_* by taking into account only the
        % private dimensions for Y. Now, based on the shared dimensions of
        % that, we select the closest (in a NN manner) X from the training
        % data.
      %   w = s1(sharedDims)+s2(sharedDims);
      %   [ind, distInd] = nn_class2(model.X(:,sharedDims), x_star(i,sharedDims), numberOfNN, 'weighted', w);
          [ind, distInd] = nn_class(model.X(:,sharedDims), x_star(i,sharedDims), numberOfNN, 'euclidean');
        indsAll(i) = ind(1);

        ZpredMu = zeros(length(ind), size(model.comp{infMod}.y,2));
        ZpredSigma = zeros(length(ind), size(model.comp{infMod}.y,2));

        % Find p(y_*|x_*) for every x_* found from the NN
        for k=1:numberOfNN
            %fprintf('.');
            x_cur = model.X(ind(k),:);
            % Make the shared dimensions of the current NN the same as the
            % ones of the latent point for y_star
           %  x_cur(sharedDims) = x_star(i,sharedDims); %%% OPTIONAL!!!

            % Make the shared dimensions for the current NN the same as the
            % ones of the closest NN to y_star but from the Z dataset
            % x_cur(sharedDims) = model.X(ind(1),sharedDims); %%%%% OPTIONAL #2 !!!!

            %[ZpredMu(k,:), ZpredSigma(k,:)] = vargplvmPosteriorMeanVar(model.comp{infMod}, model.X(ind(k),:));
            %if ~testOnTraining
            ZpredMu(k,:) = vargplvmPosteriorMeanVar(model.comp{infMod}, x_cur);
            %else
            %    ZpredMu(k,:) = model.comp{infMod}(ind(k),:);
            %end
        end
        ZpredAll(i,:) = ZpredMu(1,:);

        %-- Plots
        if exist('makePlots') & ~makePlots
            continue
        end
        % Open a big figure (first 2 args control the position, last 2 control
        % the size)
        figure('Position',[scrsz(3)/100.86 scrsz(4)/6.666 scrsz(3)/1.0457 scrsz(4)/1.0682],...
            'Name',['Fig: ' num2str(i) ' (Exp: ' num2str(experimentNo) ')'],'NumberTitle','off')
        numRows = 3;

        if testOnTraining
            numCols = ceil((numberOfNN+1)/numRows)*2;
            plotCounter = 2;
        else
            % For the real test image!
            numCols = ceil((numberOfNN+2)/numRows);
            plotCounter = 2;
        end

        if ~testOnTraining
            sil_star = Yim_test(curInd,:);
            pose_star = Z_test(curInd,:);
            subplot(numRows, numCols, 1)
            imagesc(reshape(sil_star,height,width))
            if obsMod == 1
                title(['Given (image #' num2str(curInd) ')']), colormap('gray')
            else
                title('Corresponding');
            end
            subplot(numRows, numCols, 2)
            handle = xyzankurVisualise2(pose_star);
            if obsMod == 1
                title(['Given (image #' num2str(curInd) ')']), colormap('gray')
            else
                title('Corresponding');
            end

            for k=1:numberOfNN
                subplot(numRows, numCols, k+plotCounter)
                if infMod == 1 &&  exist('imageSilhouette') && imageSilhouette
                    imagesc(reshape( ZpredMu(k,:),height,width)), title(['NN #' num2str(k)]), colormap('gray')
                else
                    handle = xyzankurVisualise2(ZpredMu(k,:)); title(['NN #' num2str(k)])
                end
            end
            %mserror(i) = mean(abs(ZpredMu(1,:) - Z_test(i,:)));
        else
            subplot(numRows, numCols, 1)
            imagesc(reshape(Yim(curInd,:),height,width)), title(['Original y (image #' num2str(curInd) ')']), colormap('gray')
            subplot(numRows, numCols, 2)
            handle = xyzankurVisualise2(model.comp{2}.y(curInd,:));
            % Start plotting from 2, the first is always the same as
            for k=2:numberOfNN
                % Start from k=2, the first NN we know it's the same as the
                % given point.
                subplot(numRows, numCols, k+plotCounter-1)
                imagesc(reshape(Yim(ind(k),:),height,width)), title(['NN #' num2str(k)]), colormap('gray')
                subplot(numRows, numCols, k+plotCounter)
                handle = xyzankurVisualise2(model.comp{2}.y(ind(k),:)); title(['NN #' num2str(k)])
                plotCounter = plotCounter+1;
            end
        end
        %     pause
        %     close
    end
    if ~testOnTraining
        meanPose = repmat(mean(Ytr{2}),size(Y_test,1),1);
        errors.meanPose = xyzankurError(meanPose, Z_test);
        errors.NNYspace = xyzankurError(Ytr{2}(mini,:), Z_test);
        errors.NNXspace = xyzankurError(Ytr{2}(indsAll,:),Z_test);
        errors.svargplvm = xyzankurError(ZpredAll, Z_test);
        fprintf('# Mean Pose Error: %d\n', errors.meanPose)
        fprintf('# NN in the Y space Error: %d\n',errors.NNYspace)
        fprintf('# NN in the X space Error: %d\n',errors.NNXspace)
        fprintf('# Svargplvm Error: %d\n', errors.svargplvm)
        xyzankurAnimCompareMultipleTEMP(Z_test, {ZpredAll, Ytr{2}(mini,:), Ytr{2}(indsAll,:)},-1,{'Gr. Truth', 'Svargplvm', 'NN_Y','NN_X'});
    end
else
    %%% Dynamics

    model.dynamics.t_star = timeStampsTest;
    model.comp{1}.dynamics.t_star = timeStampsTest;
    model.comp{2}.dynamics.t_star = timeStampsTest;

    for i=1:size(Y_test,1)
        % initialize the latent points using the nearest neighbour from the training data
        dst = dist2(Y_test(i,indexPresent), Ytr{1}(:, indexPresent));
        [mind, mini(i)] = min(dst);
    end

    indexMissingData = startInd:size(Y_test,1);

    vardistx = vardistCreate(model.dynamics.vardist.means(mini,:), model.q, 'gaussian');
    vardistx.covars = 0.2*ones(size(vardistx.covars));
    model.vardistx = vardistx;
    model.comp{1}.vardistx = vardistx;
    model.comp{2}.vardistx = vardistx;

    iters = 4000;

    [x, varx] = vargplvmOptimiseSeqDyn(model.comp{obsMod}, vardistx, Y_test, 1, iters);
    % keep the optimized variational parameters
    barmu = x;
    lambda = varx;

    % Get the variational means and variacnes for the new test sequcen and
    % update the model to be prepared for prediction
    [x_star, varx_star, modelUpdated] = vargplvmDynamicsUpdateModelTestVar(model.comp{obsMod}, barmu, lambda, Y_test);

    % Latent variables corresponding to the data  with missing dimensions
    Testmeans = x(indexMissingData, :);
    Testcovars = varx(indexMissingData, :);

    modelOrig = model;
    model.comp{obsMod} = modelUpdated;
    model.vardist = modelUpdated.vardist;
    model.comp{1}.vardist = modelUpdated.vardist;
    model.comp{2}.vardist = modelUpdated.vardist;
    model.dynamics = modelUpdated.dynamics;
    model.comp{1}.dynamics = modelUpdated.dynamics;
    model.comp{2}.dynamics = modelUpdated.dynamics;
    model.X = model.vardist.means;
    model.comp{1}.X = model.X;
    model.comp{2}.X = model.X;

    numberOfNN = 2;

    fprintf('# Finding the %d NN of X_* with the training X based only on the shared dims.\n', numberOfNN);
    ZpredAll = zeros(size(Z_test));
    indsAll = zeros(size(Z_test,1),1);
    indsAllOrig = zeros(size(Z_test,1),1);
    fprintf('# Predicting images from the NN of X_* ');
    for i=1:size(Z_test,1)
        [ind2,distInd] = nn_class(modelOrig.X(:,sharedDims), x_star(i,sharedDims), numberOfNN, 'euclidean');
        indsAllOrig(i) = ind2(1);

      % Actually this is a weighted NN.
      %  w = s1(sharedDims)+s2(sharedDims);
      %  [ind, distInd] = nn_class2(model.X(:,sharedDims), x_star(i,sharedDims), numberOfNN, 'weighted', w);
        [ind, distInd] = nn_class(model.X(:,sharedDims), x_star(i,sharedDims), numberOfNN, 'euclidean');
        indsAll(i) = ind(1);
        ZpredMu = zeros(length(ind), size(model.comp{infMod}.y,2));
        ZpredSigma = zeros(length(ind), size(model.comp{infMod}.y,2));
        % Find p(y_*|x_*) for every x_* found from the NN
        for k=1:numberOfNN
            x_cur = model.X(ind(k),:);%%%%  % modelOrig.X(ind2(k));
            x_cur(sharedDims) = x_star(i,sharedDims); %%% OPTIONAL!!!

                        %---- OPTIONAL 3
          %  xcurOrig  = x_cur(sharedDims);
          %  s1new = s1/sum(s1);
          %  x_cur(sharedDims) = s1new(sharedDims).*x_star(i,sharedDims) + (1-s1new(sharedDims)).*xcurOrig;
            %----

            ZpredMu(k,:) = vargplvmPosteriorMeanVar(model.comp{infMod}, x_cur, varx_star(i,:)); % varx_star needed?
        end
        ZpredAll(i,:) = ZpredMu(1,:);
    end
    fprintf('\n');
    meanPose = repmat(mean(Ytr{2}),size(Z_test,1),1);
    errors.meanPose = xyzankurError(meanPose, Z_test);
    errors.NNYspace = xyzankurError(Ytr{2}(mini,:), Z_test);
    errors.NNXspace = xyzankurError(Ytr{2}(indsAllOrig,:),Z_test);
    errors.svargplvm = xyzankurError(ZpredAll, Z_test);
    fprintf('# Mean Pose Error: %d\n', errors.meanPose)
    fprintf('# NN in the Y space Error: %d\n',errors.NNYspace)
    fprintf('# NN in the X space Error: %d\n',errors.NNXspace)
    fprintf('# Svargplvm Error: %d\n', errors.svargplvm)
    prunedModelUpdated = vargplvmPruneModel(modelUpdated);
   % save(['demHumanPoseSvargplvm' num2str(experimentNo) '.mat'],'barmu','lambda','prunedModelUpdated','errors');

   % xyzankurAnimCompareMultipleTEMP(Z_test, {ZpredAll, Ytr{2}(mini,:), Ytr{2}(indsAllOrig,:)},-1,{'Gr. Truth', 'Svargplvm', 'NN_Y','NN_X'});
end

author: SheffieldML

repo: vargplvm


text:
function paramlist = mtrFilenames2Paramlist(filenames,paramnames)
%% Assume each name is *kSmooth_#smooth_kLength_#len_kMidSD_#mid* and stored in
%% paramData.  Grid will be returned in the following format:
%%
%% row 1: #smooth, #len, #mid
%% row 2: #smooth, #len, #mid
%%
%% ex: paramnames = {'kLength' 'kSmooth' 'kMidSD'};

paramlist = zeros(length(filenames),length(paramnames));
for nn = 1:length(filenames)
    for pp = 1:length(paramnames)
        paramlist(nn,pp) = mtrFilename2Param(filenames{nn},paramnames{pp});
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [nodes, edges2] = grRemoveEdge(nodes, edges, edge)
%GRREMOVEEDGE Remove an edge in a graph.
%
%   [NODES2 EDGES2] = grRemoveEdge(NODES, EDGES, EDGE2REMOVE)
%   Remove an edge in the edges list, and return the modified graph.
%   The NODES array is not modified.
%
%   -----
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 13/08/2003.
%

%   HISTORY
%   10/02/2004 doc

dim = size(edges);
edges2 = zeros(dim(1)-1, 2);
edges2(1:edge-1, :) = edges(1:edge-1, :);
edges2(edge:dim(1)-1, :) = edges(edge+1:dim(1), :);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% test for sparsity-based deblurring

path(path, 'toolbox/');

n = 128;
name = 'lena';
sigma = 0.04; % noise level

% load the image
M0 = rescale( crop( load_image(name), n), 0.05 );

% options for the wavelet transform
Jmin = 4;
options.Jmin = Jmin;
options.wavelet_type = 'biorthogonal_swapped';
options.wavelet_vm = 4;

% width of the filter
options.eta = 3;
% observation
Y = perform_blurring( M0, options.eta ) + sigma*randn(n);

% do inversion by iterative thresholding
lambda_min = 0;
lambda_max = 0.05;
nbr = 10;
lambda_list = linspace(lambda_max, lambda_min, nbr);

options.niter = 12; % iteration for inversion
options.mu = 1;
X = zeros(n^2,1); M = {}; perr = [];
for lambda = lambda_list
    options.lambda = lambda;
    options.X = X;
    X = perform_iterative_thresholding( @callback_blurring, Y(:), options );
    M{end+1} = perform_wavelet_transform(reshape(X,n,n), Jmin, -1, options);
    perr(end+1) = psnr(M0, M{end});
end

% plot error along regularization path
clf;
plot(lambda_list,perr);
xlabel('\lambda'); ylabel('PSNR');

% find best regularization lambda for inversion (this is cheating ...)
[pval,i] = max(perr);
clf;
imageplot(M0, 'Original', 1,3,1);
imageplot(Y, ['Observed, PSNR=' num2str(psnr(Y,M0))], 1,3,2);
imageplot(M{i}, ['Recovered, PSNR=' num2str(pval)], 1,3,3);

author: gpeyre

repo: matlab-toolboxes


text:
function [ff,f]=lpccc2ff(cc,np,nc,c0)
%LPCCC2FF Convert complex cepstrum to complex spectrum FF=(CC,NP,NC)
%
%  Inputs: cc(nf,n)     Complex ceptral coefficients excluding c(0), one frame per row
%          np           Size of output spectrum is np+1 [n]
%                       Alternatively, a vector of output frequencies in the range 0 to 0.5
%          nc           Number of cepstral coefficients to use [np or, if np is a vector, n]
%                       Set nc=-1 to use n coefficients
%          c0(nf,1)     Cepstral coefficient c(0) [0]
%
% Outputs: ff(nf,np+2)  Complex spectrum from DC to Nyquist
%          f(1,np+2)    Normalized frequencies (0 to 0.5)
%
% The "complex cepstral coefficients", cc(n), are the inverse discrete-time Fourier transform
% of the log of the complex-valued spectrum. The cc(n) are real-valued and, for n<0, cc(n)=0.
% The "real cepstral coeffcients", rc(n), are the inverse discrete-time Fourier transform
% of the log of the magnitude spectrum; rc(0)=cc(0) and rc(n)=0.5*cc(n) for n~=0.
% For highest speed, choose np+1 to be a power of 2.

%      Copyright (C) Mike Brookes 2014
%      Version: $Id: lpccc2ff.m 5025 2014-08-22 17:07:24Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nf,mc]=size(cc);
if nargin<2 || ~numel(np)
    if nargout
        np=mc;
    else
        np=128;
    end
end
if nargin>=3 && numel(nc)==1 && nc==-1 nc=mc; end
if nargin<4 || ~numel(c0) c0=zeros(nf,1); end
if numel(np)>1 || np(1)<1
    if nargin<3 || ~numel(nc) nc=mc; end
    f=np(:)';
    if nc==mc
        ff=exp([c0 cc]*exp(-2i*pi*(0:mc)'*f));
    else
        ff=exp([c0 lpccc2cc(cc,nc)]*exp(-2i*pi*(0:nc)'*f));
    end
else
    if nargin<3 || ~numel(nc) nc=np; end
    if nc==mc
        ff=exp(rfft([c0 cc].',2*np).');
    else
        ff=exp(rfft([c0 lpccc2cc(cc,nc)].',2*np).');
    end
    f=linspace(0,0.5,np+1);
end
if ~nargout
    subplot(2,1,2);
    plot(f,unwrap(angle(ff.')));
    xlabel('Normalized frequency f/f_s');
    ylabel('Phase (rad)');
    subplot(2,1,1);
    plot(f,db(abs(ff.')));
    xlabel('Normalized frequency f/f_s');
    ylabel('Gain (dB)');
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function ok = hlp_trycompile(varargin)
% Try to auto-compile a set of binary files in a folder, and return the status.
% OK = hlp_trycompile(Options...)
%
% This function tries to ensure that a given set of functions or classes (specified by their
% MATLAB identifier), whose source files are assumed to be located in a given directory, are
% properly compiled.
%
% The Style parameter determines how the function proceeds: Either compilation is always done
% ('force'), or only if necessary ('eager', e.g. if invalid file or changed source code). The check
% for re-compilation may be done on every call ('eager'), or once per MATLAB session ('lazy').
%
% The most common use case is specifying a given directory (and omitting the identifiers). In this
% case, all source files that have a mexFunction declaration are compiled, and all other source
% files are also supplied to the compiler as additional files. In case of compilation errors,
% hlp_trycompile also tries to omit all additional files during compilation. Both the list of
% additional files (or their file name patterns) or the considered file types can be specified. The
% list of identifiers to consider can also be specified.
%
% Another possible use case is to omit both the identifiers and the directory. In this case,
% hlp_trycompile assumes that a mex file with the same identifier (and path) as the calling function
% shall be compiled. This is would be used in .m files which directly implement some fallback code
% in case that the compilation fails (or which are just stubs to trigger the on-demand compilation).
%
%
% Since there can be many different versions of a mex binary under Linux (even with the same name),
% .mex files are by default moved into a sub-directory (named according to the hostname) after
% compilation. This does not apply to .class files, which are not platform-specific.
%
% The function supports nearly all features of the underlying MEX compiler, and can thus be used
% to compile a large variety of mex packages found in the wild (in some cases with custom defines,
% libraries, or include directories).
%
% If you are debugging code with it, it is best to set Verbose to true, so that you get compilation
% output.
%
% Additional features of this function include:
%  * Figures out whether the -largeArrayDims switch should be used.
%  * Repeated calls of this function are very fast if used in 'lazy' mode (so that it can be used in
%    an inner loop).
%  * Automatically rebuilds if the mex source has changed (does not apply to misc dependency files),
%    if used in 'eager' mode.
%  * By default uses the Mathworks versions of BLAS and LAPACK (if these libraries are pulled in).
%  * Supports both '/' and '\' in directory names.
%  * Behaves reasonably in deployed mode (i.e. gives warnings if files are missing).
%  * Also compiles .java files where appropriate.
%  * Supports test code.
%
% If this function produces errors for some mex package, the most common causes are:
%  * If the platform has never been used to compile code, an appropriate compiler may have to be
%    selected using "mex -setup". If no supported compiler is installed (e.g. on Win64), it must
%    first be doenloaded and installed (there is a free standard compiler for every platform).
%  * Some unused source files are in the directory which produce errors when they are automatically
%    pulled in.
%    --> turn on verbose output and identify & remove these (or check the supplied make file for
%        what files are actually needed)
%  * The functions require a custom define switch to work.
%    --> Check the make file, and add the switch(es) using the 'Defines' parameter.
%  * The functions use non-standard C code (e.g. // comments).
%    --> Tentatively rename the offending .c files into .cpp.
%  * The functions require a specific library to work.
%    --> Check the make file, and add the libraries using the 'Libaries' parameter.
%  * The functions require specific include directories to work.
%    --> Check the make file, and add the directories using the 'IncludeDirectories' parameter.
%  * The functions require additional files that are in a different directory.
%    --> Check the make file, and add these files using the 'SupportFiles' parameter. Wildcards are
%        allowed (in particular the special '*' string, which translates into all source files in
%        the Directory).
%  * The package assumes that mex is used with the -output option to use a custom identifier name
%    --> This type of make acrobatic is not supported by hlp_trycompile; instead, rename the source
%        file which has the mexFunction definition such that it matches the target identifier.
%  * The functions require specific library directories to work.
%    --> Check the make file, and add the directories using the 'LibraryDirectories' parameter.
%
%
% In:
%   Style : execution style, can be one of the following (default: 'lazy')
%           'force' : force compilation (regardless of whether the binaries are already there)
%           'eager' : compile only if necessary, check every time that this function is called
%           'lazy'  : compile only if necessary, and don't check again during this MATLAB session
%
%   --- target files ---
%
%   Directory : directory in which the source files are located
%               (default: directory of the calling function)
%
%   Identifiers : identifier of the target function/class, or cell array of identifiers that should
%                 be compiled (default: Calling function, if no directory given, or names of all
%                 compilable source files in the directory, if a directory is given.)
%
%   FileTypes : file type patterns to consider as sources files for the Identifiers
%               (default: {'*.f','*.c','*.cpp','*.java'})
%
%
%   --- testing conditions ---
%
%   TestCode : MATLAB code (string) which evaluates to true if the compiled code is behaving
%              correctly (and false otherwise), or alternatively a function handle which does the
%              same
%
%
%   --- additional compiler inputs ---
%
%   SupportFiles : cell array of additional / supporting source filenames to include in the compilation of all
%                  Identifiers (default: '*')
%                  Note: Any file listed here will not be considered part of the Identifiers, when
%                        all contents of a directory are to be compiled.
%                  Note: If there are support source files in sub-directories, include the full path
%                        to them.
%                  Note: If this is '*', all source files that are not mex files in the given
%                        directory are used as support files.
%
%   Libraries : names of libraries to include in the compilation
%               (default: {})
%
%   IncludeDirectories : additional directories in which to search for included source files.
%                        (default: {})
%
%   LibraryDirectories : additional directories in which to search for referenced library files.
%                        (default: {})
%
%   Defines : list of defined symbols (either 'name' or 'name=value' strings)
%             (default: {})
%
%   Renaming : cell array of {sourcefile,identifier,sourcefile,identifier, ...} indicating that
%              the MEX functions generated from the respective source files should be renamed to
%              the given identifiers. Corresponds to MEX's -output option; does not apply to Java files.
%              (default: {})
%
%   Arguments : miscellaneous compiler arguments (default: {})
%               For possible arguments, type "help mex" in the command line
%
%   DebugBuild : whether to build binaries in debug mode (default: false)
%
%
%   --- user messages ---
%
%   ErrorMessage : the detail error message to display which describes what type of functionality
%                  will not be available (if any).
%                  Note: If you have a MATLAB fallback, mention this in the error message.
%
%   PreparationMessage : the message that will be displayed before compilation begins.
%                        (default: {})
%
%   Verbose : whether to display verbose compiler outputs (default: false)
%
%
%   --- misc options ---
%
%   MathworksLibs : whether to use the Mathworks versions of std. libraries instead of OS-supplied
%                   ones, if present (applies to blas and lapack) (default: true)
%
%   DebugCompile : debug the compilation process; halts in the directory prior to invoking mex
%                  (default: false)
%
%
% Examples:
%   % try to compile all mex / Java files in a given directory:
%   hlp_trycompile('Directory','/Extern/MySources');
%
%   % as before, but restrict the set of identifiers to compile to a given set
%   hlp_trycompile('Directory','/Extern/MySources','Identifiers',{'svmtrain','svmpredict'});
%
%   % try to compile mex / Java files in a given directory, and include 2 libraries in the compilation
%   hlp_trycompile('Directory','/Extern/MySources','Libraries',{'blas','lapack'});
%
%   % like before, but this time include additional source files from two other directories
%   % (the single '*' means: include non-mex sources in the specified directory)
%   hlp_trycompile('Directory','/Extern/MySources','SupportFiles',{'*','../blas/*.c','../*.cpp'});
%
%   % like before, but this time add an include directory, a library directory, and some library
%   hlp_trycompile('Directory','/Extern/MySources', 'IncludeDirectories','/boost/include','LibraryDirectories','/boost/lib','Libraries','boost_date_time-1.44');
%
%   % like before, this time specifying some custom #define's
%   hlp_trycompile('Directory','/Extern/MySources','Defines',{'DEBUG','MAX=((a)>(b)?(a):(b))'});
%
% Debugging broken compiles:
%   To debug a broken compilation, you can switch on the 'Verbose' flag, which will show the output;
%   if you are debugging automatic compilation of a BCILAB dependendency, you can add Verbose=true;
%   to the env_compile.m file and reload the toolbox, or you can manually call hlp_trycompile as in:
%   hlp_trycompile('Directory',<<path-to-dependency>>,'Style','force','Verbose',true, <<any-settings-from-env_compile.m-as-name-value-pairs>>)
%   The stream of outputs can be a bit confusing and doesn't show all error messages, so you may
%   have to also put a breakpoint at the place in the code where it says PUT BREAKPOINT HERE, and
%   then manually invoke the next doeval() statement that it encounters in the console.
%
% Use cases:
%   1) In addition to a source file mysvd.c (compiling into mysvd.mex*), a stub .m file of the same
%   name can be placed in the same directory, which contains code to compile the binary when needed.
%
%   function [U,S,V] = mysvd(X)
%   if hlp_trycompile
%       [U,S,V] = mysvd(X);
%   else
%       % either display an error message or implement some fallback code.
%   end
%
%   2) In a MATLAB function which makes use of a few mex files, ensure compilation of these files.
%   function myprocessing(X,y)
%   if ~hlp_trycompile('Identifiers',{'svmtrain.c','svmpredict.c'})
%       error('Your binary files could not be compiled.');
%   else
%       m = svmtrain(X,y);
%       l = svmpredict(X,m);
%       ...
%   end
%
%   3) In a startup script.
%   hlp_trycompile('Directory','/Extern/MySources');
%
% See also:
%   mex
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-03-09

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

persistent results;           % a map of result-tag to OK value
persistent compiler_selected; % whether a compiler has been selected (true/false, or [] if uncertain)

% read options
o = hlp_varargin2struct(varargin, ...
    ... % overall behavior
    {'style','Style'}, 'lazy', ...
    ... % target files
    {'dir','Directory'}, [], ...
    {'idents','Identifiers'}, [], ...
    {'types','FileTypes'}, {'*.c','*.C','*.cpp','*.CPP','*.f','*.F','*.java','*.Java'}, ...
    ... % test condition
    {'test','TestCode'},'', ...
    ... % additional compiler inputs
    {'support','SupportFiles'}, {'*'}, ...
    {'libs','Libraries'}, {}, ...
    {'includedirs','IncludeDirectories'}, {}, ...
    {'libdirs','LibraryDirectories'}, {}, ...
    {'defines','Defines'}, {}, ...
    {'args','Arguments'}, '', ...
    {'renaming','Renaming'}, {}, ...
    {'debug','DebugBuild'}, false, ...
    ... % messages
    {'errmsg','ErrorMessage'}, {'Some BCILAB functionality will likely not be available.'}, ...
    {'prepmsg','PreparationMessage'}, {}, ...
    {'verbose','Verbose'}, false, ...
    ... % misc
    {'mwlibs','MathworksLibs'}, true, ...
    {'debugcompile','DebugCompile'}, false ...
    );

% support for parameterless calls
if isempty(o.dir)
    % if no dir given, use the calling function's directory
    [name,file] = hlp_getcaller();
    o.dir = fileparts(file);
    if isempty(file)
        error('If hlp_trycompile is called without a directory, it must be called from within a file.'); end
    % if neither idents nor dir given, use the calling function's identifier
    if isempty(o.idents)
        o.idents = name; end
end

% uniformize ident format
if isa(o.idents,'function_handle')
    o.idents = char(o.idents); end
if ischar(o.idents)
    o.idents = {o.idents}; end
if isempty(o.idents)
    o.idents = {}; end

% decide whether a re-check can be skipped based on identifiers and directory
if strcmp(o.style,'lazy') || isdeployed
    str = java.lang.String(sprintf('%s:',o.dir,o.idents{:}));
    tag = sprintf('x%.0f',str.hashCode()+3^31);
    if isfield(results,tag)
        ok = results.(tag);
        return;
    end
end

% uniformize directory format
o.dir = path_normalize(o.dir);

% verify style
if ~any(strcmp(o.style,{'force','eager','lazy'}))
    error('Unsupported style: %s',o.style); end

% uniformize test condition
if isa(o.test,'function_handle')
    o.test = char(o.test); end

% uniformize user messages
if ischar(o.errmsg)
    o.errmsg = {o.errmsg}; end
if ischar(o.prepmsg)
    o.prepmsg = {o.prepmsg}; end

% uniformize types
if ischar(o.types)
    o.types = {o.types}; end
for t=1:length(o.types)
    if o.types{t}(1) ~= '*'
        o.types{t} = ['*' o.types{t}]; end
end

% uniformize compiler inputs
if ischar(o.support)
    o.support = {o.support}; end
if ischar(o.includedirs)
    o.includedirs = {o.includedirs}; end
if ischar(o.libdirs)
    o.libdirs = {o.libdirs}; end
if ischar(o.libs)
    o.libs = {o.libs}; end
if ischar(o.defines)
    o.defines = {o.defines}; end
if ischar(o.args)
    o.args = {o.args}; end
for i=1:length(o.support)
    o.support{i} = path_normalize(o.support{i}); end
for i=1:length(o.includedirs)
    o.includedirs{i} = path_normalize(o.includedirs{i}); end
for i=1:length(o.libdirs)
    o.libdirs{i} = path_normalize(o.libdirs{i}); end

% if a support is given as '*'
starred = strcmp('*',o.support);
if any(starred)
    % list all in the given directory that are not mex files
    infos = [];
    for t = 1:length(o.types)
        if ~isempty(infos)
            infos = [infos; dir([o.dir filesep o.types{t}])];
        else
            infos = dir([o.dir filesep o.types{t}]);
        end
    end
    fnames = {infos.name};
    supportfiles = ~cellfun(@(n)is_primary([o.dir filesep n]),fnames);
    o.support = [fnames(supportfiles) o.support(~starred)];
end

% infer directory, if not given (take it from the calling function)
if isempty(o.idents) && ~isempty(o.dir)
    % get all the source files in the given direcctory
    infos = [];
    for t = 1:length(o.types)
        if ~isempty(infos)
            infos = [infos; dir([o.dir filesep o.types{t}])];
        else
            infos = dir([o.dir filesep o.types{t}]);
        end
    end
    fnames = {infos.name};
    % ... but exclude the support files
    fnames = setdiff(fnames,o.support);
    % and apply any renamings to get the corresponding identifiers
    if ~isempty(o.renaming)
        for n=1:length(fnames)
            fnames{n} = hlp_rewrite(fnames{n},o.renaming{:}); end
    end
    % ... and strip off the extensions
    for n=1:length(fnames)
        fnames{n} = hlp_getresult(2,@fileparts,fnames{n}); end
    o.idents = fnames;
end

ok = false;
missingid = []; % indices of missing identifiers (for dont-retry-next-time beacon files)
if isdeployed
    % --- deployed mode ---
    % Can not compile, but figure out whether everything needed is present. A special consideration
    % is that both the mex files calling functions are in a mounted .ctf archive.

    % check if all identifiers are present (either as mex or class)
    for i=1:length(o.idents)
        if ~any(exist(o.idents{i}) == [2 3 8])
            missingid(end+1) = i; end
    end
    ok = isempty(missingid);

    if ~isempty(missingid)
        % not all identifiers are compiled for this platform
        disp(['Note: The MEX functions/identifiers ' format_cellstr(o.idents(missingid)) ' are not included for your platform.']);
    elseif strcmp(o.style,'force')
        % in force mode, we remark that everything is already compiled
        disp_once(['The functions ' format_cellstr(o.idents) ' are properly compiled.']);
    end

else
    % --- regular mode ---
    % here, we *do* compile what needs to be compiled

    % find out a key configuration settings
    is64bit = ~isempty(strfind(computer,'64'));
    has_largearrays = is64bit && hlp_matlab_version >= 703;

    if ispc
        warning off MATLAB:FILEATTRIB:SyntaxWarning; end

    % add a few missing defines
    if hlp_matlab_version < 703
        o.defines = [o.defines {'mwIndex=int','mwSize=int','mwSignedIndex=int'}]; end

    % rewrite blas & lapack libs....
    if o.mwlibs
        % for each type of library
        for l={'blas','lapack'}
            lib = l{1};
            % note: this code is based on SeDuMi's compile script (by Michael C. Grant)
            in_use = strcmp(o.libs,lib);
            if any(in_use)
                if ispc
                    if is64bit
                        osdir = 'win64';
                    else
                        osdir = 'win32';
                    end
                    libpath = [matlabroot '\extern\lib\' osdir '\microsoft\libmw' lib '.lib'];
                    if ~exist(libpath,'file')
                        libpath = [matlabroot '\extern\lib\' osdir '\microsoft\msvc60\libmw' lib '.lib']; end
                    if exist(libpath,'file')
                        o.libs{in_use} = libpath;
                    else
                        disp_once('Note: The Mathworks library %s was assumed to be in %s, but not found.',lib,libpath);
                    end
                else
                    o.libs{in_use} = ['mw' lib];
                end
            end
        end
    end

    try
        % remember the current directory & enter the target directory
        olddir = pwd;
        if hlp_matlab_version >= 706
            go_back = onCleanup(@()cd(olddir)); end
        if ~exist(o.dir,'dir')
            error(['The target directory ' o.dir ' does not exist.']); end
        cd(o.dir);

        % expand regex patterns in o.support
        for i=length(o.support):-1:1
            if any(o.support{i} == '*')
                found = dir(o.support{i});
                if ~isempty(found)
                    % ... and splice the results in
                    basepath = fileparts(o.support{i});
                    items = cellfun(@(x)[basepath filesep x],{found.name},'UniformOutput',false);
                    o.support = [o.support(1:i-1) items o.support(i+1:end)];
                end
            end
        end

        % find all source & target files for the respective identifiers...
        % (note that there might be multiple source files for each one)
        sources = cell(1,length(o.idents)); % list of all source file names for the corresponding identifiers (indexed like idents)
        targets = cell(1,length(o.idents)); % list of all target file names for the corresponding identifiers (indexed like idents)
        for i=1:length(o.idents)
            if ~isempty(o.renaming)
                % the renaming may yield additional source file names for the given identifiers
                idx = strcmp(o.idents{i},o.renaming(2:2:end));
                if any(idx)
                    % the identifier is a renaming target: add the corresponding source file name
                    filename = o.renaming{find(idx)*2-1};
                    % if a source file with this ident & type is present
                    if exist([o.dir filesep filename],'file')
                        % remember it & derive its respective target file name
                        sources{i}{end+1} = filename;
                        targets{i}{end+1} = [o.idents{i} '.' mexext];
                    end
                end
            end
            for t=1:length(o.types)
                filename = [o.idents{i} o.types{t}(2:end)];
                % if a source file with this ident & type is present
                if exist([o.dir filesep filename],'file')
                    % remember it
                    sources{i}{end+1} = filename;
                    % and also derive its respective target file name
                    if strcmp(o.types{t},'*.java')
                        targets{i}{end+1} = [o.idents{i} '.class'];
                    else
                        targets{i}{end+1} = [o.idents{i} '.' mexext];
                    end
                end
            end
            % check whether we have all necessary source files
            if isempty(sources{i})
                error('Did not find source file for %s',o.idents{i}); end
            if isempty(targets{i})
                error('Could not determine target file for %s',o.idents{i}); end
        end

        % check for existence (either .mex* or class) of all identifiers
        % and make a list of missing & present binary files; do this in a different directory,
        % to not shadow the mex files of interest with whatever .m files live in this directory
        cd ..
        binaries = {}; % table of existing binary file paths (indexed like idents)
        for i=1:length(o.idents)
            % get current file reference to this identifier
            binaries{i} = which(o.idents{i});
            % if it doesn't point to a .mex or .class file, ignore it
            if exist(o.idents{i}) == 7
                warning('There is a directory named %s on the path, which clashes with the name of a binary to be compiled.'); end
            if ~any(exist(o.idents{i}) == [3 8])
                binaries{i} = []; end
            if ~isempty(binaries{i})
                % check whether it is correct file path
                if ~any(binaries{i} == filesep)
                    error(['Could not determine the location of the mex file for: ' o.idents{i}]); end
                % check whether the referenced file actually exists in the file system
                if isempty(dir(binaries{i}))
                    binaries{i} = []; end
            end
            % if no binary found, record it as missing
            if isempty(binaries{i})
                missingid(end+1) = i; end
        end
        cd(o.dir);

        % check which of the existing binaries need to be re-compiled (if out of date)
        outdatedid = []; % indices of identifiers (in o.idents) that need to be recompiled
        for i=1:length(binaries)
            if ~isempty(binaries{i})
                % get the date of the binary file
                bininfo = dir(binaries{i});
                % find all corresponding source files
                srcinfo = [];
                for s=1:length(sources{i})
                    srcinfo = [srcinfo; dir([o.dir filesep sources{i}{s}])]; end
                if ~isfield(bininfo,'datenum')
                    [bininfo.datenum] = celldeal(cellfun(@datenum,{bininfo.date},'UniformOutput',false)); end
                if ~isfield(srcinfo,'datenum')
                    [srcinfo.datenum] = celldeal(cellfun(@datenum,{srcinfo.date},'UniformOutput',false)); end
                % if any of the source files has been changed
                if bininfo.datenum < max([srcinfo.datenum])
                    % check if their md5 hash is still the same...
                    if exist([binaries{i} '.md5'],'file')
                        try
                            contents = load([binaries{i} '.md5'],'-mat','srchash');
                            % need to do that over all source files...
                            srchash = [];
                            sorted_sources = sort(sources{i});
                            for s=1:length(sorted_sources)
                                srchash = [srchash hlp_cryptohash([o.dir filesep sorted_sources{s}],true)]; end
                            if ~isequal(srchash,contents.srchash)
                                % hash is different: mark binary as outdated
                                outdatedid(end+1) = i;
                            end
                        catch
                            % there was a probblem: mark as outdated
                            outdatedid(end+1) = i;
                        end
                    else
                        % no md5 present: mark as outdated
                        outdatedid(end+1) = i;
                    end
                end
            end
        end

        % we try to recompile both what's missing and what's outdated
        recompileid = [missingid outdatedid];
        javainvolved = false; % for final error/help message generation
        mexinvolved = false;  % same

        if ~isempty(recompileid)
            % need to recompile something -- display a few preparatory messages...
            for l=1:length(o.prepmsg)
                disp(o.prepmsg{l}); end

            failedid = []; % list of indices of identifier that failed the build

            % for each identifier, try to compile it
            % and record whether it failed
            for i=recompileid
                success = false;
                fprintf(['Compiling the function/class ' o.idents{i} '...']);
                % for each source file mapping to that identifier
                for s=1:length(sources{i})
                    % check type of source
                    if ~isempty(strfind(sources{i}{s},'.java'))
                        % we have a Java source file: compile
                        [errcode,result] = system(['javac ' javac_options(o) sources{i}{s}]);
                        if errcode
                            javainvolved = true;
                            % problem: show display output
                            fprintf('\n');
                            disp(result);
                        else
                            success = true;
                            break;
                        end
                    else

                        % generate MEX options
                        opts = mex_options(o);
                        o.support = unique(o.support);
                        supp = sprintf(' %s',o.support{:});

                        if isempty(compiler_selected)
                            % not clear whether a compiler has been selected yet
                            if hlp_matlab_version >= 708
                                % we can find it out programmatically
                                try
                                    cconf = mex.getCompilerConfigurations; %#ok<NASGU>
                                    compiler_selected = true;
                                catch
                                    % no compiler has been selected yet...
                                    try
                                        % display a few useful hints to the user
                                        disp(' to compile this feature, you first need to select');
                                        disp('which compiler should be used on your platform.');
                                        if ispc
                                            if is64bit
                                                disp_once('As you are on 64-bit windows, you may find that no compiler is installed.');
                                            else
                                                disp_once('On 32-bit Windows, MATLAB supplies a built-in compiler (LLC), which should');
                                                disp_once('faithfully compile most C code. For broader support across C dialects (as well as C++), ');
                                                disp_once('you should make sure that a better compiler is installed on your system and selected in the following.');
                                            end
                                            disp_once('A good choice is the free Microsoft Visual Studio 2005/2008/2010 Express compiler suite');
                                            disp_once('together with the Microsoft Platform SDK (6.1 for 2008, 7.1 for 2010) for your Windows Version.');
                                            disp_once('See also: http://argus-home.coas.oregonstate.edu/forums/development/core-software-development/compiling-64-bit-mex-files');
                                            disp_once('          http://www.mathworks.com/support/compilers/R2010b/win64.html');
                                            disp_once('The installation is easier if a professional Intel or Microsoft compiler is used.');
                                        elseif isunix
                                            disp_once('On Linux/UNIX, the best choice is usually a supported version of the GCC compiler suite.');
                                        else
                                            disp_once('On Mac OS, you need to have a supported version of Xcode/GCC installed.');
                                        end
                                        % start the compiler selection tool
                                        mex -setup
                                        % verify that a compiler has been selected
                                        cconf = mex.getCompilerConfigurations; %#ok<NASGU>
                                        compiler_selected = true;
                                    catch
                                        compiler_selected = false;
                                    end
                                end
                            else
                                disp(' you may be prompted to select a compiler in the following');
                                disp('(as BCILAB cannot auto-determine whether one is selected on your platform).');
                            end
                        end

                        if ~compiler_selected
                            fprintf('skipped (no compiler selected).\n');
                        else

                            if o.verbose || isempty(compiler_selected)
                                % this variant will also be brought up if not sure whether a compiler
                                % has already been selected...
                                doeval = @eval;
                            else
                                doeval = @evalc;
                            end

                            % try to build the file
                            try
                                mexinvolved = true;
                                % check if a renaming applies...
                                idx = strcmp(sources{i}{s},o.renaming);
                                if any(idx)
                                    rename = [' -output ' o.renaming{find(idx,1)+1} ' '];
                                else
                                    rename = '';
                                end
                                if o.debugcompile
                                    % display a debug console to allow the user to debug how their file compiles
                                    fprintf('\nExecution has been paused immediately before running mex.\n');
                                    disp(['You are in the directory "' pwd '".']);
                                    disp('The mex command that will be invoked in the following is:');
                                    if has_largearrays
                                        disp([' mex ' opts ' -largeArrayDims ' rename sources{i}{s} supp]);
                                    else
                                        disp([' mex ' opts rename sources{i}{s} supp]);
                                    end
                                    fprintf('\n\nTo proceed normally, type "dbcont".\n');
                                    keyboard;
                                end
                                if has_largearrays   % PUT BREAKPOINT HERE
                                    try
                                        % -largeArrayDims enabled
                                        try
                                            doeval(['mex ' opts ' -largeArrayDims ' rename sources{i}{s} supp]); % with supporting libaries
                                        catch
                                            doeval(['mex ' opts ' -largeArrayDims ' rename sources{i}{s}]); % without supporting libraries
                                        end
                                    catch
                                        % -largeArrayDims disabled
                                        try
                                            doeval(['mex ' opts rename sources{i}{s} supp]); % with supporting libaries
                                        catch
                                            doeval(['mex ' opts rename sources{i}{s}]); % without supporting libraries
                                        end
                                    end
                                else
                                    % -largeArrayDims disabled
                                    try
                                        doeval(['mex ' opts rename sources{i}{s} supp]); % with supporting libaries
                                    catch
                                        doeval(['mex ' opts rename sources{i}{s}]); % without supporting libraries
                                    end
                                end

                                % compilation succeeded...
                                if any(i==outdatedid)
                                    % there is an outdated binary, which needs to be deleted
                                    try
                                        delete(binaries{i});
                                    catch
                                        disp(['Could not delete outdated binary ' binaries{i}]);
                                    end
                                end
                                % check whether the file is being found now
                                if exist(o.idents{i}) == 3
                                    success = true;
                                    compiler_selected = true;
                                    break;
                                end
                            catch
                                % build failed
                            end
                        end
                    end
                end

                % check if compilation of this identifier was successful
                if success
                    % if so, we sign off the binary with an md5 hash of the sources...
                    newbinary = which(o.idents{i});
                    try
                        srchash = [];
                        sorted_sources = sort(sources{i});
                        for s=1:length(sorted_sources)
                            srchash = [srchash hlp_cryptohash([o.dir filesep sorted_sources{s}],true)]; end
                        save([newbinary '.md5'],'srchash','-mat');
                        fprintf('success.\n');
                    catch
                        disp('could not create md5 hash for the source files; other than that, successful.');
                    end
                else
                    fprintf('failed.\n');
                    failedid(end+1) = i;
                end
            end

            embed_test = false;
            if isempty(failedid)
                % all worked: now run the test code - if any - to verify the correctness of the build
                if length(recompileid) > 1
                    if ~isempty(o.test)
                        fprintf('All files in %s compiled successfully; now testing the build outputs...',o.dir);
                    else
                        fprintf('All files in %s compiled successfully.\n',o.dir);
                    end
                elseif ~isempty(o.test)
                    fprintf('Now testing the build outputs...');
                end

                % test the output
                try
                    ans = true; %#ok<NOANS>
                    eval(o.test);
                catch
                    ans = false; %#ok<NOANS>
                end
                if ans %#ok<NOANS>
                    % the test was successful; now copy the files into a platform-specific directory
                    if ~isempty(o.test)
                        % only if we have a succeeding non-empty test
                        embed_test = true;
                        fprintf('success.\n');
                    end
                    retainid = recompileid;
                    eraseid = [];
                    ok = true;
                else
                    % the test was unsuccessful: remove all newly-compiled files...
                    if ~isempty(o.test)
                        fprintf('failed.\n'); end
                    disp('The code compiled correctly but failed the build tests. Reverting the build...');
                    disp('If this is unmodified BCILAB code, please consider reporting this issue.');
                    retainid = [];
                    eraseid = recompileid;
                end
            else
                if length(recompileid) > 1
                    if isempty(setdiff(recompileid,failedid))
                        fprintf('All files in %s failed to build; this indicates a problem in your build environment/settings.\n',o.dir);
                    else
                        fprintf('Some files in %s failed to build. Please make sure that you have a supported compiler; otherwise, please report this issue.\n',o.dir);
                    end
                else
                    disp('Please make sure that you have a supported compiler and that your build environment is set up correctly.');
                    disp('Also, please consider reporting this issue.');
                end
                % compilation failed; only a part of the binaries may be available...
                retainid = setdiff(recompileid,failedid);
                eraseid = [];
            end

            % move the mex files into their own directory
            moveid = retainid(cellfun(@exist,o.idents(retainid)) == 3);
            if ~isempty(moveid)
                % some files to be moved
                dest_path = [o.dir filesep 'build-' hlp_hostname filesep];
                % create a new directory
                if ~exist(dest_path,'dir')
                    if ~mkdir(o.dir,['build-' hlp_hostname])
                        error(['unable to create directory ' dest_path]); end
                    % set permissions
                    try
                        fileattrib(dest_path,'+w','a');
                    catch
                        disp(['Note: There are permission problems for the directory ' dest_path]);
                    end
                end
                % create a new env_add.m there
                try
                    filename = [dest_path 'env_add.m'];
                    fid = fopen(filename,'w+');
                    if embed_test
                        % if we had a successful test, we use this to control inclusion of the mex files
                        fprintf(fid,o.test);
                    else
                        % otherwise we check whether any one of the identifiers is recognized by
                        % MATLAB as a mex function
                        fprintf(fid,'any(cellfun(@exist,%s)==3)',hlp_tostring(o.idents(moveid)));
                    end
                    fclose(fid);
                    fileattrib(filename,'+w','a');
                catch
                    disp(['Note: There were write permission problems for the file ' filename]);
                end
                % move the targets over there...
                movefiles = unique(o.idents(moveid));
                for t = 1:length(movefiles)
                    [d,n,x] = fileparts(which(movefiles{t}));
                    movefile([d filesep n x],[dest_path n x]);
                    try
                        movefile([d filesep n x '.md5'],[dest_path n x '.md5']);
                    catch
                    end
                end
                % add the destination path
                addpath(dest_path);
                % and to be entirely sure, CD into that directory to verify that the files are being recognized...
                % (and don't get shadowed by whatever is in the directory below)
                cd(dest_path);
                all_ok = all(strncmp(dest_path,cellfun(@which,o.idents(moveid),'UniformOutput',false),length(dest_path)));
                cd(o.dir);
                % make sure that they are still being found...
                if ~all_ok
                	error('It could not be verified that the MEX file records in %s were successfully updated to their new sub-directories.',o.dir); end
            end

            % move the java class files into their own directory
            infos = dir([o.dir filesep '*.class']);
            movefiles = {infos.name};
            moveid = retainid(cellfun(@exist,o.idents(retainid)) ~= 3);
            if ~isempty(movefiles)
                % some files to be moved
                dest_path = [o.dir filesep 'build-javaclasses' filesep];
                % create a new directory
                if ~exist(dest_path,'dir')
                    if ~mkdir(o.dir,'build-javaclasses')
                        error(['unable to create directory ' dest_path]); end
                    % set permissions
                    try
                        fileattrib(dest_path,'+w','a');
                    catch
                        disp(['Note: There are permission problems for the directory ' dest_path]);
                    end
                end
                % create a new env_add.m there
                try
                    filename = [dest_path 'env_add.m'];
                    fid = fopen(filename,'w+'); fclose(fid);
                    fileattrib(filename,'+w','a');
                catch
                    disp(['Note: There were write permission problems for the file ' filename]);
                end
                % move the targets over there...
                for t = 1:length(movefiles)
                    movefile([o.dir filesep movefiles{t}],[dest_path movefiles{t}]);
                    try
                        movefile([o.dir filesep movefiles{t} '.md5'],[dest_path movefiles{t} '.md5']);
                    catch
                    end
                end
                % add the destination path
                if isdeployed
                    warning off MATLAB:javaclasspath:jarAlreadySpecified; end
                javaaddpath(dest_path);

                % check whether the class is found
                if ~all(cellfun(@exist,o.idents(moveid)) == 8)
                    disp_once('Not all Java binaries in %s could be recognized by MATLAB.',dest_path); end
            end

            if ~isempty(eraseid)
                % some files need to be erased...
                for k=eraseid
                    for t=1:length(targets{k})
                        if exist([o.dir filesep targets{k}{t}])
                            try
                                delete(targets{k}{t});
                            catch
                                disp(['Could not delete broken binary ' binaries{i}{s}]);
                            end
                        end
                    end
                end
            end
        else
            % nothing to recompile
            ok = true;
            if strcmp(o.style,'eager') && ~isempty(o.idents) && o.verbose
                disp_once(['The functions ' format_cellstr(o.idents) ' are already compiled.']); end
        end

        % go back to the old directory
        if hlp_matlab_version < 706
            cd(olddir); end
    catch e
        ok = false; %#ok<NASGU>
        % go back to the old path in case of an error
        if hlp_matlab_version < 706
            cd(olddir); end
        rethrow(e);
    end
end

% store the OK flag in the results
if strcmp(o.style,'lazy') || isdeployed
    results.(tag) = ok; end

if ~ok
    if ~isdeployed
        % regular error summary
        if mexinvolved
            disp_once('\nIn case you need to use a better / fully supported compiler, please have a look at:');
            try
                v=version;
                releasename = v(find(v=='(')+1 : find(v==')')-1);
                if length(releasename) > 3 && releasename(1) == 'R'
                    releasename = releasename(2:end); end
                site = ['http://www.google.com/search?q=matlab+supported+compilers+' releasename];
            catch
                site = 'http://www.google.com/search?q=matlab+supported+compilers';
            end
            disp_once('  <a href="%s">%s</a>\n',site,site);
            if ispc
                if is64bit
                    disp_once('On 64-bit Windows, MATLAB comes with no built-in compiler, so you need to have one installed.');
                else
                    disp_once('On 32-bit Windows, MATLAB supplies a built-in compiler (LLC), which is, however, not very good.');
                end
                disp_once('A good choice is the free Microsoft Visual Studio 2005/2008/2010 Express compiler suite');
                disp_once('together with the Microsoft Platform SDK (6.1 for 2008, 7.1 for 2010) for your Windows Version.');
                disp_once('See also: http://argus-home.coas.oregonstate.edu/forums/development/core-software-development/compiling-64-bit-mex-files');
                disp_once('          http://www.mathworks.com/support/compilers/R2010b/win64.html');
                disp_once('The installation is easier if a professional Intel or Microsoft compiler is used.');
            elseif isunix
                disp_once('On Linux/UNIX, the best choice is usually a supported version of the GCC compiler suite.');
            else
                disp_once('On Mac OS, you need to have a supported version of Xcode/GCC installed.');
            end
        end
        if javainvolved
            disp_once('Please make sure that your system''s java configuration matches the one used by MATLAB (see "ver" command).');
        end
    end
end

% create javac options string from options struct
function opts = javac_options(o)
verbosity = hlp_rewrite(o.verbose,true,'-verbose',false,'');
if ~isempty(o.libdirs)
    cpath = ['-classpath ' sprintf('%s;',o.libdirs{:})];
    cpath(end) = [];
else
    cpath = '';
end
if ~isempty(o.includedirs)
    ipath = ['-sourcepath ' sprintf('%s;',o.includedirs{:})];
    ipath(end) = [];
else
    ipath = '';
end
debugness = hlp_rewrite(o.debug,true,'-g',false,'-g:none');
targetsource = '-target 1.6 -source 1.6';
opts = [sprintf(' %s',verbosity,cpath,ipath,debugness,targetsource) ' '];

% create mex options string from options struct
function opts = mex_options(o)
if ~isempty(o.defines)
    defs = sprintf(' -D%s',o.defines{:});
else
    defs = '';
end
debugness = hlp_rewrite(o.debug,true,'-g',false,'');
if ~isempty(o.includedirs)
    incdirs = sprintf(' -I"%s"',o.includedirs{:});
else
    incdirs = '';
end
if ~isempty(o.libs)
    if ispc
        libs = sprintf(' -l"%s"',o.libs{:});
    else
        libs = sprintf(' -l%s',o.libs{:});
    end
else
    libs = '';
end
if ~isempty(o.libdirs)
    libdirs = sprintf(' -L"%s"',o.libdirs{:});
else
    libdirs = '';
end
verbosity = hlp_rewrite(o.verbose,true,'-v',false,'');
opts = [sprintf(' %s',defs,debugness,incdirs,libs,libdirs,verbosity) ' '];

% format a non-empty cell-string array into a string
function x = format_cellstr(x)
if isempty(x)
    x = '';
else
    x = ['{' sprintf('%s, ',x{1:end-1}) x{end} '}'];
end

% check whether a given identifier is frozen in a ctf archive
function tf = in_ctf(ident) %#ok<DEFNU>
tf = isdeployed && strncmp(ctfroot,which(ident),length(ctfroot));

% normalize a directory path
function dir = path_normalize(dir)
if filesep == '\';
    dir(dir == '/') = filesep;
else
    dir(dir == '\') = filesep;
end
if dir(end) == filesep
    dir = dir(1:end-1); end

% determine if a given file is a mex source file or a java source file
% (and compiles into an identifier that is seen by MATLAB)
function tf = is_primary(filename)
if length(filename)>5 && strcmp(filename(end-4:end),'.java')
    tf = true;
    return;
else
    tf = false;
end
fid = fopen(filename);
if fid ~= -1
    try
        contents = fread(fid);
        tf = ~isempty(strfind(char(contents)','mexFunction')); %#ok<FREAD>
        fclose(fid);
    catch
        fclose(fid);
    end
end

% act like deal, but with a single cell array as input
function varargout = celldeal(argin)
varargout = argin;

% for old MATLABs that can't properly move files...
function movefile(src,dst)
try
    builtin('movefile',src,dst);
catch e
    if any([src dst]=='$') && hlp_matlab_version <= 705
        if ispc
            [errcode,text] = system(sprintf('move ''%s'' ''%s''',src,dst)); %#ok<NASGU>
        else
            [errcode,text] = system(sprintf('mv ''%s'' ''%s''',src,dst)); %#ok<NASGU>
        end
        if errcode
            error('Failed to move %s to %s.',src,dst); end
    else
        rethrow(e);
    end
end

% for old MATLABs that don't handle Java classes on the dynamic path...
function res = exist(obj,type)
if nargin > 1
    res = builtin('exist',obj,type);
    if  ~res && (hlp_matlab_version <= 704) && strcmp(type,'class') && builtin('exist',[obj '.class'],'file')
        res = 8; end
else
    res = builtin('exist',obj);
    if  ~res && (hlp_matlab_version <= 704) && builtin('exist',[obj '.class'])
        res = 8; end
end

% for old MATLABs that don't handle Java classes on the dynamic path...
function res = which(ident)
res = builtin('which',ident);
if ~any(res == filesep)
    if hlp_matlab_version <= 704
        if isempty(res) || ~isempty(strfind(res,'not found'))
            res = builtin('which',[ident '.class']); end
    else
        if ~isempty(strfind(res,'Java'))
            res = builtin('which',[ident '.class']); end
    end
end

repo: REST


text:
function [ abd, rcond, z, info ] = spbco ( abd, lda, n, m )

%*****************************************************************************80
%
%% SPBCO factors a real symmetric positive definite banded matrix.
%
%  Discussion:
%
%    SPBCO also estimates the condition of the matrix.
%
%    If RCOND is not needed, SPBFA is slightly faster.
%
%    To solve A*X = B, follow SPBCO by SPBSL.
%
%    To compute inverse(A)*C, follow SPBCO by SPBSL.
%
%    To compute determinant(A), follow SPBCO by SPBDI.
%
%  Band storage:
%
%    If A is a symmetric positive definite band matrix, the following
%    program segment will set up the input.
%
%      m = (band width above diagonal)
%      do j = 1, n
%        i1 = max (1, j-m)
%        do i = i1, j
%          k = i-j+m+1
%          abd(k,j) = a(i,j)
%        end do
%      end do
%
%    This uses M + 1 rows of A, except for the M by M upper left triangle,
%    which is ignored.
%
%    For example, if the original matrix is
%
%      11 12 13  0  0  0
%      12 22 23 24  0  0
%      13 23 33 34 35  0
%       0 24 34 44 45 46
%       0  0 35 45 55 56
%       0  0  0 46 56 66
%
%    then N = 6, M = 2  and ABD should contain
%
%       *  * 13 24 35 46
%       * 12 23 34 45 56
%      11 22 33 44 55 66
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real ABD(LDA,N), the matrix to be factored.  The columns of the
%    upper triangle are stored in the columns of ABD and the diagonals of
%    the upper triangle are stored in the rows of ABD.
%
%    Input, integer LDA, the leading dimension of the array ABD.
%    M+1 <= LDA is required.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer M, the number of diagonals above the main diagonal.
%
%    Output, real ABD(LDA,N), an upper triangular matrix R, stored in band form,
%    so that A = R'*R.  If INFO /= 0, the factorization is not complete.
%
%    Output, real RCOND, an estimate of the reciprocal condition
%    of A.  For the system A*X = B, relative perturbations in A and B of size
%    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
%    If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate underflows.
%
%    Output, real Z(N), a work vector whose contents are usually
%    unimportant.  If A is singular to working precision, then Z is an
%    approximate null vector in the sense that
%      norm(A*Z) = RCOND * norm(A) * norm(Z).
%    If INFO /= 0, Z is unchanged.
%
%    Output, integer INFO, error flag.
%    0, for normal return.
%    K, signals an error condition.  The leading minor of order K is not
%    positive definite.
%

%
%  Find the norm of A.
%
  for j = 1 : n

    l = min ( j, m+1 );
    mu = max ( m+2-j, 1 );
    z(j) = sasum ( l, abd(mu:mu+l-1,j), 1 );
    k = j - l;
    for i = mu : m
      k = k + 1;
      z(k) = z(k) + abs ( abd(i,j) );
    end

  end

  anorm = max ( z(1:n) );
%
%  Factor.
%
  [ abd, info ] = spbfa ( abd, lda, n, m );

  if ( info ~= 0 )
    return
  end
%
%  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
%
%  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
%
%  The components of E are chosen to cause maximum local
%  growth in the elements of W where R'*W = E.
%
%  The vectors are frequently rescaled to avoid overflow.
%
%  Solve R' * W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  for k = 1 : n

    if ( z(k) ~= 0.0 )
      ek = - abs ( ek ) * r4_sign ( z(k) );
    end

    if ( abd(m+1,k) < abs ( ek - z(k) ) )
      s = abd(m+1,k) / abs ( ek - z(k) );
      z(1:n) = s * z(1:n);
      ek = s * ek;
    end

    wk = ek - z(k);
    wkm = -ek - z(k);
    s = abs ( wk );
    sm = abs ( wkm );
    wk = wk / abd(m+1,k);
    wkm = wkm / abd(m+1,k);
    j2 = min ( k+m, n );
    i = m + 1;

    if ( k+1 <= j2 )

      for j = k+1 : j2
        i = i - 1;
        sm = sm + abs ( z(j) + wkm * abd(i,j) );
        z(j) = z(j) + wk * abd(i,j);
        s = s + abs ( z(j) );
      end

      if ( s < sm )

        t = wkm - wk;
        wk = wkm;
        i = m + 1;

        for j = k+1 : j2
          i = i - 1;
          z(j) = z(j) + t * abd(i,j);
        end

      end

    end

    z(k) = wk;

  end

  z(1:n) = z(1:n) / sasum ( n, z(1:n), 1 );
%
%  Solve R * Y = W.
%
  for k = n : -1 : 1

    if ( abd(m+1,k) < abs ( z(k) ) )
      s = abd(m+1,k) / abs ( z(k) );
      z(1:n) = s * z(1:n);
    end

    z(k) = z(k) / abd(m+1,k);
    lm = min ( k-1, m );
    la = m + 1 - lm;
    lb = k - lm;
    t = -z(k);
    z(lb:lb+lm-1) = z(lb:lb+lm-1) + t * abd(la:la+lm-1,k)';

  end

  z(1:n) = z(1:n) / sasum ( n, z(1:n), 1 );

  ynorm = 1.0;
%
%  Solve R' * V = Y.
%
  for k = 1 : n

    lm = min ( k-1, m );
    la = m + 1 - lm;
    lb = k - lm;

    z(k) = z(k) - sdot ( lm, abd(la:la+lm-1,k), 1, z(lb:lb+lm-1), 1 );

    if ( abd(m+1,k) < abs ( z(k) ) )
      s = abd(m+1,k) / abs ( z(k) );
      z(1:n) = s * z(1:n);
      ynorm = s * ynorm;
    end

    z(k) = z(k) / abd(m+1,k);

  end

  s = 1.0 / sasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;
%
%  Solve R * Z = W.
%
  for k = n : -1 : 1

    if ( abd(m+1,k) < abs ( z(k) ) )
      s = abd(m+1,k) / abs ( z(k) );
      z(1:n) = s * z(1:n);
      ynorm = s * ynorm;
    end

    z(k) = z(k) / abd(m+1,k);
    lm = min ( k-1, m );
    la = m + 1 - lm;
    lb = k - lm;
    t = -z(k);
    z(lb:lb+lm-1) = z(lb:lb+lm-1) + t * abd(la:la+lm-1,k)';

  end
%
%  Make ZNORM = 1.0.
%
  s = 1.0 / sasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%-----------------------------------------------------------------------
% Job saved on 04-Sep-2015 11:55:24 by cfg_util (rev $Rev: 7703 $)
% spm SPM - SPM12 (12.1)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------#
matlabbatch{1}.spm.tools.beamforming.data.val = 1;
matlabbatch{1}.spm.tools.beamforming.data.gradsource = 'inv';
matlabbatch{1}.spm.tools.beamforming.data.space = 'MNI-aligned';
matlabbatch{1}.spm.tools.beamforming.data.overwrite = 0;
matlabbatch{2}.spm.tools.beamforming.sources.BF(1) = cfg_dep('Prepare data: BF.mat file', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{2}.spm.tools.beamforming.sources.reduce_rank = [2 3];
matlabbatch{2}.spm.tools.beamforming.sources.keep3d = 1;
matlabbatch{2}.spm.tools.beamforming.sources.plugin.mesh.orient = 'unoriented';
matlabbatch{2}.spm.tools.beamforming.sources.plugin.mesh.fdownsample = 1;
matlabbatch{2}.spm.tools.beamforming.sources.plugin.mesh.symmetric = 'no';
matlabbatch{2}.spm.tools.beamforming.sources.plugin.mesh.flip = false;
matlabbatch{2}.spm.tools.beamforming.sources.visualise = 1;
matlabbatch{3}.spm.tools.beamforming.features.BF(1) = cfg_dep('Define sources: BF.mat file', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{3}.spm.tools.beamforming.features.whatconditions.all = 1;
matlabbatch{3}.spm.tools.beamforming.features.woi = [-100 0
                                                     100 200];
matlabbatch{3}.spm.tools.beamforming.features.modality = {'MEG'};
matlabbatch{3}.spm.tools.beamforming.features.fuse = 'no';
matlabbatch{3}.spm.tools.beamforming.features.plugin.vbfa.nl = 5;
matlabbatch{3}.spm.tools.beamforming.features.plugin.vbfa.nem = 50;
matlabbatch{3}.spm.tools.beamforming.features.regularisation.manual.lambda = 0;
matlabbatch{3}.spm.tools.beamforming.features.bootstrap = false;
matlabbatch{4}.spm.tools.beamforming.inverse.BF(1) = cfg_dep('Covariance features: BF.mat file', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{4}.spm.tools.beamforming.inverse.plugin.champagne.nem = 100;
matlabbatch{4}.spm.tools.beamforming.inverse.plugin.champagne.vcs = 0;
matlabbatch{4}.spm.tools.beamforming.inverse.plugin.champagne.nupd = 0;
matlabbatch{5}.spm.tools.beamforming.output.BF(1) = cfg_dep('Inverse solution: BF.mat file', substruct('.','val', '{}',{4}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{5}.spm.tools.beamforming.output.plugin.image_power.whatconditions.all = 1;
matlabbatch{5}.spm.tools.beamforming.output.plugin.image_power.sametrials = false;
matlabbatch{5}.spm.tools.beamforming.output.plugin.image_power.woi = [100 200];
matlabbatch{5}.spm.tools.beamforming.output.plugin.image_power.foi = [0 Inf];
matlabbatch{5}.spm.tools.beamforming.output.plugin.image_power.contrast = 1;
matlabbatch{5}.spm.tools.beamforming.output.plugin.image_power.result = 'bycondition';
matlabbatch{5}.spm.tools.beamforming.output.plugin.image_power.scale = 0;
matlabbatch{5}.spm.tools.beamforming.output.plugin.image_power.powermethod = 'trace';
matlabbatch{5}.spm.tools.beamforming.output.plugin.image_power.modality = 'MEG';
matlabbatch{6}.spm.tools.beamforming.write.BF(1) = cfg_dep('Output: BF.mat file', substruct('.','val', '{}',{5}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{6}.spm.tools.beamforming.write.plugin.gifti.normalise = 'all';
matlabbatch{6}.spm.tools.beamforming.write.plugin.gifti.space = 'mni';
matlabbatch{6}.spm.tools.beamforming.write.plugin.gifti.visualise = 2;

author: spm


text:
function A_com = extract_patch(A,dims,patch_size,padding)

% Extractx a patch of size patch_size centered around the centroid of each
% component.
% INPUTS:
% A:            2d matrix of spatial components
% dims:         dimensions of FOV
% patch_size:   dimensions of patch
% padding:      if true components remain centered and are zero padded,
%               otherwise they are shifted (default: true)

% OUTPUT:
% A_com:            Nd matrix of patches

if ~exist('padding','var'); padding = true; end

nd = length(dims);
if nd == 2; dims(3) = 1; patch_size(3) = 1; end
K = size(A,2);
A = A/spdiags(sqrt(sum(A.^2,1))'+eps,0,K,K);      % normalize to sum 1 for each compoennt
cm = com(A,dims(1),dims(2),dims(3));
xx = -ceil(patch_size(1)/2-1):floor(patch_size(1)/2);
yy = -ceil(patch_size(2)/2-1):floor(patch_size(2)/2);
zz = -ceil(patch_size(3)/2-1):floor(patch_size(3)/2);
A_com = zeros([patch_size,K]);

for i = 1:K
    int_x = round(cm(i,1)) + xx;
    pad_pre_x = 0;
    pad_pre_y = 0;
    pad_post_x = 0;
    pad_post_y = 0;
    if int_x(1)<1
        if padding
            pad_pre_x = 1 - int_x(1);
            int_x = 1:int_x(end);
        else
            int_x= int_x + 1 - int_x(1);
        end
    end
    if int_x(end)>dims(1)
        if padding
            pad_post_x = int_x(end) - dims(1);
            int_x = int_x(1):dims(1);
        else
            int_x = int_x - (int_x(end)-dims(1));
        end
    end
    int_y = round(cm(i,2)) + yy;
    if int_y(1)<1
        if padding
            pad_pre_y = 1 - int_y(1);
            int_y = 1:int_y(end);
        else
            int_y = int_y + 1 - int_y(1);
        end
    end
    if int_y(end)>dims(2)
        if padding
            pad_post_y = int_y(end) - dims(2);
            int_y = int_y(1):dims(2);
        else
            int_y = int_y - (int_y(end)-dims(2));
        end
    end
    if nd == 3
        int_z = round(cm(i,3)) + zz;
        if int_z(1)<1
            int_z = int_z + 1 - int_z(1);
        end
        if int_z(end)>dims(3)
            int_z = int_z - (int_z(end)-dims(3));
        end
    else
        int_z = 1;
    end
    A_temp = reshape(full(A(:,i)),dims);
    A_temp = A_temp(int_x,int_y,int_z);
    if padding
        A_temp = padarray(A_temp,[pad_pre_x,pad_pre_y],0,'pre');
        A_temp = padarray(A_temp,[pad_post_x,pad_post_y],0,'post');
    end
    A_com(:,:,i) = A_temp;
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function Prepare_data_combined_all()

    % Read in in-the-wild data and Multi-PIE data and concatenate them
    root = '../prepared_data/';
    Collect_combined_data(root, '0.25');
    Collect_combined_data(root, '0.35');
    Collect_combined_data(root, '0.5');
    Collect_combined_data(root, '1');
end

function Collect_combined_data(root, scale)

    wild_locs = dir(sprintf('%s/wild_%s*.mat', root, scale));
    mpie_locs = dir(sprintf('%s/mpie_%s*.mat', root, scale));

    % For a particular loaded m_pie check if an appropriate in-the-wild
    % exists
    for i=1:numel(mpie_locs)

        load([root, mpie_locs(i).name]);

        imgs_used_pie = actual_imgs_used;
        samples_pie = all_images;
        centres_pie = centres;
        landmark_locations_pie = landmark_locations;

        wild_to_use = -1;
        for j=1:numel(wild_locs)
            load([root, wild_locs(j).name], 'centres');
            if(isequal(centres_pie, centres))
                wild_to_use = j;
            end
        end

        % Reset the centres
        centres = centres_pie;

        if(wild_to_use ~= -1)
            load([root, wild_locs(wild_to_use).name]);
            actual_imgs_used = cat(1, actual_imgs_used, imgs_used_pie);
            all_images = cat(1, all_images, samples_pie);
            landmark_locations = cat(1, landmark_locations, landmark_locations_pie);
        end

        save(sprintf('%s/combined_%s_%d.mat', root, scale, i), 'actual_imgs_used', 'all_images', 'centres', 'landmark_locations', 'training_scale', 'visiIndex');
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
%% OPTI Toolbox AMPL Examples
%
% This file loads a number of supplied examples and shows how to solve
% them using the OPTI Toolbox. You should read and complete BasicUsage.m
% BEFORE running the below examples.
%
% The underlying parser uses Netlib's AMPL Solver Library (ASL) Interface.
%
%   Copyright (C) 2014 Jonathan Currie (IPL)

% There is also a page on the Wiki which supplements this example:
web('https://www.inverseproblem.co.nz/OPTI/index.php/File/AMPL');

%% Loading an AMPL Problem
% OPTI Toolbox is supplied with a number of example AMPL problems ranging
% from LP to NLP including continuous and discrete problems. To load an
% AMPL model it must be in .NL format, see the supplied user's guide page
% on AMPL Interfacing in order to generate the .NL file from your model. To
% load an AMPL problem, simply use the command below. Returned will be a
% optiprob structure containing the data in the file.
clc
prob = amplRead('diet.nl')

%% Example 1 - Solving a Loaded AMPL Problem
% Solving a loaded AMPL problem is simple, just pass it to the opti
% constructor and call solve. You will note that the command amplRead will
% automatically determine the problem type (LP, QP, NLP, etc) based on
% properties supplied by the ASL.
clc
Opt = opti(prob)

x = solve(Opt)

%% Loading an AMPL Nonlinear Program
% The AMPL format can easily specify nonlinear problems unlike the LP or
% MPS formats.
clc
prob = amplRead('ch3.nl')

%% Example 2 - Solving a Nonlinear Problem from AMPL
% In order to solve a nonlinear problem the model .nl file
% remains open after the amplRead function is called to allow callbacks for
% the objective, gradient, etc. Calling solve() will automatically close
% the interface for linear and quadratic problems, however you must manually
% call close for nonlinear problems. Alternatively it will be automatically
% called once Matlab is closed to clean up ASL memory.
clc
Opt = opti(prob)

x = solve(Opt)

asl('close')

%% Loading an Integer Problem
% AMPL will reorder the variables from the original model based on internal
% rules (see the supplied pdfs). Based on this the ASL indicates which
% variables have binary or integer constraints, which are passed through
% the amplRead interface:
clc
prob = amplRead('multmip1.nl')

%% Example 3 - Solving an AMPL Integer Problem
% As the optiprob structure can fully define an integer problem, no change
% is required to build and solve the MIP
clc
Opt = opti(prob)

[x,fval] = solve(Opt);
fval

%% A Note On Hessians
% As of OPTI v1.58 the Hessian of an NLP will automatically be added to the
% problem structure. However note it is the Hessian of the Lagrangian, with
% the calling form:
%
%   H = hessian(x,sigma,lambda)
%
% This means the Hessian includes second derivatives of both the objective
% and the constraints. Sigma is a scalar scaling factor on the objective
% derivatives only, and lambda are the lagrange multipliers at the current
% point x (scaling the constraint derivatives).

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [chanMapName] = make_eMouseChannelMap_3B_short(fpath,NchanTOT)
% create a channel Map file for simulated data on a section of
% an imec 3B probe to use with eMouse
% essentially identical to the 3A version

% total number of channels = 385 (in real 3A, 384 channels + digital)
chanMap = (1:NchanTOT)';

% channels to ignore in analysis and when adding data
% in real 3A data, include the reference channels and the digital channel
% replicate the refChans that are within range of the short probe to
% preserve the geometry for channel to channel correlation in noise
% generation

allRef = [37, 76, 113, 152, 189, 228, 265, 304, 341, 380];
refChan = allRef( find(allRef < max(NchanTOT)) );

connected = true(NchanTOT,1);
connected(refChan) = 0;

% copy the coordinates from MP chanmap for 3A.
halfChan = floor(NchanTOT/2);

xcoords = zeros(NchanTOT,1,'double');
ycoords = zeros(NchanTOT,1,'double');

xcoords(1:4:NchanTOT) = 43;
xcoords(2:4:NchanTOT) = 11;
xcoords(3:4:NchanTOT) = 59;
xcoords(4:4:NchanTOT) = 27;

ycoords(1:2:NchanTOT) = 20*(1:halfChan);
ycoords(2:2:NchanTOT) = 20*(1:halfChan);

% Often, multi-shank probes or tetrodes will be organized into groups of
% channels that cannot possibly share spikes with the rest of the probe. This helps
% the algorithm discard noisy templates shared across groups. In
% this case, we set kcoords to indicate which group the channel belongs to.
% In our case all channels are on the same shank in a single group so we
% assign them all to group 1.

kcoords = ones(NchanTOT,1);

% at this point in Kilosort we do data = data(connected, :), ycoords =
% ycoords(connected), xcoords = xcoords(connected) and kcoords =
% kcoords(connected) and no more channel map information is needed (in particular
% no "adjacency graphs" like in KlustaKwik).
% Now we can save our channel map for the eMouse.

% would be good to also save the sampling frequency here
fs = 30000;

chanMapName = sprintf('chanMap_3A_%dsites.mat', NchanTOT);

save(fullfile(fpath, chanMapName), 'chanMap', 'connected', 'xcoords', 'ycoords', 'kcoords', 'fs', 'NchanTOT' )

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function mapping = calc_projection_from_library(endmembers, reduced_dim)
if 1 % Make the same size for all endmembers and the mean lie at the center
    num_spectra_per_endm = zeros(1,length(endmembers));
    for i = 1:length(endmembers)
        num_spectra_per_endm(i) = size(endmembers{i},1);
    end
    num_spectra = round(mean(num_spectra_per_endm));
    spectra_train = [];
    for i = 1:length(endmembers)
        inds = linspace(1, num_spectra_per_endm(i), num_spectra);
        inds = round(inds);
        spectra_train = cat(1, spectra_train, endmembers{i}(inds,:));
    end
else % The size difference is usually not big for spectra libraries.
    spectra_train = cell2mat(endmembers');
end
[~,mapping] = pca(spectra_train, reduced_dim);

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [dat] = read_yokogawa_data_new(filename, hdr, begsample, endsample, chanindx)

% READ_YOKAGAWA_DATA_NEW reads continuous, epoched or averaged MEG data
% that has been generated by the Yokogawa MEG system and software
% and allows that data to be used in combination with FieldTrip.
%
% Use as
%   [dat] = read_yokogawa_data_new(filename, hdr, begsample, endsample, chanindx)
%
% This is a wrapper function around the function
% getYkgwData
%
% See also READ_YOKOGAWA_HEADER_NEW, READ_YOKOGAWA_EVENT

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

if ~ft_hastoolbox('yokogawa_meg_reader')
    ft_error('cannot determine whether Yokogawa toolbox is present');
end

% hdr = read_yokogawa_header(filename);
hdr = hdr.orig; % use the original Yokogawa header, not the FieldTrip header

% default is to select all channels
if nargin<5
  chanindx = 1:hdr.channel_count;
end

handles = definehandles;

switch hdr.acq_type
  case handles.AcqTypeEvokedAve
    % dat is returned as double
    start_sample  = begsample - 1; % samples start at 0
    sample_length = endsample - begsample + 1;
    epoch_count   = 1;
    start_epoch   = 0;
    dat = getYkgwData(filename, start_sample, sample_length);

  case handles.AcqTypeContinuousRaw
    % dat is returned as double
    start_sample  = begsample - 1; % samples start at 0
    sample_length = endsample - begsample + 1;
    epoch_count   = 1;
    start_epoch   = 0;
    dat = getYkgwData(filename, start_sample, sample_length);

  case handles.AcqTypeEvokedRaw
    % dat is returned as double
    begtrial = ceil(begsample/hdr.sample_count);
    endtrial = ceil(endsample/hdr.sample_count);
    if begtrial<1
      ft_error('cannot read before the begin of the file');
    elseif endtrial>hdr.actual_epoch_count
      ft_error('cannot read beyond the end of the file');
    end
    epoch_count = endtrial-begtrial+1;
    start_epoch = begtrial-1;
    % read all the necessary trials that contain the desired samples
    dat = getYkgwData(filename, start_epoch, epoch_count);
    if size(dat,2)~=epoch_count*hdr.sample_count
      ft_error('could not read all epochs');
    end
    rawbegsample = begsample - (begtrial-1)*hdr.sample_count;
    rawendsample = endsample - (begtrial-1)*hdr.sample_count;
    sample_length = rawendsample - rawbegsample + 1;
    % select the desired samples from the complete trials
    dat = dat(:,rawbegsample:rawendsample);

  otherwise
    ft_error('unknown data type');
end

if size(dat,1)~=hdr.channel_count
  ft_error('could not read all channels');
elseif size(dat,2)~=(endsample-begsample+1)
  ft_error('could not read all samples');
end

% select only the desired channels
dat = dat(chanindx,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this defines some usefull constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handles = definehandles
handles.output = [];
handles.sqd_load_flag = false;
handles.mri_load_flag = false;
handles.NullChannel         = 0;
handles.MagnetoMeter        = 1;
handles.AxialGradioMeter    = 2;
handles.PlannerGradioMeter  = 3;
handles.RefferenceChannelMark = hex2dec('0100');
handles.RefferenceMagnetoMeter       = bitor( handles.RefferenceChannelMark, handles.MagnetoMeter );
handles.RefferenceAxialGradioMeter   = bitor( handles.RefferenceChannelMark, handles.AxialGradioMeter );
handles.RefferencePlannerGradioMeter = bitor( handles.RefferenceChannelMark, handles.PlannerGradioMeter );
handles.TriggerChannel      = -1;
handles.EegChannel          = -2;
handles.EcgChannel          = -3;
handles.EtcChannel          = -4;
handles.NonMegChannelNameLength = 32;
handles.DefaultMagnetometerSize       = (4.0/1000.0);       % Square of 4.0mm in length
handles.DefaultAxialGradioMeterSize   = (15.5/1000.0);      % Circle of 15.5mm in diameter
handles.DefaultPlannerGradioMeterSize = (12.0/1000.0);      % Square of 12.0mm in length
handles.AcqTypeContinuousRaw = 1;
handles.AcqTypeEvokedAve     = 2;
handles.AcqTypeEvokedRaw     = 3;
handles.sqd = [];
handles.sqd.selected_start  = [];
handles.sqd.selected_end    = [];
handles.sqd.axialgradiometer_ch_no      = [];
handles.sqd.axialgradiometer_ch_info    = [];
handles.sqd.axialgradiometer_data       = [];
handles.sqd.plannergradiometer_ch_no    = [];
handles.sqd.plannergradiometer_ch_info  = [];
handles.sqd.plannergradiometer_data     = [];
handles.sqd.eegchannel_ch_no   = [];
handles.sqd.eegchannel_data    = [];
handles.sqd.nullchannel_ch_no   = [];
handles.sqd.nullchannel_data    = [];
handles.sqd.selected_time       = [];
handles.sqd.sample_rate         = [];
handles.sqd.sample_count        = [];
handles.sqd.pretrigger_length   = [];
handles.sqd.matching_info   = [];
handles.sqd.source_info     = [];
handles.sqd.mri_info        = [];
handles.mri                 = [];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip