text:
ngname1 = 'comp.windows.x';
ngname2 = 'comp.graphics';
load_20nggeneral;

author: sidaw

repo: nbsvm

save_path: github-repos/MATLAB/sidaw-nbsvm


text:
function [derivativeMatrix,r8fdx]=firstOrderDerivativeMatrix(xl,xu,n)
% firstOrderDerivativeMatrix precomputes the matrices used to implement numerical differentiation using 9 points.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

%  Compute the spatial increment
dx=(xu-xl)/(n-1);
r8fdx=1./(40320.*dx);
nm4=n-4;

mid_block = zeros(n-8,n);

first_row 	= [-109584 	+322560 	-564480 	+752640 	-705600 	+451584 	-188160 	+46080 	-5040];
second_row 	= [-5040 	-64224 		+141120 	-141120 	+117600 	-70560 		+28224 		-6720 	+720];
third_row 	= [+720 	-11520 		-38304 		+80640 		-50400 		+26880 		-10080 		+2304 	-240];
fourth_row 	= [-240 	+2880 		-20160 		-18144 		+50400 		-20160 		+6720 		-1440 	+144];

i_th_row 	= [+144 	-1536 		+8064 		-32256 		+0 			+32256 		-8064 		+1536 	-144];

n_min_3_row 	= [-144 	+1440 		-6720 		+20160 		-50400 		+18144 		+20160 		-2880 	+240];
n_min_2_row 	= [+240 	-2304 		+10080 		-26880 		+50400 		-80640 		+38304 		+11520 	-720];
n_min_1_row 	= [-720 	+6720 		-28224 		+70560 		-117600 	+141120 	-141120 	+64224 	+5040];
n_min_0_row 	= [+5040 	-46080 		+188160  	-451584 	+705600 	-752640 	+564480 	-322560 +109584];

first_block = [first_row;second_row;third_row;fourth_row];
first_block = [first_block zeros(4,n-9)];

row_index = 1;
for i=5:nm4
    mid_block(row_index,row_index:row_index+8)=i_th_row;
    row_index = row_index+1;
end

last_block = [n_min_3_row;n_min_2_row;n_min_1_row;n_min_0_row];
last_block=[zeros(4,n-9) last_block];

derivativeMatrix = [first_block;mid_block;last_block];

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function [CP2,CNTmat,TOTmat,LENmat] = zero_crossing(M,T,XX,ZIU,ZIL,Z,mu)
% Estimate out of control points and change-point using zero-crossing
% method.
%
% :Usage:
% ::
%
%     [CP2,CNTmat,TOTmat,LENmat] = zero_crossing(M,T,XX,ZIU,ZIL,Z,mu)
%
% Used in ewma5.m  See ewma5 for description of variables.

CP2 = NaN .* zeros(M,1); % Estimation via Zero-crossings
CNTmat  = zeros(M,1);
TOTmat = zeros(M,1);
LENmat = zeros(M,T);

wh = all(XX' - repmat(mean(XX'),size(XX,2),1) <= eps) | any(isnan(XX),2)';
wh = find(~wh);

ooc = ZIU - ZIL;    % pos if occ +, neg if ooc -, zero otherwise
for d=wh,

    % Calculate last 'zero-crossing'

    amiout = find(ooc(d,:));
    if isempty(amiout)
        amiout = 0;
    else
        amiout = amiout(1);
        amiout = ooc(d,amiout);  % first ooc tp, -1 or 1
    end

    switch amiout

        case 1
            dat = Z(d,:)-mu(d);
            tmp = ZIU(d,:);
            [a,b] = max((conv([1 1 1],tmp) == 3));  % must be OOC 3 consecutive points
            [a,zc] = max((dat(1:b) < 0).*(1:b));
            CP2(d) = zc;

        case -1
            dat = Z(d,:)-mu(d);
            tmp = ZIL(d,:);
            [a,b] = max((conv([1 1 1],tmp) == 3));  % must be OOC 3 consecutive points
            [a,zc] = max((dat(1:b) > 0).*(1:b));
            CP2(d) = zc;

        case 0
            CP2(d) = NaN;
            tmp = 0 .* ZIL(d,:);
        otherwise
            error('This should never happen'),keyboard
    end

    % Calculate number of out-of-control runs and their widths

    [cnt tot len] = cnt_runs(tmp);
    CNTmat(d) = cnt;
    TOTmat(d) = tot;
    LENmat(d,:) = len';

end;

CP2(find(CP2 == 1)) = NaN;

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function i4_factorial2_values_test ( )

%*****************************************************************************80
%
%% I4_FACTORIAL2_VALUES_TEST demonstrates the use of I4_FACTORIAL2_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_FACTORIAL2_VALUES_TEST:\n' );
  fprintf ( 1, '  I4_FACTORIAL2_VALUES returns values of\n' );
  fprintf ( 1, '  the double factorial function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N         N!!\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, fn ] = i4_factorial2_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %10d\n',  n, fn );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function g = multimodelLogLikeGradients(model)

% MULTIMODELLOGLIKEGRADIENTS Gradient of MULTIMODEL model log likelihood with respect to parameters.
% FORMAT
% DESC computes the gradient of the multi-task learning wrapper
% model's log likelihood with respect to the parameters.
% ARG model : model structure for which gradients are being
% computed.
% RETURN g : the returned gradients.
%
% SEEALSO multimodelCreate, multimodelLogLikelihood, modelLogLikeGradients
%
% COPYRIGHT : Neil D. Lawrence, 2007, 2008
%
% MODIFICATIONS : Mauricio Alvarez, 2009

% MLTOOLS

  g = zeros(1, model.numParams);
  if numel(model.outputDim) == 1
  endShared = model.numParams - model.numModels*model.numSep;
  endVal = endShared;
  for i = 1:length(model.comp)
    gModel = modelLogLikeGradients(model.comp{i});
    g(1:endShared) = g(1:endShared) + gModel(model.sharedIndices);
    if ~isempty(model.separateIndices)
      startVal = endVal + 1;
      endVal = endVal + model.numSep;
      %g(startVal:endVal) = g(model.separateIndices);
      g(startVal:endVal) = gModel(model.separateIndices);
    end
  end
  else
      startVal = 1;
      endVal = 0;
      for i = 1:length(model.comp)
          endVal = endVal + model.comp{i}.nParams;
          if (i==6) || (i ==14)
              a =1;
          end
          g(startVal:endVal) = modelLogLikeGradients(model.comp{i});
          startVal = endVal + 1;
      end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function OUT = mpower(varargin)

% overloaded power for GRIDobj
%
% Note that there is no matrix power defined for GRIDobj. Thus
% mpower overloads element-by-element powers (power).

funname = 'power';
narginchk(2,2)
OUT = builtincaller(funname,varargin{:});

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function nvar = p01_nvar ( option )

%*****************************************************************************80
%
%% P01_NVAR sets the number of variables for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option chosen for this problem.
%    For some problems, several options are available.  At least,
%    OPTION = 1 is always legal.
%
%    Output, integer NVAR, the number of variables.
%
  nvar = 3;

  return
end


text:
function data = read_line_by_line(FID);

% first pass, find the number of lines
cnt = 1;
while 1
    tmp = fgetl(FID);
    if tmp == -1
        break;
    else
        cnt = cnt + 1;
    end
end

% second pass, read the data
frewind(FID);
data = cell(cnt-1,1);
cnt = 1;
while 1
    tmp = fgetl(FID);
    if tmp == -1
        break;
    else
        data{cnt} = tmp;
        cnt = cnt + 1;
    end
end

author: singaxiong

repo: SignalGraph


text:
function [class_label,pose_label] = sydney_generate_class_and_pose_labels(class_name,rotation,old_style)
%%% Note:
%%% Class Label and Pose Label both start from 0.
%%%
%%% If the parameter "old_style" is given, then all the
%%% single-orientation classes will be labeled as rotation-class = 0

if(~exist('old_style','var')), old_style = 0; end;

classes =  {'4wd','building','bus','car','pedestrian','pillar','pole','traffic_lights','traffic_sign','tree','truck','trunk','ute','van'};

nrot = [18 9 9 18 1 1 1 9 18 1 18 1 18 18];
% nrot = [12 6 6 12 1 1 1 12 12 1 12 1 12 12];
if(old_style), nrot(nrot == 1) = 0; end;

%--- Find the Class Label
[valid_class,class_label] = ismember(class_name,classes);
if(~valid_class)
    class_label = -1;
    pose_label = -1;
    return;
end
class_label = class_label - 1;     %to start from 0

%--- Find the rotation label
cr = cumsum(nrot);
cr = circshift(cr,[0,1]);
cr(1) = 0;
pose_label = cr(class_label+1) + mod(rotation,nrot(class_label+1));

if(old_style)
    pose_label = pose_label+1;
    pose_label(nrot(class_label+1) == 0) = 0;
