text:
function [W,H] = nmf(X,K,alg,maxiter,speak)
%
% NMF wrapper function
% function [W,H] = nmf(X,K,alg[,maxiter,speak])
%
% INPUT:
%           'X'     Inputmatrix
%           'K'     Number of components
%           'alg'   Algorithm to use:
%                   'mm'     multiplicative updates using euclidean
%                            distance. Lee, D..D., and Seung, H.S., (2001)
%                   'cjlin'  alternative non-negative least squares using
%                            projected gradients, author: Chih-Jen Lin,
%                            National Taiwan University.
%                   'prob'   probabilistic NFM interpretating X as samples
%                            from a multinomial, author: Lars Kai Hansen,
%                            Technical University of Denmark
%                   'als'    Alternating Least Squares. Set negative
%                            elements to zero.
%                   'alsobs' Alternating Least Squares. Set negative elements
%                            to zero and adjusts the other elements acording
%                            to Optimal Brain Surgeon.
%           'maxiter'   Maximum number of iterations, default = 1000.
%           'speak'     Print information to screen unless speak = 0,
%                       default = 0
%
% OUTPUT:
% W       : N x K matrix
% H       : K x M matrix
%
% Kasper Winther Joergensen
% Informatics and Mathematical Modelling
% Technical University of Denmark
% kwj@imm.dtu.dk
% 2006/12/15

switch(nargin)
    case {0,1,2}
        error('Missing parameter. Type "help nmf" for usage.');
        return
    case 3
        maxiter = 1000;
        speak = 0;
    case 4
        speak = 0;
    case 5
        % empty
    otherwise
        error('Too many parameters. Type "help nmf" for usage.');
        return
end

% find dimensionallity of X
[D,N] = size(X);

% switch algorithm
switch lower(alg)
    case 'mm'
        if speak, disp('Using mm algorithm'),end
        [W,H]=nmf_mm(X,K,maxiter,speak);
    case 'prob'
        if speak, disp('Using prob algorithm'),end
        [W,H]=nmf_prob(X,K,maxiter,speak);
    case 'cjlin'
        if speak, disp('Using cjlin algorithm'),end
        [W,H]=nmf_cjlin(X,rand(D,K),rand(K,N),0.000001,10000,maxiter);
    case 'als'
        if speak, disp('Using als algorithm'),end
        [W,H]=nmf_als(X,K,maxiter,speak);
    case 'alsobs'
        if speak, disp('Using alsobs algorithm'),end
        [W,H]=nmf_alsobs(X,K,maxiter,speak);
    otherwise
        error('Unknown method. Type "help nmf" for usage.');
        return
end

[W,H,nrgy] = order_comp(W,H);

author: andrewssobral

repo: lrslibrary


text:
function [G]=TabulateGenGam2dft(Rksi,Rgam,nu);
% function [G]=TabulateGenGam1dft(Rksi,Rgam,nu,K);
% This function tabulates the gain functions G for MMSE estimation
% of speech DFT coefficients, under the assumption of a
% Generalized Gamma speech amplitude prior with gamma=2. For explanations,
% derivations and motivations, see:
% J.S. Erkelens, R.C. Hendriks, R. Heusdens, and J. Jensen,
% "Minimum mean-square error estimation of discrete Fourier coefficients
% with Generalized Gamma priors", IEEE Trans. Audio, Speech and Language
% Processing, August 2007.
% and
%J.S. Erkelens, R.C. Hendriks and  R. Heusdens
%"On the Estimation of Complex Speech DFT Coefficients without Assuming Independent Real and Imaginary Parts",
%IEEE signal processing letters 2008.
%
% INPUT variables:
% Rksi: Array of "a priori" SNR (SNRprior) values for which the gains are
% computed. NOTE: The values must be in dBs.
% Rgam: Array of "a posteriori" SNR (SNRpost) values for which the gains
% are computed. NOTE: The values must be in dBs.
% 'nu': parameter of the Generalized Gamma distribution.
%
% OUTPUT variables:
% G: Matrix with gain values for speech DFT estimation,
% evaluated at all combinations of a priori and a posteriori SNR in the
% input variables Rksi and Rgam. To be multiplied with the noisy complex DFT coefficient.
%
%
% Copyright 2007: Delft University of Technology, Information and
% Communication Theory Group. The software is free for non-commercial use.
% This program comes WITHOUT ANY WARRANTY.
%
% Last modified: 22-11-2007.

