text:
function [fitness,convar,xtxi,models] = testlist(stimLists,GA,HRF,varargin)
% [fitness,convar,xtxi,models] = testlist(stimLists,GA,HRF,[opt args])
%
% stimLists can have multiple lists (columns)
%
% 9/12/01 Tor Wager

clockStart = cputime;			  	% keeps track of starting time

for i = 1:length(varargin)
    if isstr(varargin{i})
        switch varargin{i}
        case 'rest', restMatrix = varargin{i+1};
        case 'svi', svi = varargin{i+1};
        case 'S', S = varargin{i+1};
        end
    end
end

% ----------------------------------------------------------------
% * vars setup
% ----------------------------------------------------------------

trans2switch = GA.trans2switch;
trans2block = GA.trans2block;
maxOrder = GA.maxOrder;
conditions = GA.conditions;
freqConditions = GA.freqConditions;
cbalColinPowerWeights = GA.cbalColinPowerWeights;
contrastweights = GA.contrastweights;
NumStimthresh = GA.NumStimthresh;
maxCbalDevthresh = GA.maxCbalDevthresh;
maxFreqDevthresh = GA.maxFreqDevthresh;
ISI = GA.ISI;
TR = GA.TR;
nonlinthreshold = GA.nonlinthreshold;
scanLength = GA.scanLength;
restevery = GA.restevery;
restlength = GA.restlength;
HPlength = GA.HPlength;
xc = GA.xc;
contrasts = GA.contrasts;

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

% work out whether to use contrasts.
if isempty(contrasts)
   nconds = sum(conditions > 0);
   %contrasts = zeros(nconds,nconds);
    %for i = 1:nconds
    %    contrasts(i,i) = 1;
    %end
    %disp('  ...no contrasts entered; using original model')
    docontrasts = 0;
 else
    %disp(' ...contrasts entered')
    docontrasts = 1;
    %contrasts
end

%if ~(size(contrasts,1) == sum(conditions > 0)),
    %disp('	...fyi: Num contrasts does not equal num conditions')
%end

% ----------------------------------------------------------------
% * 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

% 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
% ----------------------------------------------------------------

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
numsamps = ceil(numStim*ISI/TR);

% ----------------------------------------------------------------
% * get smoothing matrix and autocorrelation matrix
% ----------------------------------------------------------------
if isempty(HPlength),HPlength = 'none';,end
if ~isfield(GA,'LPsmooth'),GA.LPsmooth = 1;,end

dofilter = 1;
if strcmp(HPlength,'none') & ~GA.LPsmooth
    dofilter = 0;
end

%disp(['	...setting smoothing, hrf, and autocorrelation matrices, HPlength = ' num2str(HPlength) ', Smoothing = ' num2str(GA.LPsmooth)])

if ~(exist('S') == 1)
%use spm_filter
if strcmp(HPlength,'none') & GA.LPsmooth									% LP only
	%disp('		...LP only')
	[S,KL] = use_spm_filter(TR,numsamps,'hrf','none',[]);
elseif GA.LPsmooth															% HP and LP
	[S,KL,KH] = use_spm_filter(TR,numsamps,'hrf','specify',HPlength);
	%disp('		...HP and LP')
elseif strcmp(HPlength,'none') & ~GA.LPsmooth								% neither
	%disp('		...no filtering')
	S = eye(numsamps);
	dofilter = 0;
else [S,KL,KH] = use_spm_filter(TR,numsamps,'none','specify',HPlength);		% HP only
	%disp('		HP only')
end

end % if ~exist S

if ~(exist('svi') == 1)

if isempty(xc),Vi = eye(numsamps);,disp('Using white noise autocorrelation')
else Vi = getv('make',xc,numsamps);
end

svi = S * Vi * S';

end % if exist svi

clear KL, clear KH, clear Vi, clear myscannerV

% ----------------------------------------------------------------
% * 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.

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

% add intercept
% ------------------------------
if docontrasts, contrasts(:,end+1) = 0;, end

% ----------------------------------------------------------------
% * do it.
% ----------------------------------------------------------------

for z = 1:size(stimLists,2)

	stimList = stimLists(:,z);