end

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function [newmodel, HTABLE] = addReactionGEM(model, rxns, rxnNames, rxnFormulas, rev, lb, ub, nRxn, subSystems, grRules, rules, genes, HTABLE)
% Manually adds reactions to a specified model, may add one or more reactions at a time
%
% USAGE:
%
%    [newmodel, HTABLE] = addReactionGEM(model, rxns, rxnNames, rxnFormulas, rev, lb, ub, nRxn, subSystems, grRules, rules, genes, HTABLE)
%
% INPUTS:
%     model:          COBRA model structure
%     rxns:           Identifiers for the reactions
%     rxnNames:       List of reactions
%     rxnFormulas:    reactions' formulas
%     rev:            0 = irrev, 1 = rev
%     lb:             The lower bounds for fluxes
%     ub:             The upper bounds for fluxes
%     subSystems:     subSystem assignment for each reaction, default = ''
%     grRules:        A string representation of the GPR rules defined in a readable format, default = ''
%     rules:          GPR rules in evaluateable format for each reaction, default = ''
%     genes:          Identifiers of the genes in the model, default = ''
%     HTABLE:         hash table
%
% OUTPUTS:
%     newmodel:       changed model
%     HTABLE:         hash table
%
% EXAMPLE:
%
%    [modelLB_NH3] = addReactionSmiley(modelLB, 'NH3r', 'NH3 protonization', cellstr('1 NH3[c] + 1 H[c] <==> 1 NH4[c]'), 1, -1000, 1000, 'Others');
%
% NOTE:
%    Manually add reactions to a specified model, can either add one or
%    multiple reactions at a time
%
%    All syntax standards must comply with the specified model
%
%    For reaction formulas, use: '-->' for irreversible or '<==>' for
%    reversible
%
% .. Author:
%       - Aarash Bordbar 11/2/07 based on AddRxn
%       - IT 11-19-08

if ~exist('subSystems', 'var') || isempty(subSystems)
    clear subSystems;
    for i = 1:length(rev)
        subSystems(i,1) = {{''}};
    end
else
    if ischar([subSystems{:}])
        subSystems = cellfun(@(x) {x},subSystems,'UniformOutput',0);
    end
end
%We later assume, that this field exists, so we have to initialize it
%properly.
if ~isfield(model,'subSystems')
    model.subSystems = repmat({''},size(model.rxns));
end

if ~exist('grRules', 'var') || isempty(grRules)
    clear grRules;
    for i = 1:length(rev)
        grRules(i,1) = {''};
    end
end

%We later assume, that this field exists, so we have to initialize it
%properly.
if ~isfield(model,'grRules')
    model.grRules = repmat({''},size(model.rxns));
end

if ~exist('rules', 'var')|| isempty(rules)
    clear rules;
    for i = 1:length(rev)
        rules(i,1) = {''};
    end
end

if ~exist('genes', 'var')|| isempty(genes)
    genes = cell(0,1);
end

if ~exist('nRxn', 'var') || isempty(nRxn)
    nRxn = length(model.lb)+1;
end
orignMets = length(model.mets);
nMet = orignMets+1;
orignRxns = numel(model.rxns);

if (isfield(model,'genes'))
    nGenes = length(model.genes)+1;
end
newmodel = model;

putNecessary = false;
if exist('HTABLE', 'var') || length(rxnNames) > 2
    useHashTable = true;
    if exist('HTABLE', 'var') % no need to initialize HTABLE
        putNecessary = true; % will need it for output as well.
    else %create new hashtable
        HTABLE = java.util.Hashtable; % initialize HTABLE
        for i = 1:length(newmodel.mets)
            HTABLE.put(newmodel.mets{i}, i);
        end
    end
else
    useHashTable = false;
end

%showprogress(0, 'Adding Rxns ...');
for i = 1:length(rev)
    %IO indicates whether it is part of the reactants or part of the
    %product (ie. IO = -1 is reactant, IO = 1 products)
    IO = -1;
    %newmodel.rxns{nRxn} = char(rxns(i));
    newmodel.lb(nRxn,1) = lb(i,1);
    newmodel.ub(nRxn,1) = ub(i,1);
    newmodel.subSystems{nRxn,1} = subSystems(i,1);
    newmodel.grRules(nRxn,1) = grRules(i,1);
    newmodel.rules(nRxn,1) = rules(i,1);
    newmodel.c(nRxn,1) = 0;
    newmodel.rxnNames{nRxn,1} = char(rxnNames(i,1));
    %newmodel.S(:, i) = zeros(length(newmodel.mets), 1);

    %parses reaction formula into components of string
    [parsing{1,1},parsing{2,1}] = strtok(rxnFormulas{i});
    for j = 2:100
        [parsing{j,1},parsing{j+1,1}] = strtok(parsing{j,1});
        if isempty(parsing{j+1,1})==1
            break
        end
    end

    nComponents = length(parsing)-1;
    j = 1;
    while j <= nComponents
        if strcmp(parsing{j},'+') == 1
            j = j+1;
        elseif strcmp(parsing{j},' ') == 1
            j = j+1;

        %if its '-->' or '<==>' then switches IO to 1 which indicates the
        %next compounds are part of the product
        elseif strcmp(parsing{j},'<==>') == 1
            IO = 1;
            j = j+1;
        elseif strcmp(parsing{j},'=>') == 1
            IO = 1;
            j = j+1;
        elseif strcmp(parsing{j},'-->') == 1
            IO = 1;
            j = j+1;
        elseif str2double(parsing{j}) > 0

            j = j+1;
        else
            %if met already exists in newmodel then sets metLoc to
            %corresponding met location

            if useHashTable
                c = HTABLE.get(parsing{j});
            else
                c = strmatch(parsing{j},newmodel.mets,'exact');
                %display('slow method')
            end

            %c = strcmp(parsing{j},newmodel.mets);
            if ~isempty(c)
                metLoc = c;
                parsing(j);

            %if met doesn't exist then metLoc is set at nMet (the end of list)
            else
                metLoc = nMet;
                newmodel.mets{metLoc,1} = parsing{j};
                newmodel.S(metLoc,:) = 0;
                newmodel.b(metLoc,1) = 0;
                if putNecessary
                    HTABLE.put(parsing{j}, metLoc);
                end
                nMet = nMet+1;
            end
            %finds reaction coefficient of met
            rxnCoeff = str2double(parsing{j-1});
            if isnan(rxnCoeff)
                rxnCoeff=1;
                rxnNames(i,1);
            end
            s = size(model.S);
            if(metLoc <= s(1))
                origCoeff = newmodel.S(metLoc,nRxn);
                if(origCoeff < 0)
                    newmodel.rxns(nRxn);
                end
            else
                origCoeff = 0;
            end
            newmodel.S(metLoc,nRxn) = rxnCoeff*IO + origCoeff;

            j=j+1;
        end
    end

    clear parsing

    [parsing{1,1},parsing{2,1}] = strtok(grRules{i});
    if ~isempty(parsing{2,1}) %length(parsing{2,1}) ~= 0
        for j = 2:100
            [parsing{j,1},parsing{j+1,1}] = strtok(parsing{j,1});
            if isempty(parsing{j+1,1})==1
                break
            end
        end
    end
    nRxn = nRxn + 1;

    clear parsing

  %  showprogress(i/length(rev));
end

for i = 1:length(genes)
    if ~isempty(genes{i}) %length(genes{i}) ~= 0
        newmodel.genes(nGenes,1) = genes(i,1);
    end
end

%Set/Reset Gene Names etc.
newmodel.genes = genes;
if isfield(newmodel,'proteins')
    if ~all(size(newmodel.proteins) == size(newmodel.genes))
        newmodel = rmfield(model,'proteins');
        if isfield(model,'proteinNames')
            newmodel = rmfield(model,'proteinNames');
        end
    end
end
modelFields = fieldnames(newmodel);
%Remove all gene Associated fields which do not fit the new gene vector
%size, these are now invalid.
for i = 1:numel(modelFields)
    if strncmp(modelFields{i},'genes',4)
        if ~all(size(newmodel.(modelFields{i})) == size(newmodel.genes))
            newmodel = rmfield(newmodel,modelFields{i});
        end
    end
end

newmodel = extendModelFieldsForType(newmodel,'rxns','originalSize',orignRxns,'targetSize',numel(newmodel.rxns));
newmodel = extendModelFieldsForType(newmodel,'mets','originalSize',orignMets,'targetSize',numel(newmodel.mets));

author: opencobra

repo: cobratoolbox


text:
function [ radius ] = atomic_radius( Z )
%ATOMIC_RADIUS Get atomic radius of element.
%   r = ATOMIC_RADIUS(Z) returns the atomic radius of element number Z in
%   Angstroms, as given in J. Chem. Phys. 47, 1300 (1967) doi:
%   10.1063/1.1712084. Z can also be a string the the element symbol. For
%   unrecognized symbols/element numbers, a a value of 1.0 is returned.
%
%   See also CHEMSYM2NUMBER, NUMBER2CHEMSYM

    if ischar(Z)
        Z = chemsym2number(Z);
    end

    default_r = 1.0;

    r = [ ...
     0.53 ... H
     0.31 ... He
     1.67 ... Li
     1.12 ... Be
     0.87 ... B
     0.67 ... C
     0.56 ... N
     0.48 ... O
     0.42 ... F
     0.38 ... Ne
     1.90 ... Na
     1.45 ... Mg
     1.18 ... Al
     1.11 ... Si
     0.98 ... P
     0.88 ... S
     0.79 ... Cl
     0.71 ... Ar
     2.43 ... K
     1.94 ... Ca
     1.84 ... Sc
     1.76 ... Ti
     1.71 ... V
     1.66 ... Cr
     1.61 ... Mn
     1.56 ... Fe
     1.52 ... Co
     1.49 ... Ni
     1.45 ... Cu
     1.42 ... Zn
     1.36 ... Ga
     1.25 ... Ge
     1.14 ... As
     1.03 ... Se
     0.94 ... Br
     0.88 ... Kr
     2.65 ... Rb
     2.19 ... Sr
     2.12 ... Y
     2.06 ... Zr
     1.98 ... Nb
     1.90 ... Mo
     1.83 ... Tc
     1.78 ... Ru
     1.73 ... Rh
     1.69 ... Pd
     1.65 ... Ag
     1.61 ... Cd
     1.56 ... In
     1.45 ... Sn
     1.33 ... Sb
     1.23 ... Te
     1.15 ... I
     1.08 ... Xe
     2.98 ... Cs
     2.53 ... Ba
     default_r ... La (missing value)
     default_r ... Ce (missing value)
     2.47 ... Pr
     2.06 ... Nd
     2.05 ... Pm
     2.38 ... Sm
     2.31 ... Eu
     2.33 ... Gd
     2.25 ... Tb
     2.28 ... Dy
     2.26 ... Ho
     2.26 ... Er
     2.22 ... Tm
     2.22 ... Yb
     2.17 ... Lu
%{
     ... Hf
      ... Ta
      ... W
      ... Re
      ... Os
      ... Ir
      ... Pt
      ... Au
      ... Hg
      ... Tl
      ... Pb
      ... Bi
      ... Po
      ... At
      ... Rn
      ... Fr
      ... Ra
      ... Ac
      ... Th
      ... Pa
      ... U
      ... Np
      ... Pu
      ... Am
      ... Cm
      ... Bk
      ... Cf
      ... Es
      ... Fm
      ... Md
      ... No
      ... Lr
      ... Rf
      ... Db
      ... Sg
      ... Bh
      ... Hs
      ... Mt
%}
     ];
     if Z <= numel(r)
        radius = r(Z);
     else
         radius = default_r; % default radius
     end
