text:
function [Model cfg] = est_fitMVAR(varargin)
%
% Fit (adaptive) multivariate autoregressive model to EEG data. See [1] for
% details on VAR model fitting and implementations.
%
%
% Output:
%
%   Model structure with
%       .Model          (numvars x coeffs) matrix of VAR coefficients
%       .PE             (numvars x coeffs) prediction error (noise covariance) coefficients
%       .algorithm      string denoting algorithm used for estimation
%       .modelclass     string denoting model class (here, 'mvar')
%
% See Also: pop_est_fitMVAR(), pop_pre_prepData()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Chapters 3,6.
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift/
%
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

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

% pre-cache some frequently used sub-arguments
persistent subargs;
if isempty(subargs)
    subargs.algos = hlp_getMVARalgorithms('defaultNameOnly');
    subargs.all_algos = hlp_getMVARalgorithms;
    subargs.sift_domain = hlp_buildMVARHelpText;
end

verb = arg_extract(varargin,{'verb','VerbosityLevel'},[],2);

haveSigProc = hlp_isToolboxInstalled('Signal Processing Toolbox');
if haveSigProc
    taperFcns = {'rectwin','hamming','hann','bartlett','barthannwin',   ...
                 'blackman','blackmanharris','bohmanwin','chebwin',     ...
                 'flattopwin','gausswin','kaiser','nuttallwin',         ...
                 'parzenwin','taylorwin','tukeywin','triang'};
else
    taperFcns = {'rectwin'};
end

g = arg_define([0 1],varargin, ...
    arg_norep({'EEG','ALLEEG'},mandatory,[],'EEGLAB dataset'), ...
    arg_subswitch({'algorithm','Algorithm'},subargs.algos,subargs.all_algos,{'Select a model fitting algorithm.',subargs.sift_domain},'cat','Modeling Parameters','suppress',{'ModelOrder','OrderSelector','InitialState'}), ...
    arg({'morder','ModelOrder','modelOrder'},10,[1 Inf],'VAR model order.','cat','Modeling Parameters'), ...
    arg_nogui({'winStartIdx','WindowStartIndices'},[],[],'Starting indices for windows. This is a vector of sample points (start of windows) at which to estimate windowed VAR model','cat','Modeling Parameters'), ...
    arg({'winlen','WindowLength'},0.5,[eps Inf],'Sliding window length (sec)','cat','Modeling Parameters'), ...
    arg({'winstep','WindowStepSize'},0.03,[eps Inf],'Window step size (sec)','cat','Modeling Parameters'), ...
    arg({'taperfcn','TaperFunction'},'rectwin',taperFcns,{'Data tapering (windowing) function', ...
                                                            sprintf(['\n' ...
                                                                     'Each data segment (e.g. all data in a sliding window) will be multipled by the data taper to smooth endpoints towards zero.', ...
                                                                     '\n' ...
                                                                     'Available windows are:\n' ...
                                                                          'rectwin \t - Rectangular window. This is equivalent to no taper.\n', ...
                                                                          'hamming \t - Hamming window.\n', ...
                                                                          'hann \t - Hann window.\n', ...
                                                                          'bartlett \t - Bartlett window.\n', ...
                                                                          'barthannwin \t - Modified Bartlett-Hanning window.\n', ...
                                                                          'blackman \t - Blackman window.\n', ...
                                                                          'blackmanharris- Minimum 4-term Blackman-Harris window.\n', ...
                                                                          'bohmanwin \t - Bohman window.\n', ...
                                                                          'chebwin \t - Chebyshev window.\n', ...
                                                                          'flattopwin \t - Flat Top window.\n', ...
                                                                          'gausswin \t - Gaussian window.\n', ...
                                                                          'kaiser \t - Kaiser window.\n', ...
                                                                          'nuttallwin \t - Nuttall defined minimum 4-term Blackman-Harris window.\n', ...
                                                                          'parzenwin \t - Parzen (de la Valle-Poussin) window.\n', ...
                                                                          'taylorwin \t - Taylor window.\n', ...
                                                                          'tukeywin \t - Tukey window.\n', ...
                                                                          'triang \t - Triangular window.\n' ...
                                                                        ])},'cat','Modeling Parameters'), ...
    arg({'epochTimeLims','EpochTimeLimits'},[],[],'Sub-epoch time limits (sec). This is relative to event time (e.g. [-1 2]). Default is the full epoch time range','cat','Modeling Parameters'), ...
    arg({'prctWinToSample','WindowSamplePercent'},100,[1 100],'Percent of windows to sample','cat','Modeling Parameters'), ...
    arg_subtoggle({'normalize','NormalizeData'},[],@pre_normData,'Z-normalize data within windows. Note this is not recommended for short windows','cat','Window Preprocessing'), ...
    arg_subtoggle({'detrend','Detrend'},{}, ...
    {arg({'method','DetrendingMethod'},'constant',{'linear','constant'},{'Detrend data within each window.', ...
    sprintf(['\n' ...
    'Linear: removes the least-squares fit of a straight line.\n' ...
    'Constant: removes the mean from each trial (centering)' ...
    ]) ...
    } ...
    )},'Detrend or center each time window','cat','Window Preprocessing'), ...
    arg({'timer','Timer'},false,[],'Activate timer. Times are stored in EEG.CAT.Model.timeelapsed'), ...
    arg({'setArgDirectMode','SetArgDirectMode'},true,[],'Set arg_direct mode to true. Can improve speed when number of windows is large. Disable if calling this function repeatedly in a tight loop.'), ...
    arg({'verb','VerbosityLevel'},verb,{int32(0) int32(1) int32(2)},'Verbosity level. 0 = no output, 1 = text, 2 = graphical') ...
    );

% commit EEG variable to workspace
[data g] = hlp_splitstruct(g,{'EEG'});
arg_toworkspace(data);
clear data;

% do some error-checking
if isempty(g.epochTimeLims)
    g.epochTimeLims = [EEG.xmin EEG.xmax]; end
if ~(all(g.epochTimeLims>=EEG.xmin) && all(g.epochTimeLims<=EEG.xmax))
    error('Epoch time limits must be within the range [%0.3g %0.3g]',EEG.xmin,EEG.xmax); end
if isempty(g.morder) || length(g.morder)>1
    error('invalid entry for field ''morder.'' Make sure the model order is a single integer.'); end

if nargout > 1, cfg = g; end

%% do some adjustments to parameters

% ensure we are using the right model order
g.algorithm.morder = g.morder;

if rem(g.winstep,1/EEG.srate)
    if g.winstep<1/EEG.srate
        g.winstep = 1/EEG.srate;
    else
        % adjust step size to nearest multiple of sampling interval
        g.winstep = g.winstep-rem(g.winstep,1/EEG.srate);
    end

    if g.verb,
        fprintf('Adjusting window step size to nearest multiple of sampling interval\n');
        fprintf('\tstep size is now %0.5g sec\n',g.winstep);
    end
end

if rem(g.winlen,1/EEG.srate)
    if g.winlen<1/EEG.srate
        g.winlen = 1/EEG.srate;
    else
        % adjust window length to nearest multiple of sampling interval
        g.winlen = g.winlen-rem(g.winlen,1/EEG.srate);
    end

    if g.verb,
        fprintf('Adjusting window length to nearest multiple of sampling interval\n');
        fprintf('\twindow length is now %0.5g sec\n',g.winlen);
    end
end
if g.winlen > EEG.xmax-EEG.xmin
    g.winlen = EEG.xmax-EEG.xmin;
    if g.verb,
        fprintf('Window length exceeds epoch length. Adjusting window length to match epoch length\n');
        fprintf('\twindow length is now %0.5g sec\n',g.winlen);
    end
end
tidx = getindex(EEG.CAT.times,g.epochTimeLims*1000);
if ~all(isequal(EEG.CAT.times(tidx),g.epochTimeLims*1000))

    g.epochTimeLims = EEG.CAT.times(tidx)/1000;

    if g.verb
        fprintf('Adjusting epoch time limits to match sampling interval\n');
        fprintf('\tepoch limits are now [%0.5g, %0.5g] sec\n',g.epochTimeLims(1),g.epochTimeLims(2));
    end
end

winLenPnts  = round(g.winlen*EEG.srate); % window size in points
winStepPnts = round(g.winstep*EEG.srate);

