text:
function sparse_grid_laguerre_test ( )

%*****************************************************************************80
%
%% SPARSE_GRID_LAGUERRE_TEST tests the SPARSE_GRID_LAGUERRE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_LAGUERRE_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the SPARSE_GRID_LAGUERRE library.\n' );
%
%  Count number of points in sparse rule from DIM_MIN to DIM_MAX, LEVEL_MAX_MAX.
%
  dim_min = 1;
  dim_max = 5;
  level_max_min = 0;
  level_max_max = 10;

  sparse_grid_laguerre_test01 ( dim_min, dim_max, level_max_min, level_max_max );

  dim_min = 6;
  dim_max = 10;
  level_max_min = 0;
  level_max_max = 10;

  sparse_grid_laguerre_test01 ( dim_min, dim_max, level_max_min, level_max_max );

  dim_min = 100;
  dim_max = 100;
  level_max_min = 0;
  level_max_max = 2;

  sparse_grid_laguerre_test01 ( dim_min, dim_max, level_max_min, level_max_max );
%
%  Compute abstract grid indices of sparse grid points as selected from product grid
%  for DIMENSION, LEVEL_MAX.
%
  dim_num = 2;
  level_max = 3;
  sparse_grid_laguerre_test02 ( dim_num, level_max );

  dim_num = 2;
  level_max = 4;
  sparse_grid_laguerre_test02 ( dim_num, level_max );

  dim_num = 3;
  level_max = 0;
  sparse_grid_laguerre_test02 ( dim_num, level_max );

  dim_num = 3;
  level_max = 2;
  sparse_grid_laguerre_test02 ( dim_num, level_max );

  dim_num = 6;
  level_max = 2;
  sparse_grid_laguerre_test02 ( dim_num, level_max );
%
%  Compute sparse rule for DIMENSION, LEVEL_MAX.
%
  dim_num = 2;
  level_max = 0;
  sparse_grid_laguerre_test03 ( dim_num, level_max );

  dim_num = 2;
  level_max = 3;
  sparse_grid_laguerre_test03 ( dim_num, level_max );

  dim_num = 2;
  level_max = 4;
  sparse_grid_laguerre_test03 ( dim_num, level_max );

  dim_num = 3;
  level_max = 0;
  sparse_grid_laguerre_test03 ( dim_num, level_max );

  dim_num = 3;
  level_max = 2;
  sparse_grid_laguerre_test03 ( dim_num, level_max );
%
%  Test sum of weights for DIMENSION, LEVEL_MAX.
%
  dim_num = 2;
  level_max = 4;
  sparse_grid_laguerre_test04 ( dim_num, level_max );

  dim_num = 3;
  level_max = 0;
  sparse_grid_laguerre_test04 ( dim_num, level_max );

  dim_num = 3;
  level_max = 1;
  sparse_grid_laguerre_test04 ( dim_num, level_max );

  dim_num = 3;
  level_max = 6;
  sparse_grid_laguerre_test04 ( dim_num, level_max );

  dim_num = 10;
  level_max = 3;
  sparse_grid_laguerre_test04 ( dim_num, level_max );
%
%  Test monomial exactness for DIMENSION, LEVEL_MAX, DEGREE_MAX.
%
  dim_num = 2;
  level_max = 0;
  degree_max = 3;
  sparse_grid_laguerre_test05 ( dim_num, level_max, degree_max );

  dim_num = 2;
  level_max = 1;
  degree_max = 5;
  sparse_grid_laguerre_test05 ( dim_num, level_max, degree_max );

  dim_num = 2;
  level_max = 2;
  degree_max = 7;
  sparse_grid_laguerre_test05 ( dim_num, level_max, degree_max );

  dim_num = 2;
  level_max = 3;
  degree_max = 9;
  sparse_grid_laguerre_test05 ( dim_num, level_max, degree_max );

  dim_num = 2;
  level_max = 4;
  degree_max = 11;
  sparse_grid_laguerre_test05 ( dim_num, level_max, degree_max );

  dim_num = 2;
  level_max = 5;
  degree_max = 13;
  sparse_grid_laguerre_test05 ( dim_num, level_max, degree_max );

  dim_num = 3;
  level_max = 0;
  degree_max = 2;
  sparse_grid_laguerre_test05 ( dim_num, level_max, degree_max );

  dim_num = 3;
  level_max = 1;
  degree_max = 4;
  sparse_grid_laguerre_test05 ( dim_num, level_max, degree_max );

  dim_num = 3;
  level_max = 2;
  degree_max = 6;
  sparse_grid_laguerre_test05 ( dim_num, level_max, degree_max );

  dim_num = 3;
  level_max = 3;
  degree_max = 8;
  sparse_grid_laguerre_test05 ( dim_num, level_max, degree_max );
%
%  Show how to write a rule to a file.
%
  dim_num = 2;
  level_max = 3;

  sparse_grid_laguerre_test06 ( dim_num, level_max );
%
%  All done.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_LAGUERRE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
%%
% Test for OT color transfer.

addpath('../toolbox/');
addpath('../colors_functions/');
addpath('../image_blur/');
addpath('../blur_functions//');
addpath('../convolutional_wasserstein/');
addpath('../../data/images/colors/'); % low-res images
addpath('../../data/images/colors-big/'); % high-res images

% #bins
N = 160;
N = 140;
N = 200;

%%
% helpers

mmin = @(x)min(x(:));
mmax = @(x)max(x(:));
normalize = @(h)h/sum(h(:));
setfigname = @(name)set(gcf, 'Name', name, 'NumberTitle','off');
% plot histograms
delta = .1/N^2;
dispHist = @(x)-log(x+delta);
dispCell = @(H)cellfun(dispHist, H, 'UniformOutput', false);
Entropy = @(x)-sum(x(x>0).*log(x(x>0)));

%%
% Load images

if not(exist('names'))
    names = {'street-3', 'street-5'};
    names = {'nature-3' 'art'};
    names = {'street-4' 'street-6'};
    %
    names = {'street-2' 'street-7'};
    names = {'street-8' 'street-1'};
    %
    names = {'nature-7' 'street-4'};
    names = {'nature-7' 'street-2'};
    %
    names = {'yellow-2' 'purple-1'};
    %
    names = {'street-11' 'street-10'};
    %
    names = {'red-9' 'blue-7'};
end

p = length(names);

rep = ['../results/color-transport/' names{1} '-' names{2} '/'];
if not(exist(rep))
    mkdir(rep);
end

f = {}; fC = {};
for i=1:p
    f{i} = rescale( load_image(names{i}) );
    fC{i} = colorspace('RGB->LAB', f{i});
end

%%
% Compute range of ab components

arange = [];
brange = [];
for i=1:p
    arange(:,end+1) = range(fC{i}(:,:,2));
    brange(:,end+1) = range(fC{i}(:,:,3));
end
arange = [min(arange(1,:)); max(arange(2,:))];
brange = [min(brange(1,:)); max(brange(2,:))];

%%
% Blur kernel

metric_type = 'percep';
metric_type = 'unif';
switch metric_type
    case 'unif'
        mu = N/55;
        mu = N/45;
        mu = N/40;
        blur = load_filtering('imgaussian', N);
        K = @(x)blur(x,mu);
    case 'percep'
        % load the metric
        L = 50; % base lighting
        eta = .001;
        DiffDir = @(ab,h)deltaE2000([L ab],[L ab]+eta*[0 h])/eta;
        M = fit_metric_field(DiffDir, arange,brange,N);
        [e1,e2,l1,l2] = perform_tensor_decomp(M);
        M = perform_tensor_recomp(e1,e2,1./sqrt(l1),1./sqrt(l2));
        Id = zeros(N,N,2,2); Id(:,:,1,1) = 1; Id(:,:,2,2) = 1;
        M = M + 0*Id;
        % M = Id;
        % load the diffusion operator
        [blur, Delta,Grad] = blurAnisotropic(M);
        mu = .5; filtIter = 3;
        mu = .2; filtIter = 5;
        mu = .5/40; filtIter = 60;
        mu = .5/200; filtIter = 80;
        mu = .5/20; filtIter = 40;
        K = @(u)blur(u,mu,filtIter);
    otherwise
        error('Unknown metric type.');
end
Kv = @(x)apply_2d_func(K,x);

%%
% Compute histograms.

H = {}; fCi = {};
for i=1:p
    [H{i},ai,bi] = compute_histogram_2d(fC{i}(:,:,2),fC{i}(:,:,3),N, arange, brange);
    fCi{i} = cat(3,ai,bi);
