text:
%DEMEV1	Demonstrate Bayesian regression for the MLP.
%
%	Description
%	The problem consists an input variable X which sampled from a
%	Gaussian distribution, and a target variable T generated by computing
%	SIN(2*PI*X) and adding Gaussian noise. A 2-layer network with linear
%	outputs is trained by minimizing a sum-of-squares error function with
%	isotropic Gaussian regularizer, using the scaled conjugate gradient
%	optimizer. The hyperparameters ALPHA and BETA are re-estimated using
%	the function EVIDENCE. A graph  is plotted of the original function,
%	the training data, the trained network function, and the error bars.
%
%	See also
%	EVIDENCE, MLP, SCG, DEMARD, DEMMLP1
%

%	Copyright (c) Ian T Nabney (1996-2001)

clc;
disp('This demonstration illustrates the application of Bayesian')
disp('re-estimation to determine the hyperparameters in a simple regression')
disp('problem. It is based on a local quadratic approximation to a mode of')
disp('the posterior distribution and the evidence maximization framework of')
disp('MacKay.')
disp(' ')
disp('First, we generate a synthetic data set consisting of a single input')
disp('variable x sampled from a Gaussian distribution, and a target variable')
disp('t obtained by evaluating sin(2*pi*x) and adding Gaussian noise.')
disp(' ')
disp('Press any key to see a plot of the data together with the sine function.')
pause;

% Generate the matrix of inputs x and targets t.

ndata = 16;			% Number of data points.
noise = 0.1;			% Standard deviation of noise distribution.
randn('state', 0);
x = 0.25 + 0.07*randn(ndata, 1);
t = sin(2*pi*x) + noise*randn(size(x));

% Plot the data and the original sine function.
h = figure;
nplot = 200;
plotvals = linspace(0, 1, nplot)';
plot(x, t, 'ok')
xlabel('Input')
ylabel('Target')
hold on
axis([0 1 -1.5 1.5])
fplot('sin(2*pi*x)', [0 1], '-g')
legend('data', 'function');

disp(' ')
disp('Press any key to continue')
pause; clc;

disp('Next we create a two-layer MLP network having 3 hidden units and one')
disp('linear output. The model assumes Gaussian target noise governed by an')
disp('inverse variance hyperparmeter beta, and uses a simple Gaussian prior')
disp('distribution governed by an inverse variance hyperparameter alpha.')
disp(' ');
disp('The network weights and the hyperparameters are initialised and then')
disp('the weights are optimized with the scaled conjugate gradient')
disp('algorithm using the SCG function, with the hyperparameters kept')
disp('fixed. After a maximum of 500 iterations, the hyperparameters are')
disp('re-estimated using the EVIDENCE function. The process of optimizing')
disp('the weights with fixed hyperparameters and then re-estimating the')
disp('hyperparameters is repeated for a total of 3 cycles.')
disp(' ')
disp('Press any key to train the network and determine the hyperparameters.')
pause;

% Set up network parameters.
nin = 1;		% Number of inputs.
nhidden = 3;		% Number of hidden units.
nout = 1;		% Number of outputs.
alpha = 0.01;		% Initial prior hyperparameter.
beta_init = 50.0;	% Initial noise hyperparameter.

% Create and initialize network weight vector.
net = mlp(nin, nhidden, nout, 'linear', alpha, beta_init);

% Set up vector of options for the optimiser.
nouter = 3;			% Number of outer loops.
ninner = 1;			% Number of innter loops.
options = zeros(1,18);		% Default options vector.
options(1) = 1;			% This provides display of error values.
options(2) = 1.0e-7;		% Absolute precision for weights.
options(3) = 1.0e-7;		% Precision for objective function.
options(14) = 500;		% Number of training cycles in inner loop.

% Train using scaled conjugate gradients, re-estimating alpha and beta.
for k = 1:nouter
  net = netopt(net, options, x, t, 'scg');
  [net, gamma] = evidence(net, x, t, ninner);
  fprintf(1, '\nRe-estimation cycle %d:\n', k);
  fprintf(1, '  alpha =  %8.5f\n', net.alpha);
  fprintf(1, '  beta  =  %8.5f\n', net.beta);
  fprintf(1, '  gamma =  %8.5f\n\n', gamma);
  disp(' ')
  disp('Press any key to continue.')
  pause;
end

fprintf(1, 'true beta: %f\n', 1/(noise*noise));

disp(' ')
disp('Network training and hyperparameter re-estimation are now complete.')
disp('Compare the final value for the hyperparameter beta with the true')
disp('value.')
disp(' ')
disp('Notice that the final error value is close to the number of data')
disp(['points (', num2str(ndata),') divided by two.'])
disp(' ')
disp('Press any key to continue.')
pause; clc;
disp('We can now plot the function represented by the trained network. This')
disp('corresponds to the mean of the predictive distribution. We can also')
disp('plot ''error bars'' representing one standard deviation of the')
disp('predictive distribution around the mean.')
disp(' ')
disp('Press any key to add the network function and error bars to the plot.')
pause;

% Evaluate error bars.
[y, sig2] = netevfwd(mlppak(net), net, x, t, plotvals);
sig = sqrt(sig2);

% Plot the data, the original function, and the trained network function.
[y, z] = mlpfwd(net, plotvals);
figure(h); hold on;
plot(plotvals, y, '-r')
xlabel('Input')
ylabel('Target')
plot(plotvals, y + sig, '-b');
plot(plotvals, y - sig, '-b');
legend('data', 'function', 'network', 'error bars');

disp(' ')
disp('Notice how the confidence interval spanned by the ''error bars'' is')
disp('smaller in the region of input space where the data density is high,')
disp('and becomes larger in regions away from the data.')
disp(' ')
disp('Press any key to end.')
pause; clc; close(h);
%clear all

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
clear;
x = uiuc_sample;
options.J = 5;
options.M = 2;

Wop = wavelet_factory_3d_pyramid(options, options, options);

Sx = scat(x, Wop);
%%
Sx_renorm = renorm_scat_spatial(Sx);
ssx = mean(mean(format_scat(Sx),2),3);
ssx_rn = mean(mean(format_scat(Sx_renorm),2),3);
plot([ssx(2:end),ssx_rn(2:end)])

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
%Number of decision makers k, i.e.
k=4;
%Number of criterias to be evaluated n i.e.
n=4;
%Information about is the given criteria cost or benefit criteria:
%Benefit=1, Cost=2, should be given as a vector of length n. i.e.
criteria=[1 1 1 2];

%Number of attributes to be ranked (suppliers in example case). i.e.
m=6;

%Criteria for selecting Fuzzy Positive Ideal Solution FPIS and Fuzzy Negative Ideal Solution FNIS.
%Three possible option are now implemented. See [1.] for more information. i.e.

ideal=3;
%Chosen similarity measure: Integer number within {1,2,3,4}
simi=4;

lingvar

WD1={VH H H H};
WD2={H VH VH H};
WD3={MH H H MH};
WD4={M M MH MH};
WD={WD1;WD2;WD3;WD4};

FDM1={G G G F; MG G MG G; F F G VG; F P G G; MG MP MG MG; G MP F G};
FDM2={MG G MG G; G MG G G; F F G VG; MG MP MG MG; F MP F MG; MG P F VG};
FDM3={G MG MG G; F MG G F; MG P F G; MG MP MG G; F MP F G; MG P MG VG};
FDM4={G MG G G; MG G G MG; G F MG G; F P G G; MG MP MG MG; MG MP F G};
FDM={FDM1;FDM2;FDM3;FDM4};

[CCS,Sstar,Sneg,Order]=Stopsis(WD,FDM,k,m,n,ideal,criteria,simi);

disp('Similarities w.r.t. attribute and FPIS:')
Sstar
disp('Similarities w.r.t. attribute and FNIS:')
Sneg

disp('Order of the attributes:')
Order
disp('Closeness values with similarity for choosing attributes:')
CCS

author: Sable

repo: mcbench-benchmarks


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% Transformations of the time variable - discrete time signals

%upsampling-downsampling
 x=[1,2,3,4,5,6]
 a=2;
 xds=downsample(x,a)
 xds=x(1:a:end)

 a=1/2;
 xups=upsample(x,1/a)
 xups=zeros(1,1/a*length(x))
 xups(1:1/a:end)=x

 %all transformations
 n=-20:20;
 p=( (n>=-10)&(n<=10));
 x=(0.9.^n).*p;
 stem(n,x);
 legend('x[n]') ;

 figure
 stem(n+10,x)
 legend('x[n-10]')

 figure
 stem(n-10,x)
 legend('x[n+10]')

 figure
 a=2;
 xd=downsample(x,a)
 nd=-10:10;
 stem(nd,xd);
 legend('x[2n]')

 figure
 a=1/3
 xup=upsample(x,1/a)
 nup=-61:61
 stem(nup,xup)
 legend('x[(1/3) n]')

 figure
 stem(-n,x)
 legend('x[-n]')

 figure
 n=-2:4;
 n1=-fliplr(n)
 x=0.9.^n
 x1=fliplr(x)
 subplot(121);
 stem(n,x);
 title('x[n]=0.9^n');
 subplot(122);
 stem(n1,x1);
 title('x[-n]');


text:
function imgs = modcrop(imgs, modulo)

for i = 1:numel(imgs)
    sz = size(imgs{i});
    sz = sz - mod(sz, modulo);
    imgs{i} = imgs{i}(1:sz(1), 1:sz(2));
end

author: thomas-koehler

repo: SupER


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Goal: Fix a binary image by removing bad connectivities
% inObject: filename string of a binary object in .MAT format
% info = varargin{1}: connectivity
% info can be one of four values between 1 and 4 - '1=(6+,18), 2=(18,6+),
% 3=(6,26), 4=(26,6)'
% epsilon = varargin{2}: hole size, holes bigger than epsilon won't be
% filled
%
% 04/15/2002 - created by Li Shen
% Modified by Sungeun Kim (10-09-08)

function [bim, origin, vxsize, new_name] = fix_bad_topology(inObject, confs)

epsilon = confs.Epsilon;
switch deblank(char(confs.Connectivity))
    case '(6+,18)'
        conn = 1;
    case '(18,6+)'
        conn = 2;
    case '(6,26)'
        conn = 3;
    case '(26,6)'
        conn = 4;
end

% Read input object dataset
if ischar(inObject)
    [path,name,ext] = fileparts(inObject);

    if ~strcmp(ext, '.mat')
        disp('Input object data should be Matlab data file (*.mat)');
        return;
    end

    load(inObject);
    roi = bim;
else
    disp('Input Object should be the filename of a binary object (_bim.mat).');
    return;
end

tic;

infostr{1} = 'outlier or hole';
infostr{2} = 'vertex conn.';
infostr{3} = 'edge conn.';
infostr{4} = 'ring or hole';
infostr{5} = 'single ring';

% display information
disp('Make a simply connected volume by removing bad voxels');
disp(sprintf('(1) %s, (2) %s, (3) %s, (4) %s, (5) s', ...
              infostr{1}, infostr{2}, infostr{3}, infostr{4}, infostr{5}));

% Make binary image, in case reslice generates 2's
ind = find(roi>0); roi(ind) = 1;
vol = length(ind);
ind = find(roi<1); roi(ind) = 0;

% routine fix
roi = routine_fix(roi);
[vertices, faces] = gen_surf_data(roi,origin,vxsize);

% fill 3d holes
for k=1:3
    if (size(vertices,1)-size(faces,1)==2)
        disp('===> No 3D holes, skip');
        break;
    end
    disp('===> fixing 3D holes ...');
    roi = fill3dholes(roi,conn,epsilon,name);
    disp('===> routine fix again ...');
    roi = routine_fix(roi);
    [vertices, faces] = gen_surf_data(roi,origin,vxsize);
end

fvn = length(find(roi~=bim));
bim = roi;

disp(sprintf('fix_vox_num (%d) / vol (%d) = %f%%, save ...', ...
                      fvn, vol, fvn*100/vol));

