text:
 function [C, ugibb] = Csparse(arg0, arg1, arg2, arg3)
%function C = Csparse('b2info', kappa, mask, chat)
%function C = Csparse('center', Cscale, mask, chat)
%function C = Csparse('maskleak', mask, chat)
%function C = Csparse('maskleak3', [beta_xy beta_z], mask, chat)
%function C = Csparse('maskbar', mask, chat)
%function C = Csparse('test')
%function Csparse('plot', C, mask)
% create 2D (or 1D or 3D) penalty matrix C
% R = beta * C'C is the penalty Hessian (in quadratic case)
% in
%	kappa	[nx,ny]		as in :srp paper, i.e.
%				sqrt( (A.^2)'*nder2(:) ./ sum(A.^2)' )
%	mask	[nx,ny[,nz]]	support
% out
%	C	[nc,np]		np = # of pixels in mask
%
% maskleak: penalize across mask boundary, like in ASPIRE (C*1 != 0)
%	(matches "-" penalty in ASPIRE)
% maskbar: barrier across mask boundary (C*1 = 0)
% b2info: as in 'b2info' penalty in ASPIRE
% center: as in 'center' penalty in ASPIRE

if nargin == 1 && streq(arg0, 'test'), Csparse_test, return, end
if nargin < 1, ir_usage, end

warning 'Csparse is obsolete.  C2sparse is recommended instead!'

betas = 1;

if strcmp(arg0, 'maskleak3')
	Cscale = 1;
	mask = arg2 ~= 0;
	betas = arg1;
	chat = 0;		if 1==exist('arg3'), chat = arg3; end

elseif strcmp(arg0, 'maskleak') || strcmp(arg0, 'maskbar')
	Cscale = 1;
	mask = arg1 ~= 0;
	chat = 0;		if 1==exist('arg2'), chat = arg2; end

elseif strcmp(arg0, 'center')
	Cscale = arg1;
	mask = arg2 ~= 0;
	chat = 0;		if 1==exist('arg3'), chat = arg3; end
	arg0 = 'maskleak';	% trick!

elseif strcmp(arg0, 'b2info')
	kappa = arg1;
	mask = kappa~=0;	if 1==exist('arg2'), mask = arg2; end
	chat = 0;		if 1==exist('arg3'), chat = arg3; end
	if ~islogical(mask), error 'mask must be logical', end
	if any(mask(:,1)) || any(mask(:,end)) || any(mask(1,:)) || any(mask(end,:))
		error 'mask must not include outer row/col'
	end

%
% test by calculating penalty for an image here and in aspire
%
elseif streq(arg0, 'test')
	nx = 20; ny = 12;	nb = nx+1; na = ny;
	mask = zeros(nx,ny);	mask(4:14,4:8) = ones(11,5);
	mask = mask ~= 0;

	A = [spdiag(mask(:), 'nowarn'); zeros(ny,nx*ny)];
	rng(0)
	yy = rand(nb,na);
	xx = mask .* rand(nx,ny);
	dd = 1 + 1 * rand(nb,na);
	kappa = mask+0;
	jb = mask(:);
	kappa(jb) = sqrt( (A(:,jb).^2)'*dd(:) ./ sum(A(:,jb).^2)' );

whos
error('hmm')
	delete('t,A.wtf')
	wtf_write('t,A.wtf', A, nx, ny, nb, na)
	save('t,d', 'dd')
	save('t,y', 'yy')
	save('t,x', 'xx')

	err = yy(:)-A*xx(:);