end
imageplot( dispCell(H) );

%%
% Export colored histograms

for i=1:p
    H1 = render_lab_histogram(H{i},arange,brange);
    imwrite(rescale(H1), [rep 'density-' num2str(i) '.png'], 'png');
end

%%
% Compute transport coupling
%   pi = diag(w1)*K*diag(w0)
% and
%   pi*1 = w1.*K(w0) = p1
% and
%   pi'*1 = w1.*K(w0) = p0

options.tol = 1e-9;
options.tol = 0;
options.niter = 250;
options.verb = 2;
options.disp = @(w0,w1)imageplot( dispHist(w0.*K(w1)) );

clf;
[distances,w0,w1] = convolutionalDistance(H{1}, H{2}, [], K,[], options);

H1 = {}; f1 = {}; fCeq = {};
[f1{1},fCeq{1},H1{1}] = perform_color_transfer(K,w0, w1, fCi{1}, fC{1}(:,:,1),fC{2}(:,:,1), arange,brange);
[f1{2},fCeq{2},H1{2}] = perform_color_transfer(K,w1, w0, fCi{2}, fC{2}(:,:,1),fC{1}(:,:,1), arange,brange);

%%
% Save image

for i=1:p
    imwrite(rescale(f{i}), [rep names{i} '-original.jpg'], 'jpg');
%    imwrite(rescale(f1{i}), [rep names{i} '-' metric_type '-equalized.jpg'], 'jpg');
end

% display a comparison of histograms
figure(1); setfigname('Histograms');
clf;
imageplot(dispCell({H{1} H1{2}, H{2} H1{1}}), '', 2,2);
saveas(gcf, [rep names{1} '-' names{2} '-histograms.jpg'], 'jpg');

figure(2); setfigname('Images');
clf;
for i=1:2
    imageplot(f{i}, 'Original', 2,2,1+2*(i-1));
    imageplot(f1{i}, 'Equalized', 2,2,2+2*(i-1));
end
saveas(gcf, [rep names{1} '-' names{2} '-' metric_type '-images.jpg'], 'jpg');

%%
% Compute displacemet interpolation

options.disp = @(x)imageplot(dispHist(reshape(x, [N N])));

Q = 9;
tlist = linspace(0,1,Q);
% Q = 1; tlist = [1/2];

sharpening = 'none';
sharpening = 'min';
sharpening = 'arith';
sharpening = 'geom';

for i=1:Q
    t = tlist(i);
    w = [t 1-t];
    % select entropy bound
    switch sharpening
        case 'none'
            entropyLimit = [];
        case 'min'
            entropyLimit = min(Entropy(H{1}),Entropy(H{2}));
        case 'arith'
            entropyLimit = t*Entropy(H{1}) + (1-t)*Entropy(H{2});
        case 'geom'
            entropyLimit = exp( t*log(Entropy(H{1})) + (1-t)*log(Entropy(H{2})) );
        otherwise
            error('Unknown sharpening type');
    end
    % do the computation
    clf;
    Hv = reshape( cell2mat(H), [N*N p]); % as matrix
    [B,u] = convolutionalBarycenter( Hv, w, [], Kv, [],entropyLimit, options);
    B = reshape(B, [N N]);
    u = mat2cell( reshape(u, [N N p]), N, N, ones(p,1) );
    % save histogram
    H1 = render_lab_histogram(B,arange,brange);
    imwrite(rescale(H1), [rep 'barycenter-' metric_type '-' num2str(i) '.png'], 'png');
    % pi = diag(v{i})*K*diag(u{i})
    %   pi'*1 = u{i}.*K(v{i}) = bary
    %   pi*1 = v{i}.*K(u{i}) = H{i} ==> v{i} = H{i}./K(u{j})
    v = { H{1}./K(u{1}), H{2}./K(u{2}) };
    % compute average histogram
    opt.histinterp = 1-t;
    L = perform_histogram_equalization( fC{1}(:,:,1),fC{2}(:,:,1), opt );
    % perform both equalization
    H1 = {}; f1 = {}; fCeq = {};
    [f1{1},fCeq{1},H1{1}] = perform_color_transfer(K, u{1}, v{1}, fCi{1}, fC{1}(:,:,1),L, arange,brange);
    [f1{2},fCeq{2},H1{2}] = perform_color_transfer(K, u{2}, v{2}, fCi{2}, fC{2}(:,:,1),L, arange,brange);
    % write
    imwrite(rescale(f1{1}), [rep names{1} '-barycenter-' metric_type '-' num2str(i) '.jpg'], 'jpg');
    imwrite(rescale(f1{2}), [rep names{2} '-barycenter-' metric_type '-' num2str(i) '.jpg'], 'jpg');
end

author: gpeyre


text:
%MDL_QUADCOPTER Dynamic parameters for a quadrotor.
%
% MDL_QUADCOPTER is a script creates the workspace variable quad which
% describes the dynamic characterstics of a quadrotor flying robot.
%
% Properties::
%
% This is a structure with the following elements:
%
% nrotors   Number of rotors (1x1)
% J         Flyer rotational inertia matrix (3x3)
% h         Height of rotors above CoG (1x1)
% d         Length of flyer arms (1x1)
% nb        Number of blades per rotor (1x1)
% r         Rotor radius (1x1)
% c         Blade chord (1x1)
% e         Flapping hinge offset (1x1)
% Mb        Rotor blade mass (1x1)
% Mc        Estimated hub clamp mass (1x1)
% ec        Blade root clamp displacement (1x1)
% Ib        Rotor blade rotational inertia (1x1)
% Ic        Estimated root clamp inertia (1x1)
% mb        Static blade moment (1x1)
% Ir        Total rotor inertia (1x1)
% Ct        Non-dim. thrust coefficient (1x1)
% Cq        Non-dim. torque coefficient (1x1)
% sigma     Rotor solidity ratio (1x1)
% thetat    Blade tip angle (1x1)
% theta0    Blade root angle (1x1)
% theta1    Blade twist angle (1x1)
% theta75   3/4 blade angle (1x1)
% thetai    Blade ideal root approximation (1x1)
% a         Lift slope gradient (1x1)
% A         Rotor disc area (1x1)
% gamma     Lock number (1x1)
%
%
% Notes::
% - SI units are used.
%
% References::
% - Design, Construction and Control of a Large Quadrotor micro air vehicle.
%   P.Pounds, PhD thesis,
%   Australian National University, 2007.
%   http://www.eng.yale.edu/pep5/P_Pounds_Thesis_2008.pdf
% - This is a heavy lift quadrotor
%
% See also sl_quadrotor.

% MODEL: quadrotor

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

quad.nrotors = 4;                %   4 rotors
quad.g = 9.81;                   %   g       Gravity                             1x1
quad.rho = 1.184;                %   rho     Density of air                      1x1
quad.muv = 1.5e-5;               %   muv     Viscosity of air                    1x1

% Airframe
quad.M = 4;                      %   M       Mass                                1x1
Ixx = 0.082;
Iyy = 0.082;
Izz = 0.149;%0.160;
quad.J = diag([Ixx Iyy Izz]);    %   I       Flyer rotational inertia matrix     3x3

quad.h = -0.007;                 %   h       Height of rotors above CoG          1x1
quad.d = 0.315;                  %   d       Length of flyer arms                1x1

%Rotor
quad.nb = 2;                      %   b       Number of blades per rotor          1x1
quad.r = 0.165;                  %   r       Rotor radius                        1x1

quad.c = 0.018;                  %   c       Blade chord                         1x1

quad.e = 0.0;                    %   e       Flapping hinge offset               1x1
quad.Mb = 0.005;                 %   Mb      Rotor blade mass                    1x1
quad.Mc = 0.010;                 %   Mc      Estimated hub clamp mass            1x1
quad.ec = 0.004;                 %   ec      Blade root clamp displacement       1x1
quad.Ib = quad.Mb*(quad.r-quad.ec)^2/4 ;        %   Ib      Rotor blade rotational inertia      1x1
quad.Ic = quad.Mc*(quad.ec)^2/4;           %   Ic      Estimated root clamp inertia        1x1
quad.mb = quad.g*(quad.Mc*quad.ec/2+quad.Mb*quad.r/2);    %   mb      Static blade moment                 1x1
quad.Ir = quad.nb*(quad.Ib+quad.Ic);             %   Ir      Total rotor inertia                 1x1