Rgam=10.^(Rgam(:)'/10);%Rksi and Rgam are in dBs.
Rksi=10.^(Rksi/10);
G=zeros(length(Rksi),length(Rgam));
for k=1:length(Rgam)
    [g]=gam2gains(nu,Rgam(k),Rksi);
    G(:,k)=g(:);
end

function [G]=gam2gains(nu,gamm,ksi);

Q=ksi./(nu+ksi);
Teller=ConflHyperGeomFun(nu+1,2,Q.*gamm);%evaluate the confluent hypergeometric functions
Noemer=ConflHyperGeomFun(nu,1,Q.*gamm);
G=(nu.*Q).*Teller./Noemer;
% asymptotic result for Q.*gamm>700 (A&S, 13.5.1)  For very high values of nu and a posteriori SNR a value lower than 700 might be more appropriate.

I=find(Q.*gamm>700);
G(I)=Q(I);
%limit G for problematic cases
I=find(~isfinite(G));
G(I)=1;
I=find(G>1e4);
G(I)=1e4;

author: Sable

repo: mcbench-benchmarks


text:
%DTC Verzakov Tree - Trainable Decision Tree Classifier
%
%   W = DTC(A, CRIT, CHISQSTOPVAL, PRUNE, T, NFEAT)
%
% INPUT
%   A       Training dataset.
%           Object weights and class priors:
%             It is possible to assign individual weights to dataset
%             objects by using 'weights' identifier:
%
%             A = SETIDENT(A, WEIGHTS, 'weights')
%
%             If weights are not defined they assumed to be equal to 1 for
%             all objects. The actual weights used by the training routine
%             are computed by using supplied weights along with the class
%             priors (if priors are not set then the apparent priors are used):
%
%             ACTUAL_WEIGHTS(I) = M*(WEIGHTS(I)/CW(C))*(PRIOR(C)/SUM(PRIOR))
%
%             where M is the total amount of (labelled) objects in A,
%             CW is the vector weights sums for each class, and C is
%             the class of the object I. The sum of all actual weights is M.
%
%           Feature types:
%             Features are treated diffrently based on their domain information.
%             If the feature domain is empty or is the interval/collection of intervals
%             then this feature is considered to be the continuous one and
%             branches are created by splitting at the threshold value.
%             Otherwise (feature domain specifies a set of values or set of
%             names), the feature is considered to be the nominal one and
%             branches corresponding to all present feature values are
%             created.
%
%           Unknown and 'non-applicable' values:
%           If the feature value is NaN then it is assumed that its value
%           is unknown. In such a situation the object with unknown value
%           is split into fractions and sent down all branches.
%           The concept of the 'non-applicable' feature value is different
%           from the concept of the unknown (missing) value.
%           ...
%           The 'non-applicable' values for the continues features are
%           encoded as INF. If feature domain is the (set of) interval(s),
%           then INF value has to be explicitly added to the domain
%           defintion. The 'non-applicable' value of nominal features
%           does not have predefined encoding. If it is necessary user
%           have to include such value into domain definition.
%
%   CRIT   Splitting citerion name.
%          'igr' Information Gain Ratio (default):
%          As defined by Quinlan. The penalty on the number of the distinct
%          values of the continues feature is used. If the gain is zero or
%          negative due to such penalization, the split is not performed.
%          This leads to smaller trees and may give non-zero training error.
%          This criterion does not use costs. (Costs are used only at the classification step).
%
%          'gini' Gini impurity index. More precisely, the change in this
%          index. GINI index can be interpreted as a misclassification rate
%          for the stochastic prior based classifier, so costs are
%          naturally embedded. If the change in the (absolute) error less
%          or equal to 0.1 (change in the cost less or equal to 0.1 of minimal
%          absolute value of non-zero costs) the split is not performed.
%          This leads to smaller trees and may give non-zero training error.
%
%          'miscls' Classification error criterion.
%          To be used only for educational puposes because
%          it gives rather inferior results. Costs are naturally embedded.
%
%   CHISQSTOPVAL Early stopping crtitical value for the chi-squared test
%          on the difference between the original node class distribution and
%          branches class distributions. CHISQSTOPVAL is 0 by default.
%          Which means that branches will be discarded only if they bring no
%          change in class distribution.
%
%   PRUNE  Pruning type name
%          'prunep' - pessimistic (top-down) pruning as defined by Quinlan.
%          Pessimistic pruning be perfromed if cost matrix is defined.
%          'prunet' - test set (bottom-up) pruning using the (required)
%          dataset T.
%          These implementations of both pruning algorithms may be not
%          exactly correct if there are unknown values in datastes.
%
%   T      Test test for the test set pruning
%
% OUTPUT
%   W      Classifier mapping
%
% DESCRIPTION
%    If (full grown) branches of (sub)tree do not improve classification error
%    (misclassification cost) they are immediatley discarded.
%    This may happen because we use regularized posteriors. As a result
%    the algorithm is more stable, trees are smaller, but split on
%    the training set may be not perfect.
%
% REFERENCES
% [1] J.R. Quinlan, Simplifying Decision Trees, International Journal of
%     Man-Machine Studies, 27(3), pp. 221-234, 1987.
% [2] J.R. Quinlan, Improved use of continuous attributes in C4.5. Journal
%     of AI Research, 4(96), pp. 77-90, 1996.
%
% see also DATASETS, MAPPINGS, TREEC

% Copyright: S. Verzakov, s.verzakov@gmail.com
% Based on prtools' treec, tree_map, and their subroutines
% by Guido te Brake and R.P.W Duin

function w = dtc(a, crit, chisqstopval, prune, t, nfeat)

		% When no input data is given, an empty tree is defined:
	if (nargin == 0) || isempty(a)

    if nargin < 2 || isempty(crit)
			crit = 'igr';
    end

    if nargin < 3 || isempty(chisqstopval)
      chisqstopval = 0;
    end

    if nargin < 4 || isempty(prune)
      prune = '';
    end

    if nargin < 5
      t = [];
    end

    if nargin < 6
      nfeat = [];
    end

    w = prmapping('dtc', {crit, chisqstopval, prune, t, nfeat});
    w = setname(w, ['DecTree' upper(crit)]);

  elseif (nargin == 2) && ismapping(crit)
    % Execution
    w = crit;
    tree = +w;

    % B contains posteriors.
    % We do not need to convert posteriors to costs.
    % PRTools will do it automatically

    b = mapdt(tree, +a);
    b = setdat(a,b,w);
%    b = setfeatlab(b, getlabels(w));

    w = b;

  else
    % Training

    if nargin < 2 || isempty(crit)
      crit = 'igr';
    end

    if nargin < 3 || isempty(chisqstopval)
      chisqstopval = 0;
    end

    if nargin < 4 || isempty(prune)
      prune = '';
    end

    if nargin < 5
      t = [];
    end

    if nargin < 6
      nfeat = [];
    end

    if ~any(strcmpi(crit, {'igr', 'gini', 'miscls'}))
      error('Unknown splitting criterion');
    end

    islabtype(a,'crisp');
    isvaldfile(a, 1, 2); % at least 1 object per class, 2 classes
    a = seldat(prdataset(a)); % get rid of all unlabelled objects

    m = size(a,1);
    % Get weights (if defined)
    weights = getident(a, 'weights');
    if isempty(weights)
      weights = ones([m 1]);
    else
     idx = (weights > 0);
     if nnz(idx) < m
       a = a(idx, :);
       weights = weights(idx);
     end
     isvaldset(a, 1, 2);
    end

    % First get some useful parameters:

    % Sizes
    [m, k, c] = getsize(a);
    cs = classsizes(a);

    % Determine if features are categorical
    featdom = getfeatdom(a);
    featdom = featdom(:)';
    if isempty(featdom)
      feattype = zeros(1, k);
    else
      feattype = cellfun(@(x) ischar(x) || (size(x,1)==1), featdom);
    end

    % Features to use
    usefeat = true([1, k]);
    if ~isempty(nfeat)
      nfeat = min(nfeat, k);
    end

    % Labelling
    nlab = getnlab(a);

    % Get priors
    prior = getprior(a);
    prior = prior/sum(prior);

    % Define weights:
    % the total sum of weights is equal to m
    % Compute class weights
    cw = zeros(size(cs));
    for i=1:c
      cw(i) = sum(weights(nlab == i));
    end
    % Rescale objects weights
    % by class weight factors derived from priors
    cwf = m*prior./cw;
    weights = weights.*(cwf(nlab)).';

    % Miscalssification cost
    cost = a.cost;

    % Now the training can really start:
    tree = makedt(+a, feattype, usefeat, nfeat, c, nlab, weights, cost, crit, chisqstopval);

    if ~isempty(prune)
      if strcmpi(prune, 'prunep')
        if ~isempty(cost)
          error('Pessimistic prunning based on misclassification costs is not implemented');
        end
        tree = prunep(tree);
      elseif strcmpi(prune, 'prunet')
        if isempty(t)
          error('Test set is not specified for the test set prunning');
        end
         tree = prunet(tree, t, cost);
      else
        error('Unkown prunning method');
      end

      tree = cleandt(tree);
    end

    % Store the results:
    w = prmapping('dtc', 'trained', {tree}, getlablist(a), k, c);
    w = setname(w, ['DecTree' upper(crit)]);
    w = setcost(w, cost);
  end

  return

%MAKEDT General tree building algorithm
%
%   TREE = MAKEDT(A, FEATTYPE, USEFEAT, NFEAT, C, NLAB, WEIGHTS, COST, CRIT, CHISQSTOPVAL)
%
% INPUT
%   A      Data matrix
%
%   FEATTYPE Row defining the feature types (0 - continues, 1 - nominal)
%
%   USEFEAT Row defining the features to be used for splitting
%
%   NFEAT  The maximum number of features for which splitting has to be
%          attempted. If the number of the available features is geater
%          than NFEAT the random subset of NFEAT features will be used for
%          splitting.
%
%  C       Total number of classes in the initial dataset
%
%  NLAB    Numeric class labels of objects in A
%
%  WEIGHTS Weights of objects in A
%
%  COST    Misclassification costs
%
%  CRIT    Name of splitting criterion to be used
%
%  CHISQSTOPVAL Crtitical value for the chi-squared test
%
% OUTPUT
%   TREE   Structure containing arrays defining decision tree.
%          .NSMP The number of samples. NSMP(J) is the sum of weights of objects
%          which reached the node J.
%
%          .POST Posterior probabilities. POST(J, :) is the class
%          distribution at the node J. Object weights are taken into
%          account. 0 and 1 probablities are avoided by perfroming Bayes
%          uniform priors regularization.
%
%          .CIDX Class index. CIDX(J) is the class corresponding to the
%          node J if it considered as a leaf.
%
%          .ERRL Leaf error. ERRL(J) is the misclassification error (cost)
%          on training samples which reached the node J if this node is
%          considered to be a leaf.
%
%          .ERRT Leaf error. ERRT(J) is the misclassification error (cost)
%          on training samples which reached the node J if this node is
%          considered to be a root of the (sub)tree.
%
%          .SIZE Tree size. SIZE(J) is the (sub)tree size with root at the
%          node J.
%
%          .FIDX Feature index. FIDX(J) is the index of the feature on
%          which node J is split. FIDX(J) == 0 means that J is a leaf.
%
%          .FVAL Feature value(s). For nominal features FVAL{J} is the set
%          of feature values observed at the node J (LENGTH(FVAL{J} == 0
%          for the leaf, otherwise it is > 1).
%          For continues features, if LENGTH(FVAL{J}) == 1 then it contains
%          threshold THR for splitting into the left (<= THR) and the
%          right (> THR) branches. If FVAL{J} == [] (and J is not a leaf)
%          then it means that split is perfromed between applicable and
%          non-applicable values.
%
%          .NIDX Branch node indices. For nominal features NIDX{J,K} is
%          the index of branch node of node J with value FVAL{J,K}
%          of feature FIDX(J). For continues features (if LENGTH(FVAL{J}) == 1)
%          NIDX{J,1} is the index of the left branch node, NIDX{J,2} is the
%          index of the right branch node. If FVAL{J} == []
%          (and J is not a leaf) then NIDX{J,1} is the index of branch node
%          containing objects with applicable values of feature FIDX(J) and
%          NIDX{J,2} is the index of branch node containing objects with
%          non-applicable values of the same feature.
%
% This is a low-level routine called by DTC.
%
% See also IGR, GINI, MISCLS

function tree = makedt(a, feattype, usefeat, nfeat, c, nlab, weights, cost, crit, chisqstopval)

  % Construct the tree:

  % Find (absolute) class frequencies
  C = zeros(1, c);
  for j=1:c
    C(j) = sum(weights(nlab == j));
  end

  NC = nnz(C);
  tree.nsmp = sum(C);

  % regularization by 'uniform' Bayesian priors;
  C0 = C;
  C = C + 1;
  p = C/sum(C);

  if isempty(cost)
    [maxpost, cidx] = max(p);
    errc = tree.nsmp*(1-maxpost);
    %errc = tree.nsmp - C0(cidx);
  else
    costp = p*cost;
    [mincost, cidx] = min(costp);
    errc = mincost;
  end

  tree.post = p;
  tree.cidx = cidx;
  tree.errl = errc;
  tree.errt = errc;
  tree.size = 1;

  if NC ~= 1 % not a pure class dataset
    % now the tree is recursively constructed further:
		% use desired split criterion
    [fidx, fval, nb, bidx, chisqval] = findsplit(+a, feattype, usefeat, nfeat, c, nlab, weights, cost, crit);

    % When the stop criterion is not reached yet, we recursively split
		% further:
    if ~isempty(fidx) && (chisqval > chisqstopval)
      tree.fidx = fidx;
      tree.fval = {fval};

      if feattype(fidx) > 0
        usefeat(fidx) = 0;
      end

      uidx = bidx == 0;
      nu = nnz(uidx);
      if nu > 0
        knsmp = tree.nsmp - sum(weights(uidx));
      end

      tree.nidx = {zeros(1, nb)};
      tree.errt(1) = 0;
      for j=1:nb
        tree.nidx{1}(j) = tree.size(1) + 1;

        J = bidx == j;
        if nu == 0
          bweights = weights(J);
        else
          J = J | uidx;
          bweights = weights(J);
          buidx = uidx(J);
          bknsmp = sum(bweights(~buidx));
          bweights(buidx) = (bknsmp/knsmp)*bweights(buidx);
        end

        branch = makedt(+a(J, :), feattype, usefeat, nfeat, c, nlab(J), bweights, cost, crit, chisqstopval);

        branch.nidx = cellfun(@(x) x + tree.size(1)*(x>0), branch.nidx, 'UniformOutput', false);
        tree.errt(1) = tree.errt(1) + branch.errt(1);
        tree.size(1) = tree.size(1) + size(branch.nidx, 1);

        tree.nsmp = [tree.nsmp; branch.nsmp];
        tree.post = [tree.post; branch.post];
        tree.cidx = [tree.cidx; branch.cidx];
        tree.errl = [tree.errl; branch.errl];
        tree.errt = [tree.errt; branch.errt];
        tree.size = [tree.size; branch.size];
        tree.fidx = [tree.fidx; branch.fidx];
        tree.fval = [tree.fval; branch.fval];
        tree.nidx = [tree.nidx; branch.nidx];
      end
    end
  end

  % no improvement in error (cost), rollback
  if (tree.size(1) > 1) && (tree.errt(1) >= tree.errl(1))
    tree.nsmp = tree.nsmp(1);
    tree.post = tree.post(1,:);
    tree.cidx = tree.cidx(1);
    tree.errl = tree.errl(1);
    tree.errt = tree.errl(1); % sic!
    tree.size = 1;
  end

  if tree.size(1) == 1
    % We reached the stop criterion or no further split is possible
    % so we make a leaf node:
    tree.fidx = 0;
    tree.fval = {[]};
    tree.nidx = {[]};
  end

	return

%MAPDT Tree mapping and node statistic calculation
%
% 	[P, S] = MAPDT(TREE, A, NLAB, WEIGHTS)
%

function [p, s] = mapdt(tree, a, nlab, weights)

  persistent dt st

  if isstruct(tree)
    dt = tree;

    if nargin < 4
      weights = [];
    end

    if nargin < 3
      nlab = [];
    end

    m = size(a, 1);
    [n, c] = size(dt.post);
    p = zeros([m, c]);

    if (nargout < 2) || isempty(nlab)
      st = [];

      for i=1:m
        p(i, :) = mapdt(1, +a(i, :));
      end

    else
      st = zeros([n, c]);
      nlab (nlab > c) = 0;
      if isempty(weights)
        weights = ones(m, 1);
      end

      for i=1:m
        p(i, :) = mapdt(1, +a(i, :), nlab(i), weights(i));
      end
    end

    s = st;
    clear dt st

  else
    j = tree;

    while j > 0
      if ~isempty(st) && (nlab > 0)
        st(j, nlab) = st(j, nlab) + weights;
      end

      k = 0;
      fidx = dt.fidx(j);

      if fidx ~= 0
        nidx = dt.nidx{j};
        aval = a(fidx);

        if ~isnan(aval)
          fval = dt.fval{j};

          if isempty(fval)
            k = nidx(2 - (aval ~= inf));
          elseif length(fval) == 1
            if aval ~= inf
              k = nidx(2 - (aval <= fval));
            elseif length(nidx) == 3
              k = nidx(3);
            end;
          else
            k = nidx(aval == fval);
            if isempty(k)
              k = 0;
            end
          end

        else
          p = zeros([1, size(dt.post, 2)]);
          if isempty(st) || (nlab <= 0)
            for b=1:length(nidx)
              k = nidx(b);
              f = (dt.nsmp(k)/dt.nsmp(j));
              p = p + f*mapdt(nidx(b), a);
            end
          else
            for b=1:length(nidx)
              k = nidx(b);
              f = (dt.nsmp(k)/dt.nsmp(j));
              p = p + f*mapdt(nidx(b), a, nlab, f*weights);
            end
          end
          k = -1;
        end
      end

      if k == 0
        p = dt.post(j, :);
      end

      j = k;
    end
  end

  return

%FINDSPLIT General routine for finding the best split
%
% 	[FIDX, FVAL, NB, BIDX, CHISQVAL, CRITVAL] = FINDSPLI(A, FEATTYPE, USEFEAT, NFEAT, C, NLAB, WEIGHTS, COST, CRIT)
%

function [fidx, fval, nb, bidx, chisqval, critval] = findsplit(a, feattype, usefeat, nfeat, c, nlab, weights, cost, crit)

  selfeatidx = find(usefeat);
  nf = length(selfeatidx);
  if ~isempty(nfeat) && (nfeat < nf)
    permidx = randperm(length(selfeatidx));
    selfeatidx = selfeatidx(permidx(1:nfeat));
    nf = nfeat;
  end

  fval = cell([1, nf]);
  chisqval = zeros([1, nf]);
  critval = nan([1, nf]);

  % repeat for all selected features
  for f=1:nf
    fidx = selfeatidx(f);
    af = a(:, fidx);

    kidx = ~isnan(af); % known values index
    if nnz(kidx) == 0
      continue
    end

    MU = sum(weights(~kidx)) + realmin;

    af = af(kidx);
    wk = weights(kidx);
    nlabk = nlab(kidx);

    switch feattype(fidx)
      case 0 % continous/ordered feature
        naidx = af == inf;
        NA = repmat(realmin, [1 c]);
        MNA = c*realmin;
        nlabna = nlabk(naidx);
        nna = length(nlabna);
        if nna > 0
          for j = 1:c
            NA(1,j) = sum(wk(nlabna == j)) + realmin;
          end
          MNA = sum(NA);
        end

        apidx = find(~naidx);
        af = af(apidx);
        wap = wk(apidx);
        nlabap = nlabk(apidx);
        [af, sortidx] = sort(af);
        wap = wap(sortidx);
        nlabap = nlabap(sortidx);

        if length(af) == 1
          uv = af;
          ns = 0;
          sli = [];
        else
          labchngcount = cumsum([1; double(diff(nlabap) ~= 0)]);
          uniquevallowidx = find([true; ((diff(af)./(0.5*abs(af(1:end-1) + af(2:end)) + realmin)) > 1e-8)]);
          uniquevalhighidx = [(uniquevallowidx(2:end) - 1); length(af)];
          % unique values
          uv = af(uniquevalhighidx);

          % split low indices in unique values
          sli = find(labchngcount(uniquevalhighidx(2:end)) - labchngcount(uniquevallowidx(1:end-1)) > 0);
          % split low indices in af
          splitlowidx = uniquevalhighidx(sli);
          ns = length(splitlowidx);
        end

        if (ns == 0) && (nna == 0)
          continue
        end

        % applicable, left, and right branch class counts
        AP = zeros(1, c);
        L = zeros(ns, c);
        R = zeros(ns, c);

        for j = 1:c
          J = find(nlabap == j);
          mj = length(J);
          AP(j) = sum(wap(J));
          if (ns > 0) && (mj > 0)
            L(:, j) = (repmat(splitlowidx, [1, mj]) >= repmat(J.', [ns, 1]))*wap(J) + realmin;
            R(:, j) = AP(j) - L(:, j) + realmin;
          end
        end
        AP = AP + 2*realmin;

        % total count of applicable
        MAP = sum(AP);

        % known
        K = AP + NA;
        MK = MNA + MAP;

        % object counts for branches
        ML = sum(L, 2);
        MR = sum(R, 2);

        [cv, i] = feval(crit, 0, MU, MK, K, MNA, NA, MAP, AP, ML, L, MR, R, cost, weights, uv, sli);
        critval(f) = cv;

        if ~isnan(cv) && (cv > -inf)
          if (ns > 0) && ~isempty(i)
            t = splitlowidx(i);
            fval{f} = 0.5*(af(t) + af(t+1));

            if nna == 0
              APL = AP*(ML(i)/MAP);
              APR = AP*(MR(i)/MAP);
              chisqval(f) = sum(((L(i, :) - APL).^2)./APL + ((R(i, :) - APR).^2)./APR);
            else
              KNA = K*(MNA/MK);
              KL = K*(ML(i)/MK);
              KR = K*(MR(i)/MK);
              chisqval(f) = sum(((NA - KNA).^2)./KNA + ((L(i, :) - KL).^2)./KL + ((R(i, :) - KR).^2)./KR);
            end
          else
            fval{f} = [];
            KNA = K*(MNA/MK);
            KAP = K*(MAP/MK);
            chisqval(f) = sum(((NA - KNA).^2)./KNA + ((AP(i, :) - KAP).^2)./KAP);
          end
        end

      case 1 % nominal feature
        vf = unique(af);
        v = length(vf);

        B = zeros(v, c);
        for j=1:c
          J = find(j == nlab);
          mj = length(J);
          if mj > 0
            B(:, j) = (repmat(vf, [1, mj]) ==  repmat(af(J).', [v 1]))*weights(J);
          end
        end

        B = B + realmin;

        % class counts for known
        K = sum(B, 1);

        % total known count
        MK = sum(K);

        % object counts for branches
        MB = sum(B, 2);

        cv = feval(crit, 1, MU, MK, K, MB, B, cost, weights);
        critval(f) = cv;

        if ~isnan(cv) && (cv > -inf)
          KB = MB.*K/MK;
          chisqval(f) = sum(sum(((B - KB).^2)./KB));
          fval{f} = bf;
        end
    end
  end

  % best criterion over all features
  testfeatidx = find(~isnan(critval) & (critval > -inf));
  if isempty(testfeatidx)
    fidx = [];
    fval = [];
    nb = [];
    bidx = [];
    chisqval = [];
    critval = [];

  else
    [critval, fidx] = max(critval(testfeatidx));
    fidx = testfeatidx(fidx);
    fval = fval{fidx};
    chisqval = chisqval(fidx);
    fidx = selfeatidx(fidx);

    af = a(:, fidx);
    m = size(af, 1);
    bidx = zeros(m, 1);

    switch feattype(fidx)
      case 0
        apidx = af < inf;

        if ~isempty(fval)
          lidx = af <= fval;
          ridx = apidx & ~lidx;
          bidx(lidx) = 1;
          bidx(ridx) = 2;
          nb = 2;
        else
          bidx(apidx) = 1;
          nb = 1;
        end

        naidx = ~isnan(af) & ~apidx;
        if nnz(naidx) > 0
          nb = nb + 1;
          bidx(naidx) = nb;
        end

      case 1
        nb = length(fval);
        for i=1:nb
          bidx(af == fval(i)) = i;
        end
    end
  end

  return

%IGR The information gain ratio
%
% 	[CRITVAL, IDX] = IGR(FEATTYPE, MU, MK, K, MNA, NA, PRMAP, AP, ML, L, MR, R, COST, WEIGHTS, UV, SLI)
%   [CRITVAL, IDX] = IGR(FEATTYPE,MU, MK, K, MB, B, COST, WEIGHTS

function [critval, idx] = igr(feattype, varargin)

  switch feattype
    case 0
      [MU, MK, K, MNA, NA, PRMAP, AP, ML, L, MR, R, cost, weights, uv, sli] = deal(varargin{:});

      M = MK + MU;

      infoold = - K * log2(K.'/MK);
      infonew = - NA * log2(NA.'/MNA);
      infosplit = - MU * log2(MU/M) - MNA * log2(MNA/M);

      if ~isempty(R)
        infolr = - ( ...
          sum(L .* log2(L./(repmat(ML, [1 size(L, 2)]))), 2) + ...
          sum(R .* log2(R./(repmat(MR, [1 size(R, 2)]))), 2) ...
        );

        % best criterion value over all thresholds
        [infolr, idx] = min(infolr);

        infonew = infonew + infolr;
        infosplit = infosplit - [ML(idx), MR(idx)] * log2([ML(idx); MR(idx)]/M);

      else
        idx = [];
        infonew = infonew - AP * log2(AP.'/MAP);
        infosplit = infosplit - MAP * log2(MAP/M);
      end

      infothresh = log2(max(length(uv) - 1, 1));
      infogain = (infoold - infonew - infothresh);

    case 1
      idx = [];
       [MU, MK, K, MB, B, cost, weights] = deal(varargin{:});

      M = MK + MU;

      infoold = - K * log2(K.'/MK);
      infosplit = - MU.' * log2(MU/M);

      infonew = - sum(sum(B .* log2(B./(repmat(MB, [1 size(B, 2)]))), 2), 1);
      infosplit = infosplit - MB.' * log2(MB/M);

      infogain = (infoold - infonew);
  end

  % infogain = infogain / M;
  % infosplit = infospit / M;

  if infogain > 0
    critval = infogain / infosplit;
  else
    critval = -inf;
  end

  return

%GINI
%
% 	[CRITVAL, IDX] = GINI(FEATTYPE, MU, MK, K, MNA, NA, PRMAP, AP, ML, L, MR, R, COST, WEIGHTS, UV, SLI)
%   [CRITVAL, IDX] = GINI(FEATTYPE,MU, MK, K, MB, B, COST, WEIGHTS%

function [critval, idx] = gini(feattype, varargin)

  switch feattype
    case 0
      [MU, MK, K, MNA, NA, PRMAP, AP, ML, L, MR, R, cost, weights, uv, sli] = deal(varargin{:});

      if isempty(cost)
        purityold = K * (K.'/MK);
        %impurityold = MK - purityold;

        if ~isempty(R)
          puritylr = ( ...
            sum(L .* (L./(repmat(ML, [1 size(L, 2)]))), 2) + ...
            sum(R .* (R./(repmat(MR, [1 size(R, 2)]))), 2) ...
          );

          % best criterion value over all thresholds
          [puritylr, idx] = max(puritylr);
          puritynew = puritylr + NA * (NA.'/MNA);

        else
          idx = [];
          puritynew = AP * (AP.'/ MAP) + NA * (NA.'/MNA);
        end

        %impuritynew = MK - puritynew;
        deltaimpurity = puritynew - purityold;

      else
        impurityold = K * cost * (K.'/MK);

        if ~isempty(R)
          impuritylr = ...
            sum((L*cost) .* (L./(repmat(ML, [1 size(L, 2)]))), 2) + ...
            sum((R*cost) .* (R./(repmat(MR, [1 size(R, 2)]))), 2);

          [impuritylr, idx] = min(impuritylr);
          impuritynew = impuritylr + NA * cost * (NA.'/MNA);

        else
          idx = [];
          impuritynew = AP * cost* (AP.'/ MAP) + NA * cost *(NA.'/MNA);
        end

        deltaimpurity = impurityold - impuritynew;
      end

    case 1
      idx = [];
      [MU, MK, K, MB, B, cost, weights] = deal(varargin{:});

      if isempty(cost)
        purityold = K * (K.'/MK);
        %impurityold = MK - purityold;
        puritynew = sum(sum(B .* (B./(repmat(MB, [1 size(B, 2)]))), 2));
        %impuritynew = MK - puritynew;
        deltaimpurity = puritynew - purityold;

      else
        impurityold = K * cost * (K.'/MK);
        impuritynew = ( ...
          sum(sum((B*cost) .* (B./(repmat(MB, [1 size(B, 2)]))), 2), 1) ...
        );

        deltaimpurity = impurityold - impuritynew;
      end
  end

  if isempty(cost)
    deltamin = 0.1;
  else
    deltamin = 0.1*min(abs(cost(abs(cost) > 0)));
  end

  if deltaimpurity > deltamin
    critval = deltaimpurity/(MU+MK);
  else
    critval = -inf;
  end

  return

%MISCLS Miscalssification error
%
% 	[CRITVAL, IDX] = MISCLS(FEATTYPE, MU, MK, K, MNA, NA, PRMAP, AP, ML, L, MR, R, COST, WEIGHTS, UV, SLI)
%   [CRITVAL, IDX] = MISCLS(FEATTYPE,MU, MK, K, MB, B, COST, WEIGHTS
%

function [critval, idx] = miscls(feattype, varargin)

  switch feattype
    case 0
      [MU, MK, K, MNA, NA, PRMAP, AP, ML, L, MR, R, cost, weights, uv, sli] = deal(varargin{:});

      if isempty(cost)
        ccold = max(K, [], 2);
        if ~isempty(R)
          ccnew = max(L, [], 2) + max(R, [], 2);
          [ccnew, idx] = max(ccnew);
          ccnew = ccnew + max(NA, [], 2);
        else
          idx = [];
          ccnew = max(AP, [], 2) + max(NA, [], 2);
        end

        deltamc = ccnew - ccold;

      else
        mcold = min(K*cost, [], 2);
        if ~isempty(R)
          mcnew = min(L*cost, [], 2) + min(R*cost, [], 2);
          [mcnew, idx] = min(mcnew);
           mcnew = mcnew + min(NA*cost, [], 2);
        else
          idx = [];
          mcnew = min(AP*cost, [], 2) + min(NA*cost, [], 2);
        end

        deltamc = mcold - mcnew;
      end

    case 1
      idx = [];
      [MU, MK, K, MB, B, cost, weights] = deal(varargin{:});

      if isempty(cost)
        ccold = max(K, [], 2);
        ccnew = sum(max(B, [], 2), 1);
        deltamc = ccnew - ccold;
      else
        mcold = min(K*cost, [], 2);
        mcnew  = sum(min(B*cost, [], 2), 1);
        deltamc = mcold - mcnew;
      end
  end

  %if isempty(cost)
  %  deltamin = min(weights(weights > 0));
  %else
  %  deltamin = min(weights(weights > 0))*min(cost(cost > 0));
  %end
  %
  %if ~isempty(deltamin)
  %  deltamin = 0.5*deltamin;
  %else
  %  deltamin = 0;
  %end
  %
  %if deltamc <= deltamin
  %  critval = -inf;
  %else
  %  critval = deltamc / (MU+MK);
  %end

  critval = deltamc / (MU+MK);

  return

%PRUNEP Pessimistic pruning of a decision tree
%
% 	TREE = PRUNEP(TREE,NODE)
%
% Pessimistic pruning defined by Quinlan.

function tree = prunep(tree, node)

  persistent pt;

  if (nargout ~= 0) || (nargin ~= 1) || ~isscalar(tree) || ~isnumeric(tree) || ~isint(tree)
    if nargin < 2 || isempty(node)
      node = 1;
    end

    pt = tree;
    prunep(node);
    tree = pt;
    clear pt;

  else
    node = tree;
    if pt.fidx(node) > 0
      tnidx = (node+1):(node + pt.size(node) - 1);

      nleaves = nnz(pt.fidx(tnidx) == 0);
      errt = pt.errt(node) + 0.5*nleaves;
      errl = pt.errl(node) + 0.5;
      nsmp = pt.nsmp(node);
      sd = sqrt(errt*(1-errt/nsmp));

      if errl < (errt + sd)
        pt.fidx(tnidx) = -1;

        pt.errt(node) = pt.errl(node);
        pt.size(node) = 1;
        pt.fidx(node) = 0;
        pt.fval(node) = {[]};
        pt.nidx(node) = {[]};

      else
        errt = 0;
        for i=1:length(pt.nidx{node})
          idx = pt.nidx{node}(i);
          prunep(idx);
          errt = errt + pt.errt(idx);
        end
        pt.errt(idx) = errt;
      end
    end
  end

  return

%PRUNET Prune tree by testset
%
% 	TREE = PRUNET(TREE,T,COST)
%
% The test set a is used to prune a decision tree.

function tree = prunet(tree, t, cost)

  persistent pt;

  if (nargout ~= 0) || (nargin ~= 1) || ~isscalar(tree) || ~isnumeric(tree) || ~isint(tree)
    if nargin < 3
      cost = [];
    end

    [m, k, c] = getsize(t);

    cs = classsizes(t);
    prior = getprior(t);
    prior = prior/sum(prior);
    weights = m*prior./cs;

    pt = tree;
    mt = size(pt.post,1);
    [p, s] = mapdt(pt, +t, getnlab(t), weights);

    % error (cost) in each node as if there were leafs
    if isempty(cost)
      idx = sub2ind([mt, c], (1:mt).', pt.cidx);
      pt.test_errl = sum(s,2) - s(idx);
    else
      pt.test_errl = sum(s .* cost(:, pt.cidx).', 2);
    end

    pt.test_errt = pt.test_errl;

    prunet(1);
    tree = pt;
    clear pt;

  else
    node = tree;
    if pt.fidx(node) > 0
      test_errt = 0;
      errt = 0;

      for i=1:length(pt.nidx{node})
        idx = pt.nidx{node}(i);
        prunet(idx);
        test_errt = test_errt + pt.test_errt(idx);
        errt = errt + pt.errt(idx);
      end

      if pt.test_errl(node) <= test_errt
        pt.fidx(pt.nidx{node}) = -1;

        pt.errt(node) = pt.errl(node);
        pt.size(node) = 1;
        pt.fidx(node) = 0;
        pt.fval(node) = {[]};
        pt.nidx(node) = {[]};

      else
        pt.test_errt(node) = test_errt;
        pt.errt(node) = errt;
      end
    end
  end;

  return

function tree = cleandt(tree)

  rnidx = tree.fidx == -1;
  if nnz(rnidx) == 0
    return
  end

  rncount = cumsum(rnidx);
  fn = fieldnames(tree);
  for i=1:length(fn)
    tree.(fn{i})(rnidx, :) = [];
  end

  tree.nidx = cellfun(@(x) x - rncount(x).', tree.nidx, 'UniformOutput', false);

  idx = (1:length(rnidx)).';
  idx(rnidx) = [];
  tree.size = tree.size - (rncount(idx) - rncount(idx + tree.size - 1));

 	return

author: marianux

repo: ecg-kit

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


text:
function t_ext2int2ext(quiet)
%T_EXT2INT2EXT  Tests EXT2INT, INT2EXT, and related functions.
%   Includes tests for GET_REORDER, SET_REORDER, E2I_DATA, I2E_DATA
%   E2I_FIELD, I2E_FIELD, EXT2INT and INT2EXT.

%   MATPOWER
%   Copyright (c) 2009-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

num_tests = 165;
t_begin(num_tests, quiet);

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;

%% for MATLAB versions prior to R2012a (v 7.14)
if ~exist('isequaln')
    eval('isequaln = @isequalwithequalnans;');
end

if have_feature('matlab', 'vnum') < 7.001
    t_skip(num_tests, 'test requires cellfun() construct not available before MATLAB 7.1');
else
    mpce = loadcase('t_case_ext');
    mpce = add_userfcn(mpce, 'ext2int', @userfcn_test_ext2int);
    mpce = add_userfcn(mpce, 'int2ext', @userfcn_test_int2ext);
    mpci = loadcase('t_case_int');

    gen_i2e = [1 2 4];

    An = mpce.xbus;
    As = mpce.strbus;

    %%-----  get_reorder  -----
    k = [3; 7; 4; 1];
    t = 'get_reorder(A, k, 1) : numeric';
    t_is(get_reorder(An, k, 1), An(k, :), 12, t);

    t = 'get_reorder(A, k, 2) : numeric';
    t_is(get_reorder(An, k, 2), An(:, k), 12, t);

    %%-----  set_reorder  -----
    k = (2:2:10)';
    t = 'set_reorder(A, B, k, 1) : numeric';
    B = An(k, :) * -1;
    got = set_reorder(An, B, k, 1);
    ex = An;
    ex(k, :) = -ex(k, :);
    t_is(got, ex, 12, t);

    t = 'set_reorder(A, B, k, 2) : numeric';
    B = An(:, k) * -1;
    got = set_reorder(An, B, k, 2);
    ex = An;
    ex(:, k) = -ex(:, k);
    t_is(got, ex, 12, t);

    t = 'set_reorder(A, Bshort, k, 1) : numeric';
    k = (2:2:8)';
    B = An(k, 1:8) * -1;
    got = set_reorder(An, B, k, 1);
    ex = An;
    ex(k, 1:8) = -ex(k, 1:8);
    t_is(got, ex, 12, t);

    t = 'set_reorder(A, Bshort, k, 2) : numeric';
    B = An(1:8, k) * -1;
    got = set_reorder(An, B, k, 2);
    ex = An;
    ex(1:8, k) = -ex(1:8, k);
    t_is(got, ex, 12, t);

    t = 'set_reorder(A, Blong, k, 1) : numeric';
    k = (2:2:10)';
    B = An(k, :) * -1;
    B = [B B];
    got = set_reorder(An, B, k, 1);
    ex = [An zeros(size(An))];
    ex(k, :) = [-An(k, :) -An(k, :)];
    t_is(got, ex, 12, t);

    t = 'set_reorder(A, Blong, k, 2) : numeric';
    B = An(:, k) * -1;
    B = [B; B];
    got = set_reorder(An, B, k, 2);
    ex = [An; zeros(size(An))];
    ex(:, k) = [-An(:, k); -An(:, k)];
    t_is(got, ex, 12, t);

    %%-----  get_reorder (cell)  -----
    k = [3; 7; 4; 1];
    t = 'get_reorder(A, k, 1) : cell';
    t_is(cellfun(@str2num, get_reorder(As, k, 1)), cellfun(@str2num, As(k, :)), 12, t);

    t = 'get_reorder(A, k, 2) : cell';
    t_is(cellfun(@str2num, get_reorder(As, k, 2)), cellfun(@str2num, As(:, k)), 12, t);

    %%-----  set_reorder (cell)  -----
    k = (2:2:10)';
    t = 'set_reorder(A, B, k, 1) : cell';
    B = cellfun(@num2str, num2cell(An(k, :) * -1), 'UniformOutput', 0);
    got = set_reorder(As, B, k, 1);
    ex = As;
    ex(k, :) = cellfun(@num2str, num2cell(-cellfun(@str2num, ex(k, :))), 'UniformOutput', 0);
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);

    t = 'set_reorder(A, B, k, 2) : cell';
    B = cellfun(@num2str, num2cell(An(:, k) * -1), 'UniformOutput', 0);
    got = set_reorder(As, B, k, 2);
    ex = As;
    ex(:, k) = cellfun(@num2str, num2cell(-cellfun(@str2num, ex(:, k))), 'UniformOutput', 0);
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);

    t = 'set_reorder(A, Bshort, k, 1) : cell';
    k = (2:2:8)';
    B = cellfun(@num2str, num2cell(An(k, 1:8) * -1), 'UniformOutput', 0);
    got = set_reorder(As, B, k, 1);
    ex = As;
    ex(k, 1:8) = cellfun(@num2str, num2cell(-cellfun(@str2num, ex(k, 1:8))), 'UniformOutput', 0);
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);

    t = 'set_reorder(A, Bshort, k, 2) : cell';
    B = cellfun(@num2str, num2cell(An(1:8, k) * -1), 'UniformOutput', 0);
    got = set_reorder(As, B, k, 2);
    ex = As;
    ex(1:8, k) = cellfun(@num2str, num2cell(-cellfun(@str2num, ex(1:8, k))), 'UniformOutput', 0);
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);

    t = 'set_reorder(A, Blong, k, 1) : cell';
    k = (2:2:10)';
    B = cellfun(@num2str, num2cell(An(k, :) * -1), 'UniformOutput', 0);
    B = [B B];
    got = set_reorder(As, B, k, 1);
    ex = [As cell(size(As))];
    ex(k, :) = cellfun(@num2str, num2cell([-An(k, :) -An(k, :)]), 'UniformOutput', 0);
    for i = 1:size(got, 1)      %% replace [] with '-999' to make str2num happy
        for j = 1:size(got, 2)
            if isempty(got{i, j})
                got{i, j} = '-999';
            end
            if isempty(ex{i, j})
                ex{i, j} = '-999';
            end
        end
    end
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);

    t = 'set_reorder(A, Blong, k, 2) : cell';
    B = cellfun(@num2str, num2cell(An(:, k) * -1), 'UniformOutput', 0);
    B = [B; B];
    got = set_reorder(As, B, k, 2);
    ex = [As; cell(size(As))];
    ex(:, k) = cellfun(@num2str, num2cell([-An(:, k); -An(:, k)]), 'UniformOutput', 0);
    for i = 1:size(got, 1)      %% replace [] with '-999' to make str2num happy
        for j = 1:size(got, 2)
            if isempty(got{i, j})
                got{i, j} = '-999';
            end
            if isempty(ex{i, j})
                ex{i, j} = '-999';
            end
        end
    end
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);

    %%-----  mpc = ext2int/int2ext(mpc)  -----
    t = 'mpc = ext2int(mpc) : ';
    mpc = ext2int(mpce);
    t_is(mpc.bus, mpci.bus, 12, [t 'bus']);
    t_is(mpc.branch, mpci.branch, 12, [t 'branch']);
    t_is(mpc.gen, mpci.gen, 12, [t 'gen']);
    t_is(mpc.gencost, mpci.gencost, 12, [t 'gencost']);
    t_ok(isequaln(mpc.bus_name, mpci.bus_name), [t 'bus_name']);
    t_ok(isequaln(mpc.gentype, mpci.gentype), [t 'gentype']);
    t_ok(isequaln(mpc.genfuel, mpci.genfuel), [t 'genfuel']);
    t_is(mpc.A, mpci.A, 12, [t 'A']);
    t_is(mpc.N, mpci.N, 12, [t 'N']);
    t = 'mpc = ext2int(mpc) - repeat : ';
    mpc = ext2int(mpc);
    t_is(mpc.bus, mpci.bus, 12, [t 'bus']);
    t_is(mpc.branch, mpci.branch, 12, [t 'branch']);
    t_is(mpc.gen, mpci.gen, 12, [t 'gen']);
    t_is(mpc.gencost, mpci.gencost, 12, [t 'gencost']);
    t_ok(isequaln(mpc.bus_name, mpci.bus_name), [t 'bus_name']);
    t_ok(isequaln(mpc.gentype, mpci.gentype), [t 'gentype']);
    t_ok(isequaln(mpc.genfuel, mpci.genfuel), [t 'genfuel']);
    t_is(mpc.A, mpci.A, 12, [t 'A']);
    t_is(mpc.N, mpci.N, 12, [t 'N']);
    t = 'mpc = int2ext(mpc) : ';
    mpc = int2ext(mpc);
    t_is(mpc.bus, mpce.bus, 12, [t 'bus']);
    t_is(mpc.branch, mpce.branch, 12, [t 'branch']);
    t_is(mpc.gen, mpce.gen, 12, [t 'gen']);
    t_is(mpc.gencost, mpce.gencost, 12, [t 'gencost']);
    t_ok(isequaln(mpc.bus_name, mpce.bus_name), [t 'bus_name']);
    t_ok(isequaln(mpc.gentype, mpce.gentype), [t 'gentype']);
    t_ok(isequaln(mpc.genfuel, mpce.genfuel), [t 'genfuel']);
    t_is(mpc.A, mpce.A, 12, [t 'A']);
    t_is(mpc.N, mpce.N, 12, [t 'N']);

    %%-----  mpc = ext2int/int2ext(mpc, mpopt)  -----
    mpopt = mpoption();
    t = 'mpc = ext2int(mpc, mpopt) : ';
    mpc = ext2int(mpce, mpopt);
    t_is(mpc.bus, mpci.bus, 12, [t 'bus']);
    t_is(mpc.branch, mpci.branch, 12, [t 'branch']);
    t_is(mpc.gen, mpci.gen, 12, [t 'gen']);
    t_is(mpc.gencost, mpci.gencost, 12, [t 'gencost']);
    t_ok(isequaln(mpc.bus_name, mpci.bus_name), [t 'bus_name']);
    t_ok(isequaln(mpc.gentype, mpci.gentype), [t 'gentype']);
    t_ok(isequaln(mpc.genfuel, mpci.genfuel), [t 'genfuel']);
    t_is(mpc.A, mpci.A, 12, [t 'A']);
    t_is(mpc.N, mpci.N, 12, [t 'N']);
    t = 'mpc = ext2int(mpc, mpopt) - repeat : ';
    mpc = ext2int(mpc, mpopt);
    t_is(mpc.bus, mpci.bus, 12, [t 'bus']);
    t_is(mpc.branch, mpci.branch, 12, [t 'branch']);
    t_is(mpc.gen, mpci.gen, 12, [t 'gen']);
    t_is(mpc.gencost, mpci.gencost, 12, [t 'gencost']);
    t_ok(isequaln(mpc.bus_name, mpci.bus_name), [t 'bus_name']);
    t_ok(isequaln(mpc.gentype, mpci.gentype), [t 'gentype']);
    t_ok(isequaln(mpc.genfuel, mpci.genfuel), [t 'genfuel']);
    t_is(mpc.A, mpci.A, 12, [t 'A']);
    t_is(mpc.N, mpci.N, 12, [t 'N']);
    t = 'mpc = int2ext(mpc, mpopt) : ';
    mpc = int2ext(mpc, mpopt);
    t_is(mpc.bus, mpce.bus, 12, [t 'bus']);
    t_is(mpc.branch, mpce.branch, 12, [t 'branch']);
    t_is(mpc.gen, mpce.gen, 12, [t 'gen']);
    t_is(mpc.gencost, mpce.gencost, 12, [t 'gencost']);
    t_ok(isequaln(mpc.bus_name, mpce.bus_name), [t 'bus_name']);
    t_ok(isequaln(mpc.gentype, mpce.gentype), [t 'gentype']);
    t_ok(isequaln(mpc.genfuel, mpce.genfuel), [t 'genfuel']);
    t_is(mpc.A, mpce.A, 12, [t 'A']);
    t_is(mpc.N, mpce.N, 12, [t 'N']);

    %%-----  val = e2i_data/i2e_data(mpc, val, ...)  -----
    t = 'val = e2i_data(mpc, val, ''bus'')';
    mpc = ext2int(mpce);
    got = e2i_data(mpc, mpce.xbus, 'bus');
    ex = mpce.xbus;
    ex(6, :) = [];
    t_is(got, ex, 12, t);
    t = 'val = i2e_data(mpc, val, oldval, ''bus'')';
    tmp = ones(size(mpce.xbus));
    tmp(6, :) = mpce.xbus(6, :);
    got = i2e_data(mpc, ex, tmp, 'bus');
    t_is(got, mpce.xbus, 12, t);

    t = 'val = e2i_data(mpc, val, ''bus'', 2)';
    got = e2i_data(mpc, mpce.xbus, 'bus', 2);
    ex = mpce.xbus;
    ex(:, 6) = [];
    t_is(got, ex, 12, t);
    t = 'val = i2e_data(mpc, val, oldval, ''bus'', 2)';
    tmp = ones(size(mpce.xbus));
    tmp(:, 6) = mpce.xbus(:, 6);
    got = i2e_data(mpc, ex, tmp, 'bus', 2);
    t_is(got, mpce.xbus, 12, t);

    t = 'val = e2i_data(mpc, val, ''gen'')';
    got = e2i_data(mpc, mpce.xgen, 'gen');
    ex = mpce.xgen(gen_i2e, :);
    t_is(got, ex, 12, t);
    t = 'val = i2e_data(mpc, val, oldval, ''gen'')';
    tmp = ones(size(mpce.xgen));
    tmp(3, :) = mpce.xgen(3, :);
    got = i2e_data(mpc, ex, tmp, 'gen');
    t_is(got, mpce.xgen, 12, t);

    t = 'val = e2i_data(mpc, val, ''gen'', 2)';
    got = e2i_data(mpc, mpce.xgen, 'gen', 2);
    ex = mpce.xgen(:, gen_i2e);
    t_is(got, ex, 12, t);
    t = 'val = i2e_data(mpc, val, oldval, ''gen'', 2)';
    tmp = ones(size(mpce.xgen));
    tmp(:, 3) = mpce.xgen(:, 3);
    got = i2e_data(mpc, ex, tmp, 'gen', 2);
    t_is(got, mpce.xgen, 12, t);

    t = 'val = e2i_data(mpc, val, ''branch'')';
    got = e2i_data(mpc, mpce.xbranch, 'branch');
    ex = mpce.xbranch;
    ex(7, :) = [];
    t_is(got, ex, 12, t);
    t = 'val = i2e_data(mpc, val, oldval, ''branch'')';
    tmp = ones(size(mpce.xbranch));
    tmp(7, :) = mpce.xbranch(7, :);
    got = i2e_data(mpc, ex, tmp, 'branch');
    t_is(got, mpce.xbranch, 12, t);

    t = 'val = e2i_data(mpc, val, ''branch'', 2)';
    got = e2i_data(mpc, mpce.xbranch, 'branch', 2);
    ex = mpce.xbranch;
    ex(:, 7) = [];
    t_is(got, ex, 12, t);
    t = 'val = i2e_data(mpc, val, oldval, ''branch'', 2)';
    tmp = ones(size(mpce.xbranch));
    tmp(:, 7) = mpce.xbranch(:, 7);
    got = i2e_data(mpc, ex, tmp, 'branch', 2);
    t_is(got, mpce.xbranch, 12, t);

    t = 'val = e2i_data(mpc, val, {''branch'', ''gen'', ''bus''})';
    got = e2i_data(mpc, mpce.xrows, {'branch', 'gen', 'bus'});
    ex = [mpce.xbranch([1:6, 8:10], 1:4); mpce.xgen(gen_i2e, :); mpce.xbus([1:5, 7:10], 1:4); -ones(2, 4)];
    t_is(got, ex, 12, t);
    t = 'val = i2e_data(mpc, val, oldval, {''branch'', ''gen'', ''bus''})';
    tmp1 = ones(size(mpce.xbranch(:, 1:4)));
    tmp1(7, 1:4) = mpce.xbranch(7, 1:4);
    tmp2 = ones(size(mpce.xgen));
    tmp2(3, :) = mpce.xgen(3, :);
    tmp3 = ones(size(mpce.xbus(:, 1:4)));
    tmp3(6, 1:4) = mpce.xbus(6, 1:4);
    tmp = [tmp1; tmp2; tmp3];
    got = i2e_data(mpc, ex, tmp, {'branch', 'gen', 'bus'});
    t_is(got, mpce.xrows, 12, t);

    t = 'val = e2i_data(mpc, val, {''branch'', ''gen'', ''bus''}, 2)';
    got = e2i_data(mpc, mpce.xcols, {'branch', 'gen', 'bus'}, 2);
    ex = [mpce.xbranch([1:6, 8:10], 1:4); mpce.xgen(gen_i2e, :); mpce.xbus([1:5, 7:10], 1:4); -ones(2, 4)]';
    t_is(got, ex, 12, t);
    t = 'val = i2e_data(mpc, val, oldval, {''branch'', ''gen'', ''bus''}, 2)';
    tmp1 = ones(size(mpce.xbranch(:, 1:4)));
    tmp1(7, 1:4) = mpce.xbranch(7, 1:4);
    tmp2 = ones(size(mpce.xgen));
    tmp2(3, :) = mpce.xgen(3, :);
    tmp3 = ones(size(mpce.xbus(:, 1:4)));
    tmp3(6, 1:4) = mpce.xbus(6, 1:4);
    tmp = [tmp1; tmp2; tmp3]';
    got = i2e_data(mpc, ex, tmp, {'branch', 'gen', 'bus'}, 2);
    t_is(got, mpce.xcols, 12, t);

    %%-----  val = e2i_data/i2e_data(mpc, cell, ...)  -----
    t = 'val = e2i_data(mpc, cell, ''bus'')';
    mpc = ext2int(mpce);
    got = e2i_data(mpc, mpce.strbus, 'bus');
    ex = mpce.strbus;
    ex(6, :) = [];
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);
    t = 'val = i2e_data(mpc, cell, oldval, ''bus'')';
    tmp = cell(size(mpce.strbus));
    tmp(6, :) = mpce.strbus(6, :);
    got = i2e_data(mpc, ex, tmp, 'bus');
    t_is(cellfun(@str2num, got), cellfun(@str2num, mpce.strbus), 12, t);

    t = 'val = e2i_data(mpc, cell, ''bus'', 2)';
    got = e2i_data(mpc, mpce.strbus, 'bus', 2);
    ex = mpce.strbus;
    ex(:, 6) = [];
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);
    t = 'val = i2e_data(mpc, cell, oldval, ''bus'', 2)';
    tmp = cell(size(mpce.strbus));
    tmp(:, 6) = mpce.strbus(:, 6);
    got = i2e_data(mpc, ex, tmp, 'bus', 2);
    t_is(cellfun(@str2num, got), cellfun(@str2num, mpce.strbus), 12, t);

    t = 'val = e2i_data(mpc, cell, ''gen'')';
    got = e2i_data(mpc, mpce.strgen, 'gen');
    ex = mpce.strgen(gen_i2e, :);
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);
    t = 'val = i2e_data(mpc, cell, oldval, ''gen'')';
    tmp = cell(size(mpce.strgen));
    tmp(3, :) = mpce.strgen(3, :);
    got = i2e_data(mpc, ex, tmp, 'gen');
    t_is(cellfun(@str2num, got), cellfun(@str2num, mpce.strgen), 12, t);

    t = 'val = e2i_data(mpc, cell, ''gen'', 2)';
    got = e2i_data(mpc, mpce.strgen, 'gen', 2);
    ex = mpce.strgen(:, gen_i2e);
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);
    t = 'val = i2e_data(mpc, cell, oldval, ''gen'', 2)';
    tmp = cell(size(mpce.strgen));
    tmp(:, 3) = mpce.strgen(:, 3);
    got = i2e_data(mpc, ex, tmp, 'gen', 2);
    t_is(cellfun(@str2num, got), cellfun(@str2num, mpce.strgen), 12, t);

    t = 'val = e2i_data(mpc, cell, ''branch'')';
    got = e2i_data(mpc, mpce.strbranch, 'branch');
    ex = mpce.strbranch;
    ex(7, :) = [];
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);
    t = 'val = i2e_data(mpc, cell, oldval, ''branch'')';
    tmp = cell(size(mpce.strbranch));
    tmp(7, :) = mpce.strbranch(7, :);
    got = i2e_data(mpc, ex, tmp, 'branch');
    t_is(cellfun(@str2num, got), cellfun(@str2num, mpce.strbranch), 12, t);

    t = 'val = e2i_data(mpc, cell, ''branch'', 2)';
    got = e2i_data(mpc, mpce.strbranch, 'branch', 2);
    ex = mpce.strbranch;
    ex(:, 7) = [];
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);
    t = 'val = i2e_data(mpc, cell, oldval, ''branch'', 2)';
    tmp = cell(size(mpce.strbranch));
    tmp(:, 7) = mpce.strbranch(:, 7);
    got = i2e_data(mpc, ex, tmp, 'branch', 2);
    t_is(cellfun(@str2num, got), cellfun(@str2num, mpce.strbranch), 12, t);

    t = 'val = e2i_data(mpc, cell, {''branch'', ''gen'', ''bus''})';
    got = e2i_data(mpc, mpce.strrows, {'branch', 'gen', 'bus'});
    ex = [mpce.strbranch([1:6, 8:10], 1:4); mpce.strgen(gen_i2e, :); mpce.strbus([1:5, 7:10], 1:4); cellfun(@num2str, num2cell(-ones(2, 4)), 'UniformOutput', 0)];
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);
    t = 'val = i2e_data(mpc, cell, oldval, {''branch'', ''gen'', ''bus''})';
    tmp1 = cell(size(mpce.strbranch(:, 1:4)));
    tmp1(7, 1:4) = mpce.strbranch(7, 1:4);
    tmp2 = cell(size(mpce.strgen));
    tmp2(3, :) = mpce.strgen(3, :);
    tmp3 = cell(size(mpce.strbus(:, 1:4)));
    tmp3(6, 1:4) = mpce.strbus(6, 1:4);
    tmp = [tmp1; tmp2; tmp3];
    got = i2e_data(mpc, ex, tmp, {'branch', 'gen', 'bus'});
    t_is(cellfun(@str2num, got), cellfun(@str2num, mpce.strrows), 12, t);

    t = 'val = e2i_data(mpc, cell, {''branch'', ''gen'', ''bus''}, 2)';
    got = e2i_data(mpc, mpce.strcols, {'branch', 'gen', 'bus'}, 2);
    ex = [mpce.strbranch([1:6, 8:10], 1:4); mpce.strgen(gen_i2e, :); mpce.strbus([1:5, 7:10], 1:4); cellfun(@num2str, num2cell(-ones(2, 4)), 'UniformOutput', 0)]';
    t_is(cellfun(@str2num, got), cellfun(@str2num, ex), 12, t);
    t = 'val = i2e_data(mpc, cell, oldval, {''branch'', ''gen'', ''bus''}, 2)';
    tmp1 = cell(size(mpce.strbranch(:, 1:4)));
    tmp1(7, 1:4) = mpce.strbranch(7, 1:4);
    tmp2 = cell(size(mpce.strgen));
    tmp2(3, :) = mpce.strgen(3, :);
    tmp3 = cell(size(mpce.strbus(:, 1:4)));
    tmp3(6, 1:4) = mpce.strbus(6, 1:4);
    tmp = [tmp1; tmp2; tmp3]';
    got = i2e_data(mpc, ex, tmp, {'branch', 'gen', 'bus'}, 2);
    t_is(cellfun(@str2num, got), cellfun(@str2num, mpce.strcols), 12, t);

    %%-----  mpc = e2i_field/i2e_field(mpc, field, ...)  -----
    t = 'mpc = e2i_field(mpc, field, ''bus'')';
    mpc = ext2int(mpce);
    ex = mpce.xbus;
    ex(6, :) = [];
    got = e2i_field(mpc, 'xbus', 'bus');
    t_is(got.xbus, ex, 12, t);
    t = 'mpc = i2e_field(mpc, field, ''bus'')';
    got = i2e_field(got, 'xbus', 'bus');
    t_is(got.xbus, mpce.xbus, 12, t);

    t = 'mpc = e2i_field(mpc, field, ''bus'', 2)';
    ex = mpce.xbus;
    ex(:, 6) = [];
    got = e2i_field(mpc, 'xbus', 'bus', 2);
    t_is(got.xbus, ex, 12, t);
    t = 'mpc = i2e_field(mpc, field, ''bus'', 2)';
    got = i2e_field(got, 'xbus', 'bus', 2);
    t_is(got.xbus, mpce.xbus, 12, t);

    t = 'mpc = e2i_field(mpc, field, ''gen'')';
    ex = mpce.xgen(gen_i2e, :);
    got = e2i_field(mpc, 'xgen', 'gen');
    t_is(got.xgen, ex, 12, t);
    t = 'mpc = i2e_field(mpc, field, ''gen'')';
    got = i2e_field(got, 'xgen', 'gen');
    t_is(got.xgen, mpce.xgen, 12, t);

    t = 'mpc = e2i_field(mpc, field, ''gen'', 2)';
    ex = mpce.xgen(:, gen_i2e);
    got = e2i_field(mpc, 'xgen', 'gen', 2);
    t_is(got.xgen, ex, 12, t);
    t = 'mpc = i2e_field(mpc, field, ''gen'', 2)';
    got = i2e_field(got, 'xgen', 'gen', 2);
    t_is(got.xgen, mpce.xgen, 12, t);

    t = 'mpc = e2i_field(mpc, field, ''branch'')';
    ex = mpce.xbranch;
    ex(7, :) = [];
    got = e2i_field(mpc, 'xbranch', 'branch');
    t_is(got.xbranch, ex, 12, t);
    t = 'mpc = i2e_field(mpc, field, ''branch'')';
    got = i2e_field(got, 'xbranch', 'branch');
    t_is(got.xbranch, mpce.xbranch, 12, t);

    t = 'mpc = e2i_field(mpc, field, ''branch'', 2)';
    ex = mpce.xbranch;
    ex(:, 7) = [];
    got = e2i_field(mpc, 'xbranch', 'branch', 2);
    t_is(got.xbranch, ex, 12, t);
    t = 'mpc = i2e_field(mpc, field, ''branch'', 2)';
    got = i2e_field(got, 'xbranch', 'branch', 2);
    t_is(got.xbranch, mpce.xbranch, 12, t);

    t = 'mpc = e2i_field(mpc, field, {''branch'', ''gen'', ''bus''})';
    ex = [mpce.xbranch([1:6, 8:10], 1:4); mpce.xgen(gen_i2e, :); mpce.xbus([1:5, 7:10], 1:4); -ones(2, 4)];
    got = e2i_field(mpc, 'xrows', {'branch', 'gen', 'bus'});
    t_is(got.xrows, ex, 12, t);
    t = 'mpc = i2e_field(mpc, field, {''branch'', ''gen'', ''bus''})';
    got = i2e_field(got, 'xrows', {'branch', 'gen', 'bus'});
    t_is(got.xrows, mpce.xrows, 12, t);

    t = 'mpc = e2i_field(mpc, field, {''branch'', ''gen'', ''bus''}, 2)';
    ex = [mpce.xbranch([1:6, 8:10], 1:4); mpce.xgen(gen_i2e, :); mpce.xbus([1:5, 7:10], 1:4); -ones(2, 4)]';
    got = e2i_field(mpc, 'xcols', {'branch', 'gen', 'bus'}, 2);
    t_is(got.xcols, ex, 12, t);
    t = 'mpc = i2e_field(mpc, field, {''branch'', ''gen'', ''bus''})';
    got = i2e_field(got, 'xcols', {'branch', 'gen', 'bus'}, 2);
    t_is(got.xcols, mpce.xcols, 12, t);

    t = 'mpc = e2i_field(mpc, {''field1'', ''field2''}, ordering)';
    ex = mpce.x.more(gen_i2e, :);
    got = e2i_field(mpc, {'x', 'more'}, 'gen');
    t_is(got.x.more, ex, 12, t);
    t = 'mpc = i2e_field(mpc, {''field1'', ''field2''}, ordering)';
    got = i2e_field(got, {'x', 'more'}, 'gen');
    t_is(got.x.more, mpce.x.more, 12, t);

    t = 'mpc = e2i_field(mpc, {''field1'', ''field2''}, ordering, 2)';
    ex = mpce.x.more(:, gen_i2e);
    got = e2i_field(mpc, {'x', 'more'}, 'gen', 2);
    t_is(got.x.more, ex, 12, t);
    t = 'mpc = i2e_field(mpc, {''field1'', ''field2''}, ordering, 2)';
    got = i2e_field(got, {'x', 'more'}, 'gen', 2);
    t_is(got.x.more, mpce.x.more, 12, t);

    %%-----  mpc = e2i_field/i2e_field(mpc, cellfield, ...)  -----
    t = 'mpc = e2i_field(mpc, cellfield, ''bus'')';
    mpc = ext2int(mpce);
    ex = mpce.strbus;
    ex(6, :) = [];
    got = e2i_field(mpc, 'strbus', 'bus');
    t_is(cellfun(@str2num, got.strbus), cellfun(@str2num, ex), 12, t);
    t = 'mpc = i2e_field(mpc, cellfield, ''bus'')';
    got = i2e_field(got, 'strbus', 'bus');
    t_is(cellfun(@str2num, got.strbus), cellfun(@str2num, mpce.strbus), 12, t);

    t = 'mpc = e2i_field(mpc, cellfield, ''bus'', 2)';
    ex = mpce.strbus;
    ex(:, 6) = [];
    got = e2i_field(mpc, 'strbus', 'bus', 2);
    t_is(cellfun(@str2num, got.strbus), cellfun(@str2num, ex), 12, t);
    t = 'mpc = i2e_field(mpc, cellfield, ''bus'', 2)';
    got = i2e_field(got, 'strbus', 'bus', 2);
    t_is(cellfun(@str2num, got.strbus), cellfun(@str2num, mpce.strbus), 12, t);

    t = 'mpc = e2i_field(mpc, cellfield, ''gen'')';
    ex = mpce.strgen(gen_i2e, :);
    got = e2i_field(mpc, 'strgen', 'gen');
    t_is(cellfun(@str2num, got.strgen), cellfun(@str2num, ex), 12, t);
    t = 'mpc = i2e_field(mpc, cellfield, ''gen'')';
    got = i2e_field(got, 'strgen', 'gen');
    t_is(cellfun(@str2num, got.strgen), cellfun(@str2num, mpce.strgen), 12, t);

    t = 'mpc = e2i_field(mpc, cellfield, ''gen'', 2)';
    ex = mpce.strgen(:, gen_i2e);
    got = e2i_field(mpc, 'strgen', 'gen', 2);
    t_is(cellfun(@str2num, got.strgen), cellfun(@str2num, ex), 12, t);
    t = 'mpc = i2e_field(mpc, cellfield, ''gen'', 2)';
    got = i2e_field(got, 'strgen', 'gen', 2);
    t_is(cellfun(@str2num, got.strgen), cellfun(@str2num, mpce.strgen), 12, t);

    t = 'mpc = e2i_field(mpc, cellfield, ''branch'')';
    ex = mpce.strbranch;
    ex(7, :) = [];
    got = e2i_field(mpc, 'strbranch', 'branch');
    t_is(cellfun(@str2num, got.strbranch), cellfun(@str2num, ex), 12, t);
    t = 'mpc = i2e_field(mpc, cellfield, ''branch'')';
    got = i2e_field(got, 'strbranch', 'branch');
    t_is(cellfun(@str2num, got.strbranch), cellfun(@str2num, mpce.strbranch), 12, t);

    t = 'mpc = e2i_field(mpc, cellfield, ''branch'', 2)';
    ex = mpce.strbranch;
    ex(:, 7) = [];
    got = e2i_field(mpc, 'strbranch', 'branch', 2);
    t_is(cellfun(@str2num, got.strbranch), cellfun(@str2num, ex), 12, t);
    t = 'mpc = i2e_field(mpc, cellfield, ''branch'', 2)';
    got = i2e_field(got, 'strbranch', 'branch', 2);
    t_is(cellfun(@str2num, got.strbranch), cellfun(@str2num, mpce.strbranch), 12, t);

    t = 'mpc = e2i_field(mpc, cellfield, {''branch'', ''gen'', ''bus''})';
    ex = [mpce.strbranch([1:6, 8:10], 1:4); mpce.strgen(gen_i2e, :); mpce.strbus([1:5, 7:10], 1:4); cellfun(@num2str, num2cell(-ones(2, 4)), 'UniformOutput', 0)];
    got = e2i_field(mpc, 'strrows', {'branch', 'gen', 'bus'});
    t_is(cellfun(@str2num, got.strrows), cellfun(@str2num, ex), 12, t);
    t = 'mpc = i2e_field(mpc, cellfield, {''branch'', ''gen'', ''bus''})';
    got = i2e_field(got, 'strrows', {'branch', 'gen', 'bus'});
    t_is(cellfun(@str2num, got.strrows), cellfun(@str2num, mpce.strrows), 12, t);

    t = 'mpc = e2i_field(mpc, cellfield, {''branch'', ''gen'', ''bus''}, 2)';
    ex = [mpce.strbranch([1:6, 8:10], 1:4); mpce.strgen(gen_i2e, :); mpce.strbus([1:5, 7:10], 1:4); cellfun(@num2str, num2cell(-ones(2, 4)), 'UniformOutput', 0)]';
    got = e2i_field(mpc, 'strcols', {'branch', 'gen', 'bus'}, 2);
    t_is(cellfun(@str2num, got.strcols), cellfun(@str2num, ex), 12, t);
    t = 'mpc = i2e_field(mpc, cellfield, {''branch'', ''gen'', ''bus''})';
    got = i2e_field(got, 'strcols', {'branch', 'gen', 'bus'}, 2);
    t_is(cellfun(@str2num, got.strcols), cellfun(@str2num, mpce.strcols), 12, t);

    %%-----  more mpc = ext2int/int2ext(mpc)  -----
    t = 'mpc = ext2int(mpc) - bus/gen/branch only : ';
    mpce = loadcase('t_case_ext');
    mpci = loadcase('t_case_int');
    mpce = rmfield(mpce, 'gencost');
    mpce = rmfield(mpce, 'A');
    mpce = rmfield(mpce, 'N');
    mpci = rmfield(mpci, 'gencost');
    mpci = rmfield(mpci, 'A');
    mpci = rmfield(mpci, 'N');
    mpc = ext2int(mpce);
    t_is(mpc.bus, mpci.bus, 12, [t 'bus']);
    t_is(mpc.branch, mpci.branch, 12, [t 'branch']);
    t_is(mpc.gen, mpci.gen, 12, [t 'gen']);

    t = 'mpc = ext2int(mpc) - Qg cost, no N : ';
    mpce = loadcase('t_case_ext');
    mpci = loadcase('t_case_int');
    mpce = rmfield(mpce, 'N');
    mpci = rmfield(mpci, 'N');
    mpce.gencost = [mpce.gencost; mpce.gencost];
    mpci.gencost = [mpci.gencost; mpci.gencost];
    mpc = ext2int(mpce);
    t_is(mpc.bus, mpci.bus, 12, [t 'bus']);
    t_is(mpc.branch, mpci.branch, 12, [t 'branch']);
    t_is(mpc.gen, mpci.gen, 12, [t 'gen']);
    t_is(mpc.gencost, mpci.gencost, 12, [t 'gencost']);
    t_is(mpc.A, mpci.A, 12, [t 'A']);

    t = 'mpc = ext2int(mpc) - A, N are DC sized : ';
    mpce = loadcase('t_case_ext');
    mpci = loadcase('t_case_int');
    eVmQgcols = [11:20 25:28]';
    iVmQgcols = [10:18 22:24]';
    mpce.A(:, eVmQgcols) = [];
    mpce.N(:, eVmQgcols) = [];
    mpci.A(:, iVmQgcols) = [];
    mpci.N(:, iVmQgcols) = [];
    mpc = ext2int(mpce);
    t_is(mpc.bus, mpci.bus, 12, [t 'bus']);
    t_is(mpc.branch, mpci.branch, 12, [t 'branch']);
    t_is(mpc.gen, mpci.gen, 12, [t 'gen']);
    t_is(mpc.gencost, mpci.gencost, 12, [t 'gencost']);
    t_is(mpc.A, mpci.A, 12, [t 'A']);
    t_is(mpc.N, mpci.N, 12, [t 'N']);

    t = 'mpc = int2ext(mpc) - A, N are DC sized : ';
    mpc = int2ext(mpc);
    t_is(mpc.bus, mpce.bus, 12, [t 'bus']);
    t_is(mpc.branch, mpce.branch, 12, [t 'branch']);
    t_is(mpc.gen, mpce.gen, 12, [t 'gen']);
    t_is(mpc.gencost, mpce.gencost, 12, [t 'gencost']);
    t_is(mpc.A, mpce.A, 12, [t 'A']);
    t_is(mpc.N, mpce.N, 12, [t 'N']);

    t = 'mpc = ext2int(mpc) - all buses isolated : ';
    mpc = loadcase('t_case_ext');
    mpc.bus(:, BUS_TYPE) = NONE;
    try
        mpci = ext2int(mpc);
        t_is(size(mpci.bus, 1), 0, 12, [t 'internal case empty']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
    end

    t = 'mpc = int2ext(mpc) - all buses isolated : ';
    try
        mpce = int2ext(mpci);
        t_is(mpc.bus, mpce.bus, 12, [t 'bus']);
        t_is(mpc.branch, mpce.branch, 12, [t 'branch']);
        t_is(mpc.gen, mpce.gen, 12, [t 'gen']);
        t_is(mpc.gencost, mpce.gencost, 12, [t 'gencost']);
        t_is(mpc.A, mpce.A, 12, [t 'A']);
        t_is(mpc.N, mpce.N, 12, [t 'N']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
        t_skip(5, [t 'unexpected fatal error']);
    end
end

t_end;

function mpc = userfcn_test_ext2int(mpc, mpopt, args)

function results = userfcn_test_int2ext(results, mpopt, args)

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%%*****************************************************************
%% HSDHKMcorr: corrector step for the HKM direction.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function [par,dX,dy,dZ,resnrm] = HSDHKMcorr(blk,At,par,rp,Rd,sigmu,hRd,...
            dX,dZ,coeff,L,X,Z);

    global printlevel
    global matfct_options solve_ok
%%
    [rhs,EinvRc] = HSDHKMrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ);
    m = length(rp); ncolU = size(coeff.mat12,2);
    rhs = [rhs; zeros(m+ncolU-length(rhs),1)];
%%
    solve_ok = 1; resnrm = norm(rhs);
    [xx,resnrm,solve_ok] = HSDbicgstab(coeff,rhs,L,[],[],printlevel);
    if (solve_ok<=0) & (printlevel)
       fprintf('\n  warning: iterative solver fails: %3.1f.',solve_ok);
    end
    if (par.printlevel>=3); fprintf('%2.0f ',length(resnrm)-1); end
%%
    [par,dX,dy,dZ] = HSDHKMdirfun(blk,At,par,Rd,EinvRc,X,xx);
%%*****************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [verbose,verboseI,debug,doPlot,maxFunEvals,maxIter,optTol,progTol,method,...
    corrections,c1,c2,LS_init,cgSolve,qnUpdate,cgUpdate,initialHessType,...
    HessianModify,Fref,useComplex,numDiff,LS_saveHessianComp,...
    Damped,HvFunc,bbType,cycle,...
    HessianIter,outputFcn,useMex,useNegCurv,precFunc,...
    LS_type,LS_interp,LS_multi,DerivativeCheck] = ...
    minFunc_processInputOptions(o)

% Constants
SD = 0;
CSD = 1;
BB = 2;
CG = 3;
PCG = 4;
LBFGS = 5;
QNEWTON = 6;
NEWTON0 = 7;
NEWTON = 8;
TENSOR = 9;

verbose = 1;
verboseI= 1;
debug = 0;
doPlot = 0;
method = LBFGS;
cgSolve = 0;

o = toUpper(o);

if isfield(o,'DISPLAY')
    switch(upper(o.DISPLAY))
        case 0
            verbose = 0;
            verboseI = 0;
        case 'FINAL'
            verboseI = 0;
        case 'OFF'
            verbose = 0;
            verboseI = 0;
        case 'NONE'
            verbose = 0;
            verboseI = 0;
        case 'FULL'
            debug = 1;
        case 'EXCESSIVE'
            debug = 1;
            doPlot = 1;
    end
end

DerivativeCheck = 0;
if isfield(o,'DERIVATIVECHECK')
    switch(upper(o.DERIVATIVECHECK))
        case 1
            DerivativeCheck = 1;
        case 'ON'
            DerivativeCheck = 1;
    end
end

LS_init = 0;
LS_type = 1;
LS_interp = 2;
LS_multi = 0;
Fref = 1;
Damped = 0;
HessianIter = 1;
c2 = 0.9;
if isfield(o,'METHOD')
    m = upper(o.METHOD);
    switch(m)
        case 'TENSOR'
            method = TENSOR;
        case 'NEWTON'
            method = NEWTON;
        case 'MNEWTON'
            method = NEWTON;
            HessianIter = 5;
        case 'PNEWTON0'
            method = NEWTON0;
            cgSolve = 1;
        case 'NEWTON0'
            method = NEWTON0;
        case 'QNEWTON'
            method = QNEWTON;
            Damped = 1;
        case 'LBFGS'
            method = LBFGS;
        case 'BB'
            method = BB;
            LS_type = 0;
            Fref = 20;
        case 'PCG'
            method = PCG;
            c2 = 0.2;
            LS_init = 2;
        case 'SCG'
            method = CG;
            c2 = 0.2;
            LS_init = 4;
        case 'CG'
            method = CG;
            c2 = 0.2;
            LS_init = 2;
        case 'CSD'
            method = CSD;
            c2 = 0.2;
            Fref = 10;
            LS_init = 2;
        case 'SD'
            method = SD;
            LS_init = 2;
    end
end

maxFunEvals = getOpt(o,'MAXFUNEVALS',1000);
maxIter = getOpt(o,'MAXITER',500);
optTol = getOpt(o,'OPTTOL',1e-5);
progTol = getOpt(o,'PROGTOL',1e-9);
corrections = getOpt(o,'CORRECTIONS',100);
corrections = getOpt(o,'CORR',corrections);
c1 = getOpt(o,'C1',1e-4);
c2 = getOpt(o,'C2',c2);
LS_init = getOpt(o,'LS_INIT',LS_init);
cgSolve = getOpt(o,'CGSOLVE',cgSolve);
qnUpdate = getOpt(o,'QNUPDATE',3);
cgUpdate = getOpt(o,'CGUPDATE',2);
initialHessType = getOpt(o,'INITIALHESSTYPE',1);
HessianModify = getOpt(o,'HESSIANMODIFY',0);
Fref = getOpt(o,'FREF',Fref);
useComplex = getOpt(o,'USECOMPLEX',0);
numDiff = getOpt(o,'NUMDIFF',0);
LS_saveHessianComp = getOpt(o,'LS_SAVEHESSIANCOMP',1);
Damped = getOpt(o,'DAMPED',Damped);
HvFunc = getOpt(o,'HVFUNC',[]);
bbType = getOpt(o,'BBTYPE',0);
cycle = getOpt(o,'CYCLE',3);
HessianIter = getOpt(o,'HESSIANITER',HessianIter);
outputFcn = getOpt(o,'OUTPUTFCN',[]);
useMex = getOpt(o,'USEMEX',1);
useNegCurv = getOpt(o,'USENEGCURV',1);
precFunc = getOpt(o,'PRECFUNC',[]);
LS_type = getOpt(o,'LS_type',LS_type);
LS_interp = getOpt(o,'LS_interp',LS_interp);
LS_multi = getOpt(o,'LS_multi',LS_multi);
end

function [v] = getOpt(options,opt,default)
if isfield(options,opt)
    if ~isempty(getfield(options,opt))
        v = getfield(options,opt);
    else
        v = default;
    end
else
    v = default;
end
end

function [o] = toUpper(o)
if ~isempty(o)
    fn = fieldnames(o);
    for i = 1:length(fn)
        o = setfield(o,upper(fn{i}),getfield(o,fn{i}));
    end
end
end

author: zouchuhang

repo: LayoutNet


text:
classdef MimZoomTool < MimTool
    % MimZoomTool. A tool for interactively zooming an image with MimViewerPanel
    %
    %     MimZoomTool is a tool class used with MimViewerPanel to allow the user
    %     to zoom an image using mouse and keyboard controls.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Zoom'
        Cursor = 'circle'
        ToolTip = 'Zoom tool'
        Tag = 'Zoom'
        ShortcutKey = ''
    end

    properties (Access = private)
        Callback
    end

    properties (Access = private)
        StartImageLimitsMin
        StartImageLimitsMax
        ImageParameters
        StartCoords
    end

    methods
        function obj = MimZoomTool(image_parameters, callback)
            obj.Callback = callback;
            obj.ImageParameters = image_parameters;
        end

        function MouseDragged(obj, screen_coords, last_coords)
            if ~isempty(obj.StartCoords) && ~isempty(obj.StartImageLimitsMin)
                min_coords = obj.StartImageLimitsMin;
                max_coords = obj.StartImageLimitsMax;
                x_range = max_coords(1) - min_coords(1);
                y_range = max_coords(2) - min_coords(2);

                coords_offset = screen_coords - obj.StartCoords;
                y_relative_movement = coords_offset(2)/y_range;

                relative_scale = y_relative_movement;
                x_range_scale = relative_scale*x_range/1.5;
                y_range_scale = relative_scale*y_range/1.5;

                x_lim = [min_coords(1) - x_range_scale, max_coords(1) + x_range_scale];

                y_lim = [min_coords(2) - y_range_scale, max_coords(2) + y_range_scale];

                if (abs(x_lim(2) - x_lim(1)) > 10) && (abs(y_lim(2) - y_lim(1)) > 10) && ...
                        (abs(x_lim(2) - x_lim(1)) < 1000) && (abs(y_lim(2) - y_lim(1)) < 800)
                    obj.Callback.SetImageLimits([x_lim(1), y_lim(1)], [x_lim(2), y_lim(2)]);
                end
            end
        end

        function MouseDown(obj, screen_coords)
            [obj.StartImageLimitsMin, obj.StartImageLimitsMax] = obj.Callback.GetImageLimits();
            obj.StartCoords = screen_coords;
        end

        function NewOrientation(obj)
            obj.StartCoords = [];
        end

        function ImageChanged(obj)
            obj.StartImageLimitsMin = [];
            obj.StartImageLimitsMax = [];
            obj.StartCoords = [];
        end

        function Enter(obj)
            obj.StartImageLimitsMin = [];
            obj.StartImageLimitsMax = [];
            obj.StartCoords = [];
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function value = fedn ( dim_num, x )

%*****************************************************************************80
%
%% FEDN(X(1:DIM_NUM)) = exp ( sum ( X(1:DIM_NUM) ) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real X(DIM_NUM), the argument.
%
%    Output, real VALUE, the value of the function at X.
%
  value = exp ( sum ( x(1:dim_num) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function option_num = p14_option_num ( )

%*****************************************************************************80
%
%% P14_OPTION_NUM returns the number of options for problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer OPTION_NUM, the number of options.
%
  option_num = 1;

  return
end


text:
function ch = i4_to_hex_digit ( i )

%*****************************************************************************80
%
%% I4_TO_HEX_DIGIT converts a (small) I4 to a hexadecimal digit.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the integer, between 0 and 15.
%
%    Output, character CH, the hexadecimal representation of the integer.
%
  if ( 0 <= i && i <= 9 )
    ch = '0' + i;
  elseif ( 10 <= i && i <= 15 )
    ch = 'A' + i - 10;
  else
    ch = '?';
  end

  return
end


text:
%EEGLAB time-frequency functions (timefreqfunc folder):
%  angtimewarp          - Given two event marker vectors, computes a...
%  bootstat             - Accumulate surrogate data to assess significance by permutation of some...
%  correct_mc           - Compute an upper limit for the number of independent...
%  correctfit           - Correct fit using observed p-values. Use this function...
%  crossf               - Returns estimates and plots event-related coherence (ERCOH)...
%  dftfilt              - Discrete Fourier filter...
%  dftfilt2             - Discrete complex wavelet filters...
%  dftfilt3             - Discrete complex wavelet filters...
%  newcrossf            - Returns estimates and plots event-related coherence (ERCOH)...
%  newtimef             - Return estimates and plots of mean event-related (log) spectral...
%  newtimefbaseln       - Remove baseline power values for newtimef. This...
%  newtimefitc          - Function to compute inter-trial coherence (phase locking...
%  newtimefpowerunit    - Find power unit for y-axis based on input structure.
%  newtimeftrialbaseln  - Remove baseline power values for single trials in...
%  pac                  - Compute phase-amplitude coupling (power of first input...
%  pac_cont             - Compute phase-amplitude coupling (power of first input...
%  rsadjust             - Adjust l-values (Ramberg-Schmeiser distribution)...
%  rsfit                - Find p value for a given value in a given distribution...
%  rsget                - Get the p-value for a given collection of l-values...
%  rspdfsolv            - Sub-function used by RSFIT to searc for optimal...
%  rspfunc              - Sub-function used by RSGET...
%  timef                - Returns estimates and plots of mean event-related spectral...
%  timefreq             - Compute time/frequency decomposition of data trials. This...
%  timewarp             - Given two event marker vectors, computes a matrix...

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [newModel,AddedExchRxn] = addExchangeRxn(model,metList,lb,ub)
%addExchangeRxn adds exchange reactions
%
% newModel = addExchangeRxn(model,metList,lb,ub)
%
%INPUTS
% model         Cobra model structure
% metList       List of metabolites
%
%OPTIONAL INPUTS
% lb            Array of lower bounds
% ub            Array of upper bounds
%
%OUTPUT
% newModel      COBRA model with added exchange reactions
%
% Ines Thiele 02/2009

if nargin < 3
    lb = ones(length(metList),1)*min(model.lb);
end
if nargin < 4
    ub = ones(length(metList),1)*max(model.ub);
end
Revs = zeros(length(metList),1);
Revs(lb<0) = 1;

newModel = model;
AddedExchRxn = '';
for i = 1 : length(metList)
    [newModel,rxnIDexists] = addReactionOri(newModel,strcat('EX_',metList{i}),metList(i),-1,Revs(i),...
        lb(i),ub(i));
    AddedExchRxn=[AddedExchRxn;strcat('EX_',metList(i))];
end

author: opencobra

repo: cobratoolbox


text:
function x_num = box_behnken_size ( dim_num )

%*****************************************************************************80
%
%% BOX_BEHNKEN_SIZE returns the size of a Box-Behnken design.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Box, Donald Behnken,
%    Some new three level designs for the study of quantitative variables,
%    Technometrics,
%    Volume 2, pages 455-475, 1960.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, integer X_NUM, the number of elements of the design.
%    X_NUM will be equal to DIM_NUM * 2**(DIM_NUM-1) + 1.
%
  if ( 1 <= dim_num )
    x_num = 1 + dim_num * 2^( dim_num - 1 );
  else
    x_num = -1;
  end

  return
end


text:
function [index, distance, twoout]=near_lonlat(x,y,x0,y0,dist);
% NEAR_LONLAT finds the indices of (lon,lat) that are closest to the point (lon0,lat0).
%        [index,distance]=near_lonlat(lon,lat,lon0,lat0) finds the closest point and
%                                     the distance(km)
%        [index,distance]=near_lonlat(lon,lat,lon0,lat0,dist) finds all points closer than
%                                     the value of dist(km)

% Alexander Crosby 2011
% Rich Signell 2012: removed double loop, speeded up 1000x
[nx,ny]=size(x);
x2=x0*ones(size(x));
y2=y0*ones(size(y));
x3=[x(:) x2(:)].';
y3=[y(:) y2(:)].';
distance=sw_dist(y3(:),x3(:),'km');
distance=reshape(distance(1:2:end),nx,ny);

if nargin > 4,
  index=find(distance<=dist);     %finds points closer than dist
  [row col] = find(distance<=dist);
else
  index=find(distance==min(min(distance)));  % finds closest point
  [row col]=find(distance==min(min(distance)));
  index=index(1);
end
distance=distance(index);
twoout = [row col];
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function show_match(Ig1, Ig2, loc1, loc2);
n = size(loc1,1);
im3 = appendimages(Ig1,Ig2);
figure('Position', [100 100 size(im3,2) size(im3,1)]);
colormap('gray');
imagesc(im3);
hold on;
cols1 = size(Ig1,2);
for i = 1: n
  %if (match(i) > 0)
    line([loc1(i,2) loc2(i,2)+cols1], ...
         [loc1(i,1) loc2(i,1)], 'Color', 'c', 'LineWidth', 2);
  %end
end
hold off;


text:
function dom = domain_pot(pot)
% DOMAIN_POT Return the domain of this cpot.
% dom = domain_pot(pot)

dom = pot.domain;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function colorbarAxis = findColorbarAxis(obj,colorbarHandle)
    if isHG2
        colorbarAxisIndex = find(arrayfun(@(x)(isequal(getappdata(x.Handle,'ColorbarPeerHandle'),colorbarHandle)),obj.State.Axis));
    else
        colorbarAxisIndex = find(arrayfun(@(x)(isequal(getappdata(x.Handle,'LegendColorbarInnerList'),colorbarHandle) + ...
            isequal(getappdata(x.Handle,'LegendColorbarOuterList'),colorbarHandle)),obj.State.Axis));
    end
    if isempty(colorbarAxisIndex)
        colorbarAxisIndex = 1;
    end
    colorbarAxis = obj.State.Axis(colorbarAxisIndex).Handle;
end

author: goGPS-Project


text:
function [net, sName] = addOneLoop_forMeanShiftGrouping(net, sName, loopIdx, GaussianBandwidth)

if ~exist('GaussianBandwidth', 'var')
    GaussianBandwidth = 0.1;
end
%%
pre_input_layer = sName;

lName = sprintf('loop%d_meanshift_S_is_XX', loopIdx);
net.addLayer(lName, ...
    meanshift_S_is_XX(), ...
    {sName}, lName);
sName = lName;

lName = sprintf('loop%d_meanshift_G_is_Gaussian', loopIdx);
net.addLayer(lName, ...
    meanshift_G_is_Gaussian('delta', GaussianBandwidth), ...
    {sName}, lName);
G_layer = lName;
sName = lName;

lName = sprintf('loop%d_meanshift_d_is_sumG', loopIdx);
net.addLayer(lName, ...
    meanshift_d_is_sumG(), ...
    {sName}, lName);
sName = lName;

lName = sprintf('loop%d_meanshift_q_is_inv_d', loopIdx);
net.addLayer(lName, ...
    meanshift_q_is_inv_d(), ...
    {sName}, lName);
sName = lName;

lName = sprintf('loop%d_meanshift_P_is_G_diag_q', loopIdx);
net.addLayer(lName, ...
    meanshift_P_is_G_diag_q(), ...
    {G_layer, sName}, lName);
sName = lName;

lName = sprintf('loop%d_meanshift_Y_is_XP', loopIdx);
net.addLayer(lName, ...
    meanshift_Y_is_XP(), ...
    {pre_input_layer, sName}, lName);
sName = lName;

author: aimerykong


text:
function [y] = cellcolselect(x, cols)

% [Y] = CELLCOLSELECT(X, COLS) outputs cell-array Y with the same dimensionality as X
% Each cell in Y only contains the columns COLS from the original corresponding cell in X
%
% X (and Y) should be linear cell-array(s) of matrices for which the number of rows
% should be the same for all cells

nx = size(x);
if ~iscell(x) || length(nx)>2 || all(nx>1)
  error('incorrect input for cellrows');
end

if ~iscell(cols)
  y = cellfun(@colc, x, repmat(mat2cell(cols(:),length(cols),1),nx), 'UniformOutput', 0);
else
  y = cellfun(@colc, x, cols, 'UniformOutput', 0);
end

function [y] = colc(x, cols)

y = x(:,cols);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function OutStruct = prtUtilClassParentRelationshipNamesOnly(varargin)
% prtUtilClassParentRelationship(baseClassName)
% xxx Need Help xxx
% Internal

I = meta.class.fromName(varargin{1});
if nargin == 1
    OutStruct.classNames = {I.Name};
%     OutStruct.propNames = cellfun(@(c)c.Name,I.Properties,'uniformOutput',false);
%     OutStruct.methodNames = cellfun(@(c)c.Name,I.Methods,'uniformOutput',false);
%     OutStruct.AllProperties = I.Properties;
%     OutStruct.AllMethods = I.Methods;

elseif nargin == 2
    OutStruct = varargin{2};
    % Add yourself to the list of class Names
    OutStruct.classNames = cat(1,OutStruct.classNames,{I.Name});

%     cSortedProperties = cell(length(OutStruct.propNames),1);
%     for iProp = 1:length(I.Properties)
%         %cInd = find(~cellfun(@isempty, strfind(OutStruct.propNames,I.Properties{iProp}.Name)),1);
%         cInd = find(cellfun(@(s)isequal(I.Properties{iProp}.Name,s),OutStruct.propNames),1);
%
%         if isempty(cInd)
%             error('Maximum property assumption invalid');
%         end
%
%         cSortedProperties{cInd} = I.Properties{iProp};
%     end
%
%     cSortedMethods = cell(length(OutStruct.methodNames),1);
%     for iMeth = 1:length(I.Methods)
%         if strcmpi(I.Methods{iMeth}.Name,varargin{1})
%             continue
%         end
%         if strcmpi(I.Methods{iMeth}.Access,'private');
%             continue
%         end
%
%         %cInd = find(~cellfun(@isempty, strfind(OutStruct.methodNames,I.Methods{iMeth}.Name)),1);
%         cInd = find(cellfun(@(s)isequal(I.Methods{iMeth}.Name,s),OutStruct.methodNames),1);
%         if isempty(cInd)
%             error('Maximum method assumption invalid');
%         end
%
%         cSortedMethods{cInd} = I.Methods{iMeth};
%     end
%
%     OutStruct.AllProperties = cat(2,OutStruct.AllProperties,cSortedProperties);
%     OutStruct.AllMethods = cat(2,OutStruct.AllMethods,cSortedMethods);
end

for iParent = 1:length(I.SuperClasses);
    OutStruct = prtUtilClassParentRelationshipNamesOnly(I.SuperClasses{iParent}.Name,OutStruct);
end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
%% Example: Carbon Ion Treatment Plan
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2017 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
% In this example we will show
% (i) how to load patient data into matRad
% (ii) how to setup a carbon ion dose calculation plan including variable RBE optimization
% (iii) how to inversely optimize the pencil beam intensities based on the
% RBE-weighted dose
% (iv) how to inversely optimize the pencil beam intensities based on the
% biological effect
% (v) how to change the tissues' radiobiological characteristics
% (vi) how to recalculated the dose considering the previously optimized pencil beam intensities
% (vii) how to compare the two results

%% Patient Data Import
% Let's begin with a clear Matlab environment and import the liver
% patient into your workspace.

matRad_rc; %If this throws an error, run it from the parent directory first to set the paths

load('LIVER.mat');

%% Treatment Plan
% The next step is to define your treatment plan labeled as 'pln'. This
% structure requires input from the treatment planner and defines the most
% important cornerstones of your treatment plan.
%%
% First of all, we need to define what kind of radiation modality we would
% like to use. Possible values are photons, protons or carbon. In this
% example we would like to use carbon ions for treatment planning. Next, we
% need to define a treatment machine to correctly load the corresponding
% base data. matRad features generic base data in the file
% 'carbon_Generic.mat'; consequently the machine has to be set accordingly
pln.radiationMode = 'carbon';
pln.machine       = 'Generic';

%%
% Define the flavor of biological optimization for treatment planning along
% with the quantity that should be used for optimization. Possible values
% are (none: physical optimization; const_RBExD: constant RBE of 1.1;
% LEMIV_effect: effect-based optimization; LEMIV_RBExD: optimization of
% RBE-weighted dose. As we use carbon ions, we decide to use base data from
% the local effect model IV and want to optimize the RBE-weighted dose.
% Therefore we set bioOptimization to LEMIV_RBExD
pln.propOpt.bioOptimization = 'LEMIV_RBExD';

%%
% The remaining plan parameters are set like in the previous example files
pln.numOfFractions        = 30;
pln.propStf.gantryAngles  = 315;
pln.propStf.couchAngles   = 0;
pln.propStf.bixelWidth    = 3;
pln.propStf.numOfBeams    = numel(pln.propStf.gantryAngles);
pln.propStf.isoCenter     = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
pln.propOpt.runDAO        = 0;
pln.propOpt.runSequencing = 0;

% dose calculation settings
pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]
pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]
pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm]

%% Generate Beam Geometry STF
stf = matRad_generateStf(ct,cst,pln);

%%
% Let's have a closer look on the stf.ray sub-structure which contains the
% actual  beam/ray geometry information. For illustration purposes we want
% to show the last ray. Besides geometrical  information about the position
% and orientation of the ray, we can also find pencil beam information. If
% the ray coincides with the target, pencil beams were defined along the
% ray from target entry to target exit.
display(stf.ray(end));

%%
% Here are the energies selected on the last ray:
display(stf.ray(end).energy);

%% Dose Calculation
dij = matRad_calcParticleDose(ct,stf,pln,cst);

%% Inverse Optimization  for IMPT based on RBE-weighted dose
% The goal of the fluence optimization is to find a set of bixel/spot
% weights which yield the best possible dose distribution according to the
% clinical objectives and constraints underlying the radiation treatment.
resultGUI = matRad_fluenceOptimization(dij,cst,pln);

%% Plot the Resulting Dose Slice
% Let's plot the transversal iso-center dose slice
slice = round(pln.propStf.isoCenter(3)./ct.resolution.z);
figure,
imagesc(resultGUI.RBExDose (:,:,slice)),colorbar, colormap(jet);

%% Inverse Optimization  for IMPT based on biological effect
% To perform a dose optimization for carbon ions we can also use the
% biological effect instead of the RBE-weighted dose. Therefore we have to
% change the optimization mode and restart the optimization
pln.propOpt.bioOptimization = 'LEMIV_effect';
resultGUI_effect = matRad_fluenceOptimization(dij,cst,pln);

%% Visualize differences
% Through optimzation based on the biological effect we obtain a slightly
% different dose distribution as visualized by the following dose
% difference map
figure;
imagesc(resultGUI.RBExDose (:,:,slice)-resultGUI_effect.RBExDose(:,:,slice));
colorbar;
colormap(jet);

%% Change Radiosensitivity
% The previous treatment plan was optimized using an photon alpha-beta
% ratio of 2 for all tissues. Now, Let's change the radiosensitivity by
% adapting alphaX. This will change the photon alpha-beta ratio
% from 2 to 10.
for i = 1:size(cst,1)
    cst{i,5}.alphaX      = 0.5;
    cst{i,5}.TissueClass = 2;
end

%% Recalculate Plan
% Let's use the existing optimized pencil beam weights and recalculate the RBE weighted dose
resultGUI_tissue = matRad_calcDoseDirect(ct,stf,pln,cst,resultGUI.w);

%% Result Comparison
% Let's compare the new recalculation against the optimization result.
plane = 3;
doseWindow = [0 max([resultGUI_effect.RBExDose(:); resultGUI_tissue.RBExDose(:)])];

figure,
matRad_plotSliceWrapper(gca,ct,cst,1,resultGUI_effect.RBExDose,plane,slice,[],[],colorcube,[],doseWindow,[]);
title('original plan')
figure,
matRad_plotSliceWrapper(gca,ct,cst,1,resultGUI_tissue.RBExDose,plane,slice,[],[],colorcube,[],doseWindow,[]);
title('manipulated plan')
%%
% At this point we would like to see the absolute difference of the original optimization and the
% recalculation.
absDiffCube = resultGUI_effect.RBExDose-resultGUI_tissue.RBExDose;
figure,
matRad_plotSliceWrapper(gca,ct,cst,1,absDiffCube,plane,slice,[],[],colorcube);
title('absolute difference')
%%
% Plot both doses with absolute difference and gamma analysis
[gammaCube,gammaPassRate,hfigure]=matRad_compareDose(resultGUI_effect.RBExDose, resultGUI.RBExDose, ct, cst,[1 1 1],'on');

repo: matRad


text:
function [q, M] = conversionMatrix2Quaternion(M, varargin)
% conversionMatrix2Quaternion - function to convert Rotation Matrix to Quaternion coeff
%
% Syntax:  q = conversionMatrix2Quaternion(M(3,3))
%          q = conversionMatrix2Quaternion(M(1,1),M(1,2),M(1,3),M(2,1),M(2,2),M(2,3),M(3,1),M(3,2),M(3,3))
%          q = conversionMatrix2Quaternion([theta, psi, phi])
%          q = conversionMatrix2Quaternion(theta, psi, phi)
%
% Inputs:
%          M - Rotation Matrix (3,3)
%     *theta - X rotation angle (deg)
%       *psi - Y rotation angle (deg)
%       *phi - Z rotation angle (deg)
%
% Outputs:
%     q - quaternion values [q0 qx qy qz]
%     M - Rotation Matrix used to compute q
%
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%
% See also: none;

% Author: Marco Borges, Ph.D. Student, Computer/Biomedical Engineer
% UFMG, PPGEE, Neurodinamica Lab, Brazil
% email address: marcoafborges@gmail.com
% Website: http://www.cpdee.ufmg.br/
% Reference: http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/
% April 2013; v2; Last revision: 2013-04-18
% Changelog: v2 - add linear velocity extraction

%------------- BEGIN CODE --------------
% M = [m00 m01 m02;
%     m10 m11 m12;
%     m20 m21 m22];

if length(M) == 1 && nargin == 8
    M = [M, varargin{1}, varargin{2}; varargin{3}, varargin{4}, varargin{5}; varargin{6}, varargin{7}, varargin{8}];

elseif length(M) == 9
    M = [M(1), M(2), M(3); M(4), M(5), M(6); M(7), M(8), M(9)];

elseif length(M) == 3 % [theta (x rot), psi (y rot), phi (z rot)]
    theta = M(1); psi = M(2); phi = M(3);
    M = [cos(deg2rad(psi))*cos(deg2rad(phi)), cos(deg2rad(psi))*sin(deg2rad(phi)), -sin(deg2rad(psi));
         (-cos(deg2rad(theta))*sin(deg2rad(phi)))+(sin(deg2rad(theta))*sin(deg2rad(psi))*cos(deg2rad(phi))), (cos(deg2rad(theta))*cos(deg2rad(phi)))+(sin(deg2rad(theta))*sin(deg2rad(psi))*sin(deg2rad(phi))), sin(deg2rad(theta))*cos(deg2rad(psi));
         (sin(deg2rad(theta))*sin(deg2rad(phi)))+(cos(deg2rad(theta))*sin(deg2rad(psi))*cos(deg2rad(phi))), (-sin(deg2rad(theta))*cos(deg2rad(phi)))+(cos(deg2rad(theta))*sin(deg2rad(psi))*sin(deg2rad(phi))), cos(deg2rad(theta))*cos(deg2rad(psi))];

elseif length(M) == 1 && nargin == 3 % (theta (x rot), psi (y rot), phi (z rot))
    theta = M; psi = varargin{1}; phi = varargin{2};
    M = [cos(deg2rad(psi))*cos(deg2rad(phi)), cos(deg2rad(psi))*sin(deg2rad(phi)), -sin(deg2rad(psi));
         (-cos(deg2rad(theta))*sin(deg2rad(phi)))+(sin(deg2rad(theta))*sin(deg2rad(psi))*cos(deg2rad(phi))), (cos(deg2rad(theta))*cos(deg2rad(phi)))+(sin(deg2rad(theta))*sin(deg2rad(psi))*sin(deg2rad(phi))), sin(deg2rad(theta))*cos(deg2rad(psi));
         (sin(deg2rad(theta))*sin(deg2rad(phi)))+(cos(deg2rad(theta))*sin(deg2rad(psi))*cos(deg2rad(phi))), (-sin(deg2rad(theta))*cos(deg2rad(phi)))+(cos(deg2rad(theta))*sin(deg2rad(psi))*sin(deg2rad(phi))), cos(deg2rad(theta))*cos(deg2rad(psi))];
end

tr = M(1,1) + M(2,2) + M(3,3);

if (tr > 0)
    S = sqrt(tr+1.0) * 2; % S=4*qw
    qw = 0.25 * S;
    qx = (M(3,2) - M(2,3)) / S;
    qy = (M(1,3) - M(3,1)) / S;
    qz = (M(2,1) - M(1,2)) / S;
elseif ((M(1,1) > M(2,2)) && (M(1,1) > M(3,3)))
    S = sqrt(1.0 + M(1,1) - M(2,2) - M(3,3)) * 2; % S=4*qx
    qw = (M(3,2) - M(2,3)) / S;
    qx = 0.25 * S;
    qy = (M(1,2) + M(2,1)) / S;
    qz = (M(1,3) + M(3,1)) / S;
elseif (M(2,2) > M(3,3))
    S = sqrt(1.0 + M(2,2) - M(1,1) - M(3,3)) * 2; % S=4*qy
    qw = (M(1,3) - M(3,1)) / S;
    qx = (M(1,2) + M(2,1)) / S;
    qy = 0.25 * S;
    qz = (M(2,3) + M(3,2)) / S;
else
    S = sqrt(1.0 + M(3,3) - M(1,1) - M(2,2)) * 2; % S=4*qz
    qw = (M(2,1) - M(1,2)) / S;
    qx = (M(1,3) + M(3,1)) / S;
    qy = (M(2,3) + M(3,2)) / S;
    qz = 0.25 * S;
end

q = [qw qx qy qz];
%-------------- END CODE ---------------