end

author: Sable

repo: mcbench-benchmarks


text:
%PARZENC Optimisation of the Parzen classifier
%
%  [W,H] = PARZENC(A,H)
%  [W,H] = A*PARZENC([],H)
%  [W,H] = A*PARZENC(H)
%
% INPUT
%  A    dataset
%  H    smoothing parameter (may be scalar, vector of per-class
%       parameters, or matrix with parameters for each class (rows) and
%       dimension (columns))
%
% OUTPUT
%  W    trained mapping
%  H    estimated smoothing (scalar value)
%
% DESCRIPTION
% Computation of the optimum smoothing parameter H for the Parzen
% classifier between the classes in the dataset A. The leave-one-out
% Lissack & Fu estimate is used for the classification error E. The
% final classifier is stored as a mapping in W. It may be converted
% into a classifier by W*CLASSC. PARZENC cannot be used for density
% estimation.
%
% In case smoothing H is specified, no learning is performed, just the
% discriminant W is produced for the given smoothing parameters H.
% Smoothing parameters may be scalar, vector of per-class parameters, or
% a matrix with individual smoothing for each class (rows) and feature
% directions (columns)
%
% EXAMPLES
% See PREX_DENSITY for densities and PREX_PARZEN for differences between
% PARZENC, PARZENDC and PARZENM.
%
% REFERENCES
% T. Lissack and K.S. Fu, Error estimation in pattern recognition via
% L-distance between posterior density functions, IEEE Trans. Inform.
% Theory, vol. 22, pp. 34-45, 1976.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, PARZEN_MAP, PARZENML, PARZENDC, PARZENM, CLASSC,
% PREX_PARZEN

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: parzenc.m,v 1.6 2008/07/03 09:11:44 duin Exp $