quad.Ct = 0.0048;                %   Ct      Non-dim. thrust coefficient         1x1
quad.Cq = quad.Ct*sqrt(quad.Ct/2);         %   Cq      Non-dim. torque coefficient         1x1

quad.sigma = quad.c*quad.nb/(pi*quad.r);         %   sigma   Rotor solidity ratio                1x1
quad.thetat = 6.8*(pi/180);      %   thetat  Blade tip angle                     1x1
quad.theta0 = 14.6*(pi/180);     %   theta0  Blade root angle                    1x1
quad.theta1 = quad.thetat - quad.theta0;   %   theta1  Blade twist angle                   1x1
quad.theta75 = quad.theta0 + 0.75*quad.theta1;%   theta76 3/4 blade angle                     1x1
quad.thetai = quad.thetat*(quad.r/quad.e);      %   thetai  Blade ideal root approximation      1x1
quad.a = 5.5;                    %   a       Lift slope gradient                 1x1

% derived constants
quad.A = pi*quad.r^2;                 %   A       Rotor disc area                     1x1
quad.gamma = quad.rho*quad.a*quad.c*quad.r^4/(quad.Ib+quad.Ic);%   gamma   Lock number                         1x1

quad.b = quad.Ct*quad.rho*quad.A*quad.r^2; % T = b w^2
quad.k = quad.Cq*quad.rho*quad.A*quad.r^3; % Q = k w^2

quad.verbose = false;

author: Allopart

repo: rbpf-gmapping


text:
function jed = epoch_to_jed_eg_civil ( )

%*****************************************************************************80
%
%% EPOCH_TO_JED_EG_CIVIL: epoch of the Egyptian Civil calendar as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the epoch.
%
  jed = 1448637.5;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = p25_exact ( dim_num )

%*****************************************************************************80
%
%% P25_EXACT returns the exact integral for problem 25.
%
%  Discussion:
%
%    The formula in the reference seems to yield a result
%    that is too small by 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%
  c = 0.0;
  c = p25_r8 ( 'G', 'C', c );

  roundoff = eps;

  value = 1.0;

  term = 1.0;
  i = 0;

  while ( 1 )

    i = i + 1;
    term = term * c / i;

    term2 = term / ( i + 1 )^dim_num;

    if ( abs ( term2 ) <= roundoff * ( 1.0 + abs ( value ) ) )
      break
    end

    value = value + term2;

  end

  return
end


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