% The Switch Hack: transform list of stimuli into list of switches/no switches for each stimtype (R or L)
      if trans2switch, stimList = transform2switches(stimList);,end
      if trans2block,
	     restMatrix = varargin{1};
	     restlength = GA.restlength;

             restList = restMatrix(:,z);
             stimList = transform2block(stimList,restList,restlength);
      end

      % ===== counterbalancing ============================================
      %if (cbalColinPowerWeights(1) > 0)    removed 3/30 to implement criterion rather than weights.
   		[cBal,dummy,maxDev] = getCounterBal(stimList, maxOrder,conditions,freqConditions);
		 % ...if optimizing this measure, add it to fitness scores.
         if (cbalColinPowerWeights(1) > 0),fitnessMatrix(1,z) = cBal;,end
	  %end

	  % ===== frequency discrepancy ========================================
	  % calculate max discrepancy between actual frequencies and requested frequencies
	  %if size(cbalColinPowerWeights,2) > 3
      	%if (cbalColinPowerWeights(4) > 0)
            for i = 1:size(freqConditions,2)
         		freqMat(i) = sum(stimList == conditions(i)) / size(stimList,1);
      	    end
            freqMat = freqMat ./ sum(freqMat);   % adjust for probe ''0''s
			maxFreqDev = max(abs(freqConditions - freqMat));
			% ...if optimizing this measure, add it to fitness scores.
      		if (cbalColinPowerWeights(4) > 0), dummy = 1 - maxFreqDev;,end
		%end
	  %end

	  % criterion measures
	  % ====================================================================
	    go = 1;
	    if docriterion
	  	[maxNumStim,maxrest] = getMaxInARow(stimList,dojitter);
	  	if maxNumStim > NumStimthresh | maxrest > maxrestthresh | maxDev > maxCbalDevthresh | maxFreqDev > maxFreqDevthresh
		  go = 0;
		  error('Criteria not met.')
          	end
	   end

		% =====power calculation =============================================
    if (cbalColinPowerWeights(2) > 0 | cbalColinPowerWeights(3) > 0) & go   % build predictor set of vectors and convolve with HRF
        model = sampleInSeconds(stimList,ISI);
        model = getPredictors(model,HRF);
        model = resample(model,1,TR*10);
	if dosaturation,model = modelSaturation(model,nonlinthreshold);,end						% saturation (nonlinear responses)
        if size(model,1) > numsamps, model = model(1:numsamps,:);,end

	% add intercept
	% ------------------------------
	model(:,end+1) = 1;

        if dofilter
		try
            model = S * model;                                                  % temporal smoothing and HP filter
        catch
            whos model
            whos S
            error('model smoothing: wRoNg Sizes!')
        end
	end
        xtxitx = pinv(model);                                       		% inv(X'S'SX)*(SX)'; pseudoinv of (S*X)

	if docontrasts,
		try
            con1 =  model * contrasts(1,:)';
        catch
            whos model
            whos contrasts
            error('contrast multiplication: wRoNg Sizes!')
        end
	end
	%figure; subplot(2,1,1)
	%plot(model(:,1))
	%title(['Model ' num2str(z) ' Predictor 1'])
	%subplot(2,1,2)
	%plot(con1)
	%title(['Model ' num2str(z) ' Contrast 1'])
	%drawnow; pause(5)

	xtxi{z} = inv(model'*model);
	models{z} = model;

        if docontrasts											% compute the variance of the contrasts
           try
              fitness(z) = 1./(contrastweights * diag(contrasts*xtxitx*svi*xtxitx'*contrasts'));
	      convar(z) = (contrasts(1,:)*xtxitx*svi*xtxitx'*contrasts(1,:)');
              % it's 1/ because we want to minimize the variance, but we maximize fitnessMatrix value.
           catch
             disp('using contrasts: wrong sizes!!!')
             contrasts
             whos xtxitx
             whos svi
             diag(xtxitx*svi*xtxitx')
             error('Exiting...')
           end

        else
           try
               fitness(z) = 1./(contrastweights * diag(xtxitx*svi*xtxitx')); % variance of chosen regressors.
		myconvar = diag(xtxitx*svi*xtxitx');
	       convar(z) = myconvar(1);
           catch
               disp(['no contrasts version: wrong sizes!!!'])
               whos contrastweights
               whos xtxitx
               whos svi
               diag(xtxitx*svi*xtxitx')
               error('Exiting...')
           end
	end
    end

end % loop through designs

%figure;for i = 1:length(models),subplot(1,length(models),i),imagesc(models{i}),
%	title(['Model ' num2str(i) ' v = ' num2str(convar(i)) ' f = ' num2str(fitness(i))])
%end

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

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function model = setCommonZeroStandardGibbsEnergyOfFormation(model, adjustedMetList)
% Sets all Alberty's cofactor metabolites to have a common thermodynamic baseline.
%
% Sets all the exceptional metabolites to have a common baseline.
% i.e. Standard transformed Gibbs energies of reactants with the baseline adjusted
% for certain paired cofactors e.g. `fad` & `fadh2`, such that the
% difference between the two is the same as in Albertys data but
% the absolute values are consistent with the group contribution data
%
% USAGE:
%
%    model = setCommonZeroStandardGibbsEnergyOfFormation(model, adjustedMetList)
%
% INPUT:
%    model:              Thermodynamic model:
%
%                          * .met(m).dGft0 - standard transformed Gibbs energy of formation(kJ/mol)
%                          * .met(m).dGft0Keq - standard transformed Gibbs energy of formation(kJ/mol)
%                          * .met(m).dGft0Source - origin of data, `Keq` or `groupContFileName.txt`
%                          * .met(m).dGft0GroupCont - group. cont. estimate of standard transformed Gibbs energy of formation(kJ/mol)
%
% OPTIONAL INPUT:
%    adjustedMetList:
% OUTPUT:
%    model:              structure with field:
%
%                          * .met(m).dGft0 - Standard transformed Gibbs energies of reactants
%                            with the baseline adjusted for certain paired
%                            cofactors e.g. `fad` & `fadh2`, such that the
%                            difference between the two is the same as in
%                            Albertys data but the absolute values are
%                            consistent with the group contribution data
%
% .. Author: - Ronan M.T. Fleming
%
% .. Here is the list of cofactors in iAF1260 that have thermodynamic
%    properties backcalculated from Equilibrium values by Alberty
%    The compartments here are important since the same reactant in different
%    compartments may have different properties.

if ~exist('adjustedMetList','var')
%This list of cofactors contains metabolites with own baselines reported
%by Alberty in his 2006 book, plus some metabolites that appear on the
%other side of a reaction from Alberty's set. It is important to keep an
%eye on the latter as they must all be on a common baseline when the
%adjustment is done.
adjustedMetList={'coa','aacoa','accoa','ppcoa','mmcoa_R','mmcoa-R',...
    'succoa','gthrd','gthox','q8h2','q8','fmn','fmnh2','nad',...
    'nadh','nmn','fad','fadh2','nadp','nadph','malcoa'};
end

numChar=1;
[allMetCompartments,uniqueCompartments]=getCompartment(model.mets,numChar);

for p=1:length(uniqueCompartments)

    n=1;
    leftAll{n}=['fad' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['fadh2' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['q8' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['q8h2' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['nad' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['nadh' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['nad' '[' uniqueCompartments{p} ']']; % it is essential that the nad nmn pair are in this position
    rightAll{n}=['nmn' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['nadp' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['nadph' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['gthox' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['gthrd' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['fmn' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['fmnh2' '[' uniqueCompartments{p} ']'];
    n=n+1;
    %here is the list of cofactors in E. coli which are covered by Albertys data
    %see:
    %Standard transformed Gibbs energies of coenzyme A derivatives as
    %functions of pH and ionic strength
    %R.A. Alberty / Biophysical Chemistry 104 (2003) 327�334
    coAset={['coa[' uniqueCompartments{p} ']'],...
            ['aacoa[' uniqueCompartments{p} ']'],...
            ['accoa[' uniqueCompartments{p} ']'],...
            ['mmcoa-R[' uniqueCompartments{p} ']'],...
            ['mmcoa-S[' uniqueCompartments{p} ']'],...
            ['ppcoa[' uniqueCompartments{p} ']'],...
            ['succoa[' uniqueCompartments{p} ']'],...
            ['malcoa[' uniqueCompartments{p} ']']};
    leftAll{n}=[];n=n+1;
    leftAll{n}=['ppcoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['succoa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['ppcoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['mmcoa-S' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['ppcoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['coa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=[];n=n+1;
    leftAll{n}=['succoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['coa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['succoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['mmcoa-R' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['malcoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['coa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['malcoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['accoa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=[];n=n+1;
    leftAll{n}=['dpcoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['coa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['dpcoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['2tpr3dpcoa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['dpcoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['pan4p' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=[];n=n+1;
    leftAll{n}=['sbzcoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['coa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['aacoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['coa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    % Note this reaction that contains three cofactors:
    % KAT1, 3-ketoacyl-CoA thiolase [c] : aacoa + coa --> (2) accoa
    leftAll{n}=['accoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['coa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    leftAll{n}=['accoa' '[' uniqueCompartments{p} ']'];
    rightAll{n}=['aacoa' '[' uniqueCompartments{p} ']'];
    n=n+1;
    %check 'mmcoa-R'
    % methylmalonyl-CoA epimerase
    % equation:	[c] : mmcoa-R <==> mmcoa-S

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [nMet,nRxn]=size(model.S);

    exceptionMetBool=false(nMet,1);
    for m=1:nMet
        abbr=model.mets{m};
        abbr=abbr(1:end-3);
        if any(strcmp(abbr,adjustedMetList))
            exceptionMetBool(m)=1;
        end
    end
    fprintf('\n%s\n%s\n','Setting certain Alberty reactants to have a common thermodynamic',...
        'baseline so that this data can be incorporated with group contribution data.')
    fprintf('%s\n','Before baseline edit of Cofactors:')
    fprintf('%s\n','Metabolite properties:')
    %print out exceptional metabolite data before edit to baseline
    fprintf('%10s\t%20s\t%s\t\t%s\t\t%s\n','abbr','albertyAbbr','dGft0Keq','dGft0gc','dGft0');
    for m=1:nMet
        if exceptionMetBool(m)
            fprintf('%10s\t%20s\t%8.4g\t\t%8.4g\t\t%8.4g\n',model.mets{m},model.met(m).albertyAbbreviation,model.met(m).dGft0Keq,model.met(m).dGft0GroupCont,model.met(m).dGft0);
        end
    end
    fprintf('\n\n')

    fprintf('\n\n');
    fprintf('%s\n','Differences between pairs of cofactors (cytoplasm only):')
    fprintf('%22s%10s%10s%10s%10s%10s%16s%16s\n','left-right','#RxnLeft','#RxnRight','#Both','#LNotR','#RNotL','dGtr0Keq','dGrt0');
    for x=1:length(leftAll)
        abbrL=leftAll{x};
        abbrL=[abbrL(1:end-3) '[c]'];
        abbrR=rightAll{x};
        abbrR=[abbrR(1:end-3) '[c]'];

        nCoFcLRxn=length(find(model.S(strcmp(model.mets,abbrL),:)~=0));
        nCoFcRRxn=length(find(model.S(strcmp(model.mets,abbrR),:)~=0));
        nCoFc2Rxn=length(intersect(find(model.S(strcmp(model.mets,abbrL),:)~=0),find(model.S(strcmp(model.mets,abbrR),:)~=0)));
        nCoLeftOnly=length(setdiff(find(model.S(strcmp(model.mets,abbrL),:)~=0),find(model.S(strcmp(model.mets,abbrR),:)~=0)));
        nCoRightOnly=length(setdiff(find(model.S(strcmp(model.mets,abbrR),:)~=0),find(model.S(strcmp(model.mets,abbrL),:)~=0)));

        %nothing to print out if no match
        if any(strcmp(model.mets,abbrL))  && any(strcmp(model.mets,abbrR))
            fprintf('%22s%10i%10i%10i%10i%10i\t%8.4g\t%8.4g\n',[abbrL '-' abbrR],nCoFcLRxn,nCoFcRRxn,nCoFc2Rxn,nCoLeftOnly,nCoRightOnly,...
                model.met(strcmp(model.mets,abbrL)).dGft0Keq-model.met(strcmp(model.mets,abbrR)).dGft0Keq,...
                model.met(strcmp(model.mets,abbrL)).dGft0-model.met(strcmp(model.mets,abbrR)).dGft0);
        end
    end

%     fprintf('%s\n','Differences between pairs of cofactors:')
%     fprintf('%20s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n','left-right','#RxnLeft','#RxnRight','#Both','#LNotR','#RNotL','dGt0Keq','dGt0gc');
%     xx=0;
%     for x=1:length(leftAll)
%         if isempty(leftAll{x})
%             %         fprintf('\n');
%         else
%             nCoFcLRxn=length(find(model.S(strcmp(model.mets,leftAll{x}),:)~=0));
%             nCoFcRRxn=length(find(model.S(strcmp(model.mets,rightAll{x}),:)~=0));
%             nCoFc2Rxn=length(intersect(find(model.S(strcmp(model.mets,leftAll{x}),:)~=0),find(model.S(strcmp(model.mets,rightAll{x}),:)~=0)));
%             nCoLeftOnly=length(setdiff(find(model.S(strcmp(model.mets,leftAll{x}),:)~=0),find(model.S(strcmp(model.mets,rightAll{x}),:)~=0)));
%             nCoRightOnly=length(setdiff(find(model.S(strcmp(model.mets,rightAll{x}),:)~=0),find(model.S(strcmp(model.mets,leftAll{x}),:)~=0)));
%
%             %nothing to print out if no match
%             if any(strcmp(model.mets,leftAll{x}))  && any(strcmp(model.mets,rightAll{x}))
%                 fprintf('%20s\t\t%i\t\t\t%i\t\t%i\t\t%i\t\t%i\t\t%8.4g\t\t%8.4g\n',[leftAll{x} '-' rightAll{x}],nCoFcLRxn,nCoFcRRxn,nCoFc2Rxn,nCoLeftOnly,nCoRightOnly,...
%                     model.met(strcmp(model.mets,leftAll{x})).dGft0Keq-model.met(strcmp(model.mets,rightAll{x})).dGft0Keq,...
%                     model.met(strcmp(model.mets,leftAll{x})).dGft0GroupCont-model.met(strcmp(model.mets,rightAll{x})).dGft0GroupCont);
%             end
%             a=1;
%         end
%     end

    %edit the baseline of a cofactor that is not always used in a pair
    %with another cofactor
    coaSorted=0;
    x=1;
    while x<=length(leftAll)
         abbrL=leftAll{x};
         abbrR=rightAll{x};

        pause(eps)
        if strcmp('coa[c]',abbrR)
            pause(eps)
        end
        %check that both left and right cofactors are always used in the network
        %together because there is no need to change if they are always used
        %together on different sides of the reaction
        if any(strcmp(abbrL,model.mets)) && any(strcmp(abbrR,model.mets))
            %check if coa is involved since coa is a special case that is not a
            %simple cofactor pair
            if ~(any(strcmp(abbrL,coAset) | any(strcmp(abbrR,coAset))))
                nCoLeftOnly=length(setdiff(find(model.S(strcmp(model.mets,abbrL),:)~=0),find(model.S(strcmp(model.mets,abbrR),:)~=0)));
                nCoRightOnly=length(setdiff(find(model.S(strcmp(model.mets,abbrR),:)~=0),find(model.S(strcmp(model.mets,abbrL),:)~=0)));
                %identify which one in each cofactor pair is synthesised/used in
                %another reaction, besides just being used in reactions with its
                %paired cofactor
                if nCoLeftOnly==0 && nCoRightOnly~=0
                    %calculate the difference between the group contribution
                    %estimate and Albertys value
                    baseLineChange=model.met(strcmp(model.mets,abbrR)).dGft0GroupCont-model.met(strcmp(model.mets,abbrR)).dGft0Keq;
                    for m=1:nMet
                        abbr=model.mets{m};
                        if any(strcmp(model.mets{m},{abbrL,abbrR}))
                            %only change if we are using Albertys data
                            if strcmp(model.met(m).dGft0Source,'Keq')
                                model.met(m).dGft0=model.met(m).dGft0+baseLineChange;
                            end
                        end
                    end
                end

                if nCoLeftOnly~=0 && nCoRightOnly==0
                    %calculate the difference between the group contribution
                    %estimate and Albertys value for the left cofactor
                    baseLineChange=model.met(strcmp(model.mets,abbrL)).dGft0GroupCont-model.met(strcmp(model.mets,abbrL)).dGft0Keq;
                    for m=1:nMet
                        if any(strcmp(model.mets{m},{abbrL,abbrR})) % change baseline for both cofactors
                            %only change if we are using Albertys data
                            if strcmp(model.met(m).dGft0Source,'Keq')
                                model.met(m).dGft0=model.met(m).dGft0+baseLineChange;
                            end
                        end
                    end
                    if strcmp(abbrL,['nad' '[' uniqueCompartments{p} ']']) && strcmp(abbrR,['nadh' '[' uniqueCompartments{p} ']']) && any(strcmp(model.mets,['nmn' '[' uniqueCompartments{p} ']'] ))
                        fprintf('%s\n',['Assuming Albertys value for nmn is set to nad baseline']);
                        for m=1:nMet
                            if any(strcmp(model.mets{m},['nmn' '[' uniqueCompartments{p} ']']))
                                %only change if we are using Albertys data
                                if strcmp(model.met(m).dGft0Source,'Keq')
                                    model.met(m).dGft0=model.met(m).dGft0+baseLineChange;
                                end
                            end
                        end
                        %ensure that nad baseline is not changed twice due to nad also
                        %appearing as a cofactor with nmn
                        %skip the next pair as it is nad[c] & nmn[c]
                        x=x+1;
                    end
                end

            else
                %for the cofactors associated with coa by Alberty, then set coA to
                %have the same baseline as the group contribution data
                if coaSorted==0
                    pause(eps)
                    for s=1:length(coAset)
                        %Changed to make suitable for multiple compartments May 24th 2011
                        baseLineChange=model.met(strcmp(model.mets,['coa[' uniqueCompartments{p} ']'])).dGft0GroupCont-model.met(strcmp(model.mets,['coa[' uniqueCompartments{p} ']'])).dGft0Keq;
                        for m=1:nMet
                            if strcmp(coAset{s},model.mets{m})
                                %only change if we are using Albertys data
                                if strcmp(model.met(m).dGft0Source,'Keq')
                                    model.met(m).dGft0=model.met(m).dGft0+baseLineChange;
                                end
                            end
                        end
                    end
                    coaSorted=1;
                end
            end
        end
        x=x+1;
    end
    %compartment by compartment
    coASorted=0;
end

fprintf('\n\n')
fprintf('%s\n','After baseline edit of Cofactors:')
fprintf('%s\n','Metabolite properties:')
fprintf('%10s\t%30s\t%s\t%s\t%s\n','abbr','albertyAbbr','dGft0Keq',' dGft0gc','   dGft0');
%print out the new changes
for m=1:nMet
    if strcmp('nad[c]',model.mets{m})
        pause(eps)
    end
    if exceptionMetBool(m)
        fprintf('%10s\t%30s\t%8.4g\t%8.4g\t%8.4g\n',model.mets{m},model.met(m).albertyAbbreviation,model.met(m).dGft0Keq,model.met(m).dGft0GroupCont,model.met(m).dGft0);
    end
end

fprintf('\n\n');
fprintf('%s\n','Differences between pairs of cofactors, after baseline edit (cytoplasm only):')
fprintf('%22s%10s%10s%10s%10s%10s%16s%16s\n','left-right','#RxnLeft','#RxnRight','#Both','#LNotR','#RNotL','dGtr0Keq','dGrt0');
for x=1:length(leftAll)
    abbrL=leftAll{x};
    abbrL=[abbrL(1:end-3) '[c]'];
    abbrR=rightAll{x};
    abbrR=[abbrR(1:end-3) '[c]'];

    nCoFcLRxn=length(find(model.S(strcmp(model.mets,abbrL),:)~=0));
    nCoFcRRxn=length(find(model.S(strcmp(model.mets,abbrR),:)~=0));
    nCoFc2Rxn=length(intersect(find(model.S(strcmp(model.mets,abbrL),:)~=0),find(model.S(strcmp(model.mets,abbrR),:)~=0)));
    nCoLeftOnly=length(setdiff(find(model.S(strcmp(model.mets,abbrL),:)~=0),find(model.S(strcmp(model.mets,abbrR),:)~=0)));
    nCoRightOnly=length(setdiff(find(model.S(strcmp(model.mets,abbrR),:)~=0),find(model.S(strcmp(model.mets,abbrL),:)~=0)));

    %nothing to print out if no match
    if any(strcmp(model.mets,abbrL))  && any(strcmp(model.mets,abbrR))
        fprintf('%22s%10i%10i%10i%10i%10i\t%8.4g\t%8.4g\n',[abbrL '-' abbrR],nCoFcLRxn,nCoFcRRxn,nCoFc2Rxn,nCoLeftOnly,nCoRightOnly,...
            model.met(strcmp(model.mets,abbrL)).dGft0Keq-model.met(strcmp(model.mets,abbrR)).dGft0Keq,...
            model.met(strcmp(model.mets,abbrL)).dGft0-model.met(strcmp(model.mets,abbrR)).dGft0);
    end
end

%OLD CODE -may still be usefull for looking at nullspace of internal
%reactions that are also involve cofactors
% if ~exist('exceptions','var')
%     exceptions={'coa','aacoa','accoa','ppcoa','mmcoa_R','mmcoa-R',...
%         'succoa','gthrd','gthox','q8h2','q8','fmn','fmnh2','nad',...
%         'nadh','fad','fadh2','nadp','nadph','malcoa'};
% end
%
%
% [nMet,nRxn]=size(model.S);
%
% if ~isfield(model,'biomassRxnAbbr')
%     if ~exist('biomassRxnAbbr')
%         fprintf('\n%s\n','...checkObjective');
%         objectiveAbbr=checkObjective(model);
%         fprintf('%s\n',['Asumming objective is ' objectiveAbbr]);
%         model.biomassRxnAbbr=objectiveAbbr;
%     else
%         model.biomassRxnAbbr=biomassRxnAbbr;
%     end
% end
%
% %finds the reactions in the model which export/import from the model
% %boundary
% %e.g. Exchange reactions
% %     Demand reactions
% %     Sink reactions
% model=findSExRxnInd(model);
% %OUTPUT
% % model.SIntRxnInd          indices of internal reactions
% % model.SExRxnInd           indices of boundary reactions
% intRxnBool=false(1,nRxn);
% intRxnBool(model.SIntRxnInd)=1;
%

%
% %extract subnetwork
% modelE=model;
% %reactions
% exceptionRxnBool=logical(sign(sum(abs(model.S(exceptionMetBool,:)),1)));
% %only keep internal reactions involving exceptional metabolites
% modelE.S=modelE.S(:,exceptionRxnBool & intRxnBool);
% modelE.rxns=modelE.rxns(exceptionRxnBool & intRxnBool);
% modelE.lb=modelE.lb(exceptionRxnBool & intRxnBool);
% modelE.ub=modelE.ub(exceptionRxnBool & intRxnBool);
% %metabolites
% keepMetBool=false(nMet,1);
% for m=1:nMet
%     if nnz(modelE.S(m,:))~=0
%         keepMetBool(m,1)=1;
%     end
% end
% modelE.S=modelE.S(keepMetBool,:);
% modelE.mets=modelE.mets(keepMetBool);
%
% exceptionMetBool=false(nMet,1);
% for m=1:nMet
%     abbr=model.mets{m};
%     if any(strcmp(abbr(end-2:end),'[c]'))
%         abbr=abbr(1:end-3);
%         if any(strcmp(abbr,exceptions))
%             exceptionMetBool(m)=1;
%         end
%     end
% end

author: opencobra

repo: cobratoolbox


text:
%PointFeature  PointCorner feature object
%
% A superclass for image corner features.
%
% Methods::
% plot         Plot feature position
% distance     Descriptor distance
% ncc          Descriptor similarity
% uv           Return feature coordinate
% display      Display value
% char         Convert value to string
%
% Properties::
% u             horizontal coordinate
% v             vertical coordinate
% strength      feature strength
% descriptor    feature descriptor (vector)
%
% Properties of a vector of PointFeature objects are returned as a vector.
% If F is a vector (Nx1) of PointFeature objects then F.u is a 2xN matrix
% with each column the corresponding point coordinate.
%
% See also ScalePointFeature, SurfPointFeature, SiftPointFeature.

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

classdef PointFeature < handle

    properties %(GetAccess=protected)%, Hidden=true)
        u_           % feature x-coordinates
        v_           % feature y-coordinates
        strength_
        descriptor_
        image_id_
    end % properties

    methods

        function f = PointFeature(u, v, strength)
        %PointFeature.PointFeature Create a point feature object
        %
        % F = PointFeature() is a point feature object with null parameters.
        %
        % F = PointFeature(U, V) is a point feature object with specified
        % coordinates.
        %
        % F = PointFeature(U, V, STRENGTH) as above but with specified strength.

            if nargin == 0
                return;
            end
            if nargin >= 2
                f.u_ = u;
                f.v_ = v;
            end
            if nargin == 3
                f.strength_ = strength;
            end
        end

        function val = uv(features)
            val = [[features.u]; [features.v]];
        end

        function plot(features, varargin)
        %PointFeature.plot Plot feature
        %
        % F.plot() overlay a white square marker at the feature position.
        %
        % F.plot(LS) as above but the optional line style arguments LS are
        % passed to plot.
        %
        % If F is a vector then each element is plotted.

            holdon = ishold;
            hold on

            if nargin == 1
                varargin = {'ws'};
            end

            for i=1:length(features)
                plot(features(i).u_, features(i).v_, varargin{:});
            end

            if ~holdon
                hold off
            end
        end % plot

        function s = distance(f1, f2)
        %PointFeature.distance Distance between feature descriptors
        %
        % D = F.distance(F1) is the distance between feature descriptors, the norm
        % of the Euclidean distance.
        %
        % If F is a vector then D is a vector whose elements are the distance between
        % the corresponding element of F and F1.
            for i=1:length(f2)
                s(i) = norm(f1.descriptor-f2(i).descriptor);
            end
        end

        function s = ncc(f1, f2)
        %PointFeature.ncc Feature descriptor similarity
        %
        % S = F.ncc(F1) is the similarty between feature descriptors which is a
        % scalar in the interval -1 to 1, where 1 is perfect match.
        %
        % If F is a vector then D is a vector whose elements are the distance between
        % the corresponding element of F and F1.
            for i=1:length(f2)
                s(i) = dot(f1.descriptor,f2(i).descriptor);
            end
        end

        function [m,corresp] = match(f1, f2, varargin)
        %PointFeature.match Match point features
        %
        % M = F.match(F2, OPTIONS) is a vector of FeatureMatch objects that
        % describe candidate matches between the two vectors of point
        % features F and F2.
        %
        % [M,C] = F.match(F2, OPTIONS) as above but returns a correspodence
        % matrix where each row contains the indices of corresponding features
        % in F and F2  respectively.
        %
        % Options::
        % 'thresh',T    Match threshold (default 0.05)
        % 'median'      Threshold at the median distance
        % 'top',N       Take top N features
        %
        % See also FeatureMatch.

        % TODO: should use a kd-tree or FLANN

            opt.bilateral = false;
            opt = tb_optparse(opt, varargin);

            if opt.bilateral
                % bilateral matching, best match of f1 in f2 must be best match back in f1
                matches = [];
                d = distance([f1.descriptor], [f2.descriptor]);
                for i=1:numrows(d)
                    % find smallest element in the row
                    [best,ii] = min(d(i,:));
                    [~,jj] = min(d(:,ii));
                    if jj == i
                        % we have a match
                        row = d(i,:);
                        row(ii) = Inf;
                        nextbest = min(row);
                        if nextbest < 0.7*best
                            matches = [matches [i; ii]];
                        end
                    end
                end
            else

                [matches,dist,dist2] = closest([f1.descriptor], [f2.descriptor]);
                matches = [1:length(f1); matches];

                % delete matches where distance of closest match is greater than
                % 0.7 of second closest match
                k = dist > 0.7 * dist2;

                matches(:,k) = [];
            dist(k) = [];
            end

            % dist is a 1xM matrix of distance between the matched features, low is good.

            % matches is a 2xM matrix, one column per match, each column
            % is the index of the matching feature in images 1 and 2

            % sort into increasing distance
            [z,k] = sort(dist, 'ascend');
            matches = matches(:,k);
            dist = dist(:,k);

            m = [];
            cor = [];

            for i=1:numcols(matches)
                k1 = matches(1,i);
                k2 = matches(2,i);
                mm = FeatureMatch(f1(k1), f2(k2), dist(i));
                m = [m mm];
                cor(:,i) = [k1 k2]';
            end

            if nargout > 1
                corresp = cor;
            end
        end

        % methods to provide convenient access to properties of object vectors
        function val = u(f)
            val = [f.u_];
        end

        function val = v(f)
            val = [f.v_];
        end

        function val = p(f)
            val = [[f.u_]; [f.v_]];
        end

        function val = strength(f)
            val = [f.strength_];
        end

        function val = descriptor(f)
            val = [f.descriptor_];
        end

        function fp = pick(f, uv)
            %PointFeature.pick Graphically select a feature
            %
            % V = F.pick() is the id of the feature closest to the point clicked
            % by the user on a plot of the image.
            %
            if nargin < 2
                [u,v] = ginput(1);
                uv = f.uv;
            end
            k = closest([u v]', uv);
            fp = f(k);
        end

        function display(f)
        %PointFeature.display Display value
        %
        % F.display() displays a compact human-readable representation of the feature.
        % If F is a vector then the elements are printed one per line.
        %
        % Notes::
        % - This method is invoked implicitly at the command line when the result
        %   of an expression is a PointFeature object and the command has no trailing
        %   semicolon.
        %
        % See also PointFeature.char.

            disp(' ');
            disp([inputname(1), ' = '])
            disp(' ');
            if length(f) > 20
                fprintf('%d features (listing suppressed)\n  Properties:', length(f));
                for property=fieldnames(f)'
                    fprintf(' %s', property{1}(1:end-1));
                end
                fprintf('\n');
            else
                disp( char(f) );
            end
        end % display()

        function ss = char(features)
        %PointFeature.char Convert to string
        %
        % S = F.char() is a compact string representation of the point feature.
        % If F is a vector then the string has multiple lines, one per element.

            ss = [];
            for i=1:length(features)
                f = features(i);
                % display the coordinate
                s = sprintf('  (%g,%g)', f.u_, f.v_);

                % display the other properties
                for property=fieldnames(f)'
                    prop = property{1}; % convert from cell array
                    switch prop
                    case {'u_', 'v_', 'descriptor_'}
                        continue;
                    otherwise
                        val = getfield(f, prop);
                        if ~isempty(val)
                            s = strcat(s, [sprintf(', %s=', prop(1:end-1)), num2str(val, ' %g')]);
                        end
                    end
                end

                % do the descriptor last
                val = getfield(f, 'descriptor_');
                if ~isempty(val)
                    if length(val) == 1
                        % only list scalars or shortish vectors
                        s = strcat(s, [', descrip=', num2str(val, ' %g')]);
                    elseif length(val) < 4
                        % only list scalars or shortish vectors
                        s = strcat(s, [', descrip=(', num2str(val', ' %g'), ')']);
                    else
                        s = strcat(s, ', descrip= ..');
                    end
                end
                ss = strvcat(ss, s);
            end
        end

    end % methods
end % classdef

author: petercorke

repo: machinevision-toolbox-matlab


text:
function test_shortest_paths

%%
msgid = 'matlab_bgl:test_shortest_paths';

%% all_shortest_paths
Dtrue = [0 1 -3 2 -4; 3 0 -4 1 -1; 7 4 0 5 3; 2 -1 -5 0 -2; 8 5 1 6 0];
load('../graphs/clr-26-1.mat');
D = all_shortest_paths(A);
if any(any(D ~= Dtrue))
    error(msgid, 'all_shortest_paths returned an incorrect distance matrix');
end
D = all_shortest_paths(A,struct('algname','johnson'));
if any(any(D ~= Dtrue))
    error(msgid, 'all_shortest_paths(johnson) returned an incorrect distance matrix');
end
D = all_shortest_paths(A,struct('algname','floyd_warshall'));
if any(any(D ~= Dtrue))
    error(msgid, 'all_shortest_paths(floyd_warshall) returned an incorrect distance matrix');
end
At = A';
D = all_shortest_paths(At,struct('istrans',1));
if any(any(D' ~= Dtrue))
    error(msgid, 'all_shortest_paths(istrans=1) returned an incorrect distance matrix');
end
% test non-reachable vertex
A = sparse([1 1; 0 1]);
Dtrue = [0 1; 5 0];
D = all_shortest_paths(A,struct('inf',5,'algname','johnson'));
if any(any(D ~= Dtrue))
    error(msgid, 'all_shortest_paths(johnson,inf=5) returned an incorrect distance matrix');
end
D = all_shortest_paths(A,struct('inf',5,'algname','floyd_warshall'));
if any(any(D ~= Dtrue))
    error(msgid, 'all_shortest_paths(floyd_warshall,inf=5) returned an incorrect distance matrix');
end
% test edge weighted graph
Dtrue = [0 1 -3 2 -4; 3 0 -4 1 -1; 7 4 0 5 3; 2 -1 -5 0 -2; 8 5 1 6 0];
load('../graphs/clr-26-1.mat');
D = all_shortest_paths(A,struct('edge_weight','matrix'));
if any(any(D ~= Dtrue))
    error(msgid, 'all_shortest_paths(weight=matrix) returned an incorrect distance matrix');
end
v=nonzeros(A');
D = all_shortest_paths(spones(A),struct('edge_weight',v));
if any(any(D ~= Dtrue))
    error(msgid, 'all_shortest_paths(weight=matrix) returned an incorrect distance matrix');
end
try
    bc = all_shortest_paths(A,struct('edge_weight',rand(2,1)));
    error(msgid, 'all_shortest_paths(weight=rand(2,1)) did not report an error');
catch
end

% test predecessor matrix
[D P] = all_shortest_paths(A,struct('algname','floyd_warshall'));
for i=1:size(A,1)
    [d p] = shortest_paths(A,i);
    if any(D(i,:)~=d'), error(msgid,'all_shortest_paths(floyd_warshall) returned incorrect distance'); end
    if any(P(i,:)~=p), error(msgid,'all_shortest_paths(floyd_warshall) returned incorrect predecessor'); end
    % the following command should always work.
    for j=1:size(A,1)
        p=[]; while j~=0, p(end+1)=j; j=P(i,j); end; p=fliplr(p);
    end
end

%% shortest_paths


text:
function [Acc,p] = largest_component(A,sym)
% LARGEST_COMPONENT Return the largest connected component of A
%
% Acc = largest_component(A) returns the largest connected component
% of the graph A.  If A is directed, this returns the largest
% strongly connected component.
%
% Acc = largest_component(A,1) returns the largest connected piece of
% a directed graph where connectivity is undirected.  Algorithmically,
% this takes A, drops the directions, then components the largest component
% and returns just this piece of the original _directed_ network.  So the
% output Acc is directed in this case.
%
% [Acc,p] = largest_component(A,...) also returns a logical vector
% indicating which vertices in A were chosen.
%
% See also SCOMPONENTS
%
% Example:
%   load_gaimc_graph('dfs_example')
%   [Acc p] = largest_component(A); % compute the largest component
%   xy2 = xy(p,:); labels2 = labels(p); % get component metadata
%   % draw original graph
%   subplot(1,2,1); graph_draw(A,xy,'labels',labels); title('Original');
%   % draw component
%   subplot(1,2,2); graph_draw(Acc,xy2,'labels',labels2); title('Component');

% David F. Gleich
% Copyright, Stanford University, 2008-2009

% History
% 2009-04-29: Initial coding

if ~exist('sym','var') || isempty(sym), sym=0; end

if sym
    As = A|A';
    [ci sizes] = scomponents(As);
else
    [ci sizes] = scomponents(A);
end
[csize cind] = max(sizes);
p = ci==cind;
Acc = A(p,p);

author: luanfujun

repo: deep-photo-styletransfer


text:
%
% Get a motor trajectory from the b-spline control points,
% using an adaptive method to choose the number of evaluations
% based on the distance along the trajectory.
%
% Input
%
%  P: [ncpt x 2] control points
%  neval: (optional) number of evaluations
%    otherwise, we choose this adaptively
%
% Output
%  stk: [m x 2] trajectory
%
function stk = get_stk_from_bspline(P,neval)
    nland = size(P,1);
    if ~exist('neval','var')

        % set the number of evaluations adaptively,
        % based on the size of the stroke
        PM = defaultps;
        neval = PM.spline_min_neval;
        s = bspline_gen_s(nland,neval);
        stk = bspline_eval(s,P);
        sumdist = sum_pair_dist(stk);
        neval = max(neval,ceil(sumdist./PM.spline_grain));
        neval = min(neval,PM.spline_max_neval);

    end

    s = bspline_gen_s(nland,neval);
    stk = bspline_eval(s,P);
end

function s = sum_pair_dist(D)
    x1 = D(1:end-1,:);
    x2 = D(2:end,:);
    z = sqrt(sum((x1-x2).^2,2));
    s = sum(z);
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function [elem2face,face,elem2faceSign] = dof3face(elem)
%% DOF3FACE dof sturctrue for faces in 3-D.
%
% [elem2face,face,elem2faceSign] = DOF3FACE(elem) constructs data structure
% associated to faces in 3-D. elem is the connectivity matrix for a 3-D
% triangulation. elem2face is the elementwise pointer from elem to face
% indices. face is the face matrix using the ascend ordering. elem2faceSign
% records the consistency of the local faces with induced orientation and
% the global face with ascend orientation.
%
% DOF3FACE is used for elements using dof associated to faces like CR and
% WG elements in 3-D. For elements with orientation such as RT and BDM
% elements, when we use ascend ordering system, no elem2facSign is needed;
% see sc3doc for details.
%
%  See also dofedge, dof3edge
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

NT = size(elem,1);
totalFace = int32([elem(:,[2 3 4]); elem(:,[1 4 3]); ...
                   elem(:,[1 2 4]); elem(:,[1 3 2])]); % induced ordering
matlabversion = version;
if str2double(matlabversion(end-5:end-2)) > 2012
   [face, tempvar, j] = unique(sort(totalFace,2),'rows','legacy');
else
   [face, tempvar, j] = unique(sort(totalFace,2),'rows');
end
elem2face = uint32(reshape(j,NT,4));
elem2faceSign = int8(reshape(sum(sign(diff(totalFace(:,[1:3,1]),1,2)),2),NT,4));

repo: Aero-matlab


text:
% RES = pointOp(IM, LUT, ORIGIN, INCREMENT, WARNINGS)
%
% Apply a point operation, specified by lookup table LUT, to image IM.
% LUT must be a row or column vector, and is assumed to contain
% (equi-spaced) samples of the function.  ORIGIN specifies the
% abscissa associated with the first sample, and INCREMENT specifies the
% spacing between samples.  Between-sample values are estimated via
% linear interpolation.  If WARNINGS is non-zero, the function prints
% a warning whenever the lookup table is extrapolated.
%
% This function is much faster than MatLab's interp1, and allows
% extrapolation beyond the lookup table domain.  The drawbacks are
% that the lookup table must be equi-spaced, and the interpolation is
% linear.

% Eero Simoncelli, 8/96.

function res = pointOp(im, lut, origin, increment, warnings)

%% NOTE: THIS CODE IS NOT ACTUALLY USED! (MEX FILE IS CALLED INSTEAD)

% fprintf(1,'WARNING: You should compile the MEX code for "pointOp", found in the MEX subdirectory.  It is MUCH faster.\n');

X = origin + increment*[0:size(lut(:),1)-1];
Y = lut(:);

res = reshape(interp1(X, Y, im(:), 'linear'),size(im));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [nx, cx, s] = pne_callback_nose(k, nx, cx, px, s, opt)
%PNE_CALLBACK_NOSE  Callback to handle NOSE events
%   [NX, CX, S] = PNE_CALLBACK_NOSE(K, NX, CX, PX, S, OPT)
%
%   Callback to handle NOSE events, triggered by event function PNE_EVENT_NOSE
%   to indicate the limit (nose point) of the continuation curve has been
%   reached.
%
%   This function sets the msg field of the event when the nose point has
%   been found, raises the S.done flag and sets S.done_msg.
%
%   See PNE_CALLBACK_DEFAULT for details of the input and output arguments.

%   MP-Opt-Model
%   Copyright (c) 2016-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Shrirang Abhyankar, Argonne National Laboratory
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%% skip if initialize, finalize or done
if k <= 0 || s.done
    return;
end

%% handle event
if ~s.rollback || nx.step == 0
    ev = pne_detected_event(s.events, 'NOSE', 1);   %% zero only
    if ~isempty(ev)
        if nx.step == 0
            msg = sprintf('Nose point eliminated by limit induced bifurcation at %d continuation steps, lambda = %.4g.', k, nx.x(end));
        else
            msg = sprintf('Reached limit in %d continuation steps, lambda = %.4g.', k, nx.x(end));
        end

        %% the following conditional is only necessary if we also allow
        %% finding the location of the nose-point without terminating
        if ischar(opt.stop_at) && strcmp(opt.stop_at, 'NOSE')
            s.done = 1;
            s.done_msg = msg;
        end
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function w = WeightTransform_Morgan(W,N,D,Ms)

% WeightTransform_Morgan    Mapping of Subband Adaptive Filters Coefficients to
%                           Full-band Filter
%
% Weight transformation by Morgan and Thi (See Section 4.4.3.1)
%
% Arguments:
% W            Tap-weights for subbands i = 0,1,...,(N/2) in each column
% N            Number of subbands
% Ms           Length of adaptive subfilters
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

w = zeros(Ms*D,1);                    % Setup fullband weight vector

Ms = Ms*2;                            % Zero padding
G = fft(W,Ms);                        % Ms-point FFT on (N/2)+1 adaptive subfilters

% Weight trasformation for 2x-oversampling

if D == N/2                           % For 2x-oversampling

% or subband from i = 1,2,...,(N/2)-1, where i is the subband index

    X = zeros(Ms/2,N/2-1);
	for i = 1:2:(N/2-1)           % For i is an odd number
        x = reshape(G(:,i+1),Ms/4,4);
        x = [x(:,2); x(:,3)];
        X(:,i) = x;
	end
	for i = 2:2:(N/2-1)           % For i is an even number
        x = reshape(G(:,i+1),Ms/4,4);
        x = [x(:,4); x(:,1)];
        X(:,i) = x;
	end
    X = X(:);

    i = 0;   x = reshape(G(:,i+1),Ms/4,4); X = [x(:,1); X];

    i = N/2;
    if mod(i,2)==0                    % For N/2 is even, eg. N = 8, 16 subbanbs
        x = reshape(G(:,i+1),Ms/4,4); X = [X; x(:,4)];
    else                              % For N/2 is odd, eg. N = 10 subbands
        x = reshape(G(:,i+1),Ms/4,4); X = [X; x(:,2)];
    end

    X = [X; 0; conj(flipud(X(2:end)))];
    x = real(ifft(X));                % Ms*D-point ifft
    w = x(1:(Ms/2)*D);

end

% Weight trasformation for critical subsampling

if D == N                             % if critical subsampling

% For subbands from i = 1,2,...,(N/2)-1, where i is the subband index

    X = zeros(Ms,N/2-1);
	for i = 1:(N/2-1)
        x = reshape(G(:,i+1),Ms/2,2);
        x = [x(:,2); x(:,1)];
        X(:,i) = x;
	end
    X = X(:);

    i = 0;   x = reshape(G(:,i+1),Ms/2,2); X = [x(:,1); X];
    i = N/2; x = reshape(G(:,i+1),Ms/2,2); X = [X; x(:,2)];

    X = [X; 0; conj(flipud(X(2:end)))];
    x = real(ifft(X));                % Ms*D-point ifft
    w = x(1:(Ms/2)*D);

end

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function varargout = atanh(varargin)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = CreateBasicOperator('increasing','callback');
        operator.derivative = @(x)((1 + x.^2).^-0.5);
        operator.inflection = [-inf -1 0 1];
        operator.domain = [-1 1];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function tP = minus(tP,xy)
% shift triple points in x/y direction
%
% Syntax
%
%   % shift in x direction
%   tP = gB + [100,0]
%
% Input
%  tP - @triplePointList
%  xy - x and y coordinates of the shift
%
% Output
%  tP - @triplePointList

tP = tP + (-xy);

author: mtex-toolbox

repo: mtex

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


text:
function [output] = F_sigmoid(input_layer)
input = input_layer.a;
output = sigmoid(input);

end

author: singaxiong

repo: SignalGraph


text:
% DEMBASISSAMPLE Do a simple demo of a basis function.

% GPSIM

lengthScale = 1;
locationParams = [2 4 6];
numPoints = 100;
numFuncs = 3;

fontName = 'vera';
randn('seed', 1e6)
rand('seed', 1e6)

minLim = min(locationParams)-2*lengthScale;
maxLim = max(locationParams)+2*lengthScale;

colors = {'r', 'g', 'b'};

figure
hold on
basisFunctions = zeros(numPoints, length(locationParams));
t = linspace(minLim, maxLim, numPoints)';
for i = 1:length(locationParams)
  tcentred = t - locationParams(i);
  basisFunctions(:, i) = exp(-tcentred.*tcentred/((lengthScale^2)))/(sqrt(2*pi)*lengthScale);
  a = plot(t, basisFunctions(:, i), colors{i});
  set(a, 'linewidth', 2);
end
pos = get(gca, 'position');
origpos = pos;
%pos(3) = pos(3)/2;
pos(4) = pos(4)/2;
set(gca, 'position', pos);
%set(gca, 'fontname', fontName);
set(gca, 'fontsize', 20);
printPlot(['demBasisSample', num2str(0)], '../tex/diagrams/', '../html/')

decays = [0.01 0.1 1]
sensitivity = [1 1 1]
mrnaBasisFunctions = zeros(numPoints, length(decays), length(locationParams));
sigma2 = lengthScale*lengthScale;
for i = 1:length(decays)
  figure
  hold on
  for j = 1:length(locationParams);
    decays2=decays(i)*decays(i);
    tcentred = t - locationParams(j);
    mrnaBasisFunctions(:, i, j) = exp(...
        (decays2*sigma2/4 - decays(i)*(t-locationParams(j))) ...
        +lnDiffCumGaussian(2/sqrt(2)*repmat( ...
            (decays(i)*sigma2 + locationParams(j))/lengthScale, ...
                           size(t, 1), size(t, 2)),...
        -2/sqrt(2)*(t - decays(i)*sigma2 - locationParams(j))/lengthScale));
    a = plot(t, mrnaBasisFunctions(:, i, j), colors{j});
    set(a, 'linewidth', 2);
  end
  pos = get(gca, 'position');
  origpos = pos;
  %pos(3) = pos(3)/2;
  pos(4) = pos(4)/2;
  set(gca, 'position', pos);
  %set(gca, 'fontname', fontName);
  set(gca, 'fontsize', 20);
  printPlot(['demBasisSample', num2str(0), '_', num2str(i)], '../tex/diagrams/', '../html/')
end

w = randn(length(locationParams), numFuncs);
for i=1:numFuncs
  figure
  f = basisFunctions*w(:, i);
  a = plot(t, f, 'b');
  set(a, 'linewidth', 2);
  pos = get(gca, 'position');
  origpos = pos;
  pos(4) = pos(4)/2;
  set(gca, 'position', pos);
  disp(w(:, i))
  set(a, 'linewidth', 2);
  %set(gca, 'fontname', fontName);
  set(gca, 'fontsize', 20);
  printPlot(['demBasisSample', num2str(i)], '../tex/diagrams/', '../html/')
end

for i=1:length(decays)
  for j=1:numFuncs
    figure
    f = squeeze(mrnaBasisFunctions(:, i, :))*w(:, j);
    a = plot(t, f, 'b');
    set(a, 'linewidth', 2);
    pos = get(gca, 'position');
    origpos = pos;
    pos(4) = pos(4)/2;
    set(gca, 'position', pos);
    set(a, 'linewidth', 2);
    %set(gca, 'fontname', fontName);
    set(gca, 'fontsize', 20);
    printPlot(['demBasisSample', num2str(j), '_', num2str(i)], '../tex/diagrams/', '../html/')
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [pt,p0,v0,t,idx]=surfinterior(node,face)
%
% [pt,p0,v0,t,idx]=surfinterior(node,face)
%
% identify a point that is enclosed by the (closed) surface
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   node: a list of node coordinates (nn x 3)
%   face: a surface mesh triangle list (ne x 3)
%
% output:
%   pt: the interior point coordinates [x y z]
%   p0: ray origin used to determine the interior point
%   v0: the vector used to determine the interior point
%   t : ray-tracing intersection distances (with signs) from p0. the
%       intersection coordinates can be expressed as p0+t(i)*v0
%   idx: index to the face elements that intersect with the ray, order
%       match that of t
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

pt=[];
len=size(face,1);
for i=1:len
   p0=mean(node(face(i,1:3),:));
   plane=surfplane(node,face(i,:));
   v0=plane(1:3);

   [t,u,v]=raytrace(p0,v0,node,face(:,1:3));

   idx=find(u>=0 & v>=0 & u+v<=1.0 & ~isinf(t));
   [ts, uidx]=unique(sort(t(idx)));
   if(~isempty(ts) && mod(length(ts),2)==0)
       ts=reshape(ts,[2 length(ts)/2]);
       tdiff=ts(2,:)-ts(1,:);
       [maxv,maxi]=max(tdiff);
       pt=p0+v0*(ts(1,maxi)+ts(2,maxi))*0.5;
       idx=idx(uidx);
       t=t(idx);
       break;
   end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function x = cg_sparse ( n, a, b, x )

%*****************************************************************************80
%
%% CG_sparse uses the conjugate gradient method for a sparse storage matrix.
%
%  Discussion:
%
%    The linear system has the form A*x=b, where A is a positive-definite
%    symmetric matrix, stored as a full storage matrix.
%
%    The method is designed to reach the solution to the linear system
%      A * x = b
%    after N computational steps.  However, roundoff may introduce
%    unacceptably large errors for some problems.  In such a case,
%    calling the routine a second time, using the current solution estimate
%    as the new starting guess, should result in improved results.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Frank Beckman,
%    The Solution of Linear Equations by the Conjugate Gradient Method,
%    in Mathematical Methods for Digital Computers,
%    edited by John Ralston, Herbert Wilf,
%    Wiley, 1967,
%    ISBN: 0471706892,
%    LC: QA76.5.R3.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, sparse real A(N,N), the matrix.
%
%    Input, real B(N), the right hand side vector.
%
%    Input/output, real X(N).
%    On input, an estimate for the solution, which may be 0.
%    On output, the approximate solution vector.
%

%
%  Initialize
%    AP = A * x,
%    R  = b - A * x,
%    P  = b - A * x.
%
  ap(1:n,1) = a * x;

  r(1:n,1) = b(1:n,1) - ap(1:n,1);
  p(1:n,1) = b(1:n,1) - ap(1:n,1);
%
%  Do the N steps of the conjugate gradient method.
%
  for it = 1 : n
%
%  Compute the matrix*vector product AP = A*P.
%
    ap(1:n,1) = a * p;
%
%  Compute the dot products
%    PAP = P*AP,
%    PR  = P*R
%  Set
%    ALPHA = PR / PAP.
%
    pap = p(1:n,1)' * ap(1:n,1);
    pr =  p(1:n,1)' * r(1:n,1);

    if ( pap == 0.0 )
      return;
    end

    alpha = pr / pap;
%
%  Set
%    X = X + ALPHA * P
%    R = R - ALPHA * AP.
%
    x(1:n,1) = x(1:n,1) + alpha * p(1:n,1);
    r(1:n,1) = r(1:n,1) - alpha * ap(1:n,1);
%
%  Compute the vector dot product
%    RAP = R*AP
%  Set
%    BETA = - RAP / PAP.
%
    rap = r(1:n,1)' * ap(1:n,1);

    beta = - rap / pap;
%
%  Update the perturbation vector
%    P = R + BETA * P.
%
    p(1:n,1) = r(1:n,1) + beta * p(1:n,1);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [energy3M,entropy3M,sumAvg3M,corr3M,invDiffMom3M,contrast3M,...
    clustShade3M,clustPromin3M,haralCorr3M] = gpuTextureByPatch(scanArray3M, nL, ...
    patchSizeV, offsetsM, flagv, hWait, minIntensity,maxIntensity,separateDirnFlag)
% function [energy3M,entropy3M,sumAvg3M,corr3M,invDiffMom3M,contrast3M,...
%     clustShade3M,clustPromin3M] = gpuTextureByPatch(q, nL, ...
%     patchSizeV, offsetsM, flagv, hWait, minIntensity,maxIntensity,separateDirnFlag)
%
% Patch-wise texture calculation.
%
% APA, 09/09/2015

% Generate flags
if ~exist('flagv','var')
    flagv = ones(1,9);
elseif exist('flagv','var') && isempty(flagv)
    flagv = ones(1,9);
end

% Flag to draw waitbar
waitbarFlag = 0;
if exist('hWait','var') && ishandle(hWait)
    waitbarFlag = 1;
end

% Get indices of non-NaN voxels
calcIndM = ~isnan(scanArray3M);
% calcIndM = q > 0;

% % Grid resolution
slcWindow = 2 * patchSizeV(3) + 1;
rowWindow = 2 * patchSizeV(1) + 1;
colWindow = 2 * patchSizeV(2) + 1;

% Build distance matrices
numColsPad = floor(colWindow/2);
numRowsPad = floor(rowWindow/2);
numSlcsPad = floor(slcWindow/2);

% Get number of voxels per slice
[numRows, numCols, numSlices] = size(scanArray3M);
% [numRows, numCols, numSlices] = size(q);
numVoxels = numRows*numCols;

% Quantize the image

% minQ = min(scanArray3M(:));
% maxQ = max(scanArray3M(:));
% dq = (maxQ-minQ)/nL/4;
% levels = linspace(minQ,maxQ,nL);
% levels = linspace(dq,maxQ-dq,nL);
% q2 = imquantize(scanArray3M, levels);
%levels = multithresh(scanArray3M, nL);
%q3 = imquantize(scanArray3M, levels);

% q = imquantize_cerr(scanArray3M,nL);

% Pad q, so that sliding window works also for the edge voxels
%scanArrayTmp3M = padarray(scanArray3M,[numRowsPad numColsPad
%numSlcsPad],NaN,'both'); % aa commented
%q = padarray(q,[numRowsPad numColsPad numSlcsPad],NaN,'both');

% Quantize the image
if exist('minIntensity','var') && exist('maxIntensity','var')
    q = imquantize_cerr(scanArray3M,nL,minIntensity,maxIntensity);
else
    q = imquantize_cerr(scanArray3M,nL);
end

clear scanArray3M

% Pad q, so that sliding window works also for the edge voxels
%scanArrayTmp3M = padarray(scanArray3M,[numRowsPad numColsPad
%numSlcsPad],NaN,'both'); % aa commented
if exist('padarray.m','file')
    q = padarray(q,[numRowsPad numColsPad numSlcsPad],NaN,'both');
else
    q = padarray_oct(q,[numRowsPad numColsPad numSlcsPad],NaN,'both');
end
%nL = 16;
% if any(isnan(q)) % aa commented
%     nL = nL-1; % aa commented
% end % aa commented
%q = imquantize_cerr(scanArrayTmp3M,nL); % aa commented
%clear scanArrayTmp3M; % aa commented
%qmax = max(q(:));
nanFlag = 0; % the quantized image is always padded with NaN. Hence, always,
             % nanFlag = 1.
if any(isnan(q(:)))
    nanFlag = 1;
    q(isnan(q)) = nL+1;
end
lq = nL + 1;

% qs=sort(unique(q));
% lq=length(qs);
% for k = 1:length(qs)
% 	q(q==qs(k)) = k;
% end

q = uint16(q); % q is the quantized image
%numSlcsWithPadding = size(q,3);

% Create indices for 2D blocks
[m,n,~] = size(q);
m = uint32(m);
n = uint32(n);
colWindow = uint32(colWindow);
rowWindow = uint32(rowWindow);
slcWindow = uint32(slcWindow);

% Index calculation adapted from
% http://stackoverflow.com/questions/25449279/efficient-implementation-of-im2col-and-col2im

%// Start indices for each block
start_ind = reshape(bsxfun(@plus,[1:m-rowWindow+1]',[0:n-colWindow]*m),[],1); %//'

%// Row indices
lin_row = permute(bsxfun(@plus,start_ind,[0:rowWindow-1])',[1 3 2]);  %//'

%// Get linear indices based on row and col indices and get desired output
% imTmpM = A(reshape(bsxfun(@plus,lin_row,[0:ncols-1]*m),nrows*ncols,[]));
indM = reshape(bsxfun(@plus,lin_row,(0:colWindow-1)*m),rowWindow*colWindow,[]);

% Directional offsets
numOffsets = size(offsetsM,1);

% Indices of last level to filter out
nanIndV = false([lq*lq,1]);
if nanFlag
    nanIndV([lq:lq:lq*lq-lq, lq*lq-lq:lq*lq]) = true;
end

% Build levels vector for mu, sig
levRowV = repmat(1:lq,[1 lq]);
levColV = repmat(1:lq,[lq 1]);
levColV = levColV(:)';

% Build list of indices for px and contrast calculation
numElems = nL*nL; % APA 7/13
indCtrstM = false(nL,numElems); % APA 7/13
indPxM = false(nL,numElems); % APA 7/13
indPxPlusYm = false(nL,numElems); % APA 7/13
% for n=0:lq-1 % APA 7/13
for n = 0:nL-1
    % indices for p(x-y), contrast
    indCtrstV = false(lq*lq,1);
    indCtrst1V = 1:lq-n;
    indCtrst2V = 1+n:lq;
    indCtrstTmpV = indCtrst1V + (indCtrst2V-1)*lq;
    indCtrstTmpV = [indCtrstTmpV indCtrst2V + (indCtrst1V-1)*lq];
    indCtrstV(indCtrstTmpV) = 1;
    indCtrstV(nanIndV) = [];
    % indCtrstC{n+1} = indCtrstV; % APA 7/13
    indCtrstM(n+1,:) = indCtrstV;

    % indices for px
    indPxV = false(lq*lq,1);
    indPxV(lq*n+1:lq*(n+1)) = true;
    indPxV(nanIndV) = [];
    % indPxC{n+1} = indPxV; % APA 7/13
    indPxM(n+1,:) = indPxV;

end
% for n=1:2*lq
for n=1:2*nL
    % indices for p(x+y)
    indPxPlusYv = false(lq*lq,1);
    indPxPlusYv(levRowV + levColV == n) = 1;
    indPxPlusYv(nanIndV) = [];
    % indPxPlusYc{n} = indPxPlusYv; % APA 7/13
    indPxPlusYm(n,:) = indPxPlusYv;
end

% Build linear indices column/row-wise for Symmetry
indRowV = zeros(1,lq*lq,'uint16');
for i=uint16(1:lq)
    indRowV((i-1)*lq+1:(i-1)*lq+lq) = i:lq:lq*lq;
end

% Filter out NaN levels
levRowV(nanIndV) = [];
levColV(nanIndV) = [];

% Upload indices to GPU
indPxM = gpuArray(single(indPxM));
indPxPlusYm = gpuArray(single(indPxPlusYm));
indCtrstM = gpuArray(single(indCtrstM));

dim = 1;
if separateDirnFlag
    dim = numOffsets;
end

% Initialize
energy3M = [];
entropy3M = [];
sumAvg3M = [];
corr3M = [];
invDiffMom3M = [];
contrast3M = [];
clustShade3M = [];
clustPromin3M = [];
haralCorr3M = [];
if flagv(1)
    energy3M = zeros([numRows, numCols, numSlices, dim],'single');
end
if flagv(2)
    entropy3M = zeros([numRows, numCols, numSlices, dim],'single');
end
if flagv(3)
    sumAvg3M = zeros([numRows, numCols, numSlices, dim],'single');
end
if flagv(4)
    corr3M = zeros([numRows, numCols, numSlices, dim],'single');
end
if flagv(5)
    invDiffMom3M = zeros([numRows, numCols, numSlices, dim],'single');
end
if flagv(6)
    contrast3M = zeros([numRows, numCols, numSlices, dim],'single');
end
if flagv(7)
    clustShade3M = zeros([numRows, numCols, numSlices, dim],'single');
end
if flagv(8)
    clustPromin3M = zeros([numRows, numCols, numSlices, dim],'single');
end
if flagv(9)
    haralCorr3M = zeros([numRows, numCols, numSlices, dim],'single');
end

tic
% Iterate over slices. compute cooccurance for all patches per slice
parfor (slcNum = 1:numSlices,4)

    disp(['--- Texture Calculation for Slice # ', num2str(slcNum), ' ----'])
    if flagv(1), energyV = zeros(dim,numVoxels,'single'); end
    if flagv(2), entropyV = zeros(dim,numVoxels,'single'); end
    if flagv(3), sumAvgV = zeros(dim,numVoxels,'single'); end
    if flagv(4), corrV = zeros(dim,numVoxels,'single'); end
    if flagv(5), invDiffMomV = zeros(dim,numVoxels,'single'); end
    if flagv(6), contrastV = zeros(dim,numVoxels,'single'); end
    if flagv(7), clustShadeV = zeros(dim,numVoxels,'single'); end
    if flagv(8), clustProminV = zeros(dim,numVoxels,'single'); end
    if flagv(9), haralCorrV = zeros(dim,numVoxels,'single'); end

    calcSlcIndV = calcIndM(:,:,slcNum);
    calcSlcIndV = calcSlcIndV(:);
    numCalcVoxs = sum(calcSlcIndV);
    %indSlcM = indM(:,calcSlcIndV); %moved to offset loop
    % List of Voxel numbers
    voxelNumsV = uint32(0:lq*lq:lq*lq*(numCalcVoxs-1));
    %voxelNumsV = uint32(0:nL*nL:nL*nL*(numCalcVoxs-1));
    for off = 1:numOffsets
        % Get voxels for this slice
        indSlcM = indM(:,calcSlcIndV);
        %cooccurPatchM = zeros(lq*lq*numCalcVoxs,1,'double'); % APA 7/13
        %cooccurPatchM = sparse(nL*nL*numCalcVoxs,1); % APA 7/13
        %cooccurPatchM = zeros(nL*nL*numCalcVoxs,1,'single'); % APA 7/14
        offset = offsetsM(off,:);
        % Choose correct neighbors for the selected offset. i.e. the
        % correct rows from indSlcM
        indNoNeighborV = [];
        if offset(1) == 1
            indNoNeighborV = [indNoNeighborV 1:rowWindow:rowWindow*colWindow];
        elseif offset(1) == -1
            indNoNeighborV = [indNoNeighborV rowWindow:rowWindow:rowWindow*colWindow];
        end
        if offset(2) == 1
            indNoNeighborV = [indNoNeighborV 1:colWindow];
        elseif offset(2) == -1
            indNoNeighborV = [indNoNeighborV rowWindow*colWindow:-1:(rowWindow*colWindow-colWindow)+1];
        end
        indSlcM(indNoNeighborV,:) = [];

        % index for various scalar features
        featDim = min(dim,off);

        slcNumV = slcNum:slcNum+slcWindow-1; % slices within the patch
        C = cell(length(slcNumV),1);
        for iSlc = 1:length(slcNumV)
            slc = slcNumV(iSlc);
            if slc-slcNum+offset(3) >= slcWindow
                continue;
            end
            slc1M = uint16(q(numRowsPad+(1:numRows),numColsPad+(1:numCols),slc));
            slc2M = uint16(q(:,:,slc+offset(3)));
            slc2M = circshift(slc2M,offset(1:2));
            slc2M(numRowsPad+(1:numRows),numColsPad+(1:numCols)) = ...
                slc2M(numRowsPad+(1:numRows),numColsPad+(1:numCols)) + (slc1M-1)*lq;
%             for colNum = 1:numCalcVoxs
%                 cooccurPatchTmpM(:,colNum) = ...
%                     accumarray(slc2M(indSlcM(:,colNum)),1,[lq*lq,1]);
%             end
            slc2M = uint32(slc2M(indSlcM));
            %slcTransp2M = uint32(indRowV(slc2M)); % for symmetry
            %nanInd2M = nanIndV(slc2M);
            slc2M = bsxfun(@plus,slc2M,voxelNumsV);
            %nanIndTransp2M = nanIndV(slcTransp2M);
            %slcTransp2M = bsxfun(@plus,slcTransp2M,voxelNumsV);
            cooccurSlcM = accumarray(slc2M(:) ...
                ,1, [lq*lq*numCalcVoxs,1],[],...
                [],true); % patch-wise cooccurance
            indToAddV = find(cooccurSlcM > 0);
            C{iSlc} = [indToAddV,cooccurSlcM(indToAddV)];
            %cooccurPatchM(indToAddV) = cooccurPatchM(indToAddV) + cooccurSlcM(indToAddV);
        end
        IJV = cell2mat(C);
        cooccurPatchM = sparse(IJV(:,1),IJV(:,1).^0,IJV(:,2),lq*lq*numCalcVoxs,1);
        cooccurPatchM = reshape(cooccurPatchM,lq*lq,numCalcVoxs);
        cooccurPatchM = cooccurPatchM + cooccurPatchM(indRowV,:); % for symmetry
        cooccurPatchM(nanIndV,:) = [];
        %cooccurPatchM = cooccurPatchM(~nanIndV,:);
        cooccurPatchM = bsxfun(@rdivide,cooccurPatchM,sum(cooccurPatchM)+1e-5);
        cooccurPatchM = gpuArray(full(cooccurPatchM));
        % Calculate scalar texture for this offset
        % Angular Second Moment (Energy)
        if flagv(1)
            energyV(featDim,calcSlcIndV) = energyV(featDim,calcSlcIndV) + ...
                gather(sum(cooccurPatchM .* cooccurPatchM)); % cooccurPatchM.^2
        end
        % Entropy
        if flagv(2)
            entropyV(featDim,calcSlcIndV) = entropyV(featDim,calcSlcIndV) - ...
                gather(sum(cooccurPatchM.*log2(cooccurPatchM+1e-10)));
        end
        % Contrast, inverse Difference Moment
%         cooccurPatchSiz = size(cooccurPatchM,2); % APA 7/13
%         px = zeros(nL,cooccurPatchSiz,'single'); %px = []; % APA 7/13
%         pXminusY = zeros(nL,cooccurPatchSiz,'single'); %pXminusY = []; % APA 7/13
%         pXplusY = zeros(nL,cooccurPatchSiz,'single'); %pXplusY = []; % APA 7/13
%         for n=0:nL-1
%             cooccurPatchM(indPxM(n+1,:),:);
%             % px
%             % px(n+1,:) = sum(cooccurPatchM(indPxC{n+1},:),1); % APA 7/13
%             px(n+1,:) = sum(cooccurPatchM(indPxM(n+1,:),:),1); % APA 7/13
%             % p(x-y)
%             % pXminusY(n+1,:) = sum(cooccurPatchM(indCtrstC{n+1},:),1); %
%             % APA 7/13
%             pXminusY(n+1,:) = sum(cooccurPatchM(indCtrstM(n+1,:),:),1); % APA 7/13
%             % Contrast
%             if flagv(6)
%                 %contrastV(featDim,calcSlcIndV) = contrastV(featDim,calcSlcIndV) + ...
%                 %    sum(n^2*cooccurPatchM(indCtrstC{n+1},:)); % APA 7/13
%                 contrastV(featDim,calcSlcIndV) = contrastV(featDim,calcSlcIndV) + ...
%                     sum(n^2*cooccurPatchM(indCtrstM(n+1,:),:)); % APA 7/13
%
%             end
%             % inv diff moment
%             if flagv(5)
%                 %invDiffMomV(featDim,calcSlcIndV) = invDiffMomV(featDim,calcSlcIndV) + ...
%                 %    sum((1/(1+n^2))*cooccurPatchM(indCtrstC{n+1},:)); %
%                 %    APA 7/13
%                 invDiffMomV(featDim,calcSlcIndV) = invDiffMomV(featDim,calcSlcIndV) + ...
%                     sum((1/(1+n^2))*cooccurPatchM(indCtrstM(n+1,:),:)); % APA 7/13
%
%             end
%         end
%
%         for n=1:2*nL
%             % p(x+y)
%             % pXplusY(n,:) = sum(cooccurPatchM(indPxPlusYc{n},:),1); % APA
%             % 7/13
%             pXplusY(n,:) = sum(cooccurPatchM(indPxPlusYm(n,:),:),1);
%             % Sum Average
%             if flagv(3)
%                 sumAvgV(featDim,calcSlcIndV) = sumAvgV(featDim,calcSlcIndV) ...
%                     + n*pXplusY(n,:);
%             end
%         end
        px = indPxM * cooccurPatchM;
        %pXminusY = indCtrstM * cooccurPatchM;
        %pXplusY = indPxPlusYm * cooccurPatchM;
        contrastV(featDim,calcSlcIndV) = contrastV(featDim,calcSlcIndV) + ...
            gather((0:nL-1).^2 * indCtrstM * cooccurPatchM);
        invDiffMomV(featDim,calcSlcIndV) = invDiffMomV(featDim,calcSlcIndV) + ...
            gather(1./(1+(0:nL-1).^2) * indCtrstM * cooccurPatchM);
        sumAvgV(featDim,calcSlcIndV) = sumAvgV(featDim,calcSlcIndV) + ...
            gather((1:2*nL) * indPxPlusYm * cooccurPatchM);

        % weighted pixel average (mu), weighted pixel variance (sig)
        mu = (1:nL) * px;
        sig = bsxfun(@minus,(1:nL)',mu);
        sig = sum(sig .*sig .* px, 1);

        if flagv(4) || flagv(7) || flagv(8)
            levIMinusMu = bsxfun(@minus,levRowV',mu);
            levJMinusMu = bsxfun(@minus,levColV',mu);
        end

        % Correlation
        if flagv(4)
            corrV(featDim,calcSlcIndV) = corrV(featDim,calcSlcIndV) + ...
                gather(sum(levIMinusMu .* levJMinusMu  .* cooccurPatchM, 1) ...
                ./ (sig + 1e-10)); % sig.^2 to match ITK results (ITK bug)
        end

        % Cluster Shade
        if flagv(7)
            clstrV = levIMinusMu + levJMinusMu;
            clustShadeV(featDim,calcSlcIndV) = clustShadeV(featDim,calcSlcIndV) + ...
                gather(sum(clstrV.*clstrV.*clstrV .* cooccurPatchM, 1));
        end
        % Cluster Prominence
        if flagv(8)
            clstrV = levIMinusMu + levJMinusMu;
            clustProminV(featDim,calcSlcIndV) = clustProminV(featDim,calcSlcIndV) + ...
                gather(sum(clstrV.*clstrV.*clstrV.*clstrV .* cooccurPatchM, 1));
        end

        % Haralick Correlation
         if flagv(9)
             % muX = mean(px,1);
             muX = 1/nL;
             % sigX = bsxfun(@minus,px,muX);
             sigX = px - muX;
             sigX = sum(sigX .*sigX, 1)/(nL);

%              % Knuth method for mean and standard deviation (like ITK)
%              muX = px(1,:);
%              muPrevX = muX;
%              sigX = muX*0;
%              for col = 2:size(px,1)
%                  muX = muPrevX + (px(col,:) - muPrevX)/col;
%                  sigX = sigX + (px(col,:)-muX).*(px(col,:)-muPrevX);
%                  muPrevX = muX;
%              end
%              sigX = sigX/nL;

             haralCorrV(featDim,calcSlcIndV) = haralCorrV(featDim,calcSlcIndV) + ...
                 gather((levRowV .* levColV * cooccurPatchM - ...
                 muX .* muX) ./ (sigX + eps));   % (levRowV-1) .* (levColV-1) to match ITK? Bug?

         end
    end

    % Average texture from all directions
    if flagv(1)
        energyV = energyV / numOffsets;
        energy3M(:,:,slcNum,:) = reshape(energyV',[numRows, numCols, 1, dim]);
    end
    if flagv(2)
        entropyV = entropyV / numOffsets;
        entropy3M(:,:,slcNum,:) = reshape(entropyV',[numRows, numCols, 1, dim]);
    end
    if flagv(3)
        sumAvgV = sumAvgV / numOffsets;
        sumAvg3M(:,:,slcNum,:) = reshape(sumAvgV',[numRows, numCols, 1, dim]);
    end
    if flagv(4)
        corrV = corrV / numOffsets;
        corr3M(:,:,slcNum,:) = reshape(corrV',[numRows, numCols, 1, dim]);
    end
    if flagv(5)
        invDiffMomV = invDiffMomV / numOffsets;
        invDiffMom3M(:,:,slcNum,:) = reshape(invDiffMomV',[numRows, numCols, 1, dim]);
    end
    if flagv(6)
        contrastV = contrastV / numOffsets;
        contrast3M(:,:,slcNum,:) = reshape(contrastV',[numRows, numCols, 1, dim]);
    end
    if flagv(7)
        clustShadeV = clustShadeV / numOffsets;
        clustShade3M(:,:,slcNum,:) = reshape(clustShadeV',[numRows, numCols, 1, dim]);
    end
    if flagv(8)
        clustProminV = clustProminV / numOffsets;
        clustPromin3M(:,:,slcNum,:) = reshape(clustProminV',[numRows, numCols, 1, dim]);
    end
    if flagv(9)
        haralCorrV = haralCorrV / numOffsets;
        haralCorr3M(:,:,slcNum,:) = reshape(haralCorrV',[numRows, numCols, 1, dim]);
    end

    if waitbarFlag
        set(hWait, 'Vertices', [[0 0 slcNum/numSlices slcNum/numSlices]' [0 1 1 0]']);
        drawnow;
    end

end
toc

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function M = compute_compressed_sensing_matrix(p,n, type, options)

% compute_compressed_sensing_matrix - compute a CS matrix
%
%   M = compute_compressed_sensing_matrix(p,n, type, options);
%
%   M is a (p,n) matrix.
%   type can be
%       'gaussian': gaussian entries
%       'randnormed': random unit norm columns
%       'bernouilli': random +-1/sqrt(p)
%       'orthoproj': random projector M'*M=Id
%       'fourier': random Fourier rows
%       'sincos': random cosine/sine rows
%
%   Copyright (c) 2008 Gabriel Peyre

if nargin<3
    type = 'gaussian';
end

options.null = 0;
switch lower(type)
    case 'gaussian'
    	M = randn(p,n) / sqrt(p);
    case {'randnormed' 'rand'}
        M = randn(p,n);
        M = M ./ repmat( sqrt(sum(M.^2)), [p 1] );
    case 'bernouilli'
        M = sign(randn(p,n)) / sqrt(p);
    case {'orthoproj' 'randproj'}
        M = randn(p,n);
        [M,R] = qr(M'); M = M(:,1:p)';
    case 'fourier'
        [X,Y] = meshgrid(0:n-1:0:n-1);
        M = exp( 2i*pi/n * X.*Y ) / sqrt(p);
        sel = randperm(n); sel = sel(1:p);
        M = M(sel,:);
    case 'sincos'
        [X,Y] = meshgrid(0:n-1,0:n/2);
        M = cos( 2*pi/n * X.*Y );
        [X,Y] = meshgrid(0:n-1,1:n/2-1);
        M = [M; sin( 2*pi/n * X.*Y )];
        M = M ./ repmat( sqrt(sum(M.^2,2)), [1 n] );
        sel = randperm(n); sel = sel(1:p);
        M = M(sel,:) * sqrt( n/p );
    case 'projoptim'
        niter = getoptions(options, 'niter_projoptim', 40);
        M = randn(p,n);
        for i=1:niter
            % unit normed
            M = M ./ repmat( sqrt(sum(M.^2)), [p 1] );
            % projector
            [U,S,V] = svd(M, 'econ');
            S = diag(diag(S)*0+1);
            M = U*S*V';
        end
    otherwise
        error('Unknown matrice type');
end

return;
%% old code

if nargin<4
    normtype = 'normalize';
end

switch lower(type)
    case 'random'
        M = randn(n,m);
    case 'bumps'
        normtype = 'none';
        % m bumps
        sigma = 0.05;
        x = linspace(0,1,m+1)'; x(end) = []; % bump centers
        M = zeros(n,m);
        t = linspace(0,1,n+1)'; t(end) = [];
        for i=1:m
            u = t-x(i); u(u>1/2) = 1-u(u>1/2);
            M(:,i) = exp( -u.^2 / (2*sigma^2)  );
        end
end

switch lower(normtype)
    case 'normalize'
        % lines should be of unit norm
        d = sqrt( sum(M.^2,1) );
        M = M ./ repmat( d, [size(M,1) 1] );
    case 'orthogonalize'
        % M*M' = Id
        M = orth(M')';
    case 'none'
    otherwise
        error('Unknown normalization');
end

author: gpeyre

repo: matlab-toolboxes


text:
function integralImage = computeIntegralImage(image, norm)
% This function computes the integralImage given an image in input and
% evntually normalize it if the flag norm is true.

% First we pad the image this is tricky because when we compute the area of
% the a patch using the size we don't stick to the corner coordinates.
%
%
%     000000000000000000000000000000
%     0#############################
%     0#           |               #
%     0#   SUM     |               #
%     0#           |               #
%     0#-----------x               #
%     0#                           #
%     0#############################
%
%

image = padarray(image, [1 1], 0, 'pre');
[rows, columns, dimensions] = size(image);
integralImage = zeros(size(image));
image=double(image);

% For each row we compute the cumulative sum
% so for example if our row is : 1 1 1
% the cumulative sum is going to be: 1 2 3
for i = 1:rows
    cumulative_row_sum = cumsum(double(image(i, :)));
    if i == 1
       integralImage(i, :) = cumulative_row_sum;
    else
        for j = 1:columns
           integralImage(i, j) = integralImage(i-1, j) + cumulative_row_sum(j) ;
        end
    end
end

% We normalize if the flag is true (this is done if we want to display the result)
if norm==true
    integralImage = 255*(integralImage - min(integralImage(:))) / (max(integralImage(:)) - min(integralImage(:)));
    integralImage = uint8(integralImage);
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function optodedistance = ft_nirs_optodedistance(datain)
% FT_NIRS_OPTODEDISTANCE computes distances between pairs of optodes.
%
% Use as
%   distance = ft_nirs_optodedistance(indata)
% where indata is nirs data.
%
% See also FT_NIRS_REFERENCECHANNELSUBTRACTION

% You are using the FieldTrip NIRS toolbox developed and maintained by
% Artinis Medical Systems (http://www.artinis.com). For more information
% on FieldTrip, see http://www.fieldtriptoolbox.org
%
% This work is licensed under a Creative Commons Attribution-ShareAlike 4.0
% International License. To view a copy of this license, visit
% http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
% Creative Commons Attribution-ShareAlike 4.0 International License:
% -----------------------------------
% You are free to:
%
%     Share - copy and redistribute the material in any medium or format
%     Adapt - remix, transform, and build upon the material
%     for any purpose, even commercially.
%
%     The licensor cannot revoke these freedoms as long as you follow the
%     license terms.
%
% Under the following terms:
%
%     Attribution - You must give appropriate credit, provide a link to
%                    the license, and indicate if changes were made. You
%                    may do so in any reasonable manner, but not in any way
%                    that suggests the licensor endorses you or your use.
%
%     ShareAlike - If you remix, transform, or build upon the material,
%                   you must distribute your contributions under the same
%                   license as the original.
%
%     No additional restrictions - You may not apply legal terms or
%                                   technological measures that legally
%                                   restrict others from doing anything the
%                                   license permits.
%
% -----------------------------------
%
% This toolbox is not to be used for medical or clinical purposes.
%
% Copyright (c) 2016 by Artinis Medical Systems.
% Contact: askforinfo@artinis.com
%
% Main programmer:
% Marc van Wanrooij, DCN, http://www.neural-code.com
% Jörn M. Horschig, Artinis Medical Systems BV, http://www.artinis.com
% $Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the initial part deals with parsing the input options and data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ensure that the input data is raw NIRS-data, this will also do
% backward-compatibility conversions of old data that for example was
% read from an old *.mat file
datain = ft_checkdata(datain, 'datatype', 'raw', 'senstype', 'nirs');

%% Relevant parameters
label	= datain.label; % transformed channel label, combination Receiver and Transmitter
olabel	= datain.opto.optolabel; % optode label
opos	= datain.opto.optopos; % optode positions

npos	= numel(label);

xf		= opos(:,1);
yf		= opos(:,2);

%% determine distance between Receiver and Transmitter optodes
optodedistance				= NaN(npos,1);

for posIdx=1:npos
  str	= label{posIdx};
  c = textscan(str, '%s%s%s', 'Delimiter', {'-', ' '});

  chanRstr = c{1};
  chanTstr = c{2};

  idxR		= match_str(olabel,chanRstr);
  idxT		= match_str(olabel,chanTstr);

  optodedistance(posIdx)	= sqrt( (xf(idxR)-xf(idxT)).^2+(yf(idxR)-yf(idxT)).^2 ); % Pythagorean theorem
end


text:
function multimodelDisplay(model, spacing)

% MULTIMODELDISPLAY Display parameters of the MULTIMODEL model.
% FORMAT
% DESC displays the parameters of the multi-task learning wrapper
% model and the model type to the console.
% ARG model : the model to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG model : the model to display.
% ARG spacing : how many spaces to indent the display of the model by.
%
% SEEALSO : multimodelCreate, modelDisplay
%
% COPYRIGHT : Neil D. Lawrence, 2007

% MLTOOLS

if nargin > 1
  spacing = repmat(32, 1, spacing);
else
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
fprintf('Multi-model:\n')
for i = 1:length(model.comp)
  fprintf('Component %d:\n', i)
  spacing = length(spacing)+2;
  modelDisplay(model.comp{i}, spacing);
end


text:
% Demonstrate posthoc model comparison for spm_mlm_bayes

clear all
close all

N=100;
d=10;
p=5;

x=randn(N,p);
W=randn(p,d);
W(4:5,:)=0;
e=2*randn(N,d);

y=x*W+e;
verbose=1;

% Input-specific shrinkage
mlm = spm_mlm_bayes (y,x,'input',verbose);

figure
imagesc(mlm.wmean);
colormap gray
colorbar
ylabel('Inputs');
xlabel('Outputs');
title('Bayes Regression Coefficients');

figure
imagesc(mlm.wml);
colorbar
ylabel('Inputs');
xlabel('Outputs');
colormap(gray);
title('ML Regression Coefficients');

disp(' ');
disp('Posthoc questions:');
disp(' ');
disp('Is regression coefficient 1,1 non zero ?');
con=zeros(p,d);
con(1,1)=1;
con_vec=con(:)';
disp('Log Evidence in favour of this hypothesis:');
logbf = spm_mlm_posthoc (mlm,con_vec)

disp('Is regression coefficient 4,6 non zero ?');
con=zeros(p,d);
con(4,6)=1;
con_vec=con(:)';
disp('Log Evidence in favour of this hypothesis:');
logbf = spm_mlm_posthoc (mlm,con_vec)

disp('Are regression coefficients 1,1 2,1 and 3,1 non zero ?');
w=zeros(p,d);
w(1,1)=1;w(2,1)=1;w(3,1)=1;
con_vec1 = spm_mlm_makecon (mlm,w);
disp('Log Evidence in favour of this hypothesis:');
logbf = spm_mlm_posthoc (mlm,con_vec1)

disp('If we define w1 = regression coefficients 1,1 2,1 and 3,1');
disp('and w2 = regression coefficients 1,2 2,2 and 3,2');
disp('Is w1 different to w2 ?');
w=zeros(p,d);
w(1,2)=1;w(2,2)=1;w(3,2)=1;
con_vec2 = spm_mlm_makecon (mlm,w);
con_vec_diff=con_vec2-con_vec1;
disp('Log Evidence in favour of this hypothesis:');
logbf = spm_mlm_posthoc (mlm,con_vec_diff)

author: spm


text:
%INVSIGM Inverse sigmoid map
%
% 	W = W*INVSIGM
% 	B = INVSIGM(ARG)
%
% INPUT
%	ARG   Mapping/Dataset
%
% OUTPUT
%	W     Mapping transforming posterior probabilities into distances.
%
% DESCRIPTION
% The inverse sigmoidal transformation to transform a classifier to a
% mapping, transforming posterior probabilities into distances.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, CLASSC, SIGM

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function w = invsigm(arg)
		if nargin == 0

	% Create an empty mapping:
	w = prmapping('invsigm','combiner');
	w = setname(w,'Inverse Sigmoidal Mapping');

elseif isa(arg,'prmapping')

	% If the mapping requested a SIGM transformation (out_conv=1), it
	% is now removed (out_conv=0):
	if arg.out_conv == 1
		w = set(arg,'out_conv',0);
		if arg.size_out == 2
			w.size_out = 1;
			w.labels = w.labels(1,:);
			data.rot = w.data.rot(:,1);
			data.offset = w.data.offset(1);
			data.lablist_in = w.data.lablist_in;
			w.data = data;
		end
	else
		w_s = setname(prmapping('invsigm','fixed'),'Inverse Sigmoidal Mapping');
		w = arg*w_s;
	end

else
	% The data is really transformed:
	if isdatafile(arg)
		w = addpostproc(arg,invsigm);
	else  % datasets and doubles
		w = log(arg+realmin) - log(1-arg+realmin);
	end
end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [ new_pop ] = selection(old_pop, new_pop)
%   Applies tournament selection to generate a child population
%   We did not vectorize this function as we want to keep this procedure
%   separate from those of the crossover and the mutation. Otherwise
%   implementation of different kind of selection strategy will not be
%   straight forward.

global nreal ;
global nbin ;
global nbits ;
global min_binvar ;
global max_binvar ;
global pcross_bin ;
global min_realvar ;
global max_realvar ;
global pcross_real ;
global eta_c ;

[popsize,~] = size(old_pop);

old_pop =  old_pop(randperm(end),:);
% old_pop =  shuffle(old_pop); % SLOW !!!
for i = 1:4:popsize
    p1i = tournament(old_pop,i,  i+1);
    p2i = tournament(old_pop,i+2,i+3);
    if(nreal > 0)
        [c1, c2] = real_cross(old_pop(p1i,1:nreal), ...
                                old_pop(p2i,1:nreal), ...
                                pcross_real, eta_c, ...
                                min_realvar, max_realvar);
    elseif(nbin > 0)
        [c1, c2] = bin_cross_two_point(old_pop(p1i, 1:sum(nbits)), ...
                                old_pop(p2i, 1:sum(nbits)), ...
                                nbits, pcross_bin);
    end
    new_pop(i,  1:nreal) = c1 ;
    new_pop(i+1,1:nreal) = c2 ;
    p1i = tournament(old_pop,i,  i+1);
    p2i = tournament(old_pop,i+2,i+3);
    if(nreal > 0)
            [c1, c2] = real_cross(old_pop(p1i,1:nreal), ...
                                    old_pop(p2i,1:nreal), ...
                                    pcross_real, eta_c, ...
                                    min_realvar, max_realvar);
    elseif(nbin > 0)
            [c1, c2] = bin_cross_two_point(old_pop(p1i,1:sum(nbits)), ...
                                    old_pop(p2i,1:sum(nbits)), ...
                                    pcross_bin);
    end
    new_pop(i+2,1:nreal) = c1 ;
    new_pop(i+3,1:nreal) = c2 ;
end
end

function [shuffled_pop] = shuffle(pop)
%   This function uses the legacy rng for shuffling,
%   but too slow.
    [popsize,cols] = size(pop);
    shuffled_pop = zeros(popsize,cols);
    for i = 1:popsize
        shuffled_pop(i,:) = pop(rnd(1,popsize),:);
    end
end

author: chudur-budur


text:
function transformed_slice_figure = transformed_sliceBrowser(transformed_slice_figure, save_location, f, highlight_point, relevant_slice, min_dist, ...
                                                clickX, clickY, point_ind, add)

% go through histology at the same time

transformed_images_folder = fullfile(save_location, ['transformations' filesep]);
transformed_images = dir([transformed_images_folder filesep '*.tif']);

try
    ud_transformed_slice.transformed_image_names = natsortfiles({transformed_images.name});
    browse_transformed_slices = true;
catch
    disp('no transformed images found.'); ud_transformed_slice.transformed_image_names = {};
    browse_transformed_slices = false;
end

if browse_transformed_slices

    try
        figure(transformed_slice_figure);
    catch
        transformed_slice_figure = figure('Name','Transformed Slice & Probe Point Viewer');
    end
    clf(transformed_slice_figure)

    ud_transformed_slice.current_plot_handles = [];
    ud_transformed_slice.save_location = save_location;

    ud = get(f, 'UserData');

    ud_transformed_slice.transformed_images = transformed_images;
    ud_transformed_slice.transformed_images_folder = transformed_images_folder;
    ud_transformed_slice.sliceAx = axes('Position', [0.05 0.05 0.9 0.9]);
    hold(ud_transformed_slice.sliceAx, 'on');
    set(ud_transformed_slice.sliceAx, 'HitTest', 'off');
    title('Probe on Slice Viewer');

    ud_transformed_slice.im = plotTVslice(zeros(ud.ref_size, 'uint8'));

    % for probe view mode
    if highlight_point
        ud_transformed_slice.all_slices_slice_num = relevant_slice;
        if min_dist < 10
            x_plot = ud.pointList{ud.currentProbe,1}(point_ind,1);
            y_plot = ud.ref_size(1)-ud.pointList{ud.currentProbe,1}(point_ind,2);
        else
            x_plot = clickX;
            y_plot = ud.ref_size(1) - clickY;
        end
            figure(transformed_slice_figure)
            ud_transformed_slice.quiver_plot{1} = scatter(ud_transformed_slice.sliceAx, x_plot, y_plot, 200, 'ro', ...
            'MarkerFaceColor',[1 1 1],'MarkerEdgeColor', [.3 .3 .3], ...
            'MarkerFaceAlpha',.5, 'MarkerFaceAlpha',.1);

        ud_transformed_slice.quiver_plot{2} = quiver( x_plot - 30 - 10, y_plot + 30 + 10, 30, -30, 1, 'color', 'white', 'linewidth',1,'MaxHeadSize',10);
    else
        ud_transformed_slice.all_slices_slice_num = ud.slice_at_shift_start+ud.slice_shift+add;
        ud_transformed_slice.quiver_plot = [];
    end

    processed_images = dir([ud_transformed_slice.save_location filesep '*.tif']);
    processed_image_names = natsortfiles({processed_images.name});
    ud_transformed_slice.total_num_files = length(processed_image_names);

    ud_transformed_slice.slice_num =  find(strcmp([processed_image_names{ud_transformed_slice.all_slices_slice_num}(1:end-4)  '_transformed.tif'],ud_transformed_slice.transformed_image_names));

    % show image
    if size(ud_transformed_slice.slice_num,2)
        processed_image_name = ud_transformed_slice.transformed_image_names{ud_transformed_slice.slice_num};
        current_slice_image = flip(imread([transformed_images_folder processed_image_name]),1);
        ud_transformed_slice.extra_text = ' (transformed)';
    else
        processed_image_name = processed_image_names{ud_transformed_slice.all_slices_slice_num};
        current_slice_image = flip(imread(fullfile(save_location, processed_image_name)),1);
        ud_transformed_slice.extra_text = ' (not transformed)';
    end

    set(ud_transformed_slice.im, 'CData', current_slice_image);
    title(['Probe on Slice Viewer ' num2str(ud_transformed_slice.all_slices_slice_num) '/' num2str(ud_transformed_slice.total_num_files) ud_transformed_slice.extra_text])

    set(transformed_slice_figure, 'KeyPressFcn', @(slice_figure, keydata)SliceAtlasHotkeyFcn(transformed_slice_figure, keydata, f, transformed_images_folder));
    set(transformed_slice_figure, 'UserData', ud_transformed_slice)
end

% ------------------------
% react to keyboard press
% ------------------------
function SliceAtlasHotkeyFcn(fig, keydata, f, transformed_images_folder)

ud = get(fig, 'UserData');
ud_atlas_viewer = get(f, 'UserData');

processed_images = dir([ud.save_location filesep '*.tif']);
processed_image_names = natsortfiles({processed_images.name});

transformed_images = dir([transformed_images_folder filesep '*.tif']);
try
    ud.transformed_image_names = natsortfiles({transformed_images.name});
catch
    ud.transformed_image_names = {};
end

ud.transformed_images = transformed_images;

if strcmp(keydata.Key,'leftarrow')
    if ud.all_slices_slice_num > 1
        ud.all_slices_slice_num = ud.all_slices_slice_num - 1;

        ud.slice_num =  find(strcmp([processed_image_names{ud.all_slices_slice_num}(1:end-4)  '_transformed.tif'],ud.transformed_image_names));

        % show next image
        if size(ud.slice_num,2)
            processed_image_name = ud.transformed_image_names{ud.slice_num};
            current_slice_image = flip(imread(fullfile(ud.transformed_images_folder, processed_image_name)),1);
            ud.extra_text = ' (transformed)';
        else
            processed_image_name = processed_image_names{ud.all_slices_slice_num};
            current_slice_image = flip(imread(fullfile(ud.save_location, processed_image_name)),1);
            ud.extra_text = ' (not transformed)';
        end
        % reduce to 3 channels at most
        color_channels = min( 3, size(image,3));
        current_slice_image = current_slice_image(:,:,1:color_channels);
        set(ud.im, 'CData', current_slice_image);

    end
elseif strcmp(keydata.Key,'rightarrow') % break
    if ud.all_slices_slice_num < ud.total_num_files
       ud.all_slices_slice_num = ud.all_slices_slice_num + 1;

        ud.slice_num =  find(strcmp([processed_image_names{ud.all_slices_slice_num}(1:end-4)  '_transformed.tif'],ud.transformed_image_names));

        % show next image
        if size(ud.slice_num,2)
            processed_image_name = ud.transformed_image_names{ud.slice_num};
            current_slice_image = flip(imread(fullfile(ud.transformed_images_folder, processed_image_name)),1);
            ud.extra_text = ' (transformed)';
        else
            processed_image_name = processed_image_names{ud.all_slices_slice_num};
            current_slice_image = flip(imread(fullfile(ud.save_location, processed_image_name)),1);
            ud.extra_text = ' (not transformed)';
        end
        % reduce to 3 channels at most
        color_channels = min( 3, size(image,3));
        current_slice_image = current_slice_image(:,:,1:color_channels);
        set(ud.im, 'CData', current_slice_image);
    end
else
% otherwise -- call function to atlas browser
    fcn = get(f, 'KeyPressFcn');
    fcn(f, keydata);
end
        title(['Probe on Slice Viewer -- Slice ' num2str(ud.all_slices_slice_num) '/' num2str(ud.total_num_files) ud.extra_text ])
        % plot probe points for that slice
        set(ud.current_plot_handles(:), 'Visible', 'off'); ud.current_plot_handles = [];

        for probe = 1:size(ud_atlas_viewer.pointList,1)
            for point = 1:size(ud_atlas_viewer.pointList{probe,1},1)
                if ud.all_slices_slice_num == ud_atlas_viewer.pointList{probe,2}(point)
                    ud.current_plot_handles(end+1) = scatter(ud.sliceAx, ud_atlas_viewer.pointList{probe,1}(point,1), ...
                        ud_atlas_viewer.ref_size(1)-ud_atlas_viewer.pointList{probe,1}(point,2), 20, 'ro', 'MarkerFaceColor',[ .1 .1 .1],...
                    'MarkerEdgeColor', ud_atlas_viewer.ProbeColors(probe, :), 'MarkerFaceAlpha',.4,'LineWidth',1.5);
                end
            end
        end

set(fig, 'UserData', ud);

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function [exponents,base]=getexponentbase(p,x)
%GETEXPONENTBASE Internal function used in SOS programs

if isempty(p)
    exponents=[];
    base=[];
else
    p_vars = getvariables(p);
    x_vars = getvariables(x);

    base = getbase(p);

    monom_table = yalmip('monomtable');
    exponents = monom_table(p_vars,x_vars);
    if any(base(:,1))%base(1)~=0
        exponents = [spalloc(1,size(exponents,2),0);exponents];
    else
        base = base(:,2:end);
    end
end