text:
%klas_isoclass 'Performs Isoclass Clustering on multispectral image data '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros las_isoclass.pane file
%
% Parameters:
% InputFile: i 'Input data object', required: 'Input data object'
% OutputFile: o 'Output cluster object', required: 'Resulting output cluster assignment data object'
% Integer: maxnumit 'Maximum number of iterations', default: 16: 'Maximum number of iterations'
% Double: clusdist 'Threshold for combining clusters', default: 3.2: 'Threshold for combining clusters'
% Double: maxclstd 'Threshold for splitting clusters', default: 4.5: 'Threshold for splitting clusters'
% Double: ssepval 'Split separation value for clusters', default: 0: 'Split separation value for clusters'
% Integer: minnum 'Minimum number of members allowed in a cluster', default: 30: 'Minimum number of members allowed in a cluster'
% Integer: maxclust 'Maximum number of clusters', default: 16: 'Maximum number of clusters'
% Double: chnthr 'Threshold for cluster chaining', default: 3.2: 'Threshold for cluster chaining'
% OutputFile: ascii 'Isoclass statistics output (ASCII)', optional: 'Isoclass statistics output (ASCII)'
%
% Example: [o, ascii] = klas_isoclass(i, {'i','';'o','';'maxnumit',16;'clusdist',3.2;'maxclstd',4.5;'ssepval',0;'minnum',30;'maxclust',16;'chnthr',3.2;'ascii',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% las_isoclass - Performs Isoclass Clustering on multispectral image data
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% For a detailed description of the algorithm, see llas_isoclass.
%
%  RESTRICTIONS
%
%  REFERENCES
% This program is modeled after the Land Analysis System (LAS)
% program "isoclass."
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = klas_isoclass(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = klas_isoclass(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output';'maxnumit', 16;'clusdist', 3.2;'maxclstd', 4.5;'ssepval', 0;'minnum', 30;'maxclust', 16;'chnthr', 3.2;'ascii', '__output'};
maxval={0,0,2147483647,1,1,1,2147483647,255,1,1};
minval={0,0,0,1,1,1,0,2,1,1};
istoggle=[0,0,1,1,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Integer','Double','Double','Double','Integer','Integer','Double','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'las_isoclass"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [fdMc, bH, fPval, fKsstat] = calc_KStest_dMc(mCat1, mCat2)
% ----------------------------------------------------------------------------------------------
%
% Incoming variable
% mCat1 : EQ catalog period 1 (Catalog to be modified)
% mCat2 : EQ catalog period 2 (Observed catalog)
%
% Outgoing variable
%
% Author: J. Woessner, woessner@seismo.ifg,.ethz.ch
% last update: 21.01.04

fBinning = 0.1;
nSample = 100;
nMethod = 6;

% Determine exact time period
fPeriod1 = max(mCat1(:,3)) - min(mCat1(:,3));
fPeriod2 = max(mCat2(:,3)) - min(mCat2(:,3));

% Initialize values
fMinMag = min([min(mCat1(:,6)) min(mCat2(:,6))]);
fMaxMag = max([max(mCat1(:,6)) max(mCat2(:,6))]);

%% Calculate model for best fitting Mc, both periods
[mResult1, fMls1, fMc1, fMu1, fSigma1, mDataPredBest1, vPredBest1, fBValue1, fAvalue1] = calc_McCdfnormal(mCat1, fBinning);
[mResult, fMls, fMc2, fMu, fSigma, mDataPredBest, vPredBest, fBValue2, fAvalue2] = calc_McCdfnormal(mCat2, fBinning);

% Using bootstrap approach
% [fMc1, fStd_Mc1, fBvalue1, fStd_B1, fAvalue1, fStd_A1, vMc1, mBvalue1] = calc_McBboot(mCat1, fBinning, nSample, nMethod)
% [fMc2, fStd_Mc2, fBvalue2, fStd_B2, fAvalue2, fStd_A2, vMc2, mBvalue2] = calc_McBboot(mCat2, fBinning, nSample, nMethod)

% Mc difference
fdMc = fMc2-fMc1;

try
    % Select part of catalog
    vSel = mCat1(:,6) >= fMc2-fBinning/2;
    mCat1tmp=mCat1(vSel,:);

    % Create Model distribution
    vMagnitudes = [fMc2:0.1:floor(max(mCat2(:,6)))+0.1];
    % Productuvity for fMc2
    nNumberEvents = length(mCat1tmp(:,1));
    %vNumbers = 10.^(log10(nNumberEvents) - fBValue1*(vMagnitudes-fMc2));
    vNumbers = 10.^(fAvalue1- fBValue1*(vMagnitudes));%-fMc2
    vNumbers = round(vNumbers);
    vNCumFMD = round(-diff(vNumbers));
    % Calculate synthetic data below Mc
    fMinMag = min(mCat2(:,6));
    vMagstep = fMinMag:0.1:fMc2-0.1;
    vProb = normcdf(vMagstep,fMu1, fSigma1);
    vProb = vProb';
    vMagstep = vMagstep';

    % Calculate number of EQs in bins
    fN_Mc = vNCumFMD(1,1);
    vN = round(vProb(:,1)*fN_Mc);
    mNonCumModel = [vN vMagstep; vNCumFMD' vMagnitudes(:,1:end-1)'];

    % vMags = mDataPredBest(:,2)';
    % FMD to be modeled (second period)
    [vObsFMD,vBin2] = hist(mCat2(:,6),roundn(min(mCat2(:,6)),-1):0.1:floor(max(mCat2(:,6))));

    % Select bins to calculate loglikelihood
    vSel = (mNonCumModel(:,2) < min(vBin2) | mNonCumModel(:,2) > max(vBin2));
    vPredFMD = mNonCumModel(~vSel,:);
    % Normalize
    vObsFMD = ceil(vObsFMD./fPeriod2);
    vPredFMD(:,1) = vPredFMD(:,1)./fPeriod1;

    if length(vObsFMD') ~= length(vPredFMD(:,1))
        disp('warning')
    end
    % Calculate the likelihood
    vProb_ = calc_log10poisspdf2(vObsFMD', vPredFMD(:,1));
    % Sum the probabilities
    fProbability = (-1) * sum(vProb_);

    nDegFree = 1; % degree of freedom
    n_samples = length(mCat2(:,6))+length(mCat1(:,6));
    %% Corrected Akaike Information Criterion (AICc)
    fAICc = -2*(-fProbability)+2*nDegFree+2*nDegFree*(nDegFree+1)/(n_samples-nDegFree-1);

catch
    fdMc = nan;
    fProbability = nan;
    fAICc = nan;
end

% % Create model
% fMinMag = min(mCat2(:,6));
% vMagstep = fMinMag:0.1:fMc2-0.1;
% vProb = normcdf(vMagstep,fMu1, fSigma1);
% vProb = vProb';
% vMagstep = vMagstep';
% mDataPredBest1(:,1) = mDataPredBest1(:,1)./fPeriod1;
%
% % Calculate number of EQs in bins
% vSel = (roundn(mDataPredBest1(:,2),-1) == roundn(fMc2,-1));
% fN_Mc = mDataPredBest1(vSel,1);
% vN = round(vProb(:,1)*fN_Mc);
%
% vSel = mDataPredBest1(:,2) >= fMc2;
% mModelFMD = [vN vMagstep; mDataPredBest1(vSel,1) mDataPredBest1(vSel,2)];

% Reconstruct vector of magnitudes from model for Period 1
vMag = [];
% % Normalize to time period
% mModelFMD(:,1) = round(mModelFMD(:,1));
mNonCumModel(:,1) = round(mNonCumModel(:,1)./fPeriod1);
vSel = (mNonCumModel(:,1) ~= 0); % Remove bins with zero frequency of zero events
mData = mNonCumModel(vSel,:);
% mModelFMD = [round(mDataPredBest1(:,1)./fPeriod1) mDataPredBest1(:,2)];
% vSel = (mModelFMD(:,1) ~= 0); % Remove bins with zero frequency of zero events
% mData = mModelFMD(vSel,:);
for nCnt=1:length(mData(:,1))
    fM = repmat(mData(nCnt,2),mData(nCnt,1),1);
    vMag = [vMag; fM];
end

% Reconstruct vector of magnitudes from model for Period 2
% from normalized FMD
vMag2 = [];
%mObsFMD = [vObsFMD' vBin2'];
mObsFMD = [round(mDataPredBest(:,1)./fPeriod2) mDataPredBest(:,2)];
vSel = (mObsFMD(:,1) ~= 0); % Remove bins with zero frequency of zero events
mData2 = mObsFMD(vSel,:);
for nCnt=1:length(mData2(:,1))
    fM = repmat(mData2(nCnt,2),mData2(nCnt,1),1);
    vMag2 = [vMag2; fM];
end

[bH,fPval,fKsstat] = kstest2(roundn(vMag2,-1),roundn(vMag,-1),0.05,0);

figure
h1=cdfplot(vMag2)
set(h1,'Color',[1 0 0]);
hold on
h2=cdfplot(vMag)
set(h2,'Color',[0 0 0]);

h3=cdfplot(mCat2(:,6))
set(h2,'Color',[0 1 0]);
h4=cdfplot(mCat1(:,6))
set(h4,'Color',[0 0 1]);

figure
histogram(vMag2,min(vMag2):0.1:max(vMag2))
figure
histogram(vMag,min(vMag):0.1:max(vMag))

author: CelsoReyes


text:
function plot3d (type)
global A cds opt holdon leg label figmain smmethod subpl

% Check if there IS data to plot
if isempty(A) errordlg ('No data present. Sure you read a file yet?','Error'); return; end

% Add or substitute plot
figure (figmain);
subplot(subpl(1),subpl(2),subpl(3),'Parent',figmain);
if holdon hold on; else 	hold off; end

% Columns
X= A(:,opt.xc(1));
Y= A(:,opt.yc(1));
M= A(:,opt.zc);

switch type
	case 'histo'
		hist(Y,10);
	case 'histo3'
		hist3([X,Y],[10,10]);
	case 'corrxy'
		[rho,pval]= corr(X,Y);
		plot(X,Y,'o');
		L= ['rho= ',num2str(rho),' p< ',num2str(pval)];
		legend(L);
	case 'corrall'
		fprintf('Col Row  Rho  Pval \n')
		output= mcorr(A,'sig');
		fprintf('  %d   %d %4.2f  %4.2f \n',output)
	case 'qqplot'
		qqplot(X,Y);
	case 'normplot'
		normplot(X);

end

author: Sable

repo: mcbench-benchmarks


text:
function sparse_grid_test013 ( dim_min, dim_max, level_max_min, ...
  level_max_max )

%*****************************************************************************80
%
%% TEST013 tests SPARSE_GRID_F2S_SIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_MIN, the minimum spatial dimension to consider.
%
%    Input, integer DIM_MAX, the maximum spatial dimension to consider.
%
%    Input, integer LEVEL_MAX_MIN, the minimum value of LEVEL_MAX to consider.
%
%    Input, integer LEVEL_MAX_MAX, the maximum value of LEVEL_MAX to consider.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST013\n' );
  fprintf ( 1, '  SPARSE_GRID_F2S_SIZE returns the number of\n' );
  fprintf ( 1, '  distinct points in a sparse grid made up of all \n' );
  fprintf ( 1, '  product grids formed from Fejer Type 2 Slow \n' );
  fprintf ( 1, '  quadrature rules.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The sparse grid is the sum of all product grids\n' );
  fprintf ( 1, '  of order LEVEL, with\n' );
  fprintf ( 1, '    0 <= LEVEL <= LEVEL_MAX.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL is the sum of the levels of the 1D rules,\n' );
  fprintf ( 1, '  the order of the 1D rule is 2^(LEVEL+1) - 1,\n' );
  fprintf ( 1, '  the region is [-1,1]^DIM_NUM.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For this kind of rule, there is complete nesting,\n' );
  fprintf ( 1, '  that is, a sparse grid of a given level includes\n' );
  fprintf ( 1, '  ALL the points on grids of lower levels.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   DIM: ' );

  for dim_num = dim_min : dim_max
    fprintf ( 1, '  %10d', dim_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   LEVEL_MAX\n' );
  fprintf ( 1, '\n' );

  for level_max = level_max_min : level_max_max
    fprintf ( 1, '    %4d', level_max );
    for dim_num = dim_min : dim_max
      point_num = sparse_grid_f2s_size ( dim_num, level_max );
      fprintf ( 1, '  %10d', point_num );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [stats,st] = lme_fit_EM(X,Zcols,y,ni,e)
% [stats,st] = lme_fit_EM(X,Zcols,y,ni,e)
%
% Linear mixed-effects estimation by expectation maximization.
%
% Input
% X: Ordered design Matrix (according to time for each subject).
% Zcols: Vector with the indices of the colums of X that will be considered
% as random effects.
% y: Ordered data vector (according to X).
% ni: Vector whose entries are the number of repeated measures for each
% subject (ordered according to X).
% e: Convergence epsilon. Default 10^-6;
%
% Output
% stats.Bhat: Estimated vector of the population regresion parameters.
% stats.CovBhat: Estimated covariance matrix of the population regresion
% parameters.
% stats.bihat: Estimated subject especific random effects.
% stats.Covbihat: Estimated covariance of the subject especific random
% coefficients.
% stats.phisqhat: Estimated within-subject variability.
% stats.SIGMA: Estimated marginal covariance matrices for each subject
% stacked in SIGMA.
% stats.W: Inverses of the estimated marginal covariance matrices for each
% subject stacked in W.
% stats.Dhat = Estimated random effects covariance matrix.
% stats.X: Design matrix.
% stats.Zcols: Same as Zcols in the input.
% stats.re: Residuals;
% stats.ni: Same as ni in the input.;
%
% $Revision: 1.2 $  $Date: 2015/01/06 17:14:57 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:14:57 $
%    $Revision: 1.2 $
% References: Bernal-Rusiel J.L., Greve D.N., Reuter M., Fischl B., Sabuncu
% M.R., 2012. Statistical Analysis of Longitudinal Neuroimage Data with Linear
% Mixed Effects Models, NeuroImage, doi:10.1016/j.neuroimage.2012.10.065.
%
tic;
if nargin < 4
    error('Too few inputs');
elseif nargin < 5
    e = 10^-6;
end;
try
    Z = X(:,Zcols);
catch Me
    error(['The colums of X specify in Zcols are not correct: ' Me.message]);
end
nit = 1000;
st = 1;
m = length(ni);
p = size(X,2);
q = length(Zcols);
n = sum(ni);
nimax = max(ni);
W = zeros(n,nimax);
bihat = zeros(q,m);
theta3 = zeros(q*q+1,3);

%Starting values
[D,phisq] = lme_fit_init(X,Zcols,y,ni);

%% Iterations
lreml0 = -10^10;
tf = true;
it = 0;
display('Starting Expectation Maximization iterations');
while tf
    it = it+1;
    %Computation of W = SIGMA^-1 and H.
    posi = 1; H = 0; Term = 0;
    scInvD = D\eye(q)*phisq;
    for i=1:m
        posf = posi+ni(i)-1;
        Zi = Z(posi:posf,:);
        Wi = (eye(ni(i))-Zi/(Zi'*Zi+scInvD)*Zi')/phisq;
        W(posi:posf,1:ni(i)) = Wi;
        Xi = X(posi:posf,:);
        Ti = Xi'*Wi;
        H = H + Ti*Xi;
        Term = Term + Ti*y(posi:posf);
        posi = posf+1;
    end;
    invH = H\eye(p);
   %Estimation
    Bhat = invH*Term;
    r = y-X*Bhat;
    Dhat = zeros(q,q);
    phisqhat = 0;
    posi = 1; lreml = 0;
    for i=1:m
        posf = posi+ni(i)-1;
        Zi = Z(posi:posf,:);
        Wi = W(posi:posf,1:ni(i));
        Xi = X(posi:posf,:);
        ri = r(posi:posf);
        bihat(:,i) = D*Zi'*Wi*ri;
        Pi = Wi - Wi*Xi*invH*Xi'*Wi;
        phisqhat = phisqhat + (ri-Zi*bihat(:,i))'*(ri-Zi*bihat(:,i))+phisq*...
                                             trace(eye(ni(i))-phisq*Pi);
        Dhat = Dhat+bihat(:,i)*bihat(:,i)'+D-D*Zi'*Pi*Zi*D;
        lreml = lreml + log(det(Wi))-ri'*Wi*ri;
        posi = posf+1;
    end;
    phisqhat = phisqhat/n;
    Dhat = Dhat/m;
    %Restricted log-likelihood
    lreml = 0.5*(lreml - log(det(H)));
    eps = abs(lreml-lreml0);
    display(['Likelihood at EM iteration ' num2str(it) ' : ' num2str(lreml)]);
    display(['Epsilon: ' num2str(eps)]);
    %Aitken acceleration
    theta3(:,1) = theta3(:,2);
    theta3(:,2) = theta3(:,3);
    theta3(:,3) = [vec(Dhat);phisqhat];
    if mod(it,10)==0
        lhat = mean((theta3(:,3)-theta3(:,2))./(theta3(:,2)-theta3(:,1)));
        if (lhat > 0) && (lhat < 1)
           theta3(:,3) = theta3(:,2) + (theta3(:,3)-theta3(:,2))/(1-lhat);
           Dhat = reshape(theta3(1:end-1,3),q,q);
           phisqhat = theta3(end,3);
        end;
    end;

    %Termination
    if (it==nit) || (eps<e)
        tf = false;
        SIGMA = zeros(n,nimax);
        Cbihat = zeros(q*m,q);
        posi = 1;
        for i=1:m
            posf = posi+ni(i)-1;
            Zi = Z(posi:posf,:);
            SIGMA(posi:posf,1:ni(i)) = Zi*D*Zi'+ eye(ni(i))*phisq;
            Wi = W(posi:posf,1:ni(i));
            Xi = X(posi:posf,:);
            Pi = Wi - Wi*Xi*invH*Xi'*Wi;
            Cbihat((i-1)*q+1:i*q,:) = D-D*Zi'*Pi*Zi*D;
            posi = posf+1;
        end;
       stats = struct('Bhat',Bhat,'CovBhat',invH,'bihat',bihat,...
             'Covbihat',Cbihat,'phisqhat',phisq,'SIGMA',SIGMA,'W',W,...
             'Dhat',D,'X',X,'Zcols',Zcols,'re',r,'ni',ni,'lreml',lreml);
        if it == nit
            st = 0;
            display(['Algorithm does not converge after ' num2str(nit)...
                                                        ' iterations!!!']);
        end;
    else
        lreml0 = lreml;
        D = Dhat;
        phisq = phisqhat;
    end;
end
toc;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
numNoisePoints = 40;
dots = rand(numNoisePoints,2);

numClus = 2;
clusCtr = rand(numClus,2);

stdMin = .04;
stdMax = .06;
clusStdX = rand(numClus)*stdMax+stdMin;
clusStdY = rand(numClus)*stdMax+stdMin;

sizeMin = 30;
sizeMax = 50;
clusSize = randi(sizeMax, numNoisePoints)+sizeMin;

for i = 1:numClus
    clusDotsX = normrnd(clusCtr(i,1),clusStdX(i),clusSize(i),1);
    clusDotsY = normrnd(clusCtr(i,2),clusStdY(i),clusSize(i),1);
    clusDots = [clusDotsX clusDotsY];
    dots = vertcat(dots, clusDots);
end

figure;
scatter(dots(:,1),dots(:,2),5,'MarkerFaceColor','k','MarkerEdgeColor','k');

%% kmeans 1
numK1 = 15;
idx_k1 = kmeans(dots,numK1);
colorsK1 = lines(numK1);
figure; hold on;
for i = 1:length(dots)
    scatter(dots(i,1),dots(i,2),5,'MarkerFaceColor',colorsK1(idx_k1(i),:),...
        'MarkerEdgeColor',colorsK1(idx_k1(i),:));
end

%% kmeans 2
numK2 = 5;
figure; hold on;
for j = 1:numK1
    dots_k1 = dots(idx_k1 == j,:);
    idx_k2 = kmeans(dots_k1,numK2);
    colorsK2 = lines(numK2*numK1);
    for i = 1:length(dots_k1)
        ci = (j-1)*numK1+idx_k2(i);
    scatter(dots_k1(i,1),dots_k1(i,2),5,'MarkerFaceColor',colorsK2(idx_k2(ci),:),...
        'MarkerEdgeColor',colorsK2(idx_k2(ci),:));
    end
end

author: xiuyechen

repo: FishExplorer


text:
function cnn_resnet_finetune(varargin)

if ~isempty(gcp('nocreate')),
    delete(gcp)
end
addpath('resnet');
opts = cnn_setup_environment();
opts.train.gpus =  [1] ;
% opts.train.gpus = [ 1 : 3 ]
opts.dataSet = 'Remat';
% opts.dataSet = 'hmdb51';
opts.aug        = 'stretch';
opts.border     = 0;
opts.pad        = 0;
opts.train.memoryMapFile = fullfile(tempdir, 'ramdisk', ['matconvnet' num2str(feature('getpid')) '.bin']) ;
addpath('network_surgery');
opts.dataDir = fullfile(opts.dataPath, opts.dataSet) ;
opts.splitDir = [opts.dataSet '_splits'];
opts.nSplit = 1 ;
opts.dropOutRatio = 0.5 ;
opts.train.cheapResize = 0 ;
opts.inputdim  = [ 224,  224, 3] ;

opts.train.batchSize = 196  ;
opts.train.numSubBatches =  ceil(8 / max(numel(opts.train.gpus),1));

opts.train.epochFactor = 10 ;
opts.train.augmentation = 'borders25';

opts.train.backpropDepth =     cell(1, 2);
opts.train.backpropDepth(:) = {'pool5'};
opts.train.learningRate =  [1e-2*ones(1, 2) 1e-2*ones(1, 3) 1e-3*ones(1, 3) 1e-4*ones(1, 3)] ;
if strcmp(opts.dataSet, 'hmdb51')
  opts.train.learningRate =  [1e-2*ones(1, 2) 1e-2*ones(1, 1) 1e-3*ones(1, 1) 1e-4*ones(1, 1)] ;
end

model = ['img-res50-' opts.train.augmentation '-bs=' num2str(opts.train.batchSize) ...
  '-split' num2str(opts.nSplit) '-dr' num2str(opts.dropOutRatio)];

if  strfind(model, 'vgg16');
  baseModel = 'imagenet-vgg-verydeep-16.mat' ;
  opts.train.learningRate =  [1e-3*ones(1, 3) 5e-4*ones(1, 5) 5e-5*ones(1,2) 5e-6*ones(1,2)]  ;
  opts.train.backpropDepth =     cell(1, 3);
  opts.train.backpropDepth(:) = {'layer37'};
  opts.train.batchSize = 128  ;
  opts.train.numSubBatches =  ceil(16 / max(numel(opts.train.gpus),1));
elseif strfind(model, 'vgg-m');
  baseModel = 'imagenet-vgg-m-2048.mat' ;
elseif strfind(model, 'res152');
  baseModel = 'imagenet-resnet-152-dag.mat' ;
elseif strfind(model, 'res101');
  baseModel = 'imagenet-resnet-101-dag.mat' ;
elseif strfind(model, 'res50');
  baseModel = 'imagenet-resnet-50-dag.mat' ;
  opts.train.numSubBatches =  ceil(32 / max(numel(opts.train.gpus),1));
else
  error('Unknown model %s', model) ;
end
opts.model = fullfile(opts.modelPath,baseModel) ;
opts.expDir = fullfile(opts.dataDir, [opts.dataSet '-' model]) ;
opts.imdbPath = fullfile(opts.dataDir, [opts.dataSet '_resnet_split' num2str(opts.nSplit) 'imdb.mat']);

opts.train.plotDiagnostics = 0 ;
opts.train.continue = 1 ;
opts.train.prefetch = 1 ;
opts.train.expDir = opts.expDir ;
opts.train.numAugments = 1;
opts.train.frameSample = 'random';
opts.train.nFramesPerVid = 1;
opts.train.uniformAugments  = false;

[opts, varargin] = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                   Database initialization
% -------------------------------------------------------------------------
if exist(opts.imdbPath)
  imdb = load(opts.imdbPath) ;
else
  imdb = cnn_Remat_resnet_setup_data('dataPath', opts.dataPath, 'flowDir',opts.flowDir, ...
    'dataSet', opts.dataSet, 'nSplit', opts.nSplit) ;
  save(opts.imdbPath, '-struct', 'imdb', '-v6') ;
end
nClasses = length(imdb.classes.name);
if ~exist(opts.model)
  fprintf('Downloading base model file: %s ...\n', baseModel);
  mkdir(fileparts(opts.model)) ;
  urlwrite(...
  ['http://www.vlfeat.org/matconvnet/models/' baseModel], ...
    opts.model) ;
end

net = load(opts.model);
if isfield(net, 'net'), net=net.net;end

if isstruct(net.layers)
  % replace 1000-way imagenet classifiers
  for p = 1 : numel(net.params)
    sz = size(net.params(p).value);
    if any(sz == 1000)
      sz(sz == 1000) = nClasses;
      fprintf('replace classifier layer of %s\n', net.params(p).name);
      if numel(sz) > 2
         net.params(p).value = 0.01 * randn(sz,  class(net.params(p).value));
      else
         net.params(p).value = zeros(sz,  class(net.params(p).value));
      end
    end
  end

  net.meta.normalization.border = [256 256] - net.meta.normalization.imageSize(1:2);
  net = dagnn.DagNN.loadobj(net);
  if strfind(model, 'bnorm')
    net = insert_bnorm_layers(net) ;
  end
else
%   net=vl_simplenn_tidy(net);
  if isfield(net, 'meta'),
    netNorm = net.meta.normalization;
  else
    netNorm = net.normalization;
  end
  if(netNorm.imageSize(3) == 3) && ~isempty(strfind(opts.model, 'imagenet'))
    netNorm.border = [240 240] - netNorm.imageSize(1:2);
    net = replace_last_layer(net, [1 2], [1 2], nClasses, opts.dropOutRatio);
  end
  if strfind(model, 'bnorm')
    net = insert_bnorm_layers(net) ;
  end

  net = dagnn.DagNN.fromSimpleNN(net) ;

end

net = dagnn.DagNN.setLrWd(net);

net.renameVar(net.vars(1).name, 'input');

if ~isnan(opts.dropOutRatio)
  dr_layers = find(arrayfun(@(x) isa(x.block,'dagnn.DropOut'), net.layers)) ;
  if ~isempty(dr_layers)
    if opts.dropOutRatio > 0
      for i=dr_layers, net.layers(i).block.rate = opts.dropOutRatio; end
    else
      net.removeLayer({net.layers(dr_layers).name});
    end
  else
    if opts.dropOutRatio > 0
      pool5_layer = find(arrayfun(@(x) isa(x.block,'dagnn.Pooling'), net.layers)) ;
      conv_layers = pool5_layer(end);
      for i=conv_layers
        block = dagnn.DropOut() ;   block.rate = opts.dropOutRatio ;
        newName = ['drop_' net.layers(i).name];

        net.addLayer(newName, ...
          block, ...
          net.layers(i).outputs, ...
          {newName}) ;

        for l = 1:numel(net.layers)-1
          for f = net.layers(i).outputs
             sel = find(strcmp(f, net.layers(l).inputs )) ;
             if ~isempty(sel)
              [net.layers(l).inputs{sel}] = deal(newName) ;
             end
          end
        end
      end
    end
  end
end
net.layers(~cellfun('isempty', strfind({net.layers(:).name}, 'err'))) = [] ;

opts.train.derOutputs = {} ;
for l=numel(net.layers):-1:1
  if isa(net.layers(l).block, 'dagnn.Loss') && isempty(strfind(net.layers(l).name, 'err'))
      opts.train.derOutputs = {opts.train.derOutputs{:}, net.layers(l).outputs{:}, 1} ;
  end
  if isa(net.layers(l).block, 'dagnn.SoftMax')
    net.removeLayer(net.layers(l).name)
    l = l - 1;
  end
end

if isempty(opts.train.derOutputs)
  net = dagnn.DagNN.insertLossLayers(net, 'numClasses', nClasses) ;
  fprintf('setting derivative for layer %s \n', net.layers(end).name);
  opts.train.derOutputs = {opts.train.derOutputs{:}, net.layers(end).outputs{:}, 1} ;
end

lossLayers = find(arrayfun(@(x) isa(x.block,'dagnn.Loss') && strcmp(x.block.loss,'softmaxlog'),net.layers));

net.addLayer('top1error', ...
             dagnn.Loss('loss', 'classerror'), ...
             net.layers(lossLayers(end)).inputs, ...
             'top1error') ;

net.addLayer('top5error', ...
             dagnn.Loss('loss', 'topkerror', 'opts', {'topK', 5}), ...
             net.layers(lossLayers(end)).inputs, ...
             'top5error') ;

net.print() ;
net.rebuild() ;

net.meta.normalization.rgbVariance = [];
net.meta.normalization.averageImage = mean(mean(net.meta.normalization.averageImage, 1), 2);
opts.train.train = find(ismember(imdb.images.set, [1])) ;
opts.train.train = repmat(opts.train.train,1,opts.train.epochFactor);
% opts.train.valmode = '250samples';

opts.train.valmode = '30samples'

opts.train.denseEval = 1;
net.conserveMemory = 1 ;
 fn = getBatchFn(opts, net.meta);
%fn = getBatchWrapper_ucf101_rgbflow(net.meta.normalization, opts.numFetchThreads, opts.train) ;
[info] = cnn_resnet_train_dag(net, imdb, fn, opts.train) ;

end

% -------------------------------------------------------------------------
function fn = getBatchFn(opts, meta)
% -------------------------------------------------------------------------
bopts.numThreads = opts.numFetchThreads ;
bopts.pad = opts.pad ;
bopts.border = opts.border ;
bopts.transformation = opts.aug ;
bopts.imageSize = meta.normalization.imageSize ;
bopts.averageImage = meta.normalization.averageImage ;
bopts.rgbVariance = meta.normalization.rgbVariance ;
% bopts.transformation = meta.augmentation.transformation ;

fn = @(x,y) getSimpleNNBatch(bopts,x,y) ;
end

% -------------------------------------------------------------------------
function [im,labels] = getSimpleNNBatch(opts, imdb, batch)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
isVal = ~isempty(batch) && imdb.images.set(batch(1)) ~= 1 ;

if ~isVal
  % training
  im = cnn_resnet_get_batch(images, opts, ...
                              'prefetch', nargout == 0) ;
else
  % validation: disable data augmentation
%   opts.border=0;
  im = cnn_resnet_get_batch(images, opts, ...
                              'prefetch', nargout == 0, ...
                              'transformation', 'none') ;
end

if nargout > 0
  labels = imdb.images.class(batch) ;
end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [volumePVC] = PVEcorrect(volume,nIter,waveletName)
% -------------------------------------------------------------------------
% function [volumePVC] = PVEcorrect(volume,nIter,waveletName)
% -------------------------------------------------------------------------
% DESCRIPTION:
% Apply partial-volume effect (PVE) correction of an input PET volume
% using the methodology developed in ref. [1].
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Boussion, N. et al. (2009). Incorporation of wavelet-based denoising
%     in iterative deconvolution for partial volume correction in
%     whole-body PET imaging. Eur J Nucl Med Mol Imaging, 36(7), 1064-1075.
% -------------------------------------------------------------------------
% INPUTS:
% - volume: 3D array representing the input PET volume to correct for PVE.
% - nIter: Number of iterations in the deconvolution process of the PVE
%          correction (see ref. [1]). If set to a string 'compute', the
%          algorithm will determine the optimal number of iterations in
%          terms of residuals up to a maximum of 10.
% - waveletName: (optional). MATLAB name of the type of wavelet used in the
%                denoising part of the PVE correction (see ref. [1]).
%                Default is 'bior3.5'.
% -------------------------------------------------------------------------
% OUTPUTS:
% - volumePVC: 3D array representing the input PET volume corrected for
%              partial volume effects.
% -------------------------------------------------------------------------
% AUTHOR(S):
% - Andre Diamant <adboustead@gmail.com>
% - Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

% INTIALIZATION
volumePVC = zeros(size(volume));
psf = nonIsotropicGaussianPSF(1,4);

if nargin == 3
    global wavelet_name
    wavelet_name = waveletName;
end

% PVE CORRECTION
if strcmp(nIter,'compute')
    [~,residuals] = deconvlucydenoiseNew(volume,psf,10);
    [~,nIter] = min(abs(residuals));
    fprintf('*** FOUND NUMBER OF ITERATIONS TO BE %.0f ***\n',nIter)
end
[volumePVC,~] = deconvlucydenoiseNew(volume,psf,nIter);

if nargin == 3
    clear global wavelet_name
end

end

author: mvallieres

repo: radiomics


text:
% This is a template script that runs the HMM on a fictional fMRI data set
% using a number of different settings, involving how many PCA components we use to
% describe the data, and how many HMM states we set up for the HMM inference.
% It assumes that there are 2 different
% conditions (for example rest and some task). It performs permutation
% testing on the HMM results for each of the HMM settings, and plots the results.
% More specifically, it tests for fractional occupancy differences and state
% switching rate differences between the conditions, using GLM-based
% permutation testing (function hmmtest.m)
%
% This script can be adapted to block-designed real data with any number of
% conditions by configuring the variable Y.
% Note that the data contain no dynamics, and this
% script is only for the purposes of demonstration
%
% It assumes that the toolbox paths are in the right place.
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2017)

%% Prepare data

N = 25; % subjects
Q = 4; % sessions per subject
ttrial = 500; % time points
nregions = 50; % regions or voxels
Y = zeros(N*Q,2); % design matrix with conditions
X = randn(Q*N*ttrial,nregions); % all data concatenated
T = ttrial * ones(N*Q,1);  % length of data for each session
Tsubject = Q*ttrial * ones(N,1);  % length of data for each subject

% Functional connectivity profile for condition 1
FC_condition1 = randn(nregions);
FC_condition1 = FC_condition1' * FC_condition1;
% Functional connectivity profile for condition 2
FC_condition2 = randn(nregions);
FC_condition2 = FC_condition2' * FC_condition2;
% Note that there is no difference in the mean activity between conditions

for j1 = 1:N
    for j2 = 1:Q
        t = (1:ttrial) + (j2-1)*ttrial + Q*(j1-1)*ttrial;
        n = (j1-1)*Q + j2;
        if rem(j2,2)
            Y(n,1) = 1;
            X(t,:) = X(t,:) * FC_condition1;
        else
            Y(n,2) = 1;
            X(t,:) = X(t,:) * FC_condition2;
        end
        for i = 1:nregions % do some smoothing
            X(t,i) = smooth(X(t,i));
        end
    end
end

options = struct();
e = explainedvar_PCA(X,T,options); % how much variance PCA explains on the data
pc1 = find(e>0.5,1); % get no. of PCA components to explain 50% of variance
pc2 = find(e>0.6,1); % get no. of PCA components to explain 60% of variance
pc3 = find(e>0.7,1); % get no. of PCA components to explain 70% of variance
pc4 = find(e>0.8,1); % get no. of PCA components to explain 80% of variance
pc5 = find(e>0.9,1); % get no. of PCA components to explain 90% of variance
pc6 = 0; % no pca
number_pca_components = [pc1 pc2 pc3 pc4 pc5 pc6];
number_states = 4:6;

template_configuration = struct();
template_configuration.order = 0;
template_configuration.dropstates = 0;
template_configuration.verbose = 0;
template_configuration.cyc = 50;
template_configuration.initcyc = 5;

i = 1;
configurations = {}; % parameters of the HMM
% Gaussian distribution with  mean and full covariance
for k = number_states
    for pca_comp = number_pca_components
        configurations{i} = template_configuration;
        configurations{i}.K = k;
        configurations{i}.pca = pca_comp;
        configurations{i}.zeromean = 0;
        configurations{i}.covtype = 'full';
        i = i + 1;
    end
end
% Gaussian distribution with just covariance
for k = number_states
    for pca_comp = number_pca_components
        configurations{i} = template_configuration;
        configurations{i}.K = k;
        configurations{i}.pca = pca_comp;
        configurations{i}.zeromean = 1;
        configurations{i}.covtype = 'full';
        i = i + 1;
    end
end
% Gaussian distribution with just mean
for k = number_states
    for pca_comp = number_pca_components
        configurations{i} = template_configuration;
        configurations{i}.K = k;
        configurations{i}.pca = pca_comp;
        configurations{i}.zeromean = 0;
        configurations{i}.covtype = 'sharedfull';
        i = i + 1;
    end
end

%% Run HMM and perform statistical testing on the HMM results vs the conditions

L = length(configurations);
Gamma = cell(length(configurations),1);
test_group = cell(length(configurations),1);

options_test = struct();
options_test.subjectlevel = 0;
options_test.Nperm = 1000;

for i = 1:length(configurations)
    [hmm,Gamma{i}] = hmmmar(X,T,configurations{i});
    t = hmmtest(Gamma{i},T,Tsubject,Y,options_test,hmm);
    test_group{i} = t.grouplevel; % only doing group-level testing
    disp([num2str(i) ' of ' num2str(L)])
end

%% Plot

L1 = length(number_pca_components) * length(number_states);
L2 = length(number_states);
L3 = length(number_pca_components);
titles = {'Gaussian','Gaussian no Mean','Gaussian no Cov'};

% fractional occupancies
figure(1);clf(1)
ii = 1;
for i = 1:3
    for j = 1:L2
        subplot(3,L2,ii)
        K = number_states(j);
        M = zeros(L3,K);
        for l = 1:L3
            ind = (i-1)*L1 + (j-1)*L3 + l;
            M(l,:) = test_group{ind}.p_fractional_occupancy';
        end
        imagesc(M,[0 0.25]); colorbar
        title([titles{i} ' K=' num2str(K)])
        ylabel('PCA components');xlabel('States')
        ii = ii + 1;
    end
end

% switching rate
figure(2);clf(2)
ii = 1;
for i = 1:3
    subplot(1,3,i)
    M = zeros(L3,L2);
    for j = 1:L2
        for l = 1:L3
            ind = (i-1)*L1 + (j-1)*L2 + l;
            M(l,j) = test_group{ind}.p_switching_rate;
        end
    end
    imagesc(M,[0 0.25]); colorbar
    title([titles{i}])
    ylabel('PCA components');xlabel('No. States')
end

author: OHBA-analysis

repo: HMM-MAR


text:
classdef ParadigmMTCSP < ParadigmDataflowSimplified
    % Experimental paradigm for all-frequency Common Spatial Patterns.
    %
    % The basic idea is to calculate CSP for each covariance matrix in the cross-spectrum,
    % and to use multi-taper spectral estimation to ensure an optimal tradeoff between
    % spectral precision and estimation noise. The default classifier is sparse logistic
    % regression with elastic-net penalty.
    %
    % This paradigm also implements a second approach in which the cross-spectrum is not spatially
    % filtered but directly submitted to the classifier (Disciplined Cross-Spectral Regression).
    %
    % Name:
    %   Multi-Taper CSP
    %
    %                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
    %                           2013-04-26

    methods

        function defaults = preprocessing_defaults(self)
            % define the default pre-processing parameters of this paradigm
            defaults = {'FIRFilter',{[0.5 1],'highpass'}, 'EpochExtraction',[0.5 3.5],'Resampling',200};
        end

        function defaults = machine_learning_defaults(self)
            defaults = {'logreg','variant',{'lars','ElasticMixing',0.25}};
        end

        function model = feature_adapt(self,varargin)
            % adapt a feature representation using the CSP algorithm
            args = arg_define(varargin, ...
                arg_norep('signal'), ...
                arg({'timewnds','TimeWindows'},[],[],'Time windows of interest. Matrix containing one row for the start and end of each time window for which CSP patterns shall be computed. Values in seconds. If both this and the freqwnds parameter are non-empty, they should have the same number of rows.'), ...
                arg({'winfunc','WindowFunction'},'rect',{'bartlett','barthann','blackman','blackmanharris','flattop','gauss','hamming','hann','kaiser','lanczos','nuttall','rect','triang'},'Type of window function. Typical choices are rect (rectangular), hann, gauss, blackman and kaiser.'),...
                arg({'winparam','WindowParameter','param'},[],[],'Parameter of the window function. This is mandatory for cheb, kaiser and tukey and optional for some others.','shape','scalar'), ...
                arg_sub({'spectral_estimation','SpectralEstimation'},{},@utl_calc_crossspec, 'Spectral estimation parameters.'), ...
                arg({'patterns','PatternPairs'},3,uint32([1 1 64 10000]),'CSP patterns per frequency (times two).'), ...
                arg({'whycsp','SkipCSP','WhyCSP'},false,[],'Classify cross-spectrum directly. This results in much higher-dimensional features, but can be approached with appropriately regularized classifiers (see ml_trainproximal).'), ...
                arg({'normalize_spectrum','NormalizeSpectrum'},false,[],'Normalize the spectrum. Recommended if using sophisticated regularized classifiers.'), ...
                arg({'logtransform','LogTransform'},false,[],'Log-transform output. Log-transformed spectra are more likely to be separable by a linear classifier.'), ...
                arg({'vectorize_features','VectorizeFeatures'},true,[],'Vectorize the features. For compatibility with basic classifiers.'));

            if args.signal.nbchan == 1
                error('Multi-taper CSP does intrinsically not support single-channel data (it is a spatial filter).'); end
            if args.signal.nbchan < args.patterns
                error('Multi-taper CSP prefers to work on at least as many channels as you request output patterns. Please reduce the number of pattern pairs.'); end
            if isempty(args.timewnds)
                args.timewnds = struct(); end
            [C,dummy,T] = size(args.signal.data); %#ok<NASGU,ASGLU>
            if args.whycsp
                % shortcut
                for w = size(args.timewnds,1):-1:1
                    time_args{w} = arg_report('vals',@flt_window,{'time',{args.timewnds(w,:),args.winfunc,args.winparam}}); end
                model = struct('spec_args',{args.spectral_estimation}, 'time_args',{time_args},'chanlocs',{args.signal.chanlocs},'vectorize_features',{args.vectorize_features},'whycsp',{args.whycsp},'normalize_spectrum',{args.normalize_spectrum});
            else
                covar = {}; mean_covar = {}; weighted_covar = {};
                % for each time window...
                for w = size(args.timewnds,1):-1:1
                    if length(unique([args.signal.event.target]))>2
                        % SPoC version
                        time_args{w} = arg_report('vals',@flt_window,{'time',{args.timewnds(w,:),args.winfunc,args.winparam}});
                        % calc weighted and average cross-spectra
                        [mean_covar{w},weighted_covar{w}] = hlp_diskcache('featuremodels',@utl_calc_crossspec,args.spectral_estimation,'sum_weights',[args.signal.epoch.target],'signal',exp_eval_optimized(flt_window('signal',args.signal,time_args{w})));
                        mean_covar{w}(~isfinite(mean_covar{w})) = 0; weighted_covar{w}(~isfinite(weighted_covar{w})) = 0;
                        % calculate spatial filters for each frequency
                        for f=size(mean_covar{w},1):-1:1
                            [V,D] = eig(reshape(weighted_covar{w}(f,:,:),C,C),reshape(mean_covar{w,1}(f,:,:),C,C)); %#ok<NASGU>
                            P = inv(V);
                            % retain k best filters/patterns at both ends of the eigenvalue spectrum
                            filters(w,f,:,:) = real(V(:,[1:args.patterns end-args.patterns+1:end]));
                            patterns(w,f,:,:) = real(P([1:args.patterns end-args.patterns+1:end],:))';
                        end
                    else
                        % CSP version
                        for k=1:2
                            subset = exp_eval_optimized(set_picktrials(args.signal,'rank',k));
                            % pre-parse arguments for flt_window and flt_spectrum (for fast subsequent online use)
                            time_args{w} = arg_report('vals',@flt_window,{'time',{args.timewnds(w,:),args.winfunc,args.winparam}});
                            % calc cross-spectrum for the given windowed data subset
                            covar{w,k} = hlp_diskcache('featuremodels',@utl_calc_crossspec,args.spectral_estimation,'signal',exp_eval_optimized(flt_window('signal',subset,time_args{w})));
                            covar{w,k}(~isfinite(covar{w,k})) = 0;
                        end
                        % solve a CSP instance for each frequency
                        for f=size(covar{w,1},1):-1:1
                            [V,D] = eig(reshape(covar{w,1}(f,:,:),C,C),reshape(covar{w,1}(f,:,:),C,C)+reshape(covar{w,2}(f,:,:),C,C)); %#ok<NASGU>
                            P = inv(V);
                            filters(w,f,:,:) = real(V(:,[1:args.patterns end-args.patterns+1:end]));
                            patterns(w,f,:,:) = real(P([1:args.patterns end-args.patterns+1:end],:))';
                        end
                    end
                end
                model = struct('filters',{filters},'patterns',{patterns},'time_args',{time_args},'spec_args',{args.spectral_estimation}, 'covar',{covar}, 'mean_covar',{mean_covar}, 'weighted_covar',{weighted_covar}, 'chanlocs',{args.signal.chanlocs},'vectorize_features',{args.vectorize_features},'whycsp',{args.whycsp},'normalize_spectrum',{args.normalize_spectrum},'logtransform',{args.logtransform});
            end
            global tracking; %#ok<TLEV>
            tracking.inspection.signal = args.signal;
            tracking.inspection.chanlocs = args.signal.chanlocs;
        end

        function features = feature_extract(self,signal,featuremodel)
            for w = length(featuremodel.time_args):-1:1
                % extract time window
                wnd = exp_eval_optimized(flt_window('signal',signal,featuremodel.time_args{w}));
                % extract cross-spectral features (note: gigantic!)
                if featuremodel.whycsp
                    % W x F x C x C x T
                    features(w,:,:,:,:) = utl_calc_crossspec(featuremodel.spec_args,'signal',wnd,'feature_filters',false);
                    if featuremodel.normalize_spectrum
                        nfreqs = size(features,2);
                        freqs = featuremodel.spec_args.freqwnd;
                        freqs = freqs(1):(freqs(2)-freqs(1))/(nfreqs-1):freqs(2);
                        features = bsxfun(@times,features,max(1,freqs));
                    end
                else
                    % F x W x P x T
                    if onl_isonline
                        features(:,w,:,:) = utl_calc_crossspec(featuremodel.spec_args,'signal',wnd,'feature_filters',squeeze(featuremodel.filters(w,:,:,:)));
                    else
                        features(:,w,:,:) = hlp_diskcache('features',@utl_calc_crossspec,featuremodel.spec_args,'signal',wnd,'feature_filters',squeeze(featuremodel.filters(w,:,:,:)));
                    end
                    if featuremodel.normalize_spectrum
                        nfreqs = size(features,1);
                        freqs = featuremodel.spec_args.freqwnd;
                        freqs = freqs(1):(freqs(2)-freqs(1))/(nfreqs-1):freqs(2);
                        features = bsxfun(@times,features,max(1,1./freqs'));
                    end
                    if featuremodel.logtransform
                        features = log(features); end
                end
            end
            % apply minimal conditioning to features
            features = real(features);
            features(~isfinite(features)) = 0;
            % do final vectorization if desired
            if featuremodel.vectorize_features
                features = reshape(features,[],signal.trials)'; end
        end

        function visualize_model(self,parent,featuremodel,predictivemodel,varargin) %#ok<*INUSD>
            % no visualization yet
        end

        function layout = dialog_layout_defaults(self)
            % define the default configuration dialog layout
            layout = {'SignalProcessing.Resampling.SamplingRate', 'SignalProcessing.EpochExtraction', '', ...
                'Prediction.FeatureExtraction.TimeWindows','Prediction.FeatureExtraction.WindowFunction', '', ...
                'Prediction.FeatureExtraction.SpectralEstimation.FrequencyRange', ...
                'Prediction.FeatureExtraction.SpectralEstimation.TimeBandwidth', ...
                'Prediction.FeatureExtraction.SpectralEstimation.SubsampleSpectrum', ...
                'Prediction.FeatureExtraction.SpectralEstimation.RobustEstimation', '' ...
                'Prediction.FeatureExtraction.PatternPairs', 'Prediction.FeatureExtraction.VectorizeFeatures', '', ...
                'Prediction.MachineLearning.Learner'};
        end

        function tf = needs_voting(self)
            tf = false;
        end
    end
end

repo: REST


text:
function [imgOut, BTMO_segments, BTMO_which_operator] = BanterleTMO(img, BTMO_segments, BTMO_rescale)
%
%
%        [imgOut, BTMO_segments, BTMO_which_operator] = BanterleTMO(img, BTMO_segments)
%
%
%       Input:
%           -img: an HDR image with calibrated data in cd/m^2.
%           Note that this algorithm was tested with values
%           from 0.015cd/m^2 to 3,000 cd/m^2
%           -BTMO_segments: a segmented image, each value in a segment is a
%           dynamic range zone; i.e. integer values in [-6,9]. If it is not
%           provided this will be computed with the function CreateSegments
%           -BTMO_rescale:
%
%       Output:
%           -imgOut: the tone mapped image using the HybridTMO
%           -BTMO_segments: the segmentation output; the input image is
%           segmented into different zones of dynamic range
%           -BTMO_which_operator: is a string which outputs if only
%           DragoTMO is used, if only Reinhard is used, or if the full
%           BanterleTMO is used.
%
%       This TMO is an hybrid operator which merges different
%       Tone Mapping Operators: DragoTMO and ReinhardTMO
%
%     Copyright (C) 2013-15  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%     The paper describing this technique is:
%     "Dynamic Range Compression by Differential Zone Mapping Based on
%     Psychophysical Experiments"
% 	  by Francesco Banterle, Alessandro Artusi, Elena Sikudova,
%     Thomas Edward William Bashford-Rogers, Patrick Ledda, Marina Bloj, Alan Chalmers
%     in ACM Symposium on Applied Perception (SAP) - August 2012
%
%

checkNegative(img);

if(~exist('BTMO_rescale', 'var'))
    BTMO_rescale = 1;
end

L = lum(img);

if(BTMO_rescale)
    L_max = MaxQuart(L, 0.999);

    if(L_max > 0.0)
        img = img / L_max * 3000.0;
        img = ClampImg(img, 0.015, 3000.0);
        L = lum(img);
    end
end

indx = find(L > 3000);
if(~isempty(indx))
    warning(['The input image has values over 3,000 cd/m^2. ', ...
        'These values were not tested in the original paper.']);
end

indx = find(L < 0.015);
if(~isempty(indx))
    warning(['The input image has values under 0.015 cd/m^2. ', ...
        'These values were not tested in the original paper.']);
end

%segmentation
if(~exist('BTMO_segments', 'var'))
    BTMO_segments = CreateSegments(img);
else
    if(isempty(BTMO_segments))
        BTMO_segments = CreateSegments(img);
    else
        BTMO_segments = round(BTMO_segments);
    end
end

%TMO look-up table for determing the best
%TMO depending on the luminance zone. These
%values were extracted from a psychophysical
%experiment.
% 0 ---> Drago et al. 2003
% 1 ---> Reinhard et al. 2002
LumZone     = [-2, -1, 0, 1, 2, 3, 4];
TMOForZone =  [ 0,  0, 1, 0, 1, 0, 0];

%mask
mask = zeros(size(BTMO_segments));

for i=1:length(LumZone)
    mask(BTMO_segments == LumZone(i)) = TMOForZone(i);
end

imwrite(mask, 'mask.png');
%check for Drago et al.'s operator
indx0 = find(mask == 1);

%check for Reinhard et al.'s operator
indx1 = find(mask == 0);

%are both TMOs used?
if(~isempty(indx0) && ~isempty(indx1)) %pyramid blending with gamma encoding
    img_dra_tmo = DragoTMO(img);
    img_rei_tmo = ReinhardTMO(img, 0.5, -1, 'local', -1);

    gamma = 2.2;
    invGamma = 1.0 / gamma;
    imgA   = img_rei_tmo.^invGamma;
    imgB   = img_dra_tmo.^invGamma;

    imwrite(imgA, 'banterle_rei.png');
    imwrite(imgB, 'banterle_dra.png');
    imwrite(mask, 'banterle_mask.png');

    %removing gamma for linear output
    imgOut = pyrBlend(imgA, imgB, mask).^gamma;
    imwrite(imgOut.^invGamma, 'banterle_mix.png');

    BTMO_which_operator = 'BanterleTMO';
else
    %Only Reinhard et al.'s operator?
    length(indx1)
    length(indx0)
    if(isempty(indx1))
        imgOut = ReinhardTMO(img, -1, -1, 'bilateral', -1);
        BTMO_which_operator = 'ReinhardTMO';
        disp('The BanterleTMO is using ReinhardTMO only');
    end

    %Only Drago et al.'s operator?
    if(isempty(indx0))
        imgOut = DragoTMO(img);
        BTMO_which_operator = 'DragoTMO';
        disp('The BanterleTMO is using DragoTMO only');
    end
end

end

author: banterle


text:
function x = NormL1NN_project(x,weights,tau)
% Projection onto the non-negative part of the L1 ball

%idx = (x > 0);
%x(idx) = NormL1_project(x(idx),weights(idx),tau);
x(x < 0) = 0;
x = NormL1_project(x,weights,tau);

author: mpf


text:
function [points,weights] = FibonacciGaussian(N,mu,P)
%%FIBONACCIGAUSSIAN Creates a grid of points based on a transformation of a
%                   Fibonacci grid as given in [1].
%
%INPUT:
% N: The number of points to sample.
% mu: The D-by-1 vector indicating the Gaussian mean.
% P: The D-by-D covariance matrix. This must be positive definite.
%
%OUTPUT:
% points: A D-by-N matrix containing sample points for the given Gaussian.
% weights: A collection of weights representing the Gaussian PDF at each
%          sample point.
%
%EXAMPLE 1: Generates a collection of sample points and their respective
%           Gaussian weights in 1D and plots them.
% mu = 5;
% P = 2;
% [points, weights] = FibonacciGaussian(100,mu,P);
% stem(points(1,:),weights,'k',"filled","MarkerEdgeColor",'b')
%
%EXAMPLE 2: Generates a collection of sample points and their respective
%         Gaussian weights in 2D and plots them.
% mu = [5;3];
% P = [2,-0.7;-0.7,1];
% [points, weights] = FibonacciGaussian(100,mu,P);
% scatter(points(1,:),points(2,:),'k',"filled","MarkerEdgeColor",'b',"CData",weights)
% hold on
% drawEllipse(mu,P\eye(2),[],'r',"LineWidth",2)
% hold off
%
%EXAMPLE 3: Generates a collection of sample points and their respective
%           Gaussian weights in 3D and plots them.
% mu = [5;3;4];
% P = [2,-0.7,0.1;-0.7,1,-0.3;0.1,-0.3,2];
% [points, weights] = FibonacciGaussian(100,mu,P);
% scatter3(points(1,:),points(2,:),points(3,:),'k',"filled","MarkerEdgeColor",'b',"CData",1e3*weights)
% hold on
% drawEllipse(mu,P\eye(3),[],"FaceAlpha",0.1,"EdgeAlpha",0.1)
% hold off
%
%REFERENCE:
%[1] D. Frisch and U. D. Hanebeck, "Deterministic gaussian sampling with
%    generalized fibonacci grids," in 2021 IEEE 24th International
%    Conference on Information Fusion (FUSION), IEEE, 2021, pp. 1-8.
%
%December 2022 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.
D = size(mu,1);
points = FibonacciGrid(N,D);
points = sqrt(2)*erfinv(2*points-1);
nu = sqrt(sum(points.^2,2)/N);
for idx = 1:D
    points(idx,:) = points(idx,:)/nu(idx);
end
[V,E] = eig(P);
points = V*sqrt(abs(E))*points+mu;
weights = GaussianD.PDF(points,mu,P);
end
%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [mustUU, pos_mustUU, mustUU_linear, pos_mustUU_linear] = findMustUUWithGAMS(model, minFluxesW, maxFluxesW, constrOpt, excludedRxns, mustSetFirstOrder, solverName, runID, outputFolder, outputFileName, printExcel, printText, printReport, keepInputs, keepGamsOutputs, verbose)
% This function runs the second step of `optForce`, that is to solve a
% bilevel mixed integer linear programming  problem to find a second order
% `MustUU` set.
%
% USAGE:
%
%    [mustUU, pos_mustUU, mustUU_linear, pos_mustUU_linear] = findMustUUWithGAMS(model, minFluxesW, maxFluxesW, varargin)
%
% INPUTS:
%    model:                      (structure) a metabolic model with at
%                                least the following fields:
%
%                                  * .rxns - Reaction IDs in the model
%                                  * .mets - Metabolite IDs in the model
%                                  * .S -    Stoichiometric matrix (sparse)
%                                  * .b -    RHS of `Sv = b` (usually zeros)
%                                  * .c -    Objective coefficients
%                                  * .lb -   Lower bounds for fluxes
%                                  * .ub -   Upper bounds for fluxes
%    minFluxesW:                 (double array of size `n_rxns x 1`) minimum
%                                fluxes for each reaction in the model for
%                                wild-type strain. This can be obtained by
%                                running the function `FVAOptForce`. E.g.:
%                                `minFluxesW = [-90; -56];`
%    maxFluxesW:                 (double array of size `n_rxns x 1`) maximum
%                                fluxes for each reaction in the model for
%                                wild-type strain. This can be obtained by
%                                running the function `FVAOptForce`. E.g.:
%                                `maxFluxesW = [90; 56];`
%
% OPTIONAL INPUTS:
%    constrOpt:                  (Structure) structure containing
%                                additional contraints. Include here only
%                                reactions whose flux is fixed, i.e.,
%                                reactions whose lower and upper bounds
%                                have the same value. Do not include here
%                                reactions whose lower and upper bounds
%                                have different values. Such contraints
%                                should be defined in the lower and upper
%                                bounds of the model. The structure has the
%                                following fields:
%
%                                  * .rxnList - Reaction list (cell array)
%                                  * .values -  Values for constrained
%                                    reactions (double array)
%                                    E.g.: `struct('rxnList', {{'EX_gluc', 'R75', 'EX_suc'}}, 'values', [-100, 0, 155.5]');`
%    excludedRxns:               (cell array) Reactions to be excluded to
%                                the `MustUU` set. This could be used to
%                                avoid finding transporters or exchange
%                                reactions in the set. Default = empty.
%    mustSetFirstOrder:          (cell array) Reactions that belong to
%                                `MustU` and `MustL` (first order sets).
%                                Default = empty.
%    solverName:                 (string) Name of the solver used in
%                                GAMS. Default = 'cplex'.
%    runID:                      (string) ID for identifying this run.
%                                Default = ['run' date hour].
%    outputFolder:               (string) name for folder in which
%                                results will be stored. Default =
%                                'OutputsFindMustUU'.
%    outputFileName:             (string) name for files in which
%                                results. will be stored Default =
%                                'MustUUSet'.
%    printExcel:                 (double) boolean to describe wheter
%                                data must be printed in an excel file or
%                                not. Default = 1
%    printText:                  (double) boolean to describe wheter
%                                data must be printed in an plaint text
%                                file or not. Default = 1
%    printReport:                (double) 1 to generate a report in a
%                                plain text file. 0 otherwise. Default = 1
%    keepInputs:                 (double) 1 to mantain folder with
%                                inputs to run `findMustUU.gms`. 0 otherwise.
%                                Default = 1
%    keepGamsOutputs:            (double) 1 to mantain files returned by
%                                `findMustUU.gms`. 0 otherwise. Default = 1
%    verbose:                    (double) 1 to print results in console.
%                                0 otherwise. Default = 0
%
% OUTPUTS:
%    mustUU:                     (cell array of size number of sets found X
%                                2) Cell array containing the reactions IDs
%                                which belong to the `MustUU` set. Each row
%                                contain a couple of reactions that must
%                                decrease their flux.
%    pos_mustUU:                 (double array of size number of sets found
%                                X 2) double array containing the positions
%                                of each reaction in `mustUU` with regard to
%                                `model.rxns`
%    mustUU_linear:              (cell array of size number of unique
%                                reactions found X 1) Cell array containing
%                                the unique reactions ID which belong to
%                                the `MustUU` Set
%    pos_mustUU_linear:          (double array of size number of unique
%                                reactions found X 1) double array
%                                containing positions for reactions in
%                                `mustUU_linear`. with regard to `model.rxns`
%    outputFileName.xls:        (file) File containing one column
%                                array with identifiers for reactions in
%                                MustUU. This file will only be generated
%                                if the user entered `printExcel = 1`. Note
%                                that the user can choose the name of this
%                                file entering the input `outputFileName = 'PutYourOwnFileNameHere';`
%    outputFileName.txt:         (file) File containing one column
%                                array with identifiers for reactions in
%                                MustUU. This file will only be generated
%                                if the user entered `printText = 1`. Note
%                                that the user can choose the name of this
%                                file entering the input `outputFileName = 'PutYourOwnFileNameHere';`
%    outputFileName_Info.xls:    (file) File containing one column
%                                array. In each row the user will find a
%                                couple of reactions. Each couple of
%                                reaction was found in one iteration of
%                                `FindMustUU.gms`. This file will only be
%                                generated if the user entered `printExcel = 1`.
%                                Note that the user can choose the name
%                                of this file entering the input
%                                `outputFileName = 'PutYourOwnFileNameHere';`
%    outputFileName_Info.txt:    (file) File containing one column
%                                array. In each row the user will find a
%                                couple of reactions. Each couple of
%                                reaction was found in one iteration of
%                                `FindMustUU.gms`. This file will only be
%                                generated if the user entered `printText = 1`.
%                                Note that the user can choose the name
%                                of this file entering the input
%                                `outputFileName = 'PutYourOwnFileNameHere';`
%    findMustUU.lst:             (file) file autogenerated by GAMS. It
%                                contains information about equations,
%                                variables, parameters as well as
%                                information about the running (values at
%                                each iteration). This file only will be
%                                saved in the output folder is the user
%                                entered `keepGamsOutputs = 1`
%    GtoMUU.gdx:                 (file) file containing values for
%                                variables, parameters, etc. which were
%                                found by GAMS when solving `findMustUU.gms`.
%                                This file only will be saved in the output
%                                folder is the user entered `keepInputs = 1`
%
% NOTE:
%
%    This function is based in the GAMS files written by Sridhar
%    Ranganathan which were provided by the research group of Costas D.
%    Maranas. For a detailed description of the `optForce` procedure, please
%    see: `Ranganathan S, Suthers PF, Maranas CD (2010) OptForce: An
%    Optimization Procedure for Identifying All Genetic Manipulations
%    Leading to Targeted Overproductions. PLOS Computational Biology 6(4):
%    e1000744`. https://doi.org/10.1371/journal.pcbi.1000744
%
% .. Author: - Sebastian Mendoza, May 30th 2017, Center for Mathematical Modeling, University of Chile, snmendoz@uc.cl

optionalParameters = {'constrOpt', 'excludedRxns', 'mustSetFirstOrder', 'solverName', 'runID', 'outputFolder', 'outputFileName',  ...
    'printExcel', 'printText', 'printReport', 'keepInputs', 'keepGamsOutputs', 'verbose'};

if (numel(varargin) > 0 && (~ischar(varargin{1}) || ~any(ismember(varargin{1},optionalParameters))))

    tempargin = cell(1,2*(numel(varargin)));
    for i = 1:numel(varargin)

        tempargin{2*(i-1)+1} = optionalParameters{i};
        tempargin{2*(i-1)+2} = varargin{i};
    end
    varargin = tempargin;

end

parser = inputParser();
parser.addRequired('model', @(x) isstruct(x) && isfield(x, 'S') && isfield(model, 'rxns')...
    && isfield(model, 'mets') && isfield(model, 'lb') && isfield(model, 'ub') && isfield(model, 'b')...
    && isfield(model, 'c'))
parser.addRequired('minFluxesW', @isnumeric)
parser.addRequired('maxFluxesW', @isnumeric)
parser.addParamValue('constrOpt', struct('rxnList', {{}}, 'values', []),@ (x) isstruct(x) && isfield(x, 'rxnList') && isfield(x, 'values') ...
    && length(x.rxnList) == length(x.values) && length(intersect(x.rxnList, model.rxns)) == length(x.rxnList))
parser.addParamValue('excludedRxns', {}, @(x) iscell(x) && length(intersect(x, model.rxns)) == length(x))
parser.addParamValue('mustSetFirstOrder', {}, @(x) iscell(x) && length(intersect(x, model.rxns)) == length(x))
solvers = checkGAMSSolvers('MIP');
if isempty(solvers)
    error('there is no GAMS solvers available to solve Mixed Integer Programming problems') ;
else
    if ismember('cplex', lower(solvers))
        defaultSolverName = 'cplex';
    else
        defaultSolverName = lower(solvers(1));
    end
end

parser.addParamValue('solverName', defaultSolverName, @(x) ischar(x))
hour = clock; defaultRunID = ['run-' date '-' num2str(hour(4)) 'h' '-' num2str(hour(5)) 'm'];
parser.addParamValue('runID', defaultRunID, @(x) ischar(x))
parser.addParamValue('outputFolder', 'OutputsFindMustUU', @(x) ischar(x))
parser.addParamValue('outputFileName', 'MustUUSet', @(x) ischar(x))
parser.addParamValue('printExcel', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('printText', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('printReport', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('keepInputs', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('keepGamsOutputs', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('verbose', 1, @(x) isnumeric(x) || islogical(x));

parser.parse(model, minFluxesW, maxFluxesW, varargin{:})
model = parser.Results.model;
minFluxesW = parser.Results.minFluxesW;
maxFluxesW = parser.Results.maxFluxesW;
constrOpt= parser.Results.constrOpt;
excludedRxns= parser.Results.excludedRxns;
mustSetFirstOrder = parser.Results.mustSetFirstOrder;
solverName = parser.Results.solverName;
runID = parser.Results.runID;
outputFolder = parser.Results.outputFolder;
outputFileName = parser.Results.outputFileName;
printExcel = parser.Results.printExcel;
printText = parser.Results.printText;
printReport = parser.Results.printReport;
keepInputs = parser.Results.keepInputs;
keepGamsOutputs = parser.Results.keepGamsOutputs;
verbose = parser.Results.verbose;

% correct size of constrOpt
if ~isempty(constrOpt.rxnList)
    if size(constrOpt.rxnList, 1) > size(constrOpt.rxnList,2); constrOpt.rxnList = constrOpt.rxnList'; end;
    if size(constrOpt.values, 1) > size(constrOpt.values,2); constrOpt.values = constrOpt.values'; end;
end

% first, verify that GAMS is installed in your system
gamsPath = which('gams');
if isempty(gamsPath); error('OptForce: GAMS is not installed in your system. Please install GAMS.'); end;

%name of the function to solve the optimization problem in GAMS
gamsMustUUFunction = 'findMustUU.gms';
%path of that function
pathGamsFunction = which(gamsMustUUFunction);
if isempty(pathGamsFunction); error(['OptForce: ' gamsMustUUFunction ' not in MATLAB path.']); end;
%current path
workingPath = pwd;
%go to the path associate to the ID for this run.
if ~isdir(runID); mkdir(runID); end; cd(runID);

% if the user wants to generate a report.
if printReport
    %create name for file.
    hour = clock;
    reportFileName = ['report-' date '-' num2str(hour(4)) 'h' '-' num2str(hour(5)) 'm.txt'];
    freport = fopen(reportFileName, 'w');
    reportClosed = 0;
    % print date of running.
    fprintf(freport, ['findMustUUWithGAMS.m executed on ' date ' at ' num2str(hour(4)) ':' num2str(hour(5)) '\n\n']);
    % print matlab version.
    fprintf(freport, ['MATLAB: Release R' version('-release') '\n']);
    % print gams version.
    fprintf(freport, ['GAMS: ' regexprep(gamsPath, '\\', '\\\') '\n']);
    % print solver used in GAMS to solve optForce.
    fprintf(freport, ['GAMS solver: ' solverName '\n']);

    %print each of the inputs used in this running.
    fprintf(freport, '\nThe following inputs were used to run OptForce: \n');
    fprintf(freport, '\n------INPUTS------\n');
    %print model.
    fprintf(freport, '\nModel:\n');
    for i = 1:length(model.rxns)
        rxn = printRxnFormula(model, model.rxns{i}, false);
        fprintf(freport, [model.rxns{i} ': ' rxn{1} '\n']);
    end
    %print lower and upper bounds, minimum and maximum values for each of
    %the reactions in wild-type and mutant strain
    fprintf(freport, '\nLB\tUB\tMin_WT\tMax_WT\n');
    for i = 1:length(model.rxns)
        fprintf(freport, '%6.4f\t%6.4f\t%6.4f\t%6.4f\n', model.lb(i), model.ub(i), minFluxesW(i), maxFluxesW(i));
    end

    %print constraints
    fprintf(freport,'\nConstrained reactions:\n');
    for i = 1:length(constrOpt.rxnList)
        fprintf(freport,'%s: fixed in %6.4f\n',constrOpt.rxnList{i},constrOpt.values(i));
    end

    fprintf(freport, '\nExcluded Reactions:\n');
    for i = 1:length(excludedRxns)
        rxn = printRxnFormula(model, excludedRxns{i}, false);
        fprintf(freport, [excludedRxns{i} ': ' rxn{1} '\n']);
    end

    fprintf(freport, '\nReactions from first order sets(MustU and MustL):\n');
    for i = 1:length(mustSetFirstOrder)
        rxn = printRxnFormula(model, mustSetFirstOrder{i}, false);
        fprintf(freport, [mustSetFirstOrder{i} ': ' rxn{1} '\n']);
    end

    fprintf(freport,'\nrunID(Main Folder): %s \n\noutputFolder: %s \n\noutputFileName: %s \n',...
        runID, outputFolder, outputFileName);

    fprintf(freport,'\nprintExcel: %1.0f \n\nprintText: %1.0f \n\nprintReport: %1.0f \n\nkeepInputs: %1.0f  \n\nkeepGamsOutputs: %1.0f \n\nverbose: %1.0f \n',...
        printExcel,printText,printReport,keepInputs,keepGamsOutputs,verbose);

end

copyfile(pathGamsFunction);

% export inputs for running the optimization problem in GAMS to find the
% MustUU Set
inputFolder = 'InputsMustUU';
exportInputsMustOrder2ToGAMS(model, 'UU', minFluxesW, maxFluxesW, constrOpt, excludedRxns, mustSetFirstOrder, inputFolder)

% create a directory to save results if this don't exist
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder);
end

%run
if verbose
    run = system(['gams ' gamsMustUUFunction ' lo=3 --myroot=' inputFolder '/ --solverName=' solverName ' gdx=GtoMUU']);
else
    run=system(['gams ' gamsMustUUFunction ' --myroot=' inputFolder '/ --solverName=' solverName ' gdx=GtoMUU']);
end

if printReport; fprintf(freport, '\n------RESULTS------\n'); end;

%if user decide not to show inputs files for findMustUU.gms
if ~keepInputs; rmdir(inputFolder,'s'); end;

%if findMustUU.gms was executed correctly "run" should be 0
if run == 0

    if printReport; fprintf(freport,'\nGAMS was executed correctly\n'); end;
    if verbose; fprintf('GAMS was executed correctly\nSummary of information exported by GAMS:\n'); end;
    %show GAMS report in MATLAB console
    if verbose; gdxWhos GtoMUU; end;
    try
        findMustUU.name = 'findMustUU';
        rgdx('GtoMUU',findMustUU); %if do not exist the variable findMustUU in GtoMUU, an error will ocurr.
        if printReport; fprintf(freport,'\nGAMS variables were read by MATLAB correctly\n'); end;
        if verbose; fprintf('GAMS variables were read by MATLAB correctly\n'); end;

        %Using GDXMRW to read solutions found by findMustUU.gms
        %extract matrix 1 found by findMustII.gms. This matrix contains the
        %first reaction in each couple of reactions
        m1.name = 'matrix1';
        m1.compress = 'true';
        m1 = rgdx('GtoMUU',m1);
        uels_m1 = m1.uels{2};

        if ~isempty(uels_m1)
            %if the uel array for m1 is not empty, at least 1 couple of reations was found.
            if printReport; fprintf(freport,'\na MustUU set was found\n'); end;
            if verbose; fprintf('a MustUU set was found\n'); end;

            %find values for matrix 1
            val_m1 = m1.val;
            m1_full = full(sparse(val_m1(:,1),val_m1(:,2:end-1),val_m1(:,3)));

            %find values for matrix 2
            m2.name = 'matrix2';
            m2.compress = 'true';
            m2 = rgdx('GtoMUU',m2);
            uels_m2 = m2.uels{2};
            val_m2 = m2.val;
            m2_full = full(sparse(val_m2(:,1),val_m2(:,2:end-1),val_m2(:,3)));

            %initialize empty array for storing
            n_mustSet = size(m1_full,1);
            mustUU = cell(n_mustSet,2);
            pos_mustUU = zeros(size(mustUU));
            mustUU_linear = {};

            %write each couple of reactions.
            for i = 1:n_mustSet
                rxn1 = uels_m1(m1_full(i,:) == 1);
                rxn2 = uels_m2(m2_full(i,:) == 1);
                mustUU(i,1) = rxn1;
                mustUU(i,2) = rxn2;
                pos_mustUU(i,1) = find(strcmp(model.rxns,rxn1));
                pos_mustUU(i,2) = find(strcmp(model.rxns,rxn2));
                mustUU_linear = union(mustUU_linear,[rxn1;rxn2]);
            end
            pos_mustUU_linear = cell2mat(arrayfun(@(x)find(strcmp(x,model.rxns)),mustUU_linear,'UniformOutput', false))';
        else
            %if the uel array for m1 is empty, no couple of reations was found.
            if printReport; fprintf(freport,'\na MustUU set was not found\n'); end;
            if verbose; fprintf('a MustUU set was not found\n'); end;

            %initialize arrays to be returned by this function
            mustUU = {};
            pos_mustUU = [];
            mustUU_linear = {};
            pos_mustUU_linear = [];
        end

        % print info into an excel file if required by the user
        if printExcel
            if  ~isempty(uels_m1)
                currentFolder = pwd;
                cd(outputFolder);
                must = cell(size(mustUU,1),1);
                for i = 1:size(mustUU,1)
                    must{i} = strjoin(mustUU(i,:),' or ');
                end
                xlswrite([outputFileName '_Info'],[{'Reactions'};must]);
                xlswrite(outputFileName,mustUU_linear);
                cd(currentFolder);
                if verbose
                    fprintf(['MustUU set was printed in ' outputFileName '.xls  \n']);
                    fprintf(['MustUU set was also printed in ' outputFileName '_Info.xls  \n']);
                end
                if printReport
                    fprintf(freport,['\nMustUU set was printed in ' outputFileName '.xls  \n']);
                    fprintf(freport,['\nMustUU set was printed in ' outputFileName '_Info.xls  \n']);
                end
            else
                if verbose; fprintf('No mustUU set was found. Therefore, no excel file was generated\n'); end;
                if printReport; fprintf(freport, '\nNo mustUU set was found. Therefore, no excel file was generated\n'); end;
            end
        end

        % print info into a plain text file if required by the user
        if printText
            if ~isempty(uels_m1)
                currentFolder = pwd;
                cd(outputFolder);
                f = fopen([outputFileName '_Info.txt'],'w');
                fprintf(f,'Reactions\n');
                for i = 1:size(mustUU,1)
                    fprintf(f,'%s or %s\n',mustUU{i,1},mustUU{i,2});
                end
                fclose(f);

                f = fopen([outputFileName '.txt'],'w');
                for i = 1:length(mustUU_linear)
                    fprintf(f,'%s\n',mustUU_linear{i});
                end
                fclose(f);

                cd(currentFolder);
                if verbose
                    fprintf(['MustUU set was printed in ' outputFileName '.txt  \n']);
                    fprintf(['MustUU set was also printed in ' outputFileName '_Info.txt  \n']);
                end
                if printReport
                    fprintf(freport,['\nMustUU set was printed in ' outputFileName '.txt  \n']);
                    fprintf(freport,['\nMustUU set was printed in ' outputFileName '_Info.txt  \n']);
                end

            else
                if verbose; fprintf('No mustUU set was found. Therefore, no plain text file was generated\n'); end;
                if printReport; fprintf(freport, '\nNo mustUU set was found. Therefore, no plain text file was generated\n'); end;
            end
        end

        %close file for saving report
        if printReport; fclose(freport); reportClosed = 1; end;
        if printReport; movefile(reportFileName,outputFolder); end;
        delete(gamsMustUUFunction);

        %remove or move additional files that were generated during running
        if keepGamsOutputs
            if ~isdir(outputFolder); mkdir(outputFolder); end;
            movefile('GtoMUU.gdx',outputFolder);
            movefile(regexprep(gamsMustUUFunction,'gms','lst'),outputFolder);
        else
            delete('GtoMUU.gdx');
            delete(regexprep(gamsMustUUFunction,'gms','lst'));
        end

        %go back to the original path
        cd(workingPath);
    catch
        %GAMS variables were not read correctly by MATLAB
        if verbose; fprintf('GAMS variables were not read by MATLAB corretly\n'); end;
        if printReport && ~reportClosed; fprintf(freport, '\nGAMS variables were not read by MATLAB corretly\n'); fclose(freport); end;
        cd(workingPath);
        error('OptForce: GAMS variables were not read by MATLAB corretly');

    end

    %if findMustUU.gms was not executed correctly "run" should be different from 0
else
    %if GAMS was not executed correcttly
    if printReport && ~reportClosed; fprintf(freport, '\nGAMS was not executed correctly\n'); fclose(freport); end;
    if verbose; fprintf('GAMS was not executed correctly\n'); end;
    cd(workingPath);
    error('OptForce: GAMS was not executed correctly');

end

end

author: opencobra

repo: cobratoolbox


text:
function [fusion,params] = linear_fuser(w,scores)
%
%  Does affine fusion of scores: It does a weighted sum of scores and adds
%                                an offset.
%
%  Inputs:
%    scores: M-by-N matrix of N scores for each of M input systems.
%    w: Optional:
%         - when supplied, the output 'fusion' is the vector of fused scores.
%         - when w=[], the output 'fusion' is a function handle, to be used
%                      for training the fuser.
%       w is a (K+1)-vector, with one weight per system, followed by the
%       offset.
%
%    fusion: if w is given, fusion is a vector of N fused scores.
%            if w is not given, fusion is a function handle, so that
%              fusion(w) = @(w) linear_fusion(scores,w).
%    w0: default values for w, to initialize training.
%
%  For training use:
%     [fuser,params] = linear_fuser(train_scores);
%     w0 = get_w0();
%     w = train_binary_classifier(fuser,...,w0,...);
%
%  For test use:
%     fused_scores = linear_fuser(test_scores,w);
%

if nargin==0
    test_this();
    return;
end

if ~exist('scores','var') || isempty(scores)
    fusion = sprintf(['linear fuser:',repmat(' %g',1,length(w))],w);
    return;
end

wsz = size(scores,1)+1;
[whead,wtail] = splitvec_fh(wsz,w);
params.get_w0 = @() zeros(wsz,1);
%params.get_w0 = @() randn(wsz,1);
params.tail = wtail;

fusion = fusion_mv2df(whead,scores);

end

function test_this()

N = 100;
dim = 2;  % number of used systems

% ----------------synthesize training data -------------------
randn('state',0);
means = randn(dim,2)*8; %signal
[tar,non] = make_data(N,means);

% ------------- create system ------------------------------

[fuser,params] = linear_fuser([],[tar,non]);

% ------------- train it ------------------------------

ntar = size(tar,2);
nnon = size(non,2);
classf = [ones(1,ntar),-ones(1,nnon)];

prior = 0.1;
maxiters = 50;
quiet = true;
objfun = [];
w0 = params.get_w0();
[w,cxe] = train_binary_classifier(fuser,classf,w0,objfun,prior,[],0,maxiters,[],[],quiet);
fprintf('train Cxe = %g\n',cxe);

% ------------- test it ------------------------------

[tar,non] = make_data(N,means);

scores = [tar,non];
tail = [1;2;3];
wbig = [w;tail];
[fused_scores,params] = linear_fuser(wbig,scores);
check_tails = [tail,params.tail],
cxe = evaluate_objective(objfun,fused_scores,classf,prior);
fprintf('test Cxe = %g\n',cxe);

plot(fused_scores);

end

function [tar,non] = make_data(N,means)
[dim,K] = size(means);
X = 5*randn(dim,K*N); % noise
ii = 1:N;
for k=1:K
    X(:,ii) = bsxfun(@plus,means(:,k),X(:,ii));
    ii = ii+N;
end
N = K*N;
tar = X(:,1:N/2);
non = X(:,N/2+(1:N/2));

end

author: nesl


text:
function d = bayesgauss(X,C,M,P)
%BAYESGAUSS Bayes classifier for Gaussian patterns.
%   D = BAYESGAUSS(X,C,M,P) computes Bayes Gaussian decision functions
%   of the n-dimensional patterns in the rows of X. C is an array of
%   size n-by-n-by-Nc containing Nc covariance matrices of size n-by-n,
%   where Nc is the number of classes. M is a matrix of size Nc-by-n,
%   whose rows are the corresponding mean vectors. A covariance matrix
%   and a mean vector must be specified for each class. X is of size
%   K-by-n, where K is the number of patterns to be classified. P is a
%   1-by-Nc vector containing the probabilities of occurrence of each
%   class. If P is not included in the argument, the classes are assumed
%   to be equally likely.
%
%   In the output, D is a column vector of length K. Its ith element is
%   the class number assigned to the ith vector in X during
%   classification.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Verify the number of inputs.
narginchk(3,4)
n = size(C,1); % Dimension of patterns.

% Protect against the possibility that the class number is included
% as an (n + 1)th element of the vectors.
X = double(X(:,1:n));
% Number of pattern classes.
Nc = size(C,3);
% Number of patterns to classify.
K = size(X,1);
if nargin == 3
   P(1:Nc) = 1/Nc; % Classes assumed equally likely.
else
   if sum(P) ~= 1
      error('Elements of P must sum to 1.');
   end
end
% Compute the determinants.
DM = zeros(1,Nc); % Preallocate memory for loop speed.
for J = 1:Nc
   DM(J) = det(C(:,:,J));
end

% Evaluate the decision functions. Note the use of function mahalanobis
% discussed in Section 14.2.
D = zeros(K,Nc); % Preallocate memory for loop speed.
for J = 1:Nc
   Cm = C(:,:,J);
   Mm = M(J,:);
   L(1:K,1) = log(P(J));
   DET(1:K,1) = 0.5*log(DM(J));
   if P(J) == 0
      D(1:K,J) = -inf;
   else
      D(:,J) = L - DET - 0.5*mahalanobis(X,Cm,Mm);
   end
end

% Find the coordinates of the maximum value in each row. The location of
% a maximum gives the class of the corresponding pattern vector.
[i,j] = find(D == max(D,[],2));

% Re-use X. Its first element in a row is the pattern number, and the
% second is the class to which the pattern was assigned.
X = [i,j];

% Eliminate multiple classifications of the same patterns. Since the
% class assignment when two or more decision functions give the same
% value is arbitrary, we need to keep only one. Function unique
% eliminates duplicates and returns the array sorted by rows. This
% re-establishes the original order of the patterns.
[~,idx] = unique(X(:,1));
X = X(idx,:);
% X is now sorted, with the 2nd column giving the class of the pattern
% number in the 1st col.;  i.e., X(j,1) refers to the jth input pattern,
% and X(j,2) is its class number.

% Output the result of classification. d is a column vector with length
% equal to the total number of input patterns. The elements of d are the
% classes into which the patterns (in their original order) were
% classified.
d = X(:,2);

author: dipum

repo: dipum-toolbox

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


text:
function x=subint(x1,x2,x3,x4,na,nb,nc)
%subint   geometrically stretched subdivision generator
%   x=subint(x1,x2,x3,x4,na,nb,nc);
%   input
%             x1      left coordinate
%             x2      limit of uniform expansion section
%             x3      limit of equal subinterval section
%             x4      right coordinate and limit of contraction section
%             na      number of uniformly expanded subintervals
%             nb      number of intermediate uniform subintervals
%             nc      number of contracting subintervals
%    output
%              x      vector of coordinates
%
%   calls function fitint
%   sets up global variables global_N, global_INTL, global_LASTDL
%   IFISS function: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      global global_N global_INTL global_LASTDL
      intla=x2-x1;
      intlb=x3-x2;
      intlc=x4-x3;
      npts=na+nb+nc+1;
      bdl=intlb/nb;
      x=zeros(npts,1);
      x(1)=x1;
%
% uniform expansion section
      if na > 0
%% fixed value here
%        ratio = default('approximate expansion ratio (default 1.2)',1.2);
         ratio = 1.2;
         global_N=na;
         global_INTL=intla;
         global_LASTDL=bdl;
         [ratio,initdl]=fitint(ratio);
         fprintf('computed stretch ratio is %10.5g \n',ratio)
         dl=initdl;
         for node=2:na+1
            x(node)=x(node-1)+dl;
            dl=dl*ratio;
         end
         if (abs(x(na+1)-x2) > 1e-6)
            fprintf('\n warning ...\n calculated coordinate is %10.5g',x(na+1))
            fprintf('              \n      input coordinate is %10.5g',x2)
            x(na+1)=x2;
         end
      end
%
% uniform subinterval section
      dl=bdl;
      for node=na+2:na+nb+1
         x(node)=x(node-1) + dl;
      end
%
% uniform contraction section
      if nc > 0
%% fixed value here
%        ratio = default('approximate contraction ratio (default 1.2)',
         ratio = 1.2;
         global_N=nc;
         global_INTL=intlc;
         global_LASTDL=bdl;
         [ratio,initdl]=fitint(ratio);
%        fprintf('computed contraction ratio is %10.5g \n',ratio)
         dl=global_LASTDL/ratio;
         for node=na+nb+2:npts
            x(node)=x(node-1)+dl;
            dl=dl/ratio;
         end
         if (abs(x(npts)-x4) > 1e-6)
            fprintf('\n warning ...\n calculated coordinate is %10.5g',x(npts))
            fprintf('              \n      input coordinate is %10.5g',x4)
         end
      end
      x(npts)=x4;
      return