optionalParameters = {'constrOpt', '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))
solvers = checkGAMSSolvers('MIP');
if isempty(solvers)
    error('there is no GAMS solvers available to solver 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', 'OutputsFindMustU', @(x) ischar(x))
parser.addParamValue('outputFileName', 'MustUSet', @(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;
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
gamsMustUFunction = 'findMustU.gms';
%path of that function
pathGamsFunction = which(gamsMustUFunction);
if isempty(pathGamsFunction); error(['optForce: ' gamsMustUFunction ' 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, ['findMustUWithGAMS.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,'\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
% MustU Set
inputFolder = 'InputsMustU';
exportInputsMustToGAMS(model, 'U', minFluxesW, maxFluxesW, constrOpt,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 ' gamsMustUFunction ' lo=3 --myroot=' inputFolder '/ --solverName=' solverName ' gdx=GtoMU']);
else
    run = system(['gams ' gamsMustUFunction ' --myroot=' inputFolder '/ --solverName=' solverName ' gdx=GtoMU']);
end

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

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

%if findMustU.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 GtoMU; end;

    %if the problem was solved correctly, a variable named findMustU should be
    %inside of GtoMU. Otherwise, the wrong file is being read.
    try
        findMustU.name = 'findMustU';
        rgdx('GtoMU', findMustU);
        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 findMustU.gms
        %extract must U set found by findMustU.gms
        must.name = 'must';
        must.compress = 'true';
        must = rgdx('GtoMU', must);
        uelsMust = must.uels{1};

        %if the set is not empty
        if ~isempty(uelsMust)
            if printReport; fprintf(freport, '\na MustU set was found\n'); end;
            if verbose; fprintf('a MustU set was found\n'); end;
            mustUSet = uelsMust';
            %find position for reactions of the MustU set in model.rxns
            posMustU = cell2mat(arrayfun(@(x)find(strcmp(x, model.rxns)), mustUSet, 'UniformOutput', false))';

            %find values of fluxes achieved by each reaction in MustU set
            minFlux = zeros(size(mustUSet));
            maxFlux = zeros(size(mustUSet));

            %extract minimum values
            vmin.name = 'vmin';
            vmin.compress = 'true';
            vmin = rgdx('GtoMU', vmin);
            uels_vmin = vmin.uels{1};
            if ~isempty(uels_vmin)
                val_vmin = vmin.val(:,2);
                for i = 1:length(uels_vmin)
                    pos = strcmp(mustUSet, uels_vmin{i});
                    minFlux(pos == 1) = val_vmin(i);
                end
            end

            %extract miximum values
            vmax.name = 'vmax';
            vmax.compress = 'true';
            vmax = rgdx('GtoMU', vmax);
            uels_vmax = vmax.uels{1};
            if ~isempty(uels_vmax)
                val_vmax = vmax.val(:,2);
                for i = 1:length(uels_vmax)
                    pos = strcmp(mustUSet, uels_vmax{i});
                    maxFlux(pos == 1) = val_vmax(i);
                end
            end
        else
            if printReport; fprintf(freport, '\na MustU set was not found\n'); end;
            if verbose; fprintf('a MustU set was not found\n'); end;
        end

        % print info into an excel file if required by the user
        if printExcel
            if ~isempty(mustUSet)
                currentFolder = pwd;
                cd(outputFolder);
                Info = [{'Reactions'},{'Min Flux in Wild-type strain'},{'Max Flux in Wild-type strain'},{'Min Flux in Mutant strain'},{'Max Flux in Mutant strain'}];
                Info = [Info; [mustUSet, num2cell(minFluxesW(posMustU)), num2cell(maxFluxesW(posMustU)), num2cell(minFlux) ,num2cell(maxFlux)]];
                xlswrite('MustU_Info', Info);
                xlswrite(outputFileName, mustUSet);
                cd(currentFolder);
                if verbose; fprintf(['MustU set was printed in ' outputFileName '.xls  \n']); end;
                if printReport; fprintf(freport, ['\nMustU set was printed in ' outputFileName '.xls  \n']); end;
            else
                if verbose; fprintf('No mustU set was found. Therefore, no excel file was generated\n'); end;
                if printReport; fprintf(freport, '\nNo mustU 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(mustUSet)
                currentFolder = pwd;
                cd(outputFolder);
                f = fopen('MustU_Info.txt','w');
                fprintf(f, 'Reactions\tMin Flux in Wild-type strain\tMax Flux in Wild-type strain\tMin Flux in Mutant strain\tMax Flux in Mutant strain\n');
                for i = 1:length(posMustU)
                    fprintf(f, '%s\t%4.4f\t%4.4f\t%4.4f\t%4.4f\n', mustUSet{i}, minFluxesW(posMustU(i)), maxFluxesW(posMustU(i)), minFlux(i), maxFlux(i));
                end
                fclose(f);
                f = fopen([outputFileName '.txt'], 'w');
                for i = 1:length(posMustU)
                    fprintf(f, '%s\n', mustUSet{i});
                end
                fclose(f);
                cd(currentFolder);
                if verbose; fprintf(['MustU set was printed in ' outputFileName '.txt  \n']); end;
                if printReport; fprintf(freport, ['\nMustU set was printed in ' outputFileName '.txt  \n']); end;
            else
                if verbose; fprintf('No mustU set was found. Therefore, no excel file was generated\n'); end;
                if printReport; fprintf(freport, '\nNo mustU set was found. Therefore, no excel 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(gamsMustUFunction);

        %remove or move additional files that were generated during running
        if keepGamsOutputs
            if ~isdir(outputFolder); mkdir(outputFolder); end;
            movefile('GtoMU.gdx',outputFolder);
            movefile(regexprep(gamsMustUFunction, 'gms', 'lst'), outputFolder);
        else
            delete('GtoMU.gdx');
            delete(regexprep(gamsMustUFunction, '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
else
    %if GAMS was not executed correcttly
    if printReport && ~reportClosed; fprintf(freport, '\nGAMS was not executed correctly\n'); fclose(freport); end;
    if verbose; fprintf('GAMS was not executed correctly\n'); end;
    cd(workingPath);
    error('OptForce: GAMS was not executed correctly');
end
end

author: opencobra

repo: cobratoolbox


text:
function [ xtab, ytab, weight ] = hexagon_unit_set ( rule, order )

%*****************************************************************************80
%
%% HEXAGON_UNIT_SET sets a quadrature rule inside the unit hexagon in 2D.
%
%  Integration region:
%
%    The definition is given in terms of THETA, the angle in degrees of the
%    vector (X,Y).  The following six conditions apply, respectively,
%    between the bracketing values of THETA of 0, 60, 120, 180, 240,
%    300, and 360.
%
%                              0 <= Y <= - SQRT(3) * X + SQRT(3)
%                              0 <= Y <=                 SQRT(3)/2
%                              0 <= Y <=   SQRT(3) * X + SQRT(3)
%      - SQRT(3) * X - SQRT(3)   <= Y <= 0
%                    - SQRT(3)/2 <= Y <= 0
%        SQRT(3) * X - SQRT(3)   <= Y <= 0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Abramowitz and Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964.
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, integer RULE, the rule desired.
%      1, 1 point,  degree 1;
%      2, 4 points, degree 3;
%      3, 7 points, degree 3;
%      4, 7 points, degree 5;
%
%    Input, integer ORDER, the order of the desired rule.
%
%    Output, real XTAB(ORDER), YTAB(ORDER), the abscissase.
%
%    Output, real WEIGHT(ORDER), the weights.
%
  if ( rule == 1 )

    xtab(1) = 0.0;
    ytab(1) = 0.0;
    weight(1) = 1.0;
%
%  Stroud rule H2:3-1.
%
  elseif ( rule == 2 )

    a = sqrt ( 5.0 / 12.0 );
    b = 1.0 / 4.0;
    z = 0.0;

    xtab(1:4) =   [  a, -a,  z,  z ];
    ytab(1:4) =   [  z,  z,  a, -a ];
    weight(1:4) = [  b,  b,  b,  b ];
%
%  Stroud rule H2:3-2.
%
  elseif ( rule == 3 )

    a = sqrt ( 3.0 ) / 2.0;
    b =  0.5;
    c =  1.0;
    d =  5.0 / 72.0;
    e = 42.0 / 72.0;
    z =  0.0;

    xtab(1:7) =   [  z,  c, -c,  b, -b,  b, -b ];
    ytab(1:7) =   [  z,  z,  z,  a,  a, -a, -a ];
    weight(1:7) = [  e,  d,  d,  d,  d,  d,  d ];
%
%  Stroud rule H2:5-1.
%
  elseif ( rule == 4 )

    a = sqrt ( 14.0 ) / 5.0;
    b = sqrt ( 14.0 ) / 10.0;
    c = sqrt ( 42.0 ) / 10.0;
    d = 125.0 / 1008.0;
    e = 258.0 / 1008.0;
    z = 0.0;

    xtab(1:7) =   [ z,  a, -a,  b, -b,  b, -b ];
    ytab(1:7) =   [ z,  z,  z,  c,  c, -c, -c ];
    weight(1:7) = [ e,  d,  d,  d,  d,  d,  d ];

  else

    xtab = [];
    ytab = [];
    weight = [];

    fprintf ( 1, '\n' );
    fprintf ( 1, 'HEXAGON_UNIT_SET - Fatal error!\n' );
    fprintf ( 1, '  Nonexistent rule number %d requested.\n', rule );
    error ( 'HEXAGON_UNIT_SET - Fatal error!' );

  end

  return
end


text:
classdef MimMarkerFirst < MimGuiPlugin
    % MimMarkerFirst. Gui Plugin for switching transparency of zero
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimMarkerFirst is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'First marker'
        SelectedText = 'First marker'
        ToolTip = 'Go to slice containing first marker'
        Category = 'Slice navigation'
        Visibility = 'Dataset'
        Mode = 'Markers'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Icon = 'rewind.png'
        Location = 1
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.GotoFirstMarker;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && isequal(gui_app.GetCurrentModeName, MimModes.MarkerMode);
        end

        function is_selected = IsSelected(gui_app)
            is_selected = false;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function linpack_c_test29 ( )

%*****************************************************************************80
%
%% TEST29 tests CSIFA and CSISL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST29\n' );
  fprintf ( 1, '  For a single precision complex (C)\n' );
  fprintf ( 1, '  symmetric matrix (SI):\n' );
  fprintf ( 1, '  CSIFA factors the matrix.\n' );
  fprintf ( 1, '  CSISL solves a linear system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix order is N = %d\n', n );
%
%  Set the values of the matrix A.
%
  seed = 123456789;

  for i = 1 : n
    [ a(i,i), seed ] = c4_uniform_01 ( seed );
    for j = i+1 : n
      [ a(i,j), seed ] = c4_uniform_01 ( seed );
      a(j,i) = a(i,j);
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix A is\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  (%8f  %8f)', real ( a(i,j) ), imag ( a(i,j) ) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Set the values of the right hand side vector B.
%
  [ x, seed ] = c4vec_uniform_01 ( n, seed );

  b(1:n) = a(1:n,1:n) * transpose ( x(1:n) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The right hand side B is\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  (%8f  %8f)\n', real ( b(i) ), imag ( b(i) ) );
  end
%
%  Factor the matrix A.
%
  [ a, ipvt, info ] = csifa ( a, lda, n );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  CSIFA returned an error flag INFO = %d\n', info );
    return
  end
%
%  Solve the system.
%
  b = csisl ( a, lda, n, ipvt, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed                     Exact\n' );
  fprintf ( 1, '  Solution                     Solution\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  (%8f  %8f)  (%8f  %8f)\n', ...
      real ( b(i) ), imag ( b(i) ), real ( x(i) ), imag ( x(i) ) );
  end

  return
end


text:
classdef TestCopyMakeBorder
    %TestCopyMakeBorder

    methods (Static)
        function test_1
            src = imread(fullfile(mexopencv.root(),'test','img001.jpg'));
            [h,w,~] = size(src);
            border = [10, 20, 30, 40];
            types = {'Default', 'Reflect', 'Reflect101', 'Replicate', 'Wrap'};
            for i=1:numel(types)
                dst = cv.copyMakeBorder(src, border, 'BorderType',types{i});

                [hh,ww,~] = size(dst);
                assert(hh == (h+border(1)+border(2)));
                assert(ww == (w+border(3)+border(4)));
                assert(ndims(src) == ndims(dst));

                src2 = dst((1:h)+border(1), (1:w)+border(3), :);
                assert(isequal(src2, src));
            end
        end

        function test_2
            src = magic(5);
            dst = cv.copyMakeBorder(src, [1 2 3 4], ...
                'BorderType','Constant', 'Value',0);
        end

        function test_3
            src = magic(5);
            dst = cv.copyMakeBorder(src, 1, 2, 3, 4);
            dst = cv.copyMakeBorder(src, [1 2 3 4]);
        end

        function test_error_argnum
            try
                cv.copyMakeBorder();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
 function X = nufft(x, st)
%function X = nufft(x, st)
%|
%| Compute d-dimensional NUFFT of signal/image x
%|
%| in
%|	x	[N1 N2 ... Nd (L)]	L input image(s) of size
%|						N1 x N2 x ... x Nd
%|	st	structure		precomputed by nufft_init()
%| out
%|	X	[M (L)]			output spectra
%|
%| Copyright 2003-5-30, Jeff Fessler, University of Michigan

% if no arguments, then run some self tests
if (nargin < 1) || (nargin == 1 && streq(x, 'test'))
	help([mfilename '.m'])

	printm('Starting nufft() self test; after a wait, shows small errors')
	if 0
		Jd = [5 4 4];
		Nd = [23 15 19];
		alpha_user = {1, 1, 1}; % default alpha, beta
		beta_user = {0.5, 0.5, 0.5};
	else
		Jd = [5 6];
		Nd = [60 75];
		alpha_user = {1, 1};
		beta_user = {0.5, 0.5};
	end
	Kd = 2 * Nd;
	gam = 2*pi ./ Kd;
	n_shift = zeros(size(Nd));

	if 1
		oldarg = {Nd(1), Nd(2), Jd(1), Jd(2), Kd(1), Kd(2)};
		printm('err alf1 %g best %g', ...
			max(col(nufft2_err_mm('all', oldarg{:}, 1))), ...
			max(col(nufft2_err_mm('all', oldarg{:}, 'best'))) )
	end
	rng(0)
	x = randn(Nd);
%	x = [[1:Nd(1)]'*ones(1,3), ones(Nd(1), Nd(2)-3)]; % test signal
%	x = zeros(N1, N2); x(1,1) = 1;
	if 0 % test with uniform frequency locations
		o1 = 2 * pi * [0:(N1-1)]' / N1;
		o2 = 2 * pi * [0:(N2-1)]' / N2;
		[o1, o2] = ndgrid(o1, o2);
		Xf = fft2(x);
	else
		if length(Nd) == 3 % nonuniform frequencies
			[o1, o2, o3] = ndgrid(linspace(0,gam(1),11), ...
				linspace(0,gam(2),13), linspace(0,gam(3),5));
			om = [	[o1(:); [0 7.2 2.6 3.3]'], ...
				[o2(:); [0 4.2 -1. 5.5]'], ...
				[o3(:); [0 1.1 -2. 3.4]'] ];
		else
			[o1, o2] = ndgrid(linspace(-3*gam(1),3*gam(1),41), ...
				linspace(-2*gam(2),gam(2),31));
			om = [	[o1(:); [0 7.2 2.6 3.3]'], ...
				[o2(:); [0 4.2 -1. 5.5]'] ];
		end
		Y.d = dtft(x, om, 'n_shift', n_shift);
	end

	s.tab = nufft_init(om, Nd, Jd, Kd, n_shift, 'table', 2^12, 'minmax:kb');
	Y.tab = nufft(x, s.tab);
	printm('table0		max%%diff = %g', max_percent_diff(Y.d, Y.tab))

	s.mmkb = nufft_init(om, Nd, Jd, Kd, n_shift, 'minmax:kb');
	Y.mmkb = nufft(x, s.mmkb);
	printm('minmax:kb	max%%diff = %g', max_percent_diff(Y.d, Y.mmkb))

	if 0 % test multiple input case
		xx = x; xx(:,:,:,2) = x; xx(:,:,:,3) = x;
		Y.tmp = nufft(xx, s.mmkb);
		Y.tmp = Y.tmp(:,1);
		printm('multi	max%%diff = %g', max_percent_diff(Y.mmkb, Y.tmp))
	return
	end

	% kaiser with minmax best alpha,m
	s.kb = nufft_init(om, Nd, Jd, Kd, n_shift, 'kaiser');
	Y.kb = nufft(x, s.kb);
	printm('kaiser		max%%diff = %g', max_percent_diff(Y.d, Y.kb))

	% kaiser with user-specified supoptimal alpha,m for comparison
	s.ku = nufft_init(om, Nd, Jd, Kd, n_shift, 'kaiser', ...
		s.kb.kb_alf + 0.1*ones(size(s.kb.kb_alf)), s.kb.kb_m);
	Y.ku = nufft(x, s.ku);
	printm('kaiser-user	max%%diff = %g', max_percent_diff(Y.d, Y.ku))

	s.mmtu = nufft_init(om, Nd, Jd, Kd, n_shift, 'minmax:tuned');
	Y.mmtu = nufft(x, s.mmtu);
	printm('minmax:tuned	max%%diff = %g', max_percent_diff(Y.d, Y.mmtu))

	s.mm = nufft_init(om, Nd, Jd, Kd, n_shift, 'minmax:user', ...
		alpha_user, beta_user);
	Y.mm = nufft(x, s.mm);
	printm('minmax:user	max%%diff = %g', max_percent_diff(Y.d, Y.mm))

	if 0 % test 'uniform' scaling
		s.un = nufft_init(om, Nd, Jd, Kd, n_shift, 'uniform')
		Y.un = nufft(x, s.un);
		printm('user-unif max%%diff = %g', max_percent_diff(Y.mm, Y.un))
	return
	end

%	s1 = nufft_init(om, Nd, Jd, Kd, n_shift, 'loop', 'kaiser');
%	printm('kb loop max %% difference = %g', max_percent_diff(sn.p,s1.p))

 if 0 % test against old 2D version
	s2 = nufft2_init_kb(om, oldarg{:}, n_shift, 'kaiser');
%	Y2 = nufft2(x, s2);
%	printm('KB 2d vs nd max%%diff = %g', max_percent_diff(Y2, Y.kb))
	printm('KB 2d vs nd max%%diff = %g', max_percent_diff(s2.p,s.kb.p))

	s2 = nufft2_init(om, oldarg{:}, n_shift, 0, 'best');
	printm('tuned 2d vs nd max%%diff = %g', max_percent_diff(s2.p,s.mmtu.p))

	s2 = nufft2_init(om, oldarg{:}, n_shift, 0); % minmax
	printm('user 2d vs nd max%%diff = %g', max_percent_diff(s2.p,s.mm.p))

 end
return
end

Nd = st.Nd;
Kd = st.Kd;

dims = size(x);
dd = length(Nd);
if ndims(x) < dd, fail 'input signal has too few dimensions', end
if any(dims(1:dd) ~= Nd), fail 'input signal has wrong size', end

% the usual case is where L=1, i.e., there is just one input signal.
if ndims(x) == dd
	x = x .* st.sn; % apply scaling factors
	Xk = col(fftn_fast(x, Kd)); % [*Kd] oversampled FFT, padded at end

% otherwise, collapse all excess dimensions into just one
else
	xx = reshape(x, [prod(Nd) prod(dims((dd+1):end))]); % [*Nd *L]
	L = size(xx, 2);
	Xk = zeros(prod(Kd), L, class(xx)); % [*Kd *L]
	for ll=1:L
		xl = reshape(xx(:,ll), [Nd 1]); % l'th signal
		xl = xl .* st.sn; % scaling factors
		Xk(:,ll) = col(fftn_fast(xl, [Kd 1]));
	end
end

% interpolate using precomputed sparse matrix
% or with tabulated interpolator
if ~isfield(st, 'interp_table')
	X = st.p * Xk; % [M *L]
else
	X = feval(st.interp_table, st, Xk);
end

if ndims(x) > dd
	X = reshape(X, [st.M dims((dd+1):end)]); % [M (L)]
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [x, f, eflag, output, lambda] = qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt)
%QPS_GLPK  Linear Program Solver based on GLPK - GNU Linear Programming Kit.
%   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
%       QPS_GLPK(H, C, A, L, U, XMIN, XMAX, X0, OPT)
%   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_GLPK(PROBLEM)
%   A wrapper function providing a standardized interface for using
%   GLPK to solve the following LP (linear programming) problem:
%
%       min C'*X
%        X
%
%   subject to
%
%       L <= A*X <= U       (linear constraints)
%       XMIN <= X <= XMAX   (variable bounds)
%
%   Inputs (all optional except H, C, A and L):
%       H : IGNORED dummy matrix of quadratic cost coefficients
%           for QP problems, which GLPK does not handle
%       C : vector of linear cost coefficients
%       A, L, U : define the optional linear constraints. Default
%           values for the elements of L and U are -Inf and Inf,
%           respectively.
%       XMIN, XMAX : optional lower and upper bounds on the
%           X variables, defaults are -Inf and Inf, respectively.
%       X0 : optional starting value of optimization vector X (NOT USED)
%       OPT : optional options structure with the following fields,
%           all of which are also optional (default values shown in
%           parentheses)
%           verbose (0) - controls level of progress output displayed
%               0 = no progress output
%               1 = some progress output
%               2 = verbose progress output
%           glpk_opt - options struct for GLPK, value in verbose
%                   overrides these options
%       PROBLEM : The inputs can alternatively be supplied in a single
%           PROBLEM struct with fields corresponding to the input arguments
%           described above: H, c, A, l, u, xmin, xmax, x0, opt
%
%   Outputs:
%       X : solution vector
%       F : final objective function value
%       EXITFLAG : exit flag, 1 - optimal, <= 0 - infeasible, unbounded or other
%       OUTPUT : output struct with the following fields:
%           errnum - GLPK errnum output arg
%           status - GKPK status output arg
%           runtime - solver run time in seconds
%       LAMBDA : struct containing the Langrange and Kuhn-Tucker
%           multipliers on the constraints, with fields:
%           mu_l - lower (left-hand) limit on linear constraints
%           mu_u - upper (right-hand) limit on linear constraints
%           lower - lower bound on optimization variables
%           upper - upper bound on optimization variables
%
%   Note the calling syntax is almost identical to that of GLPK. The main
%   difference is that the linear constraints are specified with A, L, U
%   instead of A, B, Aeq, Beq.
%
%   Calling syntax options:
%       [x, f, exitflag, output, lambda] = ...
%           qps_glpk([], c, A, l, u, xmin, xmax, x0, opt)
%
%       x = qps_glpk([], c, A, l, u)
%       x = qps_glpk([], c, A, l, u, xmin, xmax)
%       x = qps_glpk([], c, A, l, u, xmin, xmax, x0)
%       x = qps_glpk([], c, A, l, u, xmin, xmax, x0, opt)
%       x = qps_glpk(problem), where problem is a struct with fields:
%                       H, c, A, l, u, xmin, xmax, x0, opt
%                       all fields except 'c', 'A' and 'l' or 'u' are optional
%       x = qps_glpk(...)
%       [x, f] = qps_glpk(...)
%       [x, f, exitflag] = qps_glpk(...)
%       [x, f, exitflag, output] = qps_glpk(...)
%       [x, f, exitflag, output, lambda] = qps_glpk(...)
%
%   Example: (based on example from 'doc linprog')
%       c = [-5; -4; -6];
%       A = [ 1  -1   1;
%            -3  -2  -4;
%             3   2   0];
%       l = [-Inf; -42; -Inf];
%       u = [20; Inf; 30];
%       xmin = [0; 0; 0];
%       opt = struct('verbose', 2);
%       [x, f, s, out, lambda] = qps_glpk([], c, A, l, u, xmin, [], [], opt);
%
%   See also QPS_MASTER, GLPK.

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

%% check for Optimization Toolbox
% if ~have_feature('quadprog')
%     error('qps_glpk: requires the MEX interface to GLPK');
% end

%%----- input argument handling  -----
%% gather inputs
if nargin == 1 && isstruct(H)       %% problem struct
    p = H;
    if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
    if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
    if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
    if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
    if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
    if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
    if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
    if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
    if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
else                                %% individual args
    if nargin < 9
        opt = [];
        if nargin < 8
            x0 = [];
            if nargin < 7
                xmax = [];
                if nargin < 6
                    xmin = [];
                end
            end
        end
    end
end

%% define nx, set default values for missing optional inputs
if isempty(H) || ~any(any(H))
    if isempty(A) && isempty(xmin) && isempty(xmax)
        error('qps_glpk: LP problem must include constraints or variable bounds');
    else
        if ~isempty(A)
            nx = size(A, 2);
        elseif ~isempty(xmin)
            nx = length(xmin);
        else    % if ~isempty(xmax)
            nx = length(xmax);
        end
    end
else
    error('qps_glpk: GLPK handles only LP problems, not QP problems');
    nx = size(H, 1);
end
if isempty(c)
    c = zeros(nx, 1);
end
if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
                                 (isempty(u) || all(u == Inf)))
    A = sparse(0,nx);           %% no limits => no linear constraints
end
nA = size(A, 1);                %% number of original linear constraints
if isempty(u)                   %% By default, linear inequalities are ...
    u = Inf(nA, 1);             %% ... unbounded above and ...
end
if isempty(l)
    l = -Inf(nA, 1);            %% ... unbounded below.
end
if isempty(xmin)                %% By default, optimization variables are ...
    xmin = -Inf(nx, 1);         %% ... unbounded below and ...
end
if isempty(xmax)
    xmax = Inf(nx, 1);          %% ... unbounded above.
end
if isempty(x0)
    x0 = zeros(nx, 1);
end

%% default options
if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
    verbose = opt.verbose;
else
    verbose = 0;
end

%% split up linear constraints
ieq = find( abs(u-l) <= eps );          %% equality
igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
AA = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
bb = [ u(ieq);    u(ilt);    -l(igt);    u(ibx);    -l(ibx)];

%% grab some dimensions
nlt = length(ilt);      %% number of upper bounded linear inequalities
ngt = length(igt);      %% number of lower bounded linear inequalities
nbx = length(ibx);      %% number of doubly bounded linear inequalities
neq = length(ieq);      %% number of equalities
nie = nlt+ngt+2*nbx;    %% number of inequalities

ctype = [repmat('S', neq, 1); repmat('U', nlt+ngt+2*nbx, 1)];
vtype = repmat('C', nx, 1);

%% set options struct for GLPK
if ~isempty(opt) && isfield(opt, 'glpk_opt') && ~isempty(opt.glpk_opt)
    glpk_opt = glpk_options(opt.glpk_opt);
else
    glpk_opt = glpk_options;
end
glpk_opt.msglev = verbose;

%% call the solver
t0 = tic;
[x, f, errnum, extra] = ...
    glpk(c, AA, bb, xmin, xmax, ctype, vtype, 1, glpk_opt);
output.runtime = toc(t0);

%% set exit flag
if isfield(extra, 'status')             %% status found in extra.status
    output.errnum = errnum;
    output.status = extra.status;
    eflag = -errnum;
    if eflag == 0 && extra.status == 5
        eflag = 1;
    end
else                                    %% status found in errnum
    output.errnum = [];
    output.status = errnum;
    if have_feature('octave')
        if errnum == 180 || errnum == 151 || errnum == 171
            eflag = 1;
        else
            eflag = 0;
        end
    else
        if errnum == 5
            eflag = 1;
        else
            eflag = 0;
        end
    end
end

%% repackage lambdas
if isempty(extra) || ~isfield(extra, 'lambda') || isempty(extra.lambda)
    lambda = struct( ...
        'mu_l', zeros(nA, 1), ...
        'mu_u', zeros(nA, 1), ...
        'lower', zeros(nx, 1), ...
        'upper', zeros(nx, 1) ...
    );
else
    lam.eqlin = extra.lambda(1:neq);
    lam.ineqlin = extra.lambda(neq+(1:nie));
    lam.lower = extra.redcosts;
    lam.upper = -extra.redcosts;
    lam.lower(lam.lower < 0) = 0;
    lam.upper(lam.upper < 0) = 0;
    kl = find(lam.eqlin > 0);   %% lower bound binding
    ku = find(lam.eqlin < 0);   %% upper bound binding

    mu_l = zeros(nA, 1);
    mu_l(ieq(kl)) = lam.eqlin(kl);
    mu_l(igt) = -lam.ineqlin(nlt+(1:ngt));
    mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx));

    mu_u = zeros(nA, 1);
    mu_u(ieq(ku)) = -lam.eqlin(ku);
    mu_u(ilt) = -lam.ineqlin(1:nlt);
    mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx));

    lambda = struct( ...
        'mu_l', mu_l, ...
        'mu_u', mu_u, ...
        'lower', lam.lower(1:nx), ...
        'upper', lam.upper(1:nx) ...
    );
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [x, infos] = admm_seq_conv_nmf(V, rank, t, in_options)
% ADMM based convolutive non-negative matrix factorization (ADMM-Conv-NMF).
%
% The problem of interest is defined as
%
%
% Given a non-negative matrix V, factorized non-negative matrices {W, H} are calculated.
%
%
% Inputs:
%       matrix      V
%       rank        rank
%
% Output:
%       w           solution of w
%       infos       information
%
% References:
%
%
% This file is part of NMFLibrary
%
% This file has been ported from
% convNMF_MM1.m and convNMF_MM2.m at https://github.com/lyn202206/ADMM-Convolutive-NMF
% by Yinan Li.
%
% Ported by H.Kasai on June 29, 2022
%
% Change log:
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.metric_type = 'beta-div';
    local_options.d_beta = 2;
    local_options.rho = 1;
    local_options.flag = 1;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    if ~isfield(options, 'x_init')
        W = zeros(m, rank, t);
        for i = 1 : t
            [init_factors, ~] = generate_init_factors(V, rank, init_options);
            W(:, :, i) = init_factors.W;
        end
        H = init_factors.H;
    else
        W = init_options.x_init.W;
        H = init_options.x_init.H;
    end

    % initialize
    epoch = 0;
    grad_calc_count = 0;

    options = check_divergence(options);
    sub_mode = sprintf('beta=%.1f', options.d_beta);
    if ~strcmp(options.metric_type, 'beta-div')
        sub_mode = options.metric_type;
    end
    method_name = sprintf('ADMM-Seq-Conv (%s)', sub_mode);

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % initialize for this algorithm
    X = zeros(m, n, t);
    for i=0:t-1
        tW = W(:, :, i+1);
        tH = shift_t(H, i);
        % initial X
        X(:, :, i+1) = tW * tH;
    end

    % some variable for test
    Xplus = X;
    U = zeros(m, n,t);

    alphaX = zeros(m, n,t);
    alphaH = zeros(rank,n);
    alphaW = zeros(m,rank,t);

    Wplus = W;
    Hplus = H;

    % store initial info
    clear infos;

    [Wcon, Hcon] = reconstruct_wh(W, H, t);
    [infos, f_val, optgap] = store_nmf_info(V, Wcon, Hcon, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('ADMM-Seq-Conv (%s): Epoch = 0000, cost = %.16e, optgap = %.4e\n', sub_mode, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % update H by accumulation
        P1 = zeros(rank,rank);
        P2 = zeros(rank,n);
        P3 = zeros(rank,n);
        for i = 0 : t-1
            tW = W(:, :, i+1);
            tX = shift_t(X(:, :, i+1), -i);
            talphaX = shift_t(alphaX(:, :, i+1), -i);
            P1 = P1 + tW'*tW;
            P2 = P2 + tW'*tX;
            P3 = P3 + tW'*talphaX;
        end
        H = (P1 + eye(rank)) \ (P2 + Hplus + 1/options.rho * (P3-alphaH));

        % update Hplus
        Hplus = max(H + 1/options.rho * alphaH, 0);

        % update alphaH
        alphaH = alphaH + options.rho * (H - Hplus);

        % update parameters sequentially
        for i = 0 : t-1
            % update Xplus in each time slice
            tW = Wplus(:, :, i+1);
            tH = shift_t(Hplus, i);
            % calculate Xplus
            Xplus(:, :, i+1) = tW * tH;

            % splite V in each time slice which result in U
            U(:, :, i+1) = (tW * tH + eps).*V./(sum(Xplus, 3) + t*eps);

            % update X in each time slice
            if options.d_beta ==2

                % update for Euclidean distance
                tW = W(:, :, i+1);
                tH = shift_t(H, i);
                tV = U(:, :, i+1);

                % update X
                X(:, :, i+1) = (options.rho*tW * tH + tV - alphaX(:, :, i+1)) / (1 + options.rho);

            elseif options.d_beta == 1

                % update for Kullback-Leibler divergence
                tW = W(:, :, i+1);
                tH = shift_t(H, i);
                tV = U(:, :, i+1);

                % update X
                b = options.rho * tW * tH - alphaX(:, :, i+1) - 1;
                X(:, :, i+1) = (b + sqrt(b.^2 + 4 * options.rho * tV)) / (2 * options.rho);

            elseif options.d_beta == 0

                % update for Itakura-Saito divergence
                tW = W(:, :, i+1);
                tH = shift_t(H, i);
                tV = U(:, :, i+1);

                % parameters for update
                A = alphaX(:, :, i+1) / options.rho - tW * tH;
                B = 1/(3 * options.rho) - A.^2/9;
                C = - A.^3/27 + A/(6 * options.rho) + tV / (2*options.rho);
                D = B.^3 + C.^2;

                % update X
                tX = X(:, :, i+1);
                tX(D>=0) = nthroot(C(D>=0)+sqrt(D(D>=0)),3) + ...
                nthroot(C(D>=0)-sqrt(D(D>=0)),3) - ...
                A(D>=0)/3;

                phi = acos(C(D<0) ./ ((-B(D<0)).^1.5));
                tX(D<0) = 2*sqrt(-B(D<0)).*cos(phi/3) - A(D<0)/3;
                X(:, :, i+1) = tX;
            else
                error('The options.d_beta you specified is not currently supported.')
            end

            % update alphaX
            alphaX(:, :, i+1) = alphaX(:, :, i+1) +  options.rho * (X(:, :, i+1) - tW * tH);

            % update W
            tH = shift_t(H, i);
            tX = X(:, :, i+1);
            P = tH * tH' + eye(rank);
            Q = tH * tX' + Wplus(:, :, i+1)' + 1/options.rho*(tH * alphaX(:, :, i+1)' - alphaW(:, :, i+1)');
            W(:, :, i+1) = (P \ Q)';

            % update Wplus
            Wplus(:, :, i+1) = max(W(:, :, i+1) + 1/options.rho * alphaW(:, :, i+1), 0);

            % udpate alphaW
            alphaW(:, :, i+1) = alphaW(:, :, i+1) + options.rho * (W(:, :, i+1) - Wplus(:, :, i+1));
        end

        %  if to update X before the update of U
        if options.flag == 1

            %  update Xplus in each time slice
            for i=0:t-1
                tW = Wplus(:, :, i+1);
                tH = shift_t(Hplus, i);
                % initial X
                Xplus(:, :, i+1) = tW * tH;
            end

            % splite V in each time slice which result in U
            for i=0:t-1
                tW = Wplus(:, :, i+1);
                tH = shift_t(Hplus, i);
                U(:, :, i+1) = (tW * tH + eps) .* V ./ (sum(Xplus, 3) + t*eps);
            end

            % update X in each time slice after the update of the whole W
            if options.d_beta == 2

                % update for Euclidean distance
                for i=0:t-1
                    tW = W(:, :, i+1);
                    tH = shift_t(H, i);
                    tV = U(:, :, i+1);
                    % update X
                    X(:, :, i+1) = (options.rho*tW * tH + tV-alphaX(:, :, i+1)) / (1 + options.rho);
                end

            elseif options.d_beta == 1

                % update for Kullback-Leibler divergence
                for i=0:t-1
                    tW = W(:, :, i+1);
                    tH = shift_t(H, i);
                    tV = U(:, :, i+1);
                    % update X
                    b = options.rho * tW * tH - alphaX(:, :, i+1) - 1;
                    X(:, :, i+1) = (b + sqrt(b.^2 + 4 * options.rho * tV)) / (2 * options.rho);
                end

            elseif options.d_beta == 0

                % update for Itakura-Saito divergence
                for i=0:t-1
                    tW = W(:, :, i+1);
                    tH = shift_t(H, i);
                    tV = U(:, :, i+1);

                    % parameters for update
                    A = alphaX(:, :, i+1) / options.rho - tW * tH;
                    B = 1 / (3 * options.rho) - A.^2/9;
                    C = - A.^3/27 + A / (6 * options.rho) + tV/(2 * options.rho);
                    D = B.^3 + C.^2;

                    % update X
                    tX = X(:, :, i+1);
                    tX(D>=0) = nthroot(C(D>=0)+sqrt(D(D>=0)),3) + ...
                    nthroot(C(D>=0)-sqrt(D(D>=0)),3) - ...
                    A(D>=0)/3;

                    phi = acos(C(D<0) ./ ((-B(D<0)).^1.5));
                    tX(D<0) = 2*sqrt(-B(D<0)).*cos(phi/3) - A(D<0)/3;
                    X(:, :, i+1) = tX;
                end
            else
                error('The options.d_beta you specified is not currently supported.')
            end

           % update for dual variables
            for i=0:t-1
                tW = W(:, :, i+1);
                tH = shift_t(H, i);
                % update alphaX
                alphaX(:, :, i+1) = alphaX(:, :, i+1) +  options.rho*(X(:, :, i+1)-tW * tH);
            end
        end

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        [Wcon, Hcon] = reconstruct_wh(W, H, t);
        infos = store_nmf_info(V, Wcon, Hcon, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    W = Wplus;
    H = Hplus;
    [W, H] = renormalize_convNMF(W, H);

    x.W = W;
    x.H = H;

end

function [W_concat, H_concat] = reconstruct_wh(W, H, t)

    W_concat = [];
    H_concat = [];
    for j = 1 : t
        W_concat = [W_concat W(:, :,j)];
        H_concat = [H_concat; shift_t(H, j-1)];
    end
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function scores= crossValidate(data,labels,obj,binVals,bins,params)

k = params.CVfolds;

boostIterations = 100;
% Learn classifier with all the data.

numEx = size(data,1);
posEx = labels == 1;
posCum = cumsum(posEx)/nnz(posEx);
negEx = ~posEx;
negCum = cumsum(negEx)/nnz(negEx);
numPos = sum(posEx);
numNeg = sum(negEx);

scores = zeros(1,numEx);
if numPos<1 || numNeg<1,
  obj.setStatus('Training set has only one kind of labels');
  return;
end

modLabels = sign( (labels==1)-0.5);

% rr = 1:numEx;
% splitPt = randi(numEx);
% rr = [rr(splitPt:end) rr(1:splitPt-1)];
% bStarts = round(linspace(1,numEx+1,k+1));

for bno = 1:k
%   curBlock = rr(bStarts(bno):bStarts(bno+1)-1);
%   curTest = ismember(1:numEx,curBlock);
%   curTrain = ~curTest;

  curPos = (posCum <= (bno/k)) & (posCum > ( (bno-1)/k)) & posEx;
  curNeg = (negCum <= (bno/k)) & (negCum > ( (bno-1)/k)) & negEx;

  curTest = curPos | curNeg ;
  curTrain = ~curTest;

  curTrainLabels = modLabels(curTrain);

  wt = getWeights(curTrainLabels);
  tt = tic;
  [~,curModel] = loglossboostLearnMod(data(curTrain,:),curTrainLabels,...
    boostIterations,wt,binVals,bins(:,curTrain),params);
  tScores = myBoostClassify(data(curTest,:),curModel);
  scores(curTest) = tScores;

  etime = toc(tt);
  obj.SetStatus('%d%% cross validation done. Time Remaining:%ds ',...
    round( bno/k*100), round((k-bno)*etime));
  drawnow();
end

author: kristinbranson

repo: JAABA


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Loading 3D Datasets for Manifold Learning  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Authors:
% Alberto Arrighi   (alberto.arrighi@epfl.ch)
% Ilaria Lauzana    (ilaria.lauzana@epfl.ch)
% Ludovico Novelli  (ludovico.novelli@epfl.ch)
%
% From Spring Semester Advanced Machine Learning Mini-Project
% Ecole Polytechnique Federale de Lausanne

%% GENERATE DIFFERENT DATASETS FOR MANIFOLD LEARNING EXAMPLES
%% 1) Swissroll

% Options
dataset_options                    = [];
dataset_options.numberOfPoints     = 2000;
dataset_options.name               = 'swissroll';
dataset_options.plot               = false;

% Generate dataset
[X, labels, updated_N, cmap] = ml_generate_manifold_dataset(dataset_options);

% Plot original data
plot_options            = [];
plot_options.is_eig     = false;
plot_options.points_size = 30;
plot_options.cmap       = cmap;
plot_options.title      = 'Swiss Roll';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1 = ml_plot_data(X,plot_options);
axis equal

%% 2) Swissroll with hole

% Options
dataset_options                    = [];
dataset_options.numberOfPoints     = 2000;
dataset_options.name               = 'swissroll_hole';
dataset_options.plot               = false;

% Generate dataset
[X, labels, updated_N, cmap] = ml_generate_manifold_dataset(dataset_options);

% Plot original data
plot_options            = [];
plot_options.is_eig     = false;
plot_options.points_size = 30;
plot_options.cmap       = cmap;
plot_options.title      = 'Swiss Roll w/Hole';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1 = ml_plot_data(X,plot_options);
axis equal

%% 3) S-curve

% Options
dataset_options                    = [];
dataset_options.numberOfPoints     = 2000;
dataset_options.name               = 'scurve';
dataset_options.plot               = false;

% Generate dataset
[X, labels, updated_N, cmap] = ml_generate_manifold_dataset(dataset_options);

% Plot original data
plot_options            = [];
plot_options.is_eig     = false;
plot_options.points_size = 30;
plot_options.cmap       = cmap;
plot_options.title      = 'S-Curve';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1 = ml_plot_data(X,plot_options);

%% 4) S-curve non-homogeneous density

% Options
dataset_options                    = [];
dataset_options.numberOfPoints     = 2000;
dataset_options.name               = 'scurve_mixed_density';
dataset_options.plot               = false;

% Generate dataset
[X, labels, updated_N, cmap] = ml_generate_manifold_dataset(dataset_options);

% Plot original data
plot_options             = [];
plot_options.is_eig      = false;
plot_options.points_size = 30;
plot_options.cmap        = cmap;
plot_options.title       = 'S-Curve Mixed';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1 = ml_plot_data(X,plot_options);

%% 5) Half sphere

% Options
dataset_options                    = [];
dataset_options.numberOfPoints     = 2000;
dataset_options.name               = 'sphere_hole';
dataset_options.plot               = false;

% Generate dataset
[X, labels, updated_N, cmap] = ml_generate_manifold_dataset(dataset_options);

% Plot original data
plot_options             = [];
plot_options.is_eig      = false;
plot_options.points_size = 30;
plot_options.cmap        = cmap;
plot_options.title       = 'Half-Sphere';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1 = ml_plot_data(X,plot_options);

%% 6) Broken Swissroll

% Options
options                 = [];
options.numberOfPoints  = 500;
options.name            = 'swissroll_broken';
options.plot            = false;

% Generate dataset
[X,labels,updated_N,cmap] = ml_generate_manifold_dataset(options);

% Plot original data
plot_options             = [];
plot_options.is_eig      = false;
plot_options.points_size = 30;
plot_options.labels      = labels;
plot_options.title       = 'Sphere w/Hole';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1 = ml_plot_data(X,plot_options);

author: epfl-lasa


text:
function [U,S,V] = spm_svd(X,U)
% Computationally efficient SVD (that can handle sparse arguments)
% FORMAT [U,S,V] = spm_svd(X,u)
% X    - (m x n) matrix
% u    - threshold (1 > u > 0) for normalized eigenvalues (default = 1e-6)
%      - a value of zero induces u = 64*eps
%
% U    - {m x p} singular vectors
% V    - {m x p} singular variates
% S    - {p x p} singular values
%__________________________________________________________________________
% Copyright (C) 1994-2011 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_svd.m 6110 2014-07-21 09:36:13Z karl $

% default thresholds - preclude singular vectors with small singular values
%--------------------------------------------------------------------------
if nargin < 2, U = 1e-6; end
if U >= 1; U = U - 1e-6; end
if U <= 0; U = 64*eps;   end

% deal with sparse matrices
%--------------------------------------------------------------------------
[M,N] = size(X);
p     = find(any(X,2));
q     = find(any(X,1));
X     = X(p,q);

% SVD
%--------------------------------------------------------------------------
[i, j, s] = find(X);
[m, n]    = size(X);
if any(i - j)

    % off-leading diagonal elements - full SVD
    %----------------------------------------------------------------------
    X     = full(X);
    if m > n

        [v, S, v] = svd(X'*X,0);
        S         = sparse(S);
        s         = diag(S);
        j         = find(s*length(s)/sum(s) > U);
        v         = v(:,j);
        u         = spm_en(X*v);
        S         = sqrt(S(j,j));

    elseif m < n

        [u, S, u] = svd(X*X',0);
        S         = sparse(S);
        s         = diag(S);
        j         = find(s*length(s)/sum(s) > U);
        u         = u(:,j);
        v         = spm_en(X'*u);
        S         = sqrt(S(j,j));

    else

        [u, S, v] = svd(X,0);
        S         = sparse(S);
        s         = diag(S).^2;
        j         = find(s*length(s)/sum(s) > U);
        v         = v(:,j);
        u         = u(:,j);
        S         = S(j,j);
    end

else
    S             = sparse(1:n,1:n,s,m,n);
    u             = speye(m,n);
    v             = speye(m,n);
    [i, j]        = sort(-s);
    S             = S(j,j);
    v             = v(:,j);
    u             = u(:,j);
    s             = diag(S).^2;
    j             = find(s*length(s)/sum(s) > U);
    v             = v(:,j);
    u             = u(:,j);
    S             = S(j,j);

end

% replace in full matrices
%--------------------------------------------------------------------------
j      = length(j);
U      = sparse(M,j);
V      = sparse(N,j);
if j
    U(p,:) = u;
    V(q,:) = v;
end

author: spm


text:
function view=loadResStdMap(view)
%
% view=loadMeanMap(view)
%
% Checks view.mapName field of the view structure.
% If view.mapName is not "resStdMap", then load it,
% and call setParameterMap.
%
% If you change this function make parallel changes in:
%   loadCorAnal, loadResStdMap, loadStdMap, loadMeanMap
%
% djh, 7/16/99, modified from loadCorAnal

if ~strcmp(view.mapName,'resStdMap')
   pathStr=fullfile(dataDir(view), 'resStdMap.mat');
   if ~exist(pathStr,'file')
     warning(['No ',pathStr,' file.  Run compute resStd map from Analysis menu.']);
     return
   else
      disp(['loading resStdMap from ',pathStr]);
      load(pathStr);
      view=setParameterMap(view,map,'resStdMap');
   end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft