text:
function pp = splinefit(varargin)
%SPLINEFIT Fit a spline to noisy data.
%   PP = SPLINEFIT(X,Y,BREAKS) fits a piecewise cubic spline with breaks
%   (knots) BREAKS to the noisy data (X,Y). X is a vector and Y is a vector
%   or an ND array. If Y is an ND array, then X(j) and Y(:,...,:,j) are
%   matched. Use PPVAL to evaluate PP.
%
%   PP = SPLINEFIT(X,Y,P) where P is a positive integer interpolates the
%   breaks linearly from the sorted locations of X. P is the number of
%   spline pieces and P+1 is the number of breaks.
%
%   OPTIONAL INPUT
%   Argument places 4 to 8 are reserved for optional input.
%   These optional arguments can be given in any order:
%
%   PP = SPLINEFIT(...,'p') applies periodic boundary conditions to
%   the spline. The period length is MAX(BREAKS)-MIN(BREAKS).
%
%   PP = SPLINEFIT(...,'r') uses robust fitting to reduce the influence
%   from outlying data points. Three iterations of weighted least squares
%   are performed. Weights are computed from previous residuals.
%
%   PP = SPLINEFIT(...,BETA), where 0 < BETA < 1, sets the robust fitting
%   parameter BETA and activates robust fitting ('r' can be omitted).
%   Default is BETA = 1/2. BETA close to 0 gives all data equal weighting.
%   Increase BETA to reduce the influence from outlying data. BETA close
%   to 1 may cause instability or rank deficiency.
%
%   PP = SPLINEFIT(...,N) sets the spline order to N. Default is a cubic
%   spline with order N = 4. A spline with P pieces has P+N-1 degrees of
%   freedom. With periodic boundary conditions the degrees of freedom are
%   reduced to P.
%
%   PP = SPLINEFIT(...,CON) applies linear constraints to the spline.
%   CON is a structure with fields 'xc', 'yc' and 'cc':
%       'xc', x-locations (vector)
%       'yc', y-values (vector or ND array)
%       'cc', coefficients (matrix).
%
%   Constraints are linear combinations of derivatives of order 0 to N-2
%   according to
%
%     cc(1,j)*y(x) + cc(2,j)*y'(x) + ... = yc(:,...,:,j),  x = xc(j).
%
%   The maximum number of rows for 'cc' is N-1. If omitted or empty 'cc'
%   defaults to a single row of ones. Default for 'yc' is a zero array.
%
%   EXAMPLES
%
%       % Noisy data
%       x = linspace(0,2*pi,100);
%       y = sin(x) + 0.1*randn(size(x));
%       % Breaks
%       breaks = [0:5,2*pi];
%
%       % Fit a spline of order 5
%       pp = splinefit(x,y,breaks,5);
%
%       % Fit a spline of order 3 with periodic boundary conditions
%       pp = splinefit(x,y,breaks,3,'p');
%
%       % Constraints: y(0) = 0, y'(0) = 1 and y(3) + y"(3) = 0
%       xc = [0 0 3];
%       yc = [0 1 0];
%       cc = [1 0 1; 0 1 0; 0 0 1];
%       con = struct('xc',xc,'yc',yc,'cc',cc);
%
%       % Fit a cubic spline with 8 pieces and constraints
%       pp = splinefit(x,y,8,con);
%
%       % Fit a spline of order 6 with constraints and periodicity
%       pp = splinefit(x,y,breaks,con,6,'p');
%
%   See also SPLINE, PPVAL, PPDIFF, PPINT

%   Author: Jonas Lundgren <splinefit@gmail.com> 2010

%   2009-05-06  Original SPLINEFIT.
%   2010-06-23  New version of SPLINEFIT based on B-splines.
%   2010-09-01  Robust fitting scheme added.
%   2010-09-01  Support for data containing NaNs.
%   2011-07-01  Robust fitting parameter added.

% Copyright (c) 2010, Jonas Lundgren
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% Check number of arguments
error(nargchk(3,7,nargin));

% Check arguments
[x,y,dim,breaks,n,periodic,beta,constr] = arguments(varargin{:});

% Evaluate B-splines
base = splinebase(breaks,n);
pieces = base.pieces;
A = ppval(base,x);

% Bin data
[junk,ibin] = histc(x,[-inf,breaks(2:end-1),inf]); %#ok

% Sparse system matrix
mx = numel(x);
ii = [ibin; ones(n-1,mx)];
ii = cumsum(ii,1);
jj = repmat(1:mx,n,1);
if periodic
    ii = mod(ii-1,pieces) + 1;
    A = sparse(ii,jj,A,pieces,mx);
else
    A = sparse(ii,jj,A,pieces+n-1,mx);
end

% Don't use the sparse solver for small problems
if pieces < 20*n/log(1.7*n)
    A = full(A);
end

% Solve
if isempty(constr)
    % Solve Min norm(u*A-y)
    u = lsqsolve(A,y,beta);
else
    % Evaluate constraints
    B = evalcon(base,constr,periodic);
    % Solve constraints
    [Z,u0] = solvecon(B,constr);
    % Solve Min norm(u*A-y), subject to u*B = yc
    y = y - u0*A;
    A = Z*A;
    v = lsqsolve(A,y,beta);
    u = u0 + v*Z;
end

% Periodic expansion of solution
if periodic
    jj = mod(0:pieces+n-2,pieces) + 1;
    u = u(:,jj);
end

% Compute polynomial coefficients
ii = [repmat(1:pieces,1,n); ones(n-1,n*pieces)];
ii = cumsum(ii,1);
jj = repmat(1:n*pieces,n,1);
C = sparse(ii,jj,base.coefs,pieces+n-1,n*pieces);
coefs = u*C;
coefs = reshape(coefs,[],n);

% Make piecewise polynomial
pp = mkpp(breaks,coefs,dim);

%--------------------------------------------------------------------------
function [x,y,dim,breaks,n,periodic,beta,constr] = arguments(varargin)
%ARGUMENTS Lengthy input checking
%   x           Noisy data x-locations (1 x mx)
%   y           Noisy data y-values (prod(dim) x mx)
%   dim         Leading dimensions of y
%   breaks      Breaks (1 x (pieces+1))
%   n           Spline order
%   periodic    True if periodic boundary conditions
%   beta        Robust fitting parameter, no robust fitting if beta = 0
%   constr      Constraint structure
%   constr.xc   x-locations (1 x nx)
%   constr.yc   y-values (prod(dim) x nx)
%   constr.cc   Coefficients (?? x nx)

% Reshape x-data
x = varargin{1};
mx = numel(x);
x = reshape(x,1,mx);

% Remove trailing singleton dimensions from y
y = varargin{2};
dim = size(y);
while numel(dim) > 1 && dim(end) == 1
    dim(end) = [];
end
my = dim(end);

% Leading dimensions of y
if numel(dim) > 1
    dim(end) = [];
else
    dim = 1;
end

% Reshape y-data
pdim = prod(dim);
y = reshape(y,pdim,my);

% Check data size
if mx ~= my
    mess = 'Last dimension of array y must equal length of vector x.';
    error('arguments:datasize',mess)
end

% Treat NaNs in x-data
inan = find(isnan(x));
if ~isempty(inan)
    x(inan) = [];
    y(:,inan) = [];
    mess = 'All data points with NaN as x-location will be ignored.';
    warning('arguments:nanx',mess)
end

% Treat NaNs in y-data
inan = find(any(isnan(y),1));
if ~isempty(inan)
    x(inan) = [];
    y(:,inan) = [];
    mess = 'All data points with NaN in their y-value will be ignored.';
    warning('arguments:nany',mess)
end

% Check number of data points
mx = numel(x);
if mx == 0
    error('arguments:nodata','There must be at least one data point.')
end

% Sort data
if any(diff(x) < 0)
    [x,isort] = sort(x);
    y = y(:,isort);
end

% Breaks
if isscalar(varargin{3})
    % Number of pieces
    p = varargin{3};
    if ~isreal(p) || ~isfinite(p) || p < 1 || fix(p) < p
        mess = 'Argument #3 must be a vector or a positive integer.';
        error('arguments:breaks1',mess)
    end
    if x(1) < x(end)
        % Interpolate breaks linearly from x-data
        dx = diff(x);
        ibreaks = linspace(1,mx,p+1);
        [junk,ibin] = histc(ibreaks,[0,2:mx-1,mx+1]); %#ok
        breaks = x(ibin) + dx(ibin).*(ibreaks-ibin);
    else
        breaks = x(1) + linspace(0,1,p+1);
    end
else
    % Vector of breaks
    breaks = reshape(varargin{3},1,[]);
    if isempty(breaks) || min(breaks) == max(breaks)
        mess = 'At least two unique breaks are required.';
        error('arguments:breaks2',mess);
    end
end

% Unique breaks
if any(diff(breaks) <= 0)
    breaks = unique(breaks);
end

% Optional input defaults
n = 4;                      % Cubic splines
periodic = false;           % No periodic boundaries
robust = false;             % No robust fitting scheme
beta = 0.5;                 % Robust fitting parameter
constr = [];                % No constraints

% Loop over optional arguments
for k = 4:nargin
    a = varargin{k};
    if ischar(a) && isscalar(a) && lower(a) == 'p'
        % Periodic conditions
        periodic = true;
    elseif ischar(a) && isscalar(a) && lower(a) == 'r'
        % Robust fitting scheme
        robust = true;
    elseif isreal(a) && isscalar(a) && isfinite(a) && a > 0 && a < 1
        % Robust fitting parameter
        beta = a;
        robust = true;
    elseif isreal(a) && isscalar(a) && isfinite(a) && a > 0 && fix(a) == a
        % Spline order
        n = a;
    elseif isstruct(a) && isscalar(a)
        % Constraint structure
        constr = a;
    else
        error('arguments:nonsense','Failed to interpret argument #%d.',k)
    end
end

% No robust fitting
if ~robust
    beta = 0;
end

% Check exterior data
h = diff(breaks);
xlim1 = breaks(1) - 0.01*h(1);
xlim2 = breaks(end) + 0.01*h(end);
if x(1) < xlim1 || x(end) > xlim2
    if periodic
        % Move data inside domain
        P = breaks(end) - breaks(1);
        x = mod(x-breaks(1),P) + breaks(1);
        % Sort
        [x,isort] = sort(x);
        y = y(:,isort);
    else
        mess = 'Some data points are outside the spline domain.';
        warning('arguments:exteriordata',mess)
    end
end

% Return
if isempty(constr)
    return
end

% Unpack constraints
xc = [];
yc = [];
cc = [];
names = fieldnames(constr);
for k = 1:numel(names)
    switch names{k}
        case {'xc'}
            xc = constr.xc;
        case {'yc'}
            yc = constr.yc;
        case {'cc'}
            cc = constr.cc;
        otherwise
            mess = 'Unknown field ''%s'' in constraint structure.';
            warning('arguments:unknownfield',mess,names{k})
    end
end

% Check xc
if isempty(xc)
    mess = 'Constraints contains no x-locations.';
    error('arguments:emptyxc',mess)
else
    nx = numel(xc);
    xc = reshape(xc,1,nx);
end

% Check yc
if isempty(yc)
    % Zero array
    yc = zeros(pdim,nx);
elseif numel(yc) == 1
    % Constant array
    yc = zeros(pdim,nx) + yc;
elseif numel(yc) ~= pdim*nx
    % Malformed array
    error('arguments:ycsize','Cannot reshape yc to size %dx%d.',pdim,nx)
else
    % Reshape array
    yc = reshape(yc,pdim,nx);
end

% Check cc
if isempty(cc)
    cc = ones(size(xc));
elseif numel(size(cc)) ~= 2
    error('arguments:ccsize1','Constraint coefficients cc must be 2D.')
elseif size(cc,2) ~= nx
    mess = 'Last dimension of cc must equal length of xc.';
    error('arguments:ccsize2',mess)
end

% Check high order derivatives
if size(cc,1) >= n
    if any(any(cc(n:end,:)))
        mess = 'Constraints involve derivatives of order %d or larger.';
        error('arguments:difforder',mess,n-1)
    end
    cc = cc(1:n-1,:);
end

% Check exterior constraints
if min(xc) < xlim1 || max(xc) > xlim2
    if periodic
        % Move constraints inside domain
        P = breaks(end) - breaks(1);
        xc = mod(xc-breaks(1),P) + breaks(1);
    else
        mess = 'Some constraints are outside the spline domain.';
        warning('arguments:exteriorconstr',mess)
    end
end

% Pack constraints
constr = struct('xc',xc,'yc',yc,'cc',cc);

%--------------------------------------------------------------------------
function pp = splinebase(breaks,n)
%SPLINEBASE Generate B-spline base PP of order N for breaks BREAKS

breaks = breaks(:);     % Breaks
breaks0 = breaks';      % Initial breaks
h = diff(breaks);       % Spacing
pieces = numel(h);      % Number of pieces
deg = n - 1;            % Polynomial degree

% Extend breaks periodically
if deg > 0
    if deg <= pieces
        hcopy = h;
    else
        hcopy = repmat(h,ceil(deg/pieces),1);
    end
    % to the left
    hl = hcopy(end:-1:end-deg+1);
    bl = breaks(1) - cumsum(hl);
    % and to the right
    hr = hcopy(1:deg);
    br = breaks(end) + cumsum(hr);
    % Add breaks
    breaks = [bl(deg:-1:1); breaks; br];
    h = diff(breaks);
    pieces = numel(h);
end

% Initiate polynomial coefficients
coefs = zeros(n*pieces,n);
coefs(1:n:end,1) = 1;

% Expand h
ii = [1:pieces; ones(deg,pieces)];
ii = cumsum(ii,1);
ii = min(ii,pieces);
H = h(ii(:));

% Recursive generation of B-splines
for k = 2:n
    % Antiderivatives of splines
    for j = 1:k-1
        coefs(:,j) = coefs(:,j).*H/(k-j);
    end
    Q = sum(coefs,2);
    Q = reshape(Q,n,pieces);
    Q = cumsum(Q,1);
    c0 = [zeros(1,pieces); Q(1:deg,:)];
    coefs(:,k) = c0(:);
    % Normalize antiderivatives by max value
    fmax = repmat(Q(n,:),n,1);
    fmax = fmax(:);
    for j = 1:k
        coefs(:,j) = coefs(:,j)./fmax;
    end
    % Diff of adjacent antiderivatives
    coefs(1:end-deg,1:k) = coefs(1:end-deg,1:k) - coefs(n:end,1:k);
    coefs(1:n:end,k) = 0;
end

% Scale coefficients
scale = ones(size(H));
for k = 1:n-1
    scale = scale./H;
    coefs(:,n-k) = scale.*coefs(:,n-k);
end

% Reduce number of pieces
pieces = pieces - 2*deg;

% Sort coefficients by interval number
ii = [n*(1:pieces); deg*ones(deg,pieces)];
ii = cumsum(ii,1);
coefs = coefs(ii(:),:);

% Make piecewise polynomial
pp = mkpp(breaks0,coefs,n);

%--------------------------------------------------------------------------
function B = evalcon(base,constr,periodic)
%EVALCON Evaluate linear constraints

% Unpack structures
breaks = base.breaks;
pieces = base.pieces;
n = base.order;
xc = constr.xc;
cc = constr.cc;

% Bin data
[junk,ibin] = histc(xc,[-inf,breaks(2:end-1),inf]); %#ok

% Evaluate constraints
nx = numel(xc);
B0 = zeros(n,nx);
for k = 1:size(cc,1)
    if any(cc(k,:))
        B0 = B0 + repmat(cc(k,:),n,1).*ppval(base,xc);
    end
    % Differentiate base
    coefs = base.coefs(:,1:n-k);
    for j = 1:n-k-1
        coefs(:,j) = (n-k-j+1)*coefs(:,j);
    end
    base.coefs = coefs;
    base.order = n-k;
end

% Sparse output
ii = [ibin; ones(n-1,nx)];
ii = cumsum(ii,1);
jj = repmat(1:nx,n,1);
if periodic
    ii = mod(ii-1,pieces) + 1;
    B = sparse(ii,jj,B0,pieces,nx);
else
    B = sparse(ii,jj,B0,pieces+n-1,nx);
end

%--------------------------------------------------------------------------
function [Z,u0] = solvecon(B,constr)
%SOLVECON Find a particular solution u0 and null space Z (Z*B = 0)
%         for constraint equation u*B = yc.

yc = constr.yc;
tol = 1000*eps;

% Remove blank rows
ii = any(B,2);
B2 = full(B(ii,:));

% Null space of B2
if isempty(B2)
    Z2 = [];
else
    % QR decomposition with column permutation
    [Q,R,dummy] = qr(B2); %#ok
    R = abs(R);
    jj = all(R < R(1)*tol, 2);
    Z2 = Q(:,jj)';
end

% Sizes
[m,ncon] = size(B);
m2 = size(B2,1);
nz = size(Z2,1);

% Sparse null space of B
Z = sparse(nz+1:nz+m-m2,find(~ii),1,nz+m-m2,m);
Z(1:nz,ii) = Z2;

% Warning rank deficient
if nz + ncon > m2
	mess = 'Rank deficient constraints, rank = %d.';
	warning('solvecon:deficient',mess,m2-nz);
end

% Particular solution
u0 = zeros(size(yc,1),m);
if any(yc(:))
    % Non-homogeneous case
	u0(:,ii) = yc/B2;
    % Check solution
	if norm(u0*B - yc,'fro') > norm(yc,'fro')*tol
        mess = 'Inconsistent constraints. No solution within tolerance.';
        error('solvecon:inconsistent',mess)
	end
end

%--------------------------------------------------------------------------
function u = lsqsolve(A,y,beta)
%LSQSOLVE Solve Min norm(u*A-y)

% Avoid sparse-complex limitations
if issparse(A) && ~isreal(y)
    A = full(A);
end

% Solution
u = y/A;