if isempty(g.winStartIdx)
    % starting point of each window (points)
    g.winStartIdx  = tidx(1):winStepPnts:(tidx(2)-winLenPnts)+1;
    %g.winStartIdx  =  round((double(g.epochTimeLims(1):g.winstep:g.epochTimeLims(2)-g.winlen)*EEG.srate)+1;
end

if g.prctWinToSample<100
    % randomly select percentage of windows to work with
    randwin = randperm(length(g.winStartIdx));
    randwin = sort(randwin(1:ceil(length(g.winStartIdx)*g.prctWinToSample/100)));
    g.winStartIdx = g.winStartIdx(randwin);
end

numWins   = length(g.winStartIdx);

% construct the data taper
g.taper = window_func(g.taperfcn,winLenPnts).';

%% Prepare data for model fitting

% initialize results arrays
[AR PE RC mu th lambdaOpt]  = deal(cell(1,numWins));

if g.verb==2
    waitbarTitle = sprintf('Fitting VAR[%d] model (%s)...', ...
                        g.morder, ...
                        num2str(g.algorithm.arg_selection));

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

if g.detrend.arg_selection
    % detrend each window separately
    if g.verb, fprintf('%s detrending each window...\n',firstcaps(g.detrend.method)); end
    for ch=1:EEG.CAT.nbchan
        EEG.CAT.srcdata(ch,:,:) = locdetrend_siftmod(squeeze(EEG.CAT.srcdata(ch,:,:)), ...
            EEG.srate,[g.winlen g.winstep],g.detrend.method);
    end
    if g.verb, fprintf('done.\n'); end
end

if g.timer
    timeElapsed = nan(1,numWins);
else
    timeElapsed = [];
end

%% Main loop: fit MVAR models to each window

algFcnName = hlp_nanocache('algos',10,@hlp_getMVARalgorithms,'mfileNameOnly',g.algorithm.arg_selection);

if g.setArgDirectMode && ~strcmp(algFcnName,'mvar_glADMM') % not necessary for glADMM
    % set the arg_direct flag
    % to improve speed
   g = arg_setdirect(g,true);
end

for t=1:numWins

    if g.timer, tic; end

    % get data for current window
    data = EEG.CAT.srcdata(:,g.winStartIdx(t):g.winStartIdx(t)+winLenPnts-1,:);
    if g.normalize.arg_selection
        % normalize the data window
        data = pre_normData('data',data,'method',g.normalize.method,'verb',0,'arg_direct',true);
    end

    % execute the model-fitting algorithm
    switch nargout(algFcnName)
        case 2
            [AR{t} PE{t}] = feval(algFcnName, ...
                'data',bsxfun(@times,g.taper,data), ...
                g.algorithm,'arg_direct',true);
        case 3
            [AR{t} PE{t} argsout] = feval(algFcnName, ...
                                   'data',bsxfun(@times,g.taper,data),...
                                    g.algorithm,'arg_direct',true);
            if isstruct(argsout)
                % store contents of argsout fields in cell array at index t
                % e.g. fieldname{t} = argsout.(fieldname)
                fnames = fieldnames(argsout);
                for k=1:length(fnames)
                    eval([fnames{k} '{t}=argsout.(''' fnames{k} ''');']);
                end
            end
        otherwise
            error('SIFT:est_fitMVAR:badAlgArgs','%s must output either 2 or 3 arguments',algFcnName);
    end

    if g.verb==2
        drawnow;
        % graphical waitbar
        cancel = multiWaitbar(waitbarTitle,t/numWins);
        if cancel
            if strcmpi('yes',questdlg2( ...
                            'Are you sure you want to cancel?', ...
                            'Model Fitting','Yes','No','No'));
                Model = [];
                multiWaitbar(waitbarTitle,'Close');
                return;
            else
                multiWaitbar(waitbarTitle,'ResetCancel',true);
            end
        end
    end

    if g.timer, timeElapsed(t) = toc; end
end

%% Do some cleanup
clear('-regexp','mvar_*')
if g.verb==2
    multiWaitbar(waitbarTitle,'Close');
end

%% Construct Model object
Model = hlp_sift_emptymodel;

Model.AR = AR;
Model.PE = PE;
Model.RC = RC;
Model.mu = mu;
Model.th = th;
Model.lambdaOpt = lambdaOpt;
Model.winStartTimes = (g.winStartIdx-1)/EEG.srate;
Model.morder        = g.morder;
Model.winstep       = g.winstep;
Model.winlen        = g.winlen;
Model.algorithm     = g.algorithm.arg_selection;
Model.modelclass    = 'mvar';
Model.timeelapsed   = timeElapsed;
Model.normalize     = g.normalize;
Model.modelapproach = 'Segmentation VAR';
Model.taperFcn      = g.taperfcn;

switch lower(g.algorithm.arg_selection)
    case 'group lasso dal/scsa'
        %     Model.ww = ww;
        Model.lambda = g.algorithm.dal_args.lambda;
    case 'group lasso (admm)'
        Model.lambda = g.algorithm.admm_args.lambda;
        Model.rho    = g.algorithm.admm_args.rho;
        Model.alpha  = g.algorithm.admm_args.alpha;
end


text:
function estimateCombinedRF_batchHN(pathRF,nSplit,nBoot,testSplit,testCost,seed,matlabPATH)

startpath = pwd;

% INITIALIZATIONS
cd(pathRF)
load('clinicalPARAM') % parameters gets out of there
training = load('training'); training = struct2cell(training); training = training{1};
nameOutcomes = fieldnames(training.outcomes); nOutcomes = numel(nameOutcomes);
fSetNames = fieldnames(training.textures.(nameOutcomes{1})); nFset = numel(fSetNames);
comb = cell(1,nOutcomes);
for o = 1:nOutcomes
    comb{o} = {parameters.clinical.(nameOutcomes{o})};
end
time = 60;
mkdir('batchLog_RF'), cd('batchLog_RF'), pathBatch = pwd;

% PRODUCE BATCH COMPUTATIONS FOR RF ESTIMATION
batch = 0;
save('workspace','pathRF','nameOutcomes','fSetNames','comb','testCost','nSplit','testSplit','nBoot','seed'), pause(2);
for o = 1:nOutcomes
    for f = 1:nFset
        batch = batch + 1;
        nameScript = ['batch',num2str(batch),'_script.m'];
        fid = fopen(nameScript,'w');
        fprintf(fid,'load(''workspace'')\n');
        fprintf(fid,['[indClinic,bestCost] = estimateCombinedRF(pathRF,nameOutcomes{',num2str(o),'},fSetNames{',num2str(f),'},comb{',num2str(o),'},testCost,nSplit,testSplit,nBoot,seed);\n']);
        fprintf(fid,['save(''result_batch',num2str(batch),''',''indClinic'',''bestCost'')\n']);
        fprintf(fid,['system(''touch batch',num2str(batch),'_end'');\n']);
        fprintf(fid,'clear all');
        fclose(fid);
        system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
    end
end

% WAITING LOOP
waitBatch(pathBatch,time,nOutcomes*nFset)
delete('workspace.mat')

% GROUP RESULTS IN THE 'training' STRUCT
cd(pathBatch)
batch = 0;
for o = 1:nOutcomes
    for f = 1:nFset
        batch = batch + 1;
        load(['result_batch',num2str(batch)]) % Variables 'indClinic' and 'bestCost' gets out of there
        training.clinical.bestAdd.(nameOutcomes{o}).(fSetNames{f}) = indClinic;
        training.cost.(nameOutcomes{o}).(fSetNames{f}) = bestCost;
        delete(['result_batch',num2str(batch),'.mat'])
        clear indClinic bestCost
    end
end
cd(pathRF), save('training','training')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [edges,didSucceed]=graph4VertexDeg(d)
%%GRAPH4VERTEXDEG Given a list of the degrees of vertices in a
%           non-directional graph, get a set of edges representing a graph
%           that satisfies the degree profile, or indicate if no such graph
%           can exist. Note that a set of edges for a degree profile is not
%           unique. This just returns one set.
%
%INPUTS: d An nX1 or 1Xn vector of the degrees of the n vertices. These are
%          integer values >=0. They represent the number of edges touching
%          the nodes.
%
%OUTPUTS: edges A 2XnumEdge set of edges that satisfies the degree profile.
%               An empty matrix is returned if all edges have zero degree
%               or if no set ofvedges exists that can result in the given
%               degree profile. An edge [i;j] goes between vertex i and
%               vertex j and adds to the degree count for both vertices.
%    didSucceed This is true if an edges exists that can satisfy d.
%               Otherwise, this is false and edges is an empty matrix.
%
%This function implements Algorithm H in Chapter 7 of [1].
%
%EXAMPLE:
% [edges,didSucceed]=graph4VertexDeg([2;1;6;3;2;1;1])
%The algorithm in the example will succeed and the edge set is
%edges=[7,6,2,5,5,1,1,4;
%       3,3,3,4,3,4,3,3];
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4A: Combinatorial
%    Algorithms, Part I, Boston: Addison-Wesley, 2011.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

[d,idx]=sort([d(:);0],'descend');

%Allocate space
maxEdges=ceil(sum(d)/2);
edges=zeros(2,maxEdges);
c=zeros(d(1),1);

numEdges=0;

%Step H1
k=d(1);
j=0;
while(k>0)
    j=j+1;
    while(k>d(j+1))
        c(k)=j;
        k=k-1;
    end
end
%If all of the d's are zero.
if(j==0)
    edges=[];
    didSucceed=true;
    return;
end

while(1)
    %Step H2, find n.
    n=c(1);
    if(n==0)
        %Shrink to fit and restore the original indexation.
        edges=idx(edges(:,1:numEdges));
        %Restore the original indexation.
        didSucceed=true;
        return;
    elseif(d(1)>=n)
        edges=[];
        didSucceed=false;
        return;
    end

    %Step H3
    i=1;
    t=d(1);
    r=c(t);
    j=d(n);

    %Step H4, generate a new edge.
    while(1)
        c(j)=c(j)-1;
        m=c(t);
        numEdges=numEdges+1;
        edges(:,numEdges)=[n;m];
        d(m)=d(m)-1;
        c(t)=m-1;
        j=j-1;
        if(j==0)
            break;
        elseif(m==i)
            i=r+1;
            t=d(i);
            r=c(t);
        end
    end
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [pC] = spm_dcm_symm(pV,pE)
% locks ECD orientations by introducing prior correlations
% FORMAT [pC] = spm_dcm_symm(pV,pE)
%__________________________________________________________________________
%
% pE   - prior expectation
% pV   - prior variance
% pC   - prior covariance
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_symm.m 5376 2013-04-02 09:59:01Z karl $

% Distance between homolgous sources (16mm)
%--------------------------------------------------------------------------
V     = 16;

% symmetry constraints (based on Euclidean distance from mirror image)
%==========================================================================

% diagonalise feilds
%--------------------------------------------------------------------------
feilds = fieldnames(pV);
for  i = 1:length(feilds)
    pF = getfield(pV,feilds{i});
    pV = setfield(pV,feilds{i},spm_diag(spm_vec(pF)));
end

% impose correlations between orientations (L)
%==========================================================================
n         = size(pE.Lpos,2);
Rpos      = pE.Lpos;
Rpos(1,:) = -Rpos(1,:);
D         = 128*ones(n);

% find symmetrical sources in each hemisphere
%--------------------------------------------------------------------------
for i = 1:n
    for j = 1:n
        if sign(pE.Lpos(1,i)) == sign(Rpos(1,j))
            D(i,j) = sqrt(sum(pE.Lpos(:,i) - Rpos(:,j)).^2);
        end
    end
end
D     = (D + D')/2;
DD    = zeros(n);
for i = 1:n
    [M, I] = min(D(i,:));
    if M < V
        DD(i,I) = 1;
    end
end

% reduce rank of prior covariance matrix of positions
%--------------------------------------------------------------------------
try
    pV.L = pV.L + kron(DD,diag(pV.L(1)*[-1 1 1]));
end

% and concatenate
%--------------------------------------------------------------------------
for  i = 1:length(feilds)
    pF      = getfield(pV,feilds{i});
    if ~isempty(pF)
        pC{i,i} = pF;
    end
end
pC    = spm_cat(pC);

author: spm


text:
function tapas_BLT_runStimAnswer(params)
% tapas_BLT_runStimAnswer is showing the stimulus question / answer screen

% Show stimulus question
Screen('DrawTexture', params.screen.window, params.text.imgstimanswer, [], params.screen.fit, 0);

% Flip and log
Screen('Flip', params.screen.window);
tapas_BLT_logEvent('stim_ans_on',params);

tic()
ellapsed = 0;
responded = 0;

while ellapsed < params.dur.stimAnsTimeout
    ellapsed = toc();
    keyCode = detectkey(params);

    if keyCode == params.keys.one
        responded = 1;
        rt      = ellapsed;

        % Code the answer
        if params.answertype == 1
            answer  = 1; % Yes
        elseif params.answertype == 2
            answer  = 0; % No
        end

        % Draw and flip the screen
        Screen('DrawTexture', params.screen.window, params.text.imgstimanswerSL, [], params.screen.fit, 0);
        Screen('Flip', params.screen.window);

    elseif keyCode == params.keys.two
        responded = 1;
        rt      = ellapsed;

        % Code the answer
        if params.answertype == 1
            answer  = 0; % No
        elseif params.answertype == 2
            answer  = 1; % Yes
        end

        % Draw and flip the screen
        Screen('DrawTexture', params.screen.window, params.text.imgstimanswerSR, [], params.screen.fit, 0);
        Screen('Flip', params.screen.window);
    end
end

% Display warning if too slow
if responded == 0
    Screen('DrawTexture', params.screen.window, params.text.imgpredtimeout , [], params.screen.fit, 0);
    Screen('Flip', params.screen.window);
    wait2_escapeOption(params.dur.showpredTimeoutScreen, params);
    answer = NaN;
    rt = NaN;
end

% Log the data
tapas_BLT_logData('stim_answer_rt', params, rt);
tapas_BLT_logData('stim_answer', params, answer);
tapas_BLT_logEvent('stim_ans_off',params);

end

author: translationalneuromodeling

repo: tapas


text:
function [bool, score] = score_find_in_cache(cache,j,ps,scoring_fn)
% [bool, score] = score_find_in_cache(cache,j,ps,scoring_fn)
%
% francois.olivier.c.h@gmail.com

%tic
L=size(cache,1);
N=size(cache,2)-3;

if N<1
  bool=0;
  score=0;
  return
end

parents=zeros(1,N+1);
parents(ps)=1;parents(N+1)=j;

switch scoring_fn
  case 'bic',
    fn=1;
  case 'bayesian',
    fn=2;
  otherwise,
    fn=3;
    %error(['unrecognized scoring fn ' scoring_fn]);
end

%parent = str2num(num2str(parents,'%1d'));
%[tmp y]=find(cache(:,N+3)==fn);
%if ~isempty(tmp)
%  [tmp2 y]=find(str2num(num2str(cache(tmp,1:N+1),'%1d'))==parent);
%  candidats=tmp(tmp2);
%else
%  candidats=[];
%end

[tmp y]=find(cache(2:L,N+3)==fn);
tmp=tmp+1;
[tmp2 y]=find(cache(tmp,N+1)==j);
candidats=tmp(tmp2);
if ~isempty(candidats)
  for i=1:N      % N=size(cache,2)-3;
    if ~isempty(candidats)
      [tmp2 y]=find(cache(tmp,i)==parents(i));
      candidats=intersect(candidats,tmp(tmp2));
    end
  end
end

%Tpre=toc

if ~isempty(candidats)
  bool=1;
else
  bool=0;
end

if bool
  score=cache(candidats(1),N+2);
else
  score=0;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%%
% Test for W1 on graph using linprog.

addpath('../toolbox/');
rep = 'results/w1-graph/';
[~,~] = mkdir(rep);

% helpers
normalize = @(x)x/sum(x(:));

% generate a planar graph
N = 80;
X = randn(2,N);
X(2,:) = X(2,:)*.75;
T = delaunay(X(1,:),X(2,:));

I = [T(:,1);T(:,2);T(:,2);T(:,3);T(:,3);T(:,1)];
J = [T(:,2);T(:,1);T(:,3);T(:,2);T(:,1);T(:,3)];
[I,J] = deal(I(I<J), J(I<J));
P = length(I);

% gradient operator
grad = sparse( [1:P 1:P], [I;J]', [ones(P,1) -ones(P,1)], P,N);

clf;
triplot(T,X(1,:),X(2,:));

% edge length
w = sqrt( (X(1,I)-X(1,J)).^2 + (X(2,I)-X(2,J)).^2 );
w = w(:);

% source/sink
[~,s] = min(X(1,:));
d = sqrt( (X(1,s)-X(1,:)).^2 + (X(2,s)-X(2,:)).^2 );
[~,t] = max(d);

% target measures
test = 'diracs';
test = 'spread'; Q = 14; % Q-NN

switch test
    case 'diracs'
        mu = zeros(N,1); mu(s)=1;
        nu = zeros(N,1); nu(t)=1;
    case 'spread'
        v = [s t]; Mu = {};
        for i=1:2
            % distance to target
            x = X(:,v(i));
            d = sqrt( (x(1)-X(1,:)).^2 + (x(2)-X(2,:)).^2 );
            [d,Is] = sort(d, 'ascend'); d = d(1:Q); Is = Is(1:Q);
            Mu{i} = zeros(N,1);
            Mu{i}(Is) = normalize(sqrt( 1-d/d(end) ));
        end
        mu = Mu{1}; nu = Mu{2};
end

%%
% display the measures.

clf; hold on;
triplot(T,X(1,:),X(2,:), 'k');
% vertex according to mass
for i=1:N
    plot(X(1,i),X(2,i), '.', 'MarkerSize', .01 + 50*mu(i)/max(mu), 'color', 'r');
    plot(X(1,i),X(2,i), '.', 'MarkerSize', .01 + 50*nu(i)/max(mu), 'color', 'b');
end
axis equal; axis off;
saveas(gcf, [rep 'inputs.png'], 'png');

%%
% Solve Primal and Dual

% max <u,mu-nu> s.t.  |grad(u)| <= w
cvx_solver sdpt3 % SeDuMi %
cvx_begin quiet % sdp quiet
cvx_precision high;
variable u(N,1); %  real;
% norm( grad*u, Inf ) <= 1;
abs( grad*u ) <= w;
maximize( sum( u .* (mu-nu) ) );
cvx_end

% min |w.*f|_1 s.t.   div(f)=mu-nu
cvx_solver sdpt3 % SeDuMi %
cvx_begin quiet % sdp quiet
cvx_precision high;
variable f(P,1); %  real;
(grad'*f) - (mu-nu)  == 0;
minimize( sum( abs(f).*w ) );
% maximize( sum( abs(f(:)).*w(:) ) );
cvx_end

% plot flow
K = find(abs(f)>1e-3);
clf; hold on;
triplot(T,X(1,:),X(2,:), 'k-');
for k=K'
    plot(X(1,[I(k) J(k)]),X(2,[I(k) J(k)]), 'k', 'LineWidth', 4);
end
% vertex according to mass
for i=1:N
    plot(X(1,i),X(2,i), '.', 'MarkerSize', .01 + 50*mu(i)/max(mu), 'color', 'r');
    plot(X(1,i),X(2,i), '.', 'MarkerSize', .01 + 50*nu(i)/max(mu), 'color', 'b');
end
axis equal; axis off;
saveas(gcf, [rep 'flow.png'], 'png');
saveas(gcf, [rep 'flow.eps'], 'epsc');

% plot distance function
K = find(abs(f)>1e-3);
D = u-min(u(:)); D = D/max(D(:));
h = min(abs( (grad*u) ./ w ), 1);
clf; hold on;
% triplot(T,X(1,:),X(2,:), 'k');
% edge according to saturation index
for k=1:P
    col = h(k) * [0 1 0] + (1-h(k)) * [0 0 0];
    plot(X(1,[I(k) J(k)]),X(2,[I(k) J(k)]), 'color',  col,'LineWidth', 2);
end
% shortest path
%for k=K'
%    plot(X(1,[I(k) J(k)]),X(2,[I(k) J(k)]), 'k:', 'LineWidth', 3);
%end
% vertex according to distance u
for i=1:N
    plot(X(1,i),X(2,i), '.', 'MarkerSize', 30, 'color', [D(i) 0 1-D(i)]);
end
plot(X(1,s),X(2,s), 'r.', 'MarkerSize', 40);
plot(X(1,t),X(2,t), 'b.', 'MarkerSize', 40);
axis equal; axis off;
saveas(gcf, [rep 'potential.png'], 'png');
saveas(gcf, [rep 'potential.epsc'], 'epsc');

author: optimaltransport


text:
function test_bug921

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_selectdata ft_selectdata_old ft_selectdata_new

% See also bug 798 that was reported by Yoni, from which I am reusing the data to test this bug

datadir = (dccnpath('/home/common/matlab/fieldtrip/data/test/bug798'));
load(fullfile(datadir,'t2_subj1'));
load(fullfile(datadir,'t2_subj1_null'));
load(fullfile(datadir,'t2_subj2'));
load(fullfile(datadir,'t2_subj2_null'));

argin{1} = ft_checkdata(t2_subj1, 'datatype', 'freq');
argin{2} = ft_checkdata(t2_subj1_null, 'datatype', 'freq');
argin{3} = ft_checkdata(t2_subj2, 'datatype', 'freq');
argin{4} = ft_checkdata(t2_subj2_null, 'datatype', 'freq');

% the following should select and concatenate only the powspctrm field
% but not the prob and stat fields
%data = ft_selectdata(argin{:}, 'param', 'powspctrm');

%assert(~isfield(data, 'stat'));
%assert(~isfield(data, 'prob'));
%assert(strcmp(data.dimord, 'rpt_chan_freq'));

data = cell(1,numel(argin));
[data{:}] = ft_selectdata([],argin{:});

cfg = [];
cfg.parameter = 'powspctrm';
cfg.appenddim = 'rpt';
data = ft_appendfreq(cfg, data{:});

assert(size(data.powspctrm,1)==4);
assert(size(data.powspctrm,2)==274);
assert(size(data.powspctrm,3)==1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function cx = getgain(X)

cx = X.gain;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function FreeDataCache(obj,ndataadd)

while true,

  if obj.ndatacached + ndataadd <= obj.maxdatacached || obj.ndatacached == 0,
    break;
  end
  if isempty(obj.perframehistory),
    break;
  end
  [~,j] = min(cell2mat(obj.perframehistory(:,2)));
  fn = obj.perframehistory{j,1};
  ClearDataCache(obj,fn);
%   for i = 1:obj.nexpdirs,
%     if isfield(obj.datacached{i},fn),
%       ndatacurr = 0;
%       for k = 1:numel(obj.datacached{i}),
%         ndatacurr = ndatacurr + numel(obj.datacached{i}(k).(fn));
%       end
%       obj.datacached{i} = rmfield(obj.datacached{i},fn);
%       obj.ndatacachedperexp(i) = obj.ndatacachedperexp(i) - ndatacurr;
%       obj.ndatacached = obj.ndatacached - ndatacurr;
%     end
%   end

  obj.perframehistory(j,:) = [];

end

author: kristinbranson

repo: JAABA


text:
function test_resampledesign

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY resampledesign

% since the function to test is in a private directory, we explicitely have to cd into that directory
[ftver, ftpath] = ft_version;
cd(fullfile(ftpath, 'private'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check that there are fac(N) permutations

for ntrial=1:10
  clear design
  design(1,:) = ones(1,ntrial);
  design(2,:) = 1:ntrial;

  cfg = [];
  cfg.ivar = 1;
  cfg.resampling = 'permutation';
  cfg.numrandomization = 'all';
  res = resampledesign(cfg, design);

  if size(res,1)~=max(cumprod(1:size(design,2)))
    error('incorrect number of permutations')
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the following is for Irina, who has a blocked fMRI experiment and wants
% to randomize the two conditions over the blocks, keeping the trials
% within a block together

clear design
design = [
  1 1 2 2 1 1 2 2 % condition number
  1 1 2 2 3 3 4 4 % block number
  ];

cfg = [];
cfg.ivar = 1; % condition number
cfg.wvar = 2;
cfg.resampling = 'permutation';
cfg.numrandomization = 'all';
res = resampledesign(cfg, design);

if any((res(:,2:2:end) - res(:,1:2:end)) ~= 1)
  error('trials within a block were not kept together');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear design
design = [
  1 1 2 2
  1 2 3 4 % block number
  ];

cfg = [];
cfg.ivar = 1; % condition number
cfg.wvar = 2;
cfg.resampling = 'permutation';
cfg.numrandomization = 'all';
res = resampledesign(cfg, design);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

design = [
  1 2 1 2 1 2 1 2
  1 2 3 4 5 6 7 8
  ];
cfg = [];
cfg.ivar = 1;
cfg.wvar = 2;
cfg.numrandomization = 'all';
cfg.resampling = 'permutation';
res = resampledesign(cfg, design);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

design = [
  1 2 3 4 5 6 7 8
  ];

cfg = [];
cfg.resampling = 'bootstrap';
cfg.numrandomization = 100;
res = resampledesign(cfg, design);

if ~all(size(res)==[100 8])
  error('incorrect bootstrap resampling');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

design = [
  1 2 2   % condition number
  1 2 3   % trial number
  ];
cfg = [];
cfg.ivar = 1;
cfg.numrandomization = 'all';
cfg.resampling = 'permutation';
cfg.efficient = 'no';  % NOTE THIS ONE
res = resampledesign(cfg, design);
if size(res,1)~=6
  error('incorrect number of permutations');
end

cfg.efficient = 'yes';  % NOTE THIS ONE
res = resampledesign(cfg, design);
if size(res,1)~=3
  error('incorrect number of permutations');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% here it is invalid to swap the condition for blocks, because the actual
% replications in the data would get scattered over blocks.

design = [
  1 1 1 2 2 1 1 1 2 2   % condition/ivar
  1 1 1 2 2 3 3 3 4 4   % block/wvar
  ];
cfg = [];
cfg.ivar = 1;
cfg.wvar = 2;
cfg.numrandomization = 'all';
cfg.resampling = 'permutation';
try
  % this should result in an error
  res = resampledesign(cfg, design);
  error_detected = false;
catch
  % this is ok
  error_detected = true;
end
if ~error_detected
  error('unequal block length with wvar should be caught as error');
end


text:
function [W, info] = MTFLC_ADMM(X, y, lambda1, lambda2, opts)
%
% Multi-Task Feature Learning with Calibration - ADMM
%
% OBJECTIVE
%    min_W { sum_i^m ||Xi wi - yi|| + lambda1 ||W||_{1,2} + lambda2/2 ||W||_F^2 }
%
% We use ADMM to solve the problem.
%
% INPUT
%  X - cell array of {n_i by d matrices} by m
%  y - cell array of {n_i by 1 vectors}  by m
%  lambda1 - regularization parameter of the l2,1 norm penalty
%  lambda2 - regularization parameter of the Fro norm penalty
%
% OUTPUT
%  W - task weights: d by t.
%  funcVal - the funcion value.
%
% Author: Jiayu Zhou

%% Initialization
if(nargin<5), opts = []; end

opts = setOptsDefault( opts, 'verbose', 1);
opts = setOptsDefault( opts, 'maxIter', 10000);
opts = setOptsDefault( opts, 'tol',     1e-4);
verbose = opts.verbose;

innerOpts = [];
innerOpts = setOptsDefault( innerOpts, 'maxIter', 5000);
innerOpts = setOptsDefault( innerOpts, 'tol',     1e-8);
innerOpts = setOptsDefault( innerOpts, 'tFlag',   1);

if verbose > 0
    fprintf('ADMM Config: [MaxIter %u][Tol %.4g]\n', opts.maxIter, opts.tol);
    fprintf('ADMM Subsolver Config: [MaxIter %u][Tol %.4g]\n', innerOpts.maxIter, innerOpts.tol);
end

m = length(X); % task number
d = size(X{1}, 2); % dimension.

funcVal  = zeros(opts.maxIter, 1);

% init for variables.
Th = cell(m, 1); % each element is n_i by 1
for tt = 1:m
    Th{tt} = L2proj(randn(size(X{tt}, 1), 1));
end
Z = Th; % this initialization indicates that w_init = 0.
W = zeros(d, m);

rho = 1; rhoInc = 1;

%% Computation
if verbose == 1; fprintf('Iteration:     '); end
for iter = 1: opts.maxIter

    % Update W
    innerOpts.init = W;
    [W, innerInfo] = MTFLC_ADMM_WSolver(X, y, Th, Z, rho, lambda1, lambda2, innerOpts); %#ok

    % Update Z (the same size as Tht)
    for t = 1: m
        vt = Th{t} / rho + y{t} - X{t} * W(:, t);
        vtNrm = sqrt(sum(vt.^2));
        Z{t} = max(0, 1- 1/(rho * vtNrm)) * vt;
    end

    % Update theta
    for t = 1: m
        Th{t} = Th{t} + rho * (y{t} - Z{t} - X{t} * W(:, t));
    end

    %if verbose >=2, fprintf('Th: %.4f\n', augLagObjective(W, Th, Z)); end

    % Update rho
    rho = rho * rhoInc;

    funcVal(iter) = augLagObjective(W, Th, Z);
    if verbose == 1; fprintf('\b\b\b\b\b%5i',iter); end

    % check stop criteria
    if (iter >1)
        diffW = sum(sum((W - W_old).^2));
        diffZ = 0; diffTh = 0;
        for t = 1: m
            diffZ  = diffZ  + sum((Z{t}  - Z_old{t}) .^2);
            diffTh = diffTh + sum((Th{t} - Th_old{t}).^2);
        end

        if verbose>1
            fprintf('Iter: Fv: %.4f Fv+L: %.4f dW %.4g, dZ %.4g, dTh %.4g\n', ...
                primalObjective(W), funcVal(iter), diffW, diffZ, diffTh);
        end

        if( max(diffW, max(diffZ, diffTh)) < opts.tol)
            break;
        end
    end

    W_old  = W; Z_old  = Z; Th_old = Th;
end
if verbose == 1; fprintf('\n'); end

%% Output
info.funcVal = funcVal(1: iter);
% use the last to show the 'real' objective without aug Lagrange.
info.funcVal(end + 1) = primalObjective(W);
% NOTE: the primal should be the same as augmented when converged.

%% Nested Functions
    function fvP = primalObjective(W)
        % primal objective (without augmented terms)
        %  P(W)  sum_i^m ||Xi wi - yi|| + lambda1 ||W||_{1,2} + lambda2/2 ||W||_F^2

        fvP = lambda1 * L21norm(W) + lambda2 /2 * sum(sum(W.^2));
        for i = 1: m
            fvP = fvP + sqrt(sum((X{i} * W(:, i) - y{i}).^2));
        end
    end

    function fvP = augLagObjective(W, Th, Z)
        % primal objective with augmented Lagrange terms.
        %  P(W)  sum_i^m ||Xi wi - yi|| + lambda1 ||W||_{1,2} + lambda2/2 ||W||_F^2
        %         + sum_i^m {Th_i' (y_i - z_i - X_i w_i) + 2 ||y_i - z_i - X_i w_i||^2/rho }

        fvP = lambda1 * L21norm(W) + lambda2 /2 * sum(sum(W.^2));
        for i = 1: m
            fvP = fvP + sqrt(sum((X{i} * W(:, i) - y{i}).^2)); % loss
            ti = y{i} - Z{i} - X{i} * W(:, i);
            fvP = fvP + Th{i}' * ti;             % dual
            fvP = fvP + 2 * sum(ti.^2) / rho;    % augment term.
        end
    end

end

function [Xnrm] = L21norm(X)
% ||X||_{1,2} = sum_i ||X^i||_2
Xnrm = sum(sqrt(sum(X.^2, 2)));
end

function [x] = L2proj(x)
nrm = sqrt(sum(x.^2));
x = x ./ max(1, nrm) ;
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function [favar]=loadsignres_favar(pref,favar,endo)

%%%%% Routine for information variables X
% preliminary tasks
% identify the number of endo
favar.numendo=size(endo,1);

% initiate the cells signrestable, signresperiods and signreslabels %%%% I think this step is not necessary
% favar.signrestable=cell(favar.nsignres,favar.numendo); %number of information variables to be restricted, one for the moment
% favar.signresperiods=cell(favar.nsignres,favar.numendo); %number of information variables to be restricted, one for the moment
% favar.signreslabels=favar.informationvariablestrings{1,1}; %number of information variables to be restricted, one for the moment
% load the data from Excel
% sign restrictions values
[~,~,strngs1]=xlsread(pref.excelFile,'favar.sign res values');
[~,~,strngs2]=xlsread(pref.excelFile,'favar.sign res periods');
% replace NaN entries by blanks
strngs1(cellfun(@(x) any(isnan(x)),strngs1))={[]};
strngs2(cellfun(@(x) any(isnan(x)),strngs2))={[]};
% convert all numeric entries into strings
strngs1(cellfun(@isnumeric,strngs1))=cellfun(@num2str,strngs1(cellfun(@isnumeric,strngs1)),'UniformOutput',0);
strngs2(cellfun(@isnumeric,strngs2))=cellfun(@num2str,strngs2(cellfun(@isnumeric,strngs2)),'UniformOutput',0);
% identify the non-empty entries (pairs of rows and columns)
[nerows1,neclmns1]=find(~cellfun('isempty',strngs1));
[nerows2,neclmns2]=find(~cellfun('isempty',strngs2));
% count the number of such entries
neentries1=size(nerows1,1);
neentries2=size(nerows2,1);
% all these entries contrain strings: fix them to correct potential user formatting errors
% loop over entries (value table)
for ii=1:neentries1
strngs1{nerows1(ii,1),neclmns1(ii,1)}=bear.utils.fixstring(strngs1{nerows1(ii,1),neclmns1(ii,1)});
end
% loop over entries (period table)
for ii=1:neentries2
strngs2{nerows2(ii,1),neclmns2(ii,1)}=bear.utils.fixstring(strngs2{nerows2(ii,1),neclmns2(ii,1)});
end

% information variables in X that are restricted
neclmns1index=neclmns1==1;
nerows1index=nerows1(neclmns1index,1);

favar.signresX=strngs1(nerows1index,1);
%number of restricted information variables
favar.nsignresX=size(favar.signresX,1);

%create indices for restricted information variables (advantage here: ordering in the sign res table is irrelevant)
    for jj=1:favar.nsignresX
        for ii=1:favar.nfactorvar
        favar.signresX_indexlogical{jj,1}(ii,1)=strcmp(favar.informationvariablestrings(1,ii),favar.signresX{jj,1}); % do we need this output in favar.?
        end
    end
    for jj=1:favar.nsignresX
        favar.signresX_index(jj,1)=find(favar.signresX_indexlogical{jj,1}==1);
    end

% sign restriction values

% loop over endogenous variables
for ii=1:favar.numendo
    [~,c]=find(strcmp(strngs1,endo{ii,1}));
for ll=1:favar.nsignresX
[r,~]=find(strcmp(favar.signresX,favar.informationvariablestrings{ll,1}));% first entry number of information variables to be restricted, one for the moment
   % if it is not possible to find two entries, return an error
%    if size(r,1)<2
%    message=['favar Sign restriction error: information variable ' favar.informationvariablestrings{ll,1} ' cannot be found in both rows and columns of the table. Please verify that the ''sign res values'' sheet of the Excel data file is properly filled.'];
%    msgbox(message);
%    error('programme termination: sign restriction error');
%    end
% otherwise, the greatest number in r corresponds to the row of the column labels: record it
rows(favar.nsignresX,1)=max(r);
% the greatest number in c corresponds to the column of the row labels: record it
clmns(ii,1)=max(c);
end
end
% now recover the values for the cell signrestable
% loop over endogenous (rows)
for ii=1:favar.nsignresX % loop over restricted information variables
   for jj=1:favar.numendo % loop over endogenous (columns)
   favar.signrestable{ii,jj}=strngs1{rows(ii,1),clmns(jj,1)};
   end
end
%favar sign res table can be integrated to the normal one. take the
%signrestable definition (endpoints) from the normal one, and just start
%the favar table after that

% check whether the restriction table is valid: return an error if the the column requirement for zero restrictions is not satisfied
% loop over columns
for ii=1:favar.numendo
% count the number of zero restrictions in this column
numzerores=sum(strcmp(favar.signrestable(:,ii),'0'));
   % if there are too many zero restrictions for the column, return an error
   if numzerores>favar.numendo-ii
   temp=['Zero restriction issue: you have requested ' num2str(numzerores) ' zero restrictions in column ' num2str(ii) ' of the restriction matrix, but at most ' num2str(favar.numendo-ii) ' such restrictions can be implemented.'];
   msgbox(temp);
   error('Zero restriction error');
   end
end

% sign restriction periods

% recover the rows and columns of each endogenous variable
% loop over endogenous variables
for ii=1:favar.numendo
% for each variable, there should be two entries in the table corresponding to its name
% one is the column lable, the other is the row label
[~,c]=find(strcmp(strngs1,endo{ii,1}));
[r,~]=find(strcmp(strngs1,favar.informationvariablestrings{ll,1}));% first entry number of information variables to be restricted, one for the moment
   % if it is not possible to find two entries, return an error
%    if size(r,1)<2
%    message=['Sign restriction error: endogenous variable ' favar.informationvariablestrings{ii,1} ' cannot be found in both rows and columns of the table. Please verify that the ''sign res periods'' sheet of the Excel data file is properly filled.'];
%    msgbox(message);
%    error('programme termination: sign restriction error');
%    end
% the greatest number in r corresponds to the row of the column labels: record it
rows(favar.nsignresX,1)=max(r);
% the greatest number in c corresponds to the column of the row labels: record it
clmns(ii,1)=max(c);
end
% now recover the values for the cell signresperiods
% loop over endogenous (rows)
for ii=1:1 %number of information variables to be restricted, one for the moment
   % loop over endogenous (columns)
   for jj=1:favar.numendo
   % record the value
   favar.signresperiods{ii,jj}=str2num(strngs2{rows(ii,1),clmns(jj,1)});
   end
end

% recover the labels, if any
% initiate
favar.signreslabels=cell(favar.numendo,1);
% loop over endogenous (columns)
for ii=1:favar.numendo
% the label for shock ii is found in strngs2, row 3 and column 'clmns(ii,1)'
temp=strngs2{1,clmns(ii,1)};
   % if empty, give the generic name 'shock ii'
   if isempty(temp)
   favar.signreslabels{ii,1}=['shock ' num2str(ii)];
   % else, record the name
   else
   favar.signreslabels{ii,1}=temp;
   end
end

% finally, record on Excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),strngs1,'favar.sign res values','B2');
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),strngs2,'favar.sign res periods','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
function out=atanh(x)

out=(log(1+x)-log(1-x))/2;

author: opencobra

repo: cobratoolbox


text:
function xoverKids  = crossoverNcK(parents,options,GenomeLength,FitnessFcn,unused,thisPopulation)
% Oren Rosen
% The MathWorks
% 8/29/2007
%
% This custom crossover function is written to work on a population of
% vectors of zeros and ones with the same amount of ones in each vector.
% The children that are produced from 2 parents will have the same genes
% for every element they agree on, and random choices of zerps and ones for
% the elements they don't agree on. The number of each is set so that all
% children have the same number of ones as their parents. All children
% automatically satisfy this constraint so there is no need to impose these
% constraints.

% How many children to produce?
nKids = length(parents)/2;

% Allocate space for the kids
xoverKids = zeros(nKids,GenomeLength);

% To move through the parents twice as fast as thekids are
% being produced, a separate index for the parents is needed
index = 1;

% *** Initialize ***
% Assumes all members of thisPopulation have the same number of ones.
num1s = sum(thisPopulation(1,:));
indexVec = 1:GenomeLength;

% for each kid...
for i=1:nKids

    % *** Get parents ***
    r1 = parents(index);
    index = index + 1;
    r2 = parents(index);
    index = index + 1;

    p1 = thisPopulation(r1,:);
    p2 = thisPopulation(r2,:);

    % *** Find Matching 1's and 0's ***
    % Ex: If           p1 == [ 1 0 1 0 0 1 1 0 0 0 ]
    %                  p2 == [ 1 0 0 1 0 0 1 0 1 0 ]
    %     Then matching1s == [ 1 0 0 0 0 0 1 0 0 0 ]
    %     Then matching0s == [ 0 1 0 0 1 0 0 1 0 1 ]
    matching1s = ~xor(p1,p2) & (p1 == 1);
    matching0s = ~xor(p1,p2) & (p1 == 0);

    % *** Find Matching Indices ***
    %     If       matching1s == [ 1 0 0 0 0 0 1 0 0 0 ]
    %              matching0s == [ 0 1 0 0 1 0 0 1 0 1 ]
    %     Then matching1sIndx == [ 1 7 ]
    %          matching0sIndx == [ 2 5 8 10 ]
    %         nonmatchingIndx == [ 3 4 6 9 ]
    matching1sIndx = indexVec(matching1s);
    matching0sIndx = indexVec(matching0s);
    nonmatchingIndx = setdiff(indexVec,[matching0sIndx,matching1sIndx]);

    % *** Create Child ***
    % Ex: If   num1s == 4
    %          matching1sIndx == [ 1 7 ]
    %          nonmatchingIndx == [ 3 4 6 9 ]
    %     Then numMatching1s == 2
    %     num1sToFill == 2
    %     Indx1sToFill == 2 random choices from [ 3 4 6 9 ]
    numMatching1s = numel(matching1sIndx);
    num1sToFill = num1s - numMatching1s;
    Indx1sToFill = randsample(nonmatchingIndx,num1sToFill);

    % *** Fill in 1s ***
    % Ex: If               p1 == [ 1 0 1 0 0 1 1 0 0 0 ]
    %                      p2 == [ 1 0 0 1 0 0 1 0 1 0 ]
    %     Then xoverKids(i,:) == [ 1 0 ? ? 0 ? 1 0 ? 0 ]
    %     With exactly 2 of the '?' equal to 1, the rest 0.
    xoverKids(i,matching1sIndx) = 1;
    xoverKids(i,Indx1sToFill) = 1;
end

author: Sable

repo: mcbench-benchmarks


text:
function batch = StartBatch(mfile,bfile,varargin)

%StartBatch - Start a new batch job.
%
% Batch processing is useful if you need to repeatedly perform a given analysis
% on different data sets. Using batch processing is easier and more reliable
% than writing a simple 'for' loop:
%
%  * You must only supply the code to run the analysis on one data set,
%    and a text file listing the data sets, including specific parameters
%
%  * Running the batch on a subset of the data is as easy as commenting
%    out the appropriate lines in this file
%
%  * Batch processing will not be interrupted when the analysis of one
%    particular data set generates an error
%
%  * Errors can be logged to disk
%
%  * Processing can be delayed to a later time (e.g. during the night)
%
% Results can either be kept in memory using the batch output variable, or
% stored in a central database. The latter is the recommended solution, as it
% is far more powerful and reliable. See EXAMPLES below.
%
% During batch processing, figures are typically hidden so they do not keep
% popping up as the user is performing other tasks as the batch is running in
% the background (although this is configurable). In order to keep figures
% hidden, batch functions should not use figure(h) or axes(a) which make the
% figure visible. Two helper functions, scf(h) and sca(h), can be used instead
% to set the current figure or axes without making the figure visible.
%
%  USAGE
%
%    batch = StartBatch(mfile,bfile,<options>)
%
%    mfile          batch function (M-file name or function handle)
%    bfile          batch file listing the parameters for each iteration
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'delay'       delay (in minutes) before execution starts
%     'hide'        by default, new figures are hidden to speed up batch
%                   processing (use 'off' to cancel this behavior)
%     'log'         error log file (default = none) (see NOTE below)
%    =========================================================================
%
%  NOTE
%
%    Log file names can include wildcards to indicate current date and time:
%
%      %y    year
%      %m    month
%      %d    day
%      %t    time
%
%  EXAMPLE 1 - RESULTS KEPT IN BATCH OUTPUT VARIABLE
%
%    This is the simpler solution, but it is not as reliable (e.g. the results
%    can be lost if Matlab crashes) and powerful (e.g. the results cannot be
%    shared between users) as using a database. See EXAMPLE 2 below.
%
%    Suppose you wish to count the number of 'ripples' in a set of sleep
%    sessions recorded from different rats. You could create the file
%    'batch.txt', where you would list the session files to process, as
%    well as the LFP channel to use in each case to detect ripples:
%
%        Rat-01-20100915-01-sleep     5
%        Rat-01-20100920-01-sleep     5
%        Rat-02-20101001-01-sleep     3
%        Rat-01-20101005-04-rest      1
%        ...
%
%    The function used to process the data would look like:
%
%        function n = CountRipples(session,channel)
%
%        SetCurrentSession(session);
%        lfp = GetLFP(channel);
%        fil = FilterLFP(lfp,'bandpass','ripples');
%        ripples = FindRipples(fil);
%        n = size(ripples,1);
%
%    To start processing the data in one hour, you would then call:
%
%        b = StartBatch(@CountRipples,'batch.txt','delay',60);
%
%    To get the results:
%
%        n = GetBatch(b);
%
%    To log errors:
%
%        b = StartBatch(@CountRipples,'batch.txt','log','Errors-%y%m%d.log');
%
%  NOTE
%
%    For technical reasons, the batch function must have a fixed number of output
%    parameters, e.g. it cannot have 'varargout' as one of its output parameters.
%    It is possible to circumvent this limitation if you know in advance the number
%    of output parameters that the batch function would actually return in your
%    particular case. Assuming this function is called 'VariableProcess', takes two
%    inputs and will always return three outputs for your data, you could define:
%
%        function [x,y,z] = FixedProcess(u,v)
%
%        [x,y,z] = VariableProcess(u,v);
%
%    and then use this as the batch function.
%
%  EXAMPLE 2 - RESULTS STORED IN CENTRAL DATABASE
%
%    Let us rewrite the previous example using a central database. First, we need
%    to create the database:
%
%        % This assumes your login/password are stored
%        % in a configuration file in your home directory
%        DBConnect;
%        DBCreate('Ripples_%y%m%d');
%
%    The function used to process the data would now look like:
%
%        function CountRipples(session,channel)
%
%        SetCurrentSession(session);
%        lfp = GetLFP(channel);
%        fil = FilterLFP(lfp,'bandpass','ripples');
%        ripples = FindRipples(fil);
%        n = size(ripples,1);
%        eid = [session '-' channel];
%        DBAddVariable(n,eid,'Ripple Count','','',{'CountRipples'});
%
%    This would store the number of ripples for each session, as well as useful
%    information such as who stored these data, when, using which function (the
%    actual code, not just the function name, is stored in the database).
%
%    To get the results:
%
%        n = DBGetValues;
%
%  SEE
%
%    See also Database, GetBatch, BatchInfo, CancelBatch, CleanBatches, Store,
%    Recall, Hide, sca, scf.
%

% Copyright (C) 2007-2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Default values
delay = 0;
hide = 'off';
log = '';

% Check number of parameters
if nargin < 2,
	error('Incorrect number of parameters (type ''help StartBatch'' for details).');
end

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

% Compatibility with previous versions: reverse parameter order if necessary
if isa(bfile,'function_handle') || (isstring_FMAT(mfile) && isempty(which(mfile))),
	tmp = bfile;
	bfile = mfile;
	mfile = tmp;
end

% Batch function name
if isa(mfile,'function_handle'),
	mfileName = func2str(mfile);
else
	mfileName = mfile;
end

% Check batch file and function are valid
if ~isstring_FMAT(bfile) || ~exist(bfile,'file'),
	error('Batch file not found (type ''help StartBatch'' for details).');
end
if isempty(which(mfileName)),
	error('Batch function not found (type ''help StartBatch'' for details).');
end

% Make sure M-file has a defined number of output parameters
if nargout(mfile) == -1,
	error(['''' mfileName ''' has a variable number of output arguments (type ''help StartBatch'' for details). ']);
end

% Open batch file
f = fopen(bfile,'r');
if f == -1, error(['Could not open file ''' bfile '''.']); end

% Parse batch file
b = ParseBatch(bfile);

% Open log file
b.log = -1;
if ~isempty(log),
	% Insert date
	log = InsertDate(log);
	b.log = fopen(log,'w');
	if b.log == -1, error(['Could not open log file ''' log '''.']); end
end

% Set batch function
[~,mfile] = fileparts(mfileName);
b.mfile = mfile;

b.hide = strcmp(hide,'on');
if b.hide,
	warning('Figures will be hidden (type ''help Hide'' for details).');
end

% Start timer
batch = timer('TimerFcn',{@RunBatch,b},'StartDelay',delay*60,'Tag','BatchJob');
start(batch);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function err = getHdiverrorBDM1(node,elem,divSigma,sigmah,markedElem)
%% GETHDIVERRORBDM1 Hdiv norm of the approximation error for BDM1.
%
% The input divSigma is a function bundle, and sigma_h is a vector array
% whose component including two parts, the first part (i.e., sigmah(1:NE))
% is the line integral of flux in norm direction, i.e.
% sigmah_i = \int_{e_i} sigma\cdot n.
% the second part (i.e., sigmah(NE+1:2*NE)) is the dual basis written as
% sigmah_i = 3\int_{e_i} (\lambda1-\lambda2) sigma\cdot n  (ei = < 1 , 2 >)
%
% err = getHdiverrorBDM1(node,elem,divSigma,sigmah)
%
% Example
%
%     maxIt = 5;
%     node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0]; % nodes
%     elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];    % elements
%     bdEdge = setboundary(node,elem,'Dirichlet');
%     pde = mixBCdata;
%     err = zeros(maxIt,1); N = zeros(maxIt,1);
%     for i =1:maxIt
%         [node,elem,bdEdge] = uniformrefine(node,elem,bdEdge);
%         [u,sigma,NULL] = PoissonBDM1(node,elem,pde,bdEdge);
%         err(i) = getHdiverrorBDM1(node,elem,pde.f,-sigma,[]);
%         N(i) = size(u,1);
%     end
%     r1 = showrate(N,err,2);
%     legend('||\sigma - \sigma_h||_{H(div)}',['N^{' num2str(r1) '}'],...
%            'LOCATION','Best');
%
% See also getHdiverror3RT0, getL2errorRT0.
%
% Created by Ming Wang at Jan 17, 2011, M-lint modified at May 15, 2011.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Construct Data Structure
[elem2dof,~,elem2edgeSign] = dofedge(elem);
NT = size(elem,1);
% compute div phi
[Dlambda,area] = gradbasis(node,elem);
% div phi = 2*(Dlambda_i,Rot_j);
rotMat = [0 -1; 1 0]; % rotation matrix for computing rotLambda.
divPhi(:,3) = 2*dot(Dlambda(:,:,1),Dlambda(:,:,2)*rotMat,2);
divPhi(:,1) = 2*dot(Dlambda(:,:,2),Dlambda(:,:,3)*rotMat,2);
divPhi(:,2) = 2*dot(Dlambda(:,:,3),Dlambda(:,:,1)*rotMat,2);
divSigmahp = zeros(NT,1);
for k = 1:3
    divSigmahp = divSigmahp + elem2edgeSign(:,k).*sigmah(elem2dof(:,k)).*divPhi(:,k);
end

%% compute Hdiv error element-wise
[lambda,w] = quadpts(2);
nQuad = size(lambda,1);
err = zeros(NT,1);
for p = 1:nQuad
    % quadrature points in the x-y-z coordinate
    pxy = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:);
    divSigmap = divSigma(pxy);
    % compute divSigmahp at quadrature points
    err = err + w(p)*sum((divSigmap - divSigmahp).^2,2);
end
err = err.*area;
% modify the error
err(isnan(err)) = 0; % remove the singular part
if (nargin == 5) && ~isempty(markedElem)
    err = err(markedElem); % error on some marked region
end
err = sqrt(sum(err));

repo: ifem


text:
%--- help for ts/chebyshev_box ---
%
%  Constructs chebyshev boxes for multivariate-multiperiods densities
%
%  ::
%
%    mvcb=chebyshev_box(this,gam)
%
%  Args:
%
%     this (ts | rts) : time series with many pages (number of simulations)
%       and potentially many columns (number of variables)
%
%     gam (scalar | vector) : percentile(s)
%
%  Returns:
%     :
%
%     - **mvcb** [struct] : structure with time series in its fields. Each
%       field represents a particular variable
%
%     - **gam_** [scalar|vector] : sorted percentile(s)
%
%     - **my** [ts] : mean across simulations
%
%

author: jmaih


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

% Doing some broadcasting (https://www.gnu.org/software/octave/doc/interpreter/Broadcasting.html)
%warning ("off", "Octave:broadcast");

% calculate all Euclidian distances betweek all X and all centroids (a^2 + b^2 = c^2)
%distances = sqrt(((X(:,1) - centroids(:,1)').^2) + ((X(:,2) - centroids(:,2)').^2));

% get the min value per row and corresponding index (which should be the centroid id)
%[minValues,idx]=min(distances,[],2);

% Restore broadcasting
%warning ("on", "Octave:broadcast");

for i = 1:length(idx)
    distance = zeros(K, 1);
    for j = 1:K
        % TODO(SaveTheRbtz@): Can be vectorized as diff * diff'
        distance(j) = sum(sum((X(i, :) - centroids(j, :)) .^ 2 ));
    endfor
    [value, idx(i)] = min(distance);
endfor

% =============================================================

end

author: vugsus

repo: coursera-machine-learning


text:
function fx1 = p10_fx1 ( x )

%*****************************************************************************80
%
%% P10_FX1 evaluates the derivative of the function for problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX1, the first derivative of the function at X.
%
  fx1 = ( 1.0 - 10.0 * x ) * ( 1.0 + 10.0 * x ) / ( 100.0 * x * x + 1.0 ).^2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
load('data');
[Rtilt,R] = rectify(cat(3,X,Y,Z));

figure,
XYZnew = Rtilt*[X(:),Y(:),Z(:)]';
scatter3(XYZnew(1,1:10:end),XYZnew(2,1:10:end),XYZnew(3,1:10:end),1,rgb(1:10:end,:));
axis equal;

figure,
XYZnew = R*[X(:),Y(:),Z(:)]';
scatter3(XYZnew(1,1:10:end),XYZnew(2,1:10:end),XYZnew(3,1:10:end),1,rgb(1:10:end,:));
axis equal;

author: jianxiongxiao

repo: ProfXkit


text:
function [F,failure,cause,ALREADY_MODELLED] = expandmodel(F,h,options,w)

% FIX : Current code experimental, complex, conservative, has issues with
% nonlinearities and is slow...
%
% If it wasn't for polynomials and sigmonials, it would be trivial, but the
% code is extremely messy since we want to have this functionality too

global LUbounds
global ALREADY_MODELLED
global MARKER_VARIABLES
global DUDE_ITS_A_GP
global ALREADY_MODELLED_INDEX
global REMOVE_THESE_IN_THE_END
global OPERATOR_IN_POLYNOM
global CONSTRAINTCUTSTATE

% All extended variables in the problem. It is expensive to extract this
% one so we will keep it and pass it along in the recursion
extendedvariables = yalmip('extvariables');

integers = yalmip('intvariables');
binaries = yalmip('binvariables');
if ~isempty(F)
    declareI = find(is(F,'integer'));
    declareB = find(is(F,'binary'));
    if ~isempty(declareI)
        integers = [integers depends(F(declareI))];
    end
    if ~isempty(declareB)
        binaries = [binaries depends(F(declareB))];
    end
end
yalmip('settempintvariables',integers);
yalmip('settempbinvariables',binaries);

% We keep track of all auxilliary variables introduced by YALMIP
nInitial = yalmip('nvars');

boundsAlreadySet = 0;
% Meta constraints are expanded first (iff, implies, alldifferent etc)
if ~isempty(F)
    meta = find(is(F,'meta'));
    if ~isempty(meta)
        LUbounds=setupBounds(F+extractglobalboundsmeta(F),options,extendedvariables);
        boundsAlreadySet = 1;
        F = expandmeta(F);
    end
end

if isa(F,'constraint')
    F = lmi(F);
end

if nargin < 3
    options = sdpsettings;
end

if isempty(options)
    options = sdpsettings;
end

if isfield(options,'reusemodel')

else

    ALREADY_MODELLED = [];

    % Temporary hack to deal with a bug in CPLEX. For the implies operator (and
    % some more) YALMIP creates a dummy variable x with (x==1). Cplex fails
    % to solve problem with these stupid variables kept, hence we need to
    % remove these variables and constraints...
    MARKER_VARIABLES = [];

    % Temporary hack to deal with geometric programs. GPs are messy here,
    % becasue we can by mistake claim nonconvexity, since we may have no
    % sigmonial terms but indefinite quadratic term, but the whole problem is
    % meant to be solved using a GP solver. YES, globals suck, but this is
    % only temporary...hrm.
    DUDE_ITS_A_GP = 0;

    % Keep track of expressions that already have been modelled. Note that if a
    % graph-model already has been constructed but we now require a milp, for
    % numerical reasons, we should remove the old graph descriptions (important
    % for MPT models in particular)
    % FIX: Pre-parse the whole problem etc (solves the issues with GP also)
    ALREADY_MODELLED = {};
    ALREADY_MODELLED_INDEX = [];
    REMOVE_THESE_IN_THE_END = [];

    % Nonlinear operator variables are not allowed to be used in polynomial
    % expressions, except if they are exactly modelled, i.e. modelled using
    % MILP models. We will expand the model and collect variables that are in
    % polynomials, and check in the end if they are exaclty modelled
    OPERATOR_IN_POLYNOM = [];
end
% Assume success
failure = 0;
cause = '';

% Early bail out
if isempty(extendedvariables)
    return
end

if nargin < 4
    w = [];
end

if isempty(F) & isempty(h)
    return
end

% Check if it already has ben expanded in robustify or similar
F_alreadyexpanded = [];
if ~isempty(F)
    F_alreadyexpanded = [];
    already_expanded = expanded(F);
    if all(already_expanded)
        if isempty(setdiff(getvariables(h),expanded(h)))
            return
        end
    elseif any(already_expanded)
        F_alreadyexpanded = F(find(already_expanded));
        F = F(find(~already_expanded));
    end
end

if ~isempty(F)
    % Extract all simple bounds from the model, and update the internal bounds
    % in YALMIP. This is done in order to get tighter big-M models

    if boundsAlreadySet == 0;
        LUbounds = setupBounds(F,options,extendedvariables);
    end

    % Expand equalities first, since these might generate nonconvex models,
    % thus making it unnecessaryu to generate epigraphs etc
    equalities = is(F,'equality');
    if any(equalities)
        F = [F(find(equalities));F(find(~equalities))];
    end
end

% All variable indicies used in the problem
v1 = getvariables(F);
v2 = depends(F);
v3 = getvariables(h);
v4 = depends(h);

% HACK: Performance for LARGE-scale dualizations
if isequal(v3,v4) & isequal(v1,v2)
    variables = uniquestripped([v1 v3]);
else
    variables = uniquestripped([v1 v2 v3 v4]);
end

% Index to variables modeling operators
extended = find(ismembcYALMIP(variables,extendedvariables));

if nargin < 3
    options = sdpsettings;
end

% This is a tweak to allow epxansion of bilinear terms in robust problems,
% is expression such as abs(x*w) < 1 for all -1 < w < 1
% This field is set to 1 in robustify and tells YALMIP to skip checking for
% polynomial nonconvexity in the propagation
if ~isfield(options,'expandbilinear')
    options.expandbilinear = 0;
end

% Monomial information. Expensive to retrieve, so we pass this along
[monomtable,variabletype] = yalmip('monomtable');

% Is this trivially a GP, or meant to be at least?
if strcmpi(options.solver,'gpposy') | strcmpi(options.solver,'fmincon-geometric') | strcmpi(options.solver,'mosek-geometric')
    DUDE_ITS_A_GP = 1;
else
    if ~isequal(options.solver,'fmincon') & ~isequal(options.solver,'snopt') & ~isequal(options.solver,'ipopt') & ~isequal(options.solver,'') &  ~isequal(options.solver,'mosek')
        % User has specified some other solver, which does not
        % support GPs, hence it cannot be intended to be a GP
        DUDE_ITS_A_GP = 0;
    else
        % Check to see if there are any sigmonial terms on top-level
        DUDE_ITS_A_GP = ~isempty(find(variabletype(variables) == 4));
    end
end

% Constraints generated during recursive process to model operators
F_expand = ([]);

if isempty(F)
    F = ([]);
end

% First, check the objective
variables = uniquestripped([depends(h) getvariables(h)]);
monomtable = monomtable(:,extendedvariables);

% However, some of the variables are already expanded (expand can be called
% sequentially from solvemp and solverobust)
variables = setdiff1D(variables,expanded(h));

% Determine if we should aim for MILP/EXACT model directly
if options.allowmilp == 2
    method = 'exact';
else
    method = 'graph';
end

if DUDE_ITS_A_GP == 1
    options.allowmilp = 0;
    method = 'graph';
end

% Test for very common special case with only norm expression
ExtendedMap = yalmip('extendedmap');
fail = 0;
if  0%length(ExtendedMap) > 0 &&  all(strcmp('norm',{ExtendedMap.fcn}))
    for i = 1:length(ExtendedMap)
        if ~isequal(ExtendedMap(i).arg{2},2)
            fail = 1;
            break;
        end
        if ~isreal(ExtendedMap(i).arg{1})
            fail = 1;
            break;
        end
        if any(ismembcYALMIP(getvariables(ExtendedMap(i).arg{1}),extendedvariables))
             fail = 1;
            break;
        end
    end
    for i = 1:length(ExtendedMap)
        F_expand = [F_expand, cone(ExtendedMap(i).arg{1},ExtendedMap(i).var)];
    end
    F = F + lifted(F_expand,1);
    return
end
% *************************************************************************
% OK, looks good. Apply recursive expansion on the objective
% *************************************************************************
index_in_extended = find(ismembcYALMIP(variables,extendedvariables));
allExtStructs = yalmip('extstruct');
if ~isempty(index_in_extended)
    [F_expand,failure,cause] = expand(index_in_extended,variables,h,F_expand,extendedvariables,monomtable,variabletype,'objective',0,options,method,[],allExtStructs,w);
end

% *************************************************************************
% Continue with constraints
% *************************************************************************
constraint = 1;
all_extstruct = yalmip('extstruct');
while constraint <=length(F) & ~failure

    if ~already_expanded(constraint)
        Fconstraint = F(constraint);
        variables = uniquestripped([depends(Fconstraint) getvariables(Fconstraint)]);

        % If constraint is a cut, all generated constraints must be marked
        % as cuts too
        CONSTRAINTCUTSTATE =  getcutflag(Fconstraint);
        [ix,jx,kx] = find(monomtable(variables,:));
        if ~isempty(jx) % Bug in 6.1
            if any(kx>1)
                OPERATOR_IN_POLYNOM = [OPERATOR_IN_POLYNOM extendedvariables(jx(find(kx>1)))];
            end
        end

        index_in_extended = find(ismembcYALMIP(variables,extendedvariables));
        if ~isempty(index_in_extended)
            if is(Fconstraint,'equality')
                if options.allowmilp | options.allownonconvex
                    [F_expand,failure,cause] = expand(index_in_extended,variables,-sdpvar(Fconstraint),F_expand,extendedvariables,monomtable,variabletype,['constraint #' num2str(constraint)],0,options,'exact',[],allExtStructs,w);
                else
                    failure = 1;
                    cause = ['integer model required for equality in constraint #' num2str(constraint)];
                end
            else
                [F_expand,failure,cause] = expand(index_in_extended,variables,-sdpvar(Fconstraint),F_expand,extendedvariables,monomtable,variabletype,['constraint #' num2str(constraint)],0,options,method,[],allExtStructs,w);
            end
        end
    end
    constraint = constraint+1;
end

CONSTRAINTCUTSTATE = 0;

% *************************************************************************
% Temporary hack to fix the implies operator (cplex has some problem on
% these trivial models where a variable only is used in x==1
% FIX: Automatically support this type of nonlinear operators
% *************************************************************************
if ~isempty(MARKER_VARIABLES)
    MARKER_VARIABLES = sort(MARKER_VARIABLES);
    equalities = find(is(F,'equality'));
    equalities = equalities(:)';
    remove = [];
    for j = equalities
        v = getvariables(F(j));
        if length(v)==1
            if ismembcYALMIP(v,MARKER_VARIABLES)
                remove = [remove j];
            end
        end
    end
    if ~isempty(remove)
        F(remove) = [];
    end
end

nNow = yalmip('nvars');
if nNow > nInitial
    % YALMIP has introduced auxilliary variables
    % We mark these as auxilliary
    yalmip('addauxvariables',nInitial+1:nNow);
end

F_expand = lifted(F_expand,1);
% *************************************************************************
% We are done. We might have generated some stuff more than once, but
% luckily we keep track of these mistakes and remove them in the end (this
% happens if we have constraints like (max(x)<1) + (max(x)>0) where
% the first constraint would genrate a graph-model but the second set
% creates a integer model.
% *************************************************************************
if ~failure
    F = F + F_expand;
    if length(REMOVE_THESE_IN_THE_END) > 0
        F = F(find(~ismember(getlmiid(F),REMOVE_THESE_IN_THE_END)));
    end
end

% *************************************************************************
% Normally, operators are not allowed in polynomial expressions. We do
% however allow this if the variable has been modelled with an exact MILP
% model.
% *************************************************************************
if ~failure
    dummy = unique(OPERATOR_IN_POLYNOM);
    if ~isempty(dummy)
        for i = 1:length(dummy)
            aux(i,1) = find(ALREADY_MODELLED_INDEX == dummy(i));
        end
        %    Final_model = {ALREADY_MODELLED{unique(OPERATOR_IN_POLYNOM)}};
        Final_model = {ALREADY_MODELLED{aux}};
        for i = 1:length(Final_model)
            if ~(strcmp(Final_model{i}.method,'integer') | strcmp(Final_model{i}.method,'exact') | options.allownonconvex)
                failure = 1;
                cause = 'Nonlinear operator in polynomial expression.';
                return
            end
        end
    end
end

% Append the previously appended
F = F + F_alreadyexpanded;
% Declare this model as expanded
F = expanded(F,1);

function [F_expand,failure,cause] = expand(index_in_extended,variables,expression,F_expand,extendedvariables,monomtable,variabletype,where,level,options,method,extstruct,allExtStruct,w)
global DUDE_ITS_A_GP ALREADY_MODELLED ALREADY_MODELLED_INDEX REMOVE_THESE_IN_THE_END OPERATOR_IN_POLYNOM

% *************************************************************************
% Go through all parts of expression to check for convexity/concavity
% First, a small gateway function before calling the recursive stuff
% *************************************************************************
if ~DUDE_ITS_A_GP
    [ix,jx,kx] = find(monomtable(variables,:));
    if ~isempty(jx) % Bug in 6.1
        if any(kx>1)
            OPERATOR_IN_POLYNOM = [OPERATOR_IN_POLYNOM extendedvariables(jx(find(kx>1)))];
        end
    end
end

failure = 0;
j = 1;

try
    % Optimized for 6.5 and higher (location in ismember). If user has
    % older version, it will crash and go to code below
    expression_basis = getbase(expression);
    expression_vars  = getvariables(expression);
    [yesno,location] = ismember(variables(index_in_extended),expression_vars);
    ztemp = recover(variables(index_in_extended));
    while j<=length(index_in_extended) & ~failure
        % i = index_in_extended(j);
        % zi = recover(variables(i));
        zi = ztemp(j);%recover(variables(i));
        basis = expression_basis(:,1 + location(j));
        if all(basis == 0) % The nonlinear term is inside a monomial
            if options.allownonconvex
                [F_expand,failure,cause] = expandrecursive(zi,F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,'exact',[],'exact',allExtStruct,w);
            else
                failure = 1;
                cause = 'Possible nonconvexity due to operator in monomial';
            end
            %[F_expand,failure,cause] = expandrecursive(zi,F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],'exact',allExtStruct,w);
        elseif all(basis >= 0)
            [F_expand,failure,cause] = expandrecursive(zi,F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],'convex',allExtStruct,w);
        else
            [F_expand,failure,cause] = expandrecursive(zi,F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],'concave',allExtStruct,w);
        end
        j=j+1;
    end
catch
    while j<=length(index_in_extended) & ~failure
        i = index_in_extended(j);
        basis = getbasematrix(expression,variables(i));
        if all(basis >= 0)
            [F_expand,failure,cause] = expandrecursive(recover(variables(i)),F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],'convex',allExtStruct,w);
        else
            [F_expand,failure,cause] = expandrecursive(recover(variables(i)),F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],'concave',allExtStruct,w);
        end
        j=j+1;
    end
end


text:
function value = r8vec_diff_norm_l1 ( n, a, b )

%*****************************************************************************80
%
%% R8VEC_DIFF_NORM_L1 returns the L1 norm of the difference of R8VEC's.
%
%  Discussion:
%
%    The vector L1 norm is defined as:
%
%      value = sum ( 1 <= I <= N ) abs ( A(I) ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in A.
%
%    Input, real A(N), B(N), the vectors.
%
%    Output, real VALUE, the L1 norm of A - B.
%
  value = sum ( abs ( a(1:n) - b(1:n) ) );

  return
end


text:
function value = q_measure ( n, z, triangle_order, triangle_num, triangle_node )

%*****************************************************************************80
%
%% Q_MEASURE determines the triangulated pointset quality measure Q.
%
%  Discussion:
%
%    The Q measure evaluates the uniformity of the shapes of the triangles
%    defined by a triangulated pointset.
%
%    For a single triangle T, the value of Q(T) is defined as follows:
%
%      TAU_IN = radius of the inscribed circle,
%      TAU_OUT = radius of the circumscribed circle,
%
%      Q(T) = 2 * TAU_IN / TAU_OUT
%        = ( B + C - A ) * ( C + A - B ) * ( A + B - C ) / ( A * B * C )
%
%    where A, B and C are the lengths of the sides of the triangle T.
%
%    The Q measure computes the value of Q(T) for every triangle T in the
%    triangulation, and then computes the standard deviation of this
%    set of values:
%
%      Q_MEASURE = min ( all T in triangulation ) Q(T)
%
%    In an ideally regular mesh, all triangles would have the same
%    equilateral shape, for which Q = 1.  In a good mesh, 0.5 < Q.
%
%    Given the 2D coordinates of a set of N nodes, stored as Z(1:2,1:N),
%    a triangulation is a list of NT triples of node indices that form
%    triangles.  Generally, a maximal triangulation is expected, namely,
%    a triangulation whose image is a planar graph, but for which the
%    addition of any new triangle would mean the graph was no longer planar.
%    A Delaunay triangulation is a maximal triangulation which maximizes
%    the minimum angle that occurs in any triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Max Gunzburger and John Burkardt,
%    Uniformity Measures for Point Samples in Hypercubes.
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, pages 329-345, June 2004.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real Z(DIM_NUM,N), the points.
%
%    Input, integer TRIANGLE_ORDER, the order of the triangles.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(TRIANGLE_ORDER,TRIANGLE_NUM), the triangulation.
%
%    Output, real VALUE, the Q quality measure.
%
  if ( triangle_num < 1 )
    value = -1.0;
    return
  end
%
%  Compute the mean value of Q.
%
  q_min = r8_huge ( );

  for triangle = 1 : triangle_num

    a_index = triangle_node(1,triangle);
    b_index = triangle_node(2,triangle);
    c_index = triangle_node(3,triangle);

    ab_length = sqrt ( ...
        ( z(1,a_index) - z(1,b_index) )^2 ...
      + ( z(2,a_index) - z(2,b_index) )^2 );

    bc_length = sqrt ( ...
        ( z(1,b_index) - z(1,c_index) )^2 ...
      + ( z(2,b_index) - z(2,c_index) )^2 );

    ca_length = sqrt ( ...
        ( z(1,c_index) - z(1,a_index) )^2 ...
      + ( z(2,c_index) - z(2,a_index) )^2 );

    q = ( bc_length + ca_length - ab_length ) ...
      * ( ca_length + ab_length - bc_length ) ...
      * ( ab_length + bc_length - ca_length ) ...
      / ( ab_length * bc_length * ca_length );

    q_min = min ( q_min, q );

  end

  value = q_min;

  return
end


text:
function f=burger(u,der)
% if second arg. absent u^2/2, else u.
if nargin<2
	f=0.5*u.^2;
else
	f=u;
end;

repo: Aero-matlab


text:
function planC = CERRImportDICOM(study_directory, planC_save_name, defaultSelection)
%"CERRImportDICOM"
%   Calls DICOM routines to import a DICOM directory to CERR file.
%
%   By default use no parameters, which passes the user to a GUI.  To use
%   in batch mode, pass a directory path and save file name.  The directory
%   path must contain ONLY ONE DICOM study, and will have only the first
%   scan,dose,etc extracted.
%
%Last Modified:  28 Aug 03, by ES.
%                added control for proper exit when user cancelled DICOM import
%
%                7 Feb 05, JRA Added batch mode control.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Start diary and timer for import log.
startTime = now;
tmpFileName = tempname;
diary(tmpFileName);
Ver = version;
tic;

if exist('study_directory')
    [ap,ap_xmesh,ap_ymesh,ap_zmesh,ap_ct,ap_ct_xmesh,ap_ct_ymesh,ap_ct_zmesh,ap_voi,MRCell,PETCell,SPECTCell]=dicomrt_DICOMimport(study_directory, 'r', defaultSelection);
else
    [ap,ap_xmesh,ap_ymesh,ap_zmesh,ap_ct,ap_ct_xmesh,ap_ct_ymesh,ap_ct_zmesh,ap_voi,MRCell,PETCell,SPECTCell]=dicomrt_DICOMimport;
end

if isempty(ap)==1 & ...
        isempty(ap_xmesh)==1 & ...
        isempty(ap_ymesh)==1 & ...
        isempty(ap_zmesh)==1 & ...
        isempty(ap_ct)==1 & ...
        isempty(ap_ct_xmesh)==1 & ...
        isempty(ap_ct_ymesh)==1 & ...
        isempty(ap_ct_zmesh)==1 & ...
        isempty(ap_voi)==1 & ...
        isempty(MRCell) & ...
        isempty(PETCell) & ...
        isempty(SPECTCell) , ...
        return;
else
    [planC] = dicomrt_dicomrt2cerr(ap,ap_xmesh,ap_ymesh,ap_zmesh,ap_ct,...
        ap_ct_xmesh,ap_ct_ymesh,ap_ct_zmesh,ap_voi,'CERROptions.m',MRCell,PETCell,SPECTCell);
end

%Stop diary and write it to planC.
diary off;
endTime = now;
logC = file2cell(tmpFileName);
delete(tmpFileName);
indexS = planC{end};
planC{indexS.importLog}(1).importLog = logC;
planC{indexS.importLog}(1).startTime = datestr(startTime);
planC{indexS.importLog}(1).endTime = datestr(endTime);

toc;
% Save plan  -- this line moved out from dicomrt_dicomrt2cerr.  If a passed
% save name was used, save there, else prompt.
if exist('planC_save_name')
    save_planC(planC,planC{indexS.CERROptions}, 'passed', planC_save_name);
else
    save_planC(planC,planC{indexS.CERROptions});
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function bti_test02 ( )

%*****************************************************************************80
%
%% BTI_TEST02 tests BURGERS_TIME_INVISCID.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BTI_TEST02\n' );
  fprintf ( 1, '  Test BURGERS_TIME_INVISCID\n' );

  method = 2;
  nx = 81;
  nt = 200;
  t_max = 2.0;
  bc = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Method: 2, Upwind conservative.\n' );
  fprintf ( 1, '  Initial condition: expansion\n' );
  fprintf ( 1, '  Number of space nodes = %d\n', nx );
  fprintf ( 1, '  Number of time steps = %d\n', nt );
  fprintf ( 1, '  Final time T_MAX = %g\n', t_max );
  fprintf ( 1, '  Boundary condition = %d\n', bc );

  U = burgers_time_inviscid ( method, @ic_expansion, nx, nt, t_max, bc );

  x = linspace ( -1.0, +1.0, nx );

  figure ( 2 )

  plot ( x, U(1:50:(nt+1),:), 'Linewidth', 3 )
  grid on
  xlabel ( '<-- X -->' )
  ylabel ( '<-- U(X,T) -->' )
  title ( 'Burgers, inviscid, initial expansion, upwind conservative' )

  filename = 'bti_test02.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot as "%s"\n', filename );

  return
end


text:
function [ precision, recall, PR_curve ] = evaluateDetectionCNN( objects, nElems )
%EVALUATEODCNN Evaluates the result obtained by the ObjectDetector-CNN
% w.r.t. a given intersection over union (IoU) value.

    min_conf_per = 0:0.02:1;

    %% Get labels and scores for each TP sample
    [ all_TPs, all_FPs, TP_, FP_ ] = getTPs(objects, nElems);

    %% Only continue if we have any detection
    nDetections = TP_+FP_;
    if(nDetections > 0)
        %% Sort detection confidence scores
        sorted_scores = [];
        if(TP_ > 0)
            sorted_scores = [sorted_scores; all_TPs(:,4)];
        end
        if(FP_ > 0)
            sorted_scores = [sorted_scores; all_FPs(:,4)];
        end
        sorted_scores = sort(sorted_scores, 'descend');

        %% Build precision-recall curve
        rec = zeros(1,length(min_conf_per));
        prec = zeros(1,length(min_conf_per));
        count = 1;
        for conf_per = min_conf_per

            % Only evaluate then score is >= score on current percentage position
            pos_min_score = max(round(nDetections*conf_per), 1);
            min_score = sorted_scores(pos_min_score);
            [TP, FP, FN] = measuresByConfidence(objects, min_score, all_TPs, all_FPs, TP_, FP_);

            if(conf_per == 1)
                % precision = 1 - frac_NoObj;
                % recall = detection_rate;

                precision = TP / (TP + FP); % how many are correct from all that it found?
                recall = TP / (TP + FN); % how many did it find from the total?
            end

            prec(count) = TP / (TP + FP);
            rec(count) = TP / (TP + FN);
            count = count+1;
        end
    else
        precision = 0;
        recall = 0;
        rec = zeros(1,length(min_conf_per));
        prec = zeros(1,length(min_conf_per));
    end

    PR_curve.precision = prec;
    PR_curve.recall = rec;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [w,u,eqn,info] = biharmonicP2(node,elem,bdFlag,pde,option)
%%
%   [w,u] = biharmonicP1(node,elem,bdFlag,pde) produces the mixed cubic finite element
%   approximation of the biharmonic equation, where w = laplace u
%   See also biharmonicP1, biharmonicP2, biharmonicP3.
%
%   Created by Jie Zhou. Clean up is needed
%
%   Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if nargin<5, option = []; end
tic;
%% Construct Data Structure
[elem2dof,edge,bdDof] = dofP2(elem);
N = size(node,1);  NT = size(elem,1);  Ndof = N+size(edge,1);

%% Compute geometric quantities and gradient of local basis
[Dlambda,area] = gradbasis(node,elem);

%% Assemble stiffness matrix
% Since Dphi_i*Dphi_j is quadratic, numerical quadrature rule is used here
if ~isfield(option,'quadorder')
    option.quadorder = 4;   % default order
end
[lambda, weight] = quadpts(option.quadorder);
nQuad = size(lambda,1);
ii = zeros(21*NT,1); jj = zeros(21*NT,1); sA = zeros(21*NT,1);sB = zeros(21*NT,1);
index = 0;
for i = 1:6
    for j = i:6
        Bij = 0;
        Aij = 0;
        for p = 1:nQuad
                 Bij = Bij + weight(p)*dot(Dphi(p,i),Dphi(p,j),2);
                 Aij = Aij + weight(p)*dot(phi(p,i),phi(p,j),2);
        end
        Bij = Bij.*area;
        Aij = Aij.*area;

        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        sB(index+1:index+NT) = Bij;
        sA(index+1:index+NT) = Aij;
        index = index + NT;
    end
end
clear Aij Bij
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
B = sparse(ii(diagIdx),jj(diagIdx),sB(diagIdx),Ndof,Ndof);
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Ndof,Ndof);
% A = spdiags(accumarray(ii(diagIdx),sA(diagIdx),[Ndof 1]),0,Ndof,Ndof);
BU = sparse(ii(upperIdx),jj(upperIdx),sB(upperIdx),Ndof,Ndof);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Ndof,Ndof);
B = B + BU + BU';
A = A + AU + AU';

%% boundary condition
%%
    fixedDof = [];
    isFixedDof = false(Ndof,1);
    if ~isempty(bdFlag)
        elem2edge = elem2dof(:,4:6)-N;
        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        isFixedDof(edge(isDirichlet,:)) = true;
        isFixedDof(N + find(isDirichlet')) = true;
        fixedDof = find(isFixedDof);
        freeDof = find(~isFixedDof);
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction Dphi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function s = Dphi(p,i) % gradient of basis phi
    switch i
        case 1
            s = (4*lambda(p,1)-1).*Dlambda(:,:,1);
        case 2
            s = (4*lambda(p,2)-1).*Dlambda(:,:,2);
        case 3
            s = (4*lambda(p,3)-1).*Dlambda(:,:,3);
        case 4
            s = 4*(lambda(p,2)*Dlambda(:,:,3)+lambda(p,3)*Dlambda(:,:,2));
        case 5
            s = 4*(lambda(p,3)*Dlambda(:,:,1)+lambda(p,1)*Dlambda(:,:,3));
        case 6
            s = 4*(lambda(p,1)*Dlambda(:,:,2)+lambda(p,2)*Dlambda(:,:,1));
    end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction Dphi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function s = phi(p,i) % gradient of basis phi
    switch i
        case 1
            s = (2*lambda(p,1)-1).*lambda(p,1);
        case 2
            s = (2*lambda(p,2)-1).*lambda(p,2);
        case 3
            s = (2*lambda(p,3)-1).*lambda(p,3);
        case 4
            s = 4*lambda(p,3).*lambda(p,2);
        case 5
            s = 4*lambda(p,1).*lambda(p,3);
        case 6
            s = 4*lambda(p,1).*lambda(p,2);
    end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Assemble right hand side by high order quadrature rule
% To reduce the effect of the error introduced by the numerical quadrature,
% the load term is computed using the 3rd order qudrature rule.
b = zeros(Ndof,1);
% u=zeros(Ndof,1);
w = zeros(Ndof,1);

if ~isfield(option,'fquadorder')
    option.fquadorder = 3;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    % quadrature points in the barycentric coordinate
    [lambda,w] = quadpts(option.fquadorder);
    nQuad = size(lambda,1);
%     phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
%     phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
%     phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
%     phi(:,4) = 4*lambda(:,2).*lambda(:,3);
%     phi(:,5) = 4*lambda(:,3).*lambda(:,1);
%     phi(:,6) = 4*lambda(:,1).*lambda(:,2);
    bt = zeros(NT,6);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        if isfield(pde,'f') && isnumeric(pde.f)
            fp = pde.f;        % piecewise constant
        else
            fp = pde.f(pxy);   % function handle
        end
        for j = 1:6
            bt(:,j) = bt(:,j) + w(p)*phi(p,j)*fp;
        end
    end
    bt = bt.*repmat(area,1,6);
    b = accumarray(elem2dof(:),bt(:),[Ndof 1]);
end

[b1,u] = getbdP2(b);
B(:,fixedDof)=[];
Nu=size(freeDof,1);

b(fixedDof)=[];

switch solver
    case 'direct'
        t = cputime;
        bigA = [A, B; ...
                B', sparse(Nu,Nu)];
        bigF = [b1; -b];
        bigu = bigA\bigF;
        w = bigu(1:Ndof);
        u(freeDof)=bigu(Ndof+1:Ndof+Nu);
        residual = norm(bigF-bigA*bigu);
        info = struct('solverTime',cputime - t,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
end % end of four order.

%     soln = struct('u',u,'Du',Du);
    eqn = struct('A',A,'B',B','f',b1,'g',-b,'freeNode',freeNode,'Lap',A);
    info.assembleTime = assembleTime;

 function [b1,u] = getbdP2
    %% Boundary conditions for Poisson equation: P2 quadratic FEM.
    %
    % The set up of boundary condition consists of two parts:
    %

    %
    %  Modify the right hand side b. The Neumann boundary integral is added
    %  to b. For Dirichlet boundary ndoes, b(fixedDof) is the evaluation of
    % pde.g_D.
    %
    u = zeros(Ndof,1);

    %% Part 1: Find boundary edges and modify the load b
    % Neumann boundary condition
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 6;
        end
                idxN = (bdFlag(:) == 1);      % all Neumann edges in bdFlag
        Neumannidx = elem2edge(idxN ); % index of Neumann and Robin edges
        % since boundary integral is also needed for Robin edges
        Neumann   = edge(Neumannidx,:);
        b1 = zeros(Ndof,1);
        [lambdagN,weightgN] = quadpts1(option.gNquadorder);
        nQuadgN = size(lambdagN,1);
        % quadratic bases (1---3---2)
        bdphi = zeros(nQuadgN,3);
        bdphi(:,1) = (2*lambdagN(:,1)-1).*lambdagN(:,1);
        bdphi(:,2) = (2*lambdagN(:,2)-1).*lambdagN(:,2);
        bdphi(:,3) = 4*lambdagN(:,1).*lambdagN(:,2);
        % length of edge
        el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
        ge = zeros(size(Neumann,1),3);
        for pp = 1:nQuadgN
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:);
            gNu = pde.g_N(ppxy);
            ge(:,1) = ge(:,1) + weightgN(pp)*gNu*bdphi(pp,1);
            ge(:,2) = ge(:,2) + weightgN(pp)*gNu*bdphi(pp,2);
            ge(:,3) = ge(:,3) + weightgN(pp)*gNu*bdphi(pp,3); % interior bubble
        end
        % update RHS
        ge = ge.*repmat(el,1,3);
        b1(1:N) = accumarray(Neumann(:), [ge(:,1); ge(:,2)],[N,1]);
        b1(N+Neumannidx) = b1(N+Neumannidx) + ge(:,3);

    %% Part 2: Find Dirichlet boundary edges and compute the boundary value
    % Dirichlet boundary conditions

        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        % interpolation
                idx = (fixedDof > N);                         % index of edge nodes
        u(fixedDof(~idx)) = pde.g_D(node(fixedDof(~idx),:)); % bd value at vertex dofs
        bdEdgeIdx = fixedDof(idx) - N;
        bdEdgeMid = (node(edge(bdEdgeIdx,1),:) + node(edge(bdEdgeIdx,2),:))/2;
        u(fixedDof(idx)) = pde.g_D(bdEdgeMid);
        b1 = b1 - B*u;

    end % end of getbdP2

end                 % end of function biharmonicP2