text:
function M = optimizeGA(GA, varargin)
% :Usage:
% ::
%
%     M = optimizeGA(GA, [epochdur, dur_in_sec])
%
% outputs a pseudo-random list of condition codes that optimizes
% multiple fitness measures for fMRI task designs
%
% more help can be found in ga_example_script.m and other examples
% and in Genetic_Algorithm_readme.rtf
%
%
% :Inputs:
%
%   **GA:**
%        structure with input specifications; see ga_example_script.m
%
% :Optional Inputs:
%   **epochdur:**
%        Followed by epoch duration for events in sec
%        Note: you can also enter GA.epochdur with dur in sec
%
% :Outputs:
%
%   **M:**
%        Structure with model outputs
%
%         M.ga = GA inputs from script  (all input variables, plus metrics across generations)
%
%         M.ga.HRF = HRF used (default = event, but altered if you enter the 'epochdur' optional input
%         M.ga.bestCbal = counterbalancing metric for best overall design in each generation
%         M.ga.bestCVar = weighted contrast efficiency metric for best overall design in each generation
%         M.ga.bestHrf = HRF estimation (FIR) metric for best overall design in each generation
%
%         M.ga.bestFreq = Frequency deviation from targets for best overall design in each generation
%
%         M.ga.popVar = Population variability for each generation (high is good; low may reduce performance)
%         M.ga.totalGenTime = Time for each generation to run
%         M.ga.qualifyinglists = qualifyinglists;
%         M.ga.numStimEachCond = stimuli per condition; See ga_example_script.m
%         M.ga.numGenerations = number of generations;
%         M.ga.sizeGenerations = size of Generations;
%         M.ga.nonlinthreshold = Threshold for nonlinear predictor squashing
%         M.ga.listMatrix = Lists tested in latest generation
%         M.bestLists = Stimulus sequence for best list in each generation
%         M.date = date run
%
%        These fields describe the final model:
%        stimlist: your original list of conditions
%           delta: delta matrix sampled at .01 seconds
%           model: model sampled at .01 seconds
%       modelatTR: model sampled at TR
%        smoothed: temporally smoothed model at TR
%        conmodel: predictors for contrasts
%       consmooth: temporally smoothed predictors for contrasts
%           power: power (abs(fft)) matrix for smoothed model
%        conPower: power matrix for smoothed contrasts
%          energy: sum of squared deviations from mean for each predictor
%       conEnergy: sum of squared deviations from mean for each contrast
%             eff: efficiency of design.  1/trace(xtxi)
%  resampleeffinf: efficiency loss due to resampling compared to TR = 1 s.  eff(TR=1) / eff(TR)
%    smootheffinf: loss of efficiency due to smoothing; eff(nosmooth) / eff(smooth)
%         varbeta: variance of each individual beta
%          coneff: efficiency of contrast set
%      convarbeta: efficiency of each contrast beta
%             vif: variance inflation factors for predictors
%          conVif: variance inflation factors for contrasts
%            cond: condition number of model; larger is worse, > 30 is bad
%         conCond: condition number for contrast set
%           colin: correlation coefficients between individual predictors
%        conColin: correlation coefficients between contrasts
%       cbalfreqs: counterbalancing frequencies up to 3rd order
%            freq: frequencies with which each stimulus type occurs
%
% To plot output for final list separately, try this:
%
% M = modelDiagnostics2(M,'screen');

% ..
%    Tor Wager 2/10/01
%    Edited 11/17/01
%    Edited 11/8/09    - add check for designs with missing trial types;
%                          - output and error checking a bit more
%                          informative
%    Edited 8/15/10    - fixed minor bug in contrast length checking with
%                       trans2block option
%    Edited 12/3/2020  - add optional epoch length and checking for missing regressors, label typo
% ..

% -------------------------------------------------------------------------
% OPTIONAL INPUTS
% -------------------------------------------------------------------------

% This is a compact way to assign multiple variables. The input argument
% names and variable names must match, however:

epochdur = 0;
allowable_inputs = {'epochdur'};

% optional inputs with default values - each keyword entered will create a variable of the same name

if isfield(GA, 'epochdur'), epochdur = GA.epochdur; end

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case allowable_inputs

                eval([varargin{i} ' = varargin{i+1}; varargin{i+1} = [];']);

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

N = fieldnames(GA);
for i = 1:length(N)
    eval([N{i} ' = GA.' N{i} ';'])
end

% d-optimality flag - compute d rather than a-optimality
if ~exist('dflag')==1,dflag = 0; end

clockStart = cputime;			  	% keeps track of starting time

disp('OptimizeGA')
disp('===============================')
disp(['Saturation threshold = ' num2str(nonlinthreshold)])

% Set up variables and organisms
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %
disp('Setting up variables')

% ----------------------------------------------------------------
% * jitter setup
% ----------------------------------------------------------------

dojitter = 0;
if size(freqConditions,2) == size(conditions,2) + 1
    conditions = [conditions max(conditions)+1];   % insert 0 as last condition if jitter
    dojitter = 1;
elseif size(freqConditions,2) == size(conditions,2)
else error('Frequency and condition vectors must be same length, or freq must be 1 longer for jitter')
end

nconds = sum(conditions > 0);

fprintf('You entered %01d conditions  (event types) that will each be modeled with a regessor. \n ', nconds);
fprintf('Note that if you enter 0 as a condition number, those events will not be modeled\n (they are assumed to be rest/inter-stimulus interval periods)\n\n');
pause(.5)

% ----------------------------------------------------------------
% * saturation setup
% ----------------------------------------------------------------
if isempty(nonlinthreshold) || strcmp(nonlinthreshold,'none'),dosaturation = 0; nonlinthreshold = []; else dosaturation = 1; end

% ----------------------------------------------------------------
% * contrast setup
% ----------------------------------------------------------------

% work out whether to use contrasts.
if isempty(contrasts)
    disp('  ...no contrasts entered; Calculating efficiency of model betas (regression parameters)')
    docontrasts = 0;

else
    csize = size(contrasts);
    fprintf('Contrasts entered. You entered matrix of %01d contrasts, with weights across %01d conditions.\n ', csize);

    docontrasts = 1;

    disp('Contrasts:')
    disp(contrasts)

    badconlength = (~(trans2switch || trans2block)) && ~(size(contrasts,2) == nconds);
    badconlength = badconlength | (trans2block && ~(size(contrasts,2) == nconds*2));
    if badconlength
        error(sprintf('The number of columns in your contrast matrix must equal the number of conditions\n(excluding condition number entries of 0)'))
    end

    disp(' ...adding one column to contrast matrix with weights of 0, to model the intercept')

    contrasts(:, end+1) = 0;  % add 0s contrasts to account for intercept

    if ~isempty(contrastweights)

         fprintf('You entered a contrast weight vector across %01d contrasts\n ', length(contrastweights));
         if length(contrastweights) ~= size(contrasts, 1)
             error('The length of contrastweights should equal the number of contrasts (rows) in your contrast matrix.');
         end

    end
end

pause(2)

% ----------------------------------------------------------------
% * setup freqConditions and rests
% ----------------------------------------------------------------

% normalize freqConditions
%if ~sum(freqConditions) == 1,
%	freqConditions(1:end-1) = (1 - freqConditions(end)) / (size(freqConditions,2)-1);
%	disp(['	...freqConditions does not sum to 1: normalizing to ' num2str(freqConditions)])
%end
if ~sum(freqConditions) > 1, error('freqConditions must sum to 1 or less than 1.'),end

% flag for rest length.
	if ~isempty(restevery) &&  ~isempty(restlength)
		disp(['	...Using rests of length ' num2str(restlength) ' and resting every ' num2str(restevery)])
		dorests = 1;
	else disp('	...No rests specified.'),dorests = 0;
	end

if ~(mod(scanLength / ISI,1) == 0),disp('Warning: Scan length in s is not an even multiple of ISI!'),end

% ----------------------------------------------------------------
% * initial computation of list lengths and output
% ----------------------------------------------------------------

fitnessMatrix = zeros(4,sizeGenerations);
numStim = ceil(scanLength / (ISI));
if dorests
	numRestStim = (ceil(numStim/(mean(restevery)+restlength)) - 1) * restlength;
else
	numRestStim = 0;
end
numStimEachCond = ceil((numStim-numRestStim) * freqConditions);			            % row vector of stim in each cond
numBlanks = numStim - numRestStim - sum(numStimEachCond);
numsamps = ceil(numStim*ISI/TR);
%lowerLimit = round(scanLength / TR * lowerLimit);
disp('================== Parameters ======================')
disp(['     Total stimuli in scan      : ' num2str(numStim)])
disp(['     Fixed rest length in ITIs  : ' num2str(restlength)])
disp(['     Fixed Rest stimuli in scan : ' num2str(numRestStim)])
disp(['     Num blank intervals        : ' num2str(numBlanks)])
disp(['     Num of task conds          : ' num2str(sum(conditions > 0))])
disp(['     Num stimuli each cond      : ' num2str(numStimEachCond)])
disp('====================================================')

% ----------------------------------------------------------------
% * get smoothing matrix and autocorrelation matrix
% ----------------------------------------------------------------
if ~isfield(GA,'LPsmooth'),GA.LPsmooth = 1; end
[S, ~, svi] = getSmoothing(HPlength,GA.LPsmooth,TR,numsamps,xc);

% ----------------------------------------------------------------
% * HRF GAMMA FUNCTION - spm 99
% ----------------------------------------------------------------

   HRF = spm_hrf(.1);						% SPM HRF sampled at .1 s

   if epochdur
       % Epoch dur in sec - convolve

       hrfnew = conv(ones(epochdur * 10, 1), HRF);
       hrfnew = hrfnew/ max(hrfnew);
       HRF = hrfnew;

   end

   HRF = HRF/ max(HRF);

% ----------------------------------------------------------------
% * generate an unsorted list of conditions
% ----------------------------------------------------------------

disp('	...generating lists of organisms')
unsortedList = [];
for i = 1:size(conditions,2)
    unsortedList = [unsortedList; conditions(i) * ones(numStimEachCond(i),1)];
end
unsortedList = [unsortedList; zeros(numBlanks,1)];

% ----------------------------------------------------------------
% * transform to switches warnings (and generate unsorted list)
% ----------------------------------------------------------------

% This is part of the switching hack. uses 2 conditions, then figures out 4 columns from that.
% 5 is the rest condition.
if trans2switch
    if size(conditions,2) > 2
        error('Can''t use trans2switch with other than 2 conditions.')
    end
	disp(['	...USING TRANS2SWITCH > num stim each cond = ' num2str(numStimEachCond)])

end % if trans2switch

% ----------------------------------------------------------------
% * transform to block design - double number of conditions
% ----------------------------------------------------------------
if trans2block
	disp(['	...USING TRANS2BLOCK >'])
        freqConditions = [freqConditions./2 freqConditions./2];
        disp(['     freqConditions = ' num2str(freqConditions)])

        if dojitter
           maxcond = max(conditions)-1;
           conditions = [conditions(1:maxcond) maxcond+1:2*length(conditions)-1];
           disp(['     CONDITION ' num2str(conditions(end)) ' is blank jitter.'])
	else
           maxcond = max(conditions);
           conditions = [conditions maxcond+1:2*length(conditions)];
        end
        disp(['     conditions = ' num2str(conditions)])

	% ----------------------------------------------------------------
	% * setup special first trial
	% ----------------------------------------------------------------
	if isempty(dofirst)
		dofirst = 0;
    end
    if dofirst
		disp(['	...block design: adding special predictor for first trial in block.'])
		try
			% this is really kludgy, but should give a rough estimate of how many first trials to expect.
			numfirst = numStim ./ mean(restevery);
			proportionfirst = numfirst./numStim;
		catch
			error('restevery not found...must use rest periods with transform2block!')
		end
		freqConditions = freqConditions - (proportionfirst./size(freqConditions,2));
		freqConditions = [freqConditions proportionfirst];
		disp('	...first trial special condition - adjusted expected frequencies:')
		disp(['     freqConditions = ' num2str(freqConditions)])

		disp(['	...Adding extra column to contrasts to account for 1st trial predictor:'])
		contrasts = [contrasts zeros(size(contrasts,1),1)];
		contrasts
	end
end

% ----------------------------------------------------------------
% * signal to user if jitter is used.
% ----------------------------------------------------------------
if dojitter
   disp(['REST INTERVALS DETECTED: condition ' num2str(max(conditions))])
else
   disp(['REST INTERVALS: not used (to use, freqCond should be one longer than conditions'])
end

% ----------------------------------------------------------------
% * randomize initial set of organisms
% ----------------------------------------------------------------

disp('  ...randomizing organism start state')
for z = 1:sizeGenerations % generate random ordering of x conditions
  	stimList = getRandom(unsortedList);
    listMatrix(:,z) = stimList; % a row for each subject
end

% ----------------------------------------------------------------
% * transform restevery into restMatrix
% ----------------------------------------------------------------

%	restMatrix is a matrix with period length in ITIs X list (cols)
%		that specifies how long to go for each period between a rest/probe.
if dorests

	atleast = min(restevery);
	restPeriods = numStim / atleast + 5;	% make this many rows in restMatrix
	restList = [];

	for i = 1:size(restevery,2):restPeriods
		restList = [restList;restevery'];
    end

	for i = 1:sizeGenerations
		restMatrix(:,i) = getRandom(restList);
    end

	% disp('=============== Rest calculation ===================')
	disp(['	...restMatrix is ' num2str(size(restMatrix,1)) ' X ' num2str(size(restMatrix,2))])
	restMatrix(1:5,1:10)

    % ----------------------------------------------------------------
    % * insert rests, if needed, and trim list
    % ----------------------------------------------------------------
	%	There are two listMatrices: one to be tested for fitness, also called rna, and
	%		the other for crossover (called dna).
	%
	%	The first one has the rests in it, and the 2nd just has the stimuli.
	%	#1 is called TESTlistMatrix, and #2 is called listMatrix.
	%
	%	restMatrix is a matrix with period length in ITIs X list (cols)
	%		that specifies how long to go for each period between a rest/probe.
    %
    %   We need to generate TESTlistMatrix here to make sure it's the right size.
    %   later, we'll use it in the GA.
    %
	% could also use the function dna2rna, but we need only the rests at this point
	TESTlistMatrix = insert_rests(listMatrix,restMatrix,restlength,numStim);

else
	TESTlistMatrix = listMatrix;
	restMatrix = [];
end

% ----------------------------------------------------------------
% * allow user-entered listmatrix, to use non-random start.
% ----------------------------------------------------------------
if isfield(GA,'listMatrix')

	disp('	...Existing listMatrix detected.  Starting with matrix in GA.listMatrix.')
	try
		listMatrix(:,1:size(GA.listMatrix,2)) = GA.listMatrix;
	catch
		whos listMatrix
		error('Input listMatrix size does not match generated size.')
	end
end

if size(listMatrix,2) > sizeGenerations
	listMatrix = listMatrix(:,1:sizeGenerations);
end

% ----------------------------------------------------------------
% * check for other errors, and output warnings
% ----------------------------------------------------------------

% This should never happen if using rests, because insert_rests should adjust list length to numStim if necessary, but...
if size(TESTlistMatrix,1) < numStim

    disp(['Warning: # stim = ' num2str(numStim) '; list length = ' num2str(size(TESTlistMatrix,1)) ' : Padding with zeros.'])
    listMatrix(end+1:numStim,:) = 0;

elseif size(TESTlistMatrix,1) > numStim,disp(['Warning: # stim = ' num2str(numStim) '; list length = ' num2str(size(TESTlistMatrix,1)) ' : Truncating lists.'])
    listMatrix = listMatrix(1:numStim,:);
end

% initialize bestMat to size of final list to be saved (rna)
 bestMat = uint8(zeros(numStim,1));

disp('======== Num stimuli, averaged across lists (before instruction/rest stimuli) ========')
for i = 0:length(freqConditions)

	num = mean(sum(listMatrix == i));
	freq = num / numStim;
	disp(['		Condition ' num2str(i) ': ' num2str(num) ' stimuli, ' num2str(freq*100) '% of total stimuli'])

end
disp('======================================================================================')

clear TESTlistMatrix
listMatrix = uint8(listMatrix);
restMatrix = uint8(restMatrix);

% ----------------------------------------------------------------
% * criterion measures setup
% ----------------------------------------------------------------
if ~isempty(NumStimthresh) || ~isempty(maxCbalDevthresh) || ~isempty(maxFreqDevthresh)

	docriterion = 1;
	if isempty(NumStimthresh), NumStimthresh = 10000; end
 	if isempty(maxCbalDevthresh), maxCbalDevthresh = 10000; end
	if isempty(maxFreqDevthresh), maxFreqDevthresh = 10000; end

else
	docriterion = 0;
end

maxrestthresh = 2;					% max number of rests in a row.

% ----------------------------------------------------------------
% * hox gene setup
% ----------------------------------------------------------------
if numhox > 0
    hoxlist = hoxsetup(numhox,hoxrange,sizeGenerations);
end
if numhox > 0, disp('Hox range:'), hoxrange, end

% ----------------------------------------------------------------
% * Set up block contrast, if trans2block
% ----------------------------------------------------------------

if GA.trans2block && isfield(GA, 'blockcontrast') && GA.blockcontrast

    a = [ones(GA.restevery, 1); zeros(GA.restlength, 1)]; ab = [a; -a];
    blkstimlist = repmat(ab, ( ceil(2 * size(stimList, 1)) ./ length(ab) ), 1);
    blkmodel = sampleInSeconds(blkstimlist,ISI);
    blkmodel = resample(blkmodel,1,TR*10);
    blkmodel = blkmodel(1:numsamps, :);

end

% ----------------------------------------------------------------
% * save initial list matrix
% ----------------------------------------------------------------
if dorests
    save initial_designs listMatrix restMatrix contrasts svi
else
    save initial_designs listMatrix contrasts svi
end

% d-optimality: we need the whitening matrix, inv(svi), rather than svi
% i think...
if dflag, svi = pinv(svi);  end

% -------------------------------------------------------------------------------------------------
% * genetic algorithm
% ---------------------------------------------------------------------------------------------------

%try

for generation = 1:numGenerations+1

	% -------------------------------------------------------------------------------------------------
	% * set up rna (ready to model) lists and initialize generation variables
	% -------------------------------------------------------------------------------------------------

	gentime = cputime;

	rna = dna2rna(listMatrix,restMatrix,restlength,numStim,dorests,trans2switch,trans2block,dofirst);

	countnoqualify = 0;
    genMaxNumStim = 0;
    genMaxDev = 0;
    genMaxFreqDev = 0;
    genFreqMat = zeros(size(conditions));
    genmaxrest = 0;

	% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	% * WITHIN A GENERATION
 	% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	for z = 1:sizeGenerations % do this for each organism in the generation
		stimList = double(rna(:,z));

            % check that no conditions are missing
            if length(unique(stimList)) < nconds   % we have some missing ones; not a valid design
                fitnessMatrix(:, z) = 0;

                fprintf('Warning: Some designs are incomplete (too few event types). \nYou may have too many conditions for the length of time you''re using.');

                continue
            end

		% -------------------------------------------------------------------------------------------------
		% * counterbalancing
	    %  only performs counterbalancing up to length of conditions vector - so leaves out 1st trial in 'dofirst'.
      	% if (cbalColinPowerWeights(1) > 0)    removed 3/30 to implement criterion rather than weights.
		% -------------------------------------------------------------------------------------------------

		if (cbalColinPowerWeights(1) > 0) || (docriterion && maxCbalDevthresh < 10000)
   			[cBal, ~, maxDev] = getCounterBal(stimList, maxOrder,conditions,freqConditions);
    		if (cbalColinPowerWeights(1) > 0),fitnessMatrix(1,z) = cBal; end
		else
			maxDev = 0;				% satisfy criterion without testing if this criterion not specified.
		end

		% -------------------------------------------------------------------------------------------------
		% * discrepancy from desired frequencies
		% -------------------------------------------------------------------------------------------------

    	if (cbalColinPowerWeights(4) > 0) || docriterion
			[maxFreqDev,fitnessMatrix(4,z),freqMat] = calcFreqDev(stimList,conditions,freqConditions);
		elseif docriterion
			[maxFreqDev, ~, freqMat] = calcFreqDev(stimList,conditions,freqConditions);
		end

		% -------------------------------------------------------------------------------------------------
		% * hard-constraint criteria
		% -------------------------------------------------------------------------------------------------

	    go = true;

        % It's possible to get random variation that eliminates a
        % condition. Check for this and skip if so

        if length(unique(stimList(stimList ~= 0))) ~= nconds %size(contrasts, 2)
            go = false;
        end

        if docriterion
            [maxNumStim,maxrest] = getMaxInARow(stimList,dojitter);
            if maxNumStim > NumStimthresh || maxrest > maxrestthresh || maxDev > maxCbalDevthresh || maxFreqDev > maxFreqDevthresh
                go = false;
                countnoqualify = countnoqualify + 1;
                genMaxNumStim = genMaxNumStim + maxNumStim;
                genMaxDev = genMaxDev + maxDev;
                genMaxFreqDev = genMaxFreqDev + maxFreqDev;
                genFreqMat = genFreqMat + freqMat;
                genmaxrest = genmaxrest + maxrest;
            end
        end

		% -------------------------------------------------------------------------------------------------
		% * efficiency
		% -------------------------------------------------------------------------------------------------

    	if (cbalColinPowerWeights(2) > 0) && go   % build predictor set of vectors and convolve with HRF

			model = designvector2model(stimList,ISI,HRF,TR,numsamps,nonlinthreshold,S);

            % add Block effect for mixed block/ER
            % Kludgy insertion. tor: 5/6/08
            if GA.trans2block && isfield(GA, 'blockcontrast') && GA.blockcontrast

                testcontrasts = blkdiag(1, contrasts);
                testmodel = [blkmodel model];

                xtxitx = pinv(testmodel);   % a-optimality   % inv(X'S'SX)*(SX)'; pseudoinv of (S*X)
                fitnessMatrix(2,z) = calcEfficiency([1 contrastweights],testcontrasts, xtxitx, svi, dflag);

            else
                % Original code, used in most cases
                if dflag
                    xtxitx = model' * svi * model;  % d-optimality
                    fitnessMatrix(2,z) = calcEfficiency(contrastweights,contrasts,xtxitx,[],dflag);
                else
                    xtxitx = pinv(model);   % a-optimality   % inv(X'S'SX)*(SX)'; pseudoinv of (S*X)
                    fitnessMatrix(2,z) = calcEfficiency(contrastweights,contrasts,xtxitx,svi,dflag);
                end

            end

		elseif ~go    %(cbalColinPowerWeights(2) > 0 | cbalColinPowerWeights(3) > 0) & ~go

			fitnessMatrix(2,z) = -1000;
		else
			fitnessMatrix(2,z) = 0;
		end

        % -------------------------------------------------------------------------------------------------
		% * HRF shape estimation efficiency
		% -------------------------------------------------------------------------------------------------

    	if (cbalColinPowerWeights(3) > 0) && go   % build predictor set of vectors and convolve with HRF

            delta = [];
            for i = 1:max(stimList(:,1))
                delta(:,i) = (stimList == i);
            end

            % how much HRF estimation we do (up to 30 s) matters
            % we want at least as long as the expected actual response
			[model] = tor_make_deconv_mtx2(delta,round(30 / TR),TR / ISI);  % hard-coded for 30 s of HRF estimation

            % add Block effect for mixed block/ER
            % Kludgy insertion. tor: 5/6/08
            if GA.trans2block && isfield(GA, 'blockcontrast') && GA.blockcontrast
                model = [blkmodel model];
            end

            if ~isempty(S), model = S * model; end

            if dflag
                xtxitx = model'*svi*model;  % d-optimality
                fitnessMatrix(3,z) = calcEfficiency([],[],xtxitx,[],dflag);
            else
        	    xtxitx = pinv(model);   % a-optimality   % inv(X'S'SX)*(SX)'; pseudoinv of (S*X)
                fitnessMatrix(3,z) = calcEfficiency([],[],xtxitx,svi,dflag);
            end

		elseif 	~go    %(cbalColinPowerWeights(2) > 0 | cbalColinPowerWeights(3) > 0) & ~go
			fitnessMatrix(3,z) = -1000;
		else
			fitnessMatrix(3,z) = 0;
		end

   	end

   	clear model,clear maxNumStim,clear maxDev,clear maxFreqDev,clear cBal,clear dummy,clear stimList,clear go,clear xtxitx

	% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	% * AFTER EACH GENERATION
	% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	[overallFitness,mostFit,f] = getOverallFitness(cbalColinPowerWeights,fitnessMatrix);

   % Determine variability in population - rough measure
   % -------------------------------------------------------------------------------------------------
   popuv = abs(corrcoef(double(listMatrix)));
   popVar(generation) = 1 - mean(mean(popuv));
   clear popuv

   % save summary info
   % -------------------------------------------------------------------------------------------------
   %warning off
   bestCbal(generation) = f(1);
   bestVar(generation) = f(2);
   %warning on

   qualifyinglists(generation) = sizeGenerations - countnoqualify;

   if size(cbalColinPowerWeights,2) > 3, bestFreq(generation) = f(4); end
   bestHrf(generation) = f(3);
   clear f;

   % print summary of this generation to screen
   % -------------------------------------------------------------------------------------------------
   str = (['best of generation ' num2str(generation)]);
   if (cbalColinPowerWeights(1) > 0)
       str = [str ' / cbal. = ' num2str(bestCbal(generation))];
   end
   if (cbalColinPowerWeights(2) > 0)
        str = [str ' / Wtd. con eff = ' num2str(bestVar(generation))];
   end
   if (cbalColinPowerWeights(3) > 0)
        str = [str ' / Wtd. HRF eff = ' num2str(bestHrf(generation))];
   end
   if docriterion
        str = [str ' / ok lists = ' num2str(sizeGenerations - countnoqualify)];
   end
   str = [str ' / Pop var. = ' num2str(popVar(generation)) ' / time = ' num2str(cputime - gentime)];

   disp(str)

   if countnoqualify > sizeGenerations - 11
       disp('Warning: Most designs do not meet your hard constraint criteria. You are not likely to be able to')
       disp('optimize this design without removing some hard constraints!')
       disp(' ')
       fprintf('You specified a max of %01d stim of the same type in a row.  Average in designs: %3.1f', GA.NumStimthresh, genMaxNumStim/sizeGenerations)
       disp(['Avg. deviation from freqs = ' num2str(genMaxFreqDev/sizeGenerations)])
       disp(['Avg. deviation from counterbal = ' num2str(genMaxDev/sizeGenerations)])
       disp(['Avg. frequencies are = ' num2str(genFreqMat./sizeGenerations)])
       disp(['Avg. max cond 5 is = ' num2str(genmaxrest./sizeGenerations)])
       disp(' ')

   end

	% Save best of this generation
	% -------------------------------------------------------------------------------------------------
	if dorests
		myDesign = dna2rna(listMatrix(:,mostFit),restMatrix(:,mostFit),restlength,numStim,dorests,trans2switch,trans2block,dofirst);
	else
		myDesign = dna2rna(listMatrix(:,mostFit),[],restlength,numStim,dorests,trans2switch,trans2block,dofirst);
	end
	bestMat(:,generation) = myDesign(:,1);

    % -------------------------------------------------------------------------------------------------
   	% * crossbreeding
   	% -------------------------------------------------------------------------------------------------
   	if generation < numGenerations+1 && (cputime - clockStart) < maxTime
   		if dorests
		  	[listMatrix,restMatrix] = breedorgs(overallFitness,listMatrix,alph,restMatrix);
	  	else
		  	listMatrix = breedorgs(overallFitness,listMatrix,alph);
	  	end
 	end

    % ----------------------------------------------------------------
    % * check for zeros in stimlist
    % ----------------------------------------------------------------
   %  nzeros = any(sum(listMatrix == 0));

    % -------------------------------------------------------------------------------------------------
   	% * 10% random burst if pop var is low
   	% -------------------------------------------------------------------------------------------------
    if popVar(generation) < .40
        if trans2block
            listMatrix = tenPercentRandom(listMatrix,mostFit,conditions(1:floor(length(conditions)./2)));
        else
            listMatrix = tenPercentRandom(listMatrix,mostFit,conditions);
        end
    end

	% Save in-process results for recovery in case of crash
	% -------------------------------------------------------------------------------------------------

    if mod(generation,10) == 0
		save listMat_working
        % try, pack, catch, end
    end

	totalGenTime(generation) = cputime - clockStart;

	% Monitor time and break if over
	% -------------------------------------------------------------------------------------------------
	if (cputime - clockStart) > maxTime
		break
	end

	 % just test this to see if contrast variance is the same btwn modelDiagnostics and optimizeGA
	 % remove this later.  it's slow.
	 % --------------------------------------------------------------------------------------------------
	 %test = modelDiagnostics(bestMat(:,generation),'stimList','all',ISI,TR,contrasts,'noplot');
	 %disp(['	///modelDiagnostics: se    is ' num2str(test.sebeta) ' / eff is ' num2str(1./test.sebeta)])
	 %disp(['	///modelDiagnostics: conse is ' num2str(test.consebeta) ' / eff is ' num2str(1./test.consebeta)])
	 %disp(['	///modelDiagnostics: convar is ' num2str(test.consebeta.^2) ' / eff is ' num2str(1./test.consebeta.^2)])
	 %xtxitx = pinv(test.modelatTR);
	 %gamoddiagvar = xtxitx*svi*xtxitx'
	 %[fit,testlistconvarorig,xtxi] = testlist(bestMat(:,generation),GA,HRF)
	 %[fit,testlistconvarafterM,xtxi] = testlist(test.stimlist,GA,HRF)

end  % generations loop - finish ev. algorithm

%catch
%    dbstack
%    disp('Error in GA process.')
%    lasterr
%    keyboard
%end

% AFTER ALL GENERATIONS
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %

disp('saving workspace as GAresults.mat')
save GAresults

M.origlist = listMatrix(:,mostFit);
if dorests,M.restlist = restMatrix(:,mostFit); end
M.stimlist = myDesign;
M.ga = GA;
M.ga.HRF = HRF;
M.ga.bestCbal = bestCbal;
M.ga.bestCVar = bestVar;
M.ga.bestHrf = bestHrf;
if size(cbalColinPowerWeights,2) > 3,M.ga.bestFreq = bestFreq; end
M.ga.popVar = popVar;
M.ga.totalGenTime = totalGenTime;

M.ga.qualifyinglists = qualifyinglists;
M.ga.numStimEachCond = numStimEachCond;
M.ga.numGenerations = numGenerations;
M.ga.sizeGenerations = sizeGenerations;
M.ga.nonlinthreshold = nonlinthreshold;
M.ga.listMatrix = listMatrix;
M.bestLists = bestMat;
date = clock;
M.date = [num2str(date(2)) '-' num2str(date(3)) '-' num2str(date(1))];

try

	if plotFlag
    	M = modelDiagnostics2(M,'screen', 'HRF', HRF);
	else
    	M = modelDiagnostics2(M,'noplot', 'HRF', HRF);
    end

catch

	disp('Error in modelDiagnostics2')
	lasterr
	keyboard

end

% PLOT additional OUTPUT
% -------------------- %
if plotFlag

    create_figure('GA_over_generations', 6, 1);

    subplot(6,1,1)
    plot(bestCbal, 'g', 'LineWidth', 2)
    title('Counterbalancing of best design in each generation')

    subplot(6,1,2)
    plot(bestVar, 'b', 'LineWidth', 2)
    title('Contrast efficiency for best design in each generation')

    subplot(6,1,3)
    plot(bestHrf, 'r', 'LineWidth', 2)
    title('HRF efficiency for best design in each generation')

    subplot(6,1,4)
    plot(qualifyinglists, 'k')
    title('Number of lists satisfying criteria in each generation')

    subplot(6,1,5)
    plot(bestFreq, 'k')
    title('Max % discrepancy between actual and specified frequencies')

    subplot(6,1,6)
    plot(popVar, 'm')
    title('Population variability - 1 is uncorrelated columns, 0 is all identical organisms')

%     set(gcf,'Position',[232   258   560   420]); % [612    55   665   889]);
    format short g
    format compact
    disp('Elapsed Time for this run')
    elapsedTime = clock - clockStart

end

disp(['optimizeGA done: ' num2str(cputime-clockStart) ' s'])
return

% ====== Sub-functions =================================================================

function outList = transform2switches(stimList)
% assumes 2 stimtypes, 1 and 2
% classifies them into no switch (1 and 2) or switch (3 and 4)

outList(1,1) = stimList(1,1);
for i = 2:size(stimList,1)
    switch stimList(i-1,1)		% based on previous trial

        case 1
            if stimList(i,1) == 2, outList(i,1) = 4; else outList(i,1) = stimList(i,1); end
        case 2
            if stimList(i,1) == 1, outList(i,1) = 3; else outList(i,1) = stimList(i,1); end
        case 5
            if i > 2 && stimList(i,1) == 1 && stimList(i-2,1) == 2, outList(i,1) = 3;
            elseif i > 2 && stimList(i,1) == 2 && stimList(i-2,1) == 1, outList(i,1) = 4;
            else outList(i,1) = stimList(i,1);
            end
        otherwise outList(i,1) = stimList(i,1);
    end
end

return

% USE INDEPENDENT FUNCTION TRANSFORM2BLOCK
%function sl = transform2block(sl,rl,restlength)
% sl and rl are stimlist and restlist.
% now is NOT compatible with a separate jitter blank rest interval.

%maxnum = max(sl);
%if doblank, maxnum = maxnum-1;,end

%for i = 2:2:length(rl)
%     start = sum(rl(1:i-1)) + restlength * length(rl(1:i-1)) + 1;
%     stop = min(start + rl(i) - 1,length(sl));
%     sl(start:stop) = sl(start:stop) + maxnum;
%end

%return

function [max,maxrest] = getMaxInARow(stimList,dojitter)

max = 1; counter = 1;maxrest = 1;counterrest = 1;
if dojitter, restnum = max(stimList); else restnum = Inf; end
for i = 2:size(stimList,1)
    if stimList(i,1) == stimList(i-1,1)&& ~stimList(i,1) == 0
        counter = counter + 1;
        if counter > max, max = counter; end
    else counter = 1;
    end

    if dojitter
        if stimList(i,1) == restnum && stimList(i-1,1) == restnum
            counterrest = counterrest+1;
            if counterrest > maxrest,maxrest = counterrest; end
        else counterrest = 1;
        end
    end

end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function fcnn = fcnnbp(fcnn,Rsubset)
%FCNNBP Backpropagation in fully-connected neural net.
%   FCN = FCNNBP(FCNN,Rsubset) backropagation through fully-conected
%   neural network, FCNN. Rsubset is a subset of datain.R determined in
%   training based on the size of the minibatch. The following code
%   implements Step 3 in Table 14.8.
%
%   Type
%
%   >> help fcninfo
%
%   at the prompt for detailed explanations of the components of the
%   fully-connected neural net.
%
%   This function computes the backpropagation of the output error. The
%   result is fcnn(k).D, for k = L-1:-1:2. An extra backprop step is
%   taken to generate fcnn(1).D. The extra backprop step brings the
%   backprop output (going in the reverse direction) to the same
%   location as the feedforward output of the fcnn after vectorization.
%   The extra step does not affect any other fcnn functions.
%
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Number of layers.
L = numel([fcnn.NumNodes]);

% Compute fcnn(L).D.
fcnn(L).D = (fcnn(L).A - Rsubset).*fcnn(L).HprimeZ;

% Backpropagation.
for k = L-1:-1:2
   fcnn(k).D = (fcnn(k + 1).Weights'*fcnn(k + 1).D).*fcnn(k).HprimeZ;
end

% Backprop one more layer (see comments in the help section). Note that
% there is no fcnn(1).HprimeZ because this quantity is related only to
% the fully connected net. This extra step does not affect any other fcn
% functions.
fcnn(1).D = fcnn(2).Weights'*fcnn(2).D;

author: dipum

repo: dipum-toolbox

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


text:
function [F, inls] = rEG(u, th, th7, conf)

%rEG     Robust computation of epipolar geometry based on RANSAC
%
% [F, inls] = rEG(u, th, th7, conf)
% u    point correspondences (6xn), where n is the number of corrs.
% th   threshold value for the Sampson's distance (see FDs)
% th7  threshold for inliers to iterate on F (default = th)
% conf confidence level of self-termination  (default = .95)

MAX_SAM = 100000;
iter_amount = .5;

if nargin < 3
  th7 = th;
end;

if nargin < 4
  conf = .95;
end;

len = size(u,2);
ptr = 1:len;
max_i = 8;
max_sam = MAX_SAM;

no_sam = 0;
no_mod = 0;

while no_sam < max_sam
  for pos = 1:7
      idx = pos + ceil(rand * (len-pos));
      ptr([pos, idx]) = ptr([idx, pos]);
  end;

  no_sam = no_sam +1;

  aFs = fu2F7(u(:,ptr(1:7)));

  for i = 1:size(aFs,3)
      no_mod = no_mod +1;
      aF = aFs(:,:,i);
      Ds = fFDs(aF,u);
      v  = Ds < th;
      v7 = Ds < th7;
      no_i  = sum(v);

      if max_i < no_i
        inls = v;
        F = aF;
        max_i = no_i;
        max_sam = min([max_sam,nsamples(max_i, len, 7, conf)]);
      end

      if sum(v7) >= 8 + iter_amount*(max_i - 8)
	aF = u2F(u(:,v7));
        Ds = fFDs(aF,u);
        v  = Ds < th;
	no_i = sum(v);
        if max_i < no_i
          inls = v;
          F = aF;
          max_i = no_i;
	  exp_sam = nsamples(max_i, len, 7, .95);
          max_sam = min([max_sam,exp_sam]);
        end
      end
  end
end

if no_sam == MAX_SAM
  warning(sprintf('RANSAC - termination forced after %d samples expected number of samples is %d',  no_sam, exp_sam));
end;

author: strawlab

repo: MultiCamSelfCal


text:
function coeffs=quadEllipsForm2Poly(A,x0,c)
%%QUADELLIPSFORM2POLY Given an ellipse (in 2D) or an ellipsoid (in >2D)
%          given by the equation such that a point x is in or on the
%          ellipsoid if (x-x0)'*A*(x-x0)<=c, convert it into a set of
%          coefficients for a quadratic multivariate polynomial such
%          that if a point is on or in the ellipsoid, then
%          f(x1,x2,x3,...)<=0, where f is the polynomial and the x's are
%          all of the variables. These coefficients are returned in an
%          ordering such that one can use polyValMultiDim to evaluate the
%          polynomial.
%
%INPUTS: A A numDimXnumDim real matrix that specifies the shape and
%          rotation of the ellipsoid. This doesn't actually have to be
%          symmetric or positive definite. numDim>=1.
%       x0 The numDimX1 real center of the ellipsoid.
%        c The real, scalar right-hand side of the ellipsoid equation,
%          which affects the scale of the ellipsoid.
%
%OUTPUTS: coeffs A 3X3X...X3 hypermatrix of the coefficients for the
%               multivariate quadratic polynomial. These are arranged such
%               that coeffs(a1,a2,a3...an) corresponds to the coefficient
%               of an x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1) term.
%               Thus, the number of indices coeffs takes is equal to the
%               dimensionality of x (not counting singleton dimensions at
%               the end of coeffs).
%
%The solution was obtained by just multiplying out (x-x0)'*A*(x-x0)-c and
%grouping terms. The inverse of this function is polyEllipsForm2Quad, which
%returns a symmetric A.
%
%EXAMPLE:
%In this example, we generate random values for A, x0, and c and change the
%ellipse format. For a random test point xTest, we then verify that
%(xTest-x0)'*A*(xTest-x0)-c produces the same value as using coeffs and
%xTest in polyValMultiDim. The relative error is shown and it is on the
%order of finite precision limitations.
% numDim=3;
% A=randn(numDim,numDim);
% x0=randn(numDim,1);
% c=randn(1);
% coeffs=quadEllipsForm2Poly(A,x0,c);
% xTest=randn(numDim,1);
% val1=(xTest-x0)'*A*(xTest-x0)-c;
% val2=polyValMultiDim(coeffs,xTest);
% RelDiff=(val1-val2)/val1
%
%January 2023 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Make A symmetric. This doesn't change the value (x-x0)'*A*(x-x0).
A=(A+A')/2;

numDim=size(A,1);

if(numDim==1)
    coeffs=zeros(3,1);
    coeffs(1)=x0^2*A-c;
    coeffs(2)=-2*x0*A;
    coeffs(3)=A;

elseif(numDim==2)
    coeffs=zeros(3,3);

    %The constant term.
    coeffs(1,1)=x0'*A*x0-c;

    %The linear terms.
    linTerms=-2*A*x0;%This works if A=A' --which we ensured above.
    coeffs(2,1)=linTerms(1);
    coeffs(1,2)=linTerms(2);

    %The squared terms.
    coeffs(3,1)=A(1,1);
    coeffs(1,3)=A(2,2);

    %The cross term.
    coeffs(2,2)=2*A(1,2);
elseif(numDim==3)
    coeffs=zeros(3,3,3);

    %The constant term.
    coeffs(1,1,1)=x0'*A*x0-c;

    %The linear terms.
    linTerms=-2*A*x0;
    coeffs(2,1,1)=linTerms(1);
    coeffs(1,2,1)=linTerms(2);
    coeffs(1,1,2)=linTerms(3);

    %The squared terms.
    coeffs(3,1,1)=A(1,1);
    coeffs(1,3,1)=A(2,2);
    coeffs(1,1,3)=A(3,3);

    %The cross terms.
    coeffs(2,2,1)=2*A(1,2);
    coeffs(2,1,2)=2*A(1,3);
    coeffs(1,2,2)=2*A(2,3);
else
    dims=3*ones(1,numDim);
    coeffs=zeros(dims);

    %The constant term.
    coeffs(1)=x0'*A*x0-c;

    idx=repmat({1},[1,numDim]);

    %The linear terms.
    linTerms=-2*A*x0;
    for curDim=1:numDim
        idx{curDim}=2;
        coeffs(idx{:})=linTerms(curDim);
        idx{curDim}=1;
    end

    %The squared terms.
    for curDim=1:numDim
        idx{curDim}=3;
        coeffs(idx{:})=A(curDim,curDim);
        idx{curDim}=1;
    end

    %The linear terms
    for curDim1=1:(numDim-1)
        idx{curDim1}=2;
        for curDim2=(curDim1+1):numDim
            idx{curDim2}=2;
            coeffs(idx{:})=2*A(curDim1,curDim2);
            idx{curDim2}=1;
        end
        idx{curDim1}=1;
    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 this = subsasgn(this,subs,dat)
% Overloaded subsasgn function for meeg objects.
% _________________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: subsasgn.m 2716 2009-02-09 17:14:15Z vladimir $

if isempty(subs)
    return;
end;

if this.Nsamples == 0
    error('Attempt to assign to a field of an empty meeg object.');
end;

if strcmp(subs(1).type, '.')
    if ismethod(this, subs(1).subs)
        error('meeg method names cannot be used for custom fields');
    else
        if isempty(this.other)
            this.other = struct(subs(1).subs, {dat});
        else
            this.other = builtin('subsasgn', this.other, subs, dat);
        end
    end
elseif strcmp(subs(1).type, '()')
    if numel(subs)~= 1, error('Expression too complicated');end;
    this.data.y = subsasgn(this.data.y, subs, dat);
else
    error('Unsupported assignment type for meeg.');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%filename = 'CantileverSquareSYmmetricMesh';
%filename = 'CantileverSquareMedium';
%'CantileverSquareSYmmetricMesh';
%filename = 'CantileverSquare';
%filename = 'CantileverSquareSmall';
%filename = 'CantileverSquareNew';
%'CantileverSquareNewFine';
%filename = 'Cantilever_quad_coarse';
%filename = 'LshapeTriFine';
%filename = 'LshapeTri';
filename = 'LshapeTriSmall';
%filename = 'Lshape';
%filename = 'LshapeFine';
%filename = 'ArchTriFine';
%'Arch_quad_coarse';
%'Bridge_quad_coarse';
%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';
%filename = 'Bridge';

ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'PDE';
%filterType = 'P1';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-4;
constr_initial = 1e-8;

Vfrac_final = 0.3;
optimality_final = 1e-4;
constr_final = 1e-8;

stressNormExponent_initial = 2;
stressNormExponent_final = 16;
%
optimizer = 'DualNestedInPrimal';
%optimizer = 'AlternatingPrimalDual';
optimizerUnconstrained = 'PROJECTED GRADIENT';

%optimizer = 'MMA';
%optimizer = 'IPOPT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
%vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
vademecumFileName = 'SuperEllipseQOptAnalytic';

% designVariable = 'LevelSet';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';
% initial_case = 'full';% optimizerUnconstrained = 'SLERP';

% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';
% initial_case = 'full';
% rho0 = 0.3;

line_search_initiator = 'INCREASING LAST STEP';
incrementFactor = 1.95;
%

%kfrac = 2;
nsteps = 8;%17;

plotting = true;
printing = true;
monitoring = true;
monitoring_interval = 3;
maxiter = 8000;

%
% % %
% isDirichletPartX = @(x) x > -1e-12 & x < 0.1;
% isDirichletPart1 = @(y) y > 0.20 & y < 0.30;
% isDirichletPart2 = @(y) y > 0.70 & y < 0.80;
% isDirichletPartY = @(y) isDirichletPart1(y) | isDirichletPart2(y);
% isDirichletPart = @(x,y) isDirichletPartX(x) & isDirichletPartY(y);
 isNeumannPartX = @(x) x > (1-0.05) & x < (1+1e-12);
 isNeumannPartY = @(y) y > 0.35 & y < 0.4;
%
% isNeumannPart = @(x,y) isNeumannPartX(x) & isNeumannPartY(y);
% iNotOptimizable = @(coord) isDirichletPart(coord(:,1),coord(:,2)) & ~isNeumannPart(coord(:,1),coord(:,2));
%
% costDomainNotOptimizable       = {iNotOptimizable};
% constraintDomainNotOptimizable = {[]};
%
% isDesignVariableFixed.nodes  = iNotOptimizable;
% isDesignVariableFixed.values = @(x) [m1*ones(size(x,1),1);m2*ones(size(x,1),1)];

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [nrm] = norm(tt)
%Frobenius norm of the QTT-Tucker
%   [NRM]=NORM(QTT_TUCKER) Computes the Frobenius norm of the QTT-Tucker
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
d=tt.dphys;
core=tt.core;
tuck=tt.tuck;
for i=1:d
    if (isa(tuck{i}, 'tt_matrix'))
        tuck{i} = tt_tensor(tuck{i});
    end;
   [tuck{i},rm]=qr(tuck{i},'lr');
   core{i}=ten_conv(core{i},2,rm.');
end
nrm=norm(core);

return
end

author: oseledets

repo: TT-Toolbox


text:
function y = blkdiag(varargin)
% Allows char input
ind = cellfun(@ischar, varargin);
for k = find(ind)
    varargin{k} = double(varargin{k});
end
y = builtin('blkdiag', varargin{:});
if any(ind)
    y = char(y);
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function [R O]=Forward_SCARA(theta,d3,d4,a1,a2)

th1=theta(1,1);
th2=theta(1,2);
th4=theta(1,3);
alpha=th1+th2-th4;

R=[cos(alpha),sin(alpha),0;sin(alpha),-cos(alpha),0;0,0,-1];
O=[((a1*cos(th1))+(a2*cos(th1+th2)));((a1*sin(th1))+(a2*sin(th1+th2)));(-d3-d4)];

author: Sable

repo: mcbench-benchmarks


text:
function [f,G] = tt_cp_wfg(Z,W,A,normZsqr)
%TT_CP_WFG Function and gradient of CP with missing data.
%
%   [F,G] = TT_CP_WFG(Z,W,A) computes the function and gradient values of
%   the function 0.5 * || W .* (Z - ktensor(A)) ||^2. The input A is a
%   cell array containing the factor matrices. The input W is a (dense
%   or sparse) tensor containing zeros wherever data is missing. The
%   input Z is a (dense or sparse) tensor that is assumed to have
%   zeros wherever there is missing data. The output is the function F
%   and a cell array G containing the partial derivatives with respect
%   to the factor matrices.
%
%   [F,G] = TT_CP_WFG(Z,W,A,NORMZSQR) also passes in the pre-computed
%   norm of Z, which makes the computations faster.
%
%   See also TT_CP_WFUN, TT_CP_WFG_SPARSE, TT_CP_WFG_SPARSE_SETUP.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Compute B = W.*ktensor(A)
if isa(W,'sptensor')
    B = W.*ktensor(A);
else
    B = W.*full(ktensor(A));
end

%% Compute normZ
if ~exist('normZsqr','var')
    normZsqr = norm(Z)^2;
end

% function value
f = 0.5 * normZsqr - innerprod(Z,B) + 0.5 * norm(B)^2;

% gradient computation
N = ndims(Z);
G = cell(N,1);
T = Z - B;
for n = 1:N
    G{n} = zeros(size(A{n}));
    G{n} = -mttkrp(T,A,n);
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function nop(varargin)
% this is a function that does nothing and returns nothing.
end

author: kristinbranson

repo: JAABA


function errorUnknownEnvironment()
% Throw an error to indicate an unknwon environment (i.e. not MATLAB/Octave).
    error('matlab2tikz:unknownEnvironment',...
          'Unknown environment "%s". Need MATLAB(R) or Octave.', getEnvironment);
end


text:
%SF_QUAD_Q3 Bicubic conforming shape function for quadrilaterals (Q3).
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SF_QUAD_Q3( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates conforming bicubic Q3 shape functions on quadrilaterals
%   with values defined in the nodes, edges, and cell center. XI is
%   [-1..1]^2 reference coordinates.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar:  2             Number of space dimensions
%       n_vert      scalar:  4             Number of vertices per cell
%       i_dof       scalar: 1-16           Local basis function to evaluate
%       xi          [n_sdim]               Local coordinates of evaluation point
%       aInvJac     [n,n_sdim*n_sdim]      Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [4]                    Number of local degrees of freedom on
%                                          vertices, edges, faces, and cell interiors
%       xLDof       [n_sdim,n_ldof]        Local coordinates of local dofs
%       sfun        string                 Function name of called shape function
%
%   See also SFLAG3, SF_QUAD_H3

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function test_suite=test_moxunit_util_isfolder
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_moxunit_util_isfolder_basics()
    % current direcoty
    assertTrue(moxunit_util_isfolder(pwd));
    assertTrue(moxunit_util_isfolder('.'));
    assertTrue(moxunit_util_isfolder('..'));

    % this test fails almost surely
    r = rand_str();
    assertFalse(moxunit_util_isfolder(r));

    % illegal input type
    assertExceptionThrown(@()moxunit_util_isfolder(struct));

function s=rand_str()
    n = 25;
    s=char(26*rand(1,n)+65);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [K,dK] = covRQard_fast(hyp, x, z, i)

% Rational Quadratic covariance function with Automatic Relevance Determination
% (ARD) distance measure. The covariance function is parameterized as:
%
% k(x^p,x^q) = sf^2 * [1 + (x^p - x^q)'*inv(P)*(x^p - x^q)/(2*alpha)]^(-alpha)
%
% where the P matrix is diagonal with ARD parameters ell_1^2,...,ell_D^2, where
% D is the dimension of the input space, sf2 is the signal variance and alpha
% is the shape parameter for the RQ covariance. The hyperparameters are:
%
% hyp = [ log(ell_1)
%         log(ell_2)
%          ..
%         log(ell_D)
%         log(sf)
%         log(alpha) ]
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2010-08-04.
%
% See also COVFUNCTIONS.M.

if nargin<2, K = '(D+2)'; return; end              % report number of parameters
if nargin<3, z = []; end                                   % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

[n,D] = size(x);
ell = exp(hyp(1:D));
sf2 = exp(2*hyp(D+1));
alpha = exp(hyp(D+2));

% precompute squared distances
if dg                                                               % vector kxx
  D2 = zeros(size(x,1),1);
else
  if xeqz                                                 % symmetric matrix Kxx
    D2 = sq_dist(diag(1./ell)*x');
  else                                                   % cross covariances Kxz
    D2 = sq_dist(diag(1./ell)*x',diag(1./ell)*z');
  end
end

if nargout < 2

    if nargin<4                                                        % covariances
      K = sf2*(1+0.5*D2/alpha).^(-alpha);
    else                                                               % derivatives
      if i<=D                                               % length scale parameter
        if dg
          K = D2*0;
        else
          if xeqz
            K = sf2*(1+0.5*D2/alpha).^(-alpha-1).*sq_dist(x(:,i)'/ell(i));
          else
            K = sf2*(1+0.5*D2/alpha).^(-alpha-1).*sq_dist(x(:,i)'/ell(i),z(:,i)'/ell(i));
          end
        end
      elseif i==D+1                                            % magnitude parameter
        K = 2*sf2*(1+0.5*D2/alpha).^(-alpha);
      elseif i==D+2
        K = (1+0.5*D2/alpha);
        K = sf2*K.^(-alpha).*(0.5*D2./K - alpha*log(K));
      else
        error('Unknown hyperparameter')
      end
    end

else
  M = 1+0.5*D2/alpha;
  K = sf2*M.^(-alpha);                              % covariances
  dK = zeros(size(D2,1),size(D2,2),D+2);            % prepare derivatives

  if ~dg                                            % length scale parameters
    if xeqz
       C = sq_dist_fast(x,[],ell);
    else
       error('Not supported');
       % dK(:,:,i) = sf2*M.^(-alpha-1).*sq_dist(x(:,i)'/ell(i),z(:,i)'/ell(i));
       b(:,1,:) = z;
    end
    dK(:,:,1:D) = bsxfun(@times,K./M,C);
  end

  dK(:,:,D+1) = 2*K;                                            % magnitude parameter
  dK(:,:,D+2) = K.*(0.5*D2./M - alpha*log(M));    % shape parameter
end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
% The COBRAToolbox: testRemoveDuplicateRxn.m
%
% Purpose:
%     - testRemoveDuplicateRxn tests removeDuplicateRxn
%
% Authors:
%     - CI integration: Laurent Heirendt January 2017

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testRemoveDuplicateRxn'));
cd(fileDir);

% test detection and removal of duplicate reactions
model.S = [-1, 0, 0 ,0 , 0, 0, 0;
            1, -1, 0, 0, 0, 0, 0;
            0, -1, 0,-1, 0, 0, 0;
            0, 1, 0, 1, 0, 0, 0;
            0, 1, 0, 1, 0, 0, 0;
            0, 1,-1, 0, 0, 0, 0;
            0, 0, 1,-1, 1, 0, 0;
            0, 0, 0, 1,-1,-1, 0;
            0, 0, 0, 0, 1, 0, 0;
            0, 0, 0, 0,-1, 0, 0;
            0, 0, 0, 0, 0, 1, 1;
            0, 0, 0, 0, 0, 1, -1];
model.lb = [0, 0, 0, 0, 0, 0, 0];
model.ub = [20, 20, 20, 20, 20, 20, 20];
model.rxns = {'GLCt1'; 'HEX1'; 'PGI'; 'PFK'; 'FBP'; 'FBA'; 'TPI'};
model.mets = {'glc-D[e]'; 'glc-D'; 'atp'; 'H'; 'adp'; 'g6p';'f6p'; 'fdp'; 'pi'; 'h2o'; 'g3p'; 'dhap'};
mets_length = length(model.mets);
rxns_length = length(model.rxns);

sc =  [1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];

% adding a reaction to the model
%[model,rxnIDexists] = addReaction(model,rxnName,metaboliteList,stoichCoeffList,revFlag,lowerBound,upperBound,objCoeff,subSystem,grRule,geneNameList,systNameList,checkDuplicate,printLevel)
model = addReaction(model, 'GLCt1_duplicate_reverse', model.mets, sc, 0, 0, 20,0,'temp', [], [], [], 1, 0);

method = 'FR'; %will be removed as detects reverse reaction
printLevel = 1;
removeFlag = 1;
[model, removedRxn, rxnRelationship] = checkDuplicateRxn(model, method, removeFlag, printLevel);

assert(rxns_length == length(model.rxns));

% adding a reaction to the model
model = addReaction(model, 'GLCt1_duplicate_reverse', model.mets, sc, 0, 0, 20,0,'temp', [], [], [], 1, 0);

method = 'S';%will not be removed as does not detect reverse reaction
printLevel = 1;
removeFlag = 1;
[model, removedRxn, rxnRelationship] = checkDuplicateRxn(model, method, removeFlag, printLevel);

assert(rxns_length + 1 == length(model.rxns));

%return to original directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [Ioutput]= marrfilter(I)
I1=im2double(I);
[Izerosmooth]= filtering(I1);
[Ioutput]=zerocrossing(Izerosmooth);
Ioutput1=im2uint8(Ioutput);
Ioutput1= Ioutput1>105;
imshow(Ioutput1);

figure;
  subplot(2,2,1);imshow(I);title('origional image');
  subplot(2,2,2);imshow(Izerosmooth);title('Zerosmooth image');
   subplot(2,2,3);imshow(Ioutput);title('output image');
   subplot(2,2,4);imshow(Ioutput1);title('output image with threshold');


text:
function [SOI,SOIObjX] = SelectSOI(Population,PopObjX,N,z)
% Select the solutions of interest

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

    %% Calculate angles between solutions and ASF
    Angle = acos(1-pdist2(PopObjX,PopObjX,'cosine'));
    Angle(logical(eye(length(Angle)))) = inf;
    ASF = max((PopObjX-z).*sum(PopObjX,2)./PopObjX,[],2);

    %% Delete solutions
    Remain = 1 : size(PopObjX,1);
    while length(Remain) > N
        [Dis,x] = min(Angle(Remain,Remain),[],2);
        [~,y]   = min(Dis);
        x = x(y);
        if ASF(Remain(x)) > ASF(Remain(y))
            Remain(x) = [];
        else
            Remain(y) = [];
        end
    end
    SOI     = Population(Remain);
    SOIObjX = PopObjX(Remain,:);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function c = midrad(a,r)
%MIDRAD       Initialization of interval by midpoint and radius
%  computed such that <a,r> is enclosed in interval  c
%
%   c = midrad(a,r)
%

% written  10/16/98     S.M. Rump
% modified 06/24/99     S.M. Rump  check sparsity, multi-dimensional arrays
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 11/16/02     S.M. Rump  sparse radius
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 09/10/07     S.M. Rump  performance
% modified 10/14/08     S.M. Rump  huge midpoint array and 0 radius
% modified 10/14/08     S.M. Rump  transform [inf,inf] into [realmax,inf] or <0,inf>
% modified 02/28/12     S.M. Rump  check input is numeric (thanks to Joshua Ross, Adelaide)
% modified 07/11/12     S.M. Rump  typo
%

  if ~isnumeric(a) | ~isnumeric(r)
    error('both bounds for infsup must be numeric')
  end

  index = any( r<0 );                   % take care of huge matrices
  if ~isreal(r) | any(index(:))
    error('invalid radius in call of midrad')
  end

  if ( prod(size(r))~=1) & ( issparse(a)~=issparse(r) )
    error('midpoint and radius must both be or none of them be sparse')
  end

  if isreal(a)                     % real interval (inf/sup representation)
    e = 1e-30;
    if 1+e==1-e                    % fast check for rounding to nearest
      rndold = 0;
    else
      rndold = getround;
    end
    [m,n] = size(a);                % already proper intervals
    if issparse(a) & isequal(r,0)
      r = sparse(0);
    end
    setround(-1)
    cinf = a - r;
    setround(1)
    csup = a + r;
    setround(rndold)                % set rounding to previous value
    c = intval(cinf,csup,'infsup');
  else                             % complex interval (mid/rad representation)
    c = intval(a,r,'midrad');
  end

repo: OpenMAS


text:
function [w, infos] = cg(problem, in_options)
% Linear conjugate gradient (CG) algorithm.
%
% The algorithm of interest is defined as
%
%           min f(x) = 1/2 * x^T * A * x - b^T * x.
%           where
%           x in R^d
%
% Inputs:
%       problem     function (cost/grad/hess)
%       options     options
% Output:
%       w           solution of w
%       infos       information
%
% Reference:
%       Jorge Nocedal and Stephen Wright,
%       "Numerical optimization,"
%       Springer Science & Business Media, 2006.
%
%       Algorithm 5.1, 5.2 in Section 5.1.
%
%       sub_mode    'PRELIM'    Algorithm 5.1
%       sub_mode    'STANDARD'  Algorithm 5.2
%       sub_mode    'PRECON'    Algorithm 5.3
%
% This file is part of GDLibrary and SGDLibrary.
%
% Created by H.Kasai on Oct. 30, 2016
% Modified by H.Kasai on Mar. 25, 2018
% Modified by H.Kasai on Oct. 22, 2020

    % set dimensions and samples
    d = problem.dim;
    n = problem.samples;

    % set local options
    local_options = [];
    local_options.algorithm = 'CG';
    local_options.sub_mode = 'STANDARD';

    % merge options
    options = mergeOptions(get_default_options(d), local_options);
    options = mergeOptions(options, in_options);

    if ~isfield(options, 'M')
        if strcmp(options.sub_mode, 'PRECON')
            h = problem.full_hess(w);
            M = diag(diag(h));
        else
            M = eye(d);
        end
    else
        M = options.M;
    end

%     ls_options.sub_mode = sub_mode;
%     if strcmp(sub_mode, 'STANDARD')
%         %
%     elseif strcmp(optionssub_mode, 'PRECON')
%         ls_options.M = M;
%     else
%         %
%     end

    % initialise
    iter = 0;
    grad_calc_count = 0;
    w = options.w_init;

   % initialize by BB step-size
    if strcmp(options.step_init_alg, 'bb_init')
        options.step_init = bb_init(problem, w);
    end

    % store first infos
    clear infos;
    [infos, f_val, optgap, grad, gnorm] = store_infos(problem, w, options, [], iter, grad_calc_count, 0);
    grad_old = [];

    % display infos
    if options.verbose
        fprintf('CG (%s): Iter = %03d, cost = %.16e, gnorm = %.4e, optgap = %.4e\n', options.sub_mode, iter, f_val, gnorm, optgap);
    end

    % set first residual, i.e., r = Ax - b;
    r_old = grad;
    y_old = [];
    w_old = w;
    grad_old = grad;
    prev_step = options.step_init;

    % initialise
    if strcmp(options.sub_mode, 'PRELIM')
        % set directoin
        p = -r_old;
    elseif strcmp(options.sub_mode, 'STANDARD')
        % set directoin
        p = -r_old;
    elseif strcmp(options.sub_mode, 'PRECON')
        % solve y
        y_old = M \r_old;
        % set directoin
        p = -y_old;
    else
        fprintf('sub_mode %s is not supported\n', sub_mode);
        return;
    end

    % set start time
    start_time = tic();

    % main loop
    while (optgap > options.tol_optgap) && (gnorm > options.tol_gnorm) && (iter < options.max_epoch)

        % Revert to steepest descent if is not direction of descent
        if (p'*grad > 0)
            p = -p;
        end

        % line search
        options.r_old = r_old;
        options.y_old = y_old;
        [step, ~] = options.linesearchfun(options.step_alg, problem, w, w_old, p, grad_old, prev_step, options);
        prev_step = step;

        % update w
        w_old = w;
        w = w + step * p;

        % proximal operator
        if ismethod(problem, 'prox')
            w = problem.prox(w, step);
        end

        % calculate gradient
        grad_old = grad;
        grad = problem.full_grad(w);

        if strcmp(options.sub_mode, 'PRELIM')
            % updata residual
            r = grad;

            % updata beta
            beta = r' * problem.A() * p / (p' * problem.A() * p);

            % update direction
            p = - grad + beta * p;

        elseif  strcmp(options.sub_mode, 'STANDARD')
            % updata residual
            r = r_old + step * problem.A() * p;

            % updata beta
            beta = r' * r / (r_old' * r_old);

            % updata beta
            p = - r + beta * p;

            % store r
            r_old = r;
        else
            % updata residual
            r = r_old + step * problem.A() * p;

            % solve y
            y = M \r;

            % updata beta
            beta = r' * y / (r_old' * y_old);

            % updata beta
            p = -y + beta * p;

            % store r
            r_old = r;
            y_old = y;
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % count gradient evaluations
        grad_calc_count = grad_calc_count + n;

        % update iter
        iter = iter + 1;

        % store infos
        [infos, f_val, optgap, grad, gnorm] = store_infos(problem, w, options, infos, iter, grad_calc_count, elapsed_time);

        % display infos
        if options.verbose
            fprintf('CG (%s): Iter = %03d, cost = %.16e, gnorm = %.4e, optgap = %.4e\n', options.sub_mode, iter, f_val, gnorm, optgap);
        end
    end

    if gnorm < options.tol_gnorm
        fprintf('Gradient norm tolerance reached: tol_gnorm = %g\n', options.tol_gnorm);
    elseif optgap < options.tol_optgap
        fprintf('Optimality gap tolerance reached: tol_optgap = %g\n', options.tol_optgap);
    elseif iter == options.max_epoch
        fprintf('Max iter reached: max_epoch = %g\n', options.max_epoch);
    end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function sobol_test ( )

%*****************************************************************************80
%
%% SOBOL_TEST tests the SOBOL library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SOBOL_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SOBOL library.\n' );

  sobol_test01 ( );
  sobol_test02 ( );
  sobol_test03 ( );
  sobol_test04 ( );
  sobol_test05 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SOBOL_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% This file is part of the following project:
% Oliver Parson, Siddhartha Ghosh, Mark Weal, Alex Rogers.
% Non-intrusive Load Monitoring using Prior Models of General Appliance Types.
% In: 26th AAAI Conference on Artificial Intelligence. Toronto, Canada. 2012.
% Code available for download: https://sites.google.com/site/oliparson/phd-work/research-files/aaai-2012-code.zip?attredirects=0
% Copyright: Oliver Parson et al., University of Southhampton, 2012.

% Modified by Romano Cicchetti, ETH Zurich, in the context of the NILM-Eval project

function [bnet, LL, engine] = my_learn_params_dbn_em(engine, evidence, varargin)
% LEARN_PARAMS_DBN Set the parameters in a DBN to their ML/MAP values using batch EM.
% [bnet, LLtrace, engine] = learn_params_dbn_em(engine, data, ...)
%
% data{l}{i,t} = value of node i in slice t of time-series l, or [] if hidden.
%   Suppose you have L time series, each of length T, in an O*T*L array D,
%   where O is the num of observed scalar nodes, and N is the total num nodes per slice.
%   Then you can create data as follows, where onodes is the index of the observable nodes:
%      data = cell(1,L);
%      for l=1:L
%        data{l} = cell(N, T);
%        data{l}(onodes,:) = num2cell(D(:,:,l));
%      end
% Of course it is possible for different sets of nodes to be observed in
% each slice/ sequence, and for each sequence to be a different length.
%
% LLtrace is the learning curve: the vector of log-likelihood scores at each iteration.
%
% Optional arguments [default]
%
% max_iter - specifies the maximum number of iterations [100]
% thresh - specifies the thresold for stopping EM [1e-3]
%   We stop when |f(t) - f(t-1)| / avg < threshold,
%   where avg = (|f(t)| + |f(t-1)|)/2 and f is log lik.
% verbose - display loglik at each iteration [1]
% anneal - 1 means do deterministic annealing (only for entropic priors) [0]
% anneal_rate - geometric cooling rate [0.8]
% init_temp - initial annealing temperature [10]
% final_temp - final annealing temperature [1e-3]
%

max_iter = 100;
thresh = 1e-3;
anneal = 0;
anneal_rate = 0.8;
init_temp = 10;
final_temp = 1e-3;
verbose = 1;

for i=1:2:length(varargin)
  switch varargin{i}
   case 'max_iter', max_iter = varargin{i+1};
   case 'thresh', thresh = varargin{i+1};
   case 'anneal', anneal = varargin{i+1};
   case 'anneal_rate', anneal_rate = varargin{i+1};
   case 'init_temp', init_temp = varargin{i+1};
   case 'final_temp', final_temp = varargin{i+1};
   case 'verbose', verbose = varargin{i+1};
   otherwise, error(['unrecognized argument' varargin{i}])
  end
end

% take 1 EM step at each temperature value, then when temp=0, run to convergence
% When using an entropic prior, Z = 1-T, so
% T=2 => Z=-1 (max entropy)
% T=1 => Z=0 (max likelihood)
% T=0 => Z=1 (min entropy / max structure)
num_iter = 1;
LL = [];
if anneal
  temperature = init_temp;
  while temperature > final_temp
    [engine, loglik, logpost] = EM_step(engine, evidence, temperature);
    if verbose
      fprintf('EM iteration %d, loglik = %8.4f, logpost = %8.4f, temp=%8.4f\n', ...
	      num_iter, loglik, logpost, temperature);
    end
    num_iter = num_iter + 1;
    LL = [LL loglik];
    temperature = temperature * anneal_rate;
  end
  temperature = 0;
  previous_loglik = loglik;
  previous_logpost = logpost;
else
  temperature = 0;
  previous_loglik = -inf;
  previous_logpost = -inf;
end

converged = 0;
while ~converged & (num_iter <= max_iter)
  [engine, loglik, logpost] = EM_step(engine, evidence, temperature);
  if verbose
    %fprintf('EM iteration %d, loglik = %8.4f, logpost = %8.4f\n', ...
    %	    num_iter, loglik, logpost);
    fprintf('EM iteration %d, loglik = %8.4f\n', num_iter, loglik);
  end
  num_iter = num_iter + 1;
  [converged, decreased] = em_converged(loglik, previous_loglik, thresh);
  %[converged, decreased] = em_converged(logpost, previous_logpost, thresh);
  previous_loglik = loglik;
  previous_logpost = logpost;
  LL = [LL loglik];
end

bnet = bnet_from_engine(engine);

%%%%%%%%%

function [engine, loglik, logpost] = EM_step(engine, cases, temp)

bnet = bnet_from_engine(engine); % engine contains the old params that are used for the E step
ss = length(bnet.intra);
CPDs = bnet.CPD; % these are the new params that get maximized
num_CPDs = length(CPDs);

% log P(theta|D) = (log P(D|theta) + log P(theta)) - log(P(D))
% where log P(D|theta) = sum_cases log P(case|theta)
% and log P(theta) = sum_CPDs log P(CPD) - only count once even if tied!
% logpost = log P(theta,D) (un-normalized)
% This should be negative, and increase at every step.

adjustable = zeros(1,num_CPDs);
logprior = zeros(1, num_CPDs);
for e=1:num_CPDs
  adjustable(e) = adjustable_CPD(CPDs{e});
end
adj = find(adjustable);

for e=adj(:)'
  logprior(e) = log_prior(CPDs{e});
  CPDs{e} = reset_ess(CPDs{e});
end

loglik = 0;
for l=1:length(cases)
  evidence = cases{l};
  if ~iscell(evidence)
    error('training data must be a cell array of cell arrays')
  end
  [engine, ll] = enter_evidence(engine, evidence);
  assert(~isnan(ll))
  loglik = loglik + ll;
  T = size(evidence, 2);

  % We unroll ns etc because in update_ess, we refer to nodes by their unrolled number
  % so that they extract evidence from the right place.
  % (The CPD should really store its own version of ns and cnodes...)
  ns = repmat(bnet.node_sizes_slice(:), [1 T]);
  cnodes = unroll_set(bnet.cnodes_slice, ss, T);

  %hidden_bitv = repmat(bnet.hidden_bitv(1:ss), [1 T]);
  hidden_bitv = zeros(ss, T);
  hidden_bitv(isemptycell(evidence))=1;
  % hidden_bitv(i) = 1 means node i is hidden.
  % We pass this in, rather than using isemptycell(evidence(dom)), because
  % isemptycell is very slow.

  t = 1;
  for i=1:ss
    e = bnet.equiv_class(i,1);
    if adjustable(e)
      fmarg = marginal_family(engine, i, t);
      CPDs{e} = update_ess(CPDs{e}, fmarg, evidence, ns(:), cnodes(:), hidden_bitv(:));
    end
  end

  for i=1:ss
    e = bnet.equiv_class(i,2);
    if adjustable(e)
      for t=2:T
	fmarg = marginal_family(engine, i, t);
	CPDs{e} = update_ess(CPDs{e}, fmarg, evidence, ns(:), cnodes(:), hidden_bitv(:));
      end
    end
  end
end

logpost = loglik + sum(logprior(:));

for e=adj(:)'
  CPDs{e} = maximize_params(CPDs{e}, temp);
end

engine = update_engine(engine, CPDs);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function [y] = precondDiag(r,D)
y = D.*r;

author: zouchuhang

repo: LayoutNet


text:
function AddReorderButtons(fighand)
% ADDREORDERBUTTONS  Add buttons to the toolbar to reorder plot objects
% Usage:
%    AddReorderButtons
%    AddReorderButtons(fighand)
% ADDREORDERBUTTONS adds two buttons to the top toolbar of the figure whose
% handle is fighand. One button, if clicked while a plot object is
% selected, moves that object forward in the plot ordering; the other moves
% the object backwards. This is useful primarily for figures with
% 2-dimensional axes containing several plots. If called with no arguments,
% ADDREORDERBUTTONS works on the current figure.
%
% Example:
%    fh = figure;
%    AddReorderButtons(fh);
%
% ADDREORDERBUTTONS may be most useful as part of a callback for the root
% object's DefaultFigureCreateFcn property.
%
% Geoffrey K. Adams, November 2007

if nargin == 0
    fighand = gcf;
end
% Get a handle for the figure toolbar:
tb = findall(fighand,'Type','uitoolbar');
if ~isempty(tb)
    % Add pushbuttons to the top toolbar:
    uipushtool(tb(1), ...
               'Tag', 'MoveForwardTool', ...
               'HandleVisibility', 'off', ...
               'Separator', 'on', ...
               'TooltipString', 'Move Current Plot Object Forward', ...
               'CData', forward_cdata, ...
               'ClickedCallback', @MoveObjForward);
    uipushtool(tb(1), ...
               'Tag', 'MoveBackwardTool', ...
               'HandleVisibility', 'off', ...
               'Separator', 'off', ...
               'TooltipString', 'Move Current Plot Object Backward', ...
               'CData', backward_cdata, ...
               'ClickedCallback', @MoveObjBackward);
end

function MoveObjForward(src,event)
% Move a plot object forward in the ordering, so that it lies on top of
% other plot objects.
obj = gco; % The object to move is the currently selected object.
parent = get(obj,'Parent');
if parent==0 % This means the object is probably the figure; we can't
    return;  % reorder figures, so just ignore.
end
siblings = get(parent,'Children'); % Get all children of the object's
objind = find(obj==siblings);      % parent, and find out this object's
if objind==1                       % index number. If it is already first,
    return;                        % it can't be moved forward any more.
end
% Swap the current object with the next one forward in the list:
siblings([objind-1, objind]) = siblings([objind, objind-1]);
set(parent,'Children',siblings);

function MoveObjBackward(src,event)
% Move a plot object backward in the ordering, so that it lies underneath
% other plot objects.
obj = gco;
parent = get(obj,'Parent');
if parent==0
    return;
end
siblings = get(parent,'Children');
objind = find(obj==siblings);
if objind==numel(siblings); % Just as before, but if the object is already
    return;                 % last, then it can't be moved backward.
end
% Swap the current object with the next one behind in the list:
siblings([objind, objind+1]) = siblings([objind+1, objind]);
set(parent,'Children',siblings);

% The remaining two functions just generate the icons for the buttons:
function c = forward_cdata
c(:,:,1) = [ ...
    NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0   0   0   0   0 NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0 .75 .75 .75 .75 .75   0 NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0 .75 .75 .75 .75 .75   0 NaN; ...
    NaN NaN   0 NaN NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN; ...
    NaN   0 NaN NaN NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN   0   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
      0 NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN   0 NaN NaN   0 NaN   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ];
c(:,:,2) = [ ...
    NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .25 .25 .25 .25 .25   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .25 .25 .25 .25 .25   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .25 .25 .25 .25 .25   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .25 .25 .25   0   0   0   0   0   0   0 NaN; ...
    NaN NaN NaN NaN   0 .25 .25 .25   0   0   0   0   0   0   0 NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0   0   0   0   0 NaN; ...
    NaN NaN   0 NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
    NaN   0 NaN NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN   0   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
      0 NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN   0 NaN NaN   0 NaN   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ];
c(:,:,3) = [ ...
    NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .75 .75 .75   0   0   0   0   0   0   0 NaN; ...
    NaN NaN NaN NaN   0 .75 .75 .75   0   0   0   0   0   0   0 NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0   0   0   0   0 NaN; ...
    NaN NaN   0 NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
    NaN   0 NaN NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN   0   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
      0 NaN NaN NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN   0 NaN NaN   0 NaN   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ];

function c = backward_cdata
c(:,:,1) = [ ...
    NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0   0   0   0   0 NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0 .75 .75 .75 .75 .75   0 NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0 .75 .75 .75 .75 .75   0 NaN; ...
      0   0   0 NaN NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN; ...
    NaN   0   0 NaN NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN; ...
      0 NaN   0 NaN NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN   0 NaN NaN NaN NaN   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN   0   0   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ];
c(:,:,2) = [ ...
    NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .25 .25 .25 .25 .25   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .25 .25 .25 .25 .25   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .25 .25 .25 .25 .25   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .25 .25 .25   0   0   0   0   0   0   0 NaN; ...
    NaN NaN NaN NaN   0 .25 .25 .25   0   0   0   0   0   0   0 NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0   0   0   0   0 NaN; ...
      0   0   0 NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
    NaN   0   0 NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN   0 NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN   0 NaN NaN NaN NaN   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN   0   0   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ];
c(:,:,3) = [ ...
    NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .75 .75 .75 .75 .75   0 NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN   0 .75 .75 .75   0   0   0   0   0   0   0 NaN; ...
    NaN NaN NaN NaN   0 .75 .75 .75   0   0   0   0   0   0   0 NaN; ...
    NaN NaN NaN NaN   0   0   0   0   0   0   0   0   0   0   0 NaN; ...
      0   0   0 NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
    NaN   0   0 NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN   0 NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN   0   0   0   0   0   0   0 NaN; ...
      0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN   0 NaN NaN NaN NaN   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN   0   0   0   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ...
    NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; ];


text:
function [ x, seed ] = r8_normal_01 ( seed )

%*****************************************************************************80
%
%% R8_NORMAL_01 returns a unit pseudonormal R8.
%
%  Discussion:
%
%    The standard normal probability distribution function (PDF) has
%    mean 0 and standard deviation 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the standard normal PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ r1, seed ] = r8_uniform_01 ( seed );
  [ r2, seed ] = r8_uniform_01 ( seed );
  x = sqrt ( -2.0 * log ( r1 ) ) * cos ( 2.0 * pi * r2 );

  return
end


text:
%% process CW vectors

% Input file
outputname = '../data/rt2k/allSStrRT2k.mat';

if ~exist(outputname, 'file')
allSStrRT2k = {};
path = '../data/rt2k/txt_sentoken/pos/*.txt';
names = dir(path);
for i = 1:length(names)
    fid = fopen(['../data/rt2k/txt_sentoken/pos/' names(i).name], 'r');
    content = textscan(fid, '%s', 'delimiter', ' ', 'bufsize', 100000);
    fclose(fid);

    allSStrRT2k{end+1} =  content{1};
end

path = '../data/rt2k/txt_sentoken/neg/*.txt';
names = dir(path);
for i = 1:length(names)
    fid = fopen(['../data/rt2k/txt_sentoken/neg/' names(i).name], 'r');
    content = textscan(fid, '%s', 'delimiter', ' ', 'bufsize', 100000);
    fclose(fid);

    allSStrRT2k{end+1} =  content{1};
end

labels = false(1,length(allSStrRT2k));
labels(1:1000) = true;
save(outputname, 'allSStrRT2k', 'labels');

else
    load(outputname);
end

%%
allSStr = allSStrRT2k;
N = length(allSStr);
themap = containers.Map;

count = 1;
allSNumBi = {};
wordsbi = {};
for i = 1:N
    if mod(i,200) == 0
        fprintf('N %d %d\n', i, count);
    end
    allSNumBi{i} = [];
    for j = 1:length(allSStr{i})

        currentword = lower(allSStr{i}{j});
        try themap(currentword);
        catch
            themap(currentword) = count;
            wordsbi{count} = currentword;
            count = count + 1;
        end
        allSNumBi{i} = [allSNumBi{i} themap(currentword)];
    end
end

for i = 1:N
    if mod(i,1000) == 0
        fprintf('N %d %d\n', i, count);
    end
    for j = 0:length(allSStr{i})
        if j == 0
            cbigram = strcat('*start*', '-', allSStr{i}{j+1});
        elseif j==length(allSStr{i})
            cbigram = strcat(allSStr{i}{j}, '-', '*end*');
        else
            cbigram = strcat(allSStr{i}{j} ,'-', allSStr{i}{j+1});
        end
        cbigram = lower(cbigram);
        try themap(cbigram);
        catch
            themap(cbigram) = count;
            wordsbi{count} = cbigram;
            count = count + 1;
        end

        allSNumBi{i} = [allSNumBi{i} themap(cbigram)];
    end

end

save('../data/rt2k/rt2kbi.mat', 'allSNumBi', 'wordsbi', 'allSStr', 'labels');

author: sidaw

repo: nbsvm

save_path: github-repos/MATLAB/sidaw-nbsvm


text:
function [n, ts_left, ts_right] = nex_int(filename, varname)
% nex_int(filename, varname): Read interval variable from a .nex file
%
% [n, ts_left, ts_right] = nex_int(filename, varname)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   varname - variable name
% OUTPUT:
%   n - number of intervals
%   ts_left - array of left ends of the intervals (in seconds)
%   ts_right - array of right ends of the intervals (in seconds)

n = 0;
ts_left = 0;
ts_right = 0;

if(nargin ~= 2)
   disp('2 input arguments are required')
   return
end

if(ischar(filename) == 0)
   disp('input arguments should be character arrays')
   return
end

if(ischar(varname) == 0)
   disp('input arguments should be character arrays')
   return
end

if(isempty(filename))
   [fname, pathname] = uigetfile('*.nex', 'Select a Nex file');
	filename = strcat(pathname, fname);
end

fid = fopen(filename, 'r');
if(fid == -1)
	disp('cannot open file');
   return
end

disp(strcat('file = ', filename));
magic = fread(fid, 1, 'int32');
version = fread(fid, 1, 'int32');
comment = fread(fid, 256, 'char');
freq = fread(fid, 1, 'double');
tbeg = fread(fid, 1, 'int32');
tend = fread(fid, 1, 'int32');
nvar = fread(fid, 1, 'int32');
fseek(fid, 260, 'cof');
name = zeros(1, 64);
found = 0;
for i=1:nvar
	type = fread(fid, 1, 'int32');
	var_version = fread(fid, 1, 'int32');
	name = fread(fid, [1 64], 'char');
	offset = fread(fid, 1, 'int32');
	n = fread(fid, 1, 'int32');
	name = char(name);
	name = deblank(name);
	k = strcmp(name, deblank(varname));
	if(k == 1)
		if type ~= 2
			disp(sprintf('%s is not an interval variable', deblank(varname)));
			return;
		end
		found = 1;
		fseek(fid, offset, 'bof');
		ts_left = fread(fid, [1 n], 'int32');
		ts_right = fread(fid, [1 n], 'int32');
		break
	end
	dummy = fread(fid, 128, 'char');
end

fclose(fid);

if found == 0
	disp('did not find variable in the file');
else
	ts_left = ts_left/freq;
	ts_right = ts_right/freq;
	disp(strcat('number of intervals = ', num2str(n)));
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function a = r8sp_indicator ( m, n, nz_num, row, col )

%*****************************************************************************80
%
%% R8SP_INDICATOR sets up a R8SP indicator matrix.
%
%  Discussion:
%
%    The R8SP storage format stores the row, column and value of each nonzero
%    entry of a sparse matrix.
%
%    It is possible that a pair of indices (I,J) may occur more than
%    once.  Presumably, in this case, the intent is that the actual value
%    of A(I,J) is the sum of all such entries.  This is not a good thing
%    to do, but I seem to have come across this in MATLAB.
%
%    The R8SP format is used by CSPARSE ("sparse triplet"), DLAP/SLAP
%    ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer NZ_NUM, the number of nonzero elements in the matrix.
%
%    Input, integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices
%    of the nonzero elements.
%
%    Output, real A(NZ_NUM), the nonzero elements of the matrix.
%
  fac = 10^( i4_log_10 ( n ) + 1 );

  for k = 1 : nz_num

    i = row(k);
    j = col(k);
    a(k) = fac * i + j;

  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%	SurfBox-MATLAB (c)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%	Yue Lu and Minh N. Do
%%
%%	Department of Electrical and Computer Engineering
%%	Coordinated Science Laboratory
%%	University of Illinois at Urbana-Champaign
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%	NDFBdec.m
%%
%%	First created: 04-18-05
%%	Last modified: 04-13-06
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Y, Recinfo] = NDFBdec(X, InD, OutD, OutF, Level, HGfname, varargin)

%   NDFB decomposition
%
%   X: N-dimensional (N >= 2) input signal
%
%   InD: "S" if the input X is a usual signal in the spatial domain; "F" if
%   we are given the Fourier transform of X.
%   We use this to get rid of the unnecessary (and time-consuming)
%   "ifftn-fftn" operations in the middle steps.
%
%   OutD: "S" if we want the output to be a usual signal in the spatial "F" if
%   we are want it to be in the Fourier domain.
%   We use this to get rid of the unnecessary (and time-consuming)
%   "ifftn-fftn" operations in the middle steps.
%
%   OutF: output data format. "B" for a multidimensional array output; "C"
%   for a cell array output, with each cell containing one subband.
%
%   Level: an N by N matrix. The k-th row vector specifies the
%   subsequent decomposition levels for the k-th hourglass subband. The
%   individual elements of Level must be great than or equal to 0, except
%   the diagonal elements, which should be -1.
%
%   HGfname: filter name for the hourglass filter bank. Currently supported
%   types are
%   - 'ritf': rotational-invariant tight frame  defined in the Fourier
%   domain.
%
%   Optional inputs: fine-tuning the NXDFB transform
%   - 'bo': the order of the checkerboard filters. Default: bo = 12
%
%   - 'msize': size of the mapping kernel. This controls the quality of the
%   hourglass filters. Default = 15;
%
%   - 'beta': beta value of the Kaiser window. This controls the quality of
%   the hourglass filters. Default = 2.5;
%
%   - 'lambda': the parameter lambda used in the hourglass filter design.
%   Default = 4;
%
%
%   Y: N by 1 cell array. The k-th (1 <= k <= N) cell is a multi-dim array,
%   containing the subbands from the k-th hourglass branch.
%
%   Recinfo: {Level, HGfname, bo, msize, beta, lambda, ...}
%
%   See also NDFBrec.m
%

%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Preparation          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Assign default values
bo = 12;
msize = 15;
beta = 2.5;
lambda = 4;

%% Read optional input arguments
k = 1; nvarargin = nargin - 6;
while k <= nvarargin
    switch varargin{k}
        case 'bo'
            if k == nvarargin
                error('Not enough input parameters.');
            else
                bo = varargin{k+1};
                k = k + 2;
            end
        case 'msize'
            if k == nvarargin
                error('Not enough input parameters.');
            else
                msize = varargin{k+1};
                k = k + 2;
            end
        case 'beta'
            if k == nvarargin
                error('Not enough input parameters.');
            else
                beta = varargin{k+1};
                k = k + 2;
            end
        case 'lambda'
            if k == nvarargin
                error('Not enough input parameters.');
            else
                lambda = varargin{k+1};
                if rem(lambda, 2) ~= 0
                    error('Lambda should be an even interger.');
                end
                k = k + 2;
            end
        otherwise
            error('Unrecognized input.');
    end
end

InD = upper(InD(1));
OutD = upper(OutD(1));
OutF = upper(OutF(1));

%% check the feasibility of the parameter "msize"
minsize = floor((min(size(X)) + 1) / 2);
if rem(minsize, 2) == 0
    minsize = minsize - 1;
end
if (msize > minsize)
    display(['Warning: msize changed from ' num2str(msize) ' to ' num2str(minsize) ' at the current scale.']);
    msize = minsize;
end

%% Reconstructor needs these parameters
Recinfo = {OutD, OutF, Level, HGfname, bo, msize, beta, lambda};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Signal Processing Part %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% First level: undecimated hourglass filter bank
%% Note: the output is given in the Fourier domain with the conjugate
%% symmetric part removed
Subs = HourGlassDec(X, InD, 'F', HGfname, msize, beta, lambda);

%% Dimension of the signal X
N = ndims(X);
%% Output array
Y = cell(N, 1);

%% get the decomposition checkerboard filters
Hflt = get_cbd_filters_load(bo, 'd');

%% Subsequent levels: IRC decomposition at pairs of dimensions
for k = 1 : N
    subband = Subs{k};

    %% free memory space
    Subs{k} = [];

    %% Iteratively Resampled Checkerboard Filter Bank
    subband = IRCdecF(subband, k, Level(k, :), Hflt);

    %% Convert the output to specified format and domain
    if OutD == 'F'
        if OutF == 'C'
            %% just change the data format from 'Box' to 'Cell' array.
            subband = box2cell(subband, Level(k, :));
        else
            %% the easiest case: do nothing.
        end
    else
        %% First change the data format from 'Box' to 'Cell' array.
        subTmp = box2cell(subband, Level(k, :));
        %% Free memory
        subband = cell(length(subTmp), 1);
        %% Perform IFFTN on each subband
        for n = 1 : length(subband)
            subband{n} = real(ifftn(ccsym(subTmp{n}, k, 'e')));
        end
        clear subTmp;

        if OutF == 'B'
            subband = cell2box(subband, Level(k, :));
        end
    end

    Y{k} = subband;

end

%%	This software is provided "as-is", without any express or implied
%%	warranty. In no event will the authors be held liable for any
%%	damages arising from the use of this software.

author: thomaskuestner