%	C = Csparse('maskleak', mask);		method = '@1@alg@0,quad,1,-'
	C = Csparse('b2info', kappa, mask, 1);	method = '@1@alg@0,quad,1,b2info'

	printm('like=%.11e penal=%.11e', ...
		err' * spdiag(dd, 'nowarn') * err / 2, norm(C * xx(jb)).^2 / 2)
	eval(['!i best,pwls2 - t,x.mat t,y.mat - t,d.mat t,A.wtf ' method])
	C = 0;
	return

elseif strcmp(arg0, 'plot')
	C = arg1;
	mask = arg2 ~= 0;
%	im(221, kappa, 'kappa')
	if im
		subplot(338), spy(C), title 'C'
	end
	if ncol(C) < 500
		im(339, C'*C, 'C''C')
	else
		im(339, embed(diag(C'*C), mask), 'Diag(C''C)'), cbar
	end
	return

else
	help Csparse
	error([arg0 ' unknown'])
end

	[nx,ny,nz] = size(mask);

%
% C for 1st-order differences
%
if 1
	dx	= sparse(filtmat('1d', [-1 1 0], nx));
	if ~strcmp(arg0, 'maskleak')
		dx(1,1) = 0;
	end
%	dx	= dx(2:end,:);
	dx	= kron(speye(ny*nz), dx);
	C	= dx;	clear dx
	if ny > 1
		dy	= sparse(filtmat('1d', [-1 1 0], ny));
		dy(1,1) = 0;
%		dy	= dy(2:end,:);
		if nz > 1
			dy	= kron(speye(nz), kron(dy, speye(nx)));
		else
			dy	= kron(dy, speye(nx));
		end
		b.xy = betas(1);
		C	= sqrt(b.xy) * [C; dy];	clear dy
	end
	if nz > 1
		dz	= sparse(filtmat('1d', [-1 1 0], nz));
%		dz	= dz(2:end,:);
		dz(1,1) = 0;
		dz	= kron(dz, speye(nx*ny));
		b.z = betas(2);
		C	= [C; sqrt(b.z)*dz];	clear dz
	end
end

if strcmp(arg0, 'maskleak') || strcmp(arg0, 'maskleak3')
	C = Cscale * C(:, mask(:));
	if chat
		Csparse('plot', C, mask);
	end
	return
end

% only rows (and cols) of C where all elements are in mask
if strcmp(arg0, 'maskbar')
	t = abs(C) * mask(:);
	t = full(t == sum(abs(C)')');
	C = C(t,:);
	C = C(:,mask(:));
	return
end

if ~strcmp(arg0, 'b2info'),	error bug, end

	if any(kappa(~mask(:))), error('kappa outside of mask'), end

	%
	% 'fix edges' of kappa, as in ASPIRE
	%
	offset = [-1 1 -nx nx nx-1 nx+1 1-nx -1-nx];
	for jj = find(mask(:))'
		bad = offset(mask(jj + offset) == 0);
		if bad % ~= []
			kappa(jj + bad) = kappa(jj) * ones(size(bad));
		end
	end

	t = abs(C) * spdiag(1:nx*ny, 'nowarn');
	t2 = max(t')';
	t1 = sum(t')' - t2;	% trick works with 1st order only!
	if any(t1 == t2 | t1 <= 0 | t2 <= 0)
		if sum(t1 == 0) ~= nx + ny, error bug_t1, end
		if sum(t2 == 0) ~= nx + ny, error bug_t2, end
		if any(t1(t1 ~= 0) == t2(t2 ~= 0)), error bug_t12, end
		tdiag = zeros(nx*ny,1);
		tdiag(t1 ~= 0) = sqrt(kappa(t1(t1 ~= 0)) .* kappa(t2(t1 ~= 0)));
		% all the above because now we keep zero rows in C...
	else
		tdiag = sqrt(kappa(t1) .* kappa(t2));
	end
	C = spdiag(tdiag, 'nowarn') * C;

	C = C(:,mask(:));

	if chat
		Csparse('plot', C, mask);
	end

	if nargout > 1
		ugibb.x = zeros(nx,ny);
		ugibb.y = zeros(nx,ny);
		ugibb.p = zeros(nx,ny);
		ugibb.n = zeros(nx,ny);

		t = (t2 == t1 + 1) & (tdiag ~= 0);
		sum(t)
		ugibb.x(t2(t)) = tdiag(t);

		t = (t2 == t1 + nx) & (tdiag ~= 0);
		sum(t)
		ugibb.y(t2(t)) = tdiag(t);
		if chat
			im([ugibb.x ugibb.y])
		end
	end

%
% self test
%
function Csparse_test
% run test with 'kappa' for b2info
x = shepplogan(64,64,1); x = x / max(x(:));
mask = true(size(x));
mask([1 end],:) = false; mask(:,[1 end]) = false;
im pl 2 2
im(1, mask+x, 'mask + x'), cbar

kappa = 2 - double(x > 0.5);
kappa = mask .* kappa;
im(2, kappa, 'kappa'), cbar

C = Csparse('b2info', kappa, mask, 0);
if im, im subplot 3, spy(C), end

tmp = embed(sum(abs(C))', mask);
im(4, tmp, '|C|''1'), cbar

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%SERIALLINK.TRCHAIN Convert to elementary transform sequence
%
% S = R.TRCHAIN(OPTIONS) is a sequence of elementary transforms that describe the
% kinematics of the serial link robot arm.  The string S comprises a number
% of tokens of the form X(ARG) where X is one of Tx, Ty, Tz, Rx, Ry, or Rz.
% ARG is a joint variable, or a constant angle or length dimension.
%
% For example:
%        >> mdl_puma560
%        >> p560.trchain
%        ans =
%        Rz(q1)Rx(90)Rz(q2)Tx(0.431800)Rz(q3)Tz(0.150050)Tx(0.020300)Rx(-90)
%        Rz(q4)Tz(0.431800)Rx(90)Rz(q5)Rx(-90)Rz(q6)
%
% Options::
% '[no]deg'    Express angles in degrees rather than radians (default deg)
% 'sym'        Replace length parameters by symbolic values L1, L2 etc.
%
% See also trchain, trotx, troty, trotz, transl, DHFactor.

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

function s = trchain(robot, varargin)

    opt.sym = false;
    opt.deg = true;
    opt = tb_optparse(opt, varargin);

    if opt.deg
        conv = 180/pi;
    else
        conv = 1;
    end

    s = '';
    varcount = 1;

    for j=1:robot.n
        L = robot.links(j);

        if robot.ismdh()
            % Method for modified DH parameters

            % Append Tx(a)
            if L.a ~= 0
                if opt.sym
                    s = append(s, 'Tx(L%d)', varcount);
                    varcount = varcount+1;
                else
                    s = append(s, 'Tx(%g)', L.a);
                end

            end

            % Append Rx(alpha)
            if L.alpha ~= 0
                s = append(s, 'Rx(%g)', (L.alpha*conv));
            end

            if L.isrevolute()

                % Append Tz(d)
                if L.d ~= 0
                    if opt.sym
                        s = append(s, 'Tz(L%d)', varcount);
                        varcount = varcount+1;
                    else
                        s = append(s, 'Tz(%g)', L.d);
                    end
                end

                % Append Rz(q)
                s = append(s, 'Rz(q%d)', j);
            else

                % Append Rz(theta)
                if L.theta ~= 0
                    s = append(s, 'Rz(%g)', (L.alpha*conv));
                end

                % Append Tz(q)
                s = append(s, 'Tz(q%d)', j);
            end

        else
            % Method for standard DH parameters

            if L.isrevolute()

                % Append Rz(q)
                s = append(s, 'Rz(q%d)', j);

                % Append Tz(d)
                if L.d ~= 0
                    if opt.sym
                        s = append(s, 'Tz(L%d)', varcount);
                        varcount = varcount+1;
                    else
                        s = append(s, 'Tz(%g)', L.d);
                    end
                end
            else

                % Append Rz(theta)
                if L.theta ~= 0
                    s = append(s, 'Rz(%g)', (L.alpha*conv));
                end

                % Append Tz(q)
                s = append(s, 'Tz(q%d)', j);
            end

            % Append Tx(a)
            if L.a ~= 0
                if opt.sym
                    s = append(s, 'Tx(L%d)', varcount);
                    varcount = varcount+1;
                else
                    s = append(s, 'Tx(%g)', L.a);
                end

            end

            % Append Rx(alpha)
            if L.alpha ~= 0
                s = append(s, 'Rx(%g)', (L.alpha*conv));
            end
        end
    end
end

    function s = append(s, fmt, j)
        s = strcat(s, sprintf(fmt, j));
    end

author: petercorke

repo: robotics-toolbox-matlab


text:
% TEST_THICL_LSHAPED_MP_G_NMNN: data function for Neumann boundary condition.

function g = test_thick_Lshaped_mp_g_nmnn(x, y, z, ind)
  switch (ind)
   case {1, 4}
     g = cos(z) .* exp(x) .* (sin(x.*y) + y .* cos(x.*y));
   case {2, 3}
     g = -x .* exp(x) .* cos(x.*y) .* cos(z);
   case {5}
     g = -cos(z) .* exp(x) .* (sin(x.*y) + y .* cos(x.*y));
   case {6}
     g = x .* exp(x) .* cos(x.*y) .* cos(z);
   case {7}
     g = exp(x) .* sin(x.*y) .* sin(z);
   case {8}
     g = - exp(x) .* sin(x.*y) .* sin(z);
   otherwise
    error('g_nmnn: error in the reference number for the boundary')
  end
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function color = get_new_color(i)

N = 8;
TWON = 2^N;

global GETNEWCOLOR___;
if isempty(GETNEWCOLOR___),
  GETNEWCOLOR___ = jet(TWON);
end
i = mod(i-1,TWON)+1;

if i == 1,
  color = GETNEWCOLOR___(1,:);
  return;
end
if i == 2,
  color = GETNEWCOLOR___(end,:);
  return;
end

total = 2;
for j = 1:8,
end

v = floor(linspace(1,256,total));
offset = i - (total - j);
l = v(2*offset);

color = GETNEWCOLOR___(l,:);

author: kristinbranson

repo: JAABA


text:
function imdb = setup_voc07_eb(varargin)
% cnn_voc07_eb_setup_data  Initialize PASCAL VOC2007 data with edge
% boxes

% Warning! boxes are in the format of ([y1 x1 y2 x2])

opts.dataDir = fullfile('data') ;
opts.proposalDir = fullfile(opts.dataDir,'EB');
opts.loadTest = 1;
opts = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                 Load selective search win
% -------------------------------------------------------------------------
%% Get selective search windows
files = {'EdgeBoxesVOC2007trainval.mat', ...
  'EdgeBoxesVOC2007test.mat'} ;

if ~exist(opts.proposalDir, 'dir')
  mkdir(opts.proposalDir) ;
end

for i=1:numel(files)
  outPath = fullfile(opts.proposalDir, files{i}) ;
  if ~exist(outPath, 'file')
    url = sprintf('http://groups.inf.ed.ac.uk/hbilen-data/data/WSDDN/%s',files{i}) ;
    fprintf('Downloading %s to %s\n', url, outPath) ;
    urlwrite(url,outPath) ;
  end
end

if ~isempty(opts.proposalDir)
  t1 = load([opts.proposalDir,filesep,files{1}]);
  if opts.loadTest
    t2 = load([opts.proposalDir,filesep,files{2}]);
    ssw.id = [str2double(t1.images) str2double(t2.images)];
    ssw.boxes = cat(2,t1.boxes,t2.boxes);
    ssw.boxScores = cat(2,t1.boxScores,t2.boxScores);
  else
    ssw.id = str2double(t1.images);
    ssw.boxes = t1.boxes;
    ssw.boxScores = t1.boxScores;
  end

  [~,si] = sort(ssw.id);
  ssw.id = ssw.id(si);
  ssw.boxes = ssw.boxes(si);
  ssw.boxScores = ssw.boxScores(si);
end

% -------------------------------------------------------------------------
%                                                  Load categories metadata
% -------------------------------------------------------------------------
cats = {'aeroplane','bicycle','bird','boat','bottle','bus','car',...
  'cat','chair','cow','diningtable','dog','horse','motorbike','person',...
  'pottedplant','sheep','sofa','train','tvmonitor'};

if ~exist(opts.dataDir,'dir')
  error('wrong data folder!');
end

% Download VOC Devkit and data
if ~exist(fullfile(opts.dataDir,'VOCdevkit'),'dir')
  files = {'VOCtest_06-Nov-2007.tar',...
           'VOCtrainval_06-Nov-2007.tar',...
           'VOCdevkit_08-Jun-2007.tar'} ;
  for i=1:numel(files)
    if ~exist(fullfile(opts.dataDir, files{i}), 'file')
      outPath = fullfile(opts.dataDir,files{i}) ;
      url = sprintf('http://host.robots.ox.ac.uk/pascal/VOC/voc2007/%s',files{i}) ;
      fprintf('Downloading %s to %s\n', url, outPath) ;
      urlwrite(url,outPath) ;
      untar(outPath,opts.dataDir);
    end
  end
end
addpath(fullfile(opts.dataDir, 'VOCdevkit', 'VOCcode'));

traindata = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main','train.txt'));
valdata = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main','val.txt'));
testdata = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main','test.txt'));

assert(numel(traindata)==2501);
assert(numel(valdata)==2510);
assert(numel(testdata)==4952);

imdb.classes.name = cats ;
imdb.classes.description = cats ;
imdb.imageDir = fullfile(opts.dataDir, fullfile('VOCdevkit','VOC2007','JPEGImages')) ;

% -------------------------------------------------------------------------
%                                                           Training images
% -------------------------------------------------------------------------%
names = cell(1,numel(traindata));
labels = zeros(numel(traindata),numel(cats));

% load image names
for t=1:numel(traindata)
  names{t} = sprintf('%06d.jpg',traindata(t));
  %   data{t} = imread(sprintf('%s/%s',imdb.imageDir,names{t}));
end

% load binary labels
for c=1:numel(cats)
  t = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main',[cats{c},'_train.txt']));
  labels(:,c) = t(:,2);
end

imdb.images.id = traindata';
imdb.images.name = names ;
imdb.images.set = ones(1, numel(names)) ;
imdb.images.label = labels' ;
% imdb.images.data = data;

% -------------------------------------------------------------------------
%                                                         Validation images
% -------------------------------------------------------------------------

names = cell(1,numel(valdata));
labels = zeros(numel(valdata),numel(cats));
% data = cell(1,numel(valdata));

% load image names
for t=1:numel(valdata)
  names{t} = sprintf('%06d.jpg',valdata(t));
  %   data{t} = imread(sprintf('%s/%s',imdb.imageDir,names{t}));
end

% load binary labels
for c=1:numel(cats)
  t = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main',[cats{c},'_val.txt']));
  labels(:,c) = t(:,2);
end

imdb.images.id = horzcat(imdb.images.id, valdata') ;
imdb.images.name = horzcat(imdb.images.name, names) ;
imdb.images.set = horzcat(imdb.images.set, 2*ones(1,numel(names))) ;
imdb.images.label = horzcat(imdb.images.label, labels') ;
% imdb.images.data = horzcat(imdb.images.data, data) ;

% % -------------------------------------------------------------------------
% %                                                               Test images
% % -------------------------------------------------------------------------
%
%
if opts.loadTest
  names = cell(1,numel(testdata));
  labels = zeros(numel(testdata),numel(cats));
  % data = cell(1,numel(testdata));

  % load image names
  for t=1:numel(testdata)
    names{t} = sprintf('%06d.jpg',testdata(t));
    %   data{t} = imread(sprintf('%s/%s',imdb.imageDir,names{t}));
  end

  % load binary labels
  for c=1:numel(cats)
    t = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main',[cats{c},'_test.txt']));
    labels(:,c) = t(:,2);
  end

  imdb.images.id = horzcat(imdb.images.id, testdata') ;
  imdb.images.name = horzcat(imdb.images.name, names) ;
  imdb.images.set = horzcat(imdb.images.set, 3 * ones(1,numel(names))) ;
  imdb.images.label = horzcat(imdb.images.label, labels') ;
  % imdb.images.data = horzcat(imdb.images.data, data) ;
end
% -------------------------------------------------------------------------
%                                                            Postprocessing
% -------------------------------------------------------------------------
[~,sorti] = sort(imdb.images.id);

imdb.images.id = imdb.images.id(sorti);
imdb.images.name = imdb.images.name(sorti) ;
imdb.images.set = imdb.images.set(sorti) ;
imdb.images.label = single(imdb.images.label(:,sorti)) ;
imdb.images.size = zeros(numel(imdb.images.name),2);

if ~isempty(opts.proposalDir)
  imdb.images.boxes = ssw.boxes;
  imdb.images.boxScores = ssw.boxScores;
  assert(all(ssw.id==imdb.images.id));
end

% this is zero as scores of selective search windows are not much
% informative
if ~isempty(opts.proposalDir)
  % imdb.images.boxScores = cell(size(imdb.images.boxes));
  for i=1:numel(imdb.images.boxes)
    imdb.images.boxes{i} = int16(imdb.images.boxes{i});
    imdb.images.boxScores{i} = single(imdb.images.boxScores{i});

    imf = imfinfo(fullfile(imdb.imageDir,imdb.images.name{i}));
    imdb.images.size(i,:) = [imf.Height,imf.Width];

    maxBoxes = max(imdb.images.boxes{i});
    if imdb.images.size(i,1)< max(maxBoxes([1,3]))
      error('Wrong box coordinates');
    end
    if imdb.images.size(i,2)< max(maxBoxes([2,4]))
      error('Wrong box coordinates');
    end

  end
end
end

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
function a = lt( x, y )

%   Disciplined convex programming information for LT (<):
%      The left-hand side of a less-than constraint must be convex. The
%      right-hand side must be concave. Of course, real constant and
%      affine expressions are both convex and concave and can be used on
%      either side as well.
%
%   Disciplined geometric programming information for LT (<):
%      The left-hand side of a less-than constraint must be log-convex,
%      including positive constants, monomials, posynomials, generalized
%      posynomials, and products thereof. The right-hand side must be
%      log-concave---including positive constants, monomials,
%      reciprocals of log-convex expressions, and products thereof.
%
%   Note that CVX does not distinguish between strict less-than (<) and
%   less-than-or-equal (<=) constraints; they are treated identically.
%   Feasible interior-point solvers tend to return points which satisfy
%   strict inequality, but not all solvers do.

evalin( 'caller', 'cvx_verify' );
b = cvx_pushcnstr( x, y, '<' );
if nargout, a = b; end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
% specify_parameters.m
%
% Run this program first to set up parameters for each data run.  The run
% parameters are stored in the file "CD_run_params.mat" in the current
% directory.  The program "extract_data.m" reads this file and adds some
% fixed radar parameters to the file.  Subsequent programs can read this
% file to obtain the parameters needed for the SAR processing.
%
% The radar data are stored on the CD supplied with the book, but they
% can be read from a disk file copied from the CD to your hard drive.
%
% The only file needed from the CD is "dat_01.001", which contains the
% radar data as well as the chirp replica and auxiliary data for each range
% line.  As none of the parameters like PRF and range gate delay change in
% this data set, and the replica is well defined by a linear FM chirp, the
% only data that are really needed from the CD is the raw radar signal data
% contained in "dat_01.001".  The radar parameters are stored in the
% program "extract_data.m".
%
% The radar data can be written in one block, or broken up into a number of
% blocks, with a set of files for each block.
%
% The first range cell and the number of range cells can be set to any
% value within the boundaries of the data array.  The first range line
% should be set to 1 plus a multiple of 8, and the number of lines per
% block set to a factor of 8.  Also, the number of lines should be set to a
% multiple of the number of lines per block.  If these rules are not
% obeyed, the "extract_data.m" program will quantize the values accordingly.

% Parameters to specify for each run:
% -----------------------------------
% input_path    : location of the radar data file 'dat_01.001'
% output_path   : location for the extracted data files
% output_prefix : prefix to the extracted data file names
%
% first_rg_cell : first range cell to extract
% Nrg_cells     : number of range cells to extract
% first_rg_line : first range line (azimuth sample) to extract
% Nrg_lines     : number of range lines to extract
% Nrg_lines_blk : number of range lines per block
% UseMATfiles   : Set to 0 if you want binary header, data, aux and replica
%                 files (.DAT) rather than just the MATLAB data files (.MAT)
% -------------------------------------------------------------------------
% Created :   Nov 29, 2004 by Kaan Ersahin
% Modified:   Nov 30, 2004 by Kaan Ersahin
% Modified:   Dec 3,  2004 by Ian Cumming
%             - moved fixed parameters to "extract_data.m"
%             - now use a fixed parameter file name for all programs
%             - added UseMATfiles option (=1 for MAT files)
% -------------------------------------------------------------------------

clear,   home,   format compact

% Define location of data files
% input_path      = 'Y:\RSAT\scene01\'; % These values are for runs at UBC
% output_path     = 'Y:\RSAT\EXTRACTED_DATA\';
input_path        = '.\scene01\';
output_path       = '.\scene01\';
CD_data_file_name = 'dat_01.001';
output_prefix     = 'raw';

%==========================================================================
%  Define area of interest of the radar data to be extracted
%  Note that the azimuth parameters are quantized in "extract_data.m"
%
%  Here are some possible values for first_rg_cell and first_rg_line:
%  For coal & ferry terminals, use range =  970, azimuth =  1801
%  For Vancouver airport,      use range = 1060, azimuth =  5561
%  For UBC and one ship,       use range =  760, azimuth =  7097
%  For Stanley Park & city,    use range = 1850, azimuth =  7657
%  For English Bay ships,      use range = 1050, azimuth =  7769
%  For Squamish & Garibaldi,   use range = 2640, azimuth = 16169
%  For Brackendale,            use range = 2800, azimuth = 17897 (max az)

%                              use range = 1350, azimuth =  3170

first_rg_cell   =  1050;    % Define the range limits
first_rg_line   =  7769;    % Define the azimuth limits  (19432 max)
Nrg_cells       =  2048;    % Suggest 2048 cells
Nrg_lines_blk   =  6*256;   % Suggest 1536, it should be larger than the
                            %        size of the azimuth match filter (705)
                            %              so that an image can be created.
Nrg_lines   =  1*Nrg_lines_blk;

UseMATfiles =  1; % (1) Use   MAT  files to save the extracted data
                  % (0) Use binary files to save the extracted data
SaveV6      =  0; % (1) Save data in MATLAB v6 format when using v7
                  % (0) Save data in MATLAB v7 format
                  % This variable is only used if you are running MATLAB 7
%==========================================================================

%  Save the data and run parameters in a matlab data file (*.mat)

vv = version;   vers = str2num(vv(1));  % Find MATLAB version
if vers == 6,   SaveV6 = 0;   end

if SaveV6
  save -v6 CD_run_params
  fprintf('\nData parameters written to:  CD_run_params.mat  in V6 format')
else
  save CD_run_params
  fprintf('\nData parameters written to:  CD_run_params.mat')
end

fprintf('\nNow run "extract_data.m"\n\n')
beep,   pause(0.4),   beep

%% This setup program is now finished, but you can add calls to subsequent
%%  data extraction and processing programs below, if desired.

% extract_data           % Extract the SAR data from the CD

% compute_azim_spectra   % Find the baseband Doppler centroid

% compress_data          % Perform SAR processing

author: joeyos

repo: SAR-imaging

save_path: github-repos/MATLAB/joeyos-SAR-imaging


text:
function b = dsisl ( a, lda, n, kpvt, b )

%*****************************************************************************80
%
%% DSISL solves a real symmetric system factored by DSIFA.
%
%  Discussion:
%
%    To compute inverse(A) * C where C is a matrix with P columns
%
%      call dsifa ( a, lda, n, kpvt, info )
%
%      if ( info == 0 ) then
%        do j = 1, p
%          call dsisl ( a, lda, n, kpvt, c(1,j) )
%        end do
%      end if
%
%    A division by zero may occur if the inverse is requested
%    and DSICO has set RCOND == 0.0D+00 or DSIFA has set INFO /= 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,N), the output from DSIFA.
%
%    Input, integer LDA, the leading dimension of the array A.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer KPVT(N), the pivot vector from DSIFA.
%
%    Input, real B(N), the right hand side.
%
%    Output, real B(N), the solution.
%

%
%  Loop backward applying the transformations and D inverse to B.
%
  k = n;

  while ( 0 < k )

    if ( 0 <= kpvt(k) )
%
%  1 x 1 pivot block.
%
      if ( k ~= 1 )

        kp = kpvt(k);
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end
%
%  Apply the transformation.
%
        b(1:k-1) = daxpy ( k-1, b(k), a(1:k-1,k), 1, b(1:k-1), 1 );

      end
%
%  Apply D inverse.
%
      b(k) = b(k) / a(k,k);
      k = k - 1;

    else
%
%  2 x 2 pivot block.
%
      if ( k ~= 2 )

        kp = abs ( kpvt(k) );
%
%  Interchange.
%
        if ( kp ~= k-1 )
          temp = b(k-1);
          b(k-1) = b(kp);
          b(kp) = temp;
        end
%
%  Apply the transformation.
%
        b(1:k-2) = daxpy ( k-2, b(k), a(1:k-2,k), 1, b(1:k-2), 1 );
        b(1:k-2) = daxpy ( k-2, b(k-1), a(1:k-2,k-1), 1, b(1:k-2), 1 );

      end
%
%  Apply D inverse.
%
      ak = a(k,k) / a(k-1,k);
      akm1 = a(k-1,k-1) / a(k-1,k);
      bk = b(k) / a(k-1,k);
      bkm1 = b(k-1) / a(k-1,k);
      denom = ak * akm1 - 1.0;
      b(k) = ( akm1 * bk - bkm1 ) / denom;
      b(k-1) = ( ak * bkm1 - bk ) / denom;
      k = k - 2;

    end

  end
%
%  Loop forward applying the transformations.
%
  k = 1;

  while ( k <= n )

    if ( 0 <= kpvt(k) )
%
%  1 x 1 pivot block.
%
      if ( k ~= 1 )
%
%  Apply the transformation.
%
        b(k) = b(k) + b(1:k-1)' * a(1:k-1,k);
        kp = kpvt(k);
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end

      end

      k = k + 1;

    else
%
%  2 x 2 pivot block.
%
      if ( k ~= 1 )
%
%  Apply the transformation.
%
        b(k) = b(k) + b(1:k-1) * a(1:k-1,k);
        b(k+1) = b(k+1) + b(1:k-1) * a(1:k-1,k+1);
        kp = abs ( kpvt(k) );
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end

      end

      k = k + 2;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef CBFDT < handle
    properties
        systemParam
        optType
        optParam
        xcurr
        timecurr
        solvertime
        xdim
        udim
        tlog = []
        xlog = []
        ulog = []
        omegalog = []
    end

    methods
        function self = CBFDT(system_param, x0, t0)
            self.systemParam = system_param;
            self.xcurr = x0;
            self.timecurr = t0;
            self.xdim = size(self.systemParam.A, 1);
            self.udim = size(self.systemParam.B, 2);
        end

        function setOpt(self, opt_type, opt_param)
            self.optType = opt_type;
            self.optParam = opt_param;
        end

        function sim(self, total_time)
            xk = self.xcurr;
            while self.timecurr <= total_time
                [feas, xopt, uopt, Jopt] = self.solve();
                if feas ~= 1
                    return;
                end
                uk = uopt(:,1);
                xk = self.systemParam.A * xk + self.systemParam.B * uk;
                self.xcurr = xk;
                self.timecurr = self.timecurr + self.systemParam.timestep;
                self.tlog = [self.tlog, self.timecurr];
                self.xlog = [self.xlog, xk];
                self.ulog = [self.ulog, uk];
            end
        end

        function [feas, xopt, uopt, Jopt] = solve(self)
            if strcmp(self.optType,'dclfdcbf')
                [feas, xopt, uopt, Jopt] = solveDCLFDCBF(self);
            elseif strcmp(self.optType,'nmpcdclfdcbf')
                [feas, xopt, uopt, Jopt] = solveNMPCDCLFDCBF(self);
            elseif strcmp(self.optType,'mpcdcbf')
                [feas, xopt, uopt, Jopt] = solveMPCDCBF(self);
            elseif strcmp(self.optType,'mpcgcbf')
                [feas, xopt, uopt, Jopt] = solveMPCGCBF(self);
            elseif strcmp(self.optType,'nmpcdcbf')
                [feas, xopt, uopt, Jopt] = solveNMPCDCBF(self);
            else
                disp('optimal control policy undefined');
            end
        end

        function [feas, xopt, uopt, Jopt] = solveDCLFDCBF(self)
            x = sdpvar(self.xdim, 2);
            u = sdpvar(self.udim, 1);
            s = sdpvar(1,1);
            cost = 0;
            constraints = [];
            % initial condition
            constraints = [constraints; x(:,1) == self.xcurr];
            % dynamics
            constraints = [constraints; x(:,2) == self.systemParam.A * x(:,1) + self.systemParam.B * u(:,1)];
            % DCLF constraint
            v = x(:,1)' * self.optParam.P * x(:,1);
            vnext = x(:,2)' * self.optParam.P * x(:,2);
            constraints = [constraints; vnext - v + self.optParam.alpha * v <= s];
            % DCBF constraint
            r = 1;
            b = x(:,1)' * x(:,1) - r^2;
            bnext = x(:,2)' * x(:,2) - r^2;
            constraints = [constraints; bnext - b + self.optParam.gamma * b >= 0];
            % input constraint
            constraints = [constraints; self.systemParam.ul <= u <= self.systemParam.uu];
            % cost
            cost = cost + u' * self.optParam.uWeight * u + s' * self.optParam.sWeight * s;
            opt_settings = sdpsettings('solver','ipopt','verbose',0);
            diagnostics = optimize(constraints, cost, opt_settings);
            if diagnostics.problem == 0
                feas = true;
                xopt = value(x);
                uopt = value(u);
                Jopt = value(cost);
            else
                feas = false;
                xopt= [];
                uopt = [];
                Jopt = value(cost);
            end
            self.solvertime = diagnostics.solvertime;
            fprintf('solver time: %f\n', diagnostics.solvertime);
        end

        function [feas, xopt, uopt, Jopt] = solveNMPCDCLFDCBF(self)
            x = sdpvar(self.xdim, self.optParam.horizon + 1);
            u = sdpvar(self.udim, self.optParam.horizon);
            s = sdpvar(1,self.optParam.horizonCLF);
            omega = sdpvar(1,self.optParam.horizonCBF);
            cost = 0;
            constraints = [];
            % initial condition
            constraints = [constraints; x(:,1) == self.xcurr];
            % MPC cost and constraints
            for i = 1:self.optParam.horizon
                constraints = [constraints;
                    self.systemParam.ul <= u(:,i) <= self.systemParam.uu;
                    x(:,i+1) == self.systemParam.A * x(:,i) + self.systemParam.B * u(:,i)];
                cost = cost + x(:,i)'*self.optParam.xWeight*x(:,i);
                cost = cost + u(:,i)'*self.optParam.uWeight*u(:,i);
            end
            cost = cost + x(:,self.optParam.horizon)' * self.optParam.P * x(:,self.optParam.horizon);
            % CLF constraints
            for i = 1:self.optParam.horizonCLF
                v = x(:,1)' * self.optParam.P * x(:,1);
                vnext = x(:,2)' * self.optParam.P * x(:,2);
                constraints = [constraints; vnext - v + self.optParam.alpha * v <= s(i)];
                cost = cost + s(i)' * self.optParam.sWeight * s(i);
            end
            % CBF constraints
            for i = 1:self.optParam.horizonCBF
                r = 1;
                b = x(:,1)' * x(:,1) - r^2;
                bnext = x(:,2)' * x(:,2) - r^2;
                constraints = [constraints; bnext >= omega(i) * (1 - self.optParam.gamma) * b];
                constraints = [constraints; omega(i) >= 0];
                assign(omega(i), 1);
                cost = cost + self.optParam.omegaWeight * (omega(i) - 1)^2;
            end
            opt_settings = sdpsettings('solver','ipopt','verbose',0);
            diagnostics = optimize(constraints, cost, opt_settings);
            if diagnostics.problem == 0
                feas = true;
                xopt = value(x);
                uopt = value(u);
                Jopt = value(cost);
            else
                feas = false;
                xopt= [];
                uopt = [];
                Jopt = value(cost);
            end
            self.solvertime = diagnostics.solvertime;
            fprintf('solver time: %f\n', diagnostics.solvertime);
        end

        function [feas, xopt, uopt, Jopt] = solveMPCDCBF(self)
            x = sdpvar(self.xdim, self.optParam.horizon + 1);
            u = sdpvar(self.udim, self.optParam.horizon);
            cost = 0;
            constraints = [];
            % initial condition
            constraints = [constraints; x(:,1) == self.xcurr];
            % MPC cost and constraints
            for i = 1:self.optParam.horizon
                constraints = [constraints;
                    self.systemParam.ul <= u(:,i) <= self.systemParam.uu;
                    x(:,i+1) == self.systemParam.A * x(:,i) + self.systemParam.B * u(:,i)];
                cost = cost + x(:,i)'*self.optParam.xWeight*x(:,i);
                cost = cost + u(:,i)'*self.optParam.uWeight*u(:,i);
            end
            cost = cost + x(:,self.optParam.horizon)' * self.optParam.P * x(:,self.optParam.horizon);
            % CBF constraints
            for i = 1:self.optParam.horizon
                r = 0;
                b = r - x(1,i);
                bnext = r - x(1,i+1);
                constraints = [constraints; bnext - b + self.optParam.gamma * b >= 0];
            end
            opt_settings = sdpsettings('solver','ipopt','verbose',0);
            diagnostics = optimize(constraints, cost, opt_settings);
            if diagnostics.problem == 0
                feas = true;
                xopt = value(x);
                uopt = value(u);
                Jopt = value(cost);
            else
                feas = false;
                xopt= [];
                uopt = [];
                Jopt = value(cost);
            end
            self.solvertime = diagnostics.solvertime;
            fprintf('solver time: %f\n', diagnostics.solvertime);
        end

        function [feas, xopt, uopt, Jopt] = solveMPCGCBF(self)
            x = sdpvar(self.xdim, self.optParam.horizon + 1);
            u = sdpvar(self.udim, self.optParam.horizon);
            cost = 0;
            constraints = [];
            % initial condition
            constraints = [constraints; x(:,1) == self.xcurr];
            % MPC cost and constraints
            for i = 1:self.optParam.horizon
                constraints = [constraints;
                    self.systemParam.ul <= u(:,i) <= self.systemParam.uu;
                    x(:,i+1) == self.systemParam.A * x(:,i) + self.systemParam.B * u(:,i)];
                cost = cost + x(:,i)'*self.optParam.xWeight*x(:,i);
                cost = cost + u(:,i)'*self.optParam.uWeight*u(:,i);
            end
            cost = cost + x(:,self.optParam.horizon)' * self.optParam.P * x(:,self.optParam.horizon);
            % GCBF constraints
            num_HO = 3;
            r = 0;
            b = r - x(1,1);
            bnext = r - x(1,1+num_HO);
            constraints = [constraints; bnext >= (1 - self.optParam.gamma)^num_HO * b];
            opt_settings = sdpsettings('solver','ipopt','verbose',0);
            diagnostics = optimize(constraints, cost, opt_settings);
            if diagnostics.problem == 0
                feas = true;
                xopt = value(x);
                uopt = value(u);
                Jopt = value(cost);
            else
                feas = false;
                xopt= [];
                uopt = [];
                Jopt = value(cost);
            end
            self.solvertime = diagnostics.solvertime;
            fprintf('solver time: %f\n', diagnostics.solvertime);
        end

        function [feas, xopt, uopt, Jopt] = solveNMPCDCBF(self)
            x = sdpvar(self.xdim, self.optParam.horizon + 1);
            u = sdpvar(self.udim, self.optParam.horizon);
            omega = sdpvar(1,self.optParam.horizonCBF);
            cost = 0;
            constraints = [];
            % initial condition
            constraints = [constraints; x(:,1) == self.xcurr];
            % MPC cost and constraints
            for i = 1:self.optParam.horizon
                constraints = [constraints;
                    self.systemParam.ul <= u(:,i) <= self.systemParam.uu;
                    x(:,i+1) == self.systemParam.A * x(:,i) + self.systemParam.B * u(:,i)];
                cost = cost + x(:,i)'*self.optParam.xWeight*x(:,i);
                cost = cost + u(:,i)'*self.optParam.uWeight*u(:,i);
            end
            cost = cost + x(:,self.optParam.horizon)' * self.optParam.P * x(:,self.optParam.horizon);
            % CBF constraints
            for i = 1:self.optParam.horizonCBF
                r = 0;
                b = r - x(1,i);
                bnext = r - x(1,i+1);
                constraints = [constraints; bnext >= omega(i) * (1 - self.optParam.gamma) * b];
                constraints = [constraints; omega(i) >= 0];
                assign(omega(i), 1);
                cost = cost + self.optParam.omegaWeight * (omega(i) - 1)^2;
            end
            opt_settings = sdpsettings('solver','ipopt','verbose',0);
            diagnostics = optimize(constraints, cost, opt_settings);
            if diagnostics.problem == 0
                feas = true;
                xopt = value(x);
                uopt = value(u);
                Jopt = value(cost);
            else
                feas = false;
                xopt= [];
                uopt = [];
                Jopt = value(cost);
            end
            self.solvertime = diagnostics.solvertime;
            fprintf('solver time: %f\n', diagnostics.solvertime);
        end
    end
end

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
function circle_test ( )

%*****************************************************************************80
%
%% CIRCLE_TEST applies continuation to the circle problem.
%
%  Discussion:
%
%    Our function is
%
%      f(x,y) = x^2 + y^2 - 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2014
%
%  Author:
%
%    John Burkardt
%
  step_max = 35;
  step_num = 0;
  xy = zeros ( 2, step_max );
%
%  A:
%  Choose a starting point,
%  Choose a starting "parameter",
%  Call Newton to try to get starting point to satisfy implicit function
%  while holding parameter fixed.
%
  n = 2;
  x0 = [ 0.5; -2.0 ];
  p0 = 1;
  tol = 1.0E-05;

  fx_norm = max ( abs ( f_circle ( n, x0 ) ) );
  fprintf ( 1, '  %2d  %14.6g  %14.6g  %8.2e\n', 0, x0(1), x0(2), fx_norm );

  [ status, x2 ] = newton ( n, x0, p0, @f_circle, @fp_circle, tol );

  step_num = step_num + 1;

  xy(1:2,step_num) = x2(1:2,1);
  fx_norm = max ( abs ( f_circle ( n, x2 ) ) );
  fprintf ( 1, '  %2d  %14.6g  %14.6g  %8.2e\n', step_num, x2(1), x2(2), fx_norm );
%
%  B:
%  X0 <= X2.
%  Take a step from the current point X0.
%
  x0 = x2;
  t0 = zeros ( n, 1 );
  h = 0.15;

  for step_num = 2 : step_max

    [ status, x2, t2, p2 ] = step ( n, x0, t0, p0, @f_circle, @fp_circle, h, tol );

    if ( status ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  STEP failed!\n' );
      break
    end

    xy(1:2,step_num) = x2(1:2,1);
    fx_norm = max ( abs ( f_circle ( n, x2 ) ) );
    fprintf ( 1, '  %2d  %14.6g  %14.6g  %8.2e\n', step_num, x2(1), x2(2), fx_norm );

    if ( p0 ~= p2 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Switching parameters from %d to %d\n', p0, p2 );
      p0 = p2;
    end

    x0 = x2;
    t0 = t2;

  end
%
%  Plot the points.
%
  plot ( xy(1,:), xy(2,:), 'r.-', 'Markersize', 25 );
  grid on
  axis equal
  xlabel ( '<--- X --->', 'Fontsize', 16 );
  ylabel ( '<--- Y --->', 'Fontsize', 16 );
  title ( 'Points on the circle, by the continuation method.', 'Fontsize', 24 );

  filename = 'circle_test.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot file saved as "%s"\n', filename );

  return
end


text:
classdef MimConnectionCacheMap < CoreBaseClass

    properties
        ConnectionCacheDictionary
    end

    methods
        function obj = MimConnectionCacheMap()
            obj.ConnectionCacheDictionary = containers.Map('KeyType', 'int32', 'ValueType', 'any');
        end

        function clear(obj)
            for connection = obj.ConnectionCacheDictionary.values
                connection{1}.clear();
            end
        end

        function addConnection(obj, connection)
            if ~obj.ConnectionCacheDictionary.isKey(connection.HashCode)
                obj.ConnectionCacheDictionary(connection.HashCode) = MimModelCache();
            end
        end

        function deleteConnection(obj, connection)
            if obj.ConnectionCacheDictionary.isKey(connection.HashCode)
                obj.ConnectionCacheDictionary.remove(connection.HashCode);
            end
        end

        function connections = getAllConnections(obj)
            connections = obj.ConnectionCacheDictionary.values;
        end

        function connectionCache = getConnection(obj, connection)
            if ~obj.ConnectionCacheDictionary.isKey(connection.HashCode)
                obj.ConnectionCacheDictionary(connection.HashCode) = MimConnectionCache();
            end
            connectionCache = obj.ConnectionCacheDictionary(connection.HashCode);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%% Machine Learning Online Class
%  Exercise 5 | Regularized Linear Regression and Bias-Variance
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     linearRegCostFunction.m
%     learningCurve.m
%     validationCurve.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% =========== Part 1: Loading and Visualizing Data =============
%  We start the exercise by first loading and visualizing the dataset.
%  The following code will load the dataset into your environment and plot
%  the data.
%

% Load Training Data
fprintf('Loading and Visualizing Data ...\n')

% Load from ex5data1:
% You will have X, y, Xval, yval, Xtest, ytest in your environment
load ('ex5data1.mat');

% m = Number of examples
m = size(X, 1);

% Plot training data
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');

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

%% =========== Part 2: Regularized Linear Regression Cost =============
%  You should now implement the cost function for regularized linear
%  regression.
%

theta = [1 ; 1];
J = linearRegCostFunction([ones(m, 1) X], y, theta, 1);

fprintf(['Cost at theta = [1 ; 1]: %f '...
         '\n(this value should be about 303.993192)\n'], J);

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

%% =========== Part 3: Regularized Linear Regression Gradient =============
%  You should now implement the gradient for regularized linear
%  regression.
%

theta = [1 ; 1];
[J, grad] = linearRegCostFunction([ones(m, 1) X], y, theta, 1);

fprintf(['Gradient at theta = [1 ; 1]:  [%f; %f] '...
         '\n(this value should be about [-15.303016; 598.250744])\n'], ...
         grad(1), grad(2));

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

%% =========== Part 4: Train Linear Regression =============
%  Once you have implemented the cost and gradient correctly, the
%  trainLinearReg function will use your cost function to train
%  regularized linear regression.
%
%  Write Up Note: The data is non-linear, so this will not give a great
%                 fit.
%

%  Train linear regression with lambda = 0
lambda = 0;
[theta] = trainLinearReg([ones(m, 1) X], y, lambda);

%  Plot fit over the data
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');
hold on;
plot(X, [ones(m, 1) X]*theta, '--', 'LineWidth', 2)
hold off;

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

%% =========== Part 5: Learning Curve for Linear Regression =============
%  Next, you should implement the learningCurve function.
%
%  Write Up Note: Since the model is underfitting the data, we expect to
%                 see a graph with "high bias" -- slide 8 in ML-advice.pdf
%

lambda = 0;
[error_train, error_val] = ...
    learningCurve([ones(m, 1) X], y, ...
                  [ones(size(Xval, 1), 1) Xval], yval, ...
                  lambda);

plot(1:m, error_train, 1:m, error_val);
title('Learning curve for linear regression')
legend('Train', 'Cross Validation')
xlabel('Number of training examples')
ylabel('Error')
axis([0 13 0 150])

fprintf('# Training Examples\tTrain Error\tCross Validation Error\n');
for i = 1:m
    fprintf('  \t%d\t\t%f\t%f\n', i, error_train(i), error_val(i));
end

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

%% =========== Part 6: Feature Mapping for Polynomial Regression =============
%  One solution to this is to use polynomial regression. You should now
%  complete polyFeatures to map each example into its powers
%

p = 8;

% Map X onto Polynomial Features and Normalize
X_poly = polyFeatures(X, p);
[X_poly, mu, sigma] = featureNormalize(X_poly);  % Normalize
X_poly = [ones(m, 1), X_poly];                   % Add Ones

% Map X_poly_test and normalize (using mu and sigma)
X_poly_test = polyFeatures(Xtest, p);
X_poly_test = bsxfun(@minus, X_poly_test, mu);
X_poly_test = bsxfun(@rdivide, X_poly_test, sigma);
X_poly_test = [ones(size(X_poly_test, 1), 1), X_poly_test];         % Add Ones

% Map X_poly_val and normalize (using mu and sigma)
X_poly_val = polyFeatures(Xval, p);
X_poly_val = bsxfun(@minus, X_poly_val, mu);
X_poly_val = bsxfun(@rdivide, X_poly_val, sigma);
X_poly_val = [ones(size(X_poly_val, 1), 1), X_poly_val];           % Add Ones

fprintf('Normalized Training Example 1:\n');
fprintf('  %f  \n', X_poly(1, :));

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

%% =========== Part 7: Learning Curve for Polynomial Regression =============
%  Now, you will get to experiment with polynomial regression with multiple
%  values of lambda. The code below runs polynomial regression with
%  lambda = 0. You should try running the code with different values of
%  lambda to see how the fit and learning curve change.
%

lambda = 0;
[theta] = trainLinearReg(X_poly, y, lambda);

% Plot training data and fit
figure(1);
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
plotFit(min(X), max(X), mu, sigma, theta, p);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');
title (sprintf('Polynomial Regression Fit (lambda = %f)', lambda));

figure(2);
[error_train, error_val] = ...
    learningCurve(X_poly, y, X_poly_val, yval, lambda);
plot(1:m, error_train, 1:m, error_val);

title(sprintf('Polynomial Regression Learning Curve (lambda = %f)', lambda));
xlabel('Number of training examples')
ylabel('Error')
axis([0 13 0 100])
legend('Train', 'Cross Validation')

fprintf('Polynomial Regression (lambda = %f)\n\n', lambda);
fprintf('# Training Examples\tTrain Error\tCross Validation Error\n');
for i = 1:m
    fprintf('  \t%d\t\t%f\t%f\n', i, error_train(i), error_val(i));
end

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

%% =========== Part 8: Validation for Selecting Lambda =============
%  You will now implement validationCurve to test various values of
%  lambda on a validation set. You will then use this to select the
%  "best" lambda value.
%

[lambda_vec, error_train, error_val] = ...
    validationCurve(X_poly, y, X_poly_val, yval);

close all;
plot(lambda_vec, error_train, lambda_vec, error_val);
legend('Train', 'Cross Validation');
xlabel('lambda');
ylabel('Error');

fprintf('lambda\t\tTrain Error\tValidation Error\n');
for i = 1:length(lambda_vec)
  fprintf(' %f\t%f\t%f\n', ...
            lambda_vec(i), error_train(i), error_val(i));
end

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

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function a = conex1_inverse ( alpha )

%*****************************************************************************80
%
%% CONEX1_INVERSE returns the inverse of the CONEX1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the scalar defining A.
%    A common value is 100.0.
%
%    Output, real A(4,4), the matrix.
%
  a = zeros ( 4, 4 );

  a(1,1) =  1.0;
  a(1,2) =  1.0 - alpha;
  a(1,3) =        alpha;
  a(1,4) =  2.0;

  a(2,1) =  0.0;
  a(2,2) =  1.0 + alpha;
  a(2,3) =      - alpha;
  a(2,4) =  0.0;

  a(3,1) =  0.0;
  a(3,2) = -1.0;
  a(3,3) =  1.0;
  a(3,4) =  1.0 / alpha;

  a(4,1) = 0.0;
  a(4,2) = 0.0;
  a(4,3) = 0.0;
  a(4,4) = 1.0 / alpha;

  return
end


text:
function pass = test_diff( )

tol = 2e2*chebfunpref().cheb2Prefs.chebfun2eps;

% Simple tests:
f = spherefun(@(lam,th) cos(lam).*sin(th));  % x
fx = diff(f, 1);
exact = @(lam,th) 1-(cos(lam).*sin(th)).^2;  % 1-x^2 or y^2 + z^2
pass(1) = SampleError(exact, fx) < tol;
fy = diff(f, 2);
exact = @(lam,th) -cos(lam).*sin(th).*sin(lam).*sin(th);  % -x*y
pass(2) = SampleError(exact, fy) < tol;
fz = diff(f, 3);
exact = @(lam,th) -cos(lam).*sin(th).*cos(th);  % -x*z
pass(3) = SampleError(exact, fz) < tol;

f = spherefun(@(lam,th) cos(th));  % z
fx = diff(f, 1);
exact = @(lam,th) -cos(lam).*sin(th).*cos(th);  % -x*z
pass(4) = SampleError(exact, fx) < tol;
fy = diff(f, 2);
exact = @(lam,th) -sin(lam).*sin(th).*cos(th);  % -y*z
pass(5) = SampleError(exact, fy) < tol;
fz = diff(f, 3);
exact = @(lam,th) 1-cos(th).^2;  % 1-z.^2  or x^2 + y^2
pass(6) = SampleError(exact, fz) < tol;

a = 1/sqrt(2);
f = spherefun(@(lam,th) sin(lam-a).*sin(th));  % Shifted y
fx = diff(f, 1);
exact = @(lam,th) 0.25*(2*cos(lam).*cos(2*th).*sin(lam-a) - sin(2*lam-a) ...
    - 3*sin(a));
pass(7) = SampleError(exact, fx) < tol;
fy = diff(f, 2);
exact = @(lam,th) 0.25*(2*sin(lam).*cos(2*th).*sin(lam-a) + cos(2*lam-a) ...
    + 3*cos(a));
pass(8) = SampleError(exact, fy) < tol;
fz = diff(f, 3);
exact = @(lam,th) -0.5*sin(lam-a).*sin(2*th);
pass(9) = SampleError(exact, fz) < tol;

% Gaussian
lam0 = pi/sqrt(2);
th0 = (sqrt(5) - 2)/2;
sig2 = 1;
r2 = @(lam,th) 2*(1 - (sin(th)*sin(th0)).*cos(lam-lam0)-cos(th)*cos(th0));
f = @(lam,th) exp(-sig2 * r2(lam, th));
fx = diff(spherefun(f), 1);
exact = @(lam,th) -0.5*sig2*exp(-sig2*r2(lam,th)).*(2*cos(lam)...
    .*cos(th0).*sin(2*th)+(cos(2*lam-lam0)-3*cos(lam0)-...
    2*cos(lam).*cos(lam-lam0).*cos(2*th))*sin(th0));
pass(10) = SampleError(exact, fx) < tol;

end

function sample_error = SampleError(h, g)
m = 6;
n = m;
[x, y] = getPoints(m, n);
[L2, T2] = meshgrid(x, y);
F = h(L2, T2);
approx = fevalm(g, x, y);
sample_error = norm(F(:) - approx(:), inf);
end

function [x, y] = getPoints(m, n)

x = trigpts(2*n,[-pi pi]);
y = linspace(0,pi,m).';

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% LCM - Variational Bayesian inference for latent component models.

% Likelihood: p(Y|A,S,C) = N(Y|A*S,C)
%
% Priors p(A) and p(S) can be arbitrary, because the user must specify
% the posterior estimation function.
%
% Covariance could be a Kronecker product of two arbitrary (or diagonal?)
% covariance matrices C_a and C_s? Missing values cause a nasty problem
% to Kronecker-covariance.. :(
%
% What if also an arbitrary diagonal covariance matrix C_d? The full
% joint covariance matrix would be:
%
% C = sqrt(C_d) * kron(C_a,C_s) * sqrt(C_d)
%
% Then missing values could be considered by using large (infinite)
% values in C_d..?

% Last modified 2010-10-12
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function Q = lcm(Y, D, varargin)

[M,N] = size(Y);

options = struct( ...
    'init',         [],               ...
    'prior',        [],               ...
    'common_tau',   true,             ...
    'common_nu',    false,            ...
    'rotate',       1,                ...
    'update_alpha', 1,               ...
    'update_beta',  1,               ...
    'update_nu',    1,               ...
    'autosavetime', nan,              ...
    'autosavefile', 'vbrfa_autosave', ...
    'maxiter',      100,              ...
    'robustness',   'independent-t', ... % none / multivariate-t / independent-t
    'user_data',    []);

% Parse arguments
[options, errmsg] = argparse( options, varargin{:} );
error(errmsg);

% Missing values
Obs = ~isnan(Y);
[ObsM, ObsN] = find(Obs);

% Helpful variables
Id = speye(D);
log2pi = log(2*pi);

logP = -inf;
loglikelihood = nan*zeros(options.maxiter,1);
start_time = cputime;
time = nan*zeros(options.maxiter,1);

lastsave = now;

for iter=1:options.maxiter

  t = cputime;

  %%%%%%%%%%%%%%%%%%%%%%
  %% Update parameters

  %%%%%%%%%%%%%%%%%
  %
  % X BLOCK
  %
  %%%%%%%%%%%%%%%%%

  %% UPDATE X

  % q(X)
  [X, XX] =
  for n=1:N
    Mmv = Obs(:,n);
% $$$     UTau = U(Mmv,n).*tau(Mmv);
% $$$     UWWTau = sum(bsxfun(@times, ...
% $$$                         reshape(U(Mmv,n),[1,1,sum(Mmv)]), ...
% $$$                         WtauW(:,:,Mmv)), ...
% $$$                  3);

    % Distribution
    UWtauW = wsum(WtauW(:,:,Mmv),U(Mmv,n),3);
    CovX(:,:,n) = inv( speye(D) + UWtauW );
    X(:,n) = CovX(:,:,n) * ...
             ( tauW(Mmv,:)'*spdiag(Y(Mmv,n)) - WtauMu(Mmv,:)' ) * U(Mmv,n);

    % Expectations
    XX(:,:,n) = CovX(:,:,n) + X(:,n)*X(:,n)';
  end

  % KL-term: <log q(X)> - <log p(X)>
  KL_X = 0;
  for n=1:N
    KL_X = KL_X - gaussian_entropy(logdet_cov(CovX(:,:,n)), D);
    KL_X = KL_X - gaussian_logpdf(trace(XX(:,:,n)), ...
                                  0, ...
                                  0, ...
                                  0, ...
                                  D);
  end

  %%%%%%%%%%%%%%%%%
  %
  % W BLOCK
  %
  %%%%%%%%%%%%%%%%%

  %
  % Update distributions and evaluate statistics
  %

  % Update q(W,Mu|tau)
  for m=1:M
    Nmv = Obs(m,:);

    UXX = wsum(XX(:,:,Nmv),U(m,Nmv),3);
    UX = X(:,Nmv)*U(m,Nmv)';
    UY = U(m,Nmv).*Y(m,Nmv);
    UYX = [UY*X(:,Nmv)', sum(UY)];

    % Distribution
    CovWMu(:,:,m) = inv(spdiag([alpha; beta]) + [UXX, UX; UX', sum(U(m,Nmv))]);
    W(m,:) = UYX*CovWMu(:,1:D,m);
    Mu(m) = UYX*CovWMu(:,end,m);
  end

  % Update q(tau)
  for m=1:M
    Nmv = Obs(m,:);
    a_tau(m) = p.a_tau + 0.5 * sum(Nmv);
    b_tau(m) = p.b_tau ...
        + 0.5 * (U(m,Nmv).*Y(m,Nmv))*Y(m,Nmv)' ...
        - 0.5 * [W(m,:),Mu(m)] * inv(CovWMu(:,:,m)) * [W(m,:),Mu(m)]';
  end
  if options.common_tau
    % Common tau
    a_tau(:) = p.a_tau + sum(a_tau-p.a_tau);
    b_tau(:) = p.b_tau + sum(b_tau-p.b_tau);
  end
  tau = a_tau./b_tau;
  logtau = psi(a_tau) - log(b_tau);
  for m=1:M
    % Expectations
    tauW(m,:) = tau(m) * W(m,:);
    tauMu(m) = tau(m) * Mu(m);
    WtauW(:,:,m) = tau(m)*W(m,:)'*W(m,:) + CovWMu(1:D,1:D,m);
    MutauMu(m) = tau(m)*Mu(m)^2 + CovWMu(end,end,m);
    WtauMu(m,:) = tau(m)*Mu(m)*W(m,:) + CovWMu(end,1:D,m);
  end

  % Update q(alpha)
  if (numel(options.update_alpha)==1 && iter >= options.update_alpha) || ...
        any(iter==options.update_alpha)
    a_alpha(:) = p.a_alpha + 0.5*M;
    b_alpha(:) = p.b_alpha + 0.5*diag(sum(WtauW, 3));
  end
  alpha = a_alpha ./ b_alpha;
  logalpha = psi(a_alpha) - log(b_alpha);

  % Update q(beta)
  if (numel(options.update_beta)==1 && iter >= options.update_beta) || ...
        any(iter==options.update_beta)
    a_beta = p.a_beta + 0.5*M;
    b_beta = p.b_beta + 0.5*sum(MutauMu);
  end
  beta = a_beta ./ b_beta;
  logbeta = psi(a_beta) - log(b_beta);

  %
  % Compute KL terms
  %

  % KL-term: <log q(W,Mu|tau)> - <log p(W|tau,alpha)> - <log p(Mu|tau,beta)>
  KL_WMu = 0;
  for m=1:M
    KL_WMu = KL_WMu - gaussian_entropy(-(D+1)*logtau(m) ...
                                       + logdet_cov(CovWMu(:,:,m)), ...
                                       D);
    KL_WMu = KL_WMu - gaussian_logpdf(diag(WtauW(:,:,m))'*alpha, ...
                                      0, ...
                                      0, ...
                                      -D*logtau(m) - sum(logalpha), ...
                                      D);
    KL_WMu = KL_WMu - gaussian_logpdf(diag(MutauMu(m))'*beta, ...
                                      0, ...
                                      0, ...
                                      -logtau(m) - sum(logbeta), ...
                                      1);
  end
% $$$   KL_W = 0;
% $$$   for m=1:M
% $$$     KL_W = KL_W - gaussian_entropy(-D*logtau(m) + logdet_cov(CovW(:,:,m)), D);
% $$$     KL_W = KL_W - gaussian_logpdf(diag(WtauW(:,:,m))'*w, ...
% $$$                                   0, ...
% $$$                                   0, ...
% $$$                                   -D*logtau(m) - sum(logw), ...
% $$$                                   D);
% $$$   end

  % KL-term: <log q(tau)> - <log p(tau)>
  KL_tau = 0;
  if options.common_tau
    KL_tau = KL_tau - gamma_entropy(a_tau(1), b_tau(1));
    KL_tau = KL_tau - gamma_logpdf(tau(1),logtau(1),p.a_tau,p.b_tau);
  else
    KL_tau = KL_tau - sum(gamma_entropy(a_tau, b_tau));
    KL_tau = KL_tau - sum(gamma_logpdf(tau,logtau,p.a_tau,p.b_tau));
  end

  % KL-term: <log q(alpha)> - <log p(alpha)>
  KL_alpha = 0;
  KL_alpha = KL_alpha - sum(gamma_entropy(a_alpha, b_alpha));
  KL_alpha = KL_alpha - sum(gamma_logpdf(alpha, logalpha, p.a_alpha, p.b_alpha));

  % KL-term: <log q(beta)> - <log p(beta)>
  KL_beta = 0;
  KL_beta = KL_beta - sum(gamma_entropy(a_beta, b_beta));
  KL_beta = KL_beta - sum(gamma_logpdf(beta, logbeta, p.a_beta, p.b_beta));

  %%%%%%%%%%%%%%%%%
  %
  % NOISE BLOCK
  %
  %%%%%%%%%%%%%%%%%

  %% UPDATE U AND NU
  switch options.robustness

   case 'none'
    %% Normal PCA / FA
    U(Obs) = 1;
    logU(Obs) = 0;
    nu(:) = inf;
    KL_U = 0;

   case 'multivariate-t'
    %% Joint Student-t for the dimensions (Archambaeu et al), except X
    %% independent of U, i.e., p(X) is Gaussian.

    % Evaluate squared errors
    for m=1:M
      Nmv = Obs(m,:);

      E2(m,Nmv) = tau(m)*Y(m,Nmv).^2 ...
          + mtraceprod(WtauW(:,:,m), XX(:,:,Nmv)) ...
          + MutauMu(m) ...
          - 2*Y(m,Nmv).*(tauW(m,:)*X(:,Nmv)) ...
          - 2*Y(m,Nmv).*tauMu(m) ...
          + 2*WtauMu(m,:)*X(:,Nmv);
    end

    % Update nu
    if (numel(options.update_nu)==1 && iter >= options.update_nu) ...
          || any(iter==options.update_nu)

      % Common nu
      nu(:) = t_ml(sum(E2,1), nu(1), sum(Obs,1));
 %     disp('Update nu');
    end

    % Update U
    a_U(1,:) = (nu(1) + sum(Obs,1)) / 2;
    b_U(1,:) = (nu(1) + sum(E2,1)) / 2;
    a_U(:,:) = repmat(a_U(1,:), [M,1]);
    b_U(:,:) = repmat(b_U(1,:), [M,1]);
    U(Obs) = a_U(Obs)./b_U(Obs);
    logU(Obs) = psi(a_U(Obs)) - log(b_U(Obs));

    % KL-term: <log q(U)> - <log p(U)>
    KL_U = 0;
    Nmv = sum(Obs,1) > 0;
    KL_U = KL_U - sum(gamma_entropy(a_U(1,Nmv),b_U(1,Nmv)));
    KL_U = KL_U - sum(gamma_logpdf(U(1,Nmv),logU(1,Nmv),nu(1)/2,nu(1)/2));

   case 'independent-t'
    %% Independent Student-t for each dimension (Luttinen et al)

    % Evaluate squared errors (this costs A LOT!)
    for m=1:M
      Nmv = Obs(m,:);

      E2(m,Nmv) = tau(m)*Y(m,Nmv).^2 ...
          + mtraceprod(WtauW(:,:,m), XX(:,:,Nmv)) ...
          + MutauMu(m) ...
          - 2*Y(m,Nmv).*(tauW(m,:)*X(:,Nmv)) ...
          - 2*Y(m,Nmv).*tauMu(m) ...
          + 2*WtauMu(m,:)*X(:,Nmv);
    end

    % Update nu
    if (numel(options.update_nu)==1 && iter >= options.update_nu) ...
          || any(iter==options.update_nu)
      if options.common_nu
        % Common nu
        nu(:) = t_ml(E2(Obs), nu(1));
      else
        for m=1:M
          % Separate nu
          if true
            % Type II ML
            nu(m) = t_ml(E2(m,Obs(m,:)), nu(m));
          else
            % EM ML
            nv = Obs(m,:);
            func = @(lognu) (-1-log(exp(lognu)/2)+psi(exp(lognu)/2)- ...
                             mean(logU(m,nv)-U(m,nv)));
            lognu = fzero(func, log(nu(m)));
            nu(m) = exp(lognu);
          end
        end
      end
%      disp('Update nu');
    end

    % Update U
    for m=1:M
      Nmv = Obs(m,:);
      a_U(m,Nmv) = (nu(m)+1) / 2;
      b_U(m,Nmv) = (nu(m) + E2(m,Nmv)) / 2;
% $$$     b_U(m,Nmv) = (nu(m) + E2(m,Nmv)*tau(m)) / 2;
    end
    U(Obs) = a_U(Obs)./b_U(Obs);
    logU(Obs) = psi(a_U(Obs)) - log(b_U(Obs));

    % KL-term: <log q(U)> - <log p(U)>
    KL_U = 0;
    KL_U = KL_U - sum(gamma_entropy(a_U(Obs),b_U(Obs)));
    for m=1:M
      Nmv = Obs(m,:);
      KL_U = KL_U - sum(gamma_logpdf(U(m,Nmv),logU(m,Nmv),nu(m)/2,nu(m)/2));
    end

  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %% Calculate lower bound of the log-likelihood

  monitor_loglikelihood = true;

  if monitor_loglikelihood

    % Cost from Y

    % We can take advantage of the linearity of gaussian_logpdf by summing
    % before the function call:
    tauYh = bsxfun(@plus, tauW*X, tauMu);
    y_invCov_y = (tau(ObsM).*U(Obs))'*(Y(Obs).^2);
    y_invCov_mu = U(Obs)'*(Y(Obs).*tauYh(Obs));
    mu_invCov_mu = 0;
    for m=1:M
      Nmv = Obs(m,:);
      mu_invCov_mu = mu_invCov_mu ...
          + traceprod(WtauW(:,:,m),wsum(XX(:,:,Nmv),U(m,Nmv),3)) ...
          + sum(U(m,Nmv)) * MutauMu(m) ...
          + 2*WtauMu(m,:)*(X(:,Nmv)*U(m,Nmv)');
    end
    logdet_Cov = - sum(logtau(ObsM)) - sum(logU(Obs));
    dim_y = sum(Obs(:));

    logpdf_Y = gaussian_logpdf(y_invCov_y, ...
                               y_invCov_mu, ...
                               mu_invCov_mu, ...
                               logdet_Cov, ...
                               dim_y);

    old_logP = logP;

    %    logpdf_Y

% $$$     KL_U
% $$$     KL_X
% $$$     KL_tau
% $$$     KL_W
% $$$     KL_w
    logP = logpdf_Y - KL_U - KL_X - KL_tau - KL_WMu - KL_alpha - KL_beta;

% $$$     % DEBUG STUFF
% $$$     kl_Y(iter) = logpdf_Y;
% $$$     kl_U(iter) = KL_U;
% $$$     kl_X(iter) = KL_X;
% $$$     kl_W(iter) = KL_W;
% $$$     kl_tau(iter) = KL_tau;
% $$$     kl_w(iter) = KL_w;

    loglikelihood(iter) = logP;
  else
    loglikelihood(iter) = nan;
  end

  time(iter) = cputime - start_time;

  % Debugging: Check that the bound really improves
  if iter > 1 && loglikelihood(iter) < loglikelihood(iter-1)
    logP_diff = (loglikelihood(iter-1)-loglikelihood(iter))/ ...
        abs(loglikelihood(iter-1));
    warmsg = sprintf(['Lower bound decreased %.2e percents. Bug or numerical ' ...
                      'inaccuracy?'], logP_diff);
    warning(warmsg);
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  %% Monitoring stuff

  % Show progress
  fprintf('Step %d: loglikelihood=%e (%.3f seconds)\n', ...
          iter, loglikelihood(iter), cputime-t);

  % Check whether to save the results
  if (now-lastsave)*3600*24 >= options.autosavetime
    fprintf('Saving to %s...', options.autosavefile);
    save(options.autosavefile, 'W','CovW','w','a_w','b_w','X','CovX', ...
         'tau', 'a_tau','b_tau','U','a_U','b_U','nu', 'loglikelihood', ...
         'time', 'iter', 'options');
    lastsave = now;
    fprintf(' done.\n');
  end

  % Transformations for speeding up
  if (numel(options.rotate)==1 && iter>=options.rotate && options.rotate>=1) ...
        || any(iter==options.rotate)
    %orthogonalize();
    rotation2();
  end

end

% $$$ figure
% $$$ kl = [kl_Y(:), kl_U(:), kl_X(:), kl_W(:), kl_tau(:), kl_w(:)];
% $$$ kl = kl(6:end,:);
% $$$ kl = eminus(kl,mean(kl,1));
% $$$ plot(kl)

% Results as a struct
Q.W = W;
Q.CovWMu = CovWMu;
Q.alpha = alpha;
Q.a_alpha = a_alpha;
Q.b_alpha = b_alpha;
Q.Mu = Mu;
Q.beta = beta;
Q.a_beta = beta;
Q.b_beta = b_beta;
Q.X = X;
Q.CovX = CovX;
Q.tau = tau;
Q.a_tau = a_tau;
Q.b_tau = b_tau;
Q.U = U;
Q.a_U = a_U;
Q.b_U = b_U;
Q.nu = nu;
Q.loglikelihood = loglikelihood(1:iter);
Q.time = time(1:iter);
Q.iter = iter;
Q.options = options;

%%%%%%%%%%%%%%%%%%%%%%
%% NESTED FUNCTIONS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function rotation2()

%alpha_before = alpha

% Rotates W and X, and moves bias from X to Mu.

%Y0 = bsxfun(@plus, W*X, Mu);

sumX = sum(X,2);
sumXX = [sum(XX,3), sumX; sumX', N]; % extended model: [X;1]*[X;1]'
sumWtauW = sum(WtauW,3);
sumMutauMu = sum(MutauMu);
sumWtauMu = sum(WtauMu,1);
sumWMutauWMu = [sumWtauW, sumWtauMu'; sumWtauMu, sumMutauMu];

%diag_wtauw = diag(sum(WtauW,3))

%
% Optimize the rotation
%

R0 = eye([D,D+1]);
%mycheckgrad(@cost, R0(:), 1e-6);
r = minimize(R0(:), @cost, 50);
%r(:) = randn(size(r)); % DEBUG !!!
R = [reshape(r,[D,D+1]); zeros(1,D), 1];
%R
%R = eye(D+1);% DEBUG!!!
invR = inv(R);

%
% Apply the rotation
%

% W and Mu
tauWMu = [tauW,tauMu]*R;
tauW = tauWMu(:,1:D);
tauMu = tauWMu(:,end);
for m=1:M
  Z = [WtauW(:,:,m), WtauMu(m,:)'; WtauMu(m,:), MutauMu(m)];
  % Expectations
  WtauW(:,:,m) = R(:,1:D)'*Z*R(:,1:D);
  MutauMu(m) = R(:,end)'*Z*R(:,end);
  WtauMu(m,:) = R(:,end)'*Z*R(:,1:D);
end

% X
X = invR(1:D,:)*[X;ones(1,N)];
for n=1:N
  XX(:,:,n) = invR(1:D,:)*[XX(:,:,n),X(:,n);X(:,n)',1]*invR(1:D,:)';
end

% alpha and beta
%diag_wtauw = diag(sum(WtauW,3))
b_alpha = p.b_alpha + 0.5*diag(sum(WtauW,3));
alpha = a_alpha ./ b_alpha;
b_beta = p.b_beta + 0.5*sum(MutauMu);
beta = a_beta ./ b_beta;

%alpha_after = alpha

%meanX = mean(X,2)
%sumXX = sum(XX,3) / N
%sumWW = sum(WtauW,3) / M

%Y1 = bsxfun(@plus, W*X, Mu);

%difference = sqrt(mean((Y1(:)-Y0(:)).^2))

%error('jou')

  function [f,df] = cost(r)

  R = [reshape(r,[D,D+1]); zeros(1,D), 1];
  R(end) = 1;

  invR = inv(R);
  invRR = inv(R*R');
  logdetR = logdet(R);

  % Transform W and Mu
  WW_R = sumWMutauWMu*R;
  R_WW_R = R'*WW_R;

  % Transform q(alpha) and q(beta)
  % (concatenate alpha and beta to one vector for simplicity)
  z = diag(R_WW_R);
  r_b_alpha = [repmat(p.b_alpha,D,1);p.b_beta] + 0.5*z;
  r_alpha = [a_alpha;a_beta] ./ r_b_alpha;
  r_logalpha = - log(r_b_alpha); % up to a constant
  grad_b = WW_R;
  grad_logb = grad_b * diag(1./r_b_alpha);
  grad_logalpha = -grad_logb;
  grad_alpha = grad_logalpha * diag(r_alpha);

  % Transform X

  % Cost from W and Mu
  %
  % TODO: There is something wrong with this!!
  logdet_Cov = 2*M*logdetR;
  x_invCov0_x = traceprod(spdiag(r_alpha),R_WW_R);
  x_invCov0_mu = 0;
  mu_invCov0_mu = 0;
  logdet_Cov0 = -M * sum(r_logalpha);
  grad_logdet_Cov = 2*M*invR';
  grad_x_invCov0_x = 2*WW_R*spdiag(r_alpha) ...
      + grad_alpha * diag(diag(R_WW_R));
  grad_x_invCov0_mu = 0;
  grad_mu_invCov0_mu = 0;
  grad_logdet_Cov0 = -M * grad_logalpha;
  [KL_W,dKL_W] = vb_rotationcost_gaussian(logdet_Cov, ...
                                          x_invCov0_x, ...
                                          x_invCov0_mu, ...
                                          mu_invCov0_mu, ...
                                          logdet_Cov0, ...
                                          grad_logdet_Cov, ...
                                          grad_x_invCov0_x, ...
                                          grad_x_invCov0_mu, ...
                                          grad_mu_invCov0_mu, ...
                                          grad_logdet_Cov0);
  %KL_W = 0;
  %dKL_W = 0;

  % Cost from alpha and beta
  prior_a = [repmat(p.a_alpha,D,1);p.a_beta];
  prior_b = [repmat(p.b_alpha,D,1);p.b_beta];
  post_a = [a_alpha;a_beta];
  A0_logB = prior_a'*log(r_b_alpha);
  A_B0_invB = prior_b'*r_alpha;
  grad_A0_logB = grad_logb * diag(prior_a);
  grad_A_B0_invB = grad_logalpha * diag(prior_b);
  [KL_alpha,dKL_alpha] = vb_rotationcost_gamma(A0_logB,      ...
                                               A_B0_invB,    ...
                                               grad_A0_logB, ...
                                               grad_A_B0_invB);
  %KL_alpha = 0;
  %dKL_alpha = 0;

  % Cost from X
  logdet_Cov = -2*N*logdetR;
  x_invCov0_x = traceprod(invR(1:D,:)'*invR(1:D,:),sumXX);
  x_invCov0_mu = 0;
  mu_invCov0_mu = 0;
  logdet_Cov0 = 0;
  grad_logdet_Cov = -2*N*invR';
  grad_x_invCov0_x = -2*invRR*sumXX*invR';
  grad_x_invCov0_mu = 0;
  grad_mu_invCov0_mu = 0;
  grad_logdet_Cov0 = 0;
  [KL_X,dKL_X] = vb_rotationcost_gaussian(logdet_Cov, ...
                                          x_invCov0_x, ...
                                          x_invCov0_mu, ...
                                          mu_invCov0_mu, ...
                                          logdet_Cov0, ...
                                          grad_logdet_Cov, ...
                                          grad_x_invCov0_x, ...
                                          grad_x_invCov0_mu, ...
                                          grad_mu_invCov0_mu, ...
                                          grad_logdet_Cov0);
  %KL_X = 0;
  %dKL_X = 0;

  % Vectorize results
  f = KL_W + KL_alpha + KL_X;
  df = dKL_W + dKL_alpha + dKL_X;
  df = df(1:D,:);
  df = df(:);

  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function orthogonalize

  %% Move bias (approximate: mean of X to zero)
% $$$   dmu = mean(X,2);
% $$$   X = X - repmat(dmu,1,N);
% $$$   tauMu = tauMu + tauW*dmu;
% $$$   WtauMu = WtauMu ???
% $$$   MutauMu = ???

  %% Rotate W and X

  Yh0 = tauW*X;

  % Whitening of X
  [Vx,D2x] = svd(sum(XX,3)/N);
  Dx = spdiag(sqrt(diag(D2x))); % convert to sparse diagonal form
  Qx = Dx \ Vx';
  Qw = full(Vx * Dx); % funny: could be sparse if size(Dx)=[1 1] ??
  tauW = tauW * Qw;
  WtauMu = WtauMu * Qw;
  for m=1:M
%    CovW(:,:,m) = Qw' * CovW(:,:,m) * Qw;
    WtauW(:,:,m) = Qw' * WtauW(:,:,m) * Qw;
% $$$     WW(:,:,m) = W(m,:)'*W(m,:) + CovW(:,:,m);
  end

  % Orthogonalization of W
  [Vw,Dw] = svd(sum(WtauW,3)/M);
% $$$   [Vw,Dw] = svd(sum(WW,3)/M);
  Qx = Vw' * Qx;
  Qw = Vw;
%  W = W * Qw;
  tauW = tauW * Qw;
  WtauMu = WtauMu * Qw;
  for m=1:M
%    CovW(:,:,m) = Qw' * CovW(:,:,m) * Qw;
    WtauW(:,:,m) = Qw' * WtauW(:,:,m) * Qw; %W(m,:)'*W(m,:) + CovW(:,:,m);
% $$$     WW(:,:,m) = W(m,:)'*W(m,:) + CovW(:,:,m);
  end

  % Apply rotations to X
  X = Qx * X;
  for n = 1:N
    CovX(:,:,n) = Qx * CovX(:,:,n) * Qx';
    XX(:,:,n) = X(:,n)*X(:,n)' + CovX(:,:,n);
  end

  Yh1 = tauW*X;

% $$$   diff_in_orth = norm(Yh1-Yh0)
% $$$
% $$$   sumWtauW = sum(WtauW,3)
% $$$   sumXX = sum(XX,3) / N

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%

function initialize()

  disp('Initialize variables..')

  % Initialize X
  X = nan * zeros(D,N);
  CovX = nan*zeros([D,D,N]); % repmat(eye(D),[1,1,N]);

  % Initialize q(alpha)
  a_alpha = nan * ones(D,1);
  b_alpha = nan * ones(D,1);
  alpha = 1e-3 * ones(D,1);
  logalpha = nan * ones(D,1);

  % Initialize q(beta)
  a_beta = nan;
  b_beta = nan;
  beta = 1e-3;
  logbeta = nan;

  % Initialize W and Mu
  W = randn(M,D);
  Mu = zeros(M,1);
  CovWMu = repmat(eye(D+1),[1,1,M]);

  % Initialize tau
  a_tau = nan * ones(M,1);
  b_tau = nan * ones(M,1);
  tau = 1e1 * ones(M,1);
  logtau = nan * ones(M,1);

  % Initialize nu and q(U)
  nu = 1 * ones(M,1);
  a_U = nan * ones(M,N);
  b_U = nan * ones(M,N);
  U = ones(M,N);
  logU = nan * zeros(M,N);

  % Initialize log likelihood
  loglikelihood = -Inf;

  % Initialize mean squared errors < (y-wx-m)^2 >
  E2 = zeros(size(Y));

  % Use given initialization (size matching is checked by using ':' )
  if isstruct(options.init)
    % X
    if isfield(options.init, 'X')
      disp('Custom initialization for X.');
      warning(['Currently, X is updated first, so given initialization has ' ...
               'no effect.']);
      X(:,:) = options.init.X;
    end
    if isfield(options.init, 'CovX')
      disp('Custom initialization for CovX.');
      warning(['Currently, X is updated first, so given initialization has ' ...
               'no effect.']);
      CovX(:,:,:) = options.init.CovX;
    end
    % W and Mu
    if isfield(options.init, 'W')
      disp('Custom initialization for W.');
      W(:,:) = options.init.W;
    end
    if isfield(options.init, 'Mu')
      disp('Custom initialization for Mu.');
      Mu(:,:) = options.init.Mu;
    end
    if isfield(options.init, 'CovWMu')
      disp('Custom initialization for CovWMu.');
      CovWMu(:,:,:) = options.init.CovWMu;
    end
    % alpha
    if isfield(options.init, 'alpha')
      disp('Custom initialization for alpha.');
      alpha(:) = options.init.alpha;
    end
% $$$     if isfield(options.init, 'a_alpha')
% $$$       disp('Custom initialization for a_alpha.');
% $$$       a_alpha(:) = options.init.a_alpha;
% $$$     end
% $$$     if isfield(options.init, 'b_alpha')
% $$$       disp('Custom initialization for b_alpha.');
% $$$       b_alpha(:) = options.init.b_alpha;
% $$$     end
    % beta
    if isfield(options.init, 'beta')
      disp('Custom initialization for beta.');
      beta(:) = options.init.beta;
    end
% $$$     if isfield(options.init, 'a_beta')
% $$$       disp('Custom initialization for a_beta.');
% $$$       a_beta(:) = options.init.a_beta;
% $$$     end
% $$$     if isfield(options.init, 'b_beta')
% $$$       disp('Custom initialization for b_beta.');
% $$$       b_beta(:) = options.init.b_beta;
% $$$     end
    % tau
    if isfield(options.init, 'tau')
      disp('Custom initialization for tau.');
      tau(:) = options.init.tau;
    end
% $$$     if isfield(options.init, 'a_tau')
% $$$       disp('Custom initialization for a_tau.');
% $$$       a_tau(:) = options.init.a_tau;
% $$$     end
% $$$     if isfield(options.init, 'b_tau')
% $$$       disp('Custom initialization for b_tau.');
% $$$       b_tau(:) = options.init.b_tau;
% $$$     end
    % nu
    if isfield(options.init, 'nu')
      disp('Custom initialization for nu.');
      nu(:) = options.init.nu;
    end
    % U
    if isfield(options.init, 'U')
      disp('Custom initialization for U.');
      U(:,:) = options.init.U;
    end
  end

  % Initialize second moments for X
  XX = zeros(D,D,N);
  for n=1:N
    XX(:,:,n) = X(:,n)*X(:,n)' + CovX(:,:,n);
  end

  % Initialize expectations for W and Mu
  tauW = zeros(M,D);
  tauMu = zeros(M,1);
  WtauW = zeros(D,D,M);
  MutauMu = zeros(M,1);
  WtauMu = zeros(M,D);
  for m=1:M
    tauW(m,:) = tau(m) * W(m,:);
    tauMu(m) = tau(m) * Mu(m);
    WtauW(:,:,m) = tau(m)*W(m,:)'*W(m,:) + CovWMu(1:D,1:D,m);
    MutauMu(m) = tau(m)*Mu(m)^2 + CovWMu(end,end,m);
    WtauMu(m,:) = tau(m)*Mu(m)*W(m,:) + CovWMu(end,1:D,m);
  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%

function set_priors()

  % Default prior
  p.a_alpha = 1e-5;
  p.b_alpha = 1e-5;
  p.a_beta = 1e-5;
  p.b_beta = 1e-5;
  p.a_tau = 1e-5;
  p.b_tau = 1e-5;

  % Custom prior
  if isstruct(options.prior)

    % Prior for alpha
    if isfield(options.prior, 'a_alpha')
      p.a_alpha(1) = options.prior.a_alpha;
      disp('Custom prior parameter a_alpha');
    end
    if isfield(options.prior, 'b_alpha')
      p.b_alpha(1) = options.prior.b_alpha;
      disp('Custom prior parameter b_alpha');
    end

    % Prior for beta
    if isfield(options.prior, 'a_beta')
      p.a_beta(1) = options.prior.a_beta;
      disp('Custom prior parameter a_beta');
    end
    if isfield(options.prior, 'b_beta')
      p.b_beta(1) = options.prior.b_beta;
      disp('Custom prior parameter b_beta');
    end

    % Prior for tau
    if isfield(options.prior, 'a_tau')
      p.a_tau(1) = options.prior.a_tau
      disp('Custom prior parameter a_tau');
    end
    if isfield(options.prior, 'b_tau')
      p.b_tau(1) = options.prior.b_tau
      disp('Custom prior parameter b_tau');
    end

  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function rotate

% I am not absolutely sure that this is bug free

% $$$ D = size(XX,1);
% $$$ N = size(XX,3);
% $$$ M = size(CC,3);

% R should be non-singular because we need to take log(abs(det(R)))

% X_X_ = sum(XX(:,:,1:(N-1)),3) + x0x0;
% CovA = 1/D * (sum(AA,3) - A'*A);
% $$$ Psi = sum(XX,3) - 2*A*sum(X_X,3) + A*X_X_*A' + traceprod(CovA,X_X_)*eye(D);

% Start with identity matrix which does not make the results worse
R = eye(D);

r = minimize(R(:),@cost,100);
R = reshape(r,[D,D]);

% Compute new statistics
for n=1:N
  XX(:,:,n) = R \ XX(:,:,n) / R';
  X_X(:,:,n) = R \ X_X(:,:,n) / R';
end
x0x0 = R \ x0x0 / R';
A = R \ A * R;
CovA = 1/D * trace(inv(R*R')) * R'*CovA*R;
for d=1:D
  AA(:,:,d) = A(d,:)'*A(d,:) + CovA;
end
C = C*R;
for m=1:M
  CC(:,:,m) = R'*CC(:,:,m)*R;
end

I = eye(D);
Rtmp = randn(D);

% Compute new hyperparameters here or after this function!

  function [f,df] = cost(x)

    R = reshape(x,[D,D]);
    [f,df] = bound(R);
    f = -f;
    df = -df(:);

  end

  function [l,dl] = bound(R)

  % Use QR to evaluate these more efficiently?
  RR = R'*R;
  invR = inv(R); % inv??
  invRR = inv(R*R'); % inv??
  logdetR = logdet(R);
  trinvRR = trace(invRR);

  Psi = sum(XX,3)*invR'/p.Cov_X - 2*inv(p.Cov_X)*sum(p.mu_X'*X,2);

  % A_invRR_A = A'*invRR*A + trinvRR*CovA;
  % diag_RA_invRR_AR = diag(R'*A_invRR_A*R);
  diagRCCR = diag(R'*sum(WtauW,3)*R); % not efficient

  % Bound terms
  lpX = -0.5 * traceprod(invRR,Psi);
  lqX = N * logdetR;
  lpA = 0; %-0.5 * D * sum(log(diag_RA_invRR_AR));
  lqA = 0; %-0.5*D*logabsdet(trinvRR*RR);
  lpC = -0.5 * M * sum(log(diagRCCR));
  lqC = -M * logdetR;

  % Bound
  l = lpX + lpA + lpC - lqX - lqA - lqC;

  % Derivative terms
  dlpX = invR * Psi * invRR;
  dlqX = N * invR;
  dlpA = 0; %-D*(diag(1./diag_RA_invRR_AR)*R'*A_invRR_A) ...
         %+ D*invR*A*R*diag(1./diag_RA_invRR_AR)*R'*A'*invRR ...
         %+ D*traceprod(diag(1./diag_RA_invRR_AR),R'*CovA*R)*invR*invRR;
  dlqA = 0; %-D * invR + D^2 * invR*invRR / trinvRR;
  dlpC = -M * (diag(1./diagRCCR)*R'*sum(CC,3));
  dlqC = -M * invR;

  % Derivative
  dl = (dlpX + dlpA + dlpC - dlqX - dlqA - dlqC)';
  end

end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function x = cvx_reshape( x, s, rmap, cperm, cperm2 )
sx = size( x );
if nargin < 5,
    cperm2 = [];
    if nargin < 4,
        cperm = [];
        if nargin < 3,
            rmap = [];
            if nargin < 2,
                s = sx;
            end
        end
    end
end
usp = cvx_use_sparse( s, nnz( x ), isreal( x ) );
isp = issparse( x );
if ~usp && isp,
    x = full( x );
end
if ~isequal( s, sx ) || ~isempty( rmap ) || ~isempty( cperm ) || ~isempty( cperm2 ),
    if usp && ( numel( x ) > 2147483647 || ~isempty( rmap ) || ~isempty( cperm ) ),
        [ ii, jj, x ] = find( x );
        if ~isempty( rmap ),
            sx( 1 ) = length( rmap );
            rmap = find( rmap );
            ii = rmap( ii );
        end
        if ~isempty( cperm ),
            temp = 1 : sx( 2 );
            temp( cperm ) = temp;
            jj = temp( :, jj )';
        end
        ij = ii + ( jj - 1 ) * sx( 1 ) - 1;
        ii = rem( ij, s( 1 ) ) + 1;
        jj = floor( ij / s( 1 ) ) + 1;
        if ~isempty( cperm2 ),
            temp = 1 : s( 2 );
            temp( cperm2 ) = temp;
            jj = temp( :, jj )';
        end
        x = sparse( ii, jj, x, s( 1 ), s( 2 ) );
        clear ii jj ij
        isp = true;
    else
        if ~isempty( cperm ),
            x = x( :, cperm );
        end
        if ~isempty( rmap ),
            x = x( max( 1, cumsum( rmap ) ), : );
            x( ~rmap, : ) = 0;
        end
        x = reshape( x, s );
        if ~isempty( cperm2 ),
            x = x( :, cperm2 );
        end
    end
end
if usp && ~isp,
    x = sparse( x );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
% /*M-FILE SCRIPT SwarmBat_demo_SO_std MMM SGALAB */ %
% /*==================================================================================================
%  Swarm Optimisation and Algorithm Laboratory Toolbox for MATLAB
%
%  Copyright 2012 The SxLAB Family - Yi Chen - leo.chen.yi@live.co.uk
% ====================================================================================================
% File description:
%
%Appendix comments:
%
%Usage:
%
%===================================================================================================
%  See Also:
%
%===================================================================================================
%===================================================================================================
%Revision -
%Date          Name     Description of Change  email
%16-May-2011   Chen Yi  Initial version        leo.chen.yi@live.co.uk
%26-Jul-2011   Chen Yi  Add eachtimeplot       leo.chen.yi@live.co.uk
%HISTORY$
%==================================================================================================*/
%==================================================================================================*/

% SwarmBat_demo_SO_std Begin

% fresh
clear ;
close ('all');
warning off
% to delete old output_*.txt
!del SwarmBat_O_*.*
% set working path

% begin to count time during calculating
home ;
tic % timer start >>

% data preparation

% open data files

%%%input data files
fid1  = fopen('SwarmBat_I_min_confines.txt' , 'r' );
fid2  = fopen('SwarmBat_I_max_confines.txt' , 'r' );
%fid3  = fopen('SwarmBat_I_soundspeed.txt' , 'r' );
fid4  = fopen('SwarmBat_I_distpara.txt' , 'r' );
fid5  = fopen('SwarmBat_I_population.txt' , 'r' );
fid6  = fopen('SwarmBat_I_randomstep.txt' , 'r' );
fid7  = fopen('SwarmBat_I_max_generation.txt' , 'r' );
fid8  = fopen('SwarmBat_I_loundness_min.txt' , 'r' );
fid9  = fopen('SwarmBat_I_PErate_min.txt','r');
%          fid10 = fopen('INPUT_deta_fitness_max.txt','r');
%          fid11 = fopen('INPUT_max_probability_crossover.txt','r');
%          fid12 = fopen('INPUT_probability_crossover_step.txt','r');
%          fid13 = fopen('INPUT_max_no_change_fitness_generation.txt','r');
%          fid22 = fopen('INPUT_tournament_size.txt','r');
fid10  = fopen('SwarmBat_I_loundness_alpha.txt' , 'r' );
fid11  = fopen('SwarmBat_I_PErate_gamma.txt','r');

fid20= fopen('SwarmBat_I_fHz.txt','r');

% set total test number
fid1000= fopen('SwarmBat_I_testnumber.txt','r');

%output data files
fid120 = fopen('SwarmBat_O_best_fitness.txt','a+');
fid130 = fopen('SwarmBat_O_maxfitness.txt','a+');
fid140 = fopen('SwarmBat_O_minfitness.txt','a+');
fid150 = fopen('SwarmBat_O_meanfitness.txt','a+');
fid160 = fopen('SwarmBat_O_best_solutions.txt','a+');
%          fid19 = fopen('OUTPUT_best_coding_space.txt','w+');
%          fid20 = fopen('OUTPUT_now_generation.txt','w+');
%          fid21 = fopen('OUTPUT_now_probability_crossover.txt','w+');

% begin to load data from file

% read data from these files
disp('/*==================================================================================================*/')
disp('/*  Swarm Optimisation and Algorithm Laboratory Toolbox 1.0.0.1 */ ')
disp('');
disp('    30-Mar-2012 Chen Yi leo.chen.yi@live.co.uk Glasgow ')
disp('/*==================================================================================================*/')
disp('>>>>')
disp(' Begin to Evaluate...Waiting Please ...');

min_confines = fscanf( fid1 , '%f' ); min_confines = min_confines' ;

max_confines = fscanf( fid2 , '%f' ); max_confines = max_confines';
%
%       probability_crossover = fscanf( fid3 , '%f' ); probability_mutation = fscanf(fid4,'%f');
%
population = fscanf( fid5 , '%f' );
%
random_step = fscanf( fid6 , '%f' );
%
max_generation = fscanf( fid7 , '%f' );
%
%soundspeed = fscanf( fid3 , '%f' );
%
beta = fscanf( fid4 , '%f' ); % random distribution parameters
%
loudness_min = fscanf( fid8 , '%f' );
%
PErate_min = fscanf( fid9,'%f' );
%
alpha = fscanf(fid10,'%f');
%
gamma = fscanf(fid11,'%f');
%

%
%       now_probability_crossover = probability_crossover;
%
fHz = fscanf(fid20,'%f'); % Hz = [fmin, fmax] is the frequency range
% %

testnumber = fscanf(fid1000,'%f');

disp(' the total test number is');
disp(testnumber);

% Step into SGALAB()
%
eachtimeplot       = 0;
statisticalplot    = [1,1,1,1,1];
%statisticalplot(1) - plot MIN,MEAN,MAX or not, 0-NO, 1-YES
%statisticalplot(2) - plot mAP or not,          0-NO, 1-YES
%statisticalplot(3) - plot mSTD or not,         0-NO, 1-YES
%statisticalplot(4) - plot mmAP or not,         0-NO, 1-YES
%statisticalplot(5) - plot mmSTD or not,        0-NO, 1-YES

options = { ...
    fHz,                  % 1, Hz range
    random_step,          % 2, bat random walk step
    [],                   % 3, sound speed, 340m/s in default
    beta,                 % 4, distribution parameters
    eachtimeplot,         % 5, plot each test, 1 - plot each time, 0 - no plots
    statisticalplot,      % 6, statistical plots, 1 - yes ,  0 - no plots
    'Roulettewheel',       % 7, selection method
    loudness_min,         % 8, loudness min
    PErate_min,           % 9, Pulse Emission rate min
    alpha,                % 10, alpha for loudness
    gamma   };            % 11, gamma for pulse emission rate

SwarmBat_Procedure_h = timebar('SwarmsLAB::SwarmFish','Total Progress...');

for idx = 1 : 1 : testnumber
    % Output

    timebar( SwarmBat_Procedure_h , idx/testnumber );

    disp('Test NO.'); disp(idx);
    disp('');

    [   fitness_data ,...
        best_decimal_space ,...
        error_status ]= SwarmBat__entry_SO_std...
        ( options,...
        min_confines ,...
        max_confines ,...
        population ,...
        max_generation,...
        [testnumber,idx] );

    disp('');

    if ( error_status ~= 0 )  return ;  end

    %write data to output files

    %   OUTPUT_bestfitness.txt
    fprintf( fid120 , '%f\n', fitness_data(1));

    %   OUTPUT_maxfitness.txt
    fprintf( fid130 , '%f\n', fitness_data(2));

    %   OUTPUT_minfitness.txt
    fprintf( fid140,  '%f\n', fitness_data(3));

    %   OUTPUT_meanfitness.txt
    fprintf(  fid150,  '%f\n', fitness_data(4));

    %   OUTPUT_best_result_space.txt
    fprintf(  fid160,  '%f\n', best_decimal_space );

    % %   OUTPUT_best_coding_space
    % fprintf(  fid19 , '%f\n' , best_coding_space );

    % %   OUTPUT_now_generation.txt
    % fprintf(  fid20, '%f\n' , now_generation );
    %
    % %   OUTPUT_now_probability_crossover.txt
    % fprintf(  fid21, '%f\n' , now_probability_crossover );

end

%close files
close( SwarmBat_Procedure_h );
status = fclose( 'all' );

disp('End SwarmFish Evaluating');
disp('');

disp(' More detail result in text files with " Swarm_O_*.txt " ' )
disp('----------------------------------------------------------------------------------------')
result_files = list_current_dir_files ('SwarmBat_O_*.*')

disp('----------------------------------------------------------------------------------------')

% timer end
toc
% SwarmBat_demo_SO_std End

author: Sable

repo: mcbench-benchmarks


text:
function value = r4_acos ( x )

%*****************************************************************************80
%
%% R4_ACOS evaluates the arc-cosine of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the arc-cosine of X.
%
  value = 0.5 * pi - r4_asin ( x );

  return
end


text:
function F = vertcat(varargin)
%VERTCAT   Vertical concatenation of CHEBFUN3 objects.
%   VERTCAT(G, H, W) is the vertical concatenation of CHEBFUN3 objects F, G
%   and H. This function returns a CHEBFUN3V object.
%
%   [G; H; W] is equivalent to VERTCAT(G, H, W).
%
%   VERTCAT(G) returns the CHEBFUN3 object G.
%
% See also CHEBFUN3V.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( nargin == 1 )
    % VERTCAT of one argument just returns the same thing back to the user:
    F = varargin{1};
elseif ( nargin > 1 && nargin <= 3)
    if ( isa(varargin{2}, 'chebfun3v') )
        f = varargin{1};
        F = varargin{2};
        if ( F.nComponents > 2 )
            error('CHEBFUN:CHEBFUN3:vertcat:tooManyComponents', ...
                'Only CHEBFUN3V objects with 2 or 3 components are valid.');
        else
            Fc = F.components;
            g = Fc{1};
            h = Fc{2};
            F = chebfun3v({f, g, h});
        end

    elseif ( isa(varargin{2}, 'chebfun3' ) )
        % Call the CHEBFUN3V constructor:
        F = chebfun3v(varargin{:});
    else
        error('CHEBFUN:CHEBFUN3:vertcat:badInputs', ...
            'Inputs must be either CHEBFUN3 or CHEBFUN3V objects.');
    end

else
    error('CHEBFUN:CHEBFUN3:vertcat:tooManyInputs', ...
        'Cannot vertically concatenate more than three objects.');
end

end


text:
function depthRefined = refineDepthICP(sequenceName,frameIDtarget, interval, renderK, renderW, renderH)

if ~exist('sequenceName','var')
    sequenceName = '2014-04-11_17-09-54_094959634447';
end

if ~exist('frameIDtarget','var')
    frameIDtarget = 50; %200;
end

if ~exist('interval','var')
    interval = 50;
end

data = loadSUN3Dv2(sequenceName);

% old Kinect
%renderK = [519.1638 0 320; 0 519.1638 240; 0 0 1];
if ~exist('renderK','var')
    renderK = [1035 0 768; 0 1035 636; 0 0 1]; %[690 0 512; 0 690 424; 0 0 1]; %data.camera.RGB.K2render; %[519 0 320; 0 519 240; 0 0 1];
end
if ~exist('renderW','var')
    renderW = 1536; %1024; %data.camera.RGB.width; %640;
end
if ~exist('renderH','var')
    renderH =  1272; %848; %data.camera.RGB.height; %480;
end

frameIDs = max(frameIDtarget-interval,1):min(data.imageTotalFrames,frameIDtarget+interval);

Rts = repmat([eye(3) zeros(3,1)],[1,1,length(frameIDs)]);

itarget = find(frameIDs==frameIDtarget);

SmartRejection = 2;
matched = false(1,length(frameIDs));

for i=[itarget:length(frameIDs) itarget:-1:1]
    frameID = frameIDs(i);
    depth = readDepth(data,frameID,false);
    XYZcamera(:,:,1)=data.camera.D.X .* depth;
    XYZcamera(:,:,2)=data.camera.D.Y .* depth;
    XYZcamera(:,:,3)=depth .* (~isnan(data.camera.D.X) & ~isnan(data.camera.D.Y));
    XYZcamera(:,:,4)=depth>0 & ~isnan(data.camera.D.X) & ~isnan(data.camera.D.Y);
    valid = logical(XYZcamera(:,:,4));
    valid = valid(:)';
    XYZ = reshape(XYZcamera,[],4)';
    XYZ = XYZ(1:3,valid);

    if i==itarget
        XYZtarget = XYZ;
    elseif i>itarget
        [R, t, ER, maxD] = icp(XYZold,XYZ,'Matching','kDtree','SmartRejection',SmartRejection);
        matched(i) = true;
        Rts(:,:,i) = concatenateRts(Rts(:,:,i-1),[R t]);

        [R, t, ER, maxD] = icp(XYZtarget,transformPointCloud(XYZ,Rts(:,:,i)),'Matching','kDtree','SmartRejection',SmartRejection);
        Rts(:,:,i) = concatenateRts([R t],Rts(:,:,i));

    elseif i<itarget
        [R, t, ER, maxD] = icp(XYZold,XYZ,'Matching','kDtree','SmartRejection',SmartRejection);
        matched(i) = true;

        Rts(:,:,i) = concatenateRts(Rts(:,:,i+1),[R t]);

        [R, t, ER, maxD] = icp(XYZtarget,transformPointCloud(XYZ,Rts(:,:,i)),'Matching','kDtree','SmartRejection',SmartRejection);
        Rts(:,:,i) = concatenateRts([R t],Rts(:,:,i));
    end

    XYZold = XYZ;
end

for i=1:length(frameIDs)
    if matched(i) || i==itarget
        frameID = frameIDs(i);
        depth = readDepth(data,frameID,false);
        XYZcamera(:,:,1)=data.camera.D.X .* depth;
        XYZcamera(:,:,2)=data.camera.D.Y .* depth;
        XYZcamera(:,:,3)=depth .* (~isnan(data.camera.D.X) & ~isnan(data.camera.D.Y));
        XYZcamera(:,:,4)=depth>0 & ~isnan(data.camera.D.X) & ~isnan(data.camera.D.Y);

        [~,undistortDepth(:,:,i)] = WarpDepthMatlab(XYZcamera, renderK, Rts(:,:,i), renderW, renderH);
    end
end

undistortDepth(undistortDepth(:)==0) = NaN;

%depthRefined = nanmean(undistortDepth,3);
fprintf('median filtering ...');
tic;
depthMedian = nanmedian(undistortDepth,3);
toc;

% to preserve details: check the value and the range of all values. if it is 25%-75%, keep the original values
fprintf('prctile filtering ...');
tic;
maxV = prctile(undistortDepth,75,3);
toc;
fprintf('prctile filtering ...');
tic;
minV = prctile(undistortDepth,25,3);
toc;
depthRaw = undistortDepth(:,:,itarget);

validCount = sum(double(~isnan(undistortDepth)),3);

selV = ((depthRaw < minV) | (depthRaw > maxV) | isnan(depthRaw)) & (validCount> max(3,0.25*(length(frameIDs)-1)));
ind = find(selV);
depthRefined = depthRaw;
depthRefined(ind) = depthMedian(ind);

depthRefined(isnan(depthRefined(:)))=0;

author: jianxiongxiao

repo: ProfXkit


text:
function gotostart(H)
%GOTOSTART Go to start of movie in movie player

feval(H.fcns.goto_start,[],[],H.hfig);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function h = contour(sF,varargin)
% spherical contour plot
%
% Syntax
%   contour(v,data)
%
% Input
%   sF - @S2Fun
%
% Options
%   contours - number of contours
%
% See also
%   S2Fun/plot S2Fun/contourf
%

% plot the function values
h = plot(sF,'contour',varargin{:});

% remove output if not required
if nargout == 0, clear h; end

end

author: mtex-toolbox

repo: mtex

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


text:
fid = fopen('G:\face-CASIA-WebFace\list.txt');
C = textscan(fid, '%s %d');
fclose(fid);
% C{1} = C{1}(1:2000);
% C{2} = C{2}(1:2000);
fid = fopen('G:\face-CASIA-WebFace\sim_label.txt');
C2 = textscan(fid, '%d');
% C2{1} = C2{1}(1:2000);
fclose(fid);
meanC = caffe('read_mean','D:\ThirdPartyLibrary\caffe\examples\siamese\mean.proto');

matcaffe_init(1,'D:\ThirdPartyLibrary\caffe\examples\siamese\mnist_siamese_deploy.prototxt','D:\ThirdPartyLibrary\caffe\examples\siamese\siamese_iter_86000.caffemodel');
num = length(C{2});
label = C2{1}(1:num/2);
% num = floor(num /80) * 2;
AllFeature = zeros(320,num);
for i = 1 : floor(num/100)
    disp([i floor(num/100)]);
    J = zeros(100,100,1,100,'single');
    for j = 1 : 100
        I = imread(C{1}{(i-1)*100+j});
%         I = I(end:-1:1,:);
        I = I';
%         I = I(end:-1:1,:);
        I = single(I) - meanC;
        J(:,:,1,j) = I / 128;
%         J(:,:,1,j) = I;
    end;
    H={J};
    f = caffe('forward',H);
    f = f{1};
    AllFeature(:,(i-1)*100+1:i*100) = reshape(f,[size(AllFeature,1),100]);
end;
J = zeros(100,100,1,100,'single');
for j = 1 : num - floor(num/100) * 100
    I = imread(C{1}{floor(num/100) * 100+j});
    I = single(I') - meanC;
    J(:,:,1,j) = I / 128;
end;
H={J};
f = caffe('forward',H);
f=f{1};
f = reshape(f,[size(AllFeature,1),100]);
AllFeature(:,floor(num/100) * 100+1:num) = f(:,1 : num - floor(num/100) * 100);

thresh1 = [];
thresh2 = [];
p1=1;
p2=1;
num = floor(size(AllFeature,2) / 2);
normX = AllFeature';
for i = 1:num
    if C{2}(i) == C{2}(num+i)
        thresh1 = [thresh1;pdist2(normX(i,:), normX(num+i,:))];
        p1=p1+1;
    else
        thresh2 = [thresh2;pdist2(normX(i,:), normX(num+i,:))];
        p2=p2+1;
    end;
end;
mean(thresh1) / 4 + mean(max(100 - thresh2, 0)) / 4

author: happynear

repo: FaceVerification