text:
function [H] = eeg_regress(data,robust,GP,labels)

% eeg_regress - plot multiple linear regression results
%
% [H] = eeg_regress(data,[robust],[GP],labels)
%
% data is NxP matrix, with N observations, P-1 predictors
% and the last column is the observed values to predict.
% So, if only 2 columns are given it does a simple linear
% regression.
%
% robust = 1, use robust least-squares regression
% robust = 0, use least-squares regression
% note, current robust is forced to 1.
%
% GP is Nx1 grouping variable
%
% labels = cell strings with variable names, keep to
% 5 characters or less for best format of figure text
%
% H is an array of handles to figures
%
% All the results are plotted to new figures, which can
% be saved/exported using print(H(i)...)
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:53 $

% Licence:  GNU GPL, no express or implied warranties
% History:  11/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('data','var'),
  fprintf('...no input data.\n\n');
  return
end

if ~exist('robust','var'),
  robust = 1;
else,
  warning('must use robust method, at present');
  robust = 1; % cheat for now, as some plotting stuff below requires it
end

if ~exist('labels','var'),
  for i = 1:size(data,2),
    if i == size(data,2),
      labels{i} = 'Y';
    else
      labels{i} = sprintf('X%d',i);
    end
  end
end

if ~exist('GP','var'),
  error('GP input variable is undefined\n');
end

plotmatrix(data); % to explore relationships
H(1) = gcf;

y = data(:,end);           % values to fit from last column
e = ones(length(data),1);  % for constant term in model
X = [e data(:,1:end-1)];   % predictor values from first columns

if ~robust,
  %beta = X\y;              % 3 x 1 matrix, const, beta1 & beta2
  [B,BCI,R,RCI,STATS] = regress(y,X,.05);
  P = X*B;      % fitted values

  % use B slope, with BCI intercepts, not strictly correct
  CI = [ [BCI(1,1); B(2:end)]  [BCI(1,2); B(2:end)] ];
  YCI = [X*CI(:,1),  X*CI(:,2) ];

  % Note if BCI(x,:) < or > 0, B(x) is significant

  % define figure rows & columns
  cols = 6;
  rows = size(X,2);

else
  [B,BCI,R,RCI,STATS] = regress(y,X,.05);
  clear B BCI R RCI;

  R2 = STATS(1);
  F  = STATS(2);
  Fp = STATS(3); % for total regression
  clear STATS;

  [B,STATS] = robustfit(X(:,2:end),y);

  P = X*B;           % predicted values
  R = STATS.resid;   % residuals

  BCI = [B - 2*STATS.se, B + 2*STATS.se]; % beta CI
  YCI = [P - 2*STATS.s,  P + 2*STATS.s ]; % predicted CI

  % define figure rows & columns
  cols = 5;
  rows = size(X,2);

end

for fig = 1:2,

  % Setup figure
  H(end+1) = figure('color',[0 0 0]);
  pos = get(gcf,'Position');
  set(gcf,'Position',[pos(1)-(pos(3)/2) pos(2)-pos(4) pos(3)*2 pos(4)*2]);  % double width/height

  colormap(prism);
  fontsize = 8;
  fontcolor = [1 1 1]; % white

  % plot observed, predicted and 95% CI
  if fig < 2,
    subplot(2,3,1);
  else
    subplot(rows,cols,1);
  end
  if max(GP),
    scatter(P,y,5,GP, 'filled'); hold on;
  else,
    scatter(P,y,5,'m','filled'); hold on;
  end
  %scatter(P,P,5,'w','filled');
  plot(P,P,'w-',...
    P,YCI(:,1),'w:',...
    P,YCI(:,2),'w:');
  %legend('predicted','95% CI',0);
  set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);

  title{1} = sprintf('y'' = %8.2f',B(1));
  for j = 2:length(B),
    title{1} = sprintf('%s + %6.2f*X%d',title{1},B(j),j-1);
  end
  Htitle = get(gca,'title');
  set(Htitle,'string',title{1}, 'fontsize',fontsize,'color',fontcolor);

  xlab = get(gca,'xlabel');
  set(xlab,'string','Y predicted','fontsize',fontsize,'color',fontcolor);
  ylab = get(gca,'ylabel');
  set(ylab,'string','Y observed', 'fontsize',fontsize,'color',fontcolor);

  % plot observed, predicted and 95% CI against Xmean
  if fig < 2,
    subplot(2,3,2);
  else
    subplot(rows,cols,2);
  end
  % calculate mean X values
  Xmean = mean(X(:,2:end),2);
  % plot the observed & predicted values
  if max(GP),
    scatter(Xmean,y,5,GP, 'filled'); hold on;
  else,
    scatter(Xmean,y,5,'m','filled'); hold on;
  end
  scatter(Xmean,P,5,'w','filled');
  set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);
  xlab = get(gca,'xlabel');
  set(xlab,'string','mean Xi', 'fontsize',fontsize,'color',fontcolor);
  ylab = get(gca,'ylabel');
  set(ylab,'string','Y observed','fontsize',fontsize,'color',fontcolor);

  % curve fitting, cubic polynomial
  [ Xsort, k ] = sort(Xmean);
  [p,S] = polyfit(Xsort,P(k),3);
  polfit = polyval(p,Xsort,S);
  % % plot observed, polyfit and CI
  plot(Xsort,polfit,'w-');
  % legend('data','polyfit',0);

  % plot residuals against predicted
  if fig < 2,
    subplot(2,3,4);
  else
    subplot(rows,cols,3);
  end
  if max(GP),
    scatter(P,R,5,GP,'filled'); hold on
  else,
    scatter(P,R,5,'m','filled'); hold on
  end
  line(P,zeros(size(R)),'linestyle','-','color','w');
  set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);
  xlab = get(gca,'xlabel');
  set(xlab,'string','Y predicted','fontsize',fontsize,'color',fontcolor);
  ylab = get(gca,'ylabel');
  set(ylab,'string','residual', 'fontsize',fontsize,'color',fontcolor);

  % plot normal QQ plot of residuals
  if fig < 2,
    subplot(2,3,5);
  else
    subplot(rows,cols,4);
  end
  qqplot(R);
  QQtitle = get(gca,'title');
  set(QQtitle,'string','QQ plot','fontsize',fontsize,'color',fontcolor);
  xlab = get(gca,'xlabel');
  set(xlab,'string','','fontsize',fontsize,'color',fontcolor);
  ylab = get(gca,'ylabel');
  set(ylab,'string','','fontsize',fontsize,'color',fontcolor);
  set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);

  % print model statistics
  if robust,

    %    [B,STATS] = ROBUSTFIT(...) also returns a STATS structure
    %     containing the following fields:
    %         stats.ols_s     sigma estimate (rmse) from least squares fit
    %         stats.robust_s  robust estimate of sigma
    %         stats.mad_s     MAD estimate of sigma; used for scaling
    %                         residuals during the iterative fitting
    %         stats.s         final estimate of sigma, the larger of robust_s
    %                         and a weighted average of ols_s and robust_s
    %         stats.se        standard error of coefficient estimates
    %         stats.t         ratio of b to stats.se
    %         stats.p         p-values for stats.t
    %         stats.coeffcorr estimated correlation of coefficient estimates
    %         stats.w         vector of weights for robust fit
    %         stats.h         vector of leverage values for least squares fit
    %         stats.dfe       degrees of freedom for error
    %         stats.R         R factor in QR decomposition of X matrix

    if fig < 2,
      subplot(2,3,3);
    else
      subplot(rows,cols,5);
    end
    axis off;
    txt{1} = sprintf('%6s %8s %8s','','Beta','SE');
    for i = 1:length(B),
      if i == 1,
        predictor = 'CONST';
      else
        if i-1 <= length(labels),
          predictor = labels{i-1};
        end
      end
      txt{i+1} = sprintf('%6s %8.4f %8.4f',predictor,B(i),STATS.se(i));
    end
    text(-.2,.5,txt,'Fontname','courier','fontsize',fontsize,'color',fontcolor)

    txt = '';
    txt{1} = sprintf('%6s %8s %6s %6s %5s','','F','DF','sig','r^2');
    txt{2} = sprintf('%6s %8.4f %6.2f %6.4f %4.2f','ALL',F,STATS.dfe,Fp,R2);
    for i = 1:length(B),
      if i == 1,
        predictor = 'CONST';
      else
        if i-1 <= length(labels),
          predictor = labels{i-1};
        end
      end
      txt{i+2} = sprintf('%6s %8.4f %6.2f %6.4f',predictor,STATS.t(i),STATS.dfe,STATS.p(i));
    end

    % output the correlation matrix
    datacorr = corrcoef(data);
    rtxt = num2str(datacorr,' %+5.2f');
    txt2{1} = '    Correlation Matrix';
    txt2{2} = '   ';
    for r = 1:length(datacorr),
      if r < length(datacorr),
        txt2{2}   = sprintf('%5s %5s',txt2{2},sprintf('X%d',r));
        txt2{r+2} = sprintf('%5s %s',labels{r},rtxt(r,:));
      else
        txt2{2}   = sprintf('%5s %5s',txt2{2},'Y');
        txt2{r+2} = sprintf('%5s %s',labels{r},rtxt(r,:));
      end
    end

    if fig < 2,
      text(-.2,.1,txt ,'Fontname','courier','fontsize',fontsize,'color',fontcolor);
      text(-.1,1 ,txt2,'Fontname','courier','fontsize',fontsize,'color',fontcolor);
    else,
      text(-.2,.1,txt ,'Fontname','courier','fontsize',fontsize,'color',fontcolor);
      text(-.1,1 ,txt2,'Fontname','courier','fontsize',fontsize,'color',fontcolor);
    end
    txt = '';
    txt2 = '';
  else

    % stats(1) = R^2, stats(2) = F, stats(3) = p; % for total regression
    if fig < 2,
      subplot(2,3,3);
    else
      subplot(rows,cols,5);
    end
    axis off;
    txt{1} = sprintf('%12s %10s  %10s  %10s\n','','R^2','F','sig');
    txt{2} = sprintf('%12s%10.4f  %10.4f  %10.4f','MODEL',STATS(1),STATS(2),STATS(3));
    text(-.2,.9,txt,'Fontname','courier','fontsize',fontsize,'color',fontcolor)

    if fig > 1,
      subplot(rows,cols,6); rcoplot(R, RCI);
    end
  end

end

% plot observed, predicted and 95% CI against individual predictors
for i = 2:size(X,2),
  subplot(rows,cols,cols*i-(cols-1));
  % plot the observed values
  if max(GP),
    scatter(X(:,i),y,5,GP, 'filled'); hold on;
  else
    scatter(X(:,i),y,5,'w','filled'); hold on;
  end
  set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);

  % Calculate predicted values for just this predictor
  beta = [B(1);B(i)];
  P = [X(:,1),X(:,i)]*beta;
  R = y - P;
  YCI = [P - 2*STATS.s,  P + 2*STATS.s ]; % predicted CI

  %scatter(X(:,i),P,5,'w','filled');
  plot(X(:,i),P,'w-',...
    X(:,i),YCI(:,1),'w:',...
    X(:,i),YCI(:,2),'w:');
  %legend('predicted','95% CI',0);

  title{1} = sprintf('y'' = %8.2f',B(1));
  title{1} = sprintf('%s + %6.2f*{\\bfX%d}',title{1},B(i),i-1);
  %Htitle = get(gca,'title');
  %set(Htitle,'string',title{1},'interpreter','tex');
  ylab = get(gca,'ylabel');
  set(ylab,'string',title{1},'fontsize',fontsize,'color',fontcolor);

  % plot residuals against predicted
  subplot(rows,cols,cols*i-(cols-3));
  if max(GP),
    scatter(X(:,i),R,5,GP, 'filled'); hold on;
  else
    scatter(X(:,i),R,5,'w','filled'); hold on;
  end
  line(X(:,i),zeros(size(R)),'linestyle','-','color','w');
  %xlab = get(gca,'xlabel');
  %set(xlab,'string','Xi','fontsize',fontsize,'color',fontcolor);
  ylab = get(gca,'ylabel');
  set(ylab,'string','residual','fontsize',fontsize,'color',fontcolor);
  set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);

  % plot normal QQ plot of residuals
  subplot(rows,cols,cols*i-(cols-4));
  qqplot(R);
  QQtitle = get(gca,'title');
  set(QQtitle,'string','QQ plot','fontsize',fontsize,'color',fontcolor);
  xlab = get(gca,'xlabel');
  set(xlab,'string','','fontsize',fontsize,'color',fontcolor);
  ylab = get(gca,'ylabel');
  set(ylab,'string','','fontsize',fontsize,'color',fontcolor);
  set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);
end

return

author: PatternRecognition

repo: OpenBMI


text:
function [ TreeMatric,Cost ] = DirectedMaximumSpanningTree( OriginalCostMatric,Root )
% MST on a directed graph
% Chu-Liu/Edmonds Algorithm:
%1 Discard the arcs entering the root if any; For each node other than the root, select the entering arc with the highest cost; Let the selected n-1 arcs be the set S.
%  If no cycle is formed, G( N,S ) is a MST. Otherwise, continue.
%2 For each cycle formed, contract the nodes in the cycle into a pseudo-node (k), and modify the cost of each arc which enters a node (j) in the cycle from some node (i)
%  outside the cycle according to the following equation.  c(i,k)=c(i,j)-(c(x(j),j)-min_{j}(c(x(j),j)) where c(x(j),j) is the cost of the arc in the cycle which enters j.
%3 For each pseudo-node, select the entering arc which has the smallest modified cost; Replace the arc which enters the same real node in S by the new selected arc.
%4 Go to step 2 with the contracted graph.

% This code is written by Lowell Guangdi, Email: lowellli121@gmail.com
Dim = size( OriginalCostMatric,2 );
for p = 1:Dim, OriginalCostMatric( p,p ) = 0;end

Min = min(min(OriginalCostMatric));
OriginalCostMatric = OriginalCostMatric - Min;

CostMatric = OriginalCostMatric;
CostMatric( :,Root )= zeros( Dim,1 );
while 1
     TreeMatric = CostMatric;
     % select out the maximal weights of arcs upon each node.
     for p = 1 : Dim
         if p ~= Root
            [ LocalMax,Index ] = max( TreeMatric( :,p ) );
            TreeMatric( :,p ) = zeros( Dim,1 );
            TreeMatric( Index,p ) = LocalMax;
         end
     end
     % test the existence of cycle
     [ CNumber, Component ] = conncomp( biograph( TreeMatric ),'Weak', true );
     if CNumber == 1,break; end
     % Cycle exists
     RootCluster = Component( Root );
     RootSharedNode = find( Component == RootCluster );
     if RootCluster == 1
        CaredCluster = 2;
     elseif RootCluster > 1
        CaredCluster = 1;
     end
     % change the weight here
     ClusterNode = find( Component == CaredCluster );
     CostMatric = SearchCycleNode(ClusterNode,TreeMatric, CostMatric,OriginalCostMatric,RootSharedNode);
end
TreeMatric = CostMatric;
Cost = sum(sum( TreeMatric ));
end

author: Sable

repo: mcbench-benchmarks


text:
% subsref() - index eegdata class
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function res = subsref(obj,s)

    if strcmpi(s(1).type, '.')
        res = builtin('subsref', obj, s);
        return;
    end;

    tmpMMO = memmapfile(obj.dataFile, 'writable', obj.writable, 'format', { 'single' obj.dimensions 'x' });

    subs = s(1).subs;
    finaldim = cellfun('length', subs);

    % one dimension input
    % -------------------
    if length(s) > 1 || ~strcmpi(s(1).type, '()')
        error('MMO can only map single array data files');
    end;

    % deal with transposed data
    % -------------------------
    if obj.transposed, s = transposeindices(obj, s); end;

    % convert : to real sizes
    % -----------------------
    lastdim = length(subs);
    if isstr(subs{end}) && ndims(obj) > lastdim
        for index = lastdim+1:ndims(obj)
            if index > length(obj.dimensions)
                subs{index} = 1;
            else
                subs{index} = [1:obj.dimensions(index)];
            end;
        end;
    end;
    for index = 1:length(subs)
        if isstr(subs{index}) % can only be ":"
            if index > length(obj.dimensions)
                subs{index} = 1;
            else
                subs{index} = [1:obj.dimensions(index)];
            end;
        end;
    end;
    finaldim = cellfun(@length, subs);
    finaldim(lastdim) = prod(finaldim(lastdim:end));
    finaldim(lastdim+1:end) = [];

    % non-transposed data
    % -------------------
    res = tmpMMO.data.x(subs{:});
    if length(finaldim) == 1, finaldim(2) = 1; end;
    res = reshape(res, finaldim);
    if obj.transposed
        if finaldim(end) == 1, finaldim(end) = []; end;
        if length(finaldim) <= 2, res = res';
        else
            res = reshape(res, [finaldim(1)*finaldim(2) finaldim(3)])';
            res = reshape(res, [finaldim(3) finaldim(1) finaldim(2)]);
        end;
    end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% fbp_fan_arc_example.m
% example of how to use fbp_fan_arc.m
% Copyright Jeff Fessler, University of Michigan

if ~isvar('sino')
	down = 4;
	ig = image_geom('nx', 512, 'ny', 504, 'fov', 500, 'down', down);
	sg = sino_geom('fan', 'nb', 888, 'na', 984, 'ds', 1.0, ...
		'dsd', 949, 'dod', 408, 'offset_s', 0.25, ...
		'orbit', 360, ...
		'down', down);

	ell = [0 0 200 200 0 1; 0 0 10 10 0 1];
	ell = [0 0 200 200 0 0; 0 0 10 10 0 1];
	ell = [20 -15 225 150 30 10; 40 0 10 15 0 1];
%	ell = [100 100 ds/4 ds/4 0 1]; % point source

	sino = ellipse_sino(sg, ell, 'oversample', 2);

	xtrue = ellipse_im(ig, ell, 'oversample', 4);

	% system object
	G = Gtomo2_dscmex(sg, ig);

	im plc 2 3
	im(1, xtrue, 'x'), cbar
	im(4, sino, 'sino'), cbar
prompt
end

if 0 % examine recon of flat fan-beam sinogram
	sino = repmat(abs(sg.s) < 0.8 * max(sg.s), [1 sg.na]);
	im(sino)
end

% fan-beam reconstruction
if ~isvar('recon')
	recon = fbp_fan_arc(sino, G, 'ramp');
end

if 1 % compare to new fbp2() method
	tmat = fbp2(sg, ig, 'type', 'std:mat');
	rmat = fbp2(sino, tmat);
	max_percent_diff(rmat, recon)
prompt
end

if 0 % compare to new fbp2() method
	tmex = fbp2(sg, ig, 'type', 'std:mex');
	rmex = fbp2(sino, tmex);
	mmask = (rmat > 0);
	max_percent_diff(rmex .* mmask, recon .* mmask)
	im(2, rmat .* mmask, 'mat'), cbar
	im(3, rmex .* mmask, 'mex'), cbar
return
end

if 0 % examine "aliasing"
	clim = [9.5 10.5];
%	clim = [0.8 1.2];
%	clim = [-0.2 0.2];
	im clf, im(recon, clim), cbar
return
end

if im
	im(2, recon, 'FBP matlab'), cbar
	im(5, recon - xtrue, 'error'), cbar
	iy = ig.ny/2; ix = 1:ig.nx;
	subplot(133)
	plot(ix, xtrue(ix,iy), '-', ix, recon(ix,iy), '--')
	axis([1 ig.nx -0.5 10.5]), legend('true', 'recon', 4)
end

if has_aspire % check consistency with aspire
	dir = test_dir;
	f.sino = [dir 'sino.fld'];
	f.image = [dir 'image.fld'];
	f.dsc = [dir 't.dsc'];
	fld_write(f.sino, sino)

	char_array_write(f.dsc, G.arg.args)
	f.win = 'boxcar,1,0,1';
	com = sprintf('echo y | i fbp dsc %s %s %s %s', ...
		f.image, f.sino, f.dsc, f.win)
%	eval(['!' com])
	disp(os_run(com))

	if 0 % compare filters
		tmp = fld_read('fft_filt.fld');
		sum(tmp) / sum(test)
%		im clf, plot([tmp test])
		max_percent_diff(test, tmp)
	return
	end

	if 0 % compare filtered projections
		% seem to match except for slight shift??
		tmp = fld_read('proj_filt.fld');
		im clf, plot([tmp(:,1) test(:,1)])
		plot(tmp(:,1)-test(:,1))
		max_percent_diff(test, tmp)
		sum(tmp(:)) / sum(test(:))
		minmax(test)
		minmax(tmp)
		minmax(test-tmp)
	return
	end

	im_asp = fld_read(f.image);
	im(233, im_asp, 'aspire'), cbar
	good = recon ~= 0;
	diff = im_asp - recon;
	diff = diff .* good;
	im(236, diff, 'aspire-matlab'), cbar
	im(234, (im_asp ~= 0) - (recon ~= 0), 'aspire-matlab support'), cbar
	max_percent_diff(recon.*good, im_asp.*good)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

% You need to return the following variables correctly.
Z = zeros(size(X, 1), K);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%
U_reduce = U(:, 1:K);
x = X';
Z = x' * U_reduce;

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

end

author: scruel

repo: Notes-ML-AndrewNg


text:
% DEMSHAREDVARGPLVM1 Run the shared variational GPLVM on various kinds of data.
% DESC
%
% COPYRIGHT: Andreas C. Damianou, Carl Henrik Ek,  2011
% SEEALSO : svargplvmPrepareData
%
% VARGPLVM

%{
%---- How to run this demo:
% Vargplvm created data
close all; initLatent='ppca'; dataSetNames = {}; toyDataCreate ='vargplvm';initVardistIters = 400; itNo = [100 200 1500 500 500 400];  demSharedVargplvm1

% Fols data (make sure to use a linard2 kernel)
close all; initLatent='pca'; dataSetNames = {}; toyDataCreate ='fols';initVardistIters = 170; itNo = 300; dataToKeep=50;  demSharedVargplvm1

% Human pose data with features extracted from silhouettes
clear ;close all; initLatent='pca3'; dataSetNames = {}; toyDataCreate='humanPose';initVardistIters = 160; itNo = [200 200 200 200 200 200 200 200]; dataToKeep=418;  demSharedVargplvm1
clear ;close all; experimentNo=1;initLatent='pca3'; dataSetNames = {};toyDataCreate='humanPose';initVardistIters = 180; itNo = [500 200 200 200 200 200 200 200 200];  demSharedVargplvm1

% Human pose data with the whole silhouette (TEMPdemSharedVargplvm is like
% demSharedVargplvm but just uses more latent dimensions)
clear ;close all; experimentNo=6; imageSilhouette=1;initLatent='ppca'; latentDimPerModel=7;dataSetNames = {};toyDataCreate='humanPose';initVardistIters = 380; itNo = [500 200 200 200 200 200 200 200 200 200];  indPoints=100;TEMPdemSharedVargplvm

% TODO: ___
clear ;close all; experimentNo=5;initLatent='ppca'; mappingKern = 'rbfardjit';dataSetNames = {};toyDataCreate='humanPose';initVardistIters = 180; itNo = [500 200 200 200 200 200 200 200 200];  demSharedVargplvm1
%}

%___

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

% Define constants (in a manner that allows other scripts to parametrize
% this one).
if ~exist('experimentNo') ,  experimentNo = 404;      end
if ~exist('itNo')         ,  itNo = [500 500];              end     % Default: 2000
if ~exist('indPoints')    ,  indPoints = 80;          end     % Default: 49
if ~exist('latentDimPerModel')    ,  latentDimPerModel = 3;          end
if ~exist('latentDim'), latentDim = 5; end
if ~exist('numSharedDims'), numSharedDims = 2; end
% Set to 1 to use dynamics or to 0 to use the standard var-GPLVM
if ~exist('dynUsed')      ,  dynUsed = 0;             end
if ~exist('initVardistIters'), initVardistIters = 100;      end
if ~exist('dynamicKern')   ,     dynamicKern = {'rbf', 'white', 'bias'}; end
if ~exist('mappingKern')   ,  mappingKern = {'rbfard2', 'white'}; end
% if ~exist('mappingKern'),  mappingKern = {'rbfard2', 'bias', 'white'}; end
% Set to 1 to tie the inducing points with the latent vars. X
if ~exist('fixInd')        ,     fixInd = 0;    end
% 0.1 gives around 0.5 init.covars. 1.3 biases towards 0.
if ~exist('vardistCovarsMult'),  vardistCovarsMult=1.3;                  end
% Set to empty value {} to work with toy data
if ~exist('dataSetNames')    ,    dataSetNames = {};    end
if ~exist('invWidthMultDyn'),    invWidthMultDyn = 100;                     end
if ~exist('invWidthMult'),       invWidthMult = 5;                     end
if ~exist('initX'),     initX ='ppca';   end % That's for the dynamics initialisation
if ~exist('initLatent'),     initLatent ='ppca';   end % That's for the whole latent space init.
if ~exist('dataToKeep'), dataToKeep = -1; end
if ~exist('toyDataCreate'), toyDataCreate = 'fols'; end
if ~exist('doPredictions'), doPredictions = 0; end
if ~exist('dataType'), dataType = 'default'; end
if ~exist('enableParallelism'), enableParallelism = 1; end

if ~exist('Yall')
    svargplvmPrepareData;
end

%%

%-- Load datasets
for i=1:numberOfDatasets
    Y = Yall{i};
    if dataToKeep ~= -1 && dataToKeep <= size(Y,1)
        Y = Y(1:dataToKeep,:);
    end
    dims{i} = size(Y,2);
    N{i} = size(Y,1);
    indTr = 1:N{i};
    indTs = setdiff(size(Y,1), indTr);
    Ytr{i} = Y(indTr,:); %Yts = Y(indTs,:);
    t{i} = linspace(0, 2*pi, size(Y, 1)+1)'; t{i} = t{i}(1:end-1, 1);
    timeStampsTraining{i} = t{i}(indTr,1); %timeStampsTest = t(indTs,1);
    d{i} = size(Ytr{i}, 2);
end

for i=2:numberOfDatasets
    if N{i} ~= N{i-1}
        error('The number of observations in each dataset must be the same!');
    end
end

%-- Options for the models
for i=1:numberOfDatasets
    % Set up models
    options{i} = vargplvmOptions('dtcvar');
    options{i}.kern = mappingKern; %{'rbfard2', 'bias', 'white'};
    %indPoints = 80; %%%%%
    options{i}.numActive = indPoints;
    options{i}.optimiser = 'scg2';

    % !!!!! Be careful to use the same type of scaling and bias for all
    % models!!!

    % scale = std(Ytr);
    % scale(find(scale==0)) = 1;
    %options.scaleVal = mean(std(Ytr));
    options{i}.scaleVal = sqrt(var(Ytr{i}(:)));
    if fixInd
        options{i}.fixInducing=1;
        options{i}.fixIndices=1:size(Ytr{i},1);
    end
end

X_init = svargplvmInitLatentSpace(Ytr, d, options, initLatent, latentDim, latentDimPerModel, numSharedDims);

latentDim = size(X_init,2);

% Free up some memory
clear('Y')

%-- Create the sub-models: Assume that for each dataset we have one model.
% This can be changed later, as long as we find a reasonable way to
% initialise the latent spaces.
for i=1:numberOfDatasets
    %---- Here put some code to assign X to the global common X which must
    % be created by doing pca in the concatenation of Y's...After this
    % point, model{i}.X will be the same for all i's. TODO...
    fprintf(1,'# Creating the model...\n');
    options{i}.initX = X_init;
    model{i} = vargplvmCreate(latentDim, d{i}, Ytr{i}, options{i});
    model{i}.X = X_init; %%%%%%%
    model{i} = vargplvmParamInit(model{i}, model{i}.m, model{i}.X);
    model{i}.X = X_init; %%%%%%%

    inpScales = invWidthMult./(((max(model{i}.X)-min(model{i}.X))).^2); % Default 5
    %inpScales(:) = max(inpScales); % Optional!!!!!
    model{i}.kern.comp{1}.inputScales = inpScales;

    if strcmp(model{i}.kern.type, 'rbfardjit')
        model{i}.kern.inputScales = model{i}.kern.comp{1}.inputScales;
    end
    params = vargplvmExtractParam(model{i});
    model{i} = vargplvmExpandParam(model{i}, params);
    model{i}.vardist.covars = 0.5*ones(size(model{i}.vardist.covars)) + 0.001*randn(size(model{i}.vardist.covars));

    %-------- Add dynamics to the model -----
    if dynUsed
        fprintf(1,'# Adding dynamics to the model...\n');
        optionsDyn{i}.type = 'vargpTime';
        optionsDyn{i}.t=timeStampsTraining{i};
        optionsDyn{i}.inverseWidth=invWidthMultDyn; % Default: 100
        optionsDyn{i}.initX = X_init; % initX; % That's for the dynamics

        kern = kernCreate(t{i}, dynamicKern); % Default: {'rbf','white','bias'}

        %-- Initialize each element of the compound kernel
        svargplvmInitDynKern

        optionsDyn{i}.kern = kern;
        optionsDyn{i}.vardistCovars = vardistCovarsMult; % 0.23 gives true vardist.covars around 0.5 (DEFAULT: 0.23) for the ocean dataset

        % Fill in with default values whatever is not already set
        optionsDyn{i} = vargplvmOptionsDyn(optionsDyn{i});
        model{i} = vargplvmAddDynamics(model{i}, 'vargpTime', optionsDyn{i}, optionsDyn{i}.t, 0, 0,optionsDyn{i}.seq);

        fprintf(1,'# Further calibration of the initial values...\n');
        model{i} = vargplvmInitDynamics(model{i},optionsDyn{i});

        %  to also not learn the last kernel's variance
        if numel(kern.comp) > 1 && exist('learnSecondVariance') && ~learnSecondVariance
            fprintf(1,'# The variance for %s in the dynamics is not learned!\n',kern.comp{end}.type)
            model{i}.dynamics.learnSecondVariance = 0;
            model{i}.dynamics.kern.comp{end}.inverseWidth = model{i}.dynamics.kern.comp{1}.inverseWidth/10; %%% TEMP
        end
    end

    model{i}.beta=1/(0.01*var(model{i}.m(:)));
    prunedModelInit{i} = vargplvmPruneModel(model{i});
    %disp(model{i}.vardist.covars)
end

if experimentNo == -1
    experimentNo = globalExperimentNumber('sharedVargplvm', 'dataType');
end

%modelInit = model;%%%TEMP

%--  Unify models into a structure
model = svargplvmModelCreate(model);
model.dataSetNames = dataSetNames;
model.initLatent = initLatent;
model.experimentNo = experimentNo;
model.dataType = dataType;
%%---
capName = dataType;
capName(1) = upper(capName(1));
modelType = model.type;
modelType(1) = upper(modelType(1));
fileToSave = ['dem' capName modelType num2str(experimentNo) '.mat'];
%%---

%-- Define what level of parallelism to use (w.r.t submodels or/and w.r.t
% datapoints).
if enableParallelism
    if model.numModels > 8
        fprintf('# Parallel computations w.r.t the submodels!\n');
        model.parallel = 1;
        model = svargplvmPropagateField(model,'parallel', 1);
    elseif model.N > 15
        fprintf('# Parallel computations w.r.t the datapoints!\n');
        model.vardist.parallel = 1;
        for i=1:model.numModels
            model.comp{i}.vardist.parallel = 1;
        end
    end
end

% Force kernel computations
params = svargplvmExtractParam(model);
model = svargplvmExpandParam(model, params);

%%

display = 1;
%%%% Optimisation
% do not learn beta and sigma_f for few iterations for intitialization
if initVardistIters ~=0
    model.initVardist = 1;
    model = svargplvmPropagateField(model,'initVardist', model.initVardist);
    fprintf(1,'# Intitiliazing the variational distribution...\n');
    model = svargplvmOptimise(model, display, initVardistIters); % Default: 20
    %fprintf(1,'1/b = %.4d\n',1/model.beta);

    modelInitVardist = model;
    model.initVardistIters=initVardistIters;
end

model.initVardist = 0;
model = svargplvmPropagateField(model,'initVardist', model.initVardist);

% Optimise the model.
model.iters = 0;
for i=1:length(itNo)
    iters = itNo(i); % default: 2000
    fprintf(1,'\n# Optimising the model for %d iterations (session %d)...\n',iters,i);
    model = svargplvmOptimise(model, display, iters);
    model.iters = model.iters + iters;
    % Save model
    fprintf(1,'# Saving %s\n',fileToSave);
    save(fileToSave, 'model', 'prunedModelInit');
end

for i=1:numberOfDatasets
    figure, bar(model.comp{i}.kern.comp{1}.inputScales);
    title(['Final scales for dataset ' num2str(i)]);
end

if toyData && strcmp(toyDataCreate,'fols')
    figure,plot(model.X(:,1)), hold on, plot(model.X(:,2),'r')
    hold on, plot(model.X(:,3),'g')
end

%%
%--------------------- PREDICTIONS --------------%
if ~doPredictions
    return
end
if model.numModels ~=2
    error('Predictions cannot be made (currently) with sharedVargplvm unless the number of submodels used is 2.');
end
%%

obsMod = 1; % one of the involved sub-models (the one for which we have the data)
infMod = setdiff(1:2, obsMod);

% Find the dimensions that are shared for obsMod and infMod
if ~exist('sharedDims')
    thresh = max(model.comp{obsMod}.kern.comp{1}.inputScales) * 0.001;
    retainedScales{obsMod} = find(model.comp{obsMod}.kern.comp{1}.inputScales > thresh);
    thresh = max(model.comp{infMod}.kern.comp{1}.inputScales) * 0.001;
    retainedScales{infMod} = find(model.comp{infMod}.kern.comp{1}.inputScales > thresh);
    sharedDims = intersect(retainedScales{obsMod}, retainedScales{infMod});
end

% Find X_* only for the shared dimensions (Xs*):
if ~exist('privateDims')
    privateDims{infMod} = setdiff(1:model.comp{obsMod}.q, sharedDims);
end

testOnTraining=0;

if testOnTraining
    i = size(model.comp{obsMod}.y,1); % last observation
    % Find p(X_* | Y_*): If Y* is not taken from the tr. data, then this step
    % must be an optimisation step of q(x*). X_* and Y_* here refer to the
    % spaces of the submodel obsMod.
    y_star = model.comp{obsMod}.y(i);
    x_star = model.comp{obsMod}.vardist.means(i,:);
    varx_star = model.comp{obsMod}.vardist.covars(i,:);
    [ind, distInd] = nn_class(model.X(:,sharedDims), x_star(:,sharedDims), 10, 'euclidean');

    ZpredMu = zeros(length(ind), size(Ytoy{infMod},2));
    ZpredSigma = zeros(length(ind), size(Ytoy{infMod},2));
    for k=1:length(ind)
        [ZpredMu(k,:), ZpredSigma(k,:)] = vargplvmPosteriorMeanVar(model.comp{infMod}, model.X(ind(k),:));
    end
else
    testInd = 25;
    Yts = Y_test(testInd,:); % Now this is a matrix
    for i=1:size(Yts,1)
        % initialize the latent points using the nearest neighbour
        % from the training data
        dst = dist2(Yts(i,:), model.comp{obsMod}.y);
        [mind, mini] = min(dst);

        Init(i,:) = model.vardist.means(mini,:);
        vardistx = vardistCreate(model.comp{obsMod}.vardist.means(mini,:), model.q, 'gaussian');
        vardistx.covars = model.comp{obsMod}.vardist.covars(mini,:);
        model.comp{obsMod}.vardistx = vardistx;
        display=1;
        iters = 100;
        % Find p(X_* | Y_*) which is approximated by q(X_*)
        [x_star, varx_star, modelUpdated] = vargplvmOptimisePoint(model.comp{obsMod}, vardistx, Yts(i,:), display, iters);%%%
        numberOfNN = 10;
        % Now we selected a datapoint X_* by taking into account only the
        % private dimensions for Y. Now, based on the shared dimensions of
        % that, we select the closest (in a NN manner) X from the training data.
        [ind, distInd] = nn_class(model.X(:,sharedDims), x_star(:,sharedDims), numberOfNN, 'euclidean');

        x = model.X(ind,:);
        % Careful:  Since x is coming from the training set, it is somewhat
        % natural that the predicted y's will also closely match the ones
        % from the training set. What makes it interesting, is to show that
        % the shared dimensions are the ones controlling the common
        % ambiguities, so, actually it really matters what indices of
        % training X's we find with the NN based on sharedDims.
        mu = vargplvmPosteriorMeanVar(model.comp{infMod}, x);

        ZpredK{i} = mu;

        ZpredMu(i,:) = mu(1);

        %ZpredSigma(i,:) = sigma;
    end

    errsumFull = sum((ZpredMu - Z_test(testInd,:)).^2);
    errorFull = mean(errsumFull);

    if toyData && strcmp(toyDataCreate,'humanPose')
        for j = 1:1:length(testInd)
            handle = xyzankurVisualise(ZpredK{j}(1,:),1);
            for k = 2:length(ind)
                %  xyzankurModify(handle,ZpredK{j}(k,:));            handle = xyzankurVisualise(ZpredK{j}(1,:),1);
                handle = xyzankurVisualise(ZpredK{j}(k,:),k);
                title(['Pose'])
                pause;
            end
        end
    end
end

author: SheffieldML

repo: vargplvm


text:
function [data,columns] = queryDatabase(db,query)
%QUERYDATABASE Query a database.
%   queryDatabase(DATABASE,QUERY) returns results of the query and the
%   column names as cell arrays.  Any errors will be thrown as
%   real errors.
%
%   Example:
%
%     db = database('foo','bar','baz');
%     query = 'select * from something';
%     [data,names] = queryDatabase(db,query)
%     s = cell2struct(data,names,2);
%     close(db);

% Matthew J. Simoneau
% Copyright 1984-2007 The MathWorks, Inc.

curs = exec(db,query);
error(curs.Message);
curs = fetch(curs);
error(curs.Message)
data = curs.Data;

cnames = columnnames(curs);
% This doesn't work if there are single quotes in the field names.
% eval(['columns = {' cnames(:) '};'])
% Use split instead.
columns = split(cnames(2:end-1),''',''');

close(curs)

if isequal(size(data),[1 1]) && isequal(data{1},'No Data')
    data = cell(0,length(columns));
end

function cellArray = split(str,delimiter)
%SPLIT divides a string into a cell array of strings based on a delimiter.
%   SPLIT(STR,DELIMITER) breaks STR into chunks sectioned by DELIMITER and
%   returns a cell array of these chunks.
%
%     >> split('fooMJSbarMJSbaz','MJS')
%
%     ans =
%
%         'foo'
%         'bar'
%         'baz'

% Matthew J. Simoneau, April 2003

dLength = length(delimiter);
matches = strfind(str,delimiter);
if isempty(matches)
    cellArray{1} = str;
else
    cellArray{1} = str(1:matches(1)-1);
    for i = 1:length(matches)-1
        cellArray{end+1,1} = str(matches(i)+dLength:matches(i+1)-1);
    end
    cellArray{end+1,1} = str(matches(end)+dLength:end);
end


text:
function a = daub4_matrix ( n )

%*****************************************************************************80
%
%% DAUB4_MATRIX returns the DAUB4 matrix.
%
%  Discussion:
%
%    The DAUB4 matrix is the Daubechies wavelet transformation matrix
%    with 4 coefficients.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 4, and a multiple of 2.
%
%    Output, real A(N,N), the matrix.
%
  c = [  0.4829629131445341E+00; ...
         0.8365163037378079E+00; ...
         0.2241438680420133E+00; ...
        -0.1294095225512603E+00 ];

  if ( n < 4 || mod ( n, 2 ) ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DAUB4_MATRIX - Fatal error!\n' );
    fprintf ( 1, '  N must be at least 4 and a multiple of 2.\n' );
    error ( 'DAUB4_MATRIX - Fatal error!' );
  end

  a = zeros ( n, n );

  for i = 1 : 2 : n - 1

    a(i,i)                  =   c(1);
    a(i,i+1)                =   c(2);
    a(i,i4_wrap(i+2,1,n))   =   c(3);
    a(i,i4_wrap(i+3,1,n))   =   c(4);

    a(i+1,i)                =   c(4);
    a(i+1,i+1)              = - c(3);
    a(i+1,i4_wrap(i+2,1,n)) =   c(2);
    a(i+1,i4_wrap(i+3,1,n)) = - c(1);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function lse = mylogsumexp(b)
% does logsumexp across columns
B = max(b,[],2);
if issparse(b)
    lse = log(sum(exp(b-repmat(B,[1 size(b,2)])),2))+B;
else
    lse = log(sum(exp(b-repmatC(B,[1 size(b,2)])),2))+B;
end
end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function [T] = cubic_cubic_intersect(CA,CB,tol)
  % CUBIC_CUBIC_INTERSECT Intersect two cubic Bezier curves.
  %
  % [T] = cubic_cubic_intersect(CA,CB,tol)
  %
  % Inputs:
  %   CA  #4 by 2 list of control point locations
  %   CB  #4 by 2 list of control point locations
  % Outputs:
  %   T  #T by 2 list of parameters between [0,1] of intersects, so that
  %     intersect i occurs at T(i,1) on curve A and at T(i,2) on curve B.
  %
  % See also: cubic_eval, cubic_is_flat, cubic_split
  %

  if nargin<3
    tol = 1e-7;
  end

  %plot_cubic(CA);
  %hold on;
  %arrayfun(@(p) set(p,'Color','r'), plot_cubic(CB));
  %hold off;
  %axis equal;
  %axis([-1 1 -1 1]);
  %drawnow

  %if ~bounding_box_overlap(CA,CB)
  % Avoid function overhead
  if (any(max(CA)<min(CB)) || any(max(CB)<min(CA)))
    T = zeros(0,2);
    return;
  end
  A_flat = cubic_is_flat(CA,tol);
  B_flat = cubic_is_flat(CB,tol);
  % lazy just wait until they're both flat
  if A_flat && B_flat
    s = lineSegmentIntersect([CA(1,:) CA(4,:)],[CB(1,:) CB(4,:)]);
    if s.intAdjacencyMatrix
      T = [s.intNormalizedDistance1To2 s.intNormalizedDistance2To1];
    else
      T = zeros(0,2);
    end
    return;
  end

  [CA1,CA2] = cubic_split(CA,0.5);
  [CB1,CB2] = cubic_split(CB,0.5);
  T11 = cubic_cubic_intersect(CA1,CB1,tol)*0.5+0.5*[0 0];
  T12 = cubic_cubic_intersect(CA1,CB2,tol)*0.5+0.5*[0 1];
  T21 = cubic_cubic_intersect(CA2,CB1,tol)*0.5+0.5*[1 0];
  T22 = cubic_cubic_intersect(CA2,CB2,tol)*0.5+0.5*[1 1];
  T = [T11;T12;T21;T22];
end

author: alecjacobson

repo: gptoolbox


text:
function [Q,T] = forward_kinematics(C,BE,P,dQ)
  % Compute absolute transformations from FK hierarchy and relative bone
  % rotations.
  %
  % Inputs:
  %   C  #C by dim list of joint positions
  %   BE  #BE by 2 list of bone edge indices
  %   P  #BE list of parent indices into BE
  %   dQ  #BE by 4 list of relative rotations
  % Outputs:
  %   Q  #BE by 4 list of absolute rotations
  %   T  #BE by 2 list of absolute translations
  %
  % Example:
  %   [Q,T] = forward_kinematics(C,BE,P,mpose);
  %   % stacked transposed 3x4 Affine matrices
  %   A = reshape(cat(2,permute(quat2mat(Q),[2 1 3]),permute(T,[2 3 1])),3,[])';
  %   n = size(V,1);
  %   M = reshape(bsxfun(@times,[V ones(n,1)],permute(W,[1 3 2])),n,[]);
  %   % Linear blend skinning:
  %   U = M*A

  function fk_helper(b)
    %function w = quatrotate(q,v)
    %  vv = [v 0];
    %  w = quatmultiply(quatmultiply(q,v),quatconj(q));
    %  w = w(1:3);
    %end
    if ~computed(b)
      p = P(b);
      if p < 1
        % base case for roots
        Q(b,:) = dQ(b,:);
        r = C(BE(b,1),:);
        T(b,:) = r - quatrotate(dQ(b,:),r);
      else
        % otherwise first compute parent's
        fk_helper(p);
        Q(b,:) = quatmultiply(Q(p,:),dQ(b,:));
        r = C(BE(b,1),:);
        T(b,:) = T(p,:) - quatrotate(Q(b,:),r) + ...
          quatrotate(Q(p,:),r);
      end
      computed(b) = true;
    end
  end

  % number of bones
  m = size(BE,1);
  assert(m == size(dQ,1));
  assert(m == size(P,1));

  % Whether we've already computed FK for this bone: dynamic programming
  computed = false(m,1);
  Q = zeros(m,4);
  T = zeros(m,3);
  for b = 1:m
    fk_helper(b);
  end

end


text:
% Extract the nearest exemplar from 3D pose library
% pool_2d: the 2D pose library for matching
% pool_3d: the 3D pose library
% pred: the prediction 2D pose
function [j_p] = NN_pose(pool_2d,pool_3d,pred)
    temp = zeros(1,28);
    temp_root = 0.5 * (pred(9,:) + pred(12,:));
    for a = 1:14
        temp(1,a*2-1) = pred(a,1) - temp_root(1);
        temp(1,a*2) = pred(a,2) - temp_root(2);
        y_c(a) = temp(1,2*a);
    end
    scale = max(y_c)-min(y_c);
    temp = temp / scale;

    [~,m_idx] = min(pdist2(temp, pool_2d));
    j_p = pool_3d(m_idx,:);
    j_p = reshape(j_p,3,14);
    j_p = double(j_p');
end

author: flyawaychase


text:
function [ fea, out ] = ex_heattransfer1( varargin )
%EX_HEATTRANSFER1 2D ceramic strip with radiation and convection.
%
%   [ FEA, OUT ] = EX_HEATTRANSFER1( VARARGIN ) 2D heat transfer of a ceramic strip with
%   both radiation and convection on the top boundary.
%
%         _                q = h*(T_inf-T) + epsilon*sigma*(T_inf^4-T^4)
%         ^            +------------------+
%         |            |                  |
%         |            |                  |
%       0.01m  T=900C  |                  |  T=900C
%         |            |                  |
%         |            |                  |
%         v            +------------------+
%                           dt/dn = 0
%                      |<---- 0.02m ----->|
%
%   The ceramic has a thermal conductivity of 3 W/mC and the sides are fixed
%   at a temperature of 900C while the bottom boundary is insulated. The surrounding
%   temperature is 50C. The top boundary is exposed to both natural convection (with
%   a film coefficient h=50W/m^2K) and radiation (with emissivity epsilon=0.7
%   and the Stefan-Boltzmann 5.669e-8 W/m^2K^4). The solution is sought at three
%   points along the vertical symmetry line.
%
%   Reference:
%
%      [1] Holman, J. P., Heat Transfer, Fifth Edition, New York: McGraw-Hill,
%          1981, page 96, Example 3-8.
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       hmax        scalar {0.001}         Grid cell size
%       igrid       scalar {0}/1/2         Cell type (0=quadrilaterals, 1=triangles,
%       solver      string fenics/{}       Use FEniCS or default solver
%       ischeme     scalar {0}             Time stepping scheme (0 = stationary)
%       sfun        string {sflag1}        Finite element shape function
%       iplot       scalar {1}/0           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { 'hmax',     0.001;
            'igrid',    0;
            'solver',   '';
            'ischeme',  0;
            'sfun',     'sflag1';
            'iplot',    1;
            'tol',      0.01;
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
if( opt.ischeme==2 && ~got.tol )
  opt.tol = 0.05;
end

% Geometry definition.
gobj = gobj_rectangle( 0, 0.02, 0, 0.01 );
fea.geom.objects = { gobj };

% Grid generation.
switch opt.igrid
  case 0
    fea.grid = rectgrid( round(0.02/opt.hmax), round(0.01/opt.hmax), [0 0.02;0 0.01] );
  case 1
    fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', opt.fid );
  case 2
    fea.grid = rectgrid( round(0.02/opt.hmax), round(0.01/opt.hmax), [0 0.02;0 0.01] );
    fea.grid = quad2tri( fea.grid, 1 );
end

% Problem definition.
fea.sdim  = { 'x', 'y' };             % Space coordinate name.
fea = addphys( fea, @heattransfer );  % Add heat transfer physics mode.
fea.phys.ht.sfun = { opt.sfun };      % Set shape function.

% Equation coefficients.
fea.phys.ht.eqn.coef{3,end} = 3;      % Thermal conductivity.

% Boundary conditions.
fea.phys.ht.bdr.sel = [3 1 4 1];
fea.phys.ht.bdr.coef{1,end}   = { [] 900+273 [] 900+273 };
fea.phys.ht.bdr.coef{4,end}{3}{2} = 50;
fea.phys.ht.bdr.coef{4,end}{3}{3} = 50+273;
fea.phys.ht.bdr.coef{4,end}{3}{4} = 0.7*5.669e-8;
fea.phys.ht.bdr.coef{4,end}{3}{5} = 50+273;

% Parse physics modes and problem struct.
fea = parsephys(fea);
fea = parseprob(fea);

% Compute solution.
if( strcmp(opt.solver,'fenics') )
  fea = fenics( fea, 'fid', opt.fid, ...
                'tstep', 0.1, 'tmax', 1, 'ischeme', opt.ischeme );
else
  if( opt.ischeme<=0 )
    fea.sol.u = solvestat( fea, 'fid', opt.fid, 'init', {'T0_ht'} );
  else
    [fea.sol.u,fea.sol.t] = solvetime( fea, 'fid', opt.fid, 'init', {'T0_ht'}, ...
                                       'tstep', 0.1, 'tmax', 1, 'ischeme', opt.ischeme );
  end
end

% Postprocessing.
if( opt.iplot>0 )
  postplot( fea, 'surfexpr', 'T', 'isoexpr', 'T' )
  title('Temperature, T')
end

% Error checking.
T2_sol = evalexpr( 'T', [0.01;0.01], fea );
T2_ref = 984;
T5_sol = evalexpr( 'T', [0.01;0.005], fea );
T5_ref = 1064;
T8_sol = evalexpr( 'T', [0.01;0], fea );
T8_ref = 1088;
out.err  = abs([T2_sol-T2_ref T5_sol-T5_ref T8_sol-T8_ref])./[T2_ref T5_ref T8_ref];
out.pass = all(out.err<opt.tol);

if( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function y = cvxaff( x )
global cvx___
y = subsref( cvx___.problems( x.problem_ ).duals, x.name_ );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function elli = fitEllipse(varargin)
%FITELLIPSE Fit an ellipse to a set of 2D points.
%
%   ELLI = fitEllipse(PTS)
%
%   Example
%     elli = [50 40  30 10 20];
%     pts = ellipseToPolygon(elli, 60) + randn(60,2);
%     figure; hold on; drawPoint(pts, 'ko');
%     axis equal; axis([0 100 0 100]);
%     ellFit = fitEllipse(pts);
%     drawEllipse(ellFit, 'b')
%
%   This is a rewrite of an original function from the authors cited below.
%   Changes from original submission include:
%   * convert angle of result ellipse to degrees (to comply with MatGeom
%       convention)
%   * update comments
%
% Authors:
% Andrew Fitzgibbon, Maurizio Pilu, Bob Fisher
% Reference: "Direct Least Squares Fitting of Ellipses", IEEE T-PAMI, 1999
%
% https://fr.mathworks.com/matlabcentral/fileexchange/3215-fit_ellipse
%
%  @Article{Fitzgibbon99,
%   author = "Fitzgibbon, A.~W.and Pilu, M. and Fisher, R.~B.",
%   title = "Direct least-squares fitting of ellipses",
%   journal = pami,
%   year = 1999,
%   volume = 21,
%   number = 5,
%   month = may,
%   pages = "476--480"
%  }
%
%   See also
%     geom2d, ellipses2d, createEllipse, equivalentEllipse
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2022-07-16, using Matlab 9.12.0.1884302 (R2022a)
% Copyright 2022 INRAE - BIA Research Unit - BIBS Platform (Nantes)

%% Process input arguments

if nargin==1
    var = varargin{1};
    X = var(:,1);
    Y = var(:,2);
else
    X = varargin{1};
    Y = varargin{2};
end

%% Normalize data

% recenter and reduce range
mx = mean(X);
my = mean(Y);
sx = (max(X) - min(X)) / 2;
sy = (max(Y) - min(Y)) / 2;
x = (X-mx) / sx;
y = (Y-my) / sy;

% Force to column vectors
x = x(:);
y = y(:);

%% Main processing

% Build design matrix
D = [ x.*x  x.*y  y.*y  x  y  ones(size(x)) ];

% Build scatter matrix
S = D' * D;

% Build 6x6 constraint matrix
C(6,6) = 0;
C(1,3) = -2;
C(2,2) = 1;
C(3,1) = -2;

% Solve eigensystem

% Break into blocks
tmpA = S(1:3,1:3);
tmpB = S(1:3,4:6);
tmpC = S(4:6,4:6);
tmpD = C(1:3,1:3);
tmpE = tmpC \ tmpB';

[evec_x, eval_x] = eig(tmpD \ (tmpA - tmpB*tmpE));

% Find the positive (as det(tmpD) < 0) eigenvalue
I = real(diag(eval_x)) < 1e-8 & ~isinf(diag(eval_x));

% Extract eigenvector corresponding to negative eigenvalue
A = real(evec_x(:,I));

% Recover the bottom half...
evec_y = -tmpE * A;
A = [A; evec_y];

% re-calibrate
par = [
    A(1)*sy*sy,   ...
    A(2)*sx*sy,   ...
    A(3)*sx*sx,   ...
    -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy,   ...
    -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy,   ...
    A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my   ...
    - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my   ...
    + A(6)*sx*sx*sy*sy   ...
    ]';

%% Identify parameters

% rotation angle
theta = 0.5 * atan2(par(2), par(1) - par(3));

% pre-comptute trigonometrics
cost = cos(theta);
sint = sin(theta);
sin2 = sint .* sint;
cos2 = cost .* cost;
cos_sin = sint .* cost;

%
Ao =  par(6);
Au =  par(4) .* cost + par(5) .* sint;
Av = -par(4) .* sint + par(5) .* cost;
Auu = par(1) .* cos2 + par(3) .* sin2 + par(2) .* cos_sin;
Avv = par(1) .* sin2 + par(3) .* cos2 - par(2) .* cos_sin;

% ROTATED = [Ao Au Av Auu Avv]

tuCentre = - Au./(2.*Auu);
tvCentre = - Av./(2.*Avv);
wCentre = Ao - Auu.*tuCentre.*tuCentre - Avv.*tvCentre.*tvCentre;

uCentre = tuCentre .* cost - tvCentre .* sint;
vCentre = tuCentre .* sint + tvCentre .* cost;

Ru = -wCentre ./ Auu;
Rv = -wCentre ./ Avv;

Ru = sqrt(abs(Ru)) .* sign(Ru);
Rv = sqrt(abs(Rv)) .* sign(Rv);

elli = [uCentre, vCentre, Ru, Rv, rad2deg(theta)];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
classdef FigureRendererEnum < matlab.mixin.SetGet
    %FigureRendererEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        OpenGL('OpenGL','opengl');
        Painters('Painters','painters');
    end

    properties
        name char = '';
        renderer char = '';
    end

    methods
        function obj = FigureRendererEnum(name,renderer)
            obj.name = name;
            obj.renderer = renderer;
        end
    end

    methods(Static)
        function listBoxStr = getListboxStr()
            m = enumeration('FigureRendererEnum');
            listBoxStr = {m.name};
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('FigureRendererEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('FigureRendererEnum');
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function a = reshape(a,new_size,old_modes)
%RESHAPE Reshape sparse tensor.
%
%  RESHAPE(X,SIZ) reshapes the sparse tensor to the given size. PROD(SIZ)
%  must be the same as PROD(SIZE(X)).
%
%  RESHAPE(X,SIZ,MODES) reshapes only the specifies modes and appends the
%  new reshaped modes to the end of the indices.
%
%  See also SPTENSOR, SPTENSOR/PERMUTE, RESHAPE.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if ~exist('old_modes','var')
    old_modes = 1:ndims(a);
    keep_modes = [];
else
    keep_modes = setdiff(1:ndims(a),old_modes);
end
old_size = a.size(old_modes);
keep_size = a.size(keep_modes);

if prod(new_size) ~= prod(old_size)
    error('prod(SIZ) must be the same size of prod(SIZE(X,MODES))');
end

inds = tt_sub2ind(old_size,a.subs(:,old_modes));
new_subs = tt_ind2sub(new_size,inds);

a.size = [keep_size new_size];
a.subs = [a.subs(:,keep_modes) new_subs];

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function ordered_textons = select_diverse_textons(tsim, nt)

next_t = zeros(nt, 1);
vals = zeros(nt, 1);
indices = (1:size(tsim, 1));
for t = 1:nt
    total_dsim = mean(tsim, 1);
    [vals(t), index] = max(total_dsim);
    next_t(t) = indices(index);
    tsim = tsim([(1:index-1) (index+1:end)], [(1:index-1) (index+1:end)]);
    indices = indices([(1:index-1) (index+1:end)]);
end

ordered_textons = next_t;

author: zouchuhang

repo: LayoutNet


text:
%% Demo: White balancing all images in a directory
%
% Copyright (c) 2018-present, Mahmoud Afifi
% York University, Canada
% mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
% All rights reserved.
%
% Please cite the following work if this program is used:
% Mahmoud Afifi, Brian Price, Scott Cohen, and Michael S. Brown,
% "When color constancy goes wrong: Correcting improperly white-balanced
% images", CVPR 2019.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%

%% input and options
in_image_dir = fullfile('..', 'example_images');
out_image_dir = fullfile('..', 'example_images_WB');
device = 'cpu'; %'cpu' or 'gpu'
gamut_mapping = 2; % use 1 for scaling, 2 for clipping (our paper's results
% reported using clipping). If the image is over-saturated, scaling is
% recommended.
upgraded_model = 0; % use 1 to load our new model that is upgraded with new
% training examples.

%%
switch lower(device)
    case 'cpu'
        if upgraded_model == 1
            load(fullfile('models','WB_model+.mat'));
        elseif upgraded_model == 0
            load(fullfile('models','WB_model.mat'));
        else
            error('Wrong upgraded_model value; please use 0 or 1');
        end
    case 'gpu'
        try
            gpuDevice();
        catch
            error('Cannot find a GPU device');
        end
        if upgraded_model == 1
            load(fullfile('models','WB_model+_gpu.mat'));
        elseif upgraded_model == 0
            load(fullfile('models','WB_model_gpu.mat'));
        else
            error('Wrong upgraded_model value; please use 0 or 1');
        end
    otherwise
        error('Wrong device; please use ''gpu'' or ''cpu''')
end
model.gamut_mapping = gamut_mapping;
if exist(out_image_dir,'dir')==0
    mkdir(out_image_dir);
end
imds = imageDatastore(in_image_dir);
files = imds.Files;
for f = 1 : length(files)
    fprintf('Processing image: %s\n',files{f});
    infileName = files{f};
    [~,name,ext] = fileparts(infileName);
    outfileName = fullfile(out_image_dir, [name, '_', 'WB', ext]);
    I_in = imread(infileName);
    I_corr = model.correctImage(I_in);
    if strcmpi(device,'gpu')
        imwrite(gather(I_corr),outfileName);
    else
        imwrite(I_corr,outfileName);
    end
    disp('Done!');
end

author: mahmoudnafifi


text:
function [ihs, x_mu, y_mu, sig] = plotreswh2(res, peval, dpixc, p, savethis, tilesonly, showbg, showpoints, shownumbers, showimage)
% function [ihs, x_mu, y_mu, sig] = plotreswh(res, peval, dpixc, p, savethis, tilesonly, showbg, showpoints, shownumbers, showimage)
% This is just modified version of plotreswh for more recent simulations
% (different in the way the scatter plot is made... +-1)
% showbg = 0;
% showpoints = 1;
% savethis = 0;
% tilesonly = 1;
if ~(showbg == 1)
    res.w = res.w(:,1:end-1);
    res.h = res.h(1:end-1,:);
end

if ~exist('shownumbers','var')
    shownumbers = 0;
end
if ~exist('showimage', 'var')
    showimage =1; %plot the figures
end

x_mu = [];
y_mu = [];
sig = [];

% mh = mean (res.h,2);
% [mhs,ihs] = sort(mh,'descend');
% sm = sum(res.w.^2,1);
% [mhs, ihs] = sort(sm, 'descend');
[mhs, ihs] = sortcomponents(res.w);
% [mhs, ihs] = sortcomponents(res.h);
% ts = testimportance(reshape(dpixc, peval.nx*peval.ny, peval.nt), res.w, res.h);
% [mhs, ihs] = sort(ts, 'descend');

hsort=res.h(ihs,:);
wsort = res.w(:,ihs);
if showbg > 1
    hsort = hsort(1:showbg,:);
    wsort = wsort(:,1:showbg);
end

wr = reshape(wsort,peval.nx,peval.ny,size(wsort,2));

sw = size(wr,3);
a = ceil(sqrt(sw));
b = ceil(sw/a);
maxd = max(dpixc,[],3);
maxd = maxd/max(maxd(:));
mwr = mean(squeeze(max(max(wr,[],1),[],2)));
for ii=1:sw
    [x_mu(ii), y_mu(ii), sig(ii), differ(ii)] = fitgauss2d(wr(:,:,ii));
end

if showimage
    ca
    if ~tilesonly
        if savethis
            mkdir ('w')
        end
        for ii=1:sw
            %     dipshow(ii,wr(:,:,ii))
            %     colormap('jet')
            dipshow(wr(:,:,ii));
            switch showpoints
                case 0
                    if savethis; SaveImageFULL(['w/w_' num2str(ii)]); end
                case 1
                    hold on
                    scatter(p.x_vec, p.y_vec,'r')
                    hold off
                    if savethis; SaveImageFULL(['w/w2_' num2str(ii)]); end
                case 2
                    %                 mwr = max(max(wr(:,:,ii)));
                    hold on
                    contour(mwr*maxd)
                    hold off
                    if savethis; SaveImageFULL(['w/w2_' num2str(ii)]); end
                case 3
                    %                 mwr = max(max(wr(:,:,ii)));
                    %                 [x_mu(ii), y_mu(ii), sig(ii), differ(ii)] = fitgauss2d(wr(:,:,ii));
                    hold on
                    if ~isempty(p)
                        %                     scatter(p.x_vec+.5, p.y_vec+.5,'r')
                        scatter(p.x_vec, p.y_vec,'r')
                    end
                    scatter(x_mu(ii)-1, y_mu(ii)-1,sig(ii)*100, 'rx')
                    hold off
                    if savethis; SaveImageFULL(['w/w2_' num2str(ii)]); end
            end
            %     figure(sw+1)
            %     subplot(a,b,ii)
            %     imagesc(wr(:,:,ihs(ii)))
            %     set(gca,'dataaspectratio',[1 1 1])
            %     set(gca,'xtick',[])
            %     set(gca,'ytick',[])
            %     xlabel(num2str(ii))
        end
    end

    figure ('name','w tiled');
    imstiled(wr,[],'gray',1)
    suptitle('w - ordered')
    if savethis; SaveImageFULL('w_tiled'); end
    % if showpoints
    %     for ii=1:size(res.w,2)-1
    %         dipshow(ii,wr(:,:,ihs(ii+1)))
    %         hold on
    %         scatter(p.x_vec-.5, p.y_vec-.5,'r')
    %         hold off
    %         SaveImageFULL(['w/w2_' num2str(ii)])
    %     end
    % end

    figure ('name','H');
    imagesc(hsort)
    if ~showbg; ylim([0.5 peval.ncomp+0.5]); end
    xlabel('slice # (time)')
    ylabel('component')
    title('H')
    if savethis; SaveImageFULL(['Himage']); end

    if tilesonly > 1
        cc = corrcoef(hsort');
        figure ('name','corrcoef(H'')');
        ims(cc)
        colorbar
        title('corrcoef(H'')')
        if savethis; SaveImageFULL('Hcorrel'); end

        figure ('name','abs(corrcoef(H''))');
        imagesc(abs(cc))
        set(gca, 'DataAspectRatio',[1 1 1])
        colorbar
        title('abs(corrcoef(H''))')
        if savethis; SaveImageFULL('Hcorrelabs'); end

        figure ('name','mean(H,2)');
        meanh = mean(hsort,2);
        bar(meanh)
        xlabel('component')
        ylabel('mean (H,2)')
        grid on
        title('sum(H,2)')
        % if ~showbg; xlim([1.5 peval.ncomp+1.5]); end
        if savethis; SaveImageFULL('Hmean_bar'); end

        figure ('name','corrcoef(W)');
        ims(corrcoef(wsort))
        colorbar
        title('corrcoef(W)')
        if savethis; SaveImageFULL('Wcorrel'); end
    end
    figure ('name','sum');
    dipshow(sum(dpixc,3),'gray')
    hold on
    if ~isempty(p)
        %     scatter(p.x_vec, p.y_vec,'r')
        scatter(p.x_vec, p.y_vec,'b')
    end
    scatter(x_mu-1, y_mu-1,sig.^2*100, 'rx')
    if shownumbers
        text(x_mu-1+.2, y_mu-1, num2cell([1:sw]) )
    end
    % scatter(x_mu, y_mu,sig.^2*100, [1:sw])
    if savethis; SaveImageFULL('Sum'); end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [u,v]=pow2cep(m,c,mode)
%CEP2POW convert cepstral means and variances to the power domain
% Inputs:
%    m: vector giving means in the power domain
%    c: covariance matrix in the power domain
% mode: 'c'  pow=exp(irdct(cep))   [default]
%       'f'  pow=exp(rsfft(cep)/n)  [fft length even]
%       'fo' pow=exp(rsfft(cep)/n)  [fft length odd]
%       'i'  pow=exp(cep)           [ no transformation ]
%
% Outputs:
%    u: row vector giving the cepstral means with u(1) the 0'th cepstral coefficient
%    v: cepstral covariance matrix

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: pow2cep.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3 mode='c'; end
if min(size(c))==1
   v=diag(c);
end
m=m(:)';        % force to be a row vector
q=log(1+c./(m'*m));
p=log(m)-0.5*diag(q)';
if any(mode=='f')
   n=2*length(m)-2;
   if any(mode=='o')
      n=n+1;
   end
   u=rsfft(p,n);
   v=rsfft(rsfft(q,n)',n);
elseif any(mode=='i')
    u=p;
    v=q;
else
   u=rdct(p);
   v=rdct(rdct(q)');
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
 function ob = ctranspose(ob)
%function ob = ctranspose(ob)
% "ctranspose" method for Fatrix class

ob.is_transpose = ~ob.is_transpose;

% [ob.index1 ob.index2] = deal(ob.index2, ob.index1);
ob.dim = fliplr(ob.dim);


text:
classdef one_dim_denoise_problem

    properties
        name;
        dim;
        samples;
        lambda;
        Lap;
        LaptLap;
        b;
        n;
        d;
        prox_flag;
    end

    methods
        function obj = one_dim_denoise_problem(b, lambda, varargin)

            obj.name = 'one_dim_denoise_problem';
            obj.b = b;
            obj.lambda = lambda;
            obj.samples = length(b);

            obj.n = obj.samples;
            obj.d = obj.samples;
            obj.dim = 1;
            obj.samples = obj.n;

            obj.Lap = zeros(obj.d-1, obj.d);
            for i = 1 : obj.d-1
                obj.Lap(i,i) = 1;
                obj.Lap(i,i+1) = -1;
            end
            obj.LaptLap = obj.Lap' * obj.Lap;

            obj.prox_flag = false;
        end

%         function v = prox_denoise(w, t)
%
%             v = soft_thresh(w, t * obj.lambda);
%
%         end

        function f = cost(obj, w)

            f = obj.loss(w) + obj.lambda * obj.reg(w);

        end

        function l = loss(obj, w)

            l = 1/2 * norm(obj.b-w)^2;
        end

        function r = reg(obj, w)

            % Lw = obj.Lap * w;
            % r = (Lw'*Lw)/2;
            % r = w'L'Lw/2
            r = w' * obj.LaptLap * w / 2;

        end

        function r = differentiable_reg(obj, w)

            r = w' * obj.LaptLap * w / 2;

        end

        function r = residual(obj, w)
            r = - (obj.b - w);
        end

        function g = full_grad(obj, w)
            g = - (obj.b - w) + obj.lambda * obj.reg_grad(w);
        end

        function g = grad(obj, w, indices)
            g = obj.full_grad(w);
        end

        function rg = reg_grad(obj, w)

            rg = obj.Lap' * obj.Lap * w;

         end

        function h = hess(obj, w, indices)
            error('Not implemted yet.');
        end

        function h = full_hess(obj, w)
            error('Not implemted yet.');
        end

        function hv = hess_vec(obj, w, v, indices)
            error('Not implemted yet.');
        end

        function w_opt = calc_solution(obj, w_init, options)

            w_opt = (eye(obj.d) + obj.lambda * obj.LaptLap) \ obj.b;

        end

    end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function f = fx3 ( n, x )

%*****************************************************************************80
%
%% FX3 is the third 1D example.
%
%  Discussion:
%
%    The function should be plotted over [-1.0,+1.0].
%
%    Internally, this range is mapped to [-3.0,+3.0].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Rick Archibald, Anne Gelb, Jungho Yoon,
%    Polynomial fitting for edge detection in irregularly sampled signals
%    and images,
%    SIAM Journal on Numerical Analysis,
%    Volume 43, Number 1, 2006, pages 259-279.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), the arguments.
%
%    Output, real F(N,1), the function values.
%

%
%  Destroy all row vectors!
%
  x = x ( : );
%
%  Map from the convenient range [-1,+1] to the physical range [-3,+3].
%
  x = ( ( 1.0 - x ) * ( -3.0 )   ...
      + ( 1.0 + x ) * ( +3.0 ) ) ...
      /   2.0;

  f = zeros ( n, 1 );

  i = find ( -2.0 <= x & x <= -1.0 );
  f(i) = 1.0;

  j = find ( -0.5 <= x & x <= 0.5 );
  f(j) = 0.5 + 4.0 * ( x(j) + 0.5 ).^2;

  k = find ( 1.25 <= x & 3.0 * x <= 7.0 );
  f(k) = 3.0 * ( 2.0 - x(k) );

  return
end


text:
function [] = anim8_DIC_image_3Dmeasure_points_2n(ImSet,DIC2DpairResults,DIC3DpairResults,faceMeasureString,varargin)
%% function for plotting 2D-DIC results imported from Ncorr in step 2
% called inside plotNcorrPairResults
% plotting the images chosen for stereo DIC (2 views) with the
% correlated points results plotted on top, colored as their correlation
% coefficient.
% on the left side the images from the reference camera (reference image and current images), and on the right side the
% images from the deformed camera
% requirements: GIBBON toolbox
%
%
% INPUT:
% * IMset - a 2nX1 cell array containing 2n grayscale images. The first n
% images are from camera A (the "reference" camera), and the last n images
% are from camera B (the "deformed" camera). The first image in the set is
% considered as the reference image, on which the reference grid of points
% is defined, and all the correlated points and consequent displacements
% and strains, are relative to this image.
% * DIC_2Dpair_results - containig the correlated points, correlation
% coefficients, faces..
% * optional: CorCoeffCutOff - - maximal correlation coefficient to plot
% points
% * optional: CorCoeffDispMax - maximal correlation coefficient in colorbar

%%
DIC2DpairResultsL=DIC2DpairResults{1};
DIC2DpairResultsR=DIC2DpairResults{2};
DIC3DpairResultsL=DIC3DpairResults{1};
DIC3DpairResultsR=DIC3DpairResults{2};

nImages=DIC2DpairResultsL.nImages;
nCam=DIC2DpairResultsL.nCamDef;

PointsL=DIC2DpairResultsL.Points(nImages+1:end);
PointsR=DIC2DpairResultsR.Points(1:nImages);
CorCoeffVecL=DIC2DpairResultsL.CorCoeffVec(nImages+1:end);
CorCoeffVecR=DIC2DpairResultsR.CorCoeffVec(1:nImages);

switch nargin
    case 4 % in case no results were entered
        optStruct=struct;
    case 5
        optStruct=varargin{1};
    otherwise
        error('wrong number of input arguments');
end

%% cut out point with large correlation coefficient
if ~isfield(optStruct,'CorCoeffCutOff')
    CorCoeffCutOff=max([cell2mat(CorCoeffVecL); cell2mat(CorCoeffVecR)]);
else
    CorCoeffCutOff=optStruct.CorCoeffCutOff;
end

for ii=1:nImages
    CorCoeffVecL{ii}(CorCoeffVecL{ii}>CorCoeffCutOff)=NaN;
    CorCoeffVecR{ii}(CorCoeffVecR{ii}>CorCoeffCutOff)=NaN;
end

%%
switch faceMeasureString
    case {'DispMgn'}
        PCL=DIC3DpairResultsL.Disp.DispMgn;
        PCR=DIC3DpairResultsR.Disp.DispMgn;
        cMap='jet';
     case {'DispX'}
         for ii=1:nImages
             PCL{ii,1}=DIC3DpairResultsL.Disp.DispVec{ii}(:,1);
             PCR{ii,1}=DIC3DpairResultsR.Disp.DispVec{ii}(:,1);
         end
         cMap='jet';
    case {'DispY'}
        for ii=1:nImages
            PCL{ii,1}=DIC3DpairResultsL.Disp.DispVec{ii}(:,2);
            PCR{ii,1}=DIC3DpairResultsR.Disp.DispVec{ii}(:,2);
        end
        cMap='jet';
    case {'DispZ'}
        for ii=1:nImages
            PCL{ii,1}=DIC3DpairResultsL.Disp.DispVec{ii}(:,3);
            PCR{ii,1}=DIC3DpairResultsR.Disp.DispVec{ii}(:,3);
        end
        cMap='jet';
    otherwise
        error('unexpected face measure string. plots not created');
end

%%

if ~isfield(optStruct,'PClimits')
    PCmin=0;
    PCmax=0;
    for ii=1:nImages
        PCmin=min([min(PCL{ii}(~isnan(CorCoeffVecL{ii}))) min(PCL{ii}(~isnan(CorCoeffVecL{ii}))) PCmin]);
        PCmax=max([max(PCR{ii}(~isnan(CorCoeffVecR{ii}))) max(PCR{ii}(~isnan(CorCoeffVecR{ii}))) PCmax]);
    end
    PClimits=[PCmin PCmax];
else
    PClimits=optStruct.PClimits;
end

%%
hf=figure; hold all;
hf.Units='normalized'; hf.OuterPosition=[.05 .05 .9 .9]; hf.Units='pixels';

ii=1;
hp1=imagesc(repmat(ImSet{ii},1,1,3)); hold on; axis ij
PL=PointsL{ii}(~isnan(CorCoeffVecL{ii}),:);
PR=PointsR{ii}(~isnan(CorCoeffVecR{ii}),:);
hp2=scatter(PL(:,1),PL(:,2),6,PCL{ii}(~isnan(CorCoeffVecL{ii})),'+');
hp3=scatter(PR(:,1),PR(:,2),6,PCR{ii}(~isnan(CorCoeffVecR{ii})),'+');
pbaspect([size(ImSet{ii},2) size(ImSet{ii},1) 1])
hs1=title(['Cam ' num2str(nCam) ' frame ' num2str(1)]);
colormap(cMap);
hc1=colorbar;
caxis(PClimits)
title(hc1, faceMeasureString);
hc1.FontSize=16;
axis off
drawnow

%Create the time vector
animStruct.Time=linspace(0,1,nImages);

for ii=1:nImages
    xNowL=PointsL{ii}(~isnan(CorCoeffVecL{ii}),1);
    yNowL=PointsL{ii}(~isnan(CorCoeffVecL{ii}),2);
    xNowR=PointsR{ii}(~isnan(CorCoeffVecR{ii}),1);
    yNowR=PointsR{ii}(~isnan(CorCoeffVecR{ii}),2);

    cNowL=PCL{ii};
    cNowR=PCR{ii};

    TitleNow=['Cam ' num2str(nCam) ' frame ' num2str(ii)];

    %Set entries in animation structure
    animStruct.Handles{ii}=[hp1,hp2,hp2,hp2,hp3,hp3,hp3,hs1]; %Handles of objects to animate
    animStruct.Props{ii}={'CData','XData','YData','CData','XData','YData','CData','String'}; %Properties of objects to animate
    animStruct.Set{ii}={repmat(ImSet{ii},1,1,3),xNowL,yNowL,cNowL,xNowR,yNowR,cNowR,TitleNow}; %Property values for to set in order to animate

end

anim8(hf,animStruct);

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function KLdiv = KL_divergence2(bnetP, bnetQ)
% KL_DIVERGENCE2 computes the Kullback-Leibler divergence between two BNET distributions
% KLdiv = KL_divergence2(bnetP, bnetQ)
%
% Output :
%   div = sum_x  P(x).log(P(x)/Q(x))
%
% Rem :
%   This version is optimized for memory use, but quite slow !!!
%     ==> if you have no memory problem, use kl_divergence instead
%
%   ONLY FOR TABULAR NODES
%   Make sure that you have done the params learning.
%
%   V1.1 : 8 oct 2004 (Ph. Leray - philippe.leray@univ-nantes.fr)

N = size(bnetP.dag,1);
N2 = size(bnetQ.dag,1);
ns= bnetP.node_sizes;
ns2= bnetQ.node_sizes;
if N~=N2, error('size of dags must be the same'), end
if ns~=ns2, error('node sizes of dags must be the same'), end
tiny = exp(-700);
KLdiv=0;

for i=1:prod(ns),
  inst = ind2subv(ns, i); % i'th instantiation
  Px=1; Qx=1;
  for i=1:N,
    ps = parents(bnetP.dag, i);
    e = bnetP.equiv_class(i);
    [tmp Pxi] = prob_node(bnetP.CPD{e}, inst(i), inst(ps)');
    Px=Px*Pxi;
    ps = parents(bnetQ.dag, i);
    e = bnetQ.equiv_class(i);
    [tmp Qxi] = prob_node(bnetQ.CPD{e}, inst(i), inst(ps)');
    Qx=Qx*Qxi;
  end
    Px = Px + (Px==0)*tiny; % replace 0s by tiny
    Qx = Qx + (Qx==0)*tiny; % replace 0s by tiny
  KLdiv = KLdiv + Px*log(Px/Qx);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef BundleAdjust < matlab.mixin.Copyable

    properties
        % Cameras and landmarks are represented by graph nodes
        %   coord(camera) -> camera pose, 6 vector
        %   vdata(camera) -> index into state vector or [] if a fixed camera
        %   coord(landmark) -> landmark position, first 3 elements of 6 vector
        %   vdata(landmark) -> index into state vector
        %
        % An edge from camera to a landmark holds the observed projection as a
        % property
        %   edata(edge) -> uv*

        g         % PGraph object describing visibility

        camera    % camera projection model

        cameras   % list of graph nodes corresponding to cameras
        points    % list of graph nodes corresponding to landmarks

        fixedcam     % logical array of cameras that are fixed
        fixedpoint   % logical array of landmarks that are fixed
    end

    properties (Dependent=true)
        ncams       % number of cameras
        nlandmarks  % number of landmark points
        ndim        % size of the linear system
    end

    methods

        function ba = BundleAdjust(camera)
            %BundleAdjust.BundleAdjust Bundle adjustment problem constructor
            %
            % Notes::
            % - A cameraModel file must exist on the path to evaluate the point projection
            %   and Jacobians.

            ba.camera = camera;  % stash the camera model
            ba.g = PGraph(6);    % initialize the graph, nodes have 6D coordinates
        end

        function n = get.ncams(ba)
            n = length(ba.cameras);
        end

        function n = get.nlandmarks(ba)
            n = length(ba.points);
        end

        function n = get.ndim(ba)
            n = 6*(ba.ncams - sum(ba.fixedcam)) + 3*(ba.nlandmarks - sum(ba.fixedpoint));
        end

        function n = cameraIdx(ba, i)
            n = 6*i - 5;
        end

        function n = landmarkIdx(ba, j)
            n = ba.ncams*6 + 3*j-2;
        end

        function n = cameraIdx2(ba, i)
            if ba.fixedcam(i)
                n = 0;
            else
                n = 6*(i - sum(ba.fixedcam(1:i))) - 5;
            end
        end

        function n = landmarkIdx2(ba, j)
            if ba.fixedpoint(j)
                n = 0;
            else
                n = 6*(ba.ncams - sum(ba.fixedcam)) + 3 * (j - sum(ba.fixedpoint(1:j))) - 2;
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % METHODS TO SUPPORT BUILDING PROBLEMS
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function vc = add_camera(ba, T, varargin)
            %BundleAdjust.add_camera Add camera to bundle adjustment problem
            %
            % VC = BA.add_camera(T, options) is the vertex id of a camera node added to
            % a bundle adjustment problem.  The camera has pose T which can be SE3 or a
            % vector (1x7) comprising translation and quaternion.
            %
            % VC = BA.add_camera() as above but the camera pose is the null pose in SE3.
            %
            % Options::
            % 'fixed'    This camera is fixed (anchored) and will not be adjusted in
            %            the optimization process.
            %
            % See also BundleAdjust.add_landmark, BundleAdjust.add_projection, PGraph.

            if nargin == 1
                T = SE3();
            end

            opt.fixed = false;
            opt = tb_optparse(opt, varargin);

            if isvec(T, 7)
                t = T(1:3); q = T(4:7);
            else
                [R,t] = tr2rt(T);
                q = UnitQuaternion(R).double;
                t = t';
            end

            if q(1) < 0
                q = -q;
            end
            x = [t q(2:4)];
            vc = ba.g.add_node(x);
            ba.cameras = [ba.cameras vc];
            ba.g.setvdata(vc, length(ba.cameras));

            ba.fixedcam = [ba.fixedcam opt.fixed];
        end

        function vp = add_landmark(ba, P, varargin)
            %BundleAdjust.add_landmark Add landmark to bundle adjustment problem
            %
            % VL = BA.add_landmark(P, options) is the vertex id of a landmark node added to
            % a bundle adjustment problem.  The landmark has position P (3x1).
            %
            % Options::
            % 'fixed'    This landmark is fixed (anchored) and will not be adjusted in
            %            the optimization process.
            %
            % See also BundleAdjust.add_camera, BundleAdjust.add_projection, PGraph.

            assert(isvec(P), 'P must be a 3-vector');

            opt.fixed = false;
            opt = tb_optparse(opt, varargin);

            x = [P(:); 0; 0; 0];
            vp = ba.g.add_node(x);
            ba.points = [ba.points vp];
            ba.g.setvdata(vp, length(ba.points));

            ba.fixedpoint = [ba.fixedpoint opt.fixed];

        end

        function e = add_projection(ba, c, p, uv)
            %BundleAdjust.add_projection Add camera to bundle adjustment problem
            %
            % EP = BA.add_projection(VC, VL, UV) is the edge id of an edge added to
            % a bundle adjustment problem.  it represents the observed projection UV (2x1) of the landmark node VL as
            % seen by the camera node VC.
            %
            % See also BundleAdjust.add_camera, BundleAdjust.add_landmark, PGraph.

            assert(isvec(uv, 2), 'uv must be a 2-vector');

            e = ba.g.add_edge(c, p);
            ba.g.setedata(e, uv(:));
        end

        function load_sba(ba, cameraFile, pointFile, calibFile)
            %BundleAdjust.load_sba Load bundle adjustment problem from files
            %
            % BA.load_sba(camfile, pointfile, calibfile) loads a bundle adjustment
            % problem from data files as distributed with the SBA package
            %
            % Example::
            %
            % To solve the 7-point bundle adjustment problem distributed with SBA:
            %
            %           ba = BundleAdjust();
            %           pth = 'sba-1.6/demo/'
            %           ba.load_sba([pth '7cams.txt'], [pth '7pts.txt'], [pth 'calib.txt']);
            %
            % Reference::
            % - Sparse Bundle Adjustment package by Manolis Lourakis
            %   http://users.ics.forth.gr/~lourakis/sba
            %
            % See also BundleAdjust.add_camera, BundleAdjust.add_landmark, BundleAdjust.add_projection, PGraph.

            % adopted from sba-1.6/matlab/eucsbademo.m

            % read in camera parameters
            [q1, q2, q3, q4, tx, ty, tz] = textread(cameraFile, '%f%f%f%f%f%f%f', 'commentstyle', 'shell');

            for i=1:length(q1)
                ba.add_camera( [tx(i) ty(i) tz(i) q1(i) q2(i) q3(i) q4(i)] );
            end

            % read points file line by line
            %
            % each line is a world point
            fid = fopen(pointFile);
            npts = 0;
            while ~feof(fid)
                line = fgets(fid);
                [A, count, errmsg, nextindex] = sscanf(line, '%f%f%f%f', [1, 4]); % read X, Y, Z, nframes
                if(size(A, 2)>0) % did we read anything?
                    npts=npts+1;

                    % create a node for this point
                    lm = ba.add_landmark(A(1:3));

                    % now find which cameras it was seen by
                    nframes = A(4);
                    for i=1:nframes % read "nframes" id, x, y triplets
                        [A, count, errmsg, j] = sscanf(line(nextindex:length(line)), '%f%f%f', [1, 3]); % read id, x, y
                        nextindex=nextindex+j; % skip the already read line prefix

                        % add a landmark projection
                        ba.add_projection( ba.cameras(A(1)+1), lm, A(2:3) );
                    end
                end
            end
            fclose(fid);

            % read in calibration parameters
            [a1, a2, a3] = textread(calibFile, '%f%f%f', 'commentstyle', 'shell');
            % a1(1)    a2(1)  a3(1)  a2(2)    a3(2)];
            % f/rho_u  skew   u0     f/rho_v  v0

            ba.camera.f = a1(1);
            ba.camera.rho = [1 a1(1)/a2(2)];
            ba.camera.pp = [a3(1) a3(2)];
        end

        function T = getcamera(ba, ii)
            assert(all(ii > 0), 'camera indices must be > 0');
            assert(all(ii <= ba.ncams), 'camera index out of range');
            if nargin == 1
                ii = [1:ba.ncams];
            end
            for k=1:length(ii)
                i = ii(k);
                node = ba.cameras(i);
                qv = ba.g.coord(node);
                T(k) = SE3(qv(1:3)') * UnitQuaternion.vec(qv(4:6)).SE3;
            end

        end

        function P = getlandmark(ba, jj)
            assert(all(jj > 0), 'landmark indices must be > 0');
            assert(all(jj <= ba.nlandmarks), 'landmark index out of range');
            if nargin == 1
                jj = [1:ba.nlandmarks];
            end
            for k=1:length(jj)
                j = jj(k);
                node = ba.points(j);
                qv = ba.g.coord(node);
                P(:,k) = qv(1:3);
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % METHODS TO SOLVE PROBLEMS
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function [ba2,e] = optimize(ba, varargin)
            %BundleAdjust.optimize  Optimize the solution
            %
            % BA2 = BA.optimize(

            opt.iterations = 1000;
            opt.animate = false;
            opt.retain = false;
            opt.lambda = 0.1;
            opt.lambdamin = 1e-8;
            opt.dxmin = 1e-8;
            opt.tol = 0.05;

            [opt,args] = tb_optparse(opt, varargin);

            %g2 = PGraph(pg.graph);  % deep copy

            if length(args) > 0 && isvec(args{1}, ba.ndim)
                X = args{1};
            else
                X = ba.getstate();
            end

            lambda = opt.lambda;

                        X0 = X;

            t0 = cputime();

            fprintf('Initial cost %g\n', ba.errors(X));

            for i=1:opt.iterations
                if opt.animate
                    if ~opt.retain
                        clf
                    end
                    g2.plot();
                    pause(0.5)
                end

                tic

                % solve for the step
                [dX,energy] = ba.solve(X, lambda);
                % update the state
                Xnew = ba.updatestate(X, dX);
                % compute new value of cost
                enew = ba.errors(Xnew);

                dt = toc;
                fprintf('  total cost %g (solved in %.2g sec)', enew, dt);

                % are we there yet?
                if enew < opt.tol
                    break
                end

                % have we stopped moving
                if  norm(dX) < opt.dxmin
                    break
                end

                % do the Levenberg-Marquadt thing, was it a good update?
                if enew < energy
                    % step is accepted
                    X = Xnew;
                    if lambda > opt.lambdamin
                        lambda = lambda/sqrt(2);
                    end
                    if opt.verbose
                        fprintf(' -- step accepted: lambda = %g', lambda);
                    end
                else
                    % step is rejected
                    lambda = lambda*4;
                    if opt.verbose
                        fprintf(' -- step rejected: lambda = %g', lambda);
                    end
                end

                fprintf('\n');
            end

            tf = cputime();
            fprintf('\n * %d iterations in %.1f seconds\n', i, tf-t0);
            fprintf(' * %.2f pixels RMS error\n', sqrt(enew/ba.g.ne));

            % copy the object and update it
            ba2 = copy(ba);
            ba2.g = copy(ba.g);
            ba2.setstate(X);

            if nargout > 1
                e = enew;
            end
        end

        function [deltax,err] = solve(ba, X, lambda)

            if nargin < 3
                lambda =0;
            end

            %ba.init_tables();

            % create the Hessian and error vector
            [H,b,e] = ba.build_linear_system(X);

            % add damping term to the diagonal
            for i=1:ba.ndim
                H(i,i) = H(i,i) + lambda;
            end

            % solve for the state update
            % - could replace this with the Schur complement trick
            deltax = H \ b;

            if nargout > 1
                err = e;
            end
        end

        % build the Hessian and measurement vector
        function [H,b,etotal,J] = build_linear_system(ba, X)

            % allocate storage
            H = sparse(ba.ndim, ba.ndim);
            b = zeros(ba.ndim,1);

            etotal = 0;

            % loop over cameras
            for i=1:ba.ncams

                k = ba.cameraIdx(i);
                x = X(k:k+5);

                t = x(1:3); qv = x(4:6);

                % loop over all points viewed from this camera
                for p=ba.g.edges( ba.cameras(i) )
                    v = ba.g.vertices(p);
                    j = ba.g.vdata(v(2));

                    k = ba.landmarkIdx(j);
                    x = X(k:k+2);

                    P = x';

                    uv = ba.g.edata(p);

                    % compute Jacobians and projection

                    [uvhat,JA,JB] = ba.camera.derivs(t, qv, P);

                    % compute reprojection error
                    e = uvhat - uv;
                    etotal = etotal + e'*e;

                    ii = ba.cameraIdx2(i); jj = ba.landmarkIdx2(j);

                    % compute the block components of H and b for this edge

                    if ~ba.fixedcam(i) & ~ba.fixedpoint(j)
                        % adjustable point and camera
                        H_ii= JA'*JA;
                        H_ij= JA'*JB;
                        H_jj= JB'*JB;

                        H(ii:ii+5,ii:ii+5) = H(ii:ii+5,ii:ii+5) + H_ii;
                        H(ii:ii+5,jj:jj+2) = H(ii:ii+5,jj:jj+2) + H_ij;
                        H(jj:jj+2,ii:ii+5) = H(jj:jj+2,ii:ii+5) + H_ij';
                        H(jj:jj+2,jj:jj+2) = H(jj:jj+2,jj:jj+2) + H_jj;

                        b_i = -JA'*e;
                        b_j = -JB'*e;

                        b(ii:ii+5) = b(ii:ii+5) + b_i;
                        b(jj:jj+2) = b(jj:jj+2) + b_j;

                    elseif ba.fixedcam(i) & ~ba.fixedpoint(j)
                        % fixed camera and adjustable point

                        H_jj= JB'*JB;

                        H(jj:jj+2,jj:jj+2) = H(jj:jj+2,jj:jj+2) + H_jj;

                        b_j = -JB'*e;

                        b(jj:jj+2) = b(jj:jj+2) + b_j;

                    elseif ~ba.fixedcam(i) & ba.fixedpoint(j)
                        % adjustable camera and fixed point

                        H_ii= JA'*JA;

                        H(ii:ii+5,ii:ii+5) = H(ii:ii+5,ii:ii+5) + H_ii;

                        b_i = -JA'*e;

                        b(ii:ii+5) = b(ii:ii+5) + b_i;
                    end

                end
            end

        end

        function spyH(ba, X)
            H = build_linear_system(ba, X);
            spy(H);
        end

        function X = getstate(ba)
            %BundleAdjust.getstate Get the state vector
            %
            % X = BA.getstate() is a row vector containing the state vector of the
            % problem: all camera poses followed by all landmark coordinates.
            %
            % Notes::
            % - The length of the vector is given by BA.ndim.

            X = [];

            i = 0;
            for c=ba.cameras  % step through camera nodes
                X = [X  ba.g.coord(c)'];
            end

            for p=ba.points  % step through landmark nodes
                P = ba.g.coord( p );
                X = [X  P(1:3)'];
            end
        end

        function setstate(ba, X)
            %BundleAdjust.getstate Get the state vector
            %
            % X = BA.getstate() is a row vector containing the state vector of the
            % problem: all camera poses followed by all landmark coordinates.
            %
            % Notes::
            % - The length of the vector is given by BA.ndim.

            for i=1:ba.ncams  % step through camera nodes
                if ~ba.fixedcam(i)
                    k = ba.cameraIdx(i);
                    c = ba.cameras(i);
                    ba.g.setcoord(c, X(k:k+5));
                end
            end

            for j=1:ba.nlandmarks
                if ~ba.fixedpoint(j)
                    k = ba.landmarkIdx(j);
                    lm = ba.points(j);
                    ba.g.setcoord(lm, [X(k:k+2) 0 0 0]);
                end
            end
        end

        function Xnew = updatestate(ba, X, dX)
            dX = dX(:)';

            % for each camera we need to compound the camera pose with the
            % incremental relative pose
            for i=1:ba.ncams
                k = ba.cameraIdx(i);
                if ba.fixedcam(i)
                    Xnew(k:k+5) = X(k:k+5);
                else
                    x = X(k:k+5);
                    t = x(1:3); qv = x(4:6);  % get current pose

                    k2 = ba.cameraIdx2(i);
                    dx = dX(k2:k2+5);
                    dt = dx(1:3); dqv = dx(4:6); % get incremental pose
                    tnew = t+dt;  % assume translation in old frame
                    % compound the quaternion vector rotations
                    % - function qvmul is symbolically generated by ba
                    %qvnew = qvmul(qv(1), qv(2), qv(3), dqv(1), dqv(2), dqv(3));
                    qvnew = UnitQuaternion.qvmul(qv, dqv);

                    Xnew(k:k+5) = [tnew qvnew];
                end
            end

            % for each landmark we add the increment to its position
            for j=1:ba.nlandmarks
                k = ba.landmarkIdx(j);
                x = X(k:k+2);
                if ba.fixedpoint(j)
                    Xnew(k:k+2) = x;
                else
                    k2 = ba.landmarkIdx2(j);
                    dx = dX(k2:k2+2);
                    Xnew(k:k+2) = x + dx;
                end
            end
        end

        % Compute total squared reprojection error
        function etotal = errors(ba, X)

            if nargin < 2
                X = ba.getstate();
            end
            r = ba.getresidual(X);

            etotal = sum(r(:));
        end

        function r = getresidual(ba, X)
            % this is the squared reprojection errors

            if nargin == 1
                X = ba.getstate();
            end

            % loop over cameras
            for i=1:ba.ncams

                k = ba.cameraIdx(i);
                x = X(k:k+5);

                t = x(1:3); qv = x(4:6);

                % loop over all points viewed from this camera
                for p=ba.g.edges( ba.cameras(i) )
                    v = ba.g.vertices(p);
                    j = ba.g.vdata(v(2));

                    k = ba.landmarkIdx(j);
                    x = X(k:k+2);

                    P = x';

                    uv = ba.g.edata(p);

                    uvhat = ba.camera.derivs(t, qv, P);

                    % compute reprojection error
                    e = uvhat - uv;
                    r(i,j) =  e'*e;
                end
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % SUPPORT METHODS
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function plot(ba, varargin)
            clf
            ba.g.plot('dims', 3, 'only', ba.points, 'EdgeColor', 0.8*[1 1 1], varargin{:});
            axis equal
            hold on

            colorOrder = get(gca, 'ColorOrder');
            for i=1:ba.ncams
                T = ba.getcamera(i);
                cam(i) = ba.camera.move(T);
                cidx = mod(i-1, numrows(colorOrder))+1;
                color = colorOrder(cidx,:);
                cam(i).plot_camera('scale', 0.2, 'color', color, 'persist')
            end
            hold off
            axis equal
            xlabel('X (m)')
            ylabel('Y (m)')
            zlabel('Z (m)')
            grid on
        end

        function display(ba)
            %BundleAdjust.display BundleAdjust parameters
            %
            % BA.display() displays the bundle adjustment parameters in compact format.
            %
            % Notes::
            % - This method is invoked implicitly at the command line when the result
            %   of an expression is a BundleAdjust object and the command has no trailing
            %   semicolon.
            %
            % See also BundleAdjust.char.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(ba) );
        end % display()

        function s = char(ba)
            %BundleAdjust.char Convert to string
            %
            % BA.char() is a string showing bundle adjustment parameters in compact format.
            %
            %
            % See also BundleAdjust.display.
            s = 'Bundle adjustment problem:';
            s = strvcat(s, sprintf('  %d cameras', ba.ncams));
            if any(ba.fixedcam)
                s = strvcat(s, sprintf('    locked cameras: %s', num2str(find(ba.fixedcam))));
            end

            s = strvcat(s, sprintf('  %d landmarks', ba.nlandmarks));

            s = strvcat(s, sprintf('  %d projections', ba.g.ne));

            s = strvcat(s, sprintf('  %d dimension linear problem', ba.ndim));
            c = ba.g.connectivity(ba.cameras);
            s = strvcat(s, sprintf('  landmarks per camera: min=%.1f, max=%.1f, avg=%.1f', min(c), max(c), mean(c)));
            c = ba.g.connectivity(ba.points);
            s = strvcat(s, sprintf('  cameras per landmark: min=%.1f, max=%.1f, avg=%.1f', min(c), max(c), mean(c)));
        end

    end % methods

end % classdef

author: petercorke

repo: machinevision-toolbox-matlab