% Save surface object to a new file
postfix = name(end-2:end);
if strcmp(postfix, 'bim') | strcmp(postfix,'fix')
    new_name = [confs.OutDirectory '/' name(1:end-3) 'fix'];
else
    new_name = [confs.OutDirectory '/' name '_fix'];
end

if exist(new_name,'file')
    prompt = {'Enter new filename:'};
    dlg_title = 'New File Name';
    num_lines = 1;
    def = {new_name};
    answer = inputdlg(prompt,dlg_title,num_lines,def);
    new_name = answer{1};
end

save(new_name, 'bim', 'origin', 'vxsize');

seconds = toc;
hours = seconds/3600;
disp(['seconds=', num2str(seconds), ' hours=', num2str(hours)]);

return;

%
% routine fix
%

function roi = routine_fix(roi)

DIM = size(roi);

fixset = [];

for j = 1:10
    [roi, fs1] = outl_hole(roi, DIM,1);
    [roi, fs2] = vertex_conn(roi, DIM,2);
    [roi, fs3] = edge_conn(roi, DIM,3);
    if (isempty(fs1) & isempty(fs2) & isempty(fs3))
        break;
    else
        if (~isempty(fs1))
            fixset(end+1:end+size(fs1,1),:) = fs1;
        end
        if (~isempty(fs2))
            fixset(end+1:end+size(fs2,1),:) = fs2;
        end
        if (~isempty(fs3))
            fixset(end+1:end+size(fs3,1),:) = fs3;
        end
    end
end

[roi, fs4] = ring_hole(roi, DIM, 4);
if (~isempty(fs4))
    fixset(end+1:end+size(fs4,1),:) = fs4;
    for j = 1:10
        [roi, fs1] = outl_hole(roi, DIM,1);
        [roi, fs2] = vertex_conn(roi, DIM,2);
        [roi, fs3] = edge_conn(roi, DIM,3);
        [roi, fs4] = ring_hole(roi, DIM,4);
        if (isempty(fs1) & isempty(fs2) & isempty(fs3) & isempty(fs4))
            break;
        else
            if (~isempty(fs1))
                fixset(end+1:end+size(fs1,1),:) = fs1;
            end
            if (~isempty(fs2))
                fixset(end+1:end+size(fs2,1),:) = fs2;
            end
            if (~isempty(fs3))
                fixset(end+1:end+size(fs3,1),:) = fs3;
            end
            if (~isempty(fs4))
                fixset(end+1:end+size(fs4,1),:) = fs4;
            end
        end
    end
end

fix_vox_num = size(fixset,1);

return;

%
% remove disconnected small components and holes
%

function [roi, fixset] = outl_hole(roi, d, infoid)

roi = reshape(roi,d);

fixset = [];

% check seperated small components (k=1)
% and then check holes
val = [0 1];
for k = 1:2
	[L, num] = bwlabeln(roi,6);

	if num>1
		% find the background
		for i = 1:num
            cnt(i) = length(find(L == i));
		end
		[bkgd_size, bkgd_ind] = max(cnt);
        sigind = find(L~=0 & L~=bkgd_ind);
        roi(sigind) = 0;
        [xs,ys,zs] = ind2sub(d, sigind);
        fixset(end+1:end+length(xs),:) = [[xs,ys,zs], xs*0+val(k), xs*0+infoid];
	end

	roi = 1-roi;
end

return;

%
% fix bad vertex connectivities
%

function [roi, fixset] = vertex_conn(roi, d, infoid)

roi = reshape(roi,d);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% work on 0 vertex connectivies
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% make a work area so that all border voxels belong to the background
w = ones(d+2);
w(2:d(1)+1,2:d(2)+1,2:d(3)+1) = roi;
wd = d+2;
ind = find(w == 0);
w = zeros(d+2);
w(2:d(1)+1,2:d(2)+1,2:d(3)+1) = roi;

[xs,ys,zs] = ind2sub(wd,ind);

