function [A , c] = MinVolEllipse(P, tolerance)

% [A , c] = MinVolEllipse(P, tolerance)

% Finds the minimum volume enclsing ellipsoid (MVEE) of a set of data

% points stored in matrix P. The following optimization problem is solved: 

%

% minimize       log(det(A))

% subject to     (P_i - c)' * A * (P_i - c) <= 1

%                

% in variables A and c, where P_i is the i-th column of the matrix P. 

% The solver is based on Khachiyan Algorithm, and the final solution 

% is different from the optimal value by the pre-spesified amount of 'tolerance'.

%

% inputs:

%---------

% P : (d x N) dimnesional matrix containing N points in R^d.

% tolerance : error in the solution with respect to the optimal value.

%

% outputs:

%---------

% A : (d x d) matrix of the ellipse equation in the 'center form': 

% (x-c)' * A * (x-c) = 1 

% c : 'd' dimensional vector as the center of the ellipse. 

% 

% example:

% --------

%      P = rand(5,100);

%      [A, c] = MinVolEllipse(P, .01)

%

%      To reduce the computation time, work with the boundary points only:

%      

%      K = convhulln(P');  

%      K = unique(K(:));  

%      Q = P(:,K);

%      [A, c] = MinVolEllipse(Q, .01)

%

%

% Nima Moshtagh (nima@seas.upenn.edu)

% University of Pennsylvania

%

% December 2005

% UPDATE: Jan 2009

%

% Copyright (c) 2009, Nima Moshtagh

%All rights reserved.

%

%Redistribution and use in source and binary forms, with or without

%modification, are permitted provided that the following conditions are

%met:

%

%    * Redistributions of source code must retain the above copyright

%      notice, this list of conditions and the following disclaimer.

%    * Redistributions in binary form must reproduce the above copyright

%      notice, this list of conditions and the following disclaimer in

%      the documentation and/or other materials provided with the distribution

%

%THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

%AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

%IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

%ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE

%LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR

%CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF

%SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS

%INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN

%CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)

%ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE

%POSSIBILITY OF SUCH DAMAGE.

%

%%%%%%%%%%%%%%%%%%%%% Solving the Dual problem%%%%%%%%%%%%%%%%%%%%%%%%%%%5

% ---------------------------------

% data points 

% -----------------------------------

[d N] = size(P);

Q = zeros(d+1,N);

Q(1:d,:) = P(1:d,1:N);

Q(d+1,:) = ones(1,N);

% initializations

% -----------------------------------

count = 1;

err = 1;

u = (1/N) * ones(N,1);          % 1st iteration

% Khachiyan Algorithm

% -----------------------------------

while err > tolerance,

    X = Q * diag(u) * Q';       % X = \sum_i ( u_i * q_i * q_i')  is a (d+1)x(d+1) matrix

    M = diag(Q' * inv(X) * Q);  % M the diagonal vector of an NxN matrix

    [maximum j] = max(M);

    step_size = (maximum - d -1)/((d+1)*(maximum-1));

    new_u = (1 - step_size)*u ;

    new_u(j) = new_u(j) + step_size;

    count = count + 1;

    err = norm(new_u - u);

    u = new_u;

end

%%%%%%%%%%%%%%%%%%% Computing the Ellipse parameters%%%%%%%%%%%%%%%%%%%%%%

% Finds the ellipse equation in the 'center form': 

% (x-c)' * A * (x-c) = 1

% It computes a dxd matrix 'A' and a d dimensional vector 'c' as the center

% of the ellipse. 

U = diag(u);

% the A matrix for the ellipse

% --------------------------------------------

A = (1/d) * inv(P * U * P' - (P * u)*(P*u)' );

% center of the ellipse 

% --------------------------------------------

c = P * u;