% Robust fitting
if beta > 0
    [m,n] = size(y);
    alpha = 0.5*beta/(1-beta)/m;
    for k = 1:3
        % Residual
        r = u*A - y;
        rr = r.*conj(r);
        rrmean = sum(rr,2)/n;
        rrmean(~rrmean) = 1;
        rrhat = (alpha./rrmean)'*rr;
        % Weights
        w = exp(-rrhat);
        spw = spdiags(w',0,n,n);
        % Solve weighted problem
        u = (y*spw)/(A*spw);
    end
end

repo: Aero-matlab


text:
%% import_ctgr
% Below is a demonstration of the features of the |import_ctgr| function

%%
clear; close all; clc;

%% Syntax
% |[outputStruct]=import_ctgr(fileName);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function c=recalculate_points(b)
% function c=recalculate_points(b)

if ~strcmpi(class(b),'bezier')
  error('Expecting a bezier object');
end

c=b;

c.x1=c.x0+c.cx/3;
c.x2=c.x1+(c.cx+c.bx)/3;
c.x3=c.x0+c.cx+c.bx+c.ax;

c.y1=c.y0+c.cy/3;
c.y2=c.y1+(c.cy+c.by)/3;
c.y3=c.y0+c.cy+c.by+c.ay;

author: Sable

repo: mcbench-benchmarks


text:
function fval = f1(x)
% Unimodal function f_1

Bound=[-100 100];

if nargin == 0
    fval = Bound;
else
    fval = sum(x.^2);
end

author: sfvsfv

repo: Mathematical-modeling


text:
function [XU] = Associate(Population,W,N)
% Associate each solution with a neuron

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Chao He

     A  = 1 : N;
     U  = 1 : N;
     XU = zeros(1,N);
     for i = 1 : N
         x = randi(length(A));
         [~,u] = min(pdist2(Population(A(x)).dec,W(U,:)));
         XU(U(u)) = A(x);
         A(x)     = [];
         U(u)     = [];
     end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function point_merge_test04 ( m, n, n_unique, tol, seed )

%*****************************************************************************80
%
%% POINT_MERGE_TEST04 tests uniqueness indexing with a tolerance.
%
%  Discussion:
%
%    POINT_RADIAL_TOL_UNIQUE_COUNT uses an algorithm that should be,
%      in general, O(N);
%    POINT_TOL_UNIQUE_COUNT uses an O(N^2) algorithm.
%
%    For this test, we just want to make sure the algorithms agree
%    in the counting.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 July 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POINT_MERGE_TEST04\n' );
  fprintf ( 1, '  To index the unique columns in an R8COL, we call\n' );
  fprintf ( 1, '  POINT_RADIAL_TOL_UNIQUE_COUNT, (with random center)\n' );
  fprintf ( 1, '  POINT_TOL_UNIQUE_COUNT, (with zero tolerance)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  M = %d\n', m );
  fprintf ( 1, '  N = %d\n', n );
  fprintf ( 1, '  TOL = %f\n', tol );
  fprintf ( 1, '  SEED = %d\n', seed );

  [ a, seed ] = r8col_duplicates ( m, n, n_unique, seed );

  r8mat_transpose_print ( m, n, a, '  Matrix with N_UNIQUE unique columns:' );
%
%  The form of the tolerance test means that if two vectors are initially
%  equal, they remain "tolerably equal" after the addition of random
%  perturbation vectors whose 2-norm is no greater than TOL/2.
%
  for j = 1 : n
    [ r, seed ] = r8vec_uniform_01 ( m, seed );
    r(1:m) = r(1:m) / sqrt ( sum ( r(1:m).^2 ) );
    a(1:m,j) = a(1:m,j) + 0.5 * tol * r(1:m,1);
  end

  r8mat_transpose_print ( m, n, a, '  Blurred matrix:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  N_UNIQUE =                      %d\n', n_unique );

  [ unique_num, undx, xdnu, seed ] = point_radial_tol_unique_index ( m, n, a, ...
    tol, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  POINT_RADIAL_TOL_UNIQUE_INDEX\n' );
  fprintf ( 1, '  Unique_num = %d\n', unique_num );

  i4vec_print ( unique_num, undx, '  UNDX:' );

  i4vec_print ( n, xdnu, '  XDNU:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  List of nonunique points P(J), represented by\n' );
  fprintf ( 1, '  point with index I(J).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  J, P(J)\n' );
  fprintf ( 1, '  I(J), P(I(J))\n' );
  fprintf ( 1, '  || P(J) - P(I(J)) || (should be <= TOL)\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    k = undx(xdnu(j));
    if ( j ~= k )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  %4d', j );
      for i = 1 : m
        fprintf ( 1, '  %14f', a(i,j) );
      end
      fprintf ( 1, '\n' );
      fprintf ( 1, '  %4d', k );
      for i = 1 : m
        fprintf ( 1, '  %14f', a(i,k) );
      end
      fprintf ( 1, '\n' );
      dist = sqrt ( sum ( ( a(1:m,j) - a(1:m,k) ).^2 ) );
      fprintf ( 1, '          %10f\n', dist );
    end
  end
%
%  The interpretation of XDNU is simpler for POINT_TOL_UNIQUE_INDEX.
%
  [ unique_num, xdnu ] = point_tol_unique_index ( m, n, a, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  POINT_TOL_UNIQUE_INDEX\n' );
  fprintf ( 1, '  Unique_num = %d\n', unique_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  List of nonunique points P(J), represented by\n' );
  fprintf ( 1, '  point with index I(J).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  J, P(J)\n' );
  fprintf ( 1, '  I(J), P(I(J))\n' );
  fprintf ( 1, '  || P(J) - P(I(J)) || (should be <= TOL)\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    k = xdnu(j);
    if ( j ~= k )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  %4d', j );
      for i = 1 : m
        fprintf ( 1, '  %14f', a(i,j) );
      end
      fprintf ( 1, '\n' );
      fprintf ( 1, '  %4d', k );
      for i = 1 : m
        fprintf ( 1, '  %14f', a(i,k) );
      end
      fprintf ( 1, '\n' );
      dist = sqrt ( sum ( ( a(1:m,j) - a(1:m,k) ).^2 ) );
      fprintf ( 1, '          %10f\n', dist );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function s = s_adjustl ( s )

%*****************************************************************************80
%
%% S_ADJUSTL flushes a string left.
%
%  Discussion:
%
%    Both blanks and tabs are treated as "white space".
%
%    This routine is similar to the FORTRAN90 ADJUSTL routine.
%
%  Example:
%
%    Input             Output
%
%    '     Hello'      'Hello     '
%    ' Hi there!  '    'Hi there!   '
%    'Fred  '          'Fred  '
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be shifted left.
%
%    Output, string S, has been shifted left.
%
  TAB = char ( 9 );
%
%  Check the length of the string to the last nonblank.
%  If nonpositive, return.
%
  len = s_len_trim ( s );

  if ( len <= 0 )
    return
  end
%
%  Find NONB, the location of the first nonblank, nontab.
%
  nonb = 0;

  for i = 1 : len

    if ( s(i) ~= ' ' && s(i) ~= TAB )
      nonb = i;
      break
    end

  end

  if ( nonb == 0 )
    s = ' ';
    return
  end
%
%  Shift the string left.
%
  if ( 1 < nonb )
    for i = 1 : len + 1 - nonb
      s(i) = s(i+nonb-1);
    end
  end
%
%  Blank out the end of the string.
%
  s(len+2-nonb:len) = ' ';

  return
end


text:
function recall = Recall(SEG, GT)
    % SEG, GT are the binary segmentation and ground truth areas, respectively.
    % recall
    recall = double(sum(uint8(SEG(:) & GT(:)))) / double(sum(uint8(GT(:))));
end

repo: cluster


text:
function c8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% C8MAT_PRINT_SOME prints some of a C8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the matrix.
%
%    Input, complex A(M,N), the matrix.
%
%    Input, integer ILO, JLO, IHI, JHI, the first row and
%    column, and the last row and column to be printed.
%
%    Input, string TITLE, a title.
%
  incx = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
%
%  Print the columns of the matrix, in strips of INCX.
%
  for j2lo = jlo : incx : jhi

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    inc = j2hi + 1 - j2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );
    for j = j2lo : j2hi
      j2 = j + 1 - j2lo;
      fprintf ( 1, '%10d          ', j );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '  ---\n' );
%
%  Determine the range of the rows in this strip.
%
    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi
%
%  Print out (up to) INCX entries in row I, that lie in the current strip.
%
      for j2 = 1 : inc

        j = j2lo - 1 + j2;

        if ( a(i,j) == 0.0 )
          fprintf ( 1, '       0.0          ' );
        elseif ( imag ( a(i,j) ) == 0.0 )
          fprintf ( 1, ' %9g          ', real ( a(i,j) ) );
        else
          fprintf ( 1, ' %9g %9g', real ( a(i,j) ), imag ( a(i,j) ) );
       end

      end

      fprintf ( 1, '\n' );

    end

  end

  fprintf ( 1, '\n' );

  return
end


text:
function chunk = onl_peek(streamname,samples_to_get,unit,channels_to_get)
% Peek into an online stream (generates an EEG-set like view into it).
% Chunk = onl_peek(StreamName,DesiredLength,LengthUnit,DesiredChannels)
%
% This function returns an EEGLAB dataset struct that hold the last k seconds of an online stream.
% An online stream is a data structure in the workspace that can be created with onl_newstream;
% data can be appended to it by onl_append or onl_read_background.
%
% In:
%   StreamName : Name of the online stream data structure in the MATLAB workspace to read from.
%                Must have previously been created with onl_newstream.
%
%   DesiredLength : length of the view that should be generated; should not be longer than the
%                   buffer capacity (default: 10)
%
%   LengthUnit : Can be one of the following options:
%                 * 'seconds': return the last PeekLength seconds (default)
%                 * 'samples': return the last PeekLength samples
%                 * 'index': return all samples newer than PeekLength (as a sample index)
%
%   DesiredChannels : range of channels to return ([] = all channels); (default: [])
%
% Out:
%   Chunk : An EEGLAB data set that contains the most recent data of a given stream
%
% Example:
%   % get the last 5 seconds of the stream
%   EEG = onl_peek('mystream',5)
%
%   % get the last 128 samples of the stream
%   EEG = onl_peek('mystream',128,'samples')
%
%   % get all samples past the 5000'th sample in the stream
%   EEG = onl_peek('mystream',5000,'index')
%
%   % get the default amount of data, for channels 1:10
%   EEG = onl_peek('mystream',[],[],1:10)
%
% See also:
%   onl_newstream, onl_append, onl_filtered
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-03

% try to get the stream from the base workspace
try
    chunk = evalin('base',streamname);
catch e
    % diagnose the error
    if nargin < 1
        error('You need to pass at least the name of a previously created stream.'); end
    if ~ischar(streamname) || isempty(streamname) || ~isvarname(streamname)
        error('The given StreamName argument must be the name of a variable name in the MATLAB workspace, but was: %s',hlp_tostring(streamname,10000)); end
    if strcmp(e.identifier,'MATLAB:badsubscript')
        error('BCILAB:onl_peek:improper_resolve','The raw data required by the predictor does not list the name of the needed source stream; this is likely a problem in onl_newstream/onl_newpredictor.');
    else
        error('BCILAB:onl_peek:stream_not_found','The stream named %s was not found in the base workspace.',hlp_tostring(streamname));
    end
end

try
    % set further default arguments
    if nargin < 4 || isempty(channels_to_get)
        channels_to_get = 1:size(chunk.buffer,1);
        if nargin < 3 || isempty(unit)
            unit = 'seconds'; end
        if nargin < 2 || isempty(samples_to_get)
            samples_to_get = 10; end
    end

    % determine the amount of data needed
    switch unit
        case 'index'
            samples_to_get = min(chunk.buffer_len,chunk.smax-samples_to_get);
        case 'samples'
            samples_to_get = min(chunk.buffer_len,samples_to_get);
        case 'seconds'
            samples_to_get = min(chunk.buffer_len,round(chunk.srate*samples_to_get));
        otherwise
            error('Unrecognized LengthUnit specified: %s',hlp_tostring(unit));
    end

    % extract the desired interval from the .buffer field and move it to .data
    range = 1+mod(chunk.smax-samples_to_get:chunk.smax-1,chunk.buffer_len);
    chunk.data = chunk.buffer(channels_to_get,range);

    % extract the markers, if there are any
    if chunk.mmax
        [ranks,sample_indices,record_indices] = find(chunk.marker_pos(:,range));
        if any(ranks)
            chunk.event = chunk.marker_buffer(1+mod(record_indices-1,chunk.marker_buffer_len));
            [chunk.event.latency] = arraydeal(sample_indices(:) + [chunk.event.latency]');
        end
    end

    % update meta-data
    [chunk.nbchan,chunk.pnts,chunk.trials] = size(chunk.data);
    chunk.chanlocs = chunk.chanlocs(channels_to_get);
    if chunk.xmin == 0
        chunk.xmax = (chunk.smax-1)/chunk.srate;
        chunk.xmin = chunk.xmax - (chunk.pnts-1)/chunk.srate;
    else
        chunk.xmax = chunk.xmin + (chunk.pnts-1)/chunk.srate;
    end

catch e
    if ~ischar(unit) || isempty(unit)
        error('The given Unit argument must be a string, but was: %s',hlp_tostring(unit,10000)); end
    if ~isnumeric(samples_to_get) || ~isscalar(samples_to_get)
        error('The given DesiredLength argument must be a numeric scalar, but was: %s',hlp_tostring(samples_to_get,10000)); end
    if ~isnumeric(channels_to_get) || min(size(channels_to_get) ~= 1)
        error('The given DesiredChannels argument must be a numeric vector, but was: %s',hlp_tostring(channels_to_get,10000)); end
    error('Failed to read from stream %s with error: %s\nPossibly the stream variable has been overwritten or you have the wrong name.',streamname,e.message);
end

repo: REST


text:
function sF = dthetadtheta(sF)
% second derivative in direction theta

s = size(sF);
sF = reshape(sF, []);

sF.fhat(1) = 0; % exclude some special cases
fhat_theta = zeros((sF.bandwidth+2)^2, length(sF));
for m = 0:sF.bandwidth+1
  if 0 <= m-1
    fhat_theta(m^2+1:(m+1)^2, :) = (m-1)*sqrt((m^2-(-m:m)'.^2)/((2*m-1)*(2*m+1))).*[zeros(1, length(sF)); sF.fhat((m-1)^2+1:m^2, :); zeros(1, length(sF))];
  end
  if m+1 <= sF.bandwidth
    fhat_theta(m^2+1:(m+1)^2, :) = fhat_theta(m^2+1:(m+1)^2, :)-(m+2)*sqrt(((m+1)^2-(-m:m)'.^2)/((2*m+1)*(2*m+3))).*sF.fhat((m+1)^2+2:(m+2)^2-1, :);
  end
end

sF_theta = S2FunHarmonic(fhat_theta);

fhat_theta_theta = zeros((sF_theta.bandwidth+2)^2, length(sF));
for m = 0:sF_theta.bandwidth+1
  if 0 <= m-1
    fhat_theta_theta(m^2+1:(m+1)^2, :) = (m-1)*sqrt((m^2-(-m:m)'.^2)/((2*m-1)*(2*m+1))).*[zeros(1, length(sF)); sF_theta.fhat((m-1)^2+1:m^2, :); zeros(1, length(sF))];
  end
  if m+1 <= sF_theta.bandwidth
    fhat_theta_theta(m^2+1:(m+1)^2, :) = fhat_theta_theta(m^2+1:(m+1)^2, :)-(m+2)*sqrt(((m+1)^2-(-m:m)'.^2)/((2*m+1)*(2*m+3))).*sF_theta.fhat((m+1)^2+2:(m+2)^2-1, :);
  end
end

sF_theta_theta = S2FunHarmonic(fhat_theta_theta);

f = @(v) (sF_theta_theta.eval(v)-cos(v.theta).*sF_theta.eval(v))./max(sin(v.theta).^2, 0.1);
sF = S2FunHarmonic.quadrature(f, 'bandwidth', sF.bandwidth);
sF = reshape(sF, s);

end

author: mtex-toolbox

repo: mtex

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


text:
% MAIN  --  Polar Point Mass
%
% This script runs a simulation of the polar point mass, with a controller
% that is working to hold the mass at a fixed distance from the origin
%

P.m = 1;
P.g = 9.81;
P.l = 1;
P.freq = 3*sqrt(P.g/P.l);   %Response frequence in controller
P.damp = 1.0;   %Damping ratio in controller
P.eNom = 2.5*P.m*P.g*P.l;

% z = [r; th; dr; dth];

z0min = [0.5*P.l; -pi; -0.2; -0.2];
z0max = [1.5*P.l; pi;  0.2;  0.2];

z0 = z0min + (z0max-z0min).*rand(4,1);

userFunc = @(t,z)dynamics(t,z,controller(z,P),P);

tSpan = [0;3];

sol = ode45(userFunc,tSpan,z0);

t = linspace(tSpan(1),tSpan(2),500);
z = deval(sol,t);
u = controller(z,P);

figure(234); clf;

subplot(3,2,1);
plot(t,z(1,:))
ylabel('r')

subplot(3,2,3);
plot(t,z(3,:))
ylabel('dr')

subplot(3,2,5);
plot(t,u)
ylabel('u')

subplot(3,2,2);
plot(t,z(2,:))
ylabel('th')

subplot(3,2,4);
plot(t,z(4,:))
ylabel('dth')

author: MatthewPeterKelly

repo: dscTutorials


text:
function dvh = matRad_calcDVH(cst,doseCube,dvhType,doseGrid)
% matRad dvh calculation
%
% call
%   dvh = matRad_calcDVH(cst,doseCube)
%   dvh = matRad_calcDVH(cst,doseCube,dvhType)
%   dvh = matRad_calcDVH(cst,doseCube,doseGrid)
%   dvh = matRad_calcDVH(cst,doseCube,dvhType,doseGrid)
%
% input
%   cst:        matRad cst struct
%   doseCube:   arbitrary doseCube (e.g. physicalDose)
%   dvhType:    (optional) string, 'cum' for cumulative, 'diff' for differential
%               dvh
%   doseGrid:   (optional) use predefined evaluation points. Useful when
%               comparing multiple realizations
%
% output
%   dose volume histogram
%
% References
%   van't Riet et. al., IJROBP, 1997 Feb 1;37(3):731-6.
%   Kataria et. al., J Med Phys. 2012 Oct-Dec; 37(4): 207�213.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2016 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('dvhType','var') || isempty(dvhType)
    dvhType = 'cum';
end

if ~exist('doseGrid', 'var') || isempty(doseGrid)
    maxDose = max(doseCube(:));
    minDose = min(doseCube(:));

    % get dvhPoints for every structure and every scenario the same
    n = 1000;
    if strcmp(dvhType, 'cum')
        doseGrid = linspace(0,maxDose*1.05,n);
    elseif strcmp(dvhType, 'diff')
        doseGrid = linspace(0.95*minDose,maxDose*1.05,n);
    end
end

numOfVois = size(cst,1);
dvh = struct;
for i = 1:numOfVois
    dvh(i).doseGrid     = doseGrid;
    dvh(i).volumePoints = getDVHPoints(cst, i, doseCube, doseGrid, dvhType);
    dvh(i).name         = cst{i,2};
end

end %eof

function dvh = getDVHPoints(cst, sIx, doseCube, dvhPoints, dvhType)
n = numel(dvhPoints);
dvh         = NaN * ones(1,n);
indices     = cst{sIx,4}{1};
numOfVoxels = numel(indices);

doseInVoi   = doseCube(indices);

switch dvhType
    case 'cum' % cummulative DVH
        for j = 1:n
            dvh(j) = sum(doseInVoi >= dvhPoints(j));
        end

    case 'diff' % differential DVH
        binning = (dvhPoints(2) - dvhPoints(1))/2;
        for j = 1:n % differential DVH
            dvh(j) = sum(dvhPoints(j) + binning > doseInVoi & doseInVoi > dvhPoints(j) - binning);
        end

end
dvh = dvh ./ numOfVoxels * 100;
end %eof getDVHPoints

repo: matRad


text:
function t_most_spuc(quiet, create_plots, create_pdfs, savedir)
%T_MOST_SPUC  Tests of single-period unit commitment optimizations
%
%   T_MOST_SPUC(QUIET, CREATE_PLOTS, CREATE_PDFS, SAVEDIR)
%   Can generate summary plots and save them as PDFs in a directory of
%   your choice.
%   E.g. t_most_spuc(0, 1, 1, '~/Downloads/spuc_plots')

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

if nargin < 4
    savedir = '.';              %% save in current working directory by default
    if nargin < 3
        create_pdfs = 0;        %% do NOT save plots to PDF files
        if nargin < 2
            create_plots = 0;   %% do NOT create summary plots of results
            if nargin < 1
                quiet = 0;      %% verbose by default
            end
        end
    end
end

solvers = {'CPLEX', 'GLPK', 'GUROBI', 'MOSEK', 'OT'};
fcn = {'cplex', 'glpk', 'gurobi', 'mosek', 'intlinprog'};
% solvers = {'OT'};
% fcn = {'intlinprog'};
% solvers = {'GUROBI'};
% fcn = {'gurobi'};
% solvers = {'GLPK'};
% fcn = {'glpk'};
ntests = 144;
t_begin(ntests*length(solvers), quiet);

if quiet
    verbose = 0;
else
    verbose = 0;
end
% verbose = 2;

casefile = 'ex_case3b';
mpopt = mpoption;
mpopt = mpoption(mpopt, 'out.gen', 1);
mpopt = mpoption(mpopt, 'verbose', verbose);
% mpopt = mpoption(mpopt, 'opf.violation', 1e-6, 'mips.gradtol', 1e-8, ...
%         'mips.comptol', 1e-8, 'mips.costtol', 1e-8);
mpopt = mpoption(mpopt, 'model', 'DC');
mpopt = mpoption(mpopt, 'sopf.force_Pc_eq_P0', 0);
mpopt = mpoption(mpopt, 'most.price_stage_warn_tol', 10);

%% solver options
if have_feature('cplex')
    %mpopt = mpoption(mpopt, 'cplex.lpmethod', 0);       %% automatic
    %mpopt = mpoption(mpopt, 'cplex.lpmethod', 1);       %% primal simplex
    mpopt = mpoption(mpopt, 'cplex.lpmethod', 2);       %% dual simplex
    %mpopt = mpoption(mpopt, 'cplex.lpmethod', 3);       %% network simplex
    %mpopt = mpoption(mpopt, 'cplex.lpmethod', 4);       %% barrier
    mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.mipgap', 0);
    mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.absmipgap', 0);
    mpopt = mpoption(mpopt, 'cplex.opts.threads', 2);
end
if have_feature('glpk')
    mpopt = mpoption(mpopt, 'glpk.opts.mipgap', 0);
    mpopt = mpoption(mpopt, 'glpk.opts.tolint', 1e-10);
    mpopt = mpoption(mpopt, 'glpk.opts.tolobj', 1e-10);
end
if have_feature('gurobi')
    %mpopt = mpoption(mpopt, 'gurobi.method', -1);       %% automatic
    %mpopt = mpoption(mpopt, 'gurobi.method', 0);        %% primal simplex
    mpopt = mpoption(mpopt, 'gurobi.method', 1);        %% dual simplex
    %mpopt = mpoption(mpopt, 'gurobi.method', 2);        %% barrier
    mpopt = mpoption(mpopt, 'gurobi.threads', 2);
    mpopt = mpoption(mpopt, 'gurobi.opts.MIPGap', 0);
    mpopt = mpoption(mpopt, 'gurobi.opts.MIPGapAbs', 0);
end
if have_feature('mosek')
    sc = mosek_symbcon;
    %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_FREE);            %% default
    %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_INTPNT);          %% interior point
    %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_PRIMAL_SIMPLEX);  %% primal simplex
    mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);     %% dual simplex
    %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_FREE_SIMPLEX);    %% automatic simplex
    %mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_X', 0);
    mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_NODE_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
    mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_ROOT_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
    mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_RELAX_INT', 1e-9);
    %mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_RELAX_INT', 0);
    mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_GAP', 0);
    mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_GAP', 0);
end
if have_feature('intlinprog')
    %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'interior-point');
    %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'active-set');
    %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'simplex');
    mpopt = mpoption(mpopt, 'linprog.Algorithm', 'dual-simplex');
    %mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'primal-simplex');
    mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'dual-simplex');
    mpopt = mpoption(mpopt, 'intlinprog.TolCon', 1e-9);
    mpopt = mpoption(mpopt, 'intlinprog.TolGapAbs', 0);
    mpopt = mpoption(mpopt, 'intlinprog.TolGapRel', 0);
    mpopt = mpoption(mpopt, 'intlinprog.TolInteger', 1e-6);
    %% next line is to work around a bug in intlinprog
    % (Technical Support Case #01841662)
    mpopt = mpoption(mpopt, 'intlinprog.LPPreprocess', 'none');
end
if ~verbose
    mpopt = mpoption(mpopt, 'out.all', 0);
end
% mpopt = mpoption(mpopt, 'out.all', -1);

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
[CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
    CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
    CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
    CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
    CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
    CT_MODCOST_X] = idx_ct;

%% load base case file
mpc = loadcase(casefile);
mpc.gencost(:, STARTUP) = 0;
mpc.gencost(:, SHUTDOWN) = 0;

%%-----  contingencies  -----
contab = loadgenericdata('ex_contab', 'array');
pp = [1-sum(contab(:,2)); contab(1,2); contab(2,2)];

xgd = loadxgendata('ex_xgd_uc', mpc);
[iwind, mpc, xgd] = addwind('ex_wind_uc', mpc, xgd);
mpc.reserves.zones = [mpc.reserves.zones 0];
mpc = scale_load(499, mpc, [], struct('scale', 'QUANTITY'));
mpc0 = mpc;
xgd0 = xgd;

%% data structures for results for plotting
if create_plots
    j = 1;
    Pg   = NaN(5, 7);
    Rp   = NaN(5, 7);
    Rm   = NaN(5, 7);
    lamP = NaN(3, 7);
    muF  = zeros(3, 7);
end

for s = 1:length(solvers)
    if ~have_feature(fcn{s})    %% check if we have the solver
        t_skip(ntests, sprintf('%s not installed', solvers{s}));
    else
        mpopt = mpoption(mpopt, 'opf.dc.solver', solvers{s});
        mpopt = mpoption(mpopt, 'most.solver', mpopt.opf.dc.solver);

%%-----  economic dispatch (no network)  -----
if verbose
    fprintf('\n\n');
    fprintf('--------------------------------------------\n');
    fprintf('-----  economic dispatch (no network)  -----\n');
    fprintf('--------------------------------------------\n');
end
t = sprintf('%s : economic dispatch (no network) : runopf ', solvers{s});
mpc = mpc0;
mpc.branch(:, RATE_A) = 0;
r = runuopf(mpc, mpopt);
t_ok(r.success, [t 'success']);
t_is(r.f, -488030, 7, [t 'f']);
t_is(r.gen(:, PG), [200; 199; 0; -499; 100], 7, [t 'Pg']);
t_is(r.gen(:, GEN_STATUS), [1; 1; 0; 1; 1], 7, [t 'u']);
t_is(r.bus(:, LAM_P), [30; 30; 30], 7, [t 'lam P']);

%% most
t = sprintf('%s : economic dispatch (no network) : most   ', solvers{s});
mpc = mpc0;
mpc.gen(1, GEN_STATUS) = 0;
mpc.gen(2, GEN_STATUS) = 0;
% xgd = xg
mpopt = mpoption(mpopt, 'most.dc_model', 0);
mdi = loadmd(mpc, [], xgd);
mdo = most(mdi, mpopt);
rr = mdo.flow(1,1,1).mpc;
t_ok(mdo.QP.exitflag > 0, [t 'success']);
t_is(mdo.QP.f, -488030, 7, [t 'f']);
t_is(rr.gen(:, PG), [200; 199; 0; -499; 100], 7, [t 'Pg']);
t_is(rr.gen(:, GEN_STATUS), [1; 1; 0; 1; 1], 7, [t 'u']);
% rr.gen(:, GEN_STATUS)
t_is(rr.bus(:, LAM_P), [30; 30; 30], 7, [t 'lam P']);
if create_plots
    Pg(:, j) = mdo.results.ExpectedDispatch;
    Rp(:, j) = 0;
    Rm(:, j) = 0;
    lamP(:, j) = rr.bus(:, LAM_P);
    j = j + 1;
end

%%-----  DC OPF  -----
if verbose
    fprintf('\n\n');
    fprintf('--------------------------------------------\n');
    fprintf('-----              DC OPF              -----\n');
    fprintf('--------------------------------------------\n');
end
t = sprintf('%s : DC OPF : runopf ', solvers{s});
mpc = mpc0;
% mpc = scale_load(500, mpc, [], struct('scale', 'QUANTITY'));
r = runuopf(mpc, mpopt);
t_ok(r.success, [t 'success']);
t_is(r.f, -486040, 7, [t 'f']);
t_is(r.gen(:, PG), [200; 0; 199; -499; 100], 7, [t 'Pg']);
t_is(r.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(r.bus(:, LAM_P), [40; 40; 40], 7, [t 'lam P']);
t_is(r.branch(:, MU_SF) + r.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow']);

%% most
t = sprintf('%s : DC OPF : most   ', solvers{s});
mpc = mpc0;
% mpc = scale_load(500, mpc, [], struct('scale', 'QUANTITY'));
mpopt = mpoption(mpopt, 'most.dc_model', 1);
mdi = loadmd(mpc, [], xgd);
mdo = most(mdi, mpopt);
rr = mdo.flow(1,1,1).mpc;
if verbose
    printpf(rr, [], mpopt);
end
t_ok(mdo.QP.exitflag > 0, [t 'success']);
t_is(mdo.QP.f, -486040, 7, [t 'f']);
t_is(rr.gen(:, PG), [200; 0; 199; -499; 100], 6, [t 'Pg']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
% rr.gen(:, GEN_STATUS)
t_is(rr.bus(:, LAM_P), [40; 40; 40], 7, [t 'lam P']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow']);
if create_plots
    Pg(:, j) = mdo.results.ExpectedDispatch;
    Rp(:, j) = 0;
    Rm(:, j) = 0;
    lamP(:, j) = rr.bus(:, LAM_P);
    muF(:, j)  = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
    j = j + 1;
end

%%-----  economic dispatch (w/reserves)  -----
if verbose
    fprintf('\n\n');
    fprintf('--------------------------------------------\n');
    fprintf('-----  economic dispatch (w/reserves)  -----\n');
    fprintf('--------------------------------------------\n');
end
t = sprintf('%s : economic dispatch (w/reserves) : runopf ', solvers{s});
mpc = mpc0;
mpc.branch(:, RATE_A) = 0;
mpc = toggle_reserves(mpc, 'on');
r = runuopf(mpc, mpopt);
t_ok(r.success, [t 'success']);
t_is(r.f, -486802, 5, [t 'f']);
t_is(r.gen(:, PG), [200; 139; 60; -499; 100], 7, [t 'Pg']);
t_is(r.gen(:, GEN_STATUS), [1; 1; 1; 1; 1], 7, [t 'u']);
t_is(r.bus(:, LAM_P), [32; 32; 32], 7, [t 'lam P']);
t_is(r.reserves.R, [0; 61; 89; 0; 0], 7, [t 'R']);
t_is(r.reserves.prc, [5; 5; 5; 0; 0], 7, [t 'reserve prc']);
t_is(r.reserves.mu.Pmax + r.gen(:, MU_PMAX), [7; 2; 0; 0; 32], 7, [t 'reserve muPmax']);

%% most
t = sprintf('%s : economic dispatch (w/reserves) : most   ', solvers{s});
mpc = mpc0;
% mpc = scale_load(350.8, mpc, [], struct('scale', 'QUANTITY'));
mpopt = mpoption(mpopt, 'most.dc_model', 0);
mdi = loadmd(mpc, [], xgd);
mdi.FixedReserves = mpc.reserves;
mdo = most(mdi, mpopt);
rr = mdo.flow(1,1,1).mpc;
if verbose
    printpf(rr, [], mpopt);
end
t_ok(mdo.QP.exitflag > 0, [t 'success']);
t_is(mdo.QP.f, -486802, 5, [t 'f']);
t_is(rr.gen(:, PG), [200; 139; 60; -499; 100], 7, [t 'Pg']);
t_is(rr.gen(:, GEN_STATUS), [1; 1; 1; 1; 1], 7, [t 'u']);
% rr.gen(:, GEN_STATUS)
t_is(rr.bus(:, LAM_P), [32; 32; 32], 7, [t 'lam P']);
t_is(rr.reserves.R, [0; 61; 89; 0; 0], 7, [t 'R']);
t_is(rr.reserves.prc, [5; 5; 5; 0; 0], 7, [t 'reserve prc']);
t_is(rr.reserves.mu.Pmax + rr.gen(:, MU_PMAX), [7; 2; 0; 0; 32], 7, [t 'reserve muPmax']);
if create_plots
    Pg(:, j) = mdo.results.ExpectedDispatch;
    Rp(:, j) = rr.reserves.R;
    Rm(:, j) = 0;
    lamP(:, j) = rr.bus(:, LAM_P);
    j = j + 1;
end

%%-----  DC OPF  -----
if verbose
    fprintf('\n\n');
    fprintf('--------------------------------------------\n');
    fprintf('-----        DC OPF (w/reserves)       -----\n');
    fprintf('--------------------------------------------\n');
end
t = sprintf('%s : DC OPF (w/reserves) : runopf ', solvers{s});
mpc = mpc0;
r = runopf_w_res(mpc, mpopt);
t_ok(r.success, [t 'success']);
t_is(r.f, -485656, 5, [t 'f']);
t_is(r.gen(:, PG), [156; 65; 178; -499; 100], 7, [t 'Pg']);
t_is(r.gen(:, GEN_STATUS), [1; 1; 1; 1; 1], 7, [t 'u']);
t_is(r.bus(:, LAM_P), [29; 40; 51], 7, [t 'lam P']);
t_is(r.branch(:, MU_SF) + r.branch(:, MU_ST), [0; 33; 0], 7, [t 'mu flow']);
t_is(r.reserves.R, [44; 100; 6; 0; 0], 7, [t 'R']);
t_is(r.reserves.prc, [5; 5; 5; 0; 0], 7, [t 'reserve prc']);
t_is(r.reserves.mu.Pmax + r.gen(:, MU_PMAX), [4; 0; 0; 0; 40], 7, [t 'reserve muPmax']);

%% most
t = sprintf('%s : DC OPF (w/reserves) : most   ', solvers{s});
mpc = mpc0;
% mpc = scale_load(350.8, mpc, [], struct('scale', 'QUANTITY'));
mpopt = mpoption(mpopt, 'most.dc_model', 1);
mdi = loadmd(mpc, [], xgd);
mdi.FixedReserves = mpc.reserves;
mdo = most(mdi, mpopt);
rr = mdo.flow(1,1,1).mpc;
t_ok(mdo.QP.exitflag > 0, [t 'success']);
t_is(mdo.QP.f, -485656, 5, [t 'f']);
t_is(rr.gen(:, PG), [156; 65; 178; -499; 100], 6, [t 'Pg']);
t_is(rr.gen(:, GEN_STATUS), [1; 1; 1; 1; 1], 7, [t 'u']);
% rr.gen(:, GEN_STATUS)
t_is(rr.bus(:, LAM_P), [29; 40; 51], 6, [t 'lam P']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 33; 0], 7, [t 'mu flow']);
t_is(rr.reserves.R, [44; 100; 6; 0; 0], 6, [t 'R']);
t_is(rr.reserves.prc, [5; 5; 5; 0; 0], 7, [t 'reserve prc']);
t_is(rr.reserves.mu.Pmax + rr.gen(:, MU_PMAX), [4; 0; 0; 0; 40], 7, [t 'reserve muPmax']);
if create_plots
    Pg(:, j) = mdo.results.ExpectedDispatch;
    Rp(:, j) = rr.reserves.R;
    Rm(:, j) = 0;
    lamP(:, j) = rr.bus(:, LAM_P);
    muF(:, j)  = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
    j = j + 1;
end

t = sprintf('%s : Secure DC OPF (w/cont,res,ramp) : c3sopf ', solvers{s});
mpc = mpc0;
mpc = scale_load(350, mpc, [], struct('scale', 'QUANTITY'));
xgd_table.colnames = {
    'PositiveActiveReservePrice', ...
            'PositiveActiveReserveQuantity', ...
                    'NegativeActiveReservePrice', ...
                            'NegativeActiveReserveQuantity', ...
                                    'PositiveActiveDeltaPrice', ...
                                            'NegativeActiveDeltaPrice', ...
};
xgd_table.data = [
    5       250     10      0       1e-9    1e-9;
    1e-8    100     2e-8    0       1e-9    1e-9;
    1.5     600     3       0       1e-9    1e-9;
    1e-8    800     2e-8    0       1e-9    1e-9;
    1e-8    200     2e-8    0       1e-9    1e-9;
];
% xgd = loadxgendata(xgd_table, mpc);
xgd.PositiveActiveReserveQuantity(2) = 100;
xgd.PositiveActiveReserveQuantity(4) = 800;
xgd.NegativeActiveReserveQuantity(:) = 0;
xgd.PositiveActiveReservePrice([1;3]) = [5; 1.5];
xgd.NegativeActiveReservePrice([1;3]) = [10; 3];

mpc1 = mpc;
mpc1.gen(2, GEN_STATUS) = 0;
r = c3sopf(mpc1, xgd_table.data, contab, mpopt);
% r
% r.opf_results
% r.opf_results.f
t_ok(r.success, [t 'success']);
t_is(r.opf_results.f, -340850, 4, [t 'f']);
rr = r.base;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [190; 0; 60; -350; 100], 6, [t 'Pg base']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [25; 25; 25], 7, [t 'lam P base']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow base']);
rr = r.cont(1);
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_ok(isempty(rr.gen), [t 'Pg 1']);
t_ok(isempty(rr.bus), [t 'lam P 1']);
t_ok(isempty(rr.branch), [t 'mu flow 1']);
rr = r.cont(2);
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [190; 0; 60; -300; 50], 5, [t 'Pg 2']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [0; 0; 40], 7, [t 'lam P 2']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 40], 7, [t 'mu flow 2']);

t = sprintf('%s : Secure DC OPF (w/cont,res,ramp) : most ', solvers{s});
mdi = loadmd(mpc, [], xgd, [], contab);
mdo = most(mdi, mpopt);
% mdo
% mdo.QP
% mdo.QP.f
t_ok(mdo.QP.exitflag > 0, [t 'success']);
t_is(mdo.QP.f, -340850, 5, [t 'f']);
rr = mdo.flow(1,1,1).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [190; 0; 60; -350; 100], 6, [t 'Pg base']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
% t_is(rr.bus(:, LAM_P), [20; 20; 20], 7, [t 'lam P base']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow base']);
if create_plots
    lamP(:, j) = rr.bus(:, LAM_P);
    muF(:, j)  = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr1 = rr;
rr = mdo.flow(1,1,2).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [190; 0; 60; -350; 100], 5, [t 'Pg 1']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr1.bus(:, LAM_P) + rr.bus(:, LAM_P), [25; 25; 25], 7, [t 'lam P base + lam P 1']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 1']);
if create_plots
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr = mdo.flow(1,1,3).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [190; 0; 60; -300; 50], 5, [t 'Pg 2']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [0; 0; 40], 7, [t 'lam P 2']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 40], 7, [t 'mu flow 2']);
if create_plots
    Pg(:, j) = mdo.results.ExpectedDispatch;
    Rp(:, j) = mdo.results.Rpp + mdo.results.Pc - mdo.results.ExpectedDispatch;
    Rm(:, j) = mdo.results.Rpm - mdo.results.Pc + mdo.results.ExpectedDispatch;
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
    j = j + 1;
end

t = sprintf('%s : Stochastic DC OPF (w/wind,res) : c3sopf ', solvers{s});
nt = 1;
nj = 3;
% transmat = transmatnormalpersistent(nt, nj);
transmat = {[0.158655253931457; 0.682689492137086; 0.158655253931457]};
profiles = getprofiles(uniformwindprofile(nt, nj), iwind);

pp = transmat{1};
%% contingency table
% label probty  type            row             column          chgtype newvalue
contab0 = [
    1   pp(1)   profiles.table  profiles.rows   profiles.col    profiles.chgtype  profiles.values(1)/profiles.values(2);
    3   pp(3)   profiles.table  profiles.rows   profiles.col    profiles.chgtype  profiles.values(3)/profiles.values(2);
];
mpc1 = mpc;
mpc1.gen(iwind, PMAX) = mpc1.gen(iwind, PMAX) * profiles.values(2);
mpc1.gen(3, GEN_STATUS) = 0;
wind = loadgenericdata('ex_wind_uc', 'struct', 'gen', 'wind');
offers = [xgd_table.data; wind.xgd_table.data(:, [3:8])];
% mpopt.verbose = 2;
% mpopt.out.all = -1;
r = c3sopf(mpc1, offers, contab0, mpopt);
t_ok(r.success, [t 'success']);
t_is(r.opf_results.f, -341928.605134, 5, [t 'f']);
rr = r.cont(1);
% rr.gen(:, PG)
% rr.bus(:, LAM_P)
% rr.branch(:, MU_SF) + rr.branch(:, MU_ST)
t_is(rr.gen(:, PG), [200; 150; 0; -350; 0], 7, [t 'Pg 1']);
t_is(rr.bus(:, LAM_P), [4.7596576; 4.7596576; 4.7596576], 7, [t 'lam P 1']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 1']);
rr = r.base;
t_is(rr.gen(:, PG), [200; 100; 0; -350; 50], 5, [t 'Pg 2']);
t_is(rr.bus(:, LAM_P), [20.4806848; 20.4806848; 20.4806848], 6, [t 'lam P 2']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2']);
rr = r.cont(2);
t_is(rr.gen(:, PG), [200; 65; 0; -350; 85], 5, [t 'Pg 3']);
t_is(rr.bus(:, LAM_P), [0; 0; 0], 6, [t 'lam P 3']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 3']);

t = sprintf('%s : Stochastic DC OPF (w/wind,res) : most ', solvers{s});
mdi = loadmd(mpc, transmat, xgd, [], [], profiles);
mdo = most(mdi, mpopt);
% mdo
% mdo.QP
% mdo.QP.f
t_ok(mdo.QP.exitflag > 0, [t 'success']);
t_is(mdo.QP.f, -341928.605134, 6, [t 'f']);
rr = mdo.flow(1,1,1).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 150; 0; -350; 0], 7, [t 'Pg 1']);
t_is(rr.gen(:, GEN_STATUS), [1; 1; 0; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [4.7596576; 4.7596576; 4.7596576], 7, [t 'lam P 1']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 1']);
if create_plots
    lamP(:, j) = rr.bus(:, LAM_P);
    muF(:, j)  = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr = mdo.flow(1,2,1).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 100; 0; -350; 50], 7, [t 'Pg 2']);
t_is(rr.gen(:, GEN_STATUS), [1; 1; 0; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [20.4806848; 20.4806848; 20.4806848], 6, [t 'lam P 2']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2']);
if create_plots
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr = mdo.flow(1,3,1).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 65; 0; -350; 85], 5, [t 'Pg 3']);
t_is(rr.gen(:, GEN_STATUS), [1; 1; 0; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [0; 0; 0], 6, [t 'lam P 3']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 3']);
if create_plots
    Pg(:, j) = mdo.results.ExpectedDispatch;
    Rp(:, j) = mdo.results.Rpp + mdo.results.Pc - mdo.results.ExpectedDispatch;
    Rm(:, j) = mdo.results.Rpm - mdo.results.Pc + mdo.results.ExpectedDispatch;
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
    j = j + 1;
end
% keyboard;

t = sprintf('%s : Secure Stochastic DC OPF (w/wind,cont,res,ramp) : most ', solvers{s});
mdi = loadmd(mpc, transmat, xgd, [], contab, profiles);
mdo = most(mdi, mpopt);
% mdo
% mdo.QP
% mdo.QP.f
t_ok(mdo.QP.exitflag > 0, [t 'success']);
t_is(mdo.QP.f, -338857.9224447, 4, [t 'f']);
% fprintf('%.7f\n', mdo.QP.f);

rr = mdo.flow(1,1,1).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 0; 150; -350; 0], 7, [t 'Pg 1 base']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
% t_is(rr.bus(:, LAM_P), [7.2115891; 7.2115891; 7.2115891], 6, [t 'lam P 1 base']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 2, [t 'mu flow 1 base']);
if create_plots
    lamP(:, j) = rr.bus(:, LAM_P);
    muF(:, j)  = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr1 = rr;
rr = mdo.flow(1,1,2).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 0; 150; -350; 0], 4, [t 'Pg 1 1']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
% t_is(rr.bus(:, LAM_P), [0.3807726; 0.3807726; 0.3807726], 7, [t 'lam P 1 1']);
t_is(rr1.bus(:, LAM_P) + rr.bus(:, LAM_P), [7.5923617; 7.5923617; 7.5923617], 7, [t 'lam P 1 1']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 1 1']);
if create_plots
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr = mdo.flow(1,1,3).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 0; 100; -300; 0], 4, [t 'Pg 1 2']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [0.2538484; 0.2538484; 6.3462102], 6, [t 'lam P 1 2']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 6.0923618], 6, [t 'mu flow 1 2']);
if create_plots
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr = mdo.flow(1,2,1).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 0; 100; -350; 50], 6, [t 'Pg 2 base']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [24.5768217; 24.5768217; 24.5768217], 6, [t 'lam P 2 base']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2 base']);
if create_plots
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr = mdo.flow(1,2,2).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 0; 100; -350; 50], 6, [t 'Pg 2 1']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [1.6384548; 1.6384548; 1.6384548], 6, [t 'lam P 2 1']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2 1']);
if create_plots
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr = mdo.flow(1,2,3).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 0; 60; -300; 40], 6, [t 'Pg 2 2']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [0; 0; 27.3075797], 5, [t 'lam P 2 2']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 27.3075797], 6, [t 'mu flow 2 2']);
if create_plots
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr = mdo.flow(1,3,1).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 0; 60; -350; 90], 6, [t 'Pg 3 base']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [0; 0; 0], 6, [t 'lam P 3 base']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 5, [t 'mu flow 3 base']);
if create_plots
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr = mdo.flow(1,3,2).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 0; 60; -350; 90], 6, [t 'Pg 3 1']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [0; 0; 0], 6, [t 'lam P 3 1']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 3 1']);
if create_plots
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
end

