Code covered by the BSD License

Euclidian projection on ellipsoid and conic

24 May 2010 (Updated )

Projecting a point on ellipsoid or conic in n-dimensional space

Pprod(pd)
```function [Pc Pi] = Pprod(pd)
% [Pc Pi] = Pprod(pd)
% Compute multiplication of n-polynomials, and also n leave-one-out
% polynomials
%
% INPUT:
%   pd (n x m) stores n polynomials (of order m-1) in each row
% OUTPUTS:
%   - Pc is product of all polynomials: PROD [ pd(j,:) for j=1,2, ..., n ]
%   - Pi is an array polynomials of product but leaving-one-out.
%     Pi(i,:) = PROD [ pd(j,:) for j#i ]
%
% Smart booking so as the complexity is O(n*log(n))
%
%
% Author: Bruno Luong <brunoluong@yahoo.com>
%   Original: 24-May-2010

n = size(pd,1);
k = nextpow2(n);

% Building the tree from bottom-to-top (leaf to root)
tree = cell(1,k+1);
% initialize leaf
tree{1} = struct('child', [], 'coef', num2cell(pd,2));

k = 1;
l = n; % size(tree{k},1);
while l>1
tk = tree{k};
l = ceil(l/2);
tkp1 = struct('child', cell(l,1), 'coef', cell(l,1));
% Take convolution by pair
for i=1:l
j = 2*i-1;
if j<size(tk,1)
tkp1(i) = struct('child', [j j+1], ...
'coef', conv(tk(j).coef,tk(j+1).coef));
else % alone trailing
tkp1(i) = struct('child', j, ...
'coef', tk(j).coef);
end
end % for-loop
k = k+1;
tree{k} = tkp1(:);
end % while

% Retreive the product of all the polynomials is here
Pc = tree{end}.coef;

% Building the leave-one-out products
Pi = leaveoneout(tree, cell(n, 1), size(tree,2), 1, 1);
% Put the result in array
Pi = cat(1, Pi{:});

end % Pprod

%%
% Go through the tree recursively , and take the product of polynomial
% by leaving one-out
function Pi = leaveoneout(tree, Pi, level, ind, cp)
% tree: is the tree, unchanged
% Pi: current cell list of *all* polynomials, updated when we hit
%     the first level (leaf)
% level: where we are on the tree?
% ind: indice of the node in the tree at the current level
% cp: current cumulative product of the branch currently visited

if level==1
Pi(ind) = {cp}; % done, return
else
% Recursive call
node = tree{level}(ind);
child = node.child;
up = level-1;
if length(child)==1
Pi = leaveoneout(tree, Pi, up, child, cp);
else % length(child)==1
Pi = leaveoneout(tree, Pi, up, child(1), ...
conv(cp, tree{up}(child(2)).coef));
Pi = leaveoneout(tree, Pi, up, child(2), ...
conv(cp, tree{up}(child(1)).coef));
end
end

end % leaveoneout
```