function [W,h] = parzenc(varargin)

	mapname = 'ParzenC';
  argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],[]);

  if mapping_task(argin,'definition')
    W = define_mapping(argin,'untrained',mapname);

  elseif mapping_task(argin,'training')			% Train a mapping.

    [a,h] = deal(argin{:});
    islabtype(a,'crisp','soft');
    isvaldfile(a,2,2); % at least 2 objects per class, 2 classes
    a = testdatasize(a);
    a = testdatasize(a,'objects');

    [m,k,c] = getsize(a);
    nlab = getnlab(a);

    if ~isempty(h)       % take user setting for smoothing parameter

      if size(h,1) == 1, h = repmat(h,c,1); end
      if size(h,2) == 1, h = repmat(h,1,k); end
      if any(size(h) ~= [c,k])
        error('Array with smoothing parameters has wrong size');
      end
      W = prmapping('parzen_map','trained',{a,h},getlablist(a),k,c);
      W = setname(W,'Parzen Classifier');
      return

    end

    % compute all object distances
    % make diagonal inf to exclude objects own contribution
    D = +distm(a) + diag(inf*ones(1,m));

    % find object frequencies
    if islabtype(a,'crisp')
      csize = classsizes(a);
      of = csize(nlab);
    else
      csize = sum(gettargets(a),1);
    end

    % find object weights q
    p = getprior(a);
    a = setprior(a,p);
    q = p(nlab)./csize(nlab);

    % initialise
    h = max(std(a)); % for sure a too high value
    L = -inf;
    Ln = 0;
    z = 0.1^(1/k); % initial step size

    % iterate

    iter = 0;
    prwaitbar(100,'parzenc: Optimizing smoothing parameter',m > 100);
    while abs(Ln-L) > 0.001 & z < 1

      % In L we store the best performance estimate found so far.
      % Ln is the actual performance (for the actual h)
      % If Ln > L we improve the bound L, and so we rest it.

      if Ln > L, L = Ln; end
      iter = iter+1;
      prwaitbar(100,100-100*exp(-iter/10));

      r = -0.5/(h^2);
      F = q(ones(1,m),:)'.*exp(D*r);           % density contributions
      FS = sum(F)*((m-1)/m); IFS = find(FS>0); % joint density distribution
      if islabtype(a,'crisp');
        G = sum(F .* (nlab(:,ones(1,m)) == nlab(:,ones(1,m))'));
        G = G.*(of-1)./of;                     % true-class densities
      else
        % here we are for soft labels (stored in targets)
        G = zeros(1,m);
        for j=1:c
          G = G + sum(F .* (a.targets(:,j) * a.targets(:,j)'));
          % to be corrected for bias?
        end
      end

      % performance estimate
      en = max(p)*ones(1,m);
      en(IFS) = (G(IFS))./FS(IFS);
      Ln = exp(sum(log(en))/m);

      if Ln < L            % compute next estimate
        z = sqrt(z);       % adjust stepsize up (recall: 0 < z < 1)
        h = h / z;         % if we don't improve, increase h (approach h_opt from below)
      else
        h = h * z;         % if we improve, decrease h (approach h_opt from above)
      end
    end
    prwaitbar(0);
    W = prmapping('parzen_map','trained',{a,repmat(h,c,k);},getlablist(a),k,c);
    W = setname(W,mapname);
    W = setcost(W,a);

  end

return

author: marianux

repo: ecg-kit

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


text:
function [p,msg] = spm_eeval(str,Type,n,m)
% Expression evaluation
% FORMAT [p,msg] = spm_eeval(str,Type,n,m)
% Str  - Expression to work with
%
% Type - type of evaluation
%      - 's'tring
%      - 'e'valuated string
%          - 'n'atural numbers
%          - 'w'hole numbers
%          - 'i'ntegers
%          - 'r'eals
%      - 'c'ondition indicator vector
%
% n ('e', 'c' & 'p' types)
%          - Size of matrix requred
%          - NaN for 'e' type implies no checking - returns input as evaluated
%          - length of n(:) specifies dimension - elements specify size
%          - Inf implies no restriction
%          - Scalar n expanded to [n,1] (i.e. a column vector)
%            (except 'x' contrast type when it's [n,np] for np
%          - E.g: [n,1] & [1,n] (scalar n) prompt for an n-vector,
%                         returned as column or row vector respectively
%                 [1,Inf] & [Inf,1] prompt for a single vector,
%                         returned as column or row vector respectively
%                 [n,Inf] & [Inf,n] prompts for any number of n-vectors,
%                         returned with row/column dimension n respectively.
%                 [a,b] prompts for an 2D matrix with row dimension a and
%                         column dimension b
%                 [a,Inf,b] prompt for a 3D matrix with row dimension a,
%                         page dimension b, and any column dimension.
%          - 'c' type can only deal with single vectors
%          - NaN for 'c' type treated as Inf
%          - Defaults (missing or empty) to NaN
%
% m ('n', 'w', 'n1', 'w1', 'bn1' & 'bw1' types)
%          - Maximum value (inclusive)
%
% m ('r' type)
%          - Maximum and minimum values (inclusive)
%
% m ('c' type)
%       - Number of unique conditions required by 'c' type
%          - Inf implies no restriction
%          - Defaults (missing or empty) to Inf - no restriction
%
% p     - Result
%
% msg   - Explanation of why it didn't work
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_eeval.m 4439 2011-08-25 17:47:07Z guillaume $

if nargin<4, m=[]; end
if nargin<3, n=[]; end
if nargin<2, Type='e'; end
if nargin<1, str=''; end
if isempty(str), p='!'; msg='empty input'; return, end
switch lower(Type)
case 's'
    p = str; msg = '';
case 'e'
    p = evalin('base',['[',str,']'],'''!''');
    if ischar(p)
        msg = 'evaluation error';
    else
        [p,msg] = sf_SzChk(p,n);
    end
case 'n'
    p = evalin('base',['[',str,']'],'''!''');
    if ischar(p)
        msg = 'evaluation error';
    elseif any(floor(p(:))~=p(:)|p(:)<1) || ~isreal(p)
        p='!'; msg='natural number(s) required';
    elseif ~isempty(m) && any(p(:)>m)
        p='!'; msg=['max value is ',num2str(m)];
    else
        [p,msg] = sf_SzChk(p,n);
    end
case 'w'
    p = evalin('base',['[',str,']'],'''!''');
    if ischar(p)
        msg = 'evaluation error';
    elseif any(floor(p(:))~=p(:)|p(:)<0) || ~isreal(p)
        p='!'; msg='whole number(s) required';
    elseif ~isempty(m) && any(p(:)>m)
        p='!'; msg=['max value is ',num2str(m)];
    else
        [p,msg] = sf_SzChk(p,n);
    end
case 'i'
    p = evalin('base',['[',str,']'],'''!''');
    if ischar(p)
        msg = 'evaluation error';
    elseif any(floor(p(:))~=p(:)) || ~isreal(p)
        p='!'; msg='integer(s) required';
    else
        [p,msg] = sf_SzChk(p,n);
    end
case 'p'
    p = evalin('base',['[',str,']'],'''!''');
    if ischar(p)
        msg = 'evaluation error';
    elseif length(setxor(p(:)',m))
        p='!'; msg='invalid permutation';
    else
        [p,msg] = sf_SzChk(p,n);
    end
case 'r'
    p = evalin('base',['[',str,']'],'''!''');
    if ischar(p)
        msg = 'evaluation error';
    elseif ~isreal(p)
        p='!'; msg='real number(s) required';
    elseif ~isempty(m) && ( max(p)>max(m) || min(p)<min(m) )
        p='!'; msg=sprintf('real(s) in [%g,%g] required',min(m),max(m));
    else
        [p,msg] = sf_SzChk(p,n);
    end
case 'c'
    if isempty(m), m=Inf; end
    [p,msg] = icond(str,n,m);

otherwise
    error('unrecognised type');
end
return;

function [i,msg] = icond(i,n,m)
if nargin<3, m=Inf; end
if nargin<2, n=NaN; end
if any(isnan(n(:)))
    n=Inf;
elseif (length(n(:))==2 & ~any(n==1)) | length(n(:))>2
    error('condition input can only do vectors')
end
if nargin<2, i=''; end
if isempty(i), varargout={[],'empty input'}; return, end
msg = ''; i=i(:)';

if ischar(i)
    if i(1)=='0' & all(ismember(unique(i(:)),char(abs('0'):abs('9'))))
        %-Leading zeros in a digit list
        msg = sprintf('%s expanded',i);
        z = min(find([diff(i=='0'),1]));
        i = [zeros(1,z), icond(i(z+1:end))'];
    else
        %-Try an eval, for functions & string #s
        i = evalin('base',['[',i,']'],'i');
    end
end

if ischar(i)
    %-Evaluation error from above: see if it's an 'abab' or 'a b a b' type:
    [c,null,i] = unique(lower(i(~isspace(i))));
    if all(ismember(c,char(abs('a'):abs('z'))))
        %-Map characters a-z to 1-26, but let 'r' be zero (rest)
        tmp = c-'a'+1; tmp(tmp=='r'-'a'+1)=0;
        i   = tmp(i);
        msg = [sprintf('[%s] mapped to [',c),...
                sprintf('%d,',tmp(1:end-1)),...
                sprintf('%d',tmp(end)),']'];
    else
        i = '!'; msg = 'evaluation error';
    end
elseif ~all(floor(i(:))==i(:))
    i = '!'; msg = 'must be integers';
elseif length(i)==1 & prod(n)>1
    msg = sprintf('%d expanded',i);
    i = floor(i./10.^[floor(log10(i)+eps):-1:0]);
    i = i-[0,10*i(1:end-1)];
end

%-Check size of i & #conditions
if ~ischar(i), [i,msg] = sf_SzChk(i,n,msg); end
if ~ischar(i) & isfinite(m) & length(unique(i))~=m
    i = '!'; msg = sprintf('%d conditions required',m);
end
return;

function str = sf_SzStr(n,unused)
%=======================================================================
%-Size info string constuction
if nargin<2, l=0; else l=1; end
if nargin<1, error('insufficient arguments'); end;
if isempty(n), n=NaN; end
n=n(:); if length(n)==1, n=[n,1]; end; dn=length(n);
if any(isnan(n)) || (prod(n)==1 && dn<=2) || (dn==2 && min(n)==1 && isinf(max(n)))
    str = ''; lstr = '';
elseif dn==2 && min(n)==1
    str = sprintf('[%d]',max(n));   lstr = [str,'-vector'];
elseif dn==2 && sum(isinf(n))==1
    str = sprintf('[%d]',min(n));   lstr = [str,'-vector(s)'];
else
    str=''; for i = 1:dn
        if isfinite(n(i)), str = sprintf('%s,%d',str,n(i));
        else str = sprintf('%s,*',str); end
    end
    str = ['[',str(2:end),']'];     lstr = [str,'-matrix'];
end
if l, str=sprintf('\t%s',lstr); else str=[str,' ']; end

function [p,msg] = sf_SzChk(p,n,msg)
%=======================================================================
%-Size checking
if nargin<3, msg=''; end
if nargin<2, n=[]; end; if isempty(n), n=NaN; else n=n(:)'; end
if nargin<1, error('insufficient arguments'), end

if ischar(p) || any(isnan(n(:))), return, end
if length(n)==1, n=[n,1]; end

dn = length(n);
sp = size(p);

if dn==2 && min(n)==1
    %-[1,1], [1,n], [n,1], [1,Inf], [Inf,1] - vector - allow transpose
    %---------------------------------------------------------------
    i = min(find(n==max(n)));
    if n(i)==1 && max(sp)>1
        p='!'; msg='scalar required';
    elseif ndims(p)~=2 || ~any(sp==1) || ( isfinite(n(i)) && max(sp)~=n(i) )
        %-error: Not2D | not vector | not right length
        if isfinite(n(i)), str=sprintf('%d-',n(i)); else str=''; end
        p='!'; msg=[str,'vector required'];
    elseif sp(i)==1 && n(i)~=1
        p=p'; msg=[msg,' (input transposed)'];
    end

elseif dn==2 && sum(isinf(n))==1
    %-[n,Inf], [Inf,n] - n vector(s) required - allow transposing
    %---------------------------------------------------------------
    i = find(isfinite(n));
    if ndims(p)~=2 || ~any(sp==n(i))
        p='!'; msg=sprintf('%d-vector(s) required',min(n));
    elseif sp(i)~=n
        p=p'; msg=[msg,' (input transposed)'];
    end

else
    %-multi-dimensional matrix required - check dimensions
    %---------------------------------------------------------------
    if ndims(p)~=dn || ~all( size(p)==n | isinf(n) )
        p = '!'; msg='';
        for i = 1:dn
            if isfinite(n(i)), msg = sprintf('%s,%d',msg,n(i));
            else msg = sprintf('%s,*',msg); end
        end
        msg = ['[',msg(2:end),']-matrix required'];
    end
end

author: spm


text:
function [QAM_BER, QAM_SER] = QAM_Simulate(SNRs, varargin)

% Run M_QAM_Model.mdl to generate Monte Carlo simulation results for
% QAM signals over AWGN channels

if nargin>1
    simLines = varargin{1};
end

open_system('M_QAM_Model')
maxNumBits = 1e7;
maxNumErrs = 100;
Ts = 1e-6;

% SNRs = -4:28;
QAM_BER = zeros(9,length(SNRs));
QAM_SER = zeros(9,length(SNRs));

S = simset('SrcWorkspace','current', 'DstWorkspace','current');
k=1;
for M = [4, 8, 16, 32, 64, 128, 256, 512, 1024]
    for EbNo = SNRs
        % Don't try to simulate BER < 1e-5 (too long!)
        tBER = berawgn(EbNo,'qam',M);
        if (tBER>1e-4)
            fprintf('Simulating %i-QAM, %idB, ', M, EbNo)
            sim('M_QAM_Model',[], S)
            QAM_BER(k,EbNo+5) = BER(1);
            QAM_SER(k,EbNo+5) = SER(1);
            fprintf('BER=%e, SER=%e\n', BER(1), SER(1))
            % Add the result to the plot if it exists
            if nargin>1
                set(simLines(k),'YData',QAM_BER(k,:))
            end
        end
    end
    drawnow
    k=k+1;
end
close_system('M_QAM_Model')


text:
function Cbn = att2Cbn(att);
%--------------------------------------------------------------------------
% calculate the direct cosine matrix c_nb from b => n-frame given att (h,p,r)
% assume 1 row or col input
% c_nb = r3(-h)*r2(-p)*r1(-r);
% Author:   Yudan Yi
%           Aug. 2004
%           Feb. 2005
%--------------------------------------------------------------------------
if (nargin<1) error('error in data'); end;
att = att(:);
if (length(att)<3) error('error in data'); end;
%H = att(1); P = att(2); R = att(3);
R = att(1); P = att(2); H = att(3);
Cbn(1,1)	= cos(H)*cos(P);
Cbn(2,1)	= sin(H)*cos(P);
Cbn(3,1)	=-	     sin(P);
Cbn(1,2)	=-sin(H)*cos(R)+cos(H)*sin(P)*sin(R);
Cbn(2,2)	= cos(H)*cos(R)+sin(H)*sin(P)*sin(R);
Cbn(3,2)	=                      cos(P)*sin(R);
Cbn(1,3)	= sin(H)*sin(R)+cos(H)*sin(P)*cos(R);
Cbn(2,3)	=-cos(H)*sin(R)+sin(H)*sin(P)*cos(R);
Cbn(3,3)	=                      cos(P)*cos(R);
%--------------------------------------------------------------------------

repo: ekfmonoslam


text:
function methods = get_method_configs()
% Get properties of all proposal methods. For examples, this includes paths to
% candidate files, intermediate files, and also information about how to sort
% candidates.

% If you want to add your own method, add it at the bottom.

colormap = [
228, 229, 97
163, 163, 163
218, 71, 56
219, 135, 45
145, 92, 146
83, 136, 173
106,61,154
225, 119, 174
142, 195, 129
51,160,44
223, 200, 51
92, 172, 158
177,89,40
177,89,40
188, 128, 189
177,89,40
251,154,153
31,120,180
177,89,40
177,89,40
177,89,40
177,89,40
177,89,40
177,89,40
177,89,40
] ./ 256;

  config = get_config();

  methods = [];

  i = numel(methods) + 1;
  methods(i).name = 'Objectness';
  methods(i).short_name = 'O';
  methods = add_paths(methods(i), config, 'objectness/');
  methods(i).order = 'descend';
  methods(i).extract = @run_objectness;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'Rahtu';
  methods(i).short_name = 'R1';
  methods(i) = add_paths(methods(i), config, 'categ_indep_detection/');
  methods(i).order = 'none';
  methods(i).extract = @run_categ_independent;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'SelectiveSearch';
  methods(i).short_name = 'SS';
  methods(i) = add_paths(methods(i), config, 'selective_search/');
  methods(i).order = 'ascend';
  methods(i).extract = @run_selective_search;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'RandomizedPrims';
  methods(i).short_name = 'RP';
  methods(i) = add_paths(methods(i), config, 'randomized_prims/');
  methods(i).gt_recall_num_candidates = 20000;
  methods(i).rerun_num_candidates = [];
  methods(i).repeatability_num_candidates = 5000;
  methods(i).order = 'none';
  methods(i).extract = @run_randomized_prims;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'Uniform';
  methods(i).short_name = 'U';
  methods(i) = add_paths(methods(i), config, 'random_uniform/');
  methods(i).order = 'random';
  methods(i).extract = @sample_uniform;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = true;

  i = numel(methods) + 1;
  methods(i).name = 'Gaussian';
  methods(i).short_name = 'G';
  methods(i) = add_paths(methods(i), config, 'random_gaussian/');
  methods(i).order = 'random';
  methods(i).extract = @sample_gaussian;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = true;

  i = numel(methods) + 1;
  methods(i).name = 'CPMC';
  methods(i).short_name = 'C';
  methods(i) = add_paths(methods(i), config, 'CPMC/');
  methods(i).order = 'descend';
  methods(i).extract = @run_cpmc;
  methods(i).num_candidates = false;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'Bing';
  methods(i).short_name = 'B';
  methods(i) = add_paths(methods(i), config, 'BING/');
  methods(i).order = 'descend';
  methods(i).extract = @run_bing;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'Endres';
  methods(i).short_name = 'E';
  methods(i) = add_paths(methods(i), config, 'endres/');
  methods(i).order = 'none';
  methods(i).extract = @run_prop;
  methods(i).num_candidates = false;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'Rantalankila';
  methods(i).short_name = 'R4';
  methods(i) = add_paths(methods(i), config, 'Rantalankila/');
  methods(i).order = 'none';
  methods(i).extract = @run_rantalankila;
  methods(i).num_candidates = true;
  methods(i).rerun_num_candidates = [100 1000 10000];
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'Superpixels';
  methods(i).short_name = 'SP';
  methods(i) = add_paths(methods(i), config, 'segmentation_baseline/');
  methods(i).order = 'none';
  methods(i).extract = @run_felsen_candidates;
  methods(i).num_candidates = false;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = true;

  i = numel(methods) + 1;
  methods(i).name = 'Sliding window';
  methods(i).short_name = 'SW';
  methods(i) = add_paths(methods(i), config, 'sliding_window/');
  methods(i).order = 'none';
  methods(i).extract = @sample_bing_windows;
  methods(i).num_candidates = true;
  methods(i).rerun_num_candidates = ceil(10 .^ (2:0.5:4));
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = true;

  i = numel(methods) + 1;
  methods(i).name = 'EdgeBoxes';
  methods(i).short_name = 'EB70';
  methods(i) = add_paths(methods(i), config, 'edge_boxes_70/');
  methods(i).order = 'descend';
  methods(i).extract = @run_edge_boxes;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'EdgeBoxes50';
  methods(i).short_name = 'EB50';
  methods(i) = add_paths(methods(i), config, 'edge_boxes_50/');
  methods(i).order = 'descend';
  methods(i).extract = @run_edge_boxes50;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;
  methods(i).line_style = ':';

  i = numel(methods) + 1;
  methods(i).name = 'MCG';
  methods(i).short_name = 'M';
  methods(i) = add_paths(methods(i), config, 'MCG2015/');
  methods(i).order = 'none';
  methods(i).extract = @run_MCG2015;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'EdgeBoxes90';
  methods(i).short_name = 'EB90';
  methods(i) = add_paths(methods(i), config, 'edge_boxes_90/');
  methods(i).order = 'descend';
  methods(i).extract = @run_edge_boxes90;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;
  methods(i).line_style = '-.';

  i = numel(methods) + 1;
  methods(i).name = 'Rigor';
  methods(i).short_name = 'Ri';
  methods(i) = add_paths(methods(i), config, 'rigor/');
  methods(i).order = 'none';
  methods(i).extract = @run_rigor;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;
  methods(i).rerun_num_candidates = [10 100 1000 2000 10000];

  i = numel(methods) + 1;
  methods(i).name = 'Geodesic';
  methods(i).short_name = 'G';
  methods(i) = add_paths(methods(i), config, 'geodesic/');
  methods(i).order = 'none';
  methods(i).extract = @run_gop;
  methods(i).num_candidates = false;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;
  methods(i).rerun_num_candidates = [10 100 1000 2000 10000];

  i = numel(methods) + 1;
  methods(i).name = 'EdgeBoxes60';
  methods(i).short_name = 'EB60';
  methods(i) = add_paths(methods(i), config, 'edge_boxes_60/');
  methods(i).order = 'descend';
  methods(i).extract = @run_edge_boxes60;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'EdgeBoxes80';
  methods(i).short_name = 'EB80';
  methods(i) = add_paths(methods(i), config, 'edge_boxes_80/');
  methods(i).order = 'descend';
  methods(i).extract = @run_edge_boxes80;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'EdgeBoxes55';
  methods(i).short_name = 'EB55';
  methods(i) = add_paths(methods(i), config, 'edge_boxes_55/');
  methods(i).order = 'descend';
  methods(i).extract = @(im,num) run_edge_boxes(im, num, 0.65, 0.60);
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'EdgeBoxes65';
  methods(i).short_name = 'EB65';
  methods(i) = add_paths(methods(i), config, 'edge_boxes_65/');
  methods(i).order = 'descend';
  methods(i).extract = @(im,num) run_edge_boxes(im, num, 0.65, 0.70);
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'EdgeBoxes75';
  methods(i).short_name = 'EB75';
  methods(i) = add_paths(methods(i), config, 'edge_boxes_75/');
  methods(i).order = 'descend';
  methods(i).extract = @(im,num) run_edge_boxes(im, num, 0.70, 0.80);
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'EdgeBoxes85';
  methods(i).short_name = 'EB85';
  methods(i) = add_paths(methods(i), config, 'edge_boxes_85/');
  methods(i).order = 'descend';
  methods(i).extract = @(im,num) run_edge_boxes(im, num, 0.80, 0.90);
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'Edge Boxes AR';
  methods(i).short_name = 'EBAR';
  methods(i) = add_paths(methods(i), config, 'edge_boxes_AR/');
  methods(i).order = 'descend';
  methods(i).extract = @run_edge_boxes_AR;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  % add your own method here:
  if false
  i = numel(methods) + 1;
  methods(i).name = 'The full name of your method';
  methods(i).short_name = 'a very short version of the name';
  methods(i) = add_paths(methods(i), config, 'ours-wip/');
  % This specifies how to order candidates so that the first n, are the best n
  % candidates. For example we run a method for 10000 candidates and then take
  % the first 10, instead of rerunning for 10 candidates. Valid orderings are:
  %   none: candidates are already sorted, do nothing
  %   ascend/descend: sort by score descending or ascending
  %   random: random order
  %   biggest/smallest: sort by size of the bounding boxes
  methods(i).order = 'descend';
  % A function pointer to a method that runs your proposal detector.
  methods(i).extract = @run_edge_boxes90;
  % If your method supports sorting this should be empty. If your method has to
  % be rerun for every number of candidates we want, specify the number of
  % candidates here:
  methods(i).rerun_num_candidates = []; % ceil(10 .^ (2:0.5:4));
  % Specifies whether or not your method takes the desired number of candidates
  % as an input.
  % TODO(hosang): Is this actually used anywhere?
  methods(i).num_candidates = true;
  % color for drawing
  methods(i).color = colormap(i,:);
  % This should be false. Is used for drawing baselines dashed.
  methods(i).is_baseline = false;
  end

  % do the sorting dance
  sort_keys = [num2cell([methods.is_baseline])', {methods.name}'];
  for i = 1:numel(methods)
    sort_keys{i,1} = sprintf('%d', sort_keys{i,1});
  end
  [~,idx] = sortrows(sort_keys);
  for i = 1:numel(methods)
    methods(idx(i)).sort_key = i;
  end
end

function method = add_paths(method, config, method_subdir)
  method.candidate_dir = fullfile(config.precomputed_candidates, method_subdir, 'mat');
  method.repeatability_candidate_dir = fullfile(config.precomputed_candidates, method_subdir, 'repeatability_mat');
  method.best_voc07_candidates_file = fullfile(config.precomputed_candidates, method_subdir, 'best_candidates.mat');
  method.best_imagenet_candidates_file = fullfile(config.precomputed_candidates, method_subdir, 'best_candidates_imagenet.mat');
  method.repeatability_matching_file = fullfile(config.precomputed_candidates, method_subdir, 'repeatability_matching.mat');
  method.candidate_dir_coco = fullfile(config.precomputed_candidates_coco, method_subdir, 'mat');
  method.best_coco14_candidates_file = fullfile(config.precomputed_candidates_coco, method_subdir, 'best_candidates_coco14.mat');
end

author: hosang

repo: detection-proposals


text:
function [ymu,ys2,fmu,fs2,dymu,dys2,dfmu,dfs2] = gpgrad(x0,gpstruct,method,dx)
%GPGRAD Gradient of Gaussian process mean prediction

if nargin < 3 || isempty(method); method = 'central'; end
if nargin < 4 || isempty(dx); dx = 1e-6; end

D = size(x0,2);
Nhyp = length(gpstruct.hyp);

% Initialize variables
ymu = zeros(Nhyp,1);
fmu = zeros(Nhyp,1);
ys2 = zeros(Nhyp,1);
fs2 = zeros(Nhyp,1);

if nargout > 4
    dymu = zeros(Nhyp,D);
    dfmu = zeros(Nhyp,D);
    dys2 = zeros(Nhyp,D);
    dfs2 = zeros(Nhyp,D);
end

% Compute predictions and derivatives for each hyperparameter sample
for i = 1:Nhyp

    try
        if nargout > 4
            [dy,y0] = fgrad(@(xi_) gpfunc(xi_,gpstruct,gpstruct.hyp(i)),x0,method,'Vectorized','Step',dx);
        else
            y0 = gpfunc(x0,gpstruct,gpstruct.hyp(i));
        end
    catch
        y0 = NaN(1,4);
        if nargout > 4; dy = NaN(4,D); end
        warning('bads:gpGradFail', 'Error in computing the GP derivative.');
    end

    ymu(i) = y0(1);
    ys2(i) = y0(2);
    fmu(i) = y0(3);
    fs2(i) = y0(4);

    if nargout > 4
        dymu(i,:) = dy(1,:);
        dfmu(i,:) = dy(2,:);
        dys2(i,:) = dy(3,:);
        dfs2(i,:) = dy(4,:);
    end

end

%--------------------------------------------------------------------------

function yy = gpfunc(xi,gpstruct,hyp)
    if gpstruct.gpmlext
        % Pre-computed Laplace approximation at MAP
        if isfield(gpstruct,'hypSigma') && ~isempty(gpstruct.hypSigma)
            hyp.Sigma = gpstruct.hypSigma;
        end

        if isfield(gpstruct,'post') && ~isempty(gpstruct.post)
        %    [ymuib,ys2ib,fmuib,fs2ib] = mygp(hyp,gpstruct.inf,gpstruct.mean,gpstruct.cov, ...
        %gpstruct.lik,gpstruct.x,gpstruct.post,xi);

            [ymui,ys2i,fmui,fs2i] = mgp(hyp,gpstruct.inf,gpstruct.mean,gpstruct.cov, ...
        gpstruct.lik,gpstruct.x,gpstruct.post,xi);

            %[fmuib(:) - fmui(:), sqrt(fs2ib(:)) - sqrt(fs2i(:))]

        else
            [ymui,ys2i,fmui,fs2i] = mgp(hyp,gpstruct.inf,gpstruct.mean,gpstruct.cov, ...
        gpstruct.lik,gpstruct.x,gpstruct.y,xi);
        end
        yy = [ymui,ys2i,fmui,fs2i];
    else

        if isfield(gpstruct,'post') && ~isempty(gpstruct.post)
            [ymui,ys2i,fmui,fs2i] = mygp(hyp,gpstruct.inf,gpstruct.mean,gpstruct.cov, ...
        gpstruct.lik,gpstruct.x,gpstruct.post,gpstruct.s,xi);
        else
            [ymui,ys2i,fmui,fs2i] = mygp(hyp,gpstruct.inf,gpstruct.mean,gpstruct.cov, ...
        gpstruct.lik,gpstruct.x,gpstruct.y,gpstruct.s,xi);
        end

        if 1
            % [~, ymui, ys2i] = feval(gpstruct.lik{:}, hyp.lik, [], fmui, fs2i);

            lik = hyp.lik; lik(end) = -Inf;
            % fmuiold = fmui; fs2iold = fs2i;
            [~, fmui, fs2i] = feval(gpstruct.lik{:}, lik, [], [], fmui, fs2i);
        end

        yy = [ymui,ys2i,fmui,fs2i];
        % [ymui,sqrt(ys2i),fmui,sqrt(fs2i),fmuiold,sqrt(fs2iold)]
    end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [cs,h]=m_contour(long,lat,data,varargin)
%  M_CONTOUR Draws contour lines on a map
%    M_CONTOUR(LONG,LAT,DATA,...) draw contours on a map. Behavior
%    is the same as for CONTOUR except that LONG and LAT vectors or
%    matrices must be specified.
%
%    [CS,H]=M_CONTOUR(...) returns a contour matrix C and a vector
%    H of handles to LINE or PATCH objects for use by CLABEL.
%
%    See also CONTOUR

% Rich Pawlowicz (rich@ocgy.ubc.ca) 17/Jan/1998
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

% 9/Dec/98 - made sure bad things don't happen if all your lat/long
%            points are out of the plot region.
% 6/Nov/00 - eliminate returned stuff if ';' neglected (thx to D Byrne)

global MAP_PROJECTION

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

if min(size(long))==1 && min(size(lat))==1
 [long,lat]=meshgrid(long,lat);
end

[X,Y]=m_ll2xy(long,lat,'clip','on');

i=isnan(X);      % For these we set the *data* to NaN...
data(i)=NaN;

                 % And then recompute positions without clipping. THis
                 % is necessary otherwise contouring fails (X/Y with NaN
                 % is a no-no.
if any(i(:)), [X,Y]=m_ll2xy(long,lat,'clip','off'); end

if any(~i(:))
 [cs,h]=contour(X,Y,data,varargin{:});
 set(h,'tag','m_contour');
else
  cs=[];h=[];
end

if nargout==0
 clear cs h
end

author: goGPS-Project


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: rotation of an US image, based on the book's version
%
% - load data                  (setup2DUSData)
% - transform                  (book_rigid2D)
% - interpolate                (linearInter)
% - and visualize              (viewImage2D)
%==============================================================================

setup2DUSData;

c   = (omega(2:2:end)+omega(1:2:end))'/2; alpha = pi/6;
rot = [cos(alpha),-sin(alpha);sin(alpha),cos(alpha)];
wc  = [alpha;(eye(2)-rot)*c];
xc  = getCellCenteredGrid(omega,m);
yc  = rigid2D(wc,xc);
Tc  = linearInter(dataT,omega,yc);
FAIRfigure(2); viewImage2D(Tc,omega,m,'colormap','gray(256)');
%==============================================================================


text:
% DETECTOR
% See also
%
% Please see acfReadme for an overview of the detection code.
%
% Aggregate channel features object detector:
%   acfDemoCal   - Demo for aggregate channel features object detector on Caltech dataset.
%   acfDemoInria - Demo for aggregate channel features object detector on Inria dataset.
%   acfDetect    - Run aggregate channel features object detector on given image(s).
%   acfModify    - Modify aggregate channel features object detector.
%   acfReadme    - Aggregate Channel Features Detector Overview.
%   acfSweeps    - Parameter sweeps for ACF pedestrian detector.
%   acfTest      - Test aggregate channel features object detector given ground truth.
%   acfTrain     - Train aggregate channel features object detector.
%
% Object bounding box utilities and labeling tools:
%   bbApply      - Functions for manipulating bounding boxes (bb).
%   bbGt         - Bounding box (bb) annotations struct, evaluation and sampling routines.
%   bbLabeler    - Bounding box or ellipse labeler for static images.
%   bbNms        - Bounding box (bb) non-maximal suppression (nms).

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


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

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

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

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

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

end

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

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

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

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

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

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

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

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

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

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

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

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

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

end

copyfile(pathGamsFunction);

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

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

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

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

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

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

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

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

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

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

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

            %initialize empty array for storing
            n_mustSet = size(m1_full,1);
            mustLL = cell(n_mustSet, 2);
            pos_mustLL = zeros(size(mustLL));
            mustLL_linear = {};

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

            %initialize arrays to be returned by this function
            mustLL = {};
            pos_mustLL = [];
            mustLL_linear = {};
            pos_mustLL_linear = [];
        end

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

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

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

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

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

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

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

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

    end

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

end

end


text:
function [v,faces] = generateUnitCells(xy,unitCell,varargin)
% generate a list of patches according to spatial coordinates and the unitCell
%
% Input
%  xy       - midpoints of the cells
%  unitCell - spatial coordinates of the unit cell
%
% Ouput
%  v     - list of vertices
%  faces - list of faces

% compute the vertices
x = reshape(bsxfun(@plus,xy(:,1),unitCell(:,1).'),[],1);
y = reshape(bsxfun(@plus,xy(:,2),unitCell(:,2).'),[],1);

% remove equal points
% in general every measurment point generates 4 or 6 vertex points
% some of them apear multiple times
% lets try to reduce them

if ~check_option(varargin,'noStripes')
  eps = abs(diff(unitCell));
  eps = eps(eps > max(eps(:))/10);
  eps = min(eps) / 12;
  [~,m,n] = unique(round([x-min(x) y-min(y)]./eps),'rows');
else
  [~,m,n] = uniquetol([x-min(x) y-min(y)],0.01/sqrt(size(xy,1)),'ByRows',true );
end

v = [x(m) y(m)];

% set faces
faces = reshape(n, [], size(unitCell,1));

% tic
% [~,m,n] = unique(round([x-min(x) y-min(y)]./eps),'rows');
% toc

%%
% tic
% [~,m,n] = uniquetol([x-min(x) y-min(y)],0.01/sqrt(size(xy,1)),'ByRows',true );
% toc

author: mtex-toolbox

repo: mtex

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


text:
function [xd,xn,option] = denoise(x,h,type,option)
%    [xd,xn,option] = denoise(x,h,type,option);
%
%    DENOISE is a generic program for wavelet based denoising.
%    The program will denoise the signal x using the 2-band wavelet
%    system described by the filter h using either the traditional
%    discrete wavelet transform (DWT) or the linear shift invariant
%    discrete wavelet transform (also known as the undecimated DWT
%    (UDWT)).
%
%    Input:
%       x         : 1D or 2D signal to be denoised
%       h         : Scaling filter to be applied
%       type      : Type of transform (Default: type = 0)
%                   0 --> Discrete wavelet transform (DWT)
%                   1 --> Undecimated DWT (UDWT)
%       option    : Default settings is marked with '*':
%                   *type = 0 --> option = [0 3.0 0 0 0 0]
%                   type = 1 --> option = [0 3.6 0 1 0 0]
%       option(1) : Whether to threshold low-pass part
%                   0 --> Don't threshold low pass component
%                   1 --> Threshold low pass component
%       option(2) : Threshold multiplier, c. The threshold is
%                   computed as:
%                     thld = c*MAD(noise_estimate)).
%                   The default values are:
%                     c = 3.0 for the DWT based denoising
%                     c = 3.6 for the UDWT based denoising
%       option(3) : Type of variance estimator
%                   0 --> MAD (mean absolute deviation)
%                   1 --> STD (classical numerical std estimate)
%       option(4) : Type of thresholding
%                   2 --> Soft thresholding
%                   1 --> Hard thresholding
%       option(5) : Number of levels, L, in wavelet decomposition. By
%                   setting this to the default value '0' a maximal
%                   decomposition is used.
%       option(6) : Actual threshold to use (setting this to
%                   anything but 0 will mean that option(3)
%                   is ignored)
%
%    Output:
%       xd     : Estimate of noise free signal
%       xn     : The estimated noise signal (x-xd)
%       option : A vector of actual parameters used by the
%                program. The vector is configured the same way as
%                the input option vector with one added element
%                option(7) = type.
%
%  HERE'S AN EASY WAY TO RUN THE EXAMPLES:
%  Cut-and-paste the example you want to run to a new file
%  called ex.m, for example. Delete out the % at the beginning
%  of each line in ex.m (Can use search-and-replace in your editor
%  to replace it with a space). Type 'ex' in matlab and hit return.
%
%    Example 1:
%       h = daubcqf(6); [s,N] = makesig('Doppler'); n = randn(1,N);
%       x = s + n/10;     % (approximately 10dB SNR)
%       figure;plot(x);hold on;plot(s,'r');
%
%       %Denoise x with the default method based on the DWT
%       [xd,xn,opt1] = denoise(x,h);
%       figure;plot(xd);hold on;plot(s,'r');
%
%       %Denoise x using the undecimated (LSI) wavelet transform
%       [yd,yn,opt2] = denoise(x,h,1);
%       figure;plot(yd);hold on;plot(s,'r');
%
% Example 2: (on an image)
%      h = daubcqf(6);  load lena;
%      noisyLena = lena + 25 * randn(size(lena));
%      figure; colormap(gray); imagesc(lena); title('Original Image');
%       figure; colormap(gray); imagesc(noisyLena); title('Noisy Image');
%       Denoise lena with the default method based on the DWT
%      [denoisedLena,xn,opt1] = denoise(noisyLena,h);
%      figure; colormap(gray); imagesc(denoisedLena); title('denoised Image');
%
%
%    See also: mdwt, midwt, mrdwt, mirdwt, SoftTh, HardTh, setopt
%
%Author: Jan Erik Odegard  <odegard@ece.rice.edu>

if(nargin < 2)
  error('You need to provide at least 2 inputs: x and h');
end;
if(nargin < 3),
  type = 0;
  option = [];
elseif(nargin < 4)
  option = [];
end;
if(isempty(type)),
  type = 0;
end;
if(type == 0),
  default_opt = [0 3.0 0 2 0 0];
elseif(type == 1),
  default_opt = [0 3.6 0 1 0 0];
else
  error(['Unknown denoising method',10,...
	  'If it is any good we need to have a serious talk :-)']);
end;
option = setopt(option,default_opt);
[mx,nx] = size(x);
dim = min(mx,nx);
if(dim == 1),
  n = max(mx,nx);
else
  n = dim;
end;
if(option(5) == 0),
  L = floor(log2(n));
else
  L = option(5);
end;
if(type == 0), 			% Denoising by DWT
  xd = mdwt(x,h,L);
  if (option(6) == 0),
    tmp = xd(floor(mx/2)+1:mx,floor(nx/2)+1:nx);
    if(option(3) == 0),
      thld = option(2)*median(abs(tmp(:)))/.67;
    elseif(option(3) == 1),
      thld = option(2)*std(tmp(:));
    else
      error('Unknown threshold estimator, Use either MAD or STD');
    end;
  else
    thld = option(6);
  end;
  if(dim == 1)
    ix = 1:n/(2^L);
    ykeep = xd(ix);
  else
    ix = 1:mx/(2^L);
    jx = 1:nx/(2^L);
    ykeep = xd(ix,jx);
  end;
  if(option(4) == 2),
    xd = SoftTh(xd,thld);
  elseif(option(4) == 1),
    xd = HardTh(xd,thld);
  else
    error('Unknown threshold rule. Use either Soft (2) or Hard (1)');
  end;
  if (option(1) == 0),
    if(dim == 1),
      xd(ix) = ykeep;
    else
      xd(ix,jx) = ykeep;
    end;
  end;
  xd = midwt(xd,h,L);
elseif(type == 1), 			% Denoising by UDWT
  [xl,xh] = mrdwt(x,h,L);
  if(dim == 1),
    c_offset = 1;
  else
    c_offset = 2*nx + 1;
  end;
  if (option(6) == 0),
    tmp = xh(:,c_offset:c_offset+nx-1);
    if(option(3) == 0),
      thld = option(2)*median(abs(tmp(:)))/.67;
    elseif(option(3) == 1),
      thld = option(2)*std(tmp(:));
    else
      error('Unknown threshold estimator, Use either MAD or STD');
    end;
  else
    thld = option(6);
  end;
  if(option(4) == 2),
    xh = SoftTh(xh,thld);
    if(option(1) == 1),
      xl = SoftTh(xl,thld);
    end;
  elseif(option(4) == 1),
    xh = HardTh(xh,thld);
    if(option(1) == 1),
      xl = HardTh(xl,thld);
    end;
  else
    error('Unknown threshold rule. Use either Soft (2) or Hard (1)');
  end;
  xd = mirdwt(xl,xh,h,L);
else 					% Denoising by unknown method
  error(['Unknown denoising method',10,...
         'If it is any good we need to have a serious talk :-)']);
end;
option(6) = thld;
option(7) = type;
xn = x - xd;

author: ricedsp


text:
function batch_ds_info=my_get_batch_ds_info(...
    opts, imdb, work_info, net_config, work_info_epoch, work_info_batch)

task_idx=work_info_batch.ref.task_idxes;
task_subidx=work_info_batch.ref.task_subidxes;

% only support using one image for one batch...
batch_img_num=length(task_subidx);
assert(batch_img_num==1);
assert(work_info_batch.ref.batch_task_num==1);

batch_data=get_batch_data(opts, imdb, task_idx);
img_size_input=size(batch_data.img_data);
img_size_input=img_size_input(1:2);
assert(size(batch_data.img_data, 4)==batch_img_num);

% do augmentation before crop

img_size_aug=[];

if work_info.ref.run_trn

    batch_data=batch_do_data_augmentation(opts, imdb, work_info_epoch, batch_data, work_info_batch, work_info);
    img_size_aug=size(batch_data.img_data);
    img_size_aug=img_size_aug(1:2);

    batch_data=batch_do_data_crop(opts, imdb, work_info_epoch, batch_data, work_info_batch, work_info);

else

    batch_data.aug_info=[];
    batch_data.crop_info=[];
end

net_input_img_scales=net_config.ref.net_input_img_scales;

[net_input_info, net_input_str]=gen_net_input_info(opts, imdb, batch_data, net_input_img_scales);
img_size=size(batch_data.img_data);
img_size=img_size(1:2);

batch_ds_info=[];
batch_ds_info.aug_info=batch_data.aug_info;
batch_ds_info.crop_info=batch_data.crop_info;
batch_ds_info.task_idx=task_idx;
batch_ds_info.img_size=img_size;
batch_ds_info.img_size_aug=img_size_aug;
batch_ds_info.img_size_input=img_size_input;
batch_ds_info.net_input_str=net_input_str;

batch_ds_info.net_input_info=net_input_info;
batch_ds_info.batch_data=batch_data;

end

function img_data=pre_process_img_net_input(imdb, img_data)

    data_norm_info=imdb.ref.data_norm_info;

    img_data=single(img_data);
    img_data=img_data-data_norm_info.image_mean;

end

function [net_input_info,net_input_str]=gen_net_input_info(opts, imdb, batch_data, net_input_img_scales)

feat_scale_num=length(net_input_img_scales);

org_img_size=size(batch_data.img_data);
org_img_size=org_img_size(1:2);

input_info_child_groups=cell(feat_scale_num, 1);
for s_idx=1:length(net_input_img_scales)
    one_scale=net_input_img_scales(s_idx);

    if one_scale==1
        one_img_data_scaled=batch_data.img_data;
    else
        scaled_img_size=round(org_img_size.*one_scale);
        one_img_data_scaled=imresize(batch_data.img_data, scaled_img_size);
    end

    one_child_input_info=[];
    one_child_input_info.is_group_data=false;
    one_child_input_info.x=one_img_data_scaled;
    input_info_child_groups{s_idx}=one_child_input_info;
end

for tmp_img_idx=1:length(input_info_child_groups)

    one_img=input_info_child_groups{tmp_img_idx}.x;
    assert(size(one_img, 4)==1);
    assert(isa(one_img, 'uint8'));
    one_img=pre_process_img_net_input(imdb, one_img);
    assert(isa(one_img, 'single'));
    input_info_child_groups{tmp_img_idx}.x=one_img;
end

net_input_info=[];
net_input_info.is_group_data=true;
net_input_info.data_child_groups=input_info_child_groups;

net_input_str=[];
for s_idx=1:feat_scale_num
    one_net_input_size=size(input_info_child_groups{s_idx}.x);
    if length(one_net_input_size)<4
        one_net_input_size(4)=1;
    end
    one_net_input_str=my_gen_array_str(one_net_input_size);
    net_input_str=cat(2, net_input_str, one_net_input_str);
end

net_input_info=my_init_input_info(net_input_info);

end

function batch_data=get_batch_data(train_opts, imdb, img_idx)

    batch_data=[];
    if train_opts.cache_data_mem

        if ~isfield(imdb.ref, 'data_cache_tasks') || isempty(imdb.ref.data_cache_tasks)
            imdb.ref.data_cache_tasks=cell(imdb.ref.max_task_idx, 1);
        end

        batch_data=imdb.ref.data_cache_tasks{img_idx};
        if isempty(batch_data)
            do_ds_cache_all_before_run(train_opts, imdb);
            batch_data=imdb.ref.data_cache_tasks{img_idx};
            assert(~isempty(batch_data));
        end
    end

    if isempty(batch_data)
        batch_data=do_load_batch_data(train_opts, imdb, img_idx);
    end

end

function do_ds_cache_all_before_run(train_opts, imdb)

    disp('cache all data in memory ...');
    task_idxes=imdb.ref.task_idxes;
    task_num=length(task_idxes);
    for t_idx=1:task_num
        one_task_idx=task_idxes(t_idx);
        batch_data=do_load_batch_data(train_opts, imdb, one_task_idx);
        imdb.ref.data_cache_tasks{one_task_idx}=batch_data;
        if mod(t_idx, 200)==1 || t_idx==task_num
            fprintf('caching item:%d/%d\n', t_idx, task_num)
        end
    end

end

function batch_data=do_load_batch_data(train_opts, imdb, img_idx)

    ds_info=imdb.ref.ds_info;
    img_data_batch=load_img_from_ds_info(ds_info, img_idx);
    assert(isa(img_data_batch, 'uint8'));

    img_size=[size(img_data_batch, 1), size(img_data_batch, 2)];
    img_size_input=img_size;

    if ~train_opts.use_dummy_gt
        mask_data_info=ds_info.class_idxes_mask_data_info;
        one_mask=load_mask_from_ds_info(mask_data_info, img_idx);
    else
    	one_mask=ones(img_size, 'uint8');
    end
    assert(all(size(one_mask)==img_size_input));

    todo_scale=1;
    short_edge=min(img_size);
    if short_edge<train_opts.input_img_short_edge_min
        todo_scale=train_opts.input_img_short_edge_min/short_edge;
        assert(todo_scale>=1);
    end

    if short_edge>train_opts.input_img_short_edge_max
        todo_scale=train_opts.input_img_short_edge_max/short_edge;
        assert(todo_scale<=1);
    end
    todo_scale=todo_scale*train_opts.input_img_scale;

    if todo_scale~=1
        img_data_batch=imresize(img_data_batch, todo_scale);
%         img_size=[size(img_data_batch, 1), size(img_data_batch, 2)];
    end

    batch_data=[];
    batch_data.img_data=img_data_batch;
    batch_data.label_data=one_mask;
    batch_data.img_idx=img_idx;
    batch_data.img_size_origin=img_size_input;

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
clear all; close all;
I=imread('tire.tif');
J=histeq(I);
figure;
subplot(121);
imshow(uint8(I));
subplot(122);
imshow(uint8(J));
figure;
subplot(121);
imhist(I, 64);
subplot(122);
imhist(J, 64);

author: HuangCongQing


text:
%% Machine Learning Online Class
%  Exercise 7 | Principle Component Analysis and K-Means Clustering
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     pca.m
%     projectData.m
%     recoverData.m
%     computeCentroids.m
%     findClosestCentroids.m
%     kMeansInitCentroids.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% ================= Part 1: Find Closest Centroids ====================
%  To help you implement K-Means, we have divided the learning algorithm
%  into two functions -- findClosestCentroids and computeCentroids. In this
%  part, you shoudl complete the code in the findClosestCentroids function.
%
fprintf('Finding closest centroids.\n\n');

% Load an example dataset that we will be using
load('ex7data2.mat');

% Select an initial set of centroids
K = 3; % 3 Centroids
initial_centroids = [3 3; 6 2; 8 5];

% Find the closest centroids for the examples using the
% initial_centroids
idx = findClosestCentroids(X, initial_centroids);

fprintf('Closest centroids for the first 3 examples: \n')
fprintf(' %d', idx(1:3));
fprintf('\n(the closest centroids should be 1, 3, 2 respectively)\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ===================== Part 2: Compute Means =========================
%  After implementing the closest centroids function, you should now
%  complete the computeCentroids function.
%
fprintf('\nComputing centroids means.\n\n');

%  Compute means based on the closest centroids found in the previous part.
centroids = computeCentroids(X, idx, K);

fprintf('Centroids computed after initial finding of closest centroids: \n')
fprintf(' %f %f \n' , centroids');
fprintf('\n(the centroids should be\n');
fprintf('   [ 2.428301 3.157924 ]\n');
fprintf('   [ 5.813503 2.633656 ]\n');
fprintf('   [ 7.119387 3.616684 ]\n\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =================== Part 3: K-Means Clustering ======================
%  After you have completed the two functions computeCentroids and
%  findClosestCentroids, you have all the necessary pieces to run the
%  kMeans algorithm. In this part, you will run the K-Means algorithm on
%  the example dataset we have provided.
%
fprintf('\nRunning K-Means clustering on example dataset.\n\n');

% Load an example dataset
load('ex7data2.mat');

% Settings for running K-Means
K = 3;
max_iters = 10;

% For consistency, here we set centroids to specific values
% but in practice you want to generate them automatically, such as by
% settings them to be random examples (as can be seen in
% kMeansInitCentroids).
initial_centroids = [3 3; 6 2; 8 5];

% Run K-Means algorithm. The 'true' at the end tells our function to plot
% the progress of K-Means
[centroids, idx] = runkMeans(X, initial_centroids, max_iters, true);
fprintf('\nK-Means Done.\n\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ============= Part 4: K-Means Clustering on Pixels ===============
%  In this exercise, you will use K-Means to compress an image. To do this,
%  you will first run K-Means on the colors of the pixels in the image and
%  then you will map each pixel on to it's closest centroid.
%
%  You should now complete the code in kMeansInitCentroids.m
%

fprintf('\nRunning K-Means clustering on pixels from an image.\n\n');

%  Load an image of a bird
A = double(imread('bird_small.png'));

% If imread does not work for you, you can try instead
%   load ('bird_small.mat');

A = A / 255; % Divide by 255 so that all values are in the range 0 - 1

% Size of the image
img_size = size(A);

% Reshape the image into an Nx3 matrix where N = number of pixels.
% Each row will contain the Red, Green and Blue pixel values
% This gives us our dataset matrix X that we will use K-Means on.
X = reshape(A, img_size(1) * img_size(2), 3);

% Run your K-Means algorithm on this data
% You should try different values of K and max_iters here
K = 16;
max_iters = 10;

% When using K-Means, it is important the initialize the centroids
% randomly.
% You should complete the code in kMeansInitCentroids.m before proceeding
initial_centroids = kMeansInitCentroids(X, K);

% Run K-Means
[centroids, idx] = runkMeans(X, initial_centroids, max_iters);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================= Part 5: Image Compression ======================
%  In this part of the exercise, you will use the clusters of K-Means to
%  compress an image. To do this, we first find the closest clusters for
%  each example. After that, we

fprintf('\nApplying K-Means to compress an image.\n\n');

% Find closest cluster members
idx = findClosestCentroids(X, centroids);

% Essentially, now we have represented the image X as in terms of the
% indices in idx.

% We can now recover the image from the indices (idx) by mapping each pixel
% (specified by it's index in idx) to the centroid value
X_recovered = centroids(idx,:);

% Reshape the recovered image into proper dimensions
X_recovered = reshape(X_recovered, img_size(1), img_size(2), 3);

% Display the original image
subplot(1, 2, 1);
imagesc(A);
title('Original');

% Display compressed image side by side
subplot(1, 2, 2);
imagesc(X_recovered)
title(sprintf('Compressed, with %d colors.', K));

fprintf('Program paused. Press enter to continue.\n');
pause;

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 8 Assignments/K-Means Clustering and PCA/mlclass-ex7/ex7.m