h = sum([w(sub2ind(wd,xs-1,ys,zs)), ...
         w(sub2ind(wd,xs+1,ys,zs)), ...
         w(sub2ind(wd,xs,ys-1,zs)), ...
         w(sub2ind(wd,xs,ys+1,zs)), ...
         w(sub2ind(wd,xs,ys,zs-1)), ...
         w(sub2ind(wd,xs,ys,zs+1))]')';

fixset = [];

% fix ring 0 voxels (single ring)
ind1 = ind(find(h==4 & ...
                w(sub2ind(wd,xs-1,ys,zs))==0 & ...
                w(sub2ind(wd,xs+1,ys,zs))==0 ...
           ));
ind2 = ind(find(h==4 & ...
                w(sub2ind(wd,xs,ys-1,zs))==0 & ...
                w(sub2ind(wd,xs,ys+1,zs))==0 ...
           ));
ind3 = ind(find(h==4 & ...
                w(sub2ind(wd,xs,ys,zs-1))==0 & ...
                w(sub2ind(wd,xs,ys,zs+1))==0 ...
           ));
ind1 = [ind1;ind2;ind3];
n4 = length(ind1);
if (n4~=0)
	w(ind1) = 1;
    [a, b, c] = ind2sub(wd,ind1);
	fixset(end+1:end+n4,:) = [[a, b, c]-1, ind1*0+1, ind1*0+5];
end

% right up forward, right up backward, right down forward, right down backward
ruf = find(w(sub2ind(wd,xs-1,ys+1,zs+1)) == 0 & ...
           sum([w(sub2ind(wd,xs-1,ys,zs)), ...
                w(sub2ind(wd,xs,ys+1,zs)), ...
                w(sub2ind(wd,xs-1,ys+1,zs)), ...
                w(sub2ind(wd,xs,ys,zs+1)), ...
                w(sub2ind(wd,xs-1,ys,zs+1)), ...
                w(sub2ind(wd,xs,ys+1,zs+1))]')' == 6);

rub = find(w(sub2ind(wd,xs-1,ys+1,zs-1)) == 0 & ...
           sum([w(sub2ind(wd,xs-1,ys,zs)), ...
                w(sub2ind(wd,xs,ys+1,zs)), ...
                w(sub2ind(wd,xs-1,ys+1,zs)), ...
                w(sub2ind(wd,xs,ys,zs-1)), ...
                w(sub2ind(wd,xs-1,ys,zs-1)), ...
                w(sub2ind(wd,xs,ys+1,zs-1))]')' == 6);

rdf = find(w(sub2ind(wd,xs+1,ys+1,zs+1)) == 0 & ...
           sum([w(sub2ind(wd,xs+1,ys,zs)), ...
                w(sub2ind(wd,xs,ys+1,zs)), ...
                w(sub2ind(wd,xs+1,ys+1,zs)), ...
                w(sub2ind(wd,xs,ys,zs+1)), ...
                w(sub2ind(wd,xs+1,ys,zs+1)), ...
                w(sub2ind(wd,xs,ys+1,zs+1))]')' == 6);

rdb = find(w(sub2ind(wd,xs+1,ys+1,zs-1)) == 0 & ...
           sum([w(sub2ind(wd,xs+1,ys,zs)), ...
                w(sub2ind(wd,xs,ys+1,zs)), ...
                w(sub2ind(wd,xs+1,ys+1,zs)), ...
                w(sub2ind(wd,xs,ys,zs-1)), ...
                w(sub2ind(wd,xs+1,ys,zs-1)), ...
                w(sub2ind(wd,xs,ys+1,zs-1))]')' == 6);

cube = [];

if (~isempty(ruf))
    cube(end+1:end+length(ruf),:) = ...
        [xs(ruf),ys(ruf),zs(ruf), ...
         xs(ruf)-1,ys(ruf)+1,zs(ruf)+1];
end

if (~isempty(rub))
    cube(end+1:end+length(rub),:) = ...
        [xs(rub),ys(rub),zs(rub), ...
         xs(rub)-1,ys(rub)+1,zs(rub)-1];
end

if (~isempty(rdf))
    cube(end+1:end+length(rdf),:) = ...
        [xs(rdf),ys(rdf),zs(rdf), ...
         xs(rdf)+1,ys(rdf)+1,zs(rdf)+1];
end

if (~isempty(rdb))
    cube(end+1:end+length(rdb),:) = ...
        [xs(rdb),ys(rdb),zs(rdb), ...
         xs(rdb)+1,ys(rdb)+1,zs(rdb)-1];
end

% adjust according to cube (3: fix bad 0 vertex connectivities)
if (~isempty(cube))
    len = size(cube,1);
    for i = 1:len
        k = cube(i,:);
        if (w(k(1),k(2),k(3))==0 & w(k(4),k(5),k(6))==0)
            if (nb_sum(w,k(1),k(2),k(3)) >= nb_sum(w,k(4),k(5),k(6)))
                w(k(1),k(2),k(3)) = 1;
                fixset(end+1,:) = [[k(1), k(2), k(3)]-1, 1, infoid];
            else
                w(k(4),k(5),k(6)) = 1;
                fixset(end+1,:) = [[k(4), k(5), k(6)]-1, 1, infoid];
            end
        end
    end
end

cnts = [length(ruf), length(rub), length(rdf), length(rdb)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% work on 1 vertex connectivies
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ind = find(w);
[xs,ys,zs] = ind2sub(wd,ind);

% right up forward, right up backward, right down forward, right down backward
ruf = find(w(sub2ind(wd,xs-1,ys+1,zs+1)) == 1 & ...
           sum([w(sub2ind(wd,xs-1,ys,zs)), ...
                w(sub2ind(wd,xs,ys+1,zs)), ...
                w(sub2ind(wd,xs-1,ys+1,zs)), ...
                w(sub2ind(wd,xs,ys,zs+1)), ...
                w(sub2ind(wd,xs-1,ys,zs+1)), ...
                w(sub2ind(wd,xs,ys+1,zs+1))]')' == 0);

rub = find(w(sub2ind(wd,xs-1,ys+1,zs-1)) == 1 & ...
           sum([w(sub2ind(wd,xs-1,ys,zs)), ...
                w(sub2ind(wd,xs,ys+1,zs)), ...
                w(sub2ind(wd,xs-1,ys+1,zs)), ...
                w(sub2ind(wd,xs,ys,zs-1)), ...
                w(sub2ind(wd,xs-1,ys,zs-1)), ...
                w(sub2ind(wd,xs,ys+1,zs-1))]')' == 0);

rdf = find(w(sub2ind(wd,xs+1,ys+1,zs+1)) == 1 & ...
           sum([w(sub2ind(wd,xs+1,ys,zs)), ...
                w(sub2ind(wd,xs,ys+1,zs)), ...
                w(sub2ind(wd,xs+1,ys+1,zs)), ...
                w(sub2ind(wd,xs,ys,zs+1)), ...
                w(sub2ind(wd,xs+1,ys,zs+1)), ...
                w(sub2ind(wd,xs,ys+1,zs+1))]')' == 0);

rdb = find(w(sub2ind(wd,xs+1,ys+1,zs-1)) == 1 & ...
           sum([w(sub2ind(wd,xs+1,ys,zs)), ...
                w(sub2ind(wd,xs,ys+1,zs)), ...
                w(sub2ind(wd,xs+1,ys+1,zs)), ...
                w(sub2ind(wd,xs,ys,zs-1)), ...
                w(sub2ind(wd,xs+1,ys,zs-1)), ...
                w(sub2ind(wd,xs,ys+1,zs-1))]')' == 0);

cube = [];

if (~isempty(ruf))
    cube(end+1:end+length(ruf),:) = ...
        [xs(ruf),ys(ruf),zs(ruf), ...
         xs(ruf)-1,ys(ruf)+1,zs(ruf)+1];
end

if (~isempty(rub))
    cube(end+1:end+length(rub),:) = ...
        [xs(rub),ys(rub),zs(rub), ...
         xs(rub)-1,ys(rub)+1,zs(rub)-1];
end

if (~isempty(rdf))
    cube(end+1:end+length(rdf),:) = ...
        [xs(rdf),ys(rdf),zs(rdf), ...
         xs(rdf)+1,ys(rdf)+1,zs(rdf)+1];
end

if (~isempty(rdb))
    cube(end+1:end+length(rdb),:) = ...
        [xs(rdb),ys(rdb),zs(rdb), ...
         xs(rdb)+1,ys(rdb)+1,zs(rdb)-1];
end

% adjust according to cube (4: remove bad 1 voxel connectivities)
if (~isempty(cube))
    len = size(cube,1);
    for i = 1:len
        k = cube(i,:);
        if (w(k(1),k(2),k(3))==1 & w(k(4),k(5),k(6))==1)
            if (nb_sum(w,k(1),k(2),k(3)) <= nb_sum(w,k(4),k(5),k(6)))
                w(k(1),k(2),k(3)) = 0;
                fixset(end+1,:) = [[k(1), k(2), k(3)]-1, 0, infoid];
            else
                w(k(4),k(5),k(6)) = 0;
                fixset(end+1,:) = [[k(4), k(5), k(6)]-1, 0, infoid];
            end
        end
    end
end

cnts = [length(ruf), length(rub), length(rdf), length(rdb)];

% make roi consistent with w
roi = w(2:d(1)+1,2:d(2)+1,2:d(3)+1);

return;

%
% fix bad edge connectivities
%

function [roi, fixset] = edge_conn(roi, d, infoid)

roi = reshape(roi,d);

% make a work area so that all border voxels belong to the background
w = zeros(d+2);
w(2:d(1)+1,2:d(2)+1,2:d(3)+1) = roi;
wd = d+2;

ind = find(w);
[xs,ys,zs] = ind2sub(wd,ind);

squa = [];

% right up and right down
ru = find(w(sub2ind(wd,xs-1,ys,zs))   == 0 & ...
          w(sub2ind(wd,xs,ys+1,zs))   == 0 & ...
          w(sub2ind(wd,xs-1,ys+1,zs)) == 1);

if (~isempty(ru))
    squa(end+1:end+length(ru),:) = ...
                [xs(ru),ys(ru),zs(ru),...
                 xs(ru)-1,ys(ru)+1,zs(ru),...
                 xs(ru)-1,ys(ru),zs(ru),...
                 xs(ru),ys(ru)+1,zs(ru)];
end

rd = find(w(sub2ind(wd,xs+1,ys,zs))   == 0 & ...
          w(sub2ind(wd,xs,ys+1,zs))   == 0 & ...
          w(sub2ind(wd,xs+1,ys+1,zs)) == 1);

if (~isempty(rd))
    squa(end+1:end+length(rd),:) = ...
                [xs(rd),ys(rd),zs(rd),...
                 xs(rd)+1,ys(rd)+1,zs(rd),...
                 xs(rd)+1,ys(rd),zs(rd),...
                 xs(rd),ys(rd)+1,zs(rd)];
end

% forward up and forward down
fu = find(w(sub2ind(wd,xs-1,ys,zs))   == 0 & ...
          w(sub2ind(wd,xs,ys,zs+1))   == 0 & ...
          w(sub2ind(wd,xs-1,ys,zs+1)) == 1);

if (~isempty(fu))
    squa(end+1:end+length(fu),:) = ...
                [xs(fu),ys(fu),zs(fu),...
                 xs(fu)-1,ys(fu),zs(fu)+1,...
                 xs(fu)-1,ys(fu),zs(fu),...
                 xs(fu),ys(fu),zs(fu)+1];
end

fd = find(w(sub2ind(wd,xs+1,ys,zs))   == 0 & ...
          w(sub2ind(wd,xs,ys,zs+1))   == 0 & ...
          w(sub2ind(wd,xs+1,ys,zs+1)) == 1);

if (~isempty(fd))
    squa(end+1:end+length(fd),:) = ...
                [xs(fd),ys(fd),zs(fd),...
                 xs(fd)+1,ys(fd),zs(fd)+1,...
                 xs(fd)+1,ys(fd),zs(fd),...
                 xs(fd),ys(fd),zs(fd)+1];
end

% left forward and right forward
lf = find(w(sub2ind(wd,xs,ys-1,zs))   == 0 & ...
          w(sub2ind(wd,xs,ys,zs+1))   == 0 & ...
          w(sub2ind(wd,xs,ys-1,zs+1)) == 1);

if (~isempty(lf))
    squa(end+1:end+length(lf),:) = ...
                [xs(lf),ys(lf),zs(lf),...
                 xs(lf),ys(lf)-1,zs(lf)+1,...
                 xs(lf),ys(lf)-1,zs(lf),...
                 xs(lf),ys(lf),zs(lf)+1];
end

rf = find(w(sub2ind(wd,xs,ys+1,zs))   == 0 & ...
          w(sub2ind(wd,xs,ys,zs+1))   == 0 & ...
          w(sub2ind(wd,xs,ys+1,zs+1)) == 1);

if (~isempty(rf))
    squa(end+1:end+length(rf),:) = ...
                [xs(rf),ys(rf),zs(rf),...
                 xs(rf),ys(rf)+1,zs(rf)+1,...
                 xs(rf),ys(rf)+1,zs(rf),...
                 xs(rf),ys(rf),zs(rf)+1];
end

% adjust according to squa (5: remove bad edge connectivities)
fixset = [];
if (~isempty(squa))
    len = size(squa,1);
    for i = 1:len
        k = squa(i,:);
        if (sum([w(k(1),k(2),k(3)) ...
                 w(k(4),k(5),k(6)) ...
                 w(k(7),k(8),k(9)) ...
                 w(k(10),k(11),k(12)) ...
                 ]==[1 1 0 0])==4)
            [start, val] = decide_squa(w,k);
            w(k(start),k(start+1),k(start+2)) = val;
            fixset(end+1,:) = [[k(start), k(start+1), k(start+2)]-1, val, infoid];
        end
    end
end

cnts = [length(ru),length(rd),length(fu),length(fd),length(lf),length(rf)];

% make roi consistent with w
roi = w(2:d(1)+1,2:d(2)+1,2:d(3)+1);

return;

%
% check ring on one slice, and hole between two slices
%

function [roi, fixset] = ring_hole(roi, d, infoid)

roi = reshape(roi,d);

% make a work area so that all border voxels belong to the background
w = zeros(d+2);
w(2:d(1)+1,2:d(2)+1,2:d(3)+1) = roi;
wd = d+2;

% exchange background and object
w = 1-w;

fixset = [];

for k = 1:3
    % label each slice
	for i = 1:wd(k)
        switch k
        case 1
            im = reshape(w(i,:,:),wd(2),wd(3));
        case 2
            im = reshape(w(:,i,:),wd(1),wd(3));
        case 3
            im = reshape(w(:,:,i),wd(1),wd(2));
        end
        [lab, num(i)] = bwlabel(im,4);
        L{i} = lab;
	end
    % original roi background should be labeled as '1'
    ind = find(num>1);
    for i = 1:length(ind);
        lab = L{ind(i)};
        for j = 1:num(ind(i))
            cnt(j) = length(find(lab == j));
        end
        [mval, mind] = max(cnt);
        if (mind > 1)
            disp('Make background labeled as 1');
            ind1 = find(lab == 1);
            ind2 = find(lab == mind);
            lab(ind2) = 1; % background
            lab(ind1) = mind;
            L{ind(i)} = lab;
        end
    end
    % check multiple component slices
    for i = 1:length(ind);
        lab = L{ind(i)};
        for j = 2:num(ind(i))
            % connect to background in the previous slice
            con_prev = find(lab == j & L{ind(i)-1} == 1);
            % connect to background in the next slice
            con_next = find(lab == j & L{ind(i)+1} == 1);
            if (~isempty(con_prev) & ~isempty(con_next))
                [xs, ys] = find(lab == j);
                switch k
                case 1
                    im = reshape(w(ind(i),:,:),wd(2),wd(3));
                    sigind = sub2ind(wd,xs*0+ind(i),xs,ys);
                case 2
                    im = reshape(w(:,ind(i),:),wd(1),wd(3));
                    sigind = sub2ind(wd,xs,xs*0+ind(i),ys);
                case 3
                    im = reshape(w(:,:,ind(i)),wd(1),wd(2));
                    sigind = sub2ind(wd,xs,ys,xs*0+ind(i));
                end
                w(sigind) = 0; % note that really it is assigned 1 to roi
                [xs,ys,zs] = ind2sub(wd,sigind);
                fixset(end+1:end+length(xs),:) = [[xs,ys,zs]-1, xs*0+1, xs*0+infoid];
            end
        end
    end
    % check hole between slices
    for i = 2:(wd(k)-2)
        overlap = L{i}*0;
        overlap_bkgd = find(L{i} == 1 & L{i+1} == 1);
        overlap(overlap_bkgd) = 1;
        [lab, m] = bwlabel(overlap,4);
        if (m>1)
            % deal with only the minimum one;
            for j = 1:m
                cnt(j) = length(find(lab == j));
            end
            [mval, mind] = min(cnt);
            % connect to background in the previous slice
            con_prev = find(lab == mind & L{i-1} == 1);
            % connect to background in the current slice
            con_curr = find(lab == mind & L{i} == 1);
            % connect to background in the next slice
            con_next = find(lab == mind & L{i+1} == 1);
            % connect to background in the next next slice
            con_next_next = find(lab == mind & L{i+2} == 1);
            if ((~isempty(con_curr) | ~isempty(con_prev)) & ...
                (~isempty(con_next) | ~isempty(con_next_next)))
                [xs, ys] = find(lab == mind);
                switch k
                case 1
                    sigind = [sub2ind(wd,xs*0+i,xs,ys);sub2ind(wd,xs*0+i+1,xs,ys)];
                case 2
                    sigind = [sub2ind(wd,xs,xs*0+i,ys);sub2ind(wd,xs,xs*0+i+1,ys)];
                case 3
                    sigind = [sub2ind(wd,xs,ys,xs*0+i);sub2ind(wd,xs,ys,xs*0+i+1)];
                end
                w(sigind) = 0; % note that really it is assigned 1 to roi
                [xs,ys,zs] = ind2sub(wd,sigind);
                fixset(end+1:end+length(xs),:) = [[xs,ys,zs]-1, xs*0+1, xs*0+infoid];
            end
        end
    end
    L = [];
    num = [];
end

% make roi consistent with w
roi = 1 - w(2:d(1)+1,2:d(2)+1,2:d(3)+1);

return;

%
% select a point to change its value
% according to maximum number of different neighbours
%

function [start, val] = decide_squa(w,k)

for i = 1:3:10
    cnt((i-1)/3+1) = nb_sum(w, k(i), k(i+1), k(i+2));
end

cnt(1:2) = 6 - cnt(1:2);

[v, in] = max(cnt);

start = (in-1)*3+1;
if (in < 3)
    val = 0;
else
    val = 1;
end

return;

%
% calculate the sum of the direct neighbour
%
function b = nb_sum(w, x, y, z)

b = sum([w(x-1,y,z), ...
         w(x+1,y,z), ...
         w(x,y-1,z), ...
         w(x,y+1,z), ...
         w(x,y,z-1), ...
         w(x,y,z+1)]);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%% Test Problems for DSDP
clc
clear

%% CSDP Options Function
clc
csdpset

a = csdpset

%% LP1 [-31.4]
clc
f = -[6 5]';
A = sparse([1,4; 6,4; 2, -5]);
b = [16;28;6];
lb = [0;0];
ub = [10;10];

opts = [];
opts.display = 2;
opts.maxiter = 15;
% opts.writeprob = 'prob.dat-s';
% opts.writesol = 'sol.dat-s';
[x,ff,e,i,X] = csdp(f,A,b,lb,ub,[],[],opts)

%% LP2 [2]
clc
f = [8,1]';
A = sparse([-1,-2;1,-4;3,-1;1,5;-1,1;-1,0;0,-1]);
b = [-4,2,21,39,3,0,0]';

[x,p,d,e] = csdp(f,A,b)

%% LP3 [-3.75]
clc
f = -[-1, 2]';
A = sparse([2, 1;-4, 4]);
b = [5, 5]';

[x,p,d,e] = csdp(f,A,b)

%% LP4 [-97.5]
clc
clear all
f = -[1 2 3]';
A = sparse([-1,1,1; 1,-3,1]);
b = [20,30]';
Aeq = sparse([1 1 1]);
A = [A;Aeq;-Aeq];
beq = 40;
b = [b;beq;-beq];
lb = [0 0 0]';
ub =[40 inf inf]';

[x,fval,e,i] = csdp(f,A,b,lb,ub)

% [x,ff] = clp(f,[A;Aeq],[-Inf(size(b));beq],[b;beq],lb,ub)

%% SDP1 [1.4142]
clc
clear
%Objective
f = 1;
%SDP [x sqrt(2); sqrt(2) x] >= 0
A = eye(2);
C = -[0 sqrt(2); sqrt(2) 0];
sdp = sparse([C(:) A(:)]);
%Options
opts.display=2;
[x,ff,e,i] = csdp(f,[],[],[],[],sdp,[],opts)

%% SDP 1b Johan [-1]
clc
clear
%Objective
f = 1;
%SDP [1 x; x 1] >= 0
A = [0 1; 1 0];
C = -eye(2);
sdp = sparse([C(:) A(:)]);
%Options
opts.display = 2;
tic
[x,ff,e,i] = csdp(f,[],[],[],[],sdp,[],opts)
toc

%% SDP2 [4]
clc
clear
%Objective
f = [1;1];
%Linear Constraints
lb = [0;0];
ub = [10;10];
%SDP Constraints [x1 2; 2 x2] >= 0
C = -[0 2; 2 0];
A0 = [1 0; 0 0];
A1 = [0 0; 0 1];
sdp = sparse([C(:) A0(:) A1(:)]);
%Setup Options
opts.display=2;
[x,ff,e,i] = csdp(f,[],[],lb,ub,sdp,[],opts)

%% SDP3 [1.2] (Note assume matrices are symmetric triu to avoid x5)
clc
clear
%Objective
f = [1 0 0 0]';
%SDP Constraint1 [x2 x3;x3 x4] <= x1*eye(2)
C = zeros(2);
A0 = eye(2);
A1 = -[1 0; 0 0];
A2 = -[0 1; 1 0];
A3 = -[0 0; 0 1];
sdp{1} = sparse([C(:) A0(:) A1(:) A2(:) A3(:)]);
%SDP Constraint2 [x2 x3; x3 x4] >= [1 0.2; 0.2 1]
C = [1 0.2; 0.2 1];
A0 = zeros(2);
A1 = [1 0; 0 0];
A2 = [0 1; 1 0];
A3 = [0 0; 0 1];
sdp{2} = sparse([C(:) A0(:) A1(:) A2(:) A3(:)]);
%Setup Options
opts.display=2;
[x,ff,e,i] = csdp(f,[],[],[],[],sdp,[],opts)

%% SDP4 [10.1787]
clc
clear
%Objective
f = [1 1 1];
%Linear Constraints
lb = [10;0;0];
ub = [1000;1000;1000];
%SDP Constraints [x1 1 2; 1 x2 3; 2 3 100] >= 0
C = -[0 1 2; 1 0 3; 2 3 100];
A0 = [1 0 0; 0 0 0; 0 0 0];
A1 = [0 0 0; 0 1 0; 0 0 0];
A2 = zeros(3);
sdp = sparse([C(:) A0(:) A1(:) A2(:)]);
%Setup Options
opts.display=2;
[x,ff,e,i] = csdp(f,[],[],lb,ub,sdp,[],opts)

%% SDP CSDP Example
clc
clear
%Objective
f = -[1 2];
%SDP Constraint 1
C = [2 1; 1 2];
A0 = -[3 1; 1 3];
A1 = -zeros(2);
sdp{1} = sparse([C(:) A0(:) A1(:)]);
%SDP Constraint 2
C = [3 0 1; 0 2 0; 1 0 3];
A0 = -zeros(3);
A1 = -[3 0 1; 0 4 0; 1 0 5];
sdp{2} = sparse([C(:) A0(:) A1(:)]);
%SDP Constraint 3
C = zeros(2);
A0 = -[1 0; 0 0];
A1 = -[0 0; 0 1];
sdp{3} = sparse([C(:) A0(:) A1(:)]);
%Setup Options
opts.display=2;
% opts.writesol='sol.dat-s';
[x,ff,e,i,X] = csdp(f,[],[],[],[],sdp,[],opts)

X{1}
X{2}
X{3}

%% SDP5
% clc
% clear
% load dsdpdebug
%
% opts.display=2;
% [y,fvals,exitflag,stats,X] = csdp(model.f,[],[],[],[],model.sdcone,[],opts)

%% SDP 6 (testing memory leaks - use task manager to view matlab ram)
% clc
% clear all
% %Number of cones
% n = 10000;
% %Objective
% f = [1;1];
% %Linear Constraints
% lb = [0;0];
% ub = [10;10];
% %SDP Constraints [x1 2; 2 x2] >= 0
% ind = triu(ones(2))==1;
% C = [0 2; 2 0];
% A0 = [1 0; 0 0];
% A1 = [0 0; 0 1];
% sdp = repmat({sparse([C(:) A0(:) A1(:)])},n,1);
% %Setup Options
% opts.display=2;
% opts.maxtime = 10;
% [x,ff,e,i,X] = csdp(f,[],[],lb,ub,sdp,[],opts);

%% YALMIP comparisons
% %%
% x = sdpvar(1,1);
% %y = sdpvar(1,1);
% z = sdpvar(1,1);
%
% X = [x 2;2 z];
%
% F = set('X>=0');
% % F = F+set('x>=1');
% F = F+set('z>=1');
% F = F+set('x<=10');
% F = F+set('z<=10');
% sol = solvesdp(F,x+z)
%
% %%
% [At,b1,c1,K]=readsdpa('complete2uub.dat-s')
%
% A1 = full(At)
%
% %%
% [At,b2,c2,K]=readsdpa('prob.dat-s')
%
% A2 = full(At)
%
%
% %%
% t = sdpvar(1,1);
% Y = sdpvar(2,2);
% F = set('Y<=t*eye(2)');
% F = F+set('Y>=[1 0.2;0.2 1]');
% sol = solvesdp(F,t)
%
% %%
% [At,b1,c1,K]=readsdpa('test.dat-s')
%
% A1 = full(At)
%
% %%
% [At,b2,c2,K]=readsdpa('prob.dat-s')
%
% A2 = full(At)
%
%
% %%
% clc
% x = sdpvar(1,1);
% y = sdpvar(1,1);
%
% con = [x + 4*y <= 16; 6*x + 4*y <= 28, 2*x - 5*y <= 6; 0 <= x <= 10; 0 <= y <= 10];
% obj = -6*x - 5*y;
%
% opts = sdpsettings('solver','csdp');
% sol = solvesdp(con,obj,opts)
%
% double(x)
% double(y)
% double(obj)
%
% %%
% [At,b1,c1,K]=readsdpa('lp1.dat-s')
%
% A1 = full(At)
%
% %%
% [At,b2,c2,K]=readsdpa('prob.dat-s')
%
% A2 = full(At)


text:
 function xs = tpl_inc(x, Gb, yi, bi, ri, R, niter, pixmax, curv, ...
		subout, enhance, gi, denom, relax0, chat)
%function xs = tpl_inc(x, Gb, yi, bi, ri, R, niter, pixmax, curv, ...
%		subout, enhance, gi, denom, relax0, chat)
%
% TRIOT: incremental SPS algorithm for transmission Poisson problem
% (ordered subsets separable paraboloidal surrogates)
% in
%	x	[np,1]		initial estimate
%	Gb			Gblock object (see tpl_os_sps_test.m)
%	yi,bi,ri [nb,na]	see tr_fbp.m (for model too)
%	R		penalty structure (see Robject.m)
%				(or empty for ML case)
%	niter		# iterations
%	pixmax		upper constraint for pixel values or [lower, upper]
%			can be scalar (e.g. 'inf') or an array the size of x
%	curv		'oc' for erdogan's optimal curvatures
%			'mc' for maximum curvature
%			'pc' for erdogan's fast precomputed curvatures,
%				which usually provides faster convergence,
%				but can be nonmonotone.
%	subout	[1]	nonzero value: return all subiterates
%			otherwise, return only iterates
%	enhance [1]	nonzero value: Hsiao's enhancement
%			otherwise, no enhancement
%	gi	[nd,1]	sum(G')'
%	denom	[np,1]		(if precomputed denominator)
%	relax0	[1] or [2]	relax0 or (relax0, relax_rate)
% out
%	xs [np,niter]	updated image vectors each iteration
%
% Copyright Mar 2000, Jeff Fessler, The University of Michigan
% 2002-2-20 update for Gblock object
% Modified for incremental version by Sangtae Ahn, from ~sangtaea/trans.

epsilon = 1e-10;

if nargin < 3, ir_usage, end

Gb = block_ob(Gb, 'ensure'); % make it a block object (if not already)
nblock = block_ob(Gb, 'n');
starts = subset_start(nblock);

if ~isvar('bi') || isempty(bi)
	bi = ones(size(yi));
end
if ~isvar('ri') || isempty(ri)
	ri = zeros(size(yi));
end

if ~isvar('R'), R = []; end
if ~isvar('niter')	|| isempty(niter),	niter = 2;	end
if ~isvar('pixmax')	|| isempty(pixmax),	pixmax = inf;	end
if length(pixmax) == 2
	pixmin = pixmax(1);
	pixmax = pixmax(2);
else
	pixmin = 0;
end
if ~isvar('curv')	|| isempty(curv),	curv = 'oc';	end
if ~isvar('subout')	|| isempty(subout),	subout = 0;	end
if ~isvar('enhance')	|| isempty(enhance),	enhance = 0;	end
if ~isvar('chat')	|| isempty(chat),	chat = false;	end
if ~isvar('relax0')	|| isempty(relax0),	relax0 = 1;	end

if streq(curv, 'mc')
	subscale = @(iset) iset;
else
	subscale = @(iset) 1;
end

if length(relax0) == 1
	relax_rate = 0;
elseif length(relax0) == 2
	relax_rate = relax0(2);
	relax0 = relax0(1);
else
	error relax
end

if length(niter) == 2
	osspsiter = niter(1);
	niter = niter(2);
else
	osspsiter = 1;
end

trl_check(yi, bi, ri);

if ~isvar('gi') || isempty(gi)
	gi = reshape(sum(Gb'), size(yi));	% g_i = sum_j g_ij
end

starts = subset_start(nblock);
[nb na] = size(yi);

% precompute denominator if needed
if (~isvar('denom') || isempty(denom)) && (streq(curv, 'pc') || streq(curv, 'mc'))
	ni = my_trl_curvature(yi, bi, ri, [], curv); % precomputed curvatures
	%
	% separate denominators for each subset
	%
	if streq(curv, 'mc')
		denom = zeros(numel(x), nblock);
		for iset=1:nblock
			iblock = starts(iset);
			ia = iblock:nblock:na;
			denom(:,iset) = Gb{iblock}' * col(gi(:,ia) .* ni(:,ia));
		end
	else
	%
	% one denominator shared by all subsets
	%
		denom = Gb' * col(gi .* ni) / nblock;
	end
end

% "precomputed curvature"
if ~streq(curv, 'pc')
	ni = my_trl_curvature(yi, bi, ri, [], 'pc');
	denom_pc = Gb' * col(gi .* ni) / nblock;
else
	denom_pc = denom;
end

if subout
	xs = zeros(length(x), (niter-1)*nblock+1);
else
	xs = zeros(length(x), niter);
end

x = max(x,pixmin);
x = min(x,pixmax);
xs(:,1) = x;
if subout, idx = 2; end
num_store = zeros(length(x), nblock);
den_store = zeros(length(x), nblock);

%
% initialization by os-sps-pc
%
for iter = 1:osspsiter
	for iset=1:nblock
		iblock = starts(iset);
		ia = iblock:nblock:na;

		li = Gb{iblock} * x;	% l=G*x "line integrals"
		li = reshape(li, nb, length(ia));
		bel = bi(:,ia) .* exp(-li);
		yb = bel + ri(:,ia);	% predicted meas. means
		dothi = (1 - yi(:,ia) ./ yb) .* bel;

		if streq(curv, 'oc')
			ni = my_trl_curvature(yi(:,ia), bi(:,ia), ri(:,ia), ...
				li, 'oc');
			denom = Gb{iblock}' * col(gi(:,ia) .* ni);
		end

		if isempty(R)
			grad = Gb{iblock}' * dothi(:);
			den_osps = max(denom_pc, epsilon);
			den_store(:, iset) = max(denom(:,subscale(iset)), ...
						epsilon);
		else
			grad = Gb{iblock}' * dothi(:) - R.cgrad(R, x)/nblock;
			den_osps = max(denom_pc + R.denom(R, x)/nblock, ...
					epsilon);
			den_store(:, iset) = max(denom(:,subscale(iset)) + ...
					R.denom(R, x)/nblock, epsilon);
		end

		if iter == osspsiter
			num_store(:, iset) = x.*den_store(:, iset) + grad;
		end
		x = x + grad ./ den_osps;
		x = max(x,pixmin);
		x = min(x,pixmax);
		if subout
			xs(:,idx) = x;
			idx = idx + 1;
		end
	end
	if ~subout, xs(:,iter+1) = x;	end
end
num = sum(num_store, 2);
den = sum(den_store, 2);

if 1,		% optional
	x = num ./ den;
	x = max(x,pixmin); x = min(x,pixmax);
	if subout
		xs(:,idx-1) = x;
	else
		xs(:,iter+1) = x;
	end
end

alphidx = 1;
for iter = (osspsiter+2):niter
	ticker(mfilename, iter, niter)

	relax = relax0 / (1 + relax_rate * (iter-2));

	%
	% loop over subsets
	%
	for iset=1:nblock
		iblock = starts(iset);
		ia = iblock:nblock:na;

		li = Gb{iblock} * x;	% l=G*x "line integrals"
		li = reshape(li, nb, length(ia));
		bel = bi(:,ia) .* exp(-li);
		yb = bel + ri(:,ia);	% predicted meas. means
		dothi = (1 - yi(:,ia) ./ yb) .* bel;

		% optimal curvatures (for ensured monotone increase)
		if streq(curv, 'oc')
			ni = my_trl_curvature(yi(:,ia), bi(:,ia), ri(:,ia), ...
				li, 'oc');
			denom = Gb{iblock}' * col(gi(:,ia) .* ni);
		end

		num = num - num_store(:, iset);
		den = den - den_store(:, iset);

		if isempty(R)
			grad = Gb{iblock}' * dothi(:);
			den_store(:, iset) = max(denom(:,subscale(iset)), ...
						epsilon);
		else
			grad = Gb{iblock}' * dothi(:) - R.cgrad(R, x)/nblock;
			den_store(:, iset) = max(denom(:,subscale(iset)) + ...
						R.denom(R, x)/nblock, epsilon);
		end

		num_store(:, iset) = x.*den_store(:, iset) + grad;

		num = num + num_store(:, iset);
		den = den + den_store(:, iset);

		xold = x;

		x = num ./ den;
		x = max(x,pixmin); x = min(x,pixmax);

		if enhance
			xos = xold + grad ./ (den_store(:,iset) - ...
				denom(:,subscale(iset)) + denom_pc);
			xos = max(xos,pixmin);
			xos = min(xos,pixmax);

			c2 = - sum((xos - x).^2 .* den) / 2;
			c1 = num'*(xos - x) - sum(x.*(xos-x).*den);
			c0 = num'*(x - xold) - sum(x.^2.*den - xold.^2.*den)/2;

			if 1
				if c0+c1+c2 > 0
					alpha = 1;
				else
					alpha = eql_root(-c2, -c1/2, 0.95*c0);
					alpha = min(1, alpha);
					alpha = max(0, alpha);
				end
			else
				alpha = 1/0.9;
				del_obj = -999;
				while (del_obj <=0) & (alpha > 0.0096)
					alpha = alpha * 0.9;
					del_obj = c2*alpha^2 + c1*alpha + c0;
				end

				if alpha <= 0.0096, alpha = 0; end
			end

			x = (1-alpha)*x + alpha*xos;
			alph(alphidx) = alpha;
			alphidx = alphidx + 1;

		end

		if subout
			xs(:,idx) = x;
			idx = idx + 1;
		end
	end

	if chat, printf('Range %g %g', min(x), max(x)), end
	if ~subout, xs(:,iter) = x; end
end

%figure, plot(alph, 'o-')

function obj = trl_obj(alpha, xinc, xos, linc, los, yi, bi, ri, R)

x = alpha*xos + (1 - alpha)*xinc;
li = alpha*los + (1 - alpha)*linc;
yp = bi .* exp(-li) + ri;
like = sum(yi .* log(yp) - yp);
penal = sum(R.penal(R, x));
obj = like - penal;
obj = -obj;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%--- help for mdd/cj ---
%
%mdd/cj is a function.
%    [log_mdd, extras] = cj(obj, ~, opts)
%

author: jmaih


text:
fid = fopen('pairs.txt');
CC = fscanf(fid,'%d %d');
n_set = CC(1);n_num=CC(2);

same_pair = cell(n_set*n_num,2);
diff_pair = cell(n_set*n_num,2);
lfw_label = zeros(n_set*n_num * 2,2);

for i=1:n_set
    for j = 1 : n_num
        CC = textscan(fid, '%s %d %d\n');
        p = CC{1};id1=CC{2};id2=CC{3};
        same_pair((i-1)*n_num + j,1) = {sprintf('%s/%s/%s_%04d.jpg',pwd,p{1},p{1},id1)};
        same_pair((i-1)*n_num + j,2) = {sprintf('%s/%s/%s_%04d.jpg',pwd,p{1},p{1},id2)};
        if exist('list','var')
            lfw_label((i-1)*n_num + j,1) = find(strcmp(list, sprintf('%s_%04d.jpg',p{1},id1)));
            lfw_label((i-1)*n_num + j,2) = find(strcmp(list, sprintf('%s_%04d.jpg',p{1},id2)));
        end;
    end;
    for j = 1 : n_num
         CC = textscan(fid, '%s %d %s %d\n');
         p1 = CC{1};id1=CC{2};p2=CC{3};id2=CC{4};
        diff_pair((i-1)*n_num + j,1) = {sprintf('%s/%s/%s_%04d.jpg',pwd,p1{1},p1{1},id1)};
        diff_pair((i-1)*n_num + j,2) = {sprintf('%s/%s/%s_%04d.jpg',pwd,p2{1},p2{1},id2)};
        if exist('list','var')
            lfw_label(n_set*n_num + (i-1)*n_num + j,1) = find(strcmp(list, sprintf('%s_%04d.jpg',p1{1},id1)));
            lfw_label(n_set*n_num + (i-1)*n_num + j,2) = find(strcmp(list, sprintf('%s_%04d.jpg',p2{1},id2)));
        end;
    end;
end;
fclose(fid);

if exist('feature','var')
    AllFeature1 = feature(:,lfw_label(:,1));
    AllFeature2 = feature(:,lfw_label(:,2));
end;

author: happynear

repo: FaceVerification


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.
%
%
for i=1:numel(X)
    for j=1:p
        X_poly(i,j)=X(i,:).^j;
    end
end

% =========================================================================

end

author: Ayatans

repo: Machine-Learning-homework


text:
%SETUSER Set the user field of a dataset
%
%   A = SETUSER(A,S,FIELD)
%
% INPUT
%   A      Dataset
%   S      Variable to be stored in the user field
%   FIELD  Desired field, default 'USER'.
%
% OUTPUT
%   A      Updated dataset
%
% DESCRIPTION
% Set or reset the subfield FIELD of the user field of A by S.
%
% Note the the USER field of datasets was originally intended for a user
% defined description of datasets. Later its usage was extended to a field
% for storing general information on datasets. For that reason 'old'
% datasets without a structure in the user field are transformed such that
% this information is stored in a subfield USER in the user field. It can
% be retrieved by GETUSER(A,'USER').
%
% Note also that for reasons of backward compatibility the parameter order
% of the SETUSER command differs from similar Matlab commands like
% SETFIELD: first field content, then field name.

author: marianux

repo: ecg-kit

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


text:
function bbox_reg = spp_train_bbox_regressor(imdb, roidb, spp_model, varargin)
% bbox_reg = spp_train_bbox_regressor(imdb, roidb, spp_model, varargin)
%   Trains a bounding box regressor on the image database imdb
%   for use with the SPP model spp_model. The regressor is trained
%   using ridge regression.
%
%   Keys that can be passed in:
%
%   min_overlap     Proposal boxes with this much overlap or more are used
%   layer           The CNN layer features to regress from (either 5, 6 or 7)
%   lambda          The regularization hyperparameter in ridge regression
%   robust          Throw away examples with loss in the top [robust]-quantile
%   binarize        Binarize features or leave as real values >= 0
%
% Adapted from spp code written by Ross Girshick
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

ip = inputParser;
ip.addRequired('imdb',       @iscell);
ip.addRequired('spp_model',  @isstruct);
ip.addParamValue('min_overlap', 0.5,   @isscalar);
ip.addParamValue('layer',       5,     @isscalar);
ip.addParamValue('lambda',      4000,  @isscalar);
ip.addParamValue('robust',      0,     @isscalar);
ip.addParamValue('binarize',    false, @islogical);

ip.parse(imdb, spp_model, varargin{:});
opts = ip.Results;
opts = rmfield(opts, 'spp_model');
opts = rmfield(opts, 'imdb');
opts.cache_name = spp_model.cache_name;

fprintf('\n\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n');
fprintf('Training options:\n');
disp(opts);
fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n');

imdbs_name = cell2mat(cellfun(@(x) x.name, imdb, 'UniformOutput', false));
conf = spp_config('sub_dir', fullfile(opts.cache_name, imdbs_name));
clss = spp_model.classes;
num_clss = length(clss);

% ------------------------------------------------------------------------
% Get the average norm of the features
opts.feat_norm_mean = spp_feature_stats(imdb, roidb, opts.layer, spp_model);
fprintf('average norm = %.3f\n', opts.feat_norm_mean);
% ------------------------------------------------------------------------

% ------------------------------------------------------------------------
% Get all positive examples
imdbs_name = cell2mat(cellfun(@(x) x.name, imdb, 'UniformOutput', false));
save_file = sprintf('./cachedir/%s/%s/bbox_regressor_XY_layer_%d_overlap_%.5g.mat', ...
                    spp_model.cache_name, imdbs_name, opts.layer, opts.min_overlap);
try
  load(save_file);
  fprintf('Loaded saved positives from ground truth boxes\n');
catch
  [X, Y, O, C] = get_examples(spp_model, imdb, roidb, opts);
  save(save_file, 'X', 'Y', 'O', 'C', '-v7.3');
end
for i = 1:num_clss
  fprintf('%14s has %6d samples\n', spp_model.classes{i}, length(find(C == i)));
end
X = spp_poolX_to_fcX(X, opts.layer, spp_model, false);
X = spp_scale_features(X, opts.feat_norm_mean);
% ------------------------------------------------------------------------

% use ridge regression solved by cholesky factorization
method = 'ridge_reg_chol';

models = cell(num_clss, 1);
for i = 1:num_clss
  fprintf('Training regressors for class %s (%d/%d)\n', ...
      spp_model.classes{i}, i, num_clss);
  I = find(O > opts.min_overlap & C == i);
  Xi = X(:, I);
  if opts.binarize
    Xi = single(Xi > 0);
  end
  Yi = Y(I,:);
  Oi = O(I);
  Ci = C(I);

  % add bias feature
  Xi = cat(1, Xi, ones(1, size(Xi,2), class(Xi)));

  % Center and decorrelate targets
  mu = mean(Yi);
  Yi = bsxfun(@minus, Yi, mu);
  S = Yi'*Yi / size(Yi,1);
  [V, D] = eig(S);
  D = diag(D);
  T = V*diag(1./sqrt(D+0.001))*V';
  T_inv = V*diag(sqrt(D+0.001))*V';
  Yi = Yi * T;

  models{i}.mu = mu;
  models{i}.T = T;
  models{i}.T_inv = T_inv;

  models{i}.Beta = [ ...
    solve_robust(Xi, Yi(:,1), opts.lambda, method, opts.robust) ...
    solve_robust(Xi, Yi(:,2), opts.lambda, method, opts.robust) ...
    solve_robust(Xi, Yi(:,3), opts.lambda, method, opts.robust) ...
    solve_robust(Xi, Yi(:,4), opts.lambda, method, opts.robust)];
end
bbox_reg.models = models;
bbox_reg.training_opts = opts;
save([conf.cache_dir 'bbox_regressor_final'], 'bbox_reg');

% ------------------------------------------------------------------------
function [X, Y, O, C] = get_examples(spp_model, imdb, roidb, opts)
% ------------------------------------------------------------------------
num_classes = length(spp_model.classes);

cls_counts = zeros(num_classes, 1);
for idb = 1:length(imdb)
    for i = 1:length(imdb{idb}.image_ids)
      tic_toc_print('%s: counting %d/%d\n', ...
                    procid(), i, length(imdb{idb}.image_ids));

      d = roidb{idb}.rois(i);
      [max_ov cls] = max(d.overlap, [], 2);
      sel_ex = find(max_ov >= 0.5);
      cls = cls(sel_ex);
      for j = 1:length(sel_ex)
        cls_counts(cls(j)) = cls_counts(cls(j)) + 1;
      end
    end
end
total = sum(cls_counts);
feat_dim = size(spp_model.cnn.layers(1).weights{1},2);
% features
X = zeros(feat_dim, total, 'single');
% target values
Y = zeros(total, 4, 'single');
% overlap amounts
O = zeros(total, 1, 'single');
% classes
C = zeros(total, 1, 'single');
cur = 1;

for idb = 1:length(imdb)
    for i = 1:length(imdb{idb}.image_ids)
      tic_toc_print('%s: pos features %d/%d\n', ...
                    procid(), i, length(imdb{idb}.image_ids));
      d = roidb{idb}.rois(i);
      d.feat = spp_load_cached_poolX_features(spp_model.spp_pooler, spp_model.feat_cache{idb}, ...
          imdb{idb}.name, imdb{idb}.image_ids{i}, d.boxes);

      sel_gt = find(d.class > 0);
      gt_boxes = d.boxes(sel_gt, :);
      gt_classes = d.class(sel_gt);

      max_ov = max(d.overlap, [], 2);
      sel_ex = find(max_ov >= opts.min_overlap);
      ex_boxes = d.boxes(sel_ex, :);

      X(:, cur+(0:length(sel_ex)-1)) = d.feat(:, sel_ex);

      for j = 1:size(ex_boxes, 1)
        ex_box = ex_boxes(j, :);
        ov = boxoverlap(gt_boxes, ex_box);
        [max_ov, assignment] = max(ov);
        gt_box = gt_boxes(assignment, :);
        cls = gt_classes(assignment);

        src_w = ex_box(3) - ex_box(1) + eps;
        src_h = ex_box(4) - ex_box(2) + eps;
        src_ctr_x = ex_box(1) + 0.5*src_w;
        src_ctr_y = ex_box(2) + 0.5*src_h;

        gt_w = gt_box(3) - gt_box(1) + eps;
        gt_h = gt_box(4) - gt_box(2) + eps;
        gt_ctr_x = gt_box(1) + 0.5*gt_w;
        gt_ctr_y = gt_box(2) + 0.5*gt_h;

        dst_ctr_x = (gt_ctr_x - src_ctr_x) * 1/src_w;
        dst_ctr_y = (gt_ctr_y - src_ctr_y) * 1/src_h;
        dst_scl_w = log(gt_w / src_w);
        dst_scl_h = log(gt_h / src_h);

        target = [dst_ctr_x dst_ctr_y dst_scl_w dst_scl_h];

        if 0
          % debugging visualizations and checks
          im = imread(imdb.image_at(i));
          showboxesc(im, gt_box, 'g', '-');
          showboxesc([], ex_box, 'r', '-');
          hold on;
          plot(gt_ctr_x, gt_ctr_y, 'gd');
          plot(src_ctr_x, src_ctr_y, 'rd');
          hold off;
          fprintf('target = [%.3f %.3f %.3f %.3f]\n', target(1), target(2), target(3), target(4));
          fprintf('cls = %s\n', spp_model.classes{cls});

          % check that we can correctly reconstruct the gt_box from the
          % gold-standard target
          pred_ctr_x = (target(1) * src_w) + src_ctr_x;
          pred_ctr_y = (target(2) * src_h) + src_ctr_y;
          pred_w = exp(target(3)) * src_w;
          pred_h = exp(target(4)) * src_h;
          pred_box = [pred_ctr_x - 0.5*pred_w, pred_ctr_y - 0.5*pred_h, ...
                      pred_ctr_x + 0.5*pred_w, pred_ctr_y + 0.5*pred_h];
          disp(pred_box);
          disp(gt_box);
          assert(sum(abs(pred_box - gt_box)) < 0.0001);

          pause;
        end

        assert(cur <= total);
        Y(cur, :) = target;
        O(cur)    = max_ov;
        C(cur)    = cls;
        cur = cur + 1;
      end
    end
end

% ------------------------------------------------------------------------
function [x, losses] = solve_robust(A, y, lambda, method, qtile)
% ------------------------------------------------------------------------
[x, losses] = solve(A, y, lambda, method);
fprintf('loss = %.3f\n', mean(losses));
if qtile > 0
  thresh = quantile(losses, 1-qtile);
  I = find(losses < thresh);
  [x, losses] = solve(A(I,:), y(I), lambda, method);
  fprintf('loss (robust) = %.3f\n', mean(losses));
end

% ------------------------------------------------------------------------
function [x, losses] = solve(A, y, lambda, method)
% ------------------------------------------------------------------------

%tic;
switch method
case 'ridge_reg_chol'
  % solve for x in min_x ||Ax - y||^2 + lambda*||x||^2
  %
  % solve (A'A + lambdaI)x = A'y for x using cholesky factorization
  % R'R = (A'A + lambdaI)
  % R'z = A'y  :  solve for z  =>  R'Rx = R'z  =>  Rx = z
  % Rx = z     :  solve for x
  R = chol(A*A' + lambda*eye(size(A,1)));
  z = R' \ (A*y);
  x = R \ z;
case 'ridge_reg_inv'
  % solve for x in min_x ||Ax - y||^2 + lambda*||x||^2
  x = inv(A*A' + lambda*eye(size(A,1)))*A*y;
case 'ls_mldivide'
  % solve for x in min_x ||Ax - y||^2
  if lambda > 0
    warning('ignoring lambda; no regularization used');
  end
  x = A'\y;
end
%toc;
losses = 0.5 * (A'*x - y).^2;

author: ShaoqingRen


text:
function linpack_s_test24 ( )

%*****************************************************************************80
%
%% TEST24 tests SSICO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 100;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST24\n' );
  fprintf ( 1, '  For a symmetric indefinite matrix,\n' );
  fprintf ( 1, '  SSICO estimates the reciprocal condition number.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Assign values to the matrix A.
%
  a(1:n,1:n) = 0.0;

  for i = 1 : n
    a(i,i) = 2.0;
    if ( i < n )
      a(i,i+1) = -1.0;
    end
  end
%
%  Estimate the condition.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimate the condition.\n' );

  [ a, ipvt, rcond, z ] = ssico ( a, lda, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimated reciprocal condition = %f\n', rcond );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% demonstrate usage of mean functions
%
% See also meanFunctions.m.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2014-12-08.
%                                      File automatically generated using noweb.
clear all, close all
n = 5; D = 2; x = randn(n,D);            % create a random data set

% set up simple mean functions
m0 = {'meanZero'};  hyp0 = [];      % no hyperparameters are needed
m1 = {'meanOne'};   hyp1 = [];      % no hyperparameters are needed
mc = {@meanConst};  hypc = 2;  % also function handles are possible
ml = {@meanLinear}; hypl = [2;3];              % m(x) = 2*x1 + 3*x2
mp = {@meanPoly,2}; hypp = [1;1;2;3];  % m(x) = x1+x2+2*x1^2+3*x2^2
mn = {@meanNN,[1,0; 0,1],[0.9,0.5]}; hypn = [];  % nearest neighbor
s = 12; hypd = randn(s,1);           % discrete mean with 12 hypers
md = {'meanDiscrete',s};
hyp.cov = [0;0]; hypg = [];                    % GP predictive mean
xt = randn(2*n,D); yt = sign(xt(:,1)-xt(:,2));      % training data
mg = {@meanGP,hyp,@infEP,@meanZero,@covSEiso,@likErf,xt,yt};
hype = [0;0; log(0.1)];             % regression GP predictive mean
xt = randn(2*n,D); yt = xt(:,1).*xt(:,2);           % training data
me = {@meanGPexact,@meanZero,@covSEiso,xt,yt};

% set up composite mean functions
msc = {'meanScale',{m1}};      hypsc = [3; hyp1];      % scale by 3
msu = {'meanSum',{m0,mc,ml}};  hypsu = [hyp0; hypc; hypl];    % sum
mpr = {@meanProd,{mc,ml}};     hyppr = [hypc; hypl];      % product
mpo = {'meanPow',3,msu};       hyppo = hypsu;         % third power
mask = [false,true];     % mask excluding all but the 2nd component
mma = {'meanMask',mask,ml};    hypma = hypl(mask);
mpf = {@meanPref,ml};          hyppf = 2;  % linear pref with slope

% 0) specify mean function
% mean = md; hyp = hypd; x = randi([1,s],n,1);
% mean = mn; hyp = hypn;
% mean = mg; hyp = hypg;
mean = me; hyp = hype;
% mean = m0;  hyp = hyp0;
% mean = msu; hyp = hypsu;
% mean = mpr; hyp = hyppr;
% mean = mpo; hyp = hyppo;
% mean = mpf; hyp = hyppf;

% 1) query the number of parameters
feval(mean{:})

% 2) evaluate the function on x
feval(mean{:},hyp,x)

% 3) compute the derivatives w.r.t. to hyperparameter i
i = 2; feval(mean{:},hyp,x,i)

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [strCell,IND,aIND] = prtUtilRemoveStrCell(strCell,str)
%[strCell,IND] = prtUtilRemoveStrCell(strCell,str);
%
%   Returns a cell array of strings NOT containing the character arrays in strCell
%   which contain the string str.  IND is the locations of the cells kept.
%
%   see also: prtUtilKeepStrCell
%
% xxx Need Help xxx

[garbage,aIND] = prtUtilKeepStrCell(strCell,str);

IND = setdiff(1:length(strCell),aIND);
strCell = {strCell{IND}};

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function obj = update(obj)
	% Biafra Ahanonu
	% Started: 2021.03.25 [22:11:25] (branched from ciatah.m)

	% Notify user if they are behind version-wise.
	ciapkg.io.updatePkg;
	uiwait(ciapkg.overloaded.msgbox('The CIAtah GitHub website will open. Click "Clone or download" button to download most recent version of CIAtah.'))
	web(obj.githubUrl);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
%CODEGENERATOR.GETPIBUGFIXSTRING Returns a string to fix PI-Bug in auto genereated functions.
%
% bfixString = cGen.getPiBugfixString() Is the string with explanation comment
% and variable declaration as described below.
%
% Notes::
%  In some versions the symbolic toolbox writes the constant $pi$ in
%  capital letters. This way autogenerated functions might not work properly.
%  To fix this issue a local variable is introduced:
%  PI = pi
%
% Authors::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.constructheaderstring, CodeGenerator.replaceheader.

% Copyright (C) 1993-2012, by Peter I. Corke
% Copyright (C) 2012-2013, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function bfixString = getpibugfixstring(CGen)
bfixString = [];
bfixString = [bfixString, sprintf('%s\n','%% Bugfix')];
bfixString = [bfixString, sprintf('%s\n','%  In some versions the symbolic toolbox writes the constant $pi$ in')];
bfixString = [bfixString, sprintf('%s\n','%  capital letters. This way autogenerated functions might not work properly.')];
bfixString = [bfixString, sprintf('%s\n','%  To fix this issue a local variable is introduced:')];
bfixString = [bfixString, sprintf('%s\n','PI = pi;')];
bfixString = [bfixString, sprintf('%s\n','   ')];

author: petercorke

repo: robotics-toolbox-matlab


text:
%
%  Sort the terms of a multivariate polynomial f according to
%  the lexicographical order, and combine the like terms
%
%    Syntax:  >> G = mvPolynSort(F)
%             >> G = mvPolynSort(F, tord)
%             >> [G,idg] = mvPolynSort(F)
%             >> [G,idg] = mvPolynSort(F, tord)
%
%    Input:   F -- (matrix) multivariate polynomial in coeff. matrix
%          tord -- (string, optional) 'lex' or 'grlex', term order type
%                   The default is 'lex' if not provided
%
%   Output:   G --- (matrix) sorted coeff. matrix of the input polynomial
%           idg --- (vector) the indices of G terms
%
%   Examples:
%
%   >> F
%
%   F =
%
%     0     0     4     3     3     0
%     0     0     1     0     0     1
%     0     0     1     0     0     2
%     8     7    -7     2     4    -5
%
%   >> [g,idg] = mvPolynSort(F)
%
%   G =
%
%     0     3     4     0
%     0     0     1     1
%     0     0     1     2
%    15     6    -7    -5
%
%   idg =
%
%     1     4    20    26
%
%   >> [G,idg] = mvPolynSort(F,'grlex')
%
%   G =
%
%     0     0     3     4
%     0     1     0     1
%     0     2     0     1
%    15    -5     6    -7
%
%
%   idg =
%
%     1    12    20    80
%

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%% Copyright (C) 2014, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defop  Method   @@symfun mldivide {(@var{f}, @var{g})}
%% @defopx Operator @@symfun {@var{f} \ @var{g}} {}
%% Backslash division of symbolic functions.
%%
%% Simple example:
%% @example
%% @group
%% syms f(x)
%% g(x) = sin(x);
%% @end group
%%
%% @group
%% h = f \ g
%%   @result{} h(x) = (symfun)
%%       sin(x)
%%       ──────
%%        f(x)
%% @end group
%% @end example
%%
%% @seealso{@@sym/mldivide}
%% @end defop

function h = mldivide(f, g)
  [vars, s1, s2] = helper_symfun_binops(f, g);
  h = symfun(s1 \ s2, vars);
end

%!test
%! syms x
%! f(x) = x^2;
%! assert( isa(f\f, 'symfun'))
%! assert( isa(f\x, 'symfun'))

repo: octsympy


text:
function varargout = WFT(varargin)
% Function:             Windowed Fourier transform (interface)
% Initially Developed:  Dr Qian Kemao (16 May 2009)
% Last modified:        Dr Qian Kemao (17 May 2009)
% Version:              1.0
% Copyrights:           All rights reserved.
% Contact:              mkmqian@ntu.edu.sg (Dr Qian Kemao)

% pushbutton_CallWFTKernel M-file for pushbutton_CallWFTKernel.fig
%      pushbutton_CallWFTKernel, by itself, creates a new pushbutton_CallWFTKernel or raises the existing
%      singleton*.
%
%      H = pushbutton_CallWFTKernel returns the handle to a new pushbutton_CallWFTKernel or the handle to
%      the existing singleton*.
%
%      pushbutton_CallWFTKernel('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in pushbutton_CallWFTKernel.M with the given input arguments.
%
%      pushbutton_CallWFTKernel('Property','Value',...) creates a new pushbutton_CallWFTKernel or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before WFT_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to WFT_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help pushbutton_CallWFTKernel

% Last Modified by GUIDE v2.5 20-May-2009 08:55:16

% Begin initialization code - DO NOT EDIT

%sturcture of g:
%f:                 original data
%filtered:          filtered version of f
%wrapped:           wrapped phase
%unwrapped:         unwrapped phase
%PathMap:           unwrapping path
%g.wx:              wx from wfr
%g.wy:              wy from wfr
%g.r:               r from wfr
%g.phase:           phase from wfr
%g.FringeType:      1~4
%g.AlgorithmType:   'wff' or 'wfr'
%g.a:               background intensity

gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @WFT_OpeningFcn, ...
                   'gui_OutputFcn',  @WFT_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before pushbutton_CallWFTKernel is made visible.
function WFT_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to pushbutton_CallWFTKernel (see VARARGIN)

% Choose default command line output for pushbutton_CallWFTKernel
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes pushbutton_CallWFTKernel wait for user response (see UIRESUME)
% uiwait(handles.figure_WFT);

% --- Outputs from this function are returned to the command line.
function varargout = WFT_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton_LoadMat.
function pushbutton_LoadMat_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_LoadMat (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename, pathname, filterindex] = uigetfile( ...
        '*.mat',  'MATLAB data (*.mat)','Pick a file');
if filterindex~=0
    load (strcat(pathname,filename));
    g.filtered=g.f;
    g.unwrapped=[];
    g.PathMap=[];
    imagesc(angle(g.f))
    colormap gray;
    save result g;
end

% --- Executes on button press in pushbutton_ReadfII.
function pushbutton_ReadfII_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ReadfII (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename, pathname, filterindex] = uigetfile( ...
        {'*.jpg;*.bmp;*.tif;', 'Image files (*.jpg, *.bmp, *.tif)'; ...
        '*.jpg',  'jpg (*.jpg)'; ...
        '*.bmp','bmp (*.bmp)'; ...
        '*.tif','tif (*.tif)'}, ...
        'Pick a file');
if filterindex~=0
    g.wrapped=double(imread(strcat(pathname,filename)));
    g.wrapped=g.wrapped(:,:,1);
    valmax=max(max(g.wrapped));
    valmin=min(min(g.wrapped));
    g.wrapped=(g.wrapped-valmin)/(valmax-valmin)*2*pi-pi;
    g.f=exp(sqrt(-1)*g.wrapped);
    g.filtered=g.f;
    g.unwrapped=[];
    g.PathMap=[];
    g.FringeType=2;
    imagesc(angle(g.f))
    colormap gray;
    save result g;
end

% --- Executes on button press in pushbutton_ReadfIII.
function pushbutton_ReadfIII_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ReadfIII (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename, pathname, filterindex] = uigetfile( ...
        {'*.jpg;*.bmp;*.tif;', 'Image files (*.jpg, *.bmp, *.tif)'; ...
        '*.jpg',  'jpg (*.jpg)'; ...
        '*.bmp','bmp (*.bmp)'; ...
        '*.tif','tif (*.tif)'}, ...
        'Pick a file');
if filterindex~=0
    g.f=double(imread(strcat(pathname,filename)));
    g.f=g.f(:,:,1);
    g.a=mean2(g.f);
    g.f=g.f-g.a;
    g.filtered=g.f;
    g.unwrapped=[];
    g.PathMap=[];
    g.FringeType=3;
    imagesc(g.f)
    colormap gray;
    save result g;
end

% --- Executes on button press in pushbutton_ReadfIV.
function pushbutton_ReadfIV_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ReadfIV (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename, pathname, filterindex] = uigetfile( ...
        {'*.jpg;*.bmp;*.tif;', 'Image files (*.jpg, *.bmp, *.tif)'; ...
        '*.jpg',  'jpg (*.jpg)'; ...
        '*.bmp','bmp (*.bmp)'; ...
        '*.tif','tif (*.tif)'}, ...
        'Pick a file');
if filterindex~=0
    g.f=double(imread(strcat(pathname,filename)));
    g.f=g.f(:,:,1);
    g.a=mean2(g.f);
    g.f=g.f-g.a;
    g.filtered=g.f;
    g.unwrapped=[];
    g.PathMap=[];
    g.FringeType=4;
    imagesc(g.f)
    colormap gray;
    save result g;
end

% --- Executes on button press in pushbutton_CallWFTKernel.
function pushbutton_CallWFTKernel_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_CallWFTKernel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
WFTKernel;

% --- Executes on button press in pushbutton_qg.
function pushbutton_qg_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_qg (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load result;
[g.unwrapped g.PathMap]=unwrapping_qg_trim(g.filtered);
imagesc(g.unwrapped)
save result g;

% --- Executes on button press in pushbutton_ViewOriginalPhase.
function pushbutton_ViewOriginalPhase_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ViewOriginalPhase (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load result;
imagesc(angle(g.f));
H=findobj('Tag','checkbox_SaveOption');
val=get(H,'Value');
if val==1 & ~isempty(g.f)
    imwritescale(angle(g.f),'OriginalPhase.jpg');
end

% --- Executes on button press in pushbutton_ViewOriginalAmplitude.
function pushbutton_ViewOriginalAmplitude_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ViewOriginalAmplitude (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load result;
imagesc(abs(g.f));
H=findobj('Tag','checkbox_SaveOption');
val=get(H,'Value');
if val==1 & ~isempty(g.f)
    imwritescale(abs(g.f),'OriginalAmplitude.jpg');
end

% --- Executes on button press in pushbutton_ViewFilteredPhase.
function pushbutton_ViewFilteredPhase_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ViewFilteredPhase (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load result
imagesc(angle(g.filtered));
H=findobj('Tag','checkbox_SaveOption');
val=get(H,'Value');
if val==1 & ~isempty(g.filtered)
    imwritescale(angle(g.filtered),'FilteredPhase.jpg');
end

% --- Executes on button press in pushbutton_ViewFilteredAmplitude.
function pushbutton_ViewFilteredAmplitude_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ViewFilteredAmplitude (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load result
imagesc(abs(g.filtered));
H=findobj('Tag','checkbox_SaveOption');
val=get(H,'Value');
if val==1 & ~isempty(g.filtered)
    imwritescale(abs(g.filtered),'FilteredAmplitude.jpg');
end

% --- Executes on button press in pushbutton_ViewUnwrappedPhase.
function pushbutton_ViewUnwrappedPhase_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ViewUnwrappedPhase (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load result
imagesc(g.unwrapped);
H=findobj('Tag','checkbox_SaveOption');
val=get(H,'Value');
if val==1 & ~isempty(g.unwrapped)
    imwritescale(g.unwrapped,'UnwrappedPhase.jpg');
end

% --- Executes on button press in pushbutton_ViewUnwrappingPath.
function pushbutton_ViewUnwrappingPath_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ViewUnwrappingPath (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load result
imagesc(g.PathMap);
H=findobj('Tag','checkbox_SaveOption');
val=get(H,'Value');
if val==1 & ~isempty(g.PathMap)
    imwritescale(g.PathMap,'UnwrappingPath.jpg');
end

% --- Executes on button press in pushbutton_ViewOriginalFringePattern.
function pushbutton_ViewOriginalFringePattern_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ViewOriginalFringePattern (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load result
imagesc(real(g.f));
H=findobj('Tag','checkbox_SaveOption');
val=get(H,'Value');
if val==1 & ~isempty(g.f)
    imwritescale(real(g.f),'OriginalFringePattern.jpg');
end

% --- Executes on button press in pushbutton_ViewRealPartOfFilteredResult.
function pushbutton_ViewRealPartOfFilteredResult_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ViewRealPartOfFilteredResult (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load result
imagesc(real(g.filtered));
H=findobj('Tag','checkbox_SaveOption');
val=get(H,'Value');
if val==1 & ~isempty(g.filtered)
    imwritescale(real(g.filtered),'RealPartOfFilteredResult.jpg');
end

% --- Executes on button press in pushbutton_ViewCosineOfFilteredPhase.
function pushbutton_ViewCosineOfFilteredPhase_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_ViewCosineOfFilteredPhase (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load result
imagesc(cos(angle(g.filtered)));
H=findobj('Tag','checkbox_SaveOption');
val=get(H,'Value');
if val==1 & ~isempty(g.filtered)
    imwritescale(cos(angle(g.filtered)),'CosineOfFilteredPhase.jpg');
end

% --- Executes on button press in checkbox_SaveOption.
function checkbox_SaveOption_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox_SaveOption (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox_SaveOption

% --- Executes on button press in checkbox_ViewColor.
function checkbox_ViewColor_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox_ViewColor (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox_ViewColor

% --- Executes on button press in checkbox_ViewInColor.
function checkbox_ViewInColor_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox_ViewInColor (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox_ViewInColor
if get(hObject,'Value')==1
    colormap('default');
else
    colormap(gray);
end


text:
function clean_tetgen_mesh(input_file,output_file)
  % CLEAN_TETGEN_MESH TetGen writes ALL triangles to the .mesh file, that is
  % all faces of the tetrahedra in the 3D mesh.
  % clean_tetgen_mesh(input_file,output_file) reads the vertices and tetrahedra
  % in input_file (ignoring triangles) and determines the surfaces triangles
  % then write the vertices, surfaces triangles and tetrahedra to output_file
  %
  % clean_tetgen_mesh(input_file,output_file)
  %
  % Input:
  %   input_file  path to .mesh file containing vertices and tetrahedra
  %   output_file  path to .mesh file to be written
  %
  [V,T,F] = readMESH(input_file);
  F = boundary_faces(T);
  % reverse orientation of triangles, tetgen seems to have backwards
  % orientation
  F = [F(:,3) F(:,2) F(:,1)];
  % rearrange vertices so that vertices on surface come before internal
  % vertices
  [V,T,F] = faces_first(V,T,F);
  writeMESH(output_file,V,T,F);
end

author: alecjacobson

repo: gptoolbox


text:
function [x, f, eflag, output, lambda] = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt)
%QPS_GUROBI  Quadratic Program Solver based on GUROBI.
%   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
%       QPS_GUROBI(H, C, A, L, U, XMIN, XMAX, X0, OPT)
%   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_GUROBI(PROBLEM)
%   A wrapper function providing a standardized interface for using
%   GUROBI to solve the following QP (quadratic programming)
%   problem:
%
%       min 1/2 X'*H*X + C'*X
%        X
%
%   subject to
%
%       L <= A*X <= U       (linear constraints)
%       XMIN <= X <= XMAX   (variable bounds)
%
%   Inputs (all optional except H, C, A and L):
%       H : matrix (possibly sparse) of quadratic cost coefficients
%       C : vector of linear cost coefficients
%       A, L, U : define the optional linear constraints. Default
%           values for the elements of L and U are -Inf and Inf,
%           respectively.
%       XMIN, XMAX : optional lower and upper bounds on the
%           X variables, defaults are -Inf and Inf, respectively.
%       X0 : optional starting value of optimization vector X
%       OPT : optional options structure with the following fields,
%           all of which are also optional (default values shown in
%           parentheses)
%           verbose (0) - controls level of progress output displayed
%               0 = no progress output
%               1 = some progress output
%               2 = verbose progress output
%               3 = even more verbose progress output
%           grb_opt - options struct for GUROBI, value in verbose
%                   overrides these options
%       PROBLEM : The inputs can alternatively be supplied in a single
%           PROBLEM struct with fields corresponding to the input arguments
%           described above: H, c, A, l, u, xmin, xmax, x0, opt
%
%   Outputs:
%       X : solution vector
%       F : final objective function value
%       EXITFLAG : GUROBI exit flag
%           1 = converged
%           0 or negative values = negative of GUROBI exit flag
%           (see GUROBI documentation for details)
%       OUTPUT : GUROBI output struct
%           (see GUROBI documentation for details)
%       LAMBDA : struct containing the Langrange and Kuhn-Tucker
%           multipliers on the constraints, with fields:
%           mu_l - lower (left-hand) limit on linear constraints
%           mu_u - upper (right-hand) limit on linear constraints
%           lower - lower bound on optimization variables
%           upper - upper bound on optimization variables
%
%   Note the calling syntax is almost identical to that of QUADPROG
%   from MathWorks' Optimization Toolbox. The main difference is that
%   the linear constraints are specified with A, L, U instead of
%   A, B, Aeq, Beq.
%
%   Calling syntax options:
%       [x, f, exitflag, output, lambda] = ...
%           qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt)
%
%       x = qps_gurobi(H, c, A, l, u)
%       x = qps_gurobi(H, c, A, l, u, xmin, xmax)
%       x = qps_gurobi(H, c, A, l, u, xmin, xmax, x0)
%       x = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt)
%       x = qps_gurobi(problem), where problem is a struct with fields:
%                       H, c, A, l, u, xmin, xmax, x0, opt
%                       all fields except 'c', 'A' and 'l' or 'u' are optional
%       x = qps_gurobi(...)
%       [x, f] = qps_gurobi(...)
%       [x, f, exitflag] = qps_gurobi(...)
%       [x, f, exitflag, output] = qps_gurobi(...)
%       [x, f, exitflag, output, lambda] = qps_gurobi(...)
%
%
%   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
%       H = [   1003.1  4.3     6.3     5.9;
%               4.3     2.2     2.1     3.9;
%               6.3     2.1     3.5     4.8;
%               5.9     3.9     4.8     10  ];
%       c = zeros(4,1);
%       A = [   1       1       1       1;
%               0.17    0.11    0.10    0.18    ];
%       l = [1; 0.10];
%       u = [1; Inf];
%       xmin = zeros(4,1);
%       x0 = [1; 0; 0; 1];
%       opt = struct('verbose', 2);
%       [x, f, s, out, lambda] = qps_gurobi(H, c, A, l, u, xmin, [], x0, opt);
%
%   See also QPS_MASTER, GUROBI.

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

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

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

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

%% set up options struct for Gurobi
if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt)
    g_opt = gurobi_options(opt.grb_opt);
else
    g_opt = gurobi_options;
end
if verbose > 1
    g_opt.LogToConsole = 1;
    g_opt.OutputFlag = 1;
    if verbose > 2
        g_opt.DisplayInterval = 1;
    else
        g_opt.DisplayInterval = 100;
    end
else
    g_opt.LogToConsole = 0;
    g_opt.OutputFlag = 0;
end

if ~issparse(A)
    A = sparse(A);
end
if issparse(c);
    c = full(c);
end

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

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

%% set up model
m.A     = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
m.rhs   = [ u(ieq);    u(ilt);    -l(igt);    u(ibx);    -l(ibx)    ];
m.sense = char([ double('=')*ones(1,neq) double('<')*ones(1,niq) ]);
m.lb = xmin;
m.ub = xmax;
m.obj = c';

%% call the solver
if isempty(H) || ~any(any(H))
    lpqp = 'LP';
else
    lpqp = 'QP';
    if ~issparse(H)
        H = sparse(H);
    end
    m.Q = 0.5 * H;
end
if verbose
    alg_names = {
        'automatic',
        'primal simplex',
        'dual simplex',
        'interior point',
        'concurrent',
        'deterministic concurrent',
        'deterministic concurrent simplex'
    };
    vn = gurobiver;
    fprintf('Gurobi Version %s -- %s %s solver\n', ...
        vn, alg_names{g_opt.Method+2}, lpqp);
end
results = gurobi(m, g_opt);
switch results.status
    case 'LOADED',          %% 1
        eflag = -1;
    case 'OPTIMAL',         %% 2, optimal solution found
        eflag = 1;
    case 'INFEASIBLE',      %% 3
        eflag = -3;
    case 'INF_OR_UNBD',     %% 4
        eflag = -4;
    case 'UNBOUNDED',       %% 5
        eflag = -5;
    case 'CUTOFF',          %% 6
        eflag = -6;
    case 'ITERATION_LIMIT', %% 7
        eflag = -7;
    case 'NODE_LIMIT',      %% 8
        eflag = -8;
    case 'TIME_LIMIT',      %% 9
        eflag = -9;
    case 'SOLUTION_LIMIT',  %% 10
        eflag = -10;
    case 'INTERRUPTED',     %% 11
        eflag = -11;
    case 'NUMERIC',         %% 12
        eflag = -12;
    case 'SUBOPTIMAL',      %% 13
        eflag = -13;
    case 'INPROGRESS',      %% 14
        eflag = -14;
    otherwise,
        eflag = 0;
end
output = results;

%% check for empty results (in case optimization failed)
if ~isfield(results, 'x') || isempty(results.x)
    x = NaN(nx, 1);
    lam.lower   = NaN(nx, 1);
    lam.upper   = NaN(nx, 1);
else
    x = results.x;
    lam.lower   = zeros(nx, 1);
    lam.upper   = zeros(nx, 1);
end
if ~isfield(results, 'objval') || isempty(results.objval)
    f = NaN;
else
    f = results.objval;
end
if ~isfield(results, 'pi') || isempty(results.pi)
    pi  = NaN(length(m.rhs), 1);
else
    pi  = results.pi;
end
if ~isfield(results, 'rc') || isempty(results.rc)
    rc  = NaN(nx, 1);
else
    rc  = results.rc;
end

kl = find(rc > 0);   %% lower bound binding
ku = find(rc < 0);   %% upper bound binding
lam.lower(kl)   =  rc(kl);
lam.upper(ku)   = -rc(ku);
lam.eqlin   = pi(1:neq);
lam.ineqlin = pi(neq+(1:niq));
mu_l        = zeros(nA, 1);
mu_u        = zeros(nA, 1);

%% repackage lambdas
kl = find(lam.eqlin > 0);   %% lower bound binding
ku = find(lam.eqlin < 0);   %% upper bound binding

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

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

lambda = struct( ...
    'mu_l', mu_l, ...
    'mu_u', mu_u, ...
    'lower', lam.lower, ...
    'upper', lam.upper ...
);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [ p ] = normalize_vec( x )
p = x / sum(x);
end

author: jacobeisenstein

repo: SAGE


text:
function val = readStringFromFID(FID, string)

% READSTRINGFROMFID Read an boolean from an FID.
% FORMAT
% DESC reads a string from a stream.
% ARG FID : stream to read from.
% ARG name : name of string.
% RETURN val : value of variable in file.
%
% SEEALSO : writeStringToFID, readIntFromFID, readBoolFromFID, readDoubleFromFID
%
% COPYRIGHT : Neil D. Lawrence, 2008

% NDLUTIL

lineStr = getline(FID);
tokens = tokenise(lineStr, '=');
if(length(tokens)~=2 | ~strcmp(tokens{1}, string))
  error('Incorrect file format.')
end
val = tokens{2};

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text: %MTIMES Datafile overload


text:
function planC = saveDVHMatrix(DVHNum, doseBinsV, volsHistV, planC)
%"saveDVHMatrix"
%   Store the doseBinsV and volsHistV with the binWidth shift removed.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.
%
%Usage:
%   function planC = saveDVHMatrix(DVHNum, doseBinsV, volsHistV, planC)
if isempty(doseBinsV) || isempty(volsHistV)
    error('Cannot save DVH for empty structures');
end
indexS = planC{end};

if length(doseBinsV) == 1
    doseBinsV(2) = 2*doseBinsV(1);
    volsHistV(2) = 0;
end

%Calculate width of bins, use to find middle of each bin.
binWidthsV = diff(doseBinsV);
lastBinWidth = binWidthsV(end);
binWidthsV(end+1) = lastBinWidth;

%Extract dose bin values, adding half to binwidth to get lower edge.
planC{indexS.DVH}(DVHNum).DVHMatrix(:,1) = doseBinsV - binWidthsV/2;
planC{indexS.DVH}(DVHNum).DVHMatrix(:,2) = volsHistV;
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR