Dear all,
Attached is the first version of gramian.m
It is not very efficient, any suggestions to improve speed is welcome.
I have tested functionalities described in the demo section. There are
a couple of options to parse, but the basics are there.
Looking forward to hear your comments.
--
M. Sc. Juan Pablo Carbajal
-----
PhD Student
University of Zürich
www.ailab.ch/carbajal
% Copyright 2011 Juan Pablo Carbajal
% [email protected]
% Thursday, April 21 2011
%
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
function [G M]= gramian(V,varargin)
%% TODO Help goes Here
%!demo
%! nV = 4;
%! v = 2*rand(2,nV)-1;
%! [G M] = gramian(v),
%! Angles = acos(G./M)*180/pi;
%! c=jet(nV);
%! for j=1:nV; plot([0 v(1,j)],[0 v(2,j)],['-;' num2str(j) ';'], ...
%! 'linewidth',2,'color',c(j,:)); hold on; end
%! axis equal
%! %---------------------------------------------------------------------------
%! The figure shows the original vectors. Compare with the values of the gramian
%! and the angles.
%!demo
%! t=linspace(0,2*pi,100)';
%! v=[sin(t), cos(t), sin(t).*cos(t).^2];
%! [G M] = gramian(v);
%! Angles = acos(G./M)*180/pi,
%! [G M] = gramian(v,'innerProduct','L2'),
%! [G M] = gramian(v,'innerProdcut','H1'),
%! [G M] = gramian(v,'innerProdcut','H2'),
%% Get dimensions
[d nV] = size(V);
if nV < 2
warning('GramianArgWarn',...
'size(V,1) == 1, Only 1 vector was provided.');
end
%% Parse options
options = {'interval','innerProduct'};
optValPos = zeros(1,numel(options));
if nargin > 1
[optGiven p]=ismember(options,varargin);
optValPos(optGiven) = p(optGiven)+1;
% Interval for interpolation and integration
interv = [0 1];
if optGiven(1)
interv = varargin{optValPos(1)};
end
%% Check nature of V
if iscell(V)
% The user is provinding vectors as function handles
if ~optGiven(1)
warning('GramianOptWarn',...
"Interval not provided, usign [0, 1].\n");
end
% TODO
elseif ismatrix(V)
% The user is providing numerical vectors
% Prepare intepolation for application of inner produc unless
% user has provided one
if optGiven(2)
optval = varargin{optValPos(2)};
if ischar(optval)
% The user selects a common space, L2,H1,H2
spaces={'L2','H1','H2'};
for i=1:numel(spaces)
if strcmpi(spaces{i},optval)
innerProd = innerSpline(i);
break;
end
end
elseif strcmpi(whos('opt').class,'function_handle')
% custom inner product
innerProd = optval;
else
error('GramianOptError',...
"Unrecongnized value for option innerProduct.\n");
end
else
% If not given use L2
innerProd = innerSpline(1);
end
% Spline interpolation
if size(interv,2) == 2
t = linspace(interv(1),interv(2),d);
else
t = interv;
end
%% Build Gramian
G = zeros(nV,nV);
M = G;
for irow = 1:nV
xrow = spline(t,V(:,irow));
for icol = 1:nV
ycol = spline(t,V(:,icol));
G(irow,icol) = innerProd(xrow,ycol,interv(1),interv(end));
if nargout > 1
M(irow,icol) = sqrt(innerProd(xrow,xrow,interv(1),interv(end))* ...
innerProd(ycol,ycol,interv(1),interv(end)));
end
end
end
end % End ismatrix
else
if ~isnumeric(V)
error('GramianArgError',...
"With no options, the argument must be a vector or a matrix.\n");
end
Dg = sumsq(V,1);
G = diag(Dg,0);
for i=1:nV
v0 = repmat(V(:,i),1,nV-1);
pos = [1:i-1 i+1:nV];
G(i,pos) = dot(v0,V(:,pos));
end
M = sqrt(kron(Dg,Dg'));
end
endfunction
function f = innerSpline(n)
switch n
case 1 %L2
f =@(x,y,a,b) quad(@(t) ppval(x,t).*ppval(y,t),a,b);
case 2 % H1, Sobolev 1st derivative
f =@(x,y,a,b) quad(@(t) ppval(x,t).*ppval(y,t) + ...
ppval(ppder(x),t).*ppval(ppder(y),t),a,b);
case 3 % H2, Soboloev 2nd derivative
f =@(x,y,a,b) quad(@(t) ppval(x,t).*ppval(y,t) + ...
ppval(ppder(x),t).*ppval(ppder(y),t) + ...
ppval(ppder(ppder(x)),t).*ppval(ppder(ppder(y)),t),a,b);
endswitch
endfunction
------------------------------------------------------------------------------
Benefiting from Server Virtualization: Beyond Initial Workload
Consolidation -- Increasing the use of server virtualization is a top
priority.Virtualization can reduce costs, simplify management, and improve
application availability and disaster protection. Learn more about boosting
the value of server virtualization. http://p.sf.net/sfu/vmware-sfdev2dev
_______________________________________________
Octave-dev mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