rr = mdo.flow(1,3,3).mpc;
% fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
% fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
% fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
t_is(rr.gen(:, PG), [200; 0; 60; -300; 40], 4, [t 'Pg 3 2']);
t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
t_is(rr.bus(:, LAM_P), [0; 0; 6.3462102], 5, [t 'lam P 3 2']);
t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 6.3462102], 6, [t 'mu flow 3 2']);
% keyboard;
end

if create_plots
    create_plots = 0;   %% don't do them again

    Pg(:, j) = mdo.results.ExpectedDispatch;
    Rp(:, j) = mdo.results.Rpp + mdo.results.Pc - mdo.results.ExpectedDispatch;
    Rm(:, j) = mdo.results.Rpm - mdo.results.Pc + mdo.results.ExpectedDispatch;
    lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
    muF(:, j)  = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
%    R(4:5, :) = NaN;
    R = Rp + Rm;

    labels = {'Economic Dispatch'; 'DC OPF'; 'Economic Dispatch w/reserves'; 'DC OPF w/reserves'; 'secure DC OPF'; 'stochastic DC OPF'; 'secure stochastic DC OPF'};

    figure(1);

    subplot(4, 1, 1);
    bar(abs(Pg([1:3 5],:)'),'stacked');
    title('Generation');
    ylabel('MW');
    h = gca;
    h.XTickLabel = labels;
    h.XTickLabelRotation = 20;

    subplot(4, 1, 2);
    bar(abs(R([1:3 5],:)'),'stacked');
    title('Reserves');
    ylabel('MW');
    h = gca;
    h.XTickLabel = labels;
    h.XTickLabelRotation = 20;

    subplot(4, 1, 3);
    plot(lamP');
    title('Nodal Prices');
    ylabel('$/MW');
    h = gca;
    h.XTickLabel = {'', labels{:}, ''}';
    h.XTickLabelRotation = 20;
    h = [0 8 0 100];
    axis(h);

    subplot(4, 1, 4);
    plot(muF');
    title('Flow Constraint Shadow Prices');
    ylabel('$/MW');
    h = gca;
    h.XTickLabel = {'', labels{:}, ''}';
    h.XTickLabelRotation = 20;
    h = [0 8 0 60];
    axis(h);

    if create_pdfs
        fname = 'single-period-uc';
        h = gcf;
        set(h, 'PaperSize', [11 8.5]);
        set(h, 'PaperPosition', [0.25 0.25 10.5 8]);
        print('-dpdf', fullfile(savedir, fname));
    end

    for j = 1:7;
        figure(j+1);
        if create_pdfs
            fname = sprintf('single-period-uc-%d', j);
        else
            fname = '';
        end
        plot_case(labels{j}, Pg([1:3 5], j), Rp([1:3 5], j), Rm([1:3 5], j), lamP(:, j), muF(:, j), 250, 100, savedir, fname);
    end
end
end

t_end;

function h = plot_case(label, Pg, Rp, Rm, lamP, muF, maxq, maxp, mypath, fname)

subplot(1, 3, 1);
h = bar([Pg-Rm Rm Rp], 'stacked');
set(h(2), 'FaceColor', [0 0.35 0.33]);
ah1 = gca;
title('Generation & Reserves');
ylabel('MW');
xlabel('Gen');
set(gca, 'FontSize', 14);

if nargin < 6
    maxq = ah1.YLim(2);
end
ah1.YLim(2) = maxq;
ah1.YLim(1) = 0;

subplot(1, 3, 2);
bar(lamP);
ah3 = gca;
title('Nodal Prices');
ylabel('$/MW');
xlabel('Bus');
set(gca, 'FontSize', 14);

subplot(1, 3, 3);
bar(muF);
ah4 = gca;
title('Flow Constraint Shadow Prices');
ylabel('$/MW');
xlabel('Line');
set(gca, 'FontSize', 14);

if nargin < 7
    maxp = max(ah3.YLim(2), ah4.YLim(2));
end
ah3.YLim(1) = 0;
ah4.YLim(1) = 0;
ah3.YLim(2) = maxp;
ah4.YLim(2) = maxp;

[ax,h] = suplabel(label, 't');
set(h, 'FontSize', 18)

if nargin > 7 && ~isempty(fname)
    h = gcf;
    set(h, 'PaperSize', [11 8.5]);
    set(h, 'PaperPosition', [0.25 0.25 10.5 8]);
    print('-dpdf', fullfile(mypath, fname));
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%% This function will generate the ODE function based on calculation result of the sybolic expression
% Last Update: 2019/04/22
% Coded By: K
%(Modified so that matrix demension match)

function Generate_ODE_RHS(Sindy_ODEs,var_num_state,var_num_control)
z=sym('z',[1,var_num_state]);
u=sym('u',[1,var_num_control]);
syms t
f= matlabFunction(Sindy_ODEs,'File','Sindy_ODE_RHS','Optimize',true,'Vars',{t,z,u});

author: dynamicslab

repo: SINDy-PI


text:
function [F,V,C]=delaunayZip(F1,V1,F2,V2,inputStruct)

% function [Fn]=delaunayZip(F1,V1,F2,V2,inputStruct)
%-------------------------------------------------------------------------
%
%
% Change log:
% 2018/05/09: Added growing of region size if error occurs
% 2018/05/09: Added stepwise removal of last triangle since it does not
% include next points.
% 2018/11/02: Added self-triangulation as an option (default on)
% 2018/11/02: Added walk-back for paths so that one will not be too long
% 2018/11/02: Fixed bug in relation to one of the paths "running out of
% points". Remainder is triangulated properly.
%-------------------------------------------------------------------------

%%

maxD=max([patchEdgeLengths(F1,V1);patchEdgeLengths(F2,V2)]);

defaultInputStruct.ind1=[];
defaultInputStruct.ind2=[];
defaultInputStruct.distLocal=2*maxD;
defaultInputStruct.startInd=[];
defaultInputStruct.plotOn=0;
defaultInputStruct.selfTriangulate=1; %Option to self-triangulate surfaces first
defaultInputStruct.angleThreshold=(60/180)*pi; %Angular threshold for self-triangulation
[inputStruct]=structComplete(inputStruct,defaultInputStruct,0); %Complement provided with default if missing or empty

ind1=inputStruct.ind1;
ind2=inputStruct.ind2;
distLocal=inputStruct.distLocal;
startInd=inputStruct.startInd;
selfTriangulate=inputStruct.selfTriangulate;
angleThreshold=inputStruct.angleThreshold;
plotOn=inputStruct.plotOn;

%% Self-triangulate if needed

if selfTriangulate==1
    if ind1(1)==ind1(end)
        isClosedLoop=1;
    else
        isClosedLoop=0;
    end
    numFacesInitial_1=size(F1,1);
    [F1,V1,ind1]=triSurfSelfTriangulateBoundary(F1,V1,ind1,angleThreshold,isClosedLoop);
    numFaces_1=size(F1,1);
    numFacesInitial_2=size(F2,1);
    [F2,V2,ind2]=triSurfSelfTriangulateBoundary(F2,V2,ind2,angleThreshold,isClosedLoop);
    numFaces_2=size(F2,1);

    L1=false(size(F1,1),1);
    if numFaces_1>numFacesInitial_1
        L1(end-(numFaces_1-numFacesInitial_1-1):end)=1;
    end
    L2=false(size(F2,1),1);
    if numFaces_2>numFacesInitial_2
        L2(end-(numFaces_2-numFacesInitial_2-1):end)=1;
    end
else
    L1=zeros(size(F1,1),1);
    L2=zeros(size(F2,1),1);
end

%%
[F,V,C]=joinElementSets({F1,F2},{V1,V2});
[~,~,Nv]=patchNormal(F,V);

%%

ind2=ind2+size(V1,1);
startInd(2)=startInd(2)+size(V1,1);

%%
%Create edges list
E=[ind1(1:end-1) ind1(2:end); ind2(1:end-1) ind2(2:end)];

%%

logicNotUsed=ismember((1:1:size(V,1))',E); %Logic to keep track of points that are used

%Initialize other parameters
Fn=[]; %Faces
Cn=[]; %"Colors"=step count for face group
c=1; %While loop counter variable

indGroup=startInd(:); %Initialize current group
numGroup=numel(indGroup); %Initialize current number of members of the current group
numGroupPrevious=0;

if plotOn==1
    h1=[]; %Initiate empty plot handle
    markerSize=25;

    hf=cFigure; hold on;
    gpatch(F1(~L1,:),V1,'r','none',0.2);
    gpatch(F1(L1,:),V1,'rw','k',1);
    gpatch(F2(~L2,:),V2,'b','none',0.2);
    gpatch(F2(L2,:),V2,'bw','k',1);
    plotV(V(ind1,:),'r.-','LineWidth',1,'MarkerSize',10);
    plotV(V(ind2,:),'b.-','LineWidth',1,'MarkerSize',10);

    axisGeom;
    camlight headlight;
    colormap(gjet(250));
    colorbar;
    drawnow;
end

%Turn off Delaunay warning as this case is handled properly
warning('off','MATLAB:delaunayTriangulation:ConsConsSplitWarnId');

%%
while 1

    numGroupStep=1;
    lastTry=0;
    while 1
        try
            if plotOn==1 %%Plot if plotting is on
                delete(h1); h1=[];
            end

            if plotOn==1 %%Plot if plotting is on
                figure(hf);
                h1(end+1)=plotV(V(startInd,:),'y.','MarkerSize',markerSize);
                drawnow
            end

            %Grow the current region
            while 1 %Loop to form local group (stuff attached to current group within a given distance)
                logicMember=any(ismember(E,indGroup),2); %Logic for all edges touching the current groupt
                E_sub=E(logicMember,:); %The subset of touching edges

                indGroup=unique([indGroup; E_sub(:)]); %Grow group with point indices in the edges that are touching

                d1=sqrt(sum((V(indGroup,:)-V(startInd(1)*ones(numel(indGroup),1),:)).^2,2));
                d2=sqrt(sum((V(indGroup,:)-V(startInd(2)*ones(numel(indGroup),1),:)).^2,2));
                D=min(d1,d2);

                logicKeep= (D<distLocal) & (logicNotUsed(indGroup));
                indGroup=indGroup(logicKeep); %Remove points that are too far
                if numGroup==numel(indGroup) %Compare current group size to previous step
                    break %break while loop if the group is no longer growing
                end
                numGroup=numel(indGroup); %Get new current group size
            end

            %Get curve start and end points
            logicMember=all(ismember(E,indGroup),2); %Logic for all edges touching the current group
            E_sub=E(logicMember,:); %The subset of touching edges

            [~,~,~,vCount]=cunique(E_sub); %Get vertex occurance counts
            indEndPoints=unique(E_sub(vCount==1));
            indEndPoints1=indEndPoints(ismember(indEndPoints,ind1));
            indEndPoints2=indEndPoints(ismember(indEndPoints,ind2));

            %Compose sub-curves
            logicSub1=all(ismember(E_sub,ind1),2);
            %             if ~any(logicSub1)
            %                 warning('None of the current edges nodes are a member of ind1');
            %             end
            E_sub1=E_sub(logicSub1,:);
            [indListSub1]=edgeListToCurve(E_sub1);
            indListSub1=indListSub1(:);

            logicSub2=all(ismember(E_sub,ind2),2);
            %             if ~any(logicSub2)
            %                 warning('None of the current edges nodes are a member of ind2');
            %             end
            E_sub2=E_sub(logicSub2,:);
            [indListSub2]=edgeListToCurve(E_sub2);
            indListSub2=indListSub2(:);

            %Alter segments so they can be merged into closed loop curve
            if isempty(Fn) %First iteration
                d=sqrt(sum((V([indListSub2(1) indListSub2(end)],:)-V(indListSub1(1)*ones(1,2),:)).^2,2));
                [~,indClosest]=min(d);
                if indClosest~=2 %If the start is closest to the other start
                    indListSub2=flipud(indListSub2);
                end
            else
                indStart1=find(ismember(indListSub1,startInd));
                if indStart1~=1
                    indListSub1=flipud(indListSub1);
                end
                indStart2=find(ismember(indListSub2,startInd));
                if indStart2==1
                    indListSub2=flipud(indListSub2);
                end
            end

            if any(ismember(indGroup,ind1))==0
                startInd
                ind1
                indGroup
            end

            %             if any(ismember(indGroup,ind1)) && any(ismember(indGroup,ind2))
            %             else
            %                 %One of the two sets contains only 1 member, triangulate
            %                 %remaining points manually
            %
            %                 %-----------------------------------
            %                 if nnz(ismember(indGroup,ind1))==0
            %                     indLeft=flipud(indListSub2);
            %                     sKeep=startInd(ismember(startInd,ind1));
            %                 else%if nnz(ismember(indGroup,ind2))==0
            %                     indLeft=indListSub1;
            %                     sKeep=startInd(ismember(startInd,ind2));
            %                 end
            %
            %                 f=[];
            %                 for qStep=1:1:numel(indLeft)-1
            %                     sVary=indLeft(1);
            %                     f=[f; [sKeep sVary indLeft(2)]];
            %                     indLeft=indLeft(2:end); %First is part of startInd
            %                 end
            %             end

            %Shorten curves if they stick out to much wrt end of other curve
            d=sqrt(sum((V(indListSub2,:)-V(indListSub1(end)*ones(1,numel(indListSub2)),:)).^2,2));
            [~,indClosest]=min(d);
            if indClosest~=1 && indClosest~=numel(indListSub2)
                indListSub2=indListSub2(indClosest:end);
            end

            d=sqrt(sum((V(indListSub1,:)-V(indListSub2(1)*ones(1,numel(indListSub1)),:)).^2,2));
            [~,indClosest]=min(d);
            if indClosest~=numel(indListSub1) && indClosest~=1
                indListSub1=indListSub1(1:indClosest);
            end

            %Create closed curve
            if nnz(ismember(indGroup,ind1))==0
                indListSub=[startInd(ismember(startInd,ind1)); indListSub2; startInd(ismember(startInd,ind1))]; %Closed loop
            elseif nnz(ismember(indGroup,ind2))==0
                indListSub=[startInd(ismember(startInd,ind2)); indListSub1; startInd(ismember(startInd,ind2))]; %Closed loop
            else
                indListSub=[indListSub1;indListSub2;indListSub1(1)]; %Closed loop
            end

            indGroup=indListSub(1:end-1);

            if plotOn==1 %%Plot if plotting is on
                figure(hf);
                if ~isempty(indEndPoints1)
                    h1(end+1)=plotV(V(indEndPoints1,:),'r.','MarkerSize',markerSize);
                end
                if ~isempty(indEndPoints2)
                    h1(end+1)=plotV(V(indEndPoints2,:),'b.','MarkerSize',markerSize);
                end
                if ~isempty(indListSub1)
                    h1(end+1)=plotV(V(indListSub1,:),'r-','LineWidth',3);
                end
                if ~isempty(indListSub2)
                    h1(end+1)=plotV(V(indListSub2,:),'b-','LineWidth',3);
                end
                if ~isempty(indListSub)
                    h1(end+1)=plotV(V(indListSub,:),'k-','LineWidth',2);
                end
                drawnow
            end

            %Rotate current point set
            V_now=V(indListSub(1:end-1),:); %Current closed curve coordinate set
            [R]=pointSetPrincipalDir(V_now); %Fit local coordinate system with 3rd direction pointing outward of local planar-ish region
            V_now_R=V_now*R; %Rotate coordinate set to prepare for 2D Delaunay based triangulation

            %Do 2D Delaunay triangulation
            DT_contraints=[(1:numel(indGroup))' ([2:numel(indGroup) 1])']; %Constraints are edges forming boundary
            DT = delaunayTriangulation(V_now_R(:,[1 2]),DT_contraints); % Initial triangulation
            f=DT.ConnectivityList; %Get faces set
            L = isInterior(DT); %Remove faces not inside region
            f=f(L,:); %Faces excluding external faces (only keep those inside constraint edges)
            f=indGroup(f); %Change indices to overall system

            %Remove first and last triangles if there are more than 3 triangles, this
            %will improve triangulation quality as these triangles were not fully
            %embedded in the point set.
            %                 if size(f,1)>2
            %                     logicKeep=~(any(ismember(f,indEndPoints1(~ismember(indEndPoints1,Fn))),2)...
            %                         & any(ismember(f,indEndPoints2(~ismember(indEndPoints2,Fn))),2))...
            %                         | any(ismember(f,startInd),2);
            %                     f=f(logicKeep,:);
            %
            %                     [~,~,IND_FF]=tesIND(f,V);
            %                     logicNeighbours=sum(IND_FF>1,2)>1;
            %                     logicNeighbours(any(ismember(f,startInd),2))=1;
            %                     f=f(logicNeighbours,:);
            %                 end

            break

        catch ME
            if lastTry==1
                rethrow(ME)
            end

            if numGroupStep==numGroupPrevious
                lastTry=1;
            end
            distLocal=distLocal+maxD;
            warning(['Failed using current step size, increasing to: ',num2str(distLocal)]);
            numGroupPrevious=numGroupStep;
        end
    end
    distLocal=inputStruct.distLocal; %reset

    if dot(mean(patchNormal(f,V),1),mean(Nv(indListSub(1:end-1),:),1))<0
        f=fliplr(f);
    end

    if plotOn==1 %%Plot if plotting is on
        figure(hf);
        V_now_mean=mean(V_now); %Mean of coordinate set
        h1(end+1)=plotV(V_now_mean,'kx','MarkerSize',markerSize);
        gpatch(f,V,c*ones(size(f,1),1),'k',1);
        drawnow
    end

    %Collect faces and color data
    Fn=[Fn;f];
    Cn=[Cn;c*ones(size(f,1),1)];

    E_Fn=sort(patchEdges(Fn,0),2);
    ind_E_Fn= reshape(sub2indn(size(V,1)*ones(1,2),E_Fn),size(Fn));

    Es=sort(E,2);
    ind_E= sub2indn(size(V,1)*ones(1,2),Es);

    logicEdgesUsed=ismember(ind_E,ind_E_Fn);

    E_sub=E(~logicEdgesUsed,:);
    indUnusedPoints=unique(E_sub);
    logicNotUsed=false(size(logicNotUsed));
    logicNotUsed(indUnusedPoints)=1;

    if any(logicNotUsed)==0
        break
    end

    %Get curve start and end points
    [~,~,~,vCount]=cunique(E_sub); %Get vertex occurance counts
    indEndPoints=unique(E_sub(vCount==1));
    indEndPoints1=indEndPoints(ismember(indEndPoints,ind1));
    indEndPoints2=indEndPoints(ismember(indEndPoints,ind2));

    e=patchEdges(Fn,1);
    edgeEndPoints=e(any(ismember(e,indEndPoints1),2)&any(ismember(e,indEndPoints2),2),:);
    logicMemberOfLast=any(ismember(edgeEndPoints,f),2);
    edgeEndPoints=edgeEndPoints(logicMemberOfLast,:);
    if ~isempty(edgeEndPoints)
        startInd=edgeEndPoints(1,:);
    else
        startInd=[indListSub1(end) indListSub2(1)];
    end

    c=c+1;
end

F=[F;Fn];
C=[C; max(C(:))+Cn];

%Turn Delaunay warning back on
warning('on','MATLAB:delaunayTriangulation:ConsConsSplitWarnId');

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
function Y = vl_nnloss(X,c,dzdy,varargin)

% --------------------------------------------------------------------
% pixel-level L2 loss
% --------------------------------------------------------------------
if nargin <= 2 || isempty(dzdy)
    t = ((X-c).^2)/2;
    Y = sum(t(:))/size(X,4); % reconstruction error per sample;
else
    Y = bsxfun(@minus,X,c).*dzdy;
end

author: cszn

repo: DnCNN

save_path: github-repos/MATLAB/cszn-DnCNN


text:
function print_same_line(string)

% This code has been taken from Les Schaffer on comp.soft-sys.matlab
% Hope he will not mind having this sudden fame.
% Send any comment to lerouxni@iro.umontreal.ca

% small function to print a string starting from the beginning of the
% actual line

%%  the ascii code for starting an escape sequence
esc_code = 27;

%% clear to end of line
%% ESC [2K
er = [ char(esc_code) '[' '2' 'K' ];

%% goto beginning of line
%% ESC [1G
bl = [ char(esc_code) '[' '1' 'G' ];

new_string = strcat(bl, er, string);
fprintf(1, new_string);


text:
clearvars;
close all;
% clc;

rng default;

M = 300;
N = 1000;
K = 60;

A = randn(M, N);
x0 = zeros(N,1);

permutation = randperm(N);
indices = permutation(1:K);
% non-zero entries
x0(indices) = randn(K,1);

b = A * x0;

% orthonormalize the rows of A
% perform QR decomposition of rows of A
[Q, R] = qr(A', 0);
% keep the orthogonal rows
A = Q';
% change b accordingly to the new coordinates
b = R'\b;

options.verbose = 0;
options.max_iterations = 200;
options.tolerance = 5e-4;
solver = spx.pursuit.single.L1_ADMM_YZ(A, options);
x = solver.solve_bp(b);

r = x - x0;
max_diff = max(abs(r));
rel_error = norm(x-x0) /norm(x0);
iterations = solver.details.iterations(1);
elapsed_time = solver.details.elapsed_times(1);
fprintf('Iterations: %d, Relative error: %e, max diff: %.4f, time: %.4f seconds\n', ...
    iterations, rel_error, max_diff, elapsed_time);

print = 1;
if print
    subplot(211);
    stem(x0, '.');
    subplot(212);
    stem(x, '.');
    figure;
    iterations = solver.details.iterations(1);
    iters = 1:iterations;
    primal_objectives = solver.details.primal_objectives(iters, 1);
    dual_objectives = solver.details.dual_objectives(iters, 1);
    plot(iters, primal_objectives);
    hold on;
    plot(iters, dual_objectives);
    xlabel('Iterations');
    ylabel('Objective Value');
    legend({'Primal Objective', 'Dual Objective'});
    grid on;
end

author: indigits

repo: sparse-plex


text:
% =========================================================================
% Test code for Super-Resolution Convolutional Neural Networks (SRCNN)
%
% Reference
%   Chao Dong, Chen Change Loy, Kaiming He, Xiaoou Tang. Learning a Deep Convolutional Network for Image Super-Resolution,
%   in Proceedings of European Conference on Computer Vision (ECCV), 2014
%
%   Chao Dong, Chen Change Loy, Kaiming He, Xiaoou Tang. Image Super-Resolution Using Deep Convolutional Networks,
%   arXiv:1501.00092
%
% Chao Dong
% IE Department, The Chinese University of Hong Kong
% For any question, send email to ndc.forward@gmail.com
% =========================================================================

close all;
clear all;

%% read ground truth image
im  = imread('Set5\butterfly_GT.bmp');
%im  = imread('Set14\zebra.bmp');

%% set parameters
up_scale = 3;
model = 'model\9-5-5(ImageNet)\x3.mat';
% up_scale = 3;
% model = 'model\9-3-5(ImageNet)\x3.mat';
% up_scale = 3;
% model = 'model\9-1-5(91 images)\x3.mat';
% up_scale = 2;
% model = 'model\9-5-5(ImageNet)\x2.mat';
% up_scale = 4;
% model = 'model\9-5-5(ImageNet)\x4.mat';

%% work on illuminance only
if size(im,3)>1
    im = rgb2ycbcr(im);
    im = im(:, :, 1);
end
im_gnd = modcrop(im, up_scale);
im_gnd = single(im_gnd)/255;

%% bicubic interpolation
im_l = imresize(im_gnd, 1/up_scale, 'bicubic');
im_b = imresize(im_l, up_scale, 'bicubic');

%% SRCNN
im_h = SRCNN(model, im_b);

%% remove border
im_h = shave(uint8(im_h * 255), [up_scale, up_scale]);
im_gnd = shave(uint8(im_gnd * 255), [up_scale, up_scale]);
im_b = shave(uint8(im_b * 255), [up_scale, up_scale]);

%% compute PSNR
psnr_bic = compute_psnr(im_gnd,im_b);
psnr_srcnn = compute_psnr(im_gnd,im_h);

%% show results
fprintf('PSNR for Bicubic Interpolation: %f dB\n', psnr_bic);
fprintf('PSNR for SRCNN Reconstruction: %f dB\n', psnr_srcnn);

figure, imshow(im_b); title('Bicubic Interpolation');
figure, imshow(im_h); title('SRCNN Reconstruction');

%imwrite(im_b, ['Bicubic Interpolation' '.bmp']);
%imwrite(im_h, ['SRCNN Reconstruction' '.bmp']);

author: thomas-koehler

repo: SupER


text:
% TENSORLAB
% Version 2.01, 2014-02-05
%
% BLOCK TERM DECOMPOSITION
% Algorithms
%    btd_minf     - BTD by unconstrained nonlinear optimization.
%    btd_nls      - BTD by nonlinear least squares.
% Initialization
%    btd_rnd      - Pseudorandom initialization for BTD.
% Utilities
%    btdgen       - Generate full tensor given a BTD.
%    btdres       - Residual of a BTD.
%
% CANONICAL POLYADIC DECOMPOSITION
% Algorithms
%    cpd          - Canonical polyadic decomposition.
%    cpd_als      - CPD by alternating least squares.
%    cpd_minf     - CPD by unconstrained nonlinear optimization.
%    cpd_nls      - CPD by nonlinear least squares.
%    cpd3_sd      - CPD by simultaneous diagonalization.
%    cpd3_sgsd    - CPD by simultaneous generalized Schur decomposition.
% Initialization
%    cpd_gevd     - CPD by a generalized eigenvalue decomposition.
%    cpd_rnd      - Pseudorandom initialization for CPD.
% Line and plane search
%    cpd_aels     - CPD approximate enhanced line search.
%    cpd_els      - CPD exact line search.
%    cpd_eps      - CPD exact plane search.
%    cpd_lsb      - CPD line search by Bro.
% Utilities
%    cpderr       - Errors between factor matrices in a CPD.
%    cpdgen       - Generate full tensor given a polyadic decomposition.
%    cpdres       - Residual of a polyadic decomposition.
%    rankest      - Estimate rank.
%
% COMPLEX OPTIMIZATION
% Nonlinear least squares
%    nls_gncgs    - Nonlinear least squares by Gauss-Newton with CG-Steihaug.
%    nls_gndl     - Nonlinear least squares by Gauss-Newton with dogleg trust region.
%    nls_lm       - Nonlinear least squares by Levenberg-Marquardt.
%    nlsb_gndl    - Bound-constrained NLS by projected Gauss-Newton dogleg TR.
% Unconstrained nonlinear optimization
%    minf_lbfgs   - Minimize a function by L-BFGS with line search.
%    minf_lbfgsdl - Minimize a function by L-BFGS with dogleg trust region.
%    minf_ncg     - Minimize a function by nonlinear conjugate gradient.
% Utilities
%    deriv        - Approximate gradient and Jacobian.
%    ls_mt        - Strong Wolfe line search by More-Thuente.
%    mpcg         - Modified preconditioned conjugate gradients method.
%
% LOW MULTILINEAR RANK APPROXIMATION
% Algorithms
%    lmlra        - Low multilinear rank approximation.
%    lmlra_hooi   - LMLRA by higher-order orthogonal iteration.
%    lmlra_minf   - LMLRA by unconstrained nonlinear optimization.
%    lmlra_nls    - LMLRA by nonlinear least squares.
%    lmlra3_dgn   - LMLRA by a differential-geometric Newton method.
%    lmlra3_rtr   - LMLRA by a Riemannian trust region method.
%    mlsvd        - (Truncated) multilinear singular value decomposition.
% Initialization
%    lmlra_aca    - LMLRA by adaptive cross-approximation.
%    lmlra_rnd    - Pseudorandom initialization for LMLRA.
% Utilities
%    lmlraerr     - Errors between factor matrices in a LMLRA.
%    lmlragen     - Generate full tensor given a core tensor and factor matrices.
%    lmlrares     - Residual of a LMLRA.
%    mlrank       - Multilinear rank.
%    mlrankest    - Estimate multilinear rank.
%
% STRUCTURED DATA FUSION
% Algorithms
%   sdf_minf      - Structured data fusion by unconstrained nonlinear optimization.
%   sdf_nls       - Structured data fusion by nonlinear least squares.
% Structure
%   struct_abs        - Absolute value.
%   struct_band       - Band matrix.
%   struct_cell2mat   - Convert the contents of a cell array into a matrix.
%   struct_conj       - Complex conjugate.
%   struct_ctranspose - Complex conjugate transpose.
%   struct_diag       - Diagonal matrix.
%   struct_gram       - Gramian matrix.
%   struct_hankel     - Hankel matrix.
%   struct_inv        - Matrix inverse.
%   struct_invsqrtm   - Matrix inverse square root.
%   struct_invtransp  - Matrix inverse transpose.
%   struct_LL1        - Structure of third factor matrix in a rank-(Lr,Lr,1) BTD.
%   struct_log        - Natural logarithm.
%   struct_matvec     - Matrix-vector and matrix-matrix product.
%   struct_nonneg     - Nonnegative array.
%   struct_normalize  - Normalize columns to unit norm.
%   struct_orth       - Rectangular matrix with orthonormal columns.
%   struct_plus       - Plus.
%   struct_poly       - Matrix with columns as polynomials.
%   struct_power      - Array power.
%   struct_rational   - Matrix with columns as rational functions.
%   struct_rbf        - Matrix with columns as sums of Gaussian RBF kernels.
%   struct_sigmoid    - Constrain array elements to an interval.
%   struct_sqrt       - Square root.
%   struct_sum        - Sum of elements.
%   struct_times      - Times.
%   struct_toeplitz   - Toeplitz matrix.
%   struct_transpose  - Transpose.
%   struct_tridiag    - Tridiagonal matrix.
%   struct_tril       - Lower triangular matrix.
%   struct_triu       - Upper triangular matrix.
%   struct_vander     - Vandermonde matrix.
%
% UTILITIES
% Clustering
%    gap          - Optimal clustering based on the gap statistic.
%    kmeans       - Cluster multivariate data using the k-means++ algorithm.
% Polynomials
%    polymin      - Minimize a polynomial.
%    polymin2     - Minimize bivariate and real polyanalytic polynomials.
%    polyval2     - Evaluate bivariate and univariate polyanalytic polynomials.
%    polysol2     - Solve a system of two bivariate polynomials.
%    ratmin       - Minimize a rational function.
%    ratmin2      - Minimize bivariate and real polyanalytic rational functions.
% Statistics
%    cum4         - Fourth-order cumulant tensor.
%    scov         - Shifted covariance matrices.
% Tensors
%    dotk         - Dot product in K-fold precision.
%    fmt          - Format data set.
%    frob         - Frobenius norm.
%    ful          - Convert formatted data set to an array.
%    kr           - Khatri-Rao product.
%    kron         - Kronecker product.
%    mat2tens     - Tensorize a matrix.
%    mtkronprod   - Compute a matricized tensor Kronecker product.
%    mtkrprod     - Compute a matricized tensor Khatri-Rao product.
%    noisy        - Generate a noisy version of a given array.
%    sumk         - Summation in K-fold precision.
%    tens2mat     - Matricize a tensor.
%    tens2vec     - Vectorize a tensor.
%    tmprod       - Mode-n tensor-matrix product.
%    vec2tens     - Tensorize a vector.
% Visualization
%    slice3       - Visualize a third-order tensor with slices.
%    spy3         - Visualize a third-order tensor's sparsity pattern.
%    surf3        - Visualize a third-order tensor with surfaces.
%    voxel3       - Visualize a third-order tensor with voxels.

author: andrewssobral

repo: lrslibrary


text:
function c = rdivide(a,b)
	checkinputs(a,b);
	c = uint64rdivide(uint64(a),uint64(b));
end