CloudRunner

... share your algorithms as a web service

Algorithm: Ellipsoidal toolbox essentials

Description:
Essential functions of the ellipsoidal toolbox for plotting in various algorithms.
See ellipsoid_toolbox_copyright.txt in the files below.
Paper: A. A. Kurzhanskiy and P. Varaiya. Ellipsoidal toolbox. Technical Report UCB/EECS-2006-46, EECS Department, University of California, Berkeley, May 2006.
Tags: toolbox
Usage: Algorithm is public.
Viewed 32255 times, called 43 times
Upload:
Empty star Empty star Empty star Empty star Empty star
0 votes
Florian Pfaff
02/10/2013 8:19 p.m. (version #1)

Run Algorithm

:
: Allow cached result?
Add varargin parameter - Remove varargin parameter
Close

Please Wait

Computation is running...

Plots

showing plot # of

Result Value(s)


		

Resource Usage

Execution of the algorithm took seconds.
Peak memory usage for running the algorithm was kilobytes (total usage: kilobytes).

Run ID / Link

ID of this run:
Link:

Result from Cache

Result has been delivered from cache and computed on (UTC).

Matlab log

(show)

		

Error

Using this algorithm in your local MATLAB environment is easy: Click here for instructions!

Usage Instructions for CloudRunner Client

  1. Download the CloudRunner Client by clicking here and place the downloaded file in your MATLAB working directory.

  2. Inside MATLAB, initialize the CloudRunner Client by calling CloudRunner:
    >> CloudRunner

    A login dialog will prompt for your CloudRunner mail address and password. For a start, you can leave the dialog empty and just click "Connect".

    Alternatively, you can provide the login credentials (or empty strings to skip login) as a parameter and hence skip the login dialog. This is useful when using CloudRunner in non-interactive scripts.
    >> CloudRunner('mail@example.com', 'password')

  3. Select this algorithm by its URL. Selecting an algorithm creates the lambda function that proxies calls to the algorithm to the server for execution:
    >> CloudRunnerSelect('http://www.cloudrunner.eu/algorithm/100/ellipsoidal-toolbox-essentials/version/1/')

    For the sake of convenience, you can also use the algorithm ID instead of its URL for public algorithms.

  4. Call functions from the algorithm like any regular local function.

Note: You can find further information on the help page.

Source Code

File:

 1 function out = ellipsoids_init(varargin)
 2 % 
 3 % ELLIPSOIDS_INIT - initializes Ellipsoidal Toolbox.
 4 % 
 5 % Any routine of Ellipsoidal Toolbox can be called with user-specified values
 6 % of different parameters. To make Ellipsoidal Toolbox as user-friendly as
 7 % possible, we provide the option to store default values of the parameters
 8 % in variable ellOptions, which is stored in MATLAB's workspace as a global
 9 % variable (i.e. it stays there unless one types 'clear all'). 
10 %
11 
12   global ellOptions;
13 
14 
15   ellOptions.version = '1.1.3';
16   ellOptions.verbose = 1; % verbosity 1 ==> ON, 0 ==> OFF
17   if ellOptions.verbose > 0
18     fprintf('Initializing Ellipsoidal Toolbox version %s ...\n', ellOptions.version);
19   end
20 
21 
22   ellOptions.abs_tol = 1e-7; % absolute tolerance
23   ellOptions.rel_tol = 1e-6; % relative tolerance
24 
25 
26   % ODE solver parameters.
27   ellOptions.time_grid          = 200;   % density of the time grid
28   ellOptions.ode_solver         = 1;     % 1 ==> RK45, 2 ==> RK23, 3 ==> Adams
29   ellOptions.norm_control       = 'on';  % on/off norm control in ODE solver
30   ellOptions.ode_solver_options = 0;     % 0 - default, 1 - user-defined
31   
32   
33   % Solver for nonlinear optimization problem with nonlinear constraints:
34   %        Minimize:   f(x)
35   %        Subject to: g(x) <= 0, h(x) = 0.
36   % Used for distance calculation.
37   % nlcp_solver = 0 - use the solver that comes with Ellipsoidal Toolbox.
38   % nlcp_solver = 1 - use the routines from MATLAB Optimization Toolbox.
39   ellOptions.nlcp_solver = 0;
40 
41 
42   ellOptions.plot2d_grid = 200; % grid density for plotting in 2D
43   ellOptions.plot3d_grid = 200; % grid density for plotting in 3D
44 
45 
46   % YALMIP settings.
47   try
48     ellOptions.sdpsettings = sdpsettings('Verbose', 0, 'warning', 0, 'cachesolvers', 1);
49   catch
50     warning('YALMIP not found, some functionality may be not accessible.');
51     ellOptions.sdpsettings = [];
52   end
53 
54   return;
 1 Ellipsoidal Toolbox Copyright
 2 
 3 
 4 Copyright (c) 2004-2008 The Regents of the University of California. All rights reserved.
 5 Permission is hereby granted, without written agreement and without license or royalty
 6 fees, to use, copy, modify, and distribute this software and its documentation for any
 7 purpose, provided that the above copyright notice and the following two paragraphs appear
 8 in all copies of this software.
 9 
10 IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
11 SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF THIS SOFTWARE AND
12 ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY
13 OF SUCH DAMAGE.
14 
15 THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
16 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
17 THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF CALIFORNIA
18 HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
 1 function h = ell_plot(x, varargin)
 2 %
 3 %
 4 % Description:
 5 % ------------
 6 %
 7 %    Wrapper for PLOT and PLOT3 functions.
 8 %    First argument must be a vector, or an array of vectors, in 1D, 2D or 3D.
 9 %    Other arguments are the same as for PLOT and PLOT3 functions.
10 %
11 %
12 % Output:
13 % -------
14 %
15 %    Plot handle.
16 %
17 %
18 % See also:
19 % ---------
20 %
21 %    PLOT, PLOT3.
22 %
23 
24 %
25 % Author:
26 % -------
27 %
28 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
29 %
30 
31   [m, n] = size(x);
32   if m > 3
33     error('ELL_PLOT: can only plot 1D, 2D and 3D vectors.');
34   end
35 
36 
37   switch m
38     case 1,
39       y = zeros(1, n);
40       h = plot(x, y, varargin{:});
41 
42     case 2,
43       xx = x(1, :);
44       yy = x(2, :);
45       h  = plot(xx, yy, varargin{:});  
46 
47     otherwise,
48       xx = x(1, :);
49       yy = x(2, :);
50       zz = x(3, :);
51       h  = plot3(xx, yy, zz, varargin{:});  
52 
53   end
54 
55   if nargout == 0
56     clear h;
57   end 
58 
59   return;
 1 function V = volume(E)
 2 %
 3 % VOLUME - returns the volume of the ellipsoid.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    V = VOLUME(E)  Computes the volume of ellipsoids in ellipsoidal array E.
10 %
11 %    The volume of ellipsoid E(q, Q) with center q and shape matrix Q is given by
12 %                  V = S sqrt(det(Q))
13 %    where S is the volume of unit ball.
14 %
15 %
16 % Output:
17 % -------
18 %
19 %    V - array of volume values, same size as E.
20 %
21 %
22 % See also:
23 % ---------
24 %
25 %    ELLIPSOID/ELLIPSOID.
26 %
27 
28 %
29 % Author:
30 % -------
31 %
32 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
33 %
34 
35   global ellOptions;
36 
37   if ~isstruct(ellOptions)
38     evalin('base', 'ellipsoids_init;');
39   end
40 
41   if ~(isa(E, 'ellipsoid'))
42     error('VOLUME: input argument must be ellipsoid.');
43   end
44 
45   [m, n] = size(E);
46   V      = [];
47   for i = 1:m
48     v = [];
49     for j = 1:n
50       Q = E(i, j).shape;
51       if isdegenerate(E(i, j))
52         S = 0;
53       else
54         N = size(Q, 1) - 1;
55         if mod(N, 2) > 0
56           k = (N + 1)/2;
57           S = (pi^k)/factorial(k);
58         else
59           k = N/2;
60           S = ((2^(2*k + 1))*(pi^k)*factorial(k))/factorial(2*k + 1);
61         end
62       end
63       v = [v S*sqrt(det(Q))];
64     end
65     V = [V; v];
66   end
67 
68   return;
 1 function I = uminus(E)
 2 %
 3 % Description:
 4 % ------------
 5 %
 6 %    Changes the sign of the center of ellipsoid.
 7 %
 8 
 9 %
10 % Author:
11 % -------
12 %
13 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
14 %
15 
16   if ~(isa(E, 'ellipsoid'))
17     error('UMINUS: input argument must be array of ellipsoids.');
18   end
19 
20   I      = E;
21   [m, n] = size(I);
22 
23   for i = 1:m
24     for j = 1:n
25       I(i, j).center = - I(i, j).center;
26     end
27   end
28 
29   return;
 1 function T = trace(E)
 2 %
 3 % TRACE - returns the trace of the ellipsoid.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    T = TRACE(E)  Computes the trace of ellipsoids in ellipsoidal array E.
10 %
11 %
12 % Output:
13 % -------
14 %
15 %    T - array of trace values, same size as E.
16 %
17 %
18 % See also:
19 % ---------
20 %
21 %    ELLIPSOID/ELLIPSOID.
22 %
23 
24 %
25 % Author:
26 % -------
27 %
28 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
29 %
30 
31   global ellOptions;
32 
33   if ~isstruct(ellOptions)
34     evalin('base', 'ellipsoids_init;');
35   end
36 
37   [m, n] = size(E);
38   T      = [];
39 
40   for i = 1:m
41     t = [];
42     for j = 1:n
43       t = [t trace(double(E(i, j)))];
44     end
45     T = [T; t];
46   end
47 
48   return;
 1 function EM = shape(E, A)
 2 %
 3 % SHAPE - modifies the shape matrix of the ellipsoid without changing its center.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    EM = SHAPE(E, A)  Modifies the shape matrices of the ellipsoids in the
10 %                      ellipsoidal array E. The centers remain untouched -
11 %                      that is the difference of the function SHAPE and
12 %                      linear transformation A*E.
13 %                      A is expected to be a scalar or a square matrix
14 %                      of suitable dimension.
15 %        
16 %
17 % Output:
18 % -------
19 %
20 %    EM - array of modified ellipsoids.
21 %
22 %
23 % See also:
24 % ---------
25 %
26 %    ELLIPSOID/ELLIPSOID.
27 %
28 
29 %
30 % Author:
31 % -------
32 %
33 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
34 %
35 
36   if ~(isa(A, 'double')) | ~(isa(E, 'ellipsoid'))
37     msg = sprintf('SHAPE: expected arguments are:\n');
38     msg = sprintf('%s       - array of ellipsoids of the same dimension,\n', msg);
39     msg = sprintf('%s       - scalar, or square matrix of the same dimension as ellipsoids.\n', msg);
40     error(msg);
41   end
42 
43   [m, n] = size(A); 
44   if m ~= n
45     error('SHAPE: only square matrices are allowed.');
46   end
47   d      = dimension(E);
48   k      = max(max(d));
49   l      = min(min(d));
50   if ((k ~= l) & (n ~= 1) & (m ~= 1)) | ((k ~= n) & (n ~= 1) & (m ~= 1))
51     error('SHAPE: dimensions do not match.');
52   end
53 
54   EM     = [];
55   [m, n] = size(E);
56   for i = 1:m
57     for j = 1:n
58      Q    = A*(E(i, j).shape)*A';
59      Q    = 0.5*(Q + Q');
60      r(j) = ellipsoid(E(i, j).center, Q);
61     end
62     EM = [EM; r];
63     clear r;
64   end
65 
66   return;
  1 function [res, x] = rho(E, L)
  2 %
  3 % RHO - computes the values of the support function for given ellipsoid
  4 %       and given direction.
  5 %
  6 %
  7 % Description:
  8 % ------------
  9 %
 10 %         RES = RHO(E, L)  Computes the support function of the ellipsoid E
 11 %                          in directions specified by the columns of matrix L.
 12 %                          Or, if E is array of ellipsoids, L is expected to be
 13 %                          a single vector.
 14 %
 15 %    [RES, X] = RHO(E, L)  Computes the support function of the ellipsoid E
 16 %                          in directions specified by the columns of matrix L,
 17 %                          and boundary points X of this ellipsoid that correspond
 18 %                          to directions in L.
 19 %                          Or, if E is array of ellipsoids, and L - single vector,
 20 %                          then support functions and corresponding boundary
 21 %                          points are computed for all the given ellipsoids in
 22 %                          the array in the specified direction L.
 23 %
 24 %    The support function is defined as
 25 %       (1)  rho(l | E) = sup { <l, x> : x belongs to E }.
 26 %    For ellipsoid E(q,Q), where q is its center and Q - shape matrix,
 27 %    it is simplified to
 28 %       (2)  rho(l | E) = <q, l> + sqrt(<l, Ql>)
 29 %    Vector x, at which the maximum at (1) is achieved is defined by
 30 %       (3)  q + Ql/sqrt(<l, Ql>)
 31 %
 32 %
 33 % Output:
 34 % -------
 35 %
 36 %    RES - the values of the support function for the specified ellipsoid E
 37 %          and directions L. Or, if E is an array of ellipsoids, and L - single
 38 %          vector, then these are values of the support function for all the
 39 %          ellipsoids in the array in the given direction.
 40 %      X - boundary points of the ellipsoid E that correspond to directions in L.
 41 %          Or, if E is an array of ellipsoids, and L - single vector,
 42 %          then these are boundary points of all the ellipsoids in the array
 43 %          in the given direction.
 44 %
 45 %
 46 % See also:
 47 % ---------
 48 %
 49 %    ELLIPSOID/ELLIPSOID.
 50 %
 51 
 52 %
 53 % Author:
 54 % -------
 55 %
 56 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 57 %
 58 
 59   if ~(isa(E, 'ellipsoid'))
 60     error('RHO: first argument must be ellipsoid.');
 61   end
 62 
 63   [m, n] = size(E);
 64   [k, d] = size(L);
 65   if (m > 1) | (n > 1)
 66     if d > 1
 67       msg = sprintf('RHO: arguments must be single ellipsoid and matrix of direction vectors,\n     or array of ellipsoids and single direction vector.');
 68       error(msg);
 69     end
 70     ea = 1;
 71   else
 72     ea = 0;
 73   end
 74 
 75   dims = dimension(E);
 76   mn   = min(min(dims));
 77   mx   = max(max(dims));
 78   if mn ~= mx
 79     error('RHO: ellipsoids in the array must be of the same dimension.');
 80   end
 81   if mn ~= k
 82     error('RHO: dimensions of the direction vector and the ellipsoid do not match.');
 83   end
 84 
 85   
 86   res = [];
 87   x   = [];
 88   if ea > 0 % multiple ellipsoids, one direction
 89     for i = 1:m
 90       r = [];
 91       for j = 1:n
 92         q  = E(i, j).center;
 93         Q  = E(i, j).shape;
 94         sr = sqrt(L'*Q*L);
 95         if sr == 0
 96           sr = eps;
 97         end
 98         r  = [r (q'*L + sr)];
 99         x  = [x (((Q*L)/sr) + q)];
100       end
101       res = [res; r];
102     end
103   else % one ellipsoid, multiple directions
104     q = E.center;
105     Q = E.shape;
106     for i = 1:d
107       l   = L(:, i);
108       sr  = sqrt(l'*Q*l);
109       if sr == 0
110         sr = eps;
111       end
112       res = [res (q'*l + sr)];
113       x   = [x (((Q*l)/sr) + q)];
114     end
115   end
116 
117   return;
 1 function EP = projection(E, B)
 2 %
 3 % PROJECTION - computes projection of the ellipsoid onto the given subspace.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    EP = PROJECTION(E, B)  Computes projection of the ellipsoid E onto a subspace,
10 %                           specified by orthogonal basis vectors B.
11 %                           E can be an array of ellipsoids of the same dimension.
12 %                           Columns of B must be orthogonal vectors.
13 %
14 %
15 % Output:
16 % -------
17 %
18 %    EP - projected ellipsoid (or array of ellipsoids), generally, of lower dimension.
19 %
20 %
21 % See also:
22 % ---------
23 %
24 %    ELLIPSOID/ELLIPSOID.
25 %
26 
27 %
28 % Author:
29 % -------
30 %
31 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
32 %
33 
34   global ellOptions;
35 
36   if ~isstruct(ellOptions)
37     evalin('base', 'ellipsoids_init;');
38   end
39 
40   if ~(isa(E, 'ellipsoid')) | ~(isa(B, 'double'))
41     error('PROJECTION: arguments must be array of ellipsoids and matrix with orthogonal columns.');
42   end
43 
44   [k, l] = size(B);
45   dims   = dimension(E);
46   m      = min(dims); m = min(m);
47   n      = max(dims); n = max(n);
48   if (m ~= n)
49     error('PROJECTION: ellipsoids in the array must be of the same dimenion.');
50   end
51   if (k ~= n)
52     error('PROJECTION: dimension of basis vectors does not dimension of ellipsoids.');
53   end
54   if (k < l)
55     msg = sprintf('PROJECTION: number of basis vectors must be less or equal to %d.', n);
56     error(msg);
57   end
58 
59   % check the orthogonality of the columns of B
60   for i = 1:(l - 1)
61     v = B(:, i);
62     for j = (i + 1):l
63       if abs(v'*B(:, j)) > ellOptions.abs_tol
64         error('PROJECTION: basis vectors must be orthogonal.');
65       end
66     end
67   end
68 
69   % normalize the basis vectors
70   for i = 1:l
71     BB(:, i) = B(:, i)/norm(B(:, i));
72   end
73 
74   % compute projection
75   [m, n] = size(E);
76   for i = 1:m
77     for j = 1:n
78       r(j) = BB'*E(i, j);
79     end
80     if i == 1
81       EP = r;
82     else
83       EP = [EP; r];
84     end
85     clear r;
86   end
87 
88   return;
 1 function LC = rm_bad_directions(Q1, Q2, L)
 2 %
 3 % RM_BAD_DIRECTIONS - remove bad directions from the given list.
 4 %                     Bad directions are those which should not be used
 5 %                     for the support function of geometric difference
 6 %                     of two ellipsoids.
 7 %
 8 
 9   LC = [];
10   T  = ell_simdiag(Q2, Q1);
11   a  = min(abs(diag(T*Q1*T')));
12   if a < 1
13     return;
14   end
15 
16   n = size(L, 2);
17   for i = 1:n
18     l = L(:, i);
19     if (sqrt(l'*Q1*l)/sqrt(l'*Q2*l)) <= a
20       LC = [LC l];
21     end
22   end
23 
24   return;
 1 function R = regularize(Q)
 2 %
 3 % REGULARIZE - regularization of singular symmetric matrix.
 4 %
 5 
 6   global ellOptions;
 7 
 8   if Q ~= Q'
 9     error('REGULARIZE: matrix must be symmetric.');
10   end
11 
12   [m, n] = size(Q);
13   r      = rank(Q);
14 
15   if r < n
16     [U S V] = svd(Q);
17     E       = ellOptions.abs_tol * eye(n - r);
18     R       = Q + (U * [zeros(r, r) zeros(r, (n-r)); zeros((n-r), r) E] * U');
19     R       = 0.5*(R + R');
20   else
21     R = Q;
22   end
23 
24   return;
 1 function res = my_color_table(ch)
 2 %
 3 % MY_COLOR_TABLE - returns the code of the color defined by single letter.
 4 %
 5 
 6   if ~(ischar(ch))
 7     res = [0 0 0];
 8     return;
 9   end
10 
11   switch ch
12     case 'r',
13       res = [1 0 0];
14 
15     case 'g',
16       res = [0 1 0];
17 
18     case 'b',
19       res = [0 0 1];
20 
21     case 'y',
22       res = [1 1 0];
23 
24     case 'c',
25       res = [0 1 1];
26 
27     case 'm',
28       res = [1 0 1];
29       
30     case 'w',
31       res = [1 1 1];
32 
33     otherwise,
34       res = [0 0 0];
35   end
36 
37   return;
 1 function x = ellbndr_3d(E)
 2 %
 3 % ELLBNDR_3D - compute the boundary of 3D ellipsoid.
 4 %
 5 
 6   global ellOptions;
 7 
 8   M   = ellOptions.plot3d_grid/2;
 9   N   = M/2;
10 
11   psy = linspace(0, pi, N);
12   phi = linspace(0, 2*pi, M);
13 
14   l   = [];
15   for i = 2:(N - 1)
16     arr = cos(psy(i))*ones(1, M);
17     l   = [l [cos(phi)*sin(psy(i)); sin(phi)*sin(psy(i)); arr]];
18   end
19 
20   [r, x] = rho(E, l);
21 
22   return;
 1 function x = ellbndr_2d(E)
 2 %
 3 % ELLBNDR_2D - compute the boundary of 2D ellipsoid.
 4 %
 5 
 6   global ellOptions;
 7 
 8   N      = ellOptions.plot2d_grid;
 9   phi    = linspace(0, 2*pi, N);
10   l      = [cos(phi); sin(phi)];
11   [r, x] = rho(E, l);
12   x      = [x x(:, 1)];
13 
14   return;
 1 function P = polar(E)
 2 %
 3 % POLAR - computes the polar ellipsoids.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    P = POLAR(E)  Computes the polar ellipsoids for those ellipsoids in E,
10 %                  for which the origin is an interior point.
11 %                  For those ellipsoids in E, for which this condition
12 %                  does not hold, an empty ellipsoid is returned.
13 %
14 %
15 %    Given ellipsoid E(q, Q) where q is its center, and Q - its shape matrix,
16 %    the polar set to E(q, Q) is defined as follows:
17 %
18 %             P = { l in R^n  | <l, q> + sqrt(<l, Q l>) <= 1 }
19 %
20 %    If the origin is an interior point of ellipsoid E(q, Q),
21 %    then its polar set P is an ellipsoid.
22 %
23 %
24 % Output:
25 % -------
26 %
27 %    P - array of polar ellipsoids.
28 %
29 %
30 % See also:
31 % ---------
32 %
33 %    ELLIPSOID/ELLIPSOID.
34 %
35 
36 %
37 % Author:
38 % -------
39 %
40 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
41 %
42 
43   if ~(isa(E, 'ellipsoid'))
44     error('POLAR: input argument must be array of ellipsoids.');
45   end
46 
47   [m, n] = size(E);
48   P      = [];
49 
50   for i = 1:m
51     PP = [];
52     for j = 1:n
53       if isdegenerate(E(i, j))
54         PP = [PP ellipsoid];
55       else
56         [q, Q] = parameters(E(i, j));
57         d      = size(Q, 2);
58         z      = zeros(d, 1);
59         chk    = (z' - q') * ell_inv(Q) * (z - q);
60         if chk < 1
61           M  = ell_inv(Q - q*q');
62           M  = 0.5*(M + M');
63           w  = -M * q;
64           W  = (1 + q'*M*q)*M;
65           PP = [PP ellipsoid(w, W)];
66         else
67           PP = [PP ellipsoid];
68 	end
69       end
70     end
71     P = [P; PP];
72   end
73 
74   return;
 1 function res = plus(X, Y)
 2 %
 3 %
 4 % Description:
 5 % ------------
 6 %
 7 %    Operation
 8 %              E + b
 9 %    where E is an ellipsoid in R^n, and b - vector in R^n.
10 %    If E(q, Q) is an ellipsoid with center q and shape matrix Q, then
11 %          E(q, Q) + b = E(q + b, Q).
12 %
13 %
14 % Output:
15 % -------
16 %
17 %    Resulting ellipsoid E(q + b, Q).
18 %
19 %
20 % See also:
21 % ---------
22 %
23 %    ELLIPSOID/ELLIPSOID.
24 %
25 
26 %
27 % Author:
28 % -------
29 %
30 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
31 %
32 
33   if isa(X, 'ellipsoid') & ~(isa(Y, 'double'))
34     error('PLUS: this operation is only permitted between ellipsoid and vector in R^n.');
35   end
36   if isa(Y, 'ellipsoid') & ~(isa(X, 'double'))
37     error('PLUS: this operation is only permitted between ellipsoid and vector in R^n.');
38   end
39 
40   if isa(X, 'ellipsoid')
41     E = X;
42     b = Y;
43   else
44     E = Y;
45     b = X;
46   end
47 
48   d = dimension(E);
49   m = max(d);
50   n = min(d);
51   if m ~= n
52     error('PLUS: all ellipsoids in the array must be of the same dimension.');
53   end
54 
55   [k, l] = size(b);
56   if (l ~= 1) | (k ~= n)
57     error('PLUS: vector dimension does not match.');
58   end
59 
60   [m, n] = size(E);
61   for i = 1:m
62     for j = 1:n
63       r(j)        = E(i, j);
64       r(j).center = E(i, j).center + b;
65     end
66     if i == 1
67       res = r;
68     else
69       res = [res; r];
70     end
71     clear r;
72   end
73 
74   return;
  1 function plot3(varargin)
  2 %
  3 % PLOT3 - plots ellipsoids in 2D or 3D.
  4 %
  5 %
  6 % Description:
  7 % ------------
  8 %
  9 % PLOT3(E, OPTIONS) plots ellipsoid E if 1 <= dimension(E) <= 3.
 10 %
 11 %                 PLOT3(E)  Plots E in default (red) color.
 12 %             PLOT3(EA, E)  Plots array of ellipsoids EA and single ellipsoid E.
 13 %  PLOT3(E1, 'g', E2, 'b')  Plots E1 in green and E2 in blue color.
 14 %       PLOT3(EA, Options)  Plots EA using options given in the Options structure.
 15 %
 16 %
 17 % Options.newfigure    - if 1, each plot command will open a new figure window.
 18 % Options.fill         - if 1, ellipsoids in 2D will be filled with color.
 19 % Options.width        - line width for 1D and 2D plots.
 20 % Options.color        - sets default colors in the form [x y z].
 21 % Options.shade = 0-1  - level of transparency (0 - transparent, 1 - opaque).
 22 %
 23 %
 24 % Output:
 25 % -------
 26 %
 27 %    None.
 28 %
 29 %
 30 % See also:
 31 % ---------
 32 %
 33 %    ELLIPSOID/ELLIPSOID, ELLIPSOID/PLOT.
 34 %
 35 
 36 % 
 37 % Author:
 38 % -------
 39 %
 40 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 41 %
 42 
 43   global ellOptions;
 44 
 45   if ~isstruct(ellOptions)
 46     evalin('base', 'ellipsoids_init;');
 47   end
 48 
 49   nai = nargin;
 50   E   = varargin{1};
 51   if ~isa(E, 'ellipsoid')
 52     error('PLOT: input argument must be an ellipsoid.');
 53   end
 54 
 55   if nai > 1
 56     if isstruct(varargin{nai})
 57       Options = varargin{nai};
 58       nai     = nai - 1;
 59     else
 60       Options = [];
 61     end
 62   else
 63     Options = [];
 64   end
 65 
 66   if ~isfield(Options, 'newfigure')
 67     Options.newfigure = 0;
 68   end
 69 
 70   ucolor    = [];
 71   vcolor    = [];
 72   ells      = [];
 73   ell_count = 0;
 74 
 75   for i = 1:nai
 76     if isa(varargin{i}, 'ellipsoid')
 77       E      = varargin{i};
 78       [m, n] = size(E);
 79       cnt    = m * n;
 80       E1     = reshape(E, 1, cnt);
 81       ells   = [ells E1];
 82       if (i < nai) & ischar(varargin{i + 1})
 83         clr = my_color_table(varargin{i + 1});
 84         val = 1;
 85       else
 86         clr = [0 0 0];
 87         val = 0;
 88       end
 89       for j = (ell_count + 1):(ell_count + cnt)
 90         ucolor(j) = val;
 91         vcolor    = [vcolor; clr];
 92       end
 93       ell_count = ell_count + cnt;
 94     end
 95   end
 96 
 97   if ~isfield(Options, 'color')
 98     % Color maps:
 99     %    hsv       - Hue-saturation-value color map.
100     %    hot       - Black-red-yellow-white color map.
101     %    gray      - Linear gray-scale color map.
102     %    bone      - Gray-scale with tinge of blue color map.
103     %    copper    - Linear copper-tone color map.
104     %    pink      - Pastel shades of pink color map.
105     %    white     - All white color map.
106     %    flag      - Alternating red, white, blue, and black color map.
107     %    lines     - Color map with the line colors.
108     %    colorcube - Enhanced color-cube color map.
109     %    vga       - Windows colormap for 16 colors.
110     %    jet       - Variant of HSV.
111     %    prism     - Prism color map.
112     %    cool      - Shades of cyan and magenta color map.
113     %    autumn    - Shades of red and yellow color map.
114     %    spring    - Shades of magenta and yellow color map.
115     %    winter    - Shades of blue and green color map.
116     %    summer    - Shades of green and yellow color map.
117     
118     auxcolors  = hsv(ell_count);
119     colors     = auxcolors;
120     multiplier = 7;
121     if mod(size(auxcolors, 1), multiplier) == 0
122       multiplier = multiplier + 1;
123     end
124     
125     for i = 1:ell_count
126       jj           = mod(i*multiplier, size(auxcolors, 1)) + 1;
127       colors(i, :) = auxcolors(jj, :);
128     end
129     colors        = flipud(colors);
130     Options.color = colors;
131   else
132     if size(Options.color, 1) ~= ell_count
133       if size(Options.color, 1) > ell_count
134         Options.color = Options.color(1:ell_count, :);
135       else
136         Options.color = repmat(Options.color, ell_count, 1);
137       end
138     end
139   end
140 
141   if ~isfield(Options, 'shade')
142     Options.shade = 0.4*ones(1, ell_count);
143   else
144     [m, n] = size(Options.shade);
145     m      = m * n;
146     if m == 1
147       Options.shade = Options.shade * ones(1, ell_count);
148     else
149       Options.shade = reshape(Options.shade, 1, m);
150       if m < ell_count
151         for i = (m + 1):ell_count
152           Options.shade = [Options.shade 0.4];
153         end
154       end
155     end
156   end
157 
158   if ~isfield(Options, 'width')
159     Options.width = ones(1, ell_count);
160   else
161     [m, n] = size(Options.width);
162     m      = m * n;
163     if m == 1
164       Options.width = Options.width * ones(1, ell_count);
165     else
166       Options.width = reshape(Options.width, 1, m);
167       if m < ell_count
168         for i = (m + 1):ell_count
169           Options.width = [Options.width 1];
170         end
171       end
172     end
173   end
174 
175   if ~isfield(Options, 'fill')
176     Options.fill = zeros(1, ell_count);
177   else
178     [m, n] = size(Options.fill);
179     m      = m * n;
180     if m == 1
181       Options.fill = Options.fill * ones(1, ell_count);
182     else
183       Options.fill = reshape(Options.fill, 1, m);
184       if m < ell_count
185         for i = (m + 1):ell_count
186           Options.fill = [Options.fill 0];
187         end
188       end
189     end
190   end
191 
192   if size(Options.color, 1) < ell_count
193     error('PLOT: not enough colors.');
194   end
195 
196   dims = dimension(ells);
197   m    = min(dims);
198   n    = max(dims);
199   if m ~= n
200     error('PLOT: ellipsoids must be of the same dimension.');
201   end
202   if (n > 3) | (n < 1)
203     error('PLOT: ellipsoid dimension can be 1, 2 or 3.');
204   end
205 
206   if ellOptions.verbose > 0
207     if ell_count == 1
208       fprintf('Plotting ellipsoid...\n');
209     else
210       fprintf('Plotting %d ellipsoids...\n', ell_count);
211     end
212   end
213 
214   ih = ishold;
215 
216   for i = 1:ell_count
217     if Options.newfigure ~= 0
218       figure;
219     else
220       newplot;
221     end
222 
223     hold on;
224 
225     E = ells(i);
226     q = E.center;
227     Q = E.shape;
228 
229     if ucolor(i) == 1
230       clr = vcolor(i, :);
231     else
232       clr = Options.color(i, :);
233     end
234       
235     switch n
236       case 2,
237         x = ellbndr_2d(E);
238         if Options.fill(i) ~= 0
239           fill(x(1, :), x(2, :), clr);
240         end
241         h = ell_plot(x);
242         set(h, 'Color', clr, 'LineWidth', Options.width(i));
243         h = ell_plot(q, '.');
244         set(h, 'Color', clr);
245 
246       case 3,
247         x    = ellbndr_3d(E);
248         chll = convhulln(x');
249         vs   = size(x, 2);
250         patch('Vertices', x', 'Faces', chll, ...
251               'FaceVertexCData', clr(ones(1, vs), :), 'FaceColor', 'flat', ...
252               'EdgeColor', clr, 'FaceAlpha', 0);
253         %shading interp;
254         lighting phong;
255         material('metal');
256         view(3);
257         %camlight('headlight','local');
258         %camlight('headlight','local');
259         %camlight('right','local');
260         %camlight('left','local');
261 
262       otherwise,
263         h = ell_plot([(q-sqrt(Q)) (q+sqrt(Q))]);
264         set(h, 'Color', clr, 'LineWidth', Options.width(i));
265         h = ell_plot(q(1, 1), '*');
266         set(h, 'Color', clr);
267 
268     end
269   end
270 
271   if ih == 0;
272     hold off;
273   end
274 
275   return;
The file is not a text file or too large to be displayed here. You can still download the file using the link below.
  1 function plot(varargin)
  2 %
  3 % PLOT - plots ellipsoids in 2D or 3D.
  4 %
  5 %
  6 % Description:
  7 % ------------
  8 %
  9 % PLOT(E, OPTIONS) plots ellipsoid E if 1 <= dimension(E) <= 3.
 10 %
 11 %                  PLOT(E)  Plots E in default (red) color.
 12 %              PLOT(EA, E)  Plots array of ellipsoids EA and single ellipsoid E.
 13 %   PLOT(E1, 'g', E2, 'b')  Plots E1 in green and E2 in blue color.
 14 %        PLOT(EA, Options)  Plots EA using options given in the Options structure.
 15 %
 16 %
 17 % Options.newfigure    - if 1, each plot command will open a new figure window.
 18 % Options.fill         - if 1, ellipsoids in 2D will be filled with color.
 19 % Options.width        - line width for 1D and 2D plots.
 20 % Options.color        - sets default colors in the form [x y z].
 21 % Options.shade = 0-1  - level of transparency (0 - transparent, 1 - opaque).
 22 %
 23 %
 24 % Output:
 25 % -------
 26 %
 27 %    None.
 28 %
 29 %
 30 % See also:
 31 % ---------
 32 %
 33 %    ELLIPSOID/ELLIPSOID.
 34 %
 35 
 36 % 
 37 % Author:
 38 % -------
 39 %
 40 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 41 %
 42 
 43   global ellOptions;
 44 
 45   if ~isstruct(ellOptions)
 46     evalin('base', 'ellipsoids_init;');
 47   end
 48 
 49   nai = nargin;
 50   E   = varargin{1};
 51   if ~isa(E, 'ellipsoid')
 52     error('PLOT: input argument must be an ellipsoid.');
 53   end
 54 
 55   if nai > 1
 56     if isstruct(varargin{nai})
 57       Options = varargin{nai};
 58       nai     = nai - 1;
 59     else
 60       Options = [];
 61     end
 62   else
 63     Options = [];
 64   end
 65 
 66   if ~isfield(Options, 'newfigure')
 67     Options.newfigure = 0;
 68   end
 69 
 70   ucolor    = [];
 71   vcolor    = [];
 72   ells      = [];
 73   ell_count = 0;
 74 
 75   for i = 1:nai
 76     if isa(varargin{i}, 'ellipsoid')
 77       E      = varargin{i};
 78       [m, n] = size(E);
 79       cnt    = m * n;
 80       E1     = reshape(E, 1, cnt);
 81       ells   = [ells E1];
 82       if (i < nai) & ischar(varargin{i + 1})
 83         clr = my_color_table(varargin{i + 1});
 84         val = 1;
 85       else
 86         clr = [0 0 0];
 87         val = 0;
 88       end
 89       for j = (ell_count + 1):(ell_count + cnt)
 90         ucolor(j) = val;
 91         vcolor    = [vcolor; clr];
 92       end
 93       ell_count = ell_count + cnt;
 94     end
 95   end
 96 
 97   if ~isfield(Options, 'color')
 98     % Color maps:
 99     %    hsv       - Hue-saturation-value color map.
100     %    hot       - Black-red-yellow-white color map.
101     %    gray      - Linear gray-scale color map.
102     %    bone      - Gray-scale with tinge of blue color map.
103     %    copper    - Linear copper-tone color map.
104     %    pink      - Pastel shades of pink color map.
105     %    white     - All white color map.
106     %    flag      - Alternating red, white, blue, and black color map.
107     %    lines     - Color map with the line colors.
108     %    colorcube - Enhanced color-cube color map.
109     %    vga       - Windows colormap for 16 colors.
110     %    jet       - Variant of HSV.
111     %    prism     - Prism color map.
112     %    cool      - Shades of cyan and magenta color map.
113     %    autumn    - Shades of red and yellow color map.
114     %    spring    - Shades of magenta and yellow color map.
115     %    winter    - Shades of blue and green color map.
116     %    summer    - Shades of green and yellow color map.
117     
118     auxcolors  = hsv(ell_count);
119     colors     = auxcolors;
120     multiplier = 7;
121     if mod(size(auxcolors, 1), multiplier) == 0
122       multiplier = multiplier + 1;
123     end
124     
125     for i = 1:ell_count
126       jj           = mod(i*multiplier, size(auxcolors, 1)) + 1;
127       colors(i, :) = auxcolors(jj, :);
128     end
129     colors        = flipud(colors);
130     Options.color = colors;
131   else
132     if size(Options.color, 1) ~= ell_count
133       if size(Options.color, 1) > ell_count
134         Options.color = Options.color(1:ell_count, :);
135       else
136         Options.color = repmat(Options.color, ell_count, 1);
137       end
138     end
139   end
140 
141   if ~isfield(Options, 'shade')
142     Options.shade = 0.4*ones(1, ell_count);
143   else
144     [m, n] = size(Options.shade);
145     m      = m * n;
146     if m == 1
147       Options.shade = Options.shade * ones(1, ell_count);
148     else
149       Options.shade = reshape(Options.shade, 1, m);
150       if m < ell_count
151         for i = (m + 1):ell_count
152           Options.shade = [Options.shade 0.4];
153         end
154       end
155     end
156   end
157 
158   if ~isfield(Options, 'width')
159     Options.width = ones(1, ell_count);
160   else
161     [m, n] = size(Options.width);
162     m      = m * n;
163     if m == 1
164       Options.width = Options.width * ones(1, ell_count);
165     else
166       Options.width = reshape(Options.width, 1, m);
167       if m < ell_count
168         for i = (m + 1):ell_count
169           Options.width = [Options.width 1];
170         end
171       end
172     end
173   end
174 
175   if ~isfield(Options, 'fill')
176     Options.fill = zeros(1, ell_count);
177   else
178     [m, n] = size(Options.fill);
179     m      = m * n;
180     if m == 1
181       Options.fill = Options.fill * ones(1, ell_count);
182     else
183       Options.fill = reshape(Options.fill, 1, m);
184       if m < ell_count
185         for i = (m + 1):ell_count
186           Options.fill = [Options.fill 0];
187         end
188       end
189     end
190   end
191 
192   if size(Options.color, 1) < ell_count
193     error('PLOT: not enough colors.');
194   end
195 
196   dims = dimension(ells);
197   m    = min(dims);
198   n    = max(dims);
199   if m ~= n
200     error('PLOT: ellipsoids must be of the same dimension.');
201   end
202   if (n > 3) | (n < 1)
203     error('PLOT: ellipsoid dimension can be 1, 2 or 3.');
204   end
205 
206   if ellOptions.verbose > 0
207     if ell_count == 1
208       fprintf('Plotting ellipsoid...\n');
209     else
210       fprintf('Plotting %d ellipsoids...\n', ell_count);
211     end
212   end
213 
214   ih = ishold;
215 
216   for i = 1:ell_count
217     if Options.newfigure ~= 0
218       figure;
219     else
220       newplot;
221     end
222 
223     hold on;
224 
225     E = ells(i);
226     q = E.center;
227     Q = E.shape;
228 
229     if ucolor(i) == 1
230       clr = vcolor(i, :);
231     else
232       clr = Options.color(i, :);
233     end
234       
235     switch n
236       case 2,
237         x = ellbndr_2d(E);
238         if Options.fill(i) ~= 0
239           fill(x(1, :), x(2, :), clr);
240         end
241         h = ell_plot(x);
242         set(h, 'Color', clr, 'LineWidth', Options.width(i));
243         h = ell_plot(q, '.');
244         set(h, 'Color', clr);
245         hasbehavior(h, 'legend', false);
246         
247       case 3,
248         x    = ellbndr_3d(E);
249         chll = convhulln(x');
250         vs   = size(x, 2);
251         patch('Vertices', x', 'Faces', chll, ...
252               'FaceVertexCData', clr(ones(1, vs), :), 'FaceColor', 'flat', ...
253               'FaceAlpha', Options.shade(1, i));
254         shading interp;
255         lighting phong;
256         material('metal');
257         view(3);
258         %camlight('headlight','local');
259         %camlight('headlight','local');
260         %camlight('right','local');
261         %camlight('left','local');
262 
263       otherwise,
264         h = ell_plot([(q-sqrt(Q)) (q+sqrt(Q))]);
265         set(h, 'Color', clr, 'LineWidth', Options.width(i));
266         h = ell_plot(q(1, 1), '*');
267         set(h, 'Color', clr);
268 
269     end
270   end
271 
272   if ih == 0;
273     hold off;
274   end
275 
276   return;
 1 function [q, Q] = parameters(E)
 2 %
 3 % PARAMETERS - returns parameters of the ellipsoid.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    [q, Q] = PARAMETERS(E)  Extracts the values of the center q and
10 %                            the shape matrix Q from the ellipsoid object E.
11 %
12 %
13 % Output:
14 % -------
15 %
16 %    q - center of the ellipsoid E.
17 %    Q - shape matrix of the ellipsoid E.
18 %
19 %
20 % See also:
21 % ---------
22 %
23 %    ELLIPSOID/ELLIPSOID, DIMENSION, ISDEGENERATE.
24 %
25 
26 %
27 % Author:
28 % -------
29 %
30 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
31 %
32 
33   [m, n] = size(E);
34   if (m > 1) | (n > 1)
35     error('PARAMETERS: the argument of this function must be single ellipsoid.');
36   end
37   
38   if nargout < 2
39     q = E.shape;
40   else
41     q = E.center;
42     Q = E.shape;
43   end
44 
45   return;
 1 function res = ne(E1, E2)
 2 %
 3 %
 4 % Description:
 5 % ------------
 6 %
 7 %    The opposite of EQ.
 8 %
 9 
10 %
11 % Author:
12 % -------
13 %
14 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
15 %
16 
17   res = ~(eq(E1, E2));
18 
19   return;
 1 function res = mtimes(A, E)
 2 %
 3 %
 4 % Description:
 5 % ------------
 6 %
 7 %    Multiplication of the ellipsoid by a matrix or a scalar.
 8 %    If E(q,Q) is an ellipsoid, and A - matrix of suitable dimensions,
 9 %    then
10 %          A E(q, Q) = E(Aq, AQA').
11 %        
12 %
13 %
14 % Output:
15 % -------
16 %
17 %    Resulting ellipsoid E(Aq, AQA').
18 %
19 %
20 % See also:
21 % ---------
22 %
23 %    ELLIPSOID/ELLIPSOID.
24 %
25 
26 %
27 % Author:
28 % -------
29 %
30 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
31 %
32 
33   if ~(isa(A, 'double')) | ~(isa(E, 'ellipsoid'))
34     msg = sprintf('MTIMES: first multiplier is expected to be a matrix or a scalar,\n        and second multiplier - an ellipsoid.');
35     error(msg);
36   end
37 
38   [m, n] = size(A); 
39   d      = dimension(E);
40   k      = max(d);
41   l      = min(d);
42   if ((k ~= l) & (n ~= 1) & (m ~= 1)) | ((k ~= n) & (n ~= 1) & (m ~= 1))
43     error('MTIMES: dimensions do not match.');
44   end
45 
46   [m, n] = size(E);
47   for i = 1:m
48     for j = 1:n
49      Q    = A*(E(i, j).shape)*A';
50      Q    = 0.5*(Q + Q');
51      r(j) = ellipsoid(A*(E(i, j).center), Q);
52     end
53     if i == 1
54       res = r;
55     else
56       res = [res; r];
57     end
58     clear r;
59   end
60 
61   return;
 1 function EO = move2origin(E)
 2 %
 3 % MOVE2ORIGIN - moves ellipsoids in the given array to the origin.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    EO = MOVE2ORIGIN(E)  Replaces the centers of ellipsoids in E with zero vectors.
10 %
11 %
12 % Output:
13 % -------
14 %
15 %    EO - array of ellipsoids with the same shapes as in E centered at the origin.
16 %
17 %
18 % See also:
19 % ---------
20 %
21 %    ELLIPSOID/ELLIPSOID.
22 %
23 
24 %
25 % Author:
26 % -------
27 %
28 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
29 %
30 
31   if ~(isa(E, 'ellipsoid'))
32     error('MOVE2ORIGIN: argument must be array of ellipsoids.');
33   end
34 
35   EO     = E;
36   [m, n] = size(EO);
37 
38   for i = 1:m
39     for j = 1:n
40       d               = dimension(EO(i, j));
41       EO(i, j).center = zeros(d, 1);
42     end
43   end
44 
45   return;
 1 function res = minus(E, b)
 2 %
 3 %
 4 % Description:
 5 % ------------
 6 %
 7 %    Operation
 8 %              E - b
 9 %    where E is an ellipsoid in R^n, and b - vector in R^n.
10 %    If E(q, Q) is an ellipsoid with center q and shape matrix Q, then
11 %          E(q, Q) - b = E(q - b, Q).
12 %
13 %
14 % Output:
15 % -------
16 %
17 %    Resulting ellipsoid E(q - b, Q).
18 %
19 %
20 % See also:
21 % ---------
22 %
23 %    ELLIPSOID/ELLIPSOID.
24 %
25 
26 %
27 % Author:
28 % -------
29 %
30 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
31 %
32 
33   if ~(isa(E, 'ellipsoid'))
34     error('MINUS: first argument must be ellipsoid.');
35   end
36   if isa(E, 'ellipsoid') & ~(isa(b, 'double'))
37     error('MINUS: this operation is only permitted between ellipsoid and vector in R^n.');
38   end
39 
40   d = dimension(E);
41   m = max(d);
42   n = min(d);
43   if m ~= n
44     error('MINUS: all ellipsoids in the array must be of the same dimension.');
45   end
46 
47   [k, l] = size(b);
48   if (l ~= 1) | (k ~= n)
49     error('MINUS: vector dimension does not match.');
50   end
51 
52   [m, n] = size(E);
53   for i = 1:m
54     for j = 1:n
55       r(j)        = E(i, j);
56       r(j).center = E(i, j).center - b;
57     end
58     if i == 1
59       res = r;
60     else
61       res = [res; r];
62     end
63     clear r;
64   end
65 
66   return;
  1 function IA = minksum_ea(E, L)
  2 %
  3 % MINKSUM_IA - computation of internal approximating ellipsoids of the geometric
  4 %              sum of ellipsoids in given directions.
  5 %
  6 %
  7 % Description:
  8 % ------------
  9 %
 10 %    IA = MINKSUM_IA(E, L)  Computes tight internal approximating ellipsoids
 11 %                           for the geometric sum of the ellipsoids in the array E
 12 %                           in directions specified by columns of L.
 13 %                           If ellipsoids in E are n-dimensional, matrix L
 14 %                           must have dimension (n x k) where k can be arbitrarily
 15 %                           chosen. In this case, the output of the function
 16 %                           will contain k ellipsoids computed for k directions
 17 %                           specified in L.
 18 %
 19 %    Let E(q1, Q1), E(q2, Q2), ..., E(qm, Qm) be ellipsoids in R^n,
 20 %    and l - some vector in R^n. Then tight internal approximating ellipsoid E(q, Q)
 21 %    for the geometric sum E(q1, Q1) + E(q2, Q2) + ... + E(qm, Qm) in direction l,
 22 %    is such that
 23 %                 rho(l | E(q, Q)) = rho(l | (E(q1, Q1) + ... + E(qm, Qm)))
 24 %    and is defined as follows:
 25 %     q = q1 + q2 + ... + qm,
 26 %     Q = (S1 Q1^(1/2) + ... + Sm Qm^(1/2))' (S1 Q1^(1/2) + ... + Sm Qm^(1/2)),
 27 %    where S1 = I (identity), and S2, ..., Sm are orthogonal matrices such that
 28 %    vectors (S1 Q1^(1/2) l), ..., (Sm Qm^(1/2) l) are parallel.
 29 %
 30 %
 31 % Output:
 32 % -------
 33 %
 34 %    IA - array internal approximating ellipsoids.
 35 %
 36 %
 37 % See also:
 38 % ---------
 39 %
 40 %    ELLIPSOID/ELLIPSOID, MINKSUM, MINKSUM_EA, MINKDIFF, MINKDIFF_EA, MINKDIFF_IA.
 41 %
 42 
 43 %
 44 % Author:
 45 % -------
 46 %
 47 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 48 %
 49 
 50   global ellOptions;
 51 
 52   if ~isstruct(ellOptions)
 53     evalin('base', 'ellipsoids_init;');
 54   end
 55 
 56   if ~(isa(E, 'ellipsoid'))
 57     error('MINKSUM_IA: first argument must be array of ellipsoids.');
 58   end
 59 
 60   dims = dimension(E);
 61   m    = min(min(dims));
 62   n    = max(max(dims));
 63   if m ~= n
 64     error('MINKSUM_IA: ellipsoids in the array must be of the same dimension.');
 65   end
 66 
 67   [k, d] = size(L);
 68   if (k ~= n)
 69     msg = sprintf('MINKSUM_IA: second argument must be vector(s) in R^%d.', n);
 70     error(msg);
 71   end
 72 
 73   [m, n] = size(E);
 74   if (m == 1) & (n == 1)
 75     IA = E;
 76     return;
 77   end
 78 
 79   IA = [];
 80   for ii = 1:d
 81     l = L(:, ii);
 82     for i = 1:m
 83       for j = 1:n
 84         Q = E(i, j).shape;
 85         if size(Q, 1) > rank(Q)
 86           if ellOptions.verbose > 0
 87             fprintf('MINKSUM_IA: Warning! Degenerate ellipsoid.\n');
 88             fprintf('            Regularizing...\n');
 89           end
 90           Q = regularize(Q);
 91         end
 92         Q = sqrtm(Q);
 93         if (i == 1) & (j == 1)
 94           q = E(i, j).center;
 95           v = Q * l;
 96           M = Q;
 97         else
 98           q = q + E(i, j).center;
 99           T = ell_valign(v, Q*l);
100           M = M + T*Q;
101         end
102       end
103     end
104     IA = [IA ellipsoid(q, M'*M)];
105   end
106 
107   return;
  1 function EA = minksum_ea(E, L)
  2 %
  3 % MINKSUM_EA - computation of external approximating ellipsoids of the geometric
  4 %              sum of ellipsoids in given directions.
  5 %
  6 %
  7 % Description:
  8 % ------------
  9 %
 10 %    EA = MINKSUM_EA(E, L)  Computes tight external approximating ellipsoids
 11 %                           for the geometric sum of the ellipsoids in the array E
 12 %                           in directions specified by columns of L.
 13 %                           If ellipsoids in E are n-dimensional, matrix L
 14 %                           must have dimension (n x k) where k can be arbitrarily
 15 %                           chosen. In this case, the output of the function
 16 %                           will contain k ellipsoids computed for k directions
 17 %                           specified in L.
 18 %
 19 %    Let E(q1, Q1), E(q2, Q2), ..., E(qm, Qm) be ellipsoids in R^n,
 20 %    and l - some vector in R^n. Then tight external approximating ellipsoid E(q, Q)
 21 %    for the geometric sum E(q1, Q1) + E(q2, Q2) + ... + E(qm, Qm) in direction l,
 22 %    is such that
 23 %                 rho(l | E(q, Q)) = rho(l | (E(q1, Q1) + ... + E(qm, Qm)))
 24 %    and is defined as follows:
 25 %                      q = q1 + q2 + ... + qm,
 26 %                      Q = (p1 + ... + pm)((1/p1)Q1 + ... + (1/pm)Qm),
 27 %    where
 28 %          p1 = sqrt(<l, Q1l>), ..., pm = sqrt(<l, Qml>).
 29 %
 30 %
 31 % Output:
 32 % -------
 33 %
 34 %    EA - array external approximating ellipsoids.
 35 %
 36 %
 37 % See also:
 38 % ---------
 39 %
 40 %    ELLIPSOID/ELLIPSOID, MINKSUM, MINKSUM_IA, MINKDIFF, MINKDIFF_EA, MINKDIFF_IA.
 41 %
 42 
 43 %
 44 % Author:
 45 % -------
 46 %
 47 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 48 %
 49 
 50   global ellOptions;
 51 
 52   if ~isstruct(ellOptions)
 53     evalin('base', 'ellipsoids_init;');
 54   end
 55 
 56   if ~(isa(E, 'ellipsoid'))
 57     error('MINKSUM_EA: first argument must be array of ellipsoids.');
 58   end
 59 
 60   dims = dimension(E);
 61   m    = min(min(dims));
 62   n    = max(max(dims));
 63   if m ~= n
 64     error('MINKSUM_EA: ellipsoids in the array must be of the same dimension.');
 65   end
 66 
 67   [k, d] = size(L);
 68   if (k ~= n)
 69     msg = sprintf('MINKSUM_EA: second argument must be vector(s) in R^%d.', n);
 70     error(msg);
 71   end
 72 
 73   [m, n] = size(E);
 74   if (m == 1) & (n == 1)
 75     EA = E;
 76     return;
 77   end
 78 
 79   EA = [];
 80   for ii = 1:d
 81     l = L(:, ii);
 82     for i = 1:m
 83       for j = 1:n
 84         Q = E(i, j).shape;
 85         if size(Q, 1) > rank(Q)
 86           if ellOptions.verbose > 0
 87             fprintf('MINKSUM_EA: Warning! Degenerate ellipsoid.\n');
 88             fprintf('            Regularizing...\n');
 89           end
 90           Q = regularize(Q);
 91         end
 92   
 93         p = sqrt(l'*Q*l);
 94 
 95         if (i == 1) & (j == 1)
 96           q = E(i, j).center;
 97           S = (1/p) * Q;
 98           P = p;
 99         else
100           q = q + E(i, j).center;
101           S = S + ((1/p) * Q);
102           P = P + p;
103         end
104       end
105     end
106     S  = 0.5*P*(S + S');
107     EA = [EA ellipsoid(q, S)];
108   end
109 
110   return;
  1 function [y, Y] = minksum(varargin)
  2 %
  3 % MINKSUM - computes geometric (Minkowski) sum of ellipsoids in 2D or 3D.
  4 %
  5 %
  6 % Description:
  7 % ------------
  8 %
  9 % MINKSUM(EA, OPTIONS)  Computes geometric sum of ellipsoids in the array EA,
 10 %                       if 1 <= min(dimension(EA)) = max(dimension(EA)) <= 3,
 11 %                       and plots it if no output arguments are specified.
 12 %
 13 %    [y, Y] = MINKSUM(EA)  Computes geometric sum of ellipsoids in EA.
 14 %                          Here y is the center, and Y - array of
 15 %                          boundary points.
 16 %             MINKSUM(EA)  Plots geometric sum of ellipsoids in EA
 17 %                          in default (red) color.
 18 %    MINKSUM(EA, Options)  Plots geometric sum of EA using options
 19 %                          given in the Options structure.
 20 %
 21 %
 22 % Options.show_all     - if 1, displays also ellipsoids in the given array EA.
 23 % Options.newfigure    - if 1, each plot command will open a new figure window.
 24 % Options.fill         - if 1, the resulting set in 2D will be filled with color.
 25 % Options.color        - sets default colors in the form [x y z].
 26 % Options.shade = 0-1  - level of transparency (0 - transparent, 1 - opaque).
 27 %
 28 %
 29 % Output:
 30 % -------
 31 %
 32 %    y - center of the resulting set.
 33 %    Y - set of boundary points (vertices) of resulting set.
 34 %
 35 %
 36 % See also:
 37 % ---------
 38 %
 39 %    ELLIPSOID/ELLIPSOID, MINKSUM_EA, MINKSUM_IA, MINKDIFF, MINKDIFF_EA, MINKDIFF_IA.
 40 %
 41 
 42 % 
 43 % Author:
 44 % -------
 45 %
 46 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 47 %
 48 
 49   global ellOptions;
 50 
 51   if ~isstruct(ellOptions)
 52     evalin('base', 'ellipsoids_init;');
 53   end
 54 
 55   nai = nargin;
 56   E   = varargin{1};
 57   if ~isa(E, 'ellipsoid')
 58     error('MINKSUM: input argument must be an array of ellipsoids.');
 59   end
 60 
 61   ells   = varargin{1};
 62   [m, n] = size(ells);
 63   cnt    = m * n;
 64   ells   = reshape(ells, 1, cnt);
 65   dims   = dimension(ells);
 66   m      = min(dims);
 67   n      = max(dims);
 68 
 69   if m ~= n
 70     error('MINKSUM: ellipsoids must be of the same dimension.');
 71   end
 72   if n > 3
 73     error('MINKSUM: ellipsoid dimension must be not higher than 3.');
 74   end
 75 
 76   if nai > 1
 77     if isstruct(varargin{nai})
 78       Options = varargin{nai};
 79       nai     = nai - 1;
 80     else
 81         Options = [];
 82     end
 83   else
 84     Options = [];
 85   end
 86 
 87   if ~isfield(Options, 'newfigure')
 88     Options.newfigure = 0;
 89   end
 90 
 91   if ~isfield(Options, 'fill')
 92     Options.fill = 0;
 93   end
 94 
 95   if ~isfield(Options, 'show_all')
 96     Options.show_all = 0;
 97   end
 98 
 99   if ~isfield(Options, 'color')
100     Options.color = [1 0 0];
101   end
102 
103   if ~isfield(Options, 'shade')
104     Options.shade = 0.4*ones(1, cnt);
105   else
106     Options.shade = Options.shade(1, 1);
107   end
108 
109   %opts.fill            = Options.fill;
110   %opts.shade(1, 1:cnt) = Options.shade;
111 
112   if nargout == 0
113     ih = ishold;
114   end
115 
116   if (Options.show_all ~= 0) & (nargout == 0)
117     plot(ells, 'b');
118     hold on;
119     if Options.newfigure ~= 0
120       figure;
121     else
122       newplot;
123     end
124   end
125 
126   if (ellOptions.verbose > 0) & (cnt > 1)
127     if nargout == 0
128       fprintf('Computing and plotting geometric sum of %d ellipsoids...\n', cnt);
129     else
130       fprintf('Computing geometric sum of %d ellipsoids...\n', cnt);
131     end
132   end
133 	
134   clr = Options.color;
135 
136   for i = 1:cnt
137     E = ells(i);
138 
139     switch n
140       case 2,
141         if i == 1
142           Y = ellbndr_2d(E);
143           y = E.center;
144         else
145           Y = Y + ellbndr_2d(E);
146           y = y + E.center;
147         end
148         if (i == cnt) & (nargout == 0)
149           if Options.fill ~= 0
150             fill(Y(1, :), Y(2, :), clr);
151             hold on;
152           end
153           h = ell_plot(Y);
154           hold on;
155           set(h, 'Color', clr, 'LineWidth', 2);
156           h = ell_plot(y, '.');
157           set(h, 'Color', clr);
158           hasbehavior(h, 'legend', false);
159         end
160 
161       case 3,
162         if i == 1
163           Y = ellbndr_3d(E);
164           y = E.center;
165         else
166           Y = Y + ellbndr_3d(E);
167           y = y + E.center;
168         end
169         if (i == cnt) & (nargout == 0)
170           chll = convhulln(Y');
171           vs   = size(Y, 2);
172           patch('Vertices', Y', 'Faces', chll, ...
173                 'FaceVertexCData', clr(ones(1, vs), :), 'FaceColor', 'flat', ...
174                 'FaceAlpha', Options.shade(1, 1));
175           hold on;
176           shading interp;
177           lighting phong;
178           material('metal');
179           view(3);
180           %camlight('headlight','local');
181           %camlight('headlight','local');
182           %camlight('right','local');
183           %camlight('left','local');
184         end
185 
186       otherwise,
187         if i == 1
188           y       = E.center;
189           Y(1, 1) = E.center - sqrt(E.shape);
190           Y(1, 2) = E.center + sqrt(E.shape);
191         else
192           y       = y + E.center;
193           Y(1, 1) = Y(1, 1) + E.center - sqrt(E.shape);
194           Y(1, 2) = Y(1, 2) + E.center + sqrt(E.shape);
195         end
196         if (i == cnt) & (nargout == 0)
197           h = ell_plot(Y);
198           hold on;
199           set(h, 'Color', clr, 'LineWidth', 2);
200           h = ell_plot(y, '*');
201           set(h, 'Color', clr);
202         end
203 
204     end
205   end
206 
207   if nargout == 0
208     if ih == 0
209       hold off;
210     end
211   end
212 
213   if nargout == 1
214     y = Y;
215     clear Y;
216   end
217   if nargout == 0
218     clear y, Y;
219   end
220 
221   return;
 1 function IA = minkpm_ea(EE, E2, L)
 2 %
 3 % MINKPM_IA - computation of internal approximating ellipsoids
 4 %             of (E1 + E2 + ... + En) - E in given directions.
 5 %
 6 %
 7 % Description:
 8 % ------------
 9 %
10 %    IA = MINKPM_IA(EE, E, L)  Computes internal approximating ellipsoids
11 %                              of (E1 + E2 + ... + En) - E,
12 %                              where E1, E2, ..., En are ellipsoids in array EE,
13 %                              in directions specified by columns of matrix L.
14 %
15 %
16 % Output:
17 % -------
18 %
19 %    IA - array of internal approximating ellipsoids
20 %         (empty, if for all specified directions approximations cannot be computed).
21 %
22 %
23 % See also:
24 % ---------
25 %
26 %    ELLIPSOID/ELLIPSOID, MINKPM, MINKPM_EA, MINKSUM_IA, MINKDIFF_IA, MINKMP_IA.
27 %
28 
29 %
30 % Author:
31 % -------
32 %
33 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
34 %
35 
36   global ellOptions;
37 
38   if ~isstruct(ellOptions)
39     evalin('base', 'ellipsoids_init;');
40   end
41 
42   if ~(isa(EE, 'ellipsoid')) | ~(isa(E2, 'ellipsoid'))
43     error('MINKPM_IA: first and second arguments must be ellipsoids.');
44   end
45 
46   [m, n] = size(E2);
47   if (m ~= 1) | (n ~= 1)
48     error('MINKPM_IA: second argument must be single ellipsoid.');
49   end
50 
51   k  = size(L, 1);
52   n  = dimension(E2);
53   mn = min(min(dimension(EE)));
54   mx = max(max(dimension(EE)));
55   if (mn ~= mx) | (mn ~= n)
56     error('MINKPM_IA: all ellipsoids must be of the same dimension.');
57   end
58   if n ~= k
59     error('MINKPM_IA: dimension of the direction vectors must be the same as dimension of ellipsoids.');
60   end
61 
62   N                  = size(L, 2);
63   IA                 = [];
64   ES                 = minksum_ia(EE, L);
65   vrb                = ellOptions.verbose;
66   ellOptions.verbose = 0;
67 
68   for i = 1:N
69     E = ES(i);
70     l = L(:, i);
71     if isbigger(E, E2)
72       if ~isbaddirection(E, E2, l)
73         IA = [IA minkdiff_ia(E, E2, l)];
74       end
75     end
76   end
77 
78   ellOptions.verbose = vrb;
79 
80   if isempty(IA)
81     if ellOptions.verbose > 0
82       fprintf('MINKPM_IA: cannot compute internal approximation for any\n');
83       fprintf('           of the specified directions.\n');
84     end
85   end
86 
87   return;
 1 function EA = minkpm_ea(EE, E2, L)
 2 %
 3 % MINKPM_EA - computation of external approximating ellipsoids
 4 %             of (E1 + E2 + ... + En) - E in given directions.
 5 %
 6 %
 7 % Description:
 8 % ------------
 9 %
10 %    EA = MINKPM_EA(EE, E, L)  Computes external approximating ellipsoids
11 %                              of (E1 + E2 + ... + En) - E,
12 %                              where E1, E2, ..., En are ellipsoids in array EE,
13 %                              in directions specified by columns of matrix L.
14 %
15 %
16 % Output:
17 % -------
18 %
19 %    EA - array of external approximating ellipsoids
20 %         (empty, if for all specified directions approximations cannot be computed).
21 %
22 %
23 % See also:
24 % ---------
25 %
26 %    ELLIPSOID/ELLIPSOID, MINKPM, MINKPM_IA, MINKSUM_EA, MINKDIFF_EA, MINKMP_EA.
27 %
28 
29 %
30 % Author:
31 % -------
32 %
33 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
34 %
35 
36   global ellOptions;
37 
38   if ~isstruct(ellOptions)
39     evalin('base', 'ellipsoids_init;');
40   end
41 
42   if ~(isa(EE, 'ellipsoid')) | ~(isa(E2, 'ellipsoid'))
43     error('MINKPM_EA: first and second arguments must be ellipsoids.');
44   end
45 
46   [m, n] = size(E2);
47   if (m ~= 1) | (n ~= 1)
48     error('MINKPM_EA: second argument must be single ellipsoid.');
49   end
50 
51   k  = size(L, 1);
52   n  = dimension(E2);
53   mn = min(min(dimension(EE)));
54   mx = max(max(dimension(EE)));
55   if (mn ~= mx) | (mn ~= n)
56     error('MINKPM_EA: all ellipsoids must be of the same dimension.');
57   end
58   if n ~= k
59     error('MINKPM_EA: dimension of the direction vectors must be the same as dimension of ellipsoids.');
60   end
61 
62   N                  = size(L, 2);
63   EA                 = [];
64   vrb                = ellOptions.verbose;
65   ellOptions.verbose = 0;
66   
67   % sanity check: the approximated set should be nonempty
68   for i = 1:N
69     [U, S, V] = svd(L(:, i));
70     ET        = minksum_ea(EE, U);
71     if min(ET > E2) < 1
72       if vrb > 0
73         fprintf('MINKPM_EA: the resulting set is empty.\n');
74       end
75       ellOptions.verbose = vrb;
76       return;
77     end
78   end
79 
80   ES = minksum_ea(EE, L);
81 
82   for i = 1:N
83     E = ES(i);
84     l = L(:, i);
85     if ~isbaddirection(E, E2, l)
86       EA = [EA minkdiff_ea(E, E2, l)];
87     end
88   end
89   
90   ellOptions.verbose = vrb;
91 
92   if isempty(EA)
93     if ellOptions.verbose > 0
94       fprintf('MINKPM_EA: cannot compute external approximation for any\n');
95       fprintf('           of the specified directions.\n');
96     end
97   end
98 
99   return;
  1 function [y, Y] = minkpm(varargin)
  2 %
  3 % MINKPM - computes and plots geometric (Minkowski) difference of the geometric
  4 %          sum of ellipsoids and a single ellipsoid in 2D or 3D:
  5 %
  6 %          (E1 + E2 + ... + En) - E
  7 %
  8 %
  9 % Description:
 10 % ------------
 11 %
 12 % MINKPM(EA, E, OPTIONS)  Computes geometric difference of the geometric sum
 13 %                         of ellipsoids in EA and ellipsoid E,
 14 %                         if 1 <= dimension(EA) = dimension(E) <= 3,
 15 %                         and plots it if no output arguments are specified.
 16 %
 17 %    [y, Y] = MINKPM(EA, E)  Computes (geometric sum of ellipsoids in EA) - E.
 18 %                            Here y is the center, and Y - array of boundary points.
 19 %             MINKPM(EA, E)  Plots (geometric sum of ellipsoids in EA) - E
 20 %                            in default (red) color.
 21 %    MINKPM(EA, E, Options)  Plots (geometric sum of ellipsoids in EA) - E
 22 %                            using options given in the Options structure.
 23 %
 24 % Options.show_all     - if 1, displays also ellipsoids E1 and E2.
 25 % Options.newfigure    - if 1, each plot command will open a new figure window.
 26 % Options.fill         - if 1, the resulting set in 2D will be filled with color.
 27 % Options.color        - sets default colors in the form [x y z].
 28 % Options.shade = 0-1  - level of transparency (0 - transparent, 1 - opaque).
 29 %
 30 %
 31 % Output:
 32 % -------
 33 %
 34 %    y - center of the resulting set.
 35 %    Y - set of boundary points (vertices) of resulting set.
 36 %
 37 %
 38 % See also:
 39 % ---------
 40 %
 41 %    ELLIPSOID/ELLIPSOID, MINKSUM, MINKDIFF, MINKMP.
 42 %
 43 
 44 % 
 45 % Author:
 46 % -------
 47 %
 48 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 49 %
 50 
 51   global ellOptions;
 52 
 53   if ~isstruct(ellOptions)
 54     evalin('base', 'ellipsoids_init;');
 55   end
 56 
 57   if nargin < 2
 58     error('MINKPM: first and second arguments must be ellipsoids.');
 59   end
 60 
 61   EE = varargin{1};
 62   E2 = varargin{2};
 63 
 64   if ~(isa(EE, 'ellipsoid')) | ~(isa(E2, 'ellipsoid'))
 65     error('MINKPM: first and second arguments must be ellipsoids.');
 66   end
 67 
 68   [m, n] = size(E2);
 69   if (m ~= 1) | (n ~= 1)
 70     error('MINKPM: second argument must be single ellipsoid.');
 71   end
 72 
 73   dims = dimension(EE);
 74   k    = min(min(dims));
 75   m    = max(max(dims));
 76   n    = dimension(E2);
 77   if (k ~= m) | (m ~= n)
 78     error('MINKPM: all ellipsoids must be of the same dimension.');
 79   end
 80   if n > 3
 81     error('MINKPM: ellipsoid dimension must be not higher than 3.');
 82   end
 83 
 84   switch n
 85     case 2,
 86       phi = linspace(0, 2*pi, ellOptions.plot2d_grid);
 87       L   = [cos(phi); sin(phi)];
 88 
 89     case 3,
 90       M   = ellOptions.plot3d_grid/2;
 91       N   = M/2;
 92       psy = linspace(0, pi, N);
 93       phi = linspace(0, 2*pi, M);
 94       L   = [];
 95       for i = 2:(N - 1)
 96         arr = cos(psy(i))*ones(1, M);
 97         L   = [L [cos(phi)*sin(psy(i)); sin(phi)*sin(psy(i)); arr]];
 98       end
 99       
100     otherwise,
101       L = [-1 1];
102 
103   end
104 
105   vrb                = ellOptions.verbose;
106   ellOptions.verbose = 0;
107   EA                 = minksum_ea(EE, L);
108   ellOptions.verbose = vrb;
109   
110   if min(EA > E2) == 0
111     switch nargout
112       case 0,
113         fprintf('The resulting set is empty.');
114         return;
115 
116       case 1,
117         y = [];
118         return;
119 
120       otherwise,
121         y = [];
122         Y = [];
123         return;
124 
125     end
126   end
127 
128   if nargin > 2
129     if isstruct(varargin{3})
130       Options = varargin{3};
131     else
132       Options = [];
133     end
134   else
135     Options = [];
136   end
137 
138   if ~isfield(Options, 'newfigure')
139     Options.newfigure = 0;
140   end
141 
142   if ~isfield(Options, 'fill')
143     Options.fill = 0;
144   end
145 
146   if ~isfield(Options, 'show_all')
147     Options.show_all = 0;
148   end
149 
150   if ~isfield(Options, 'color')
151     Options.color = [1 0 0];
152   end
153 
154   if ~isfield(Options, 'shade')
155     Options.shade = 0.4;
156   else
157     Options.shade = Options.shade(1, 1);
158   end
159 
160   clr  = Options.color;
161 
162   if nargout == 0
163     ih = ishold;
164   end
165 
166   if (Options.show_all ~= 0) & (nargout == 0)
167     plot(EE, 'b', E2, 'k');
168     hold on;
169     if Options.newfigure ~= 0
170       figure;
171     else
172       newplot;
173     end
174   end
175 
176   if ellOptions.verbose > 0
177     if nargout == 0
178       fprintf('Computing and plotting (sum(E_i) - E) ...\n');
179     else
180       fprintf('Computing (sum(E_i) - E) ...\n');
181     end
182   end
183 
184   y                  = EA(1).center - E2.center;
185   Y                  = [];
186 %  EF                 = [];
187 %  LL                 = [];
188   N                  = size(L, 2);
189   ellOptions.verbose = 0;
190 
191 %  for i = 1:N
192 %    l = L(:, i);
193 %    E = EA(i);
194 %    if ~isbaddirection(E, E2, l)
195 %      EF = [EF minkdiff_ea(E, E2, l)];
196 %      LL = [LL l];
197 %    end
198 %  end
199 %
200 %  M = size(EF, 2);
201 %  
202 %  for i = 1:N
203 %    l    = L(:, i);
204 %    mval = 0;
205 %
206 %    for j = 1:M
207 %      Q = parameters(EF(j));
208 %      v = l' * ell_inv(Q) * l;
209 %      if v > mval
210 %        mval = v;
211 %      end
212 %    end
213 %    
214 %    if mval > 0
215 %      Y = [Y ((l/sqrt(mval))+y)];
216 %    end
217 %  end
218 %
219 %  if isempty(Y)
220 %    Y = y;
221 %  end
222 
223   switch n
224     case 2,
225       EF = [];
226       LL = [];
227       for i = 1:N
228         l = L(:, i);
229         E = EA(i);
230         if ~isbaddirection(E, E2, l)
231           EF = [EF minkdiff_ea(E, E2, l)];
232           LL = [LL l];
233         end
234       end
235       M = size(EF, 2);
236       for i = 1:N
237         l    = L(:, i);
238         mval = 0;
239         for j = 1:M
240           Q = parameters(EF(j));
241           v = l' * ell_inv(Q) * l;
242           if v > mval
243             mval = v;
244           end
245         end
246         if mval > 0
247           Y = [Y ((l/sqrt(mval))+y)];
248         end
249       end
250       if isempty(Y)
251         Y = y;
252       end
253       Y = [Y Y(:, 1)];
254       if nargout == 0
255         if Options.fill ~= 0
256           fill(Y(1, :), Y(2, :), clr);
257           hold on;
258         end
259         h = ell_plot(Y);
260         hold on;
261         set(h, 'Color', clr, 'LineWidth', 2);
262         h = ell_plot(y, '.');
263         set(h, 'Color', clr);
264       end
265 
266     case 3,
267       for i = 1:N
268         l = L(:, i);
269         E = EA(i);
270         if ~isbaddirection(E, E2, l)
271           I = minksum_ia(EE, l);
272           if isbigger(I, E2);
273             if ~isbaddirection(I, E2, l)
274               [r, x] = rho(minkdiff_ea(E, E2, l), l);
275               Y      = [Y x];
276             end
277           end
278         end
279       end
280       if isempty(Y)
281         Y = y;
282       end
283       if nargout == 0
284         vs = size(Y, 2);
285         if vs > 1
286           chll = convhulln(Y');
287           patch('Vertices', Y', 'Faces', chll, ...
288                 'FaceVertexCData', clr(ones(1, vs), :), 'FaceColor', 'flat', ...
289                 'FaceAlpha', Options.shade(1, 1));
290         else
291           h = ell_plot(y, '*');
292           set(h, 'Color', clr);
293         end
294         hold on;
295         shading interp;
296         lighting phong;
297         material('metal');
298         view(3);
299         %camlight('headlight','local');
300         %camlight('headlight','local');
301         %camlight('right','local');
302         %camlight('left','local');
303       end
304 
305     otherwise,
306       Y       = [y y];
307       Y(1, 1) = EA(1).center - E2.center + sqrt(E2.shape) - sqrt(EA(1).shape);
308       Y(1, 2) = EA(1).center - E2.center + sqrt(EA(1).shape) - sqrt(E2.shape);
309       if nargout == 0
310         h = ell_plot(Y);
311         hold on;
312         set(h, 'Color', clr, 'LineWidth', 2);
313         h = ell_plot(y, '*');
314         set(h, 'Color', clr);
315       end
316 
317   end
318 
319   ellOptions.verbose = vrb;
320 
321   if nargout == 0
322     if ih == 0
323       hold off;
324     end
325   end
326 
327   if nargout == 1
328     y = Y;
329   end
330   if nargout == 0
331     clear y, Y;
332   end
333 
334   return;
 1 function IA = minkmp_ia(E1, E2, EE, L)
 2 %
 3 % MINKMP_IA - computation of internal approximating ellipsoids
 4 %             of (E0 - E) + (E1 + ... + En) in given directions.
 5 %
 6 %
 7 % Description:
 8 % ------------
 9 %
10 % IA = MINKMP_IA(E0, E, EE, L)  Computes internal approximating ellipsoids
11 %                               of (E0 - E) + (E1 + E2 + ... + En),
12 %                               where E1, E2, ..., En are ellipsoids in array EE,
13 %                               in directions specified by columns of matrix L.
14 %
15 %
16 % Output:
17 % -------
18 %
19 %    IA - array of internal approximating ellipsoids
20 %         (empty, if for all specified directions approximations cannot be computed).
21 %
22 %
23 % See also:
24 % ---------
25 %
26 %    ELLIPSOID/ELLIPSOID, MINKMP, MINKMP_EA, MINKSUM_IA, MINKDIFF_IA, MINKPM_IA.
27 %
28 
29 %
30 % Author:
31 % -------
32 %
33 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
34 %
35 
36   global ellOptions;
37 
38   if ~isstruct(ellOptions)
39     evalin('base', 'ellipsoids_init;');
40   end
41 
42   if ~(isa(EE, 'ellipsoid')) | ~(isa(E2, 'ellipsoid')) | ~(isa(E1, 'ellipsoid'))
43     error('MINKMP_IA: first, second and third arguments must be ellipsoids.');
44   end
45 
46   [k, l] = size(E1);
47   [m, n] = size(E2);
48   if (k ~= 1) | (l ~= 1) | (m ~= 1) | (n ~= 1)
49     error('MINKMP_IA: first and second arguments must be single ellipsoids.');
50   end
51 
52   k  = size(L, 1);
53   m  = dimension(E1);
54   n  = dimension(E2);
55   mn = min(min(dimension(EE)));
56   mx = max(max(dimension(EE)));
57   if (mn ~= mx) | (mn ~= n) | (m ~= n)
58     error('MINKMP_IA: all ellipsoids must be of the same dimension.');
59   end
60   if n ~= k
61     error('MINKMP_IA: dimension of the direction vectors must be the same as dimension of ellipsoids.');
62   end
63 
64   IA = [];
65 
66   if ~isbigger(E1, E2)
67     if ellOptions.verbose > 0
68       fprintf('MINKMP_IA: the resulting set is empty.\n');
69     end
70     return;
71   end
72 
73   LL                 = [];
74   N                  = size(L, 2);
75   [m, n]             = size(EE);
76   EE                 = reshape(EE, 1, m*n);
77   vrb                = ellOptions.verbose;
78   ellOptions.verbose = 0;
79 
80 
81   for i = 1:N
82     l = L(:, i);
83     if ~isbaddirection(E1, E2, l)
84       LL = [LL l];
85       IA = [IA minksum_ia([minkdiff_ia(E1, E2, l) EE], l)];
86     end
87   end
88   
89   ellOptions.verbose = vrb;
90 
91   if isempty(IA)
92     if ellOptions.verbose > 0
93       fprintf('MINKMP_IA: cannot compute external approximation for any\n');
94       fprintf('           of the specified directions.\n');
95     end
96   end
97 
98   return;
 1 function EA = minkmp_ea(E1, E2, EE, L)
 2 %
 3 % MINKMP_EA - computation of external approximating ellipsoids
 4 %             of (E0 - E) + (E1 + ... + En) in given directions.
 5 %
 6 %
 7 % Description:
 8 % ------------
 9 %
10 % EA = MINKMP_EA(E0, E, EE, L)  Computes external approximating ellipsoids
11 %                               of (E0 - E) + (E1 + E2 + ... + En),
12 %                               where E1, E2, ..., En are ellipsoids in array EE,
13 %                               in directions specified by columns of matrix L.
14 %
15 %
16 % Output:
17 % -------
18 %
19 %    EA - array of external approximating ellipsoids
20 %         (empty, if for all specified directions approximations cannot be computed).
21 %
22 %
23 % See also:
24 % ---------
25 %
26 %    ELLIPSOID/ELLIPSOID, MINKMP, MINKMP_IA, MINKSUM_EA, MINKDIFF_EA, MINKPM_EA.
27 %
28 
29 %
30 % Author:
31 % -------
32 %
33 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
34 %
35 
36   global ellOptions;
37 
38   if ~isstruct(ellOptions)
39     evalin('base', 'ellipsoids_init;');
40   end
41 
42   if ~(isa(EE, 'ellipsoid')) | ~(isa(E2, 'ellipsoid')) | ~(isa(E1, 'ellipsoid'))
43     error('MINKMP_EA: first, second and third arguments must be ellipsoids.');
44   end
45 
46   [k, l] = size(E1);
47   [m, n] = size(E2);
48   if (k ~= 1) | (l ~= 1) | (m ~= 1) | (n ~= 1)
49     error('MINKMP_EA: first and second arguments must be single ellipsoids.');
50   end
51 
52   k  = size(L, 1);
53   m  = dimension(E1);
54   n  = dimension(E2);
55   mn = min(min(dimension(EE)));
56   mx = max(max(dimension(EE)));
57   if (mn ~= mx) | (mn ~= n) | (m ~= n)
58     error('MINKMP_EA: all ellipsoids must be of the same dimension.');
59   end
60   if n ~= k
61     error('MINKMP_EA: dimension of the direction vectors must be the same as dimension of ellipsoids.');
62   end
63 
64   EA = [];
65 
66   if ~isbigger(E1, E2)
67     if ellOptions.verbose > 0
68       fprintf('MINKMP_EA: the resulting set is empty.\n');
69     end
70     return;
71   end
72 
73   LL                 = [];
74   N                  = size(L, 2);
75   [m, n]             = size(EE);
76   EE                 = reshape(EE, 1, m*n);
77   vrb                = ellOptions.verbose;
78   ellOptions.verbose = 0;
79 
80   for i = 1:N
81     l = L(:, i);
82     if ~isbaddirection(E1, E2, l)
83       LL = [LL l];
84       EA = [EA minksum_ea([minkdiff_ea(E1, E2, l) EE], l)];
85     end
86   end
87 
88   ellOptions.verbose = vrb;
89 
90   if isempty(EA)
91     if ellOptions.verbose > 0
92       fprintf('MINKMP_EA: cannot compute external approximation for any\n');
93       fprintf('           of the specified directions.\n');
94     end
95   end
96 
97   return;
  1 function [y, Y] = minkmp(varargin)
  2 %
  3 % MINKMP - computes and plots geometric (Minkowski) sum of the geometric difference
  4 %          of two ellipsoids and the geometric sum of n ellipsoids in 2D or 3D:
  5 %
  6 %          (E0 - E) + (E1 + E2 + ... + En)
  7 %
  8 %
  9 % Description:
 10 % ------------
 11 %
 12 % MINKMP(E0, E, EE, OPTIONS)  Computes geometric sum of the geometric difference
 13 %                             of two ellipsoids E0 - E and the geometric sum of
 14 %                             ellipsoids in the ellipsoidal array EE, if
 15 %                             1 <= dimension(E0) = dimension(E) = dimension(EE) <= 3,
 16 %                             and plots it if no output arguments are specified.
 17 %
 18 %    [y, Y] = MINKMP(E0, E, EE)  Computes (E0 - E) + (geometric sum of ellipsoids in EE).
 19 %                                Here y is the center, and Y - array of boundary points.
 20 %             MINKMP(E0, E, EE)  Plots (E0 - E) + (geometric sum of ellipsoids in EE)
 21 %                                in default (red) color.
 22 %    MINKMP(E0, E, EE, Options)  Plots (E0 - E) + (geometric sum of ellipsoids in EE)
 23 %                                using options given in the Options structure.
 24 %
 25 % Options.show_all     - if 1, displays also ellipsoids E1 and E2.
 26 % Options.newfigure    - if 1, each plot command will open a new figure window.
 27 % Options.fill         - if 1, the resulting set in 2D will be filled with color.
 28 % Options.color        - sets default colors in the form [x y z].
 29 % Options.shade = 0-1  - level of transparency (0 - transparent, 1 - opaque).
 30 %
 31 %
 32 % Output:
 33 % -------
 34 %
 35 %    y - center of the resulting set.
 36 %    Y - set of boundary points (vertices) of resulting set.
 37 %
 38 %
 39 % See also:
 40 % ---------
 41 %
 42 %    ELLIPSOID/ELLIPSOID, MINKSUM, MINKDIFF, MINKPM.
 43 %
 44 
 45 % 
 46 % Author:
 47 % -------
 48 %
 49 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 50 %
 51 
 52   global ellOptions;
 53 
 54   if ~isstruct(ellOptions)
 55     evalin('base', 'ellipsoids_init;');
 56   end
 57 
 58   if nargin < 3
 59     error('MINKMP: first, second and third arguments must be ellipsoids.');
 60   end
 61 
 62   E1 = varargin{1};
 63   E2 = varargin{2};
 64   EE = varargin{3};
 65 
 66   if ~(isa(E1, 'ellipsoid')) | ~(isa(E2, 'ellipsoid')) | ~(isa(EE, 'ellipsoid'))
 67     error('MINKMP: first, second and third arguments must be ellipsoids.');
 68   end
 69 
 70   [k, l] = size(E1);
 71   [m, n] = size(E2);
 72   if (k ~= 1) | (l ~= 1) | (m ~= 1) | (n ~= 1)
 73     error('MINKMP: first and second arguments must be single ellipsoid.');
 74   end
 75 
 76   m    = dimension(E1);
 77   n    = dimension(E2);
 78   dims = dimension(EE);
 79   mn   = min(min(dims));
 80   mx   = max(max(dims));
 81   if (mn ~= mx) | (mn ~= n) | (m ~= n)
 82     error('MINKMP: all ellipsoids must be of the same dimension.');
 83   end
 84   if n > 3
 85     error('MINKMP: ellipsoid dimension must be not higher than 3.');
 86   end
 87 
 88   if ~isbigger(E1, E2)
 89     switch nargout
 90       case 0,
 91         fprintf('The resulting set is empty.');
 92         return;
 93 
 94       case 1,
 95         y = [];
 96         return;
 97 
 98       otherwise,
 99         y = [];
100         Y = [];
101         return;
102 
103     end
104   end
105 
106   switch n
107     case 2,
108       phi = linspace(0, 2*pi, ellOptions.plot2d_grid);
109       L   = [cos(phi); sin(phi)];
110 
111     case 3,
112       M   = ellOptions.plot3d_grid/2;
113       N   = M/2;
114       psy = linspace(0, pi, N);
115       phi = linspace(0, 2*pi, M);
116       L   = [];
117       for i = 2:(N - 1)
118         arr = cos(psy(i))*ones(1, M);
119         L   = [L [cos(phi)*sin(psy(i)); sin(phi)*sin(psy(i)); arr]];
120       end
121       
122     otherwise,
123       L = [-1 1];
124 
125   end
126 
127   if nargin > 3
128     if isstruct(varargin{4})
129       Options = varargin{4};
130     else
131       Options = [];
132     end
133   else
134     Options = [];
135   end
136 
137   if ~isfield(Options, 'newfigure')
138     Options.newfigure = 0;
139   end
140 
141   if ~isfield(Options, 'fill')
142     Options.fill = 0;
143   end
144 
145   if ~isfield(Options, 'show_all')
146     Options.show_all = 0;
147   end
148 
149   if ~isfield(Options, 'color')
150     Options.color = [1 0 0];
151   end
152 
153   if ~isfield(Options, 'shade')
154     Options.shade = 0.4;
155   else
156     Options.shade = Options.shade(1, 1);
157   end
158 
159   clr  = Options.color;
160 
161   if nargout == 0
162     ih = ishold;
163   end
164 
165   if (Options.show_all ~= 0) & (nargout == 0)
166     plot(EE, 'b', E2, 'k', E1, 'g');
167     hold on;
168     if Options.newfigure ~= 0
169       figure;
170     else
171       newplot;
172     end
173   end
174 
175   if ellOptions.verbose > 0
176     if nargout == 0
177       fprintf('Computing and plotting (E0 - E) + sum(E_i) ...\n');
178     else
179       fprintf('Computing (E0 - E) + sum(E_i) ...\n');
180     end
181   end
182 
183   vrb                = ellOptions.verbose;
184   ellOptions.verbose = 0;
185   [x, X]             = minksum(EE);
186   ellOptions.verbose = vrb;
187   y                  = E1.center - E2.center + x;
188   Y                  = [];
189   N                  = size(L, 2);
190   bd                 = isbaddirection(E1, E2, L);
191   bd                 = find(bd == 0);
192 
193   if isempty(bd)
194     x = E1.center - E2.center;
195   end
196 
197   l       = L(:, bd(1));
198   [r, x1] = rho(E1, l);
199   [r, x2] = rho(E2, l);
200   x       = x1 - x2;
201 
202   for i = 1:N
203     l = L(:, i);
204     if ~isbaddirection(E1, E2, l)
205       [r, x1] = rho(E1, l);
206       [r, x2] = rho(E2, l);
207       x       = x1 - x2;
208     end
209     Y = [Y (x+X(:, i))];
210   end
211 
212   switch n
213     case 2,
214       Y = [Y Y(:, 1)];
215 
216       if nargout == 0
217         if Options.fill ~= 0
218           fill(Y(1, :), Y(2, :), clr);
219           hold on;
220         end
221         h = ell_plot(Y);
222         hold on;
223         set(h, 'Color', clr, 'LineWidth', 2);
224         h = ell_plot(y, '.');
225         set(h, 'Color', clr);
226       end
227 
228     case 3,
229       if nargout == 0
230         vs = size(Y, 2);
231         if vs > 1
232           chll = convhulln(Y');
233           patch('Vertices', Y', 'Faces', chll, ...
234                 'FaceVertexCData', clr(ones(1, vs), :), 'FaceColor', 'flat', ...
235                 'FaceAlpha', Options.shade(1, 1));
236         else
237           h = ell_plot(y, '*');
238           set(h, 'Color', clr);
239         end
240         hold on;
241         shading interp;
242         lighting phong;
243         material('metal');
244         view(3);
245         %camlight('headlight','local');
246         %camlight('headlight','local');
247         %camlight('right','local');
248         %camlight('left','local');
249       end
250 
251     otherwise,
252       Y       = [y y];
253       Y(1, 1) = EA(1).center - E2.center + sqrt(E2.shape) - sqrt(EA(1).shape);
254       Y(1, 2) = EA(1).center - E2.center + sqrt(EA(1).shape) - sqrt(E2.shape);
255       if nargout == 0
256         h = ell_plot(Y);
257         hold on;
258         set(h, 'Color', clr, 'LineWidth', 2);
259         h = ell_plot(y, '*');
260         set(h, 'Color', clr);
261       end
262 
263   end
264 
265   if nargout == 0
266     if ih == 0
267       hold off;
268     end
269   end
270 
271   if nargout == 1
272     y = Y;
273   end
274   if nargout == 0
275     clear y, Y;
276   end
277 
278   return;
  1 function IA = minkdiff_ia(E1, E2, L)
  2 %
  3 % MINKDIFF_IA - computation of internal approximating ellipsoids of the geometric
  4 %               difference of two ellipsoids in given directions.
  5 %
  6 %
  7 % Description:
  8 % ------------
  9 %
 10 %    IA = MINKDIFF_IA(E1, E2, L)  Computes internal approximating ellipsoids
 11 %                                 of the geometric difference of two ellipsoids E1 - E2
 12 %                                 in directions specified by columns of matrix L.
 13 %
 14 %    First condition for the approximations to be computed, is that ellipsoid
 15 %    E1 must be bigger than ellipsoid E2 in the sense that if they had the same
 16 %    center, E2 would be contained inside E1. Otherwise, the geometric difference
 17 %    E1 - E2 is an empty set.
 18 %    Second condition for the approximation in the given direction l to exist,
 19 %    is the following. Given 
 20 %                   P = sqrt(<l, Q1 l>)/sqrt(<l, Q2 l>)
 21 %    where Q1 is the shape matrix of ellipsoid E1, and Q2 - shape matrix of E2,
 22 %    and R being minimal root of the equation
 23 %                   det(Q1 - R Q2) = 0,
 24 %    parameter P should be less than R.
 25 %    If these two conditions are satisfied, then internal approximating
 26 %    ellipsoid for the geometric difference E1 - E2 in the direction l
 27 %    is defined by its shape matrix
 28 %                 Q = (1 - (1/P)) Q1 + (1 - P) Q2
 29 %    and its center
 30 %                 q = q1 - q2,
 31 %    where q1 is center of E1 and q2 - center of E2.
 32 %
 33 %
 34 % Output:
 35 % -------
 36 %
 37 %    IA - array of internal approximating ellipsoids
 38 %         (empty, if for all specified directions approximations cannot be computed).
 39 %
 40 %
 41 % See also:
 42 % ---------
 43 %
 44 %    ELLIPSOID/ELLIPSOID, MINKDIFF_EA, MINKDIFF, ISBIGGER, ISBADDIRECTION,
 45 %                         MINKSUM_EA, MINKSUM_IA.
 46 %
 47 
 48 %
 49 % Author:
 50 % -------
 51 %
 52 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 53 %
 54 
 55   global ellOptions;
 56 
 57   if ~isstruct(ellOptions)
 58     evalin('base', 'ellipsoids_init;');
 59   end
 60 
 61   if ~(isa(E1, 'ellipsoid')) | ~(isa(E2, 'ellipsoid'))
 62     error('MINKDIFF_IA: first and second arguments must be single ellipsoids.');
 63   end
 64 
 65   [k, l] = size(E1);
 66   [m, n] = size(E2);
 67   if (k ~= 1) | (l ~= 1) | (m ~= 1) | (n ~= 1)
 68     error('MINKDIFF_IA: first and second arguments must be single ellipsoids.');
 69   end
 70 
 71   IA = [];
 72 
 73   if isbigger(E1, E2) == 0
 74     if ellOptions.verbose > 0
 75       fprintf('MINKDIFF_IA: geometric difference of these two ellipsoids is empty set.\n');
 76     end
 77     return;
 78   end
 79 
 80   k = size(L, 1);
 81   n = dimension(E1);
 82   if k ~= n
 83     error('MINKDIFF_IA: dimension of the direction vectors must be the same as dimension of ellipsoids.');
 84   end
 85   q  = E1.center - E2.center;
 86   Q1 = E1.shape;
 87   if rank(Q1) < size(Q1, 1)
 88     Q1 = regularize(Q1);
 89   end
 90   Q2 = E2.shape;
 91   if rank(Q2) < size(Q2, 1)
 92     Q2 = regularize(Q2);
 93   end
 94   L  = rm_bad_directions(Q1, Q2, L);
 95   m  = size(L, 2);
 96   if m < 1
 97     if ellOptions.verbose > 0
 98       fprintf('MINKDIFF_IA: cannot compute internal approximation for any\n');
 99       fprintf('             of the specified directions.\n');
100     end
101     return;
102   end
103   for i = 1:m
104     l  = L(:, i);
105     p  = (sqrt(l'*Q1*l))/(sqrt(l'*Q2*l));
106     Q  = (1 - (1/p))*Q1 + (1 - p)*Q2;
107     IA = [IA ellipsoid(q, Q)];
108   end 
109 
110   return;
  1 function EA = minkdiff_ea(E1, E2, L)
  2 %
  3 % MINKDIFF_EA - computation of external approximating ellipsoids of the geometric
  4 %               difference of two ellipsoids in given directions.
  5 %
  6 %
  7 % Description:
  8 % ------------
  9 %
 10 %    EA = MINKDIFF_EA(E1, E2, L)  Computes external approximating ellipsoids
 11 %                                 of the geometric difference of two ellipsoids E1 - E2
 12 %                                 in directions specified by columns of matrix L.
 13 %
 14 %    First condition for the approximations to be computed, is that ellipsoid
 15 %    E1 must be bigger than ellipsoid E2 in the sense that if they had the same
 16 %    center, E2 would be contained inside E1. Otherwise, the geometric difference
 17 %    E1 - E2 is an empty set.
 18 %    Second condition for the approximation in the given direction l to exist,
 19 %    is the following. Given 
 20 %                   P = sqrt(<l, Q1 l>)/sqrt(<l, Q2 l>)
 21 %    where Q1 is the shape matrix of ellipsoid E1, and Q2 - shape matrix of E2,
 22 %    and R being minimal root of the equation
 23 %                   det(Q1 - R Q2) = 0,
 24 %    parameter P should be less than R.
 25 %    If both of these conditions are satisfied, then external approximating
 26 %    ellipsoid is defined by its shape matrix
 27 %             Q = (Q1^(1/2) + S Q2^(1/2))' (Q1^(1/2) + S Q2^(1/2)),
 28 %    where S is orthogonal matrix such that vectors Q1^(1/2) l and S Q2^(1/2) l
 29 %    are parallel, and its center
 30 %             q = q1 - q2,
 31 %    where q1 is center of ellipsoid E1 and q2 - center of E2.
 32 %
 33 %
 34 % Output:
 35 % -------
 36 %
 37 %    EA - array of external approximating ellipsoids
 38 %         (empty, if for all specified directions approximations cannot be computed).
 39 %
 40 %
 41 % See also:
 42 % ---------
 43 %
 44 %    ELLIPSOID/ELLIPSOID, MINKDIFF_IA, MINKDIFF, ISBIGGER, ISBADDIRECTION,
 45 %                         MINKSUM_EA, MINKSUM_IA.
 46 %
 47 
 48 %
 49 % Author:
 50 % -------
 51 %
 52 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 53 %
 54 
 55   global ellOptions;
 56 
 57   if ~isstruct(ellOptions)
 58     evalin('base', 'ellipsoids_init;');
 59   end
 60 
 61   if ~(isa(E1, 'ellipsoid')) | ~(isa(E2, 'ellipsoid'))
 62     error('MINKDIFF_EA: first and second arguments must be single ellipsoids.');
 63   end
 64 
 65   [k, l] = size(E1);
 66   [m, n] = size(E2);
 67   if (k ~= 1) | (l ~= 1) | (m ~= 1) | (n ~= 1)
 68     error('MINKDIFF_EA: first and second arguments must be single ellipsoids.');
 69   end
 70 
 71   EA = [];
 72 
 73   if isbigger(E1, E2) == 0
 74     if ellOptions.verbose > 0
 75       fprintf('MINKDIFF_EA: geometric difference of these two ellipsoids is empty set.\n');
 76     end
 77     return;
 78   end
 79 
 80   k = size(L, 1);
 81   n = dimension(E1);
 82   if k ~= n
 83     error('MINKDIFF_EA: dimension of the direction vectors must be the same as dimension of ellipsoids.');
 84   end
 85   q  = E1.center - E2.center;
 86   Q1 = E1.shape;
 87   Q2 = E2.shape;
 88   L  = rm_bad_directions(Q1, Q2, L);
 89   m  = size(L, 2);
 90   if m < 1
 91     if ellOptions.verbose > 0
 92       fprintf('MINKDIFF_EA: cannot compute external approximation for any\n');
 93       fprintf('             of the specified directions.\n');
 94     end
 95     return;
 96   end
 97   if rank(Q1) < size(Q1, 1)
 98     Q1 = regularize(Q1);
 99   end
100   if rank(Q2) < size(Q2, 1)
101     Q2 = regularize(Q2);
102   end
103 
104   Q1 = sqrtm(Q1);
105   Q2 = sqrtm(Q2);
106 
107   for i = 1:m
108     l  = L(:, i);
109     T  = ell_valign(Q1*l, Q2*l);
110     Q  = Q1 - T*Q2;
111     EA = [EA ellipsoid(q, Q'*Q)];
112   end 
113 
114   return;
  1 function [y, Y] = minkdiff(varargin)
  2 %
  3 % MINKDIFF - computes geometric (Minkowski) difference of two ellipsoids in 2D or 3D.
  4 %
  5 %
  6 % Description:
  7 % ------------
  8 %
  9 % MINKDIFF(E1, E2, OPTIONS)  Computes geometric difference of two ellipsoids
 10 %                            E1 - E2, if 1 <= dimension(E1) = dimension(E2) <= 3,
 11 %                            and plots it if no output arguments are specified.
 12 %
 13 %    [y, Y] = MINKDIFF(E1, E2)  Computes geometric difference of two ellipsoids
 14 %                               E1 - E2.
 15 %                               Here y is the center, and Y - array of
 16 %                               boundary points.
 17 %             MINKDIFF(E1, E2)  Plots geometric difference of two ellipsoids
 18 %                               E1 - E2 in default (red) color.
 19 %    MINKDIFF(E1, E2, Options)  Plots geometric difference E1 - E2 using options
 20 %                               given in the Options structure.
 21 %
 22 %    In order for the geometric difference to be nonempty set, ellipsoid E1
 23 %    must be bigger than E2 in the sense that if E1 and E2 had the same center,
 24 %    E2 would be contained inside E1.
 25 %
 26 %
 27 % Options.show_all     - if 1, displays also ellipsoids E1 and E2.
 28 % Options.newfigure    - if 1, each plot command will open a new figure window.
 29 % Options.fill         - if 1, the resulting set in 2D will be filled with color.
 30 % Options.color        - sets default colors in the form [x y z].
 31 % Options.shade = 0-1  - level of transparency (0 - transparent, 1 - opaque).
 32 %
 33 %
 34 % Output:
 35 % -------
 36 %
 37 %    y - center of the resulting set.
 38 %    Y - set of boundary points (vertices) of resulting set.
 39 %
 40 %
 41 % See also:
 42 % ---------
 43 %
 44 %    ELLIPSOID/ELLIPSOID, MINKDIFF_EA, MINKDIFF_IA, ISBIGGER, ISBADDIRECTION,
 45 %                         MINKSUM, MINKSUM_EA, MINKSUM_IA.
 46 %
 47 
 48 % 
 49 % Author:
 50 % -------
 51 %
 52 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 53 %
 54 
 55   global ellOptions;
 56 
 57   if ~isstruct(ellOptions)
 58     evalin('base', 'ellipsoids_init;');
 59   end
 60 
 61   if nargin < 2
 62     error('MINKDIFF: first and second arguments must be single ellipsoids.');
 63   end
 64 
 65   E1 = varargin{1};
 66   E2 = varargin{2};
 67 
 68   if ~(isa(E1, 'ellipsoid')) | ~(isa(E2, 'ellipsoid'))
 69     error('MINKDIFF: first and second arguments must be single ellipsoids.');
 70   end
 71 
 72   if isbigger(E1, E2) == 0
 73     switch nargout
 74       case 0,
 75         fprintf('Geometric difference of these two ellipsoids is empty set.');
 76         return;
 77 
 78       case 1,
 79         y = [];
 80         return;
 81 
 82       otherwise,
 83         y = [];
 84         Y = [];
 85         return;
 86 
 87     end
 88   end
 89 
 90   if nargin > 2
 91     if isstruct(varargin{3})
 92       Options = varargin{3};
 93     else
 94       Options = [];
 95     end
 96   else
 97     Options = [];
 98   end
 99 
100   if ~isfield(Options, 'newfigure')
101     Options.newfigure = 0;
102   end
103 
104   if ~isfield(Options, 'fill')
105     Options.fill = 0;
106   end
107 
108   if ~isfield(Options, 'show_all')
109     Options.show_all = 0;
110   end
111 
112   if ~isfield(Options, 'color')
113     Options.color = [1 0 0];
114   end
115 
116   if ~isfield(Options, 'shade')
117     Options.shade = 0.4;
118   else
119     Options.shade = Options.shade(1, 1);
120   end
121 
122   clr  = Options.color;
123   m    = dimension(E1);
124   n    = dimension(E2);
125   if m ~= n
126     error('MINKDIFF: ellipsoids must be of the same dimension.');
127   end
128   if n > 3
129     error('MINKDIFF: ellipsoid dimension must be not higher than 3.');
130   end
131 
132   %opts.fill          = Options.fill;
133   %opts.shade(1, 1:2) = Options.shade;
134 
135   if nargout == 0
136     ih = ishold;
137   end
138 
139   if (Options.show_all ~= 0) & (nargout == 0)
140     plot([E1 E2], 'b');
141     hold on;
142     if Options.newfigure ~= 0
143       figure;
144     else
145       newplot;
146     end
147   end
148 
149   if ellOptions.verbose > 0
150     if nargout == 0
151       fprintf('Computing and plotting geometric difference of two ellipsoids...\n');
152     else
153       fprintf('Computing geometric difference of two ellipsoids...\n');
154     end
155   end
156 	
157   Q1 = E1.shape;
158   if rank(Q1) < size(Q1, 1)
159     Q1 = regularize(Q1);
160   end
161   Q2 = E2.shape;
162   if rank(Q2) < size(Q2, 1)
163     Q2 = regularize(Q2);
164   end
165   switch n
166     case 2,
167       y      = E1.center - E2.center;
168       phi    = linspace(0, 2*pi, ellOptions.plot2d_grid);
169       l = rm_bad_directions(Q1, Q2, [cos(phi); sin(phi)]);
170       if size(l, 2) > 0
171         [r, Y] = rho(E1, l);
172         [r, X] = rho(E2, l);
173         Y      = Y - X;
174         Y      = [Y Y(:, 1)];
175       else
176         Y = y;
177       end
178       if nargout == 0
179         if Options.fill ~= 0
180           fill(Y(1, :), Y(2, :), clr);
181           hold on;
182         end
183         h = ell_plot(Y);
184         hold on;
185         set(h, 'Color', clr, 'LineWidth', 2);
186         h = ell_plot(y, '.');
187         set(h, 'Color', clr);
188       end
189 
190     case 3,
191       y   = E1.center - E2.center;
192       M   = ellOptions.plot3d_grid/2;
193       N   = M/2;
194       psy = linspace(0, pi, N);
195       phi = linspace(0, 2*pi, M);
196       l   = [];
197       for i = 2:(N - 1)
198         arr = cos(psy(i))*ones(1, M);
199         l   = [l [cos(phi)*sin(psy(i)); sin(phi)*sin(psy(i)); arr]];
200       end
201       l = rm_bad_directions(Q1, Q2, l);
202       if size(l, 2) > 0
203         [r, Y] = rho(E1, l);
204         [r, X] = rho(E2, l);
205         Y      = Y - X;
206       else
207         Y = y;
208       end
209       if nargout == 0
210         vs   = size(Y, 2);
211         if vs > 1
212           chll = convhulln(Y');
213           patch('Vertices', Y', 'Faces', chll, ...
214                 'FaceVertexCData', clr(ones(1, vs), :), 'FaceColor', 'flat', ...
215                 'FaceAlpha', Options.shade(1, 1));
216         else
217           h = ell_plot(y, '*');
218           set(h, 'Color', clr);
219         end
220         hold on;
221         shading interp;
222         lighting phong;
223         material('metal');
224         view(3);
225         %camlight('headlight','local');
226         %camlight('headlight','local');
227         %camlight('right','local');
228         %camlight('left','local');
229       end
230 
231     otherwise,
232       y       = E1.center - E2.center;
233       Y(1, 1) = E1.center - E2.center + sqrt(E2.shape) - sqrt(E1.shape);
234       Y(1, 2) = E1.center - E2.center + sqrt(E1.shape) - sqrt(E2.shape);
235       if nargout == 0
236         h = ell_plot(Y);
237         hold on;
238         set(h, 'Color', clr, 'LineWidth', 2);
239         h = ell_plot(y, '*');
240         set(h, 'Color', clr);
241       end
242 
243   end
244 
245   if nargout == 0
246     if ih == 0
247       hold off;
248     end
249   end
250 
251   if nargout == 1
252     y = Y;
253   end
254   if nargout == 0
255     clear y, Y;
256   end
257 
258   return;
 1 function M = mineig(E)
 2 %
 3 % MINEIG - return the minimal eigenvalue of the ellipsoid.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    M = MINEIG(E)  Returns the smallest eigenvalues of ellipsoids in the array E.
10 %
11 %
12 % Output:
13 % -------
14 %
15 %    M - array of minimal eigenvalues of ellipsoids in the input array E.
16 %
17 %
18 % See also:
19 % ---------
20 %
21 %    ELLIPSOID/ELLIPSOID, ISDEGENERATE, MAXEIG.
22 %
23 
24 %
25 % Author:
26 % -------
27 %
28 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
29 %
30 
31   global ellOptions;
32 
33   if ~isstruct(ellOptions)
34     evalin('base', 'ellipsoids_init;');
35   end
36 
37   if ~(isa(E, 'ellipsoid'))
38     error('MINEIG: input argument must be ellipsoid.')
39   end
40 
41   [m, n] = size(E);
42   M      = [];
43   for i = 1:m
44     mx = [];
45     for j = 1:n
46       if isdegenerate(E(i, j))
47         mx = [mx 0];
48       else
49         mx = [mx min(eig(E(i, j).shape))];
50       end
51     end
52     M = [M; mx];
53   end
54 
55   return;
 1 function M = maxeig(E)
 2 %
 3 % MAXEIG - return the maximal eigenvalue of the ellipsoid.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    M = MAXEIG(E)  Returns the largest eigenvalues of ellipsoids in the array E.
10 %
11 %
12 % Output:
13 % -------
14 %
15 %    M - array of maximal eigenvalues of ellipsoids in the input array E.
16 %
17 %
18 % See also:
19 % ---------
20 %
21 %    ELLIPSOID/ELLIPSOID, ISDEGENERATE, MINEIG.
22 %
23 
24 %
25 % Author:
26 % -------
27 %
28 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
29 %
30 
31   global ellOptions;
32 
33   if ~isstruct(ellOptions)
34     evalin('base', 'ellipsoids_init;');
35   end
36 
37   if ~(isa(E, 'ellipsoid'))
38     error('MAXEIG: input argument must be ellipsoid.')
39   end
40 
41   [m, n] = size(E);
42   M      = [];
43   for i = 1:m
44     mx = [];
45     for j = 1:n
46       mx = [mx max(eig(E(i, j).shape))];
47     end
48     M = [M; mx];
49   end
50 
51   return;
 1 function res = lt(E1, E2)
 2 %
 3 %
 4 % Description:
 5 % ------------
 6 %
 7 %    The opposite of GT.
 8 %
 9 
10 %
11 % Author:
12 % -------
13 %
14 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
15 %
16 
17   res = gt(E2, E1);
18 
19   return;
 1 function res = le(E1, E2)
 2 %
 3 %
 4 % Description:
 5 % ------------
 6 %
 7 %    Same as LT.
 8 %
 9 
10 %
11 % Author:
12 % -------
13 %
14 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
15 %
16 
17   res = lt(E1, E2);
18 
19   return;
  1 function res = isinternal(E, X, s)
  2 %
  3 % ISINTERNAL - checks if given points belong to the union or intersection
  4 %              of ellipsoids in the given array.
  5 %
  6 %
  7 % Description:
  8 % ------------
  9 %
 10 %    RES = ISINTERNAL(E, X, s)  Checks if vectors specified as columns of matrix X
 11 %                               belong to the union (s = 'u'), or
 12 %                               intersection (s = 'i') of the ellipsoids in E.
 13 %                               If E is a single ellipsoid, then this function
 14 %                               checks if points in X belong to E or not.
 15 %                               Ellipsoids in E must be of the same dimension.
 16 %                               Column size of matrix X should match the dimension
 17 %                               of ellipsoids.
 18 %
 19 %    Let E(q, Q) be an ellipsoid with center q and shape matrix Q.
 20 %    Checking if given vector x belongs to E(q, Q) is equivalent to checking
 21 %    if inequality
 22 %                    <(x - q), Q^(-1)(x - q)> <= 1
 23 %    holds.
 24 %    If x belongs to at least one of the ellipsoids in the array, then it belongs
 25 %    to the union of these ellipsoids. If x belongs to all ellipsoids in the array,
 26 %    then it belongs to the intersection of these ellipsoids.
 27 %    The default value of the specifier s = 'u'.
 28 %
 29 %    WARNING: be careful with degenerate ellipsoids.
 30 %
 31 %
 32 % Output:
 33 % -------
 34 %
 35 %    Array of 1 and/or 0.
 36 %    1 - if vector belongs to the union or intersection of ellipsoids,
 37 %    0 - otherwise.
 38 %
 39 %
 40 % See also:
 41 % ---------
 42 %
 43 %    ELLIPSOID/ELLIPSOID, DIMENSION, ISDEGENERATE, DISTANCE.
 44 %
 45 
 46 %
 47 % Author:
 48 % -------
 49 %
 50 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 51 %
 52 
 53   global ellOptions;
 54 
 55   if ~isstruct(ellOptions)
 56     evalin('base', 'ellipsoids_init;');
 57   end
 58 
 59   if ~isa(E, 'ellipsoid')
 60     error('ISINTERNAL: first argument must be an ellipsoid, or an array of ellipsoids.');
 61   end
 62 
 63   dims = dimension(E);
 64   m    = min(min(dims));
 65   n    = max(max(dims));
 66   if m ~= n
 67     error('ISINTERNAL: ellipsoids must be of the same dimension.');
 68   end
 69 
 70   if ~(isa(X, 'double'))
 71     error('ISINTERNAL: second argument must be an array of vectors.');
 72   end
 73 
 74   if (nargin < 3) | ~(ischar(s))
 75     s = 'u';
 76   end
 77 
 78   res = [];
 79 
 80   if (s ~= 'u') & (s ~= 'i')
 81     error('ISINTERNAL: third argument is expected to be either ''u'', or ''i''.');
 82   end
 83 
 84   [k, l] = size(X);
 85   if k ~= n
 86     error('ISINTERNAL: dimensions of ellipsoid and vector do not match.');
 87   end
 88 
 89   for i = 1:l
 90     res = [res isinternal_sub(E, X(:, i), s, k)];
 91   end
 92 
 93   return;
 94 
 95 
 96 
 97 %%%%%%%%
 98   
 99 function res = isinternal_sub(E, x, s, k)
100 %
101 % ISINTERNAL_SUB - compute result for single vector.
102 %
103 
104   global ellOptions;
105 
106   if s == 'u'
107     res = 0;
108   else
109     res = 1;
110   end
111 
112   [m, n] = size(E);
113   for i = 1:m
114     for j = 1:n
115       q = x - E(i, j).center;
116       Q = E(i, j).shape;
117 
118       if rank(Q) < k
119         if ellOptions.verbose > 0
120           fprintf('ISINTERNAL: Warning! There is degenerate ellipsoid in the array.\n');
121           fprintf('            Regularizing...\n');
122         end
123         Q = regularize(Q);
124       end
125  
126       r = q' * ell_inv(Q) * q;
127       if (s == 'u')
128         if (r < 1) | (abs(r - 1) < ellOptions.abs_tol)
129           res = 1;
130           return;
131         end
132       else
133         if (r > 1) & (abs(r - 1) > ellOptions.abs_tol)
134           res = 0;
135           return;
136         end
137       end
138     end
139   end
140   
141   return;
  1 function [res, status] = isinside(E1, E2, s)
  2 %
  3 % ISINSIDE - checks if the intersection of ellipsoids contains the union or
  4 %            intersection of given ellipsoids or polytopes.
  5 %
  6 %
  7 % Description:
  8 % ------------
  9 %
 10 % RES = ISINSIDE(E1, E2, s) Checks if the union (s = 'u') or intersection (s = 'i')
 11 %                           of ellipsoids in E2 lies inside the intersection of
 12 %                           ellipsoids in E1.
 13 %                           Ellipsoids in E1 and E2 must be of the same dimension.
 14 %                           s = 'u' (default) - union of ellipsoids in E2.
 15 %                           s = 'i' - intersection.
 16 %   RES = ISINSIDE(E, P, s) Checks if the union (s = 'u') or intersection (s = 'i')
 17 %                           of polytopes in P lies inside the intersection of
 18 %                           ellipsoids in E.
 19 %                           Ellipsoids in E and polytopes in P must be
 20 %                           of the same dimension.
 21 %                           s = 'u' (default) - union of polytopes in P.
 22 %                           s = 'i' - intersection.
 23 %
 24 %    To check if the union of ellipsoids E2 belongs to the intersection of
 25 %    ellipsoids E1, it is enough to check that every ellipsoid of E2 is 
 26 %    contained in every ellipsoid of E1.
 27 %    Checking if the intersection of ellipsoids in E2 is inside intersection E1
 28 %    can be formulated as quadratically constrained quadratic programming (QCQP)
 29 %    problem.
 30 %    Let E(q, Q) be an ellipsoid with center q and shape matrix Q.
 31 %    To check if this ellipsoid contains the intersection of ellipsoids
 32 %    E(q1, Q1), E(q2, Q2), ..., E(qn, Qn), we define the QCQP problem:
 33 %                      J(x) = <(x - q), Q^(-1)(x - q)> --> max
 34 %    with constraints:
 35 %                       <(x - q1), Q1^(-1)(x - q1)> <= 1   (1)
 36 %                       <(x - q2), Q2^(-1)(x - q2)> <= 1   (2)
 37 %                       ................................
 38 %                       <(x - qn), Qn^(-1)(x - qn)> <= 1   (n)
 39 %
 40 %    If this problem is feasible, i.e. inequalities (1)-(n) do not contradict,
 41 %    or, in other words, intersection of ellipsoids E(q1, Q1), E(q2, Q2), ..., E(qn, Qn)
 42 %    is nonempty, then we can find vector y such that it satisfies inequalities (1)-(n)
 43 %    and maximizes function J. If J(y) <= 1, then ellipsoid E(q, Q) contains
 44 %    the given intersection, otherwise, it does not.
 45 %
 46 %    The intersection of polytopes is a polytope, which is computed
 47 %    by the standard routine of MPT. If the vertices of this polytope
 48 %    belong to the intersection of ellipsoids, then the polytope itself
 49 %    belongs to this intersection.
 50 %    Checking if the union of polytopes belongs to the intersection
 51 %    of ellipsoids is the same as checking if its convex hull belongs 
 52 %    to this intersection.
 53 %
 54 %    We use YALMIP as interface to the optimization tools.
 55 %    (http://control.ee.ethz.ch/~joloef/yalmip.php)
 56 %
 57 %
 58 % Output:
 59 % -------
 60 %
 61 %    RES - result:
 62 %           -1 - problem is infeasible,
 63 %                for example, if s = 'i', but the intersection of ellipsoids in E2
 64 %                is an empty set;
 65 %            0 - intersection is empty;
 66 %            1 - if intersection is nonempty.
 67 %    S   - (optional) status variable returned by YALMIP.
 68 %
 69 %
 70 % See also:
 71 % ---------
 72 %
 73 %    ELLIPSOID/ELLIPSOID, INTERSECT, ISINTERNAL,
 74 %    POLYTOPE/POLYTOPE.
 75 %
 76 
 77 %
 78 % Author:
 79 % -------
 80 %
 81 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 82 %
 83 
 84   global ellOptions;
 85 
 86   if ~isstruct(ellOptions)
 87     evalin('base', 'ellipsoids_init;');
 88   end
 89 
 90   if ~(isa(E1, 'ellipsoid'))
 91     error('ISINSIDE: first input argument must be ellipsoid.');
 92   end
 93 
 94   if ~(isa(E2, 'ellipsoid')) & ~(isa(E2, 'polytope'))
 95     error('ISINSIDE: second input arguments must be ellipsoids or polytope.');
 96   end
 97 
 98   if (nargin < 3) | ~(ischar(s))
 99     s = 'u';
100   end
101 
102   status = [];
103 
104   if isa(E2, 'polytope')
105     [m, n] = size(E2);
106     if s == 'i'
107       PP = E2(1);
108       for j = 1:n
109         PP = PP & E2(j);
110       end
111       X = extreme(PP);
112     else
113       X = [];
114       for j = 1:n
115         X = [X; extreme(E2(j))];
116       end
117     end
118     if isempty(X)
119       res = -1;
120     else
121       res = min(isinternal(E1, X', 'i'));
122     end
123 
124     if nargout < 2
125       clear status;
126     end
127 
128     return;
129   end
130 
131   if s == 'u'
132     [m, n] = size(E1);
133     res    = 1;
134     for i = 1:m
135       for j = 1:n
136         if min(min(contains(E1(i, j), E2))) < 1
137           res = 0;
138           if nargout < 2
139             clear status;
140           end
141           return;
142         end
143       end
144     end
145   elseif min(size(E2) == [1 1]) == 1
146     res = 1;
147     if min(min(contains(E1, E2))) < 1
148       res = 0;
149     end
150   else
151     dims = dimension(E1);
152     m    = min(min(dims));
153     n    = max(max(dims));
154     dims = dimension(E2);
155     k    = min(min(dims));
156     l    = max(max(dims));
157     if (m ~= n) | (k ~= l) | (k ~= m)
158       error('ISINSIDE: ellipsoids must be of the same dimension.');
159     end
160     if ellOptions.verbose > 0
161       fprintf('Invoking YALMIP...\n');
162     end
163     [m, n] = size(E1);
164     res    = 1;
165     for i = 1:m
166       for j = 1:n
167         [res, status] = qcqp(E2, E1(i, j));
168         if res < 1
169           if nargout < 2
170             clear status;
171           end
172           return;
173 	end
174       end
175     end
176   end
177 
178   if nargout < 2
179     clear status;
180   end
181 
182   return;
183 
184 
185 
186 
187 
188 %%%%%%%%
189 
190 function [res, status] = qcqp(EA, E)
191 %
192 % QCQP - formulate quadratically constrained quadratic programming problem
193 %        and invoke external solver.
194 %
195 
196   global ellOptions;
197 
198   [q, Q] = parameters(E(1, 1));
199   if size(Q, 2) > rank(Q)
200     if ellOptions.verbose > 0
201       fprintf('QCQP: Warning! Degenerate ellipsoid.\n');
202       fprintf('      Regularizing...\n');
203     end
204     Q = regularize(Q);
205   end
206   Q = ell_inv(Q);
207   Q = 0.5*(Q + Q');
208 
209   % YALMIP
210   x         = sdpvar(length(Q), 1);
211   objective = x'*Q*x + 2*(-Q*q)'*x + (q'*Q*q - 1);
212   F         = set([]);
213   
214   [m, n] = size(EA);
215   for i = 1:m
216     for j = 1:n
217       [q, Q] = parameters(EA(i, j));
218       if size(Q, 2) > rank(Q)
219         Q = regularize(Q);
220       end
221       Q = ell_inv(Q);
222       Q = 0.5*(Q + Q');
223 
224       % YALMIP
225       F = F + set(x'*Q*x + 2*(-Q*q)'*x + (q'*Q*q - 1) <= 0);
226     end
227   end
228 
229   % YALMIP 
230   status = solvesdp(F, -objective, ellOptions.sdpsettings);
231   
232   if status.problem ~= 0
233     % problem is infeasible, or global minimum cannot be found
234     res = -1;
235     return;
236   end
237 
238   if double(objective) < ellOptions.abs_tol
239     res = 1;
240   else
241     res = 0;
242   end
243 
244   return;
 1 function res = isempty(E)
 2 %
 3 % ISEMPTY - checks if the ellipsoid object is empty.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    RES = ISEMPTY(E)  Given ellipsoidal array E, returns array of ones and zeros
10 %                      specifying which ellipsoids in the array are empty.
11 %
12 %
13 % Output:
14 % -------
15 %
16 %    1 - if ellipsoid is empty, 0 - otherwise.
17 %
18 %
19 % See also:
20 % ---------
21 %
22 %    ELLIPSOID/ELLIPSOID.
23 %
24 
25 %
26 % Author:
27 % -------
28 %
29 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
30 %
31 
32   global ellOptions;
33 
34   if ~isstruct(ellOptions)
35     evalin('base', 'ellipsoids_init;');
36   end
37 
38   if ~(isa(E, 'ellipsoid'))
39     error('ISEMPTY: input argument must be ellipsoid.');
40   end
41 
42   res = ~dimension(E);
43 
44   return;
 1 function res = isdegenerate(E)
 2 %
 3 % ISDEGENERATE - checks if the ellipsoid is degenerate.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %          RES = ISDEGENERATE(E)  Returns 1 if ellipsoid E is degenerate,
10 %                                 0 - otherwise.
11 %
12 %
13 % Output:
14 % -------
15 %
16 %    1 - if ellipsoid E is degenerate, 0 - otherwise.
17 %
18 %
19 % See also:
20 % ---------
21 %
22 %    ELLIPSOID/ELLIPSOID, DIMENSION.
23 %
24 
25 %
26 % Author:
27 % -------
28 %
29 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
30 %
31 
32   [m, n] = size(E);
33 
34   for i = 1:m
35     for j = 1:n
36       if rank(E(i, j).shape) < size(E(i, j).shape, 1)
37         res(i, j) = 1;
38       else
39         res(i, j) = 0;
40       end
41     end
42   end
43 
44   return;
 1 function res = isbigger(E1, E2)
 2 %
 3 % ISBIGGER - checks if one ellipsoid would contain the other if their centers
 4 %            would coincide.
 5 %
 6 %
 7 % Description:
 8 % ------------
 9 %
10 %    RES = ISBIGGER(E1, E2)  Given two single ellipsoids of the same dimension,
11 %                            E1 and E2, check if E1 would contain E2 inside if
12 %                            they were both centered at origin.
13 %
14 %
15 % Output:
16 % -------
17 %
18 %    1 - if ellipsoid E1 would contain E2 inside, 0 - otherwise.
19 %
20 %
21 % See also:
22 % ---------
23 %
24 %    ELLIPSOID/ELLIPSOID.
25 %
26 
27 %
28 % Author:
29 % -------
30 %
31 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
32 %
33 
34   global ellOptions;
35 
36   if ~isstruct(ellOptions)
37     evalin('base', 'ellipsoids_init;');
38   end
39 
40   if ~(isa(E1, 'ellipsoid')) | ~(isa(E2, 'ellipsoid'))
41     error('ISBIGGER: both arguments must be single ellipsoids.');
42   end
43 
44   [k, l] = size(E1);
45   [m, n] = size(E2);
46   if (k > 1) | (l > 1) | (m > 1) | (n > 1)
47     error('ISBIGGER: both arguments must be single ellipsoids.');
48   end
49 
50   [m, r1] = dimension(E1);
51   [n, r2] = dimension(E2);
52   if m ~= n
53     error('ISBIGGER: both ellipsoids must be of the same dimension.');
54   end
55   if r1 < r2
56     res = 0;
57     return;
58   end
59 
60   A = E1.shape;
61   B = E2.shape;
62   if r1 < m
63     if ellOptions.verbose > 0
64       fprintf('ISBIGGER: Warning! First ellipsoid is degenerate.');
65       fprintf('          Regularizing...');
66     end
67     A = regularize(A);
68   end
69 
70   T = ell_simdiag(A, B);
71   if max(abs(diag(T*B*T'))) < (1 + ellOptions.abs_tol)
72     res = 1;
73   else
74     res = 0;
75   end
76 
77   return;
 1 function res = isbaddirection(E1, E2, L)
 2 %
 3 % ISBADDIRECTION - checks if ellipsoidal approximations of geometric difference
 4 %                  of two ellipsoids can be computed for given directions.
 5 %
 6 %
 7 % Description:
 8 % ------------
 9 %
10 %    RES = ISBADDIRECTION(E1, E2, L)  Checks if it is possible to build ellipsoidal
11 %                                     approximation of the geometric difference
12 %                                     of two ellipsoids E1 - E2 in directions
13 %                                     specified by matrix L (columns of L are
14 %                                     direction vectors).
15 %
16 %    Type 'help minkdiff_ea' or 'help minkdiff_ia' for more information.
17 %
18 %
19 % Output:
20 % -------
21 %
22 %    RES - array of 1 or 0 with length being equal to the number of columns
23 %          in matrix L.
24 %          1 marks direction vector as bad - ellipsoidal approximation cannot
25 %          be computed for this direction.
26 %          0 means the opposite.
27 %
28 %
29 % See also:
30 % ---------
31 %
32 %    ELLIPSOID/ELLIPSOID, MINKDIFF, MINKDIFF_EA, MINKDIFF_IA.
33 %
34 
35 %
36 % Author:
37 % -------
38 %
39 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
40 %
41 
42   global ellOptions;
43 
44   if ~isstruct(ellOptions)
45     evalin('base', 'ellipsoids_init;');
46   end
47 
48   [k, d] = size(L);
49 
50   if isbigger(E1, E2) == 0
51     if ellOptions.verbose > 0
52       fprintf('ISBADDIRECTION: geometric difference of these two ellipsoids is empty set.\n');
53       fprintf('                All directions are bad.\n');
54     end
55     res = ones(1, d);
56     return;
57   end
58 
59   n = dimension(E1);
60 
61   if k ~= n
62     error('ISBADDIRECTION: direction vectors must be of the same dimension as ellipsoids.');
63   end
64 
65   res = [];
66   Q1  = E1.shape;
67   Q2  = E2.shape;
68   T   = ell_simdiag(Q2, Q1);
69   a   = min(abs(diag(T*(Q1)*T')));
70   for i = 1:d
71     l = L(:, i);
72     p = sqrt(l'*Q1*l)/sqrt(l'*Q2*l);
73     if a > p
74       res = [res 0];
75     else
76       res = [res 1];
77     end
78   end
79 
80   return;
 1 function I = inv(E)
 2 %
 3 % INV - inverts shape matrices of ellipsoids in the given array.
 4 %
 5 %
 6 % Description:
 7 % ------------
 8 %
 9 %    I = INV(E)  Inverts shape matrices of ellipsoids in the array E.
10 %                In case shape matrix is sigular, it is regularized before inversion.
11 %
12 %
13 % Output:
14 % -------
15 %
16 %    I - array of ellipsoids with inverted shape matrices.
17 %
18 %
19 % See also:
20 % ---------
21 %
22 %    ELLIPSOID/ELLIPSOID.
23 %
24 
25 %
26 % Author:
27 % -------
28 %
29 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
30 %
31 
32   if ~(isa(E, 'ellipsoid'))
33     error('INV: input argument must be array of ellipsoids.');
34   end
35 
36   I      = E;
37   [m, n] = size(I);
38 
39   for i = 1:m
40     for j = 1:n
41       if isdegenerate(I(i, j))
42         Q = regularize(I(i, j).shape);
43       else
44         Q = I(i, j).shape;
45       end
46       Q             = ell_inv(Q);
47       I(i, j).shape = 0.5*(Q + Q');
48     end
49   end
50 
51   return;
  1 function E = intersection_ia(E1, X)
  2 %
  3 % INTERSECTION_IA - internal ellipsoidal approximation of the intersection of
  4 %                   of ellipsoid and ellipsoid, or ellipsoid and halfspace,
  5 %                   or ellipsoid and polytope.
  6 %
  7 %
  8 % Description:
  9 % ------------
 10 %
 11 %  E = INTERSECTION_IA(E1, E2) Given two ellipsoidal arrays of equal sizes,
 12 %                              E1 and E2, or, alternatively, E1 or E2 must be
 13 %                              a single ellipsoid, comuptes the internal
 14 %                              ellipsoidal approximations of intersections of
 15 %                              two corresponding ellipsoids from E1 and from E2.
 16 %   E = INTERSECTION_IA(E1, H) Given array of ellipsoids E1 and array of
 17 %                              hyperplanes H whose sizes match, computes
 18 %                              the internal ellipsoidal approximations of
 19 %                              intersections of ellipsoids and halfspaces
 20 %                              defined by hyperplanes in H.
 21 %                              If v is normal vector of hyperplane and c - shift,
 22 %                              then this hyperplane defines halfspace
 23 %                                         <v, x> <= c.
 24 %   E = INTERSECTION_IA(E1, P) Given array of ellipsoids E1 and array of
 25 %                              polytopes P whose sizes match, computes
 26 %                              the internal ellipsoidal approximations of
 27 %                              intersections of ellipsoids E1 and polytopes P.
 28 %
 29 %
 30 % Output:
 31 % -------
 32 %
 33 %    E - array of internal approximating ellipsoids;
 34 %        entries can be empty ellipsoids if the corresponding intersection
 35 %        is empty.
 36 %
 37 %
 38 % See also:
 39 % ---------
 40 %
 41 %    ELLIPSOID/ELLIPSOID, INTERSECTION_EA, INTERSECT, DISTANCE,
 42 %    HYPERPLANE/HYPERPLANE, POLYTOPE/POLYTOPE.
 43 %
 44 
 45 %
 46 % Author:
 47 % -------
 48 %
 49 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 50 %
 51 
 52   global ellOptions;
 53 
 54   if ~isstruct(ellOptions)
 55     evalin('base', 'ellipsoids_init;');
 56   end
 57 
 58   if ~(isa(E1, 'ellipsoid'))
 59     error('INTERSECTION_IA: first input argument must be ellipsoid.');
 60   end
 61   if ~(isa(X, 'ellipsoid')) & ~(isa(X, 'hyperplane')) & ~(isa(X, 'polytope'))
 62     error('INTERSECTION_IA: second input argument must be ellipsoid, hyperplane or polytope.');
 63   end
 64 
 65   [k, l] = size(E1);
 66   [m, n] = size(X);
 67   dims1  = dimension(E1);
 68 
 69   if isa(X, 'polytope')
 70     dims2 = [];
 71     for i = 1:m
 72       dd = [];
 73       for j = 1:n
 74         dd = [dd dimension(X(j))];
 75       end
 76       dims2 = [dims2; dd];
 77     end
 78   else
 79     dims2 = dimension(X);
 80   end
 81 
 82   mn1   = min(min(dims1));
 83   mn2   = min(min(dims2));
 84   mx1   = max(max(dims1));
 85   mx2   = max(max(dims2));
 86   if (mn1 ~= mx1) | (mn2 ~= mx2) | (mx1 ~= mx2)
 87     if isa(X, 'hyperplane')
 88       error('INTERSECTION_IA: ellipsoids and hyperplanes must be of the same dimension.');
 89     elseif isa(X, 'polytope')
 90       error('INTERSECTION_IA: ellipsoids and polytopes must be of the same dimension.');
 91     else
 92       error('INTERSECTION_IA: ellipsoids must be of the same dimension.');
 93     end
 94   end
 95 
 96   t1     = k * l;
 97   t2     = m * n;
 98   if (t1 > 1) & (t2 > 1) & ((k ~= m) | (l ~= n))
 99     if isa(X, 'hyperplane')
100       error('INTERSECTION_IA: sizes of ellipsoidal and hyperplane arrays do not match.');
101     elseif isa(X, 'polytope')
102       error('INTERSECTION_IA: sizes of ellipsoidal and polytope arrays do not match.');
103     else
104       error('INTERSECTION_IA: sizes of ellipsoidal arrays do not match.');
105     end
106   end
107 
108   E = [];
109   if (t1 > 1) & (t2 > 1)
110     for i = 1:k
111       e = [];
112       for j = 1:l
113         if isa(X, 'polytope')
114           e = [e l_polyintersect(E1(i, j), X(j))];
115         else
116           e = [e l_intersection_ia(E1(i, j), X(i, j))];
117         end
118       end
119       E = [E; e];
120     end
121   elseif t1 > 0
122     for i = 1:k
123       e = [];
124       for j = 1:l
125         if isa(X, 'polytope')
126           e = [e l_polyintersect(E1(i, j), X)];
127         else
128           e = [e l_intersection_ia(E1(i, j), X)];
129         end
130       end
131       E = [E; e];
132     end
133   else
134     for i = 1:m
135       e = [];
136       for j = 1:n
137         if isa(X, 'polytope')
138           e = [e l_polyintersect(E1(i, j), X(j))];
139         else
140           e = [e l_intersection_ia(E1, X(i, j))];
141         end
142       end
143       E = [E; e];
144     end
145   end
146 
147   return;
148 
149 
150 
151 
152 
153 %%%%%%%%
154 
155 function E = l_intersection_ia(E1, E2)
156 %
157 % L_INTERSECTION_IA - computes internal ellipsoidal approximation of intersection
158 %                     of single ellipsoid with single ellipsoid or halfspace.
159 %
160 
161   global ellOptions;
162 
163   if isa(E2, 'ellipsoid')
164     if E1 == E2
165       E = E1;
166     elseif ~intersect(E1, E2)
167       E = ellipsoid;
168     else
169       E = ellintersection_ia([E1 E2]);
170     end
171     return;
172   end
173 
174   q1 = E1.center;
175   Q1 = E1.shape;
176   if rank(Q1) < size(Q1, 1)
177     Q1 = ell_inv(regularize(Q1));
178   else
179     Q1 = ell_inv(Q1);
180   end
181 
182   [v, c] = parameters(-E2);
183   c      = c/sqrt(v'*v);
184   v      = v/sqrt(v'*v);
185   if (v'*q1 > c) & ~(intersect(E1, E2))
186     E = E1;
187     return;
188   end
189   if (v'*q1 < c) & ~(intersect(E1, E2))
190     E = ellipsoid;
191     return;
192   end
193 
194   [q, Q] = parameters(hpintersection(E1, E2));
195   [r, x] = rho(E1, v);
196   h      = 2*sqrt(maxeig(E1));
197   q2     = q + h*v;
198   Q2     = (v*v')/(h^2);
199   st1    = 0;
200   st2    = 0;
201   e1     = (q1 - q)'*Q1*(q1 - q);
202   e2     = (q2 - x)'*Q2*(q2 - x);
203   a1     = (1 - e2)/(1 - e1*e2);
204   a2     = (1 - e1)/(1 - e1*e2);
205   Q      = a1*Q1 + a2*Q2;
206   Q      = 0.5*(Q + Q');
207   q      = ell_inv(Q)*(a1*Q1*q1 + a2*Q2*q2);
208   Q      = Q/(1 - (a1*q1'*Q1*q1 + a2*q2'*Q2*q2 - q'*Q*q));
209   Q      = ell_inv(Q);
210   Q      = (1-ellOptions.abs_tol)*0.5*(Q + Q');
211   E      = ellipsoid(q, Q);
212   
213   return;
214 
215   
216 
217 
218 
219 %%%%%%%%
220 
221 function EA = l_polyintersect(E, P)
222 %
223 % L_POLYINTERSECT - computes internal ellipsoidal approximation of intersection
224 %                   of single ellipsoid with single polytope.
225 %
226 
227   global ellOptions;
228 
229   EA = E;
230   HA = polytope2hyperplane(P);
231   n  = size(HA, 2);
232 
233   for i = 1:n
234     EA = intersection_ia(EA, HA(i));
235   end
236 
237   if isinside(E, P)
238     EA = getInnerEllipsoid(P);
239     return;
240   end
241 
242   return;
  1 function E = intersection_ea(E1, X)
  2 %
  3 % INTERSECTION_EA - external ellipsoidal approximation of the intersection of
  4 %                   two ellipsoids, or ellipsoid and halfspace, or ellipsoid
  5 %                   and polytope.
  6 %
  7 %
  8 % Description:
  9 % ------------
 10 %
 11 %     E = INTERSECTION_EA(E1, E2) Given two ellipsoidal arrays of equal sizes,
 12 %                                 E1 and E2, or, alternatively, E1 or E2 must be
 13 %                                 a single ellipsoid, computes the ellipsoid 
 14 %                                 that contains the intersection of two
 15 %                                 corresponding ellipsoids from E1 and from E2.
 16 %      E = INTERSECTION_EA(E1, H) Given array of ellipsoids E1 and array of
 17 %                                 hyperplanes H whose sizes match, computes
 18 %                                 the external ellipsoidal
 19 %                                 approximations of intersections of ellipsoids
 20 %                                 and halfspaces defined by hyperplanes in H.
 21 %                                 If v is normal vector of hyperplane and c - shift,
 22 %                                 then this hyperplane defines halfspace
 23 %                                         <v, x> <= c.
 24 %      E = INTERSECTION_EA(E1, P) Given array of ellipsoids E1 and array of
 25 %                                 polytopes P whose sizes match, computes
 26 %                                 the external ellipsoidal approximations
 27 %                                 of intersections of ellipsoids E1 and
 28 %                                 polytopes P.
 29 %
 30 %    The method used to compute the minimal volume overapproximating ellipsoid
 31 %    is described in "Ellipsoidal Calculus Based on Propagation and Fusion"
 32 %    by Lluis Ros, Assumpta Sabater and Federico Thomas;
 33 %    IEEE Transactions on Systems, Man and Cybernetics, Vol.32, No.4, pp.430-442, 2002.
 34 %    For more information, visit
 35 %               http://www-iri.upc.es/people/ros/ellipsoids.html
 36 %
 37 %
 38 % Output:
 39 % -------
 40 %
 41 %    E - array of external approximating ellipsoids;
 42 %        entries can be empty ellipsoids if the corresponding intersection
 43 %        is empty.
 44 %
 45 %
 46 % See also:
 47 % ---------
 48 %
 49 %    ELLIPSOID/ELLIPSOID, INTERSECTION_IA, INTERSECT, DISTANCE,
 50 %    HYPERPLANE/HYPERPLANE, POLYTOPE/POLYTOPE.
 51 %
 52 
 53 %
 54 % Author:
 55 % -------
 56 %
 57 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 58 %
 59 
 60   global ellOptions;
 61 
 62   if ~isstruct(ellOptions)
 63     evalin('base', 'ellipsoids_init;');
 64   end
 65 
 66   if ~(isa(E1, 'ellipsoid'))
 67     error('INTERSECTION_EA: first input argument must be ellipsoid.');
 68   end
 69   if ~(isa(X, 'ellipsoid')) & ~(isa(X, 'hyperplane')) & ~(isa(X, 'polytope'))
 70     error('INTERSECTION_EA: second input argument must be ellipsoid, hyperplane or polytope.');
 71   end
 72 
 73   [k, l] = size(E1);
 74   [m, n] = size(X);
 75   dims1  = dimension(E1);
 76 
 77   if isa(X, 'polytope')
 78     dims2 = [];
 79     for i = 1:m
 80       dd = [];
 81       for j = 1:n
 82         dd = [dd dimension(X(j))];
 83       end
 84       dims2 = [dims2; dd];
 85     end
 86   else
 87     dims2 = dimension(X);
 88   end
 89   mn1   = min(min(dims1));
 90   mn2   = min(min(dims2));
 91   mx1   = max(max(dims1));
 92   mx2   = max(max(dims2));
 93 
 94   if (mn1 ~= mx1) | (mn2 ~= mx2) | (mx1 ~= mx2)
 95     if isa(X, 'hyperplane')
 96       error('INTERSECTION_EA: ellipsoids and hyperplanes must be of the same dimension.');
 97     elseif isa(X, 'polytope')
 98       error('INTERSECTION_EA: ellipsoids and polytopes must be of the same dimension.');
 99     else
100       error('INTERSECTION_EA: ellipsoids must be of the same dimension.');
101     end
102   end
103 
104   t1     = k * l;
105   t2     = m * n;
106   if (t1 > 1) & (t2 > 1) & ((k ~= m) | (l ~= n))
107     if isa(X, 'hyperplane')
108       error('INTERSECTION_EA: sizes of ellipsoidal and hyperplane arrays do not match.');
109     elseif isa(X, 'polytope')
110       error('INTERSECTION_EA: sizes of ellipsoidal and polytope arrays do not match.');
111     else
112       error('INTERSECTION_EA: sizes of ellipsoidal arrays do not match.');
113     end
114   end
115 
116   E = [];
117   if (t1 > 1) & (t2 > 1)
118     for i = 1:k
119       e = [];
120       for j = 1:l
121         if isa(X, 'polytope')
122           e = [e l_polyintersect(E1, X(j))];
123         else
124           e = [e l_intersection_ea(E1(i, j), X(i, j))];
125         end
126       end
127       E = [E; e];
128     end
129   elseif t1 > 0
130     for i = 1:k
131       e = [];
132       for j = 1:l
133         if isa(X, 'polytope')
134           e = [e l_polyintersect(E1, X)];
135         else
136           e = [e l_intersection_ea(E1(i, j), X)];
137         end
138       end
139       E = [E; e];
140     end
141   else
142     for i = 1:m
143       e = [];
144       for j = 1:n
145         if isa(X, 'polytope')
146           e = [e l_polyintersect(E1, X(j))];
147         else
148           e = [e l_intersection_ea(E1, X(i, j))];
149         end
150       end
151       E = [E; e];
152     end
153   end
154 
155   return;
156 
157 
158 
159 
160 
161 %%%%%%%%
162 
163 function E = l_intersection_ea(E1, E2)
164 %
165 % L_INTERSECTION_EA - computes external ellipsoidal approximation of intersection
166 %                     of single ellipsoid with single ellipsoid or halfspace.
167 %
168 
169   global ellOptions;
170 
171   q1 = E1.center;
172   Q1 = E1.shape;
173   if rank(Q1) < size(Q1, 1)
174     Q1 = ell_inv(regularize(Q1));
175   else
176     Q1 = ell_inv(Q1);
177   end
178 
179   if isa(E2, 'hyperplane')
180     [v, c] = parameters(-E2);
181     c      = c/sqrt(v'*v);
182     v      = v/sqrt(v'*v);
183     if (v'*q1 > c) & ~(intersect(E1, E2))
184       E = E1;
185       return;
186     end
187     if (v'*q1 < c) & ~(intersect(E1, E2))
188       E = ellipsoid;
189       return;
190     end
191     h  = 2*sqrt(maxeig(E1));
192     q2 = c*v + h*v;
193     Q2 = (v*v')/(h^2);
194 
195     [q, Q] = parameters(hpintersection(E1, E2));
196     q2     = q + h*v;
197   else
198     if E1 == E2
199       E = E1;
200       return;
201     end
202     if ~intersect(E1, E2)
203       E = ellipsoid;
204       return;
205     end
206     q2 = E2.center;
207     Q2 = E2.shape;
208     if rank(Q2) < size(Q2, 1)
209       Q2 = ell_inv(regularize(Q2));
210     else
211       Q2 = ell_inv(Q2);
212     end
213   end
214 
215   a = l_get_lambda(q1, Q1, q2, Q2, isa(E2, 'hyperplane'));
216   X = a*Q1 + (1 - a)*Q2;
217   X = 0.5*(X + X');
218  % if rank(X) < size(X, 1)
219  %   X = regularize(X);
220  % end
221   Y = ell_inv(X);
222   Y = 0.5*(Y + Y');
223   k = 1 - a*(1 - a)*(q2 - q1)'*Q2*Y*Q1*(q2 - q1);
224   q = Y*(a*Q1*q1 + (1 - a)*Q2*q2);
225   Q = (1+ellOptions.abs_tol)*k*Y;
226   E = ellipsoid(q, Q); 
227   
228   return;
229 
230 
231 
232 
233 
234 %%%%%%%%
235 
236 function a = l_get_lambda(q1, Q1, q2, Q2, flag)
237 %
238 % L_GET_LAMBDA - find parameter value for minimal volume ellipsoid.
239 %
240 
241   [a, f] = fzero(@ell_fusionlambda, 0.5, [], q1, Q1, q2, Q2, size(q1, 1));
242 
243   if (a < 0) | (a > 1) 
244     if flag | (det(Q1) > det(Q2))
245       a = 1;
246     else
247       a = 0;
248     end
249   end
250 
251   return;
252 
253 
254 
255 
256 
257 %%%%%%%%
258 
259 function EA = l_polyintersect(E, P)
260 %
261 % L_POLYINTERSECT - computes external ellipsoidal approximation of intersection
262 %                   of single ellipsoid with single polytope.
263 %
264 
265   global ellOptions;
266 
267   EA = E;
268   HA = polytope2hyperplane(P);
269   n  = size(HA, 2);
270 
271   if isinside(E, P)
272     EA = getOutterEllipsoid(P);
273     return;
274   end
275 
276   for i = 1:n
277     EA = intersection_ea(EA, HA(i));
278   end
279 
280   return;
  1 function [res, status] = intersect(E, X, s)
  2 %
  3 % INTERSECT - checks if the union or intersection of ellipsoids intersects
  4 %             given ellipsoid, hyperplane or polytope.
  5 %
  6 %
  7 % Description:
  8 % ------------
  9 %
 10 % RES = INTERSECT(E, X, s)  Checks if the union (s = 'u') or intersection (s = 'i')
 11 %                           of ellipsoids in E intersects with objects in X.
 12 %                           X can be array of ellipsoids, array of hyperplanes,
 13 %                           or array of polytopes.
 14 %                           Ellipsoids, hyperplanes or polytopes in X must have
 15 %                           the same dimension as ellipsoids in E.
 16 %                           s = 'u' (default) - union of ellipsoids in E.
 17 %                           s = 'i' - intersection.
 18 %
 19 %    If we need to check the intersection of union of ellipsoids in E (s = 'u'),
 20 %    or if E is a single ellipsoid, it can be done by calling distance function
 21 %    for each of the ellipsoids in E and X, and if it returns negative value,
 22 %    the intersection is nonempty.
 23 %    Checking if the intersection of ellipsoids in E (with size of E greater than 1)
 24 %    intersects with ellipsoids or hyperplanes in X is more difficult.
 25 %    This problem can be formulated as quadratically constrained quadratic
 26 %    programming (QCQP) problem.
 27 %    Let E(q, Q) be an ellipsoid with center q and shape matrix Q.
 28 %    To check if this ellipsoid intersects (or touches) the intersection
 29 %    of ellipsoids E(q1, Q1), E(q2, Q2), ..., E(qn, Qn), we define the QCQP
 30 %    problem:
 31 %                      J(x) = <(x - q), Q^(-1)(x - q)> --> min
 32 %    with constraints:
 33 %                       <(x - q1), Q1^(-1)(x - q1)> <= 1   (1)
 34 %                       <(x - q2), Q2^(-1)(x - q2)> <= 1   (2)
 35 %                       ................................
 36 %                       <(x - qn), Qn^(-1)(x - qn)> <= 1   (n)
 37 %
 38 %    If this problem is feasible, i.e. inequalities (1)-(n) do not contradict,
 39 %    or, in other words, intersection of ellipsoids E(q1, Q1), E(q2, Q2), ..., E(qn, Qn)
 40 %    is nonempty, then we can find vector y such that it satisfies inequalities (1)-(n)
 41 %    and minimizes function J. If J(y) <= 1, then ellipsoid E(q, Q) intersects
 42 %    or touches the given intersection, otherwise, it does not.
 43 %    To check if E(q, Q) intersects the union of E(q1, Q1), E(q2, Q2), ..., E(qn, Qn),
 44 %    we compute the distances from this ellipsoids to those in the union. 
 45 %    If at least one such distance is negative, then E(q, Q) does intersect
 46 %    the union.
 47 %
 48 %    If we check the intersection of ellipsoids with hyperplane H(v, c),
 49 %    it is enough to check the feasibility of the problem
 50 %                        1'x --> min
 51 %    with constraints (1)-(n), plus
 52 %                      <v, x> - c = 0.
 53 %
 54 %    Checking the intersection of ellipsoids with polytope P(A, b) reduces
 55 %    to checking the feasibility of the problem 
 56 %                        1'x --> min
 57 %    with constraints (1)-(n), plus
 58 %                         Ax <= b.
 59 %
 60 %    We use YALMIP as interface to optimization tools.
 61 %    (http://control.ee.ethz.ch/~joloef/yalmip.php)
 62 %
 63 %
 64 % Output:
 65 % -------
 66 %
 67 %    RES - result:
 68 %           -1 - problem is infeasible,
 69 %                for example, if s = 'i', but the intersection of ellipsoids in E
 70 %                is an empty set;
 71 %            0 - intersection is empty;
 72 %            1 - if intersection is nonempty.
 73 %
 74 %     S   - (optional) status variable returned by YALMIP.
 75 %
 76 %
 77 % See also:
 78 % ---------
 79 %
 80 %    ELLIPSOID/ELLIPSOID, ISINSIDE, DISTANCE,
 81 %    HYPERPLANE/HYPERPLANE, POLYTOPE/POLYTOPE,
 82 %    YALMIP.
 83 %
 84 
 85 %
 86 % Author:
 87 % -------
 88 %
 89 %    Alex Kurzhanskiy <akurzhan@eecs.berkeley.edu>
 90 %
 91 
 92   global ellOptions;
 93 
 94   if ~isstruct(ellOptions)
 95     evalin('base', 'ellipsoids_init;');
 96   end
 97 
 98   if ~(isa(E, 'ellipsoid'))
 99     error('INTERSECT: first input argument must be ellipsoid.');
100   end
101   if ~(isa(X, 'ellipsoid')) & ~(isa(X, 'hyperplane')) & ~(isa(X, 'polytope'))
102     error('INTERSECT: second input argument must be ellipsoid, hyperplane or polytope.');
103   end
104 
105   if (nargin < 3) | ~(ischar(s))
106     s = 'u';
107   end
108 
109   if s == 'u'
110     [m, n] = size(E);
111     res    = (distance(E(1, 1), X) <= ellOptions.abs_tol);
112     for i = 1:m
113       for j = 1:n
114         if (i > 1) | (j > 1)
115           res = res | (distance(E(i, j), X) <= ellOptions.abs_tol);
116         end
117       end
118     end
119     status = [];
120   elseif min(size(E) == [1 1]) == 1
121     res    = (distance(E, X) <= ellOptions.abs_tol);
122     status = [];
123   elseif isa(X, 'ellipsoid')
124     dims = dimension(E);
125     m    = min(min(dims));
126     n    = max(max(dims));
127     dims = dimension(X);
128     k    = min(min(dims));
129     l    = max(max(dims));
130     if (m ~= n) | (k ~= l) | (k ~= m)
131       error('INTERSECT: ellipsoids must be of the same dimension.');
132     end
133     if ellOptions.verbose > 0
134       fprintf('Invoking YALMIP...\n');
135     end
136     [m, n] = size(X);
137     res    = [];
138     status = [];
139     for i = 1:m
140       r = [];
141       s = [];
142       for j = 1:n
143         [rr, ss] = qcqp(E, X(i, j));
144         r        = [r rr];
145 	s        = [s ss];
146       end
147       res    = [res; r];
148       status = [status; s];
149     end
150   elseif isa(X, 'hyperplane')
151     dims = dimension(E);
152     m    = min(min(dims));
153     n    = max(max(dims));
154     dims = dimension(X);
155     k    = min(min(dims));
156     l    = max(max(dims));
157     if (m ~= n) | (k ~= l) | (k ~= m)
158       error('INTERSECT: ellipsoids and hyperplanes must be of the same dimension.');
159     end
160     if ellOptions.verbose > 0
161       fprintf('Invoking YALMIP...\n');
162     end
163     [m, n] = size(X);
164     res    = [];
165     status = [];
166     for i = 1:m
167       r = [];
168       s = [];
169       for j = 1:n
170         [rr, ss] = lqcqp(E, X(i, j));
171         r        = [r rr];
172         s        = [s ss];
173       end
174       res    = [res; r];
175       status = [status; s];
176     end
177   else
178     [m, n] = size(X);
179     dims   = dimension(E);
180     mm     = min(min(dims));
181     nn     = max(max(dims));
182     dims   = [];
183     for i = 1:m
184       dd = [];
185       for j = 1:n
186         dd = [dd dimension(X(j))];
187       end
188       dims = [dims; dd];
189     end
190     k = min(min(dims));
191     l = max(max(dims));
192     if (mm ~= nn) | (k ~= l) | (k ~= mm)
193       error('INTERSECT: ellipsoids and polytopes must be of the same dimension.');
194     end
195     if ellOptions.verbose > 0
196       fprintf('Invoking YALMIP...\n');
197     end
198     res    = [];
199     status = [];
200     for i = 1:m
201       r = [];
202       s = [];
203       for j = 1:n
204         [rr, ss] = lqcqp2(E, X(j));
205         r        = [r rr];
206         s        = [s ss];
207       end
208       res    = [res; r];
209       status = [status; s];
210     end
211   end
212 
213   if nargout < 2
214     clear status;
215   end
216 
217   return;
218 
219 
220 
221 
222 
223 %%%%%%%%
224 
225 function [res, status] = qcqp(EA, E)
226 %
227 % QCQP - formulate quadratically constrained quadratic programming problem
228 %        and invoke external solver.
229 %
230 
231   global ellOptions;
232 
233   [q, Q] = parameters(E(1, 1));
234   if size(Q, 2) > rank(Q)
235     if ellOptions.verbose > 0
236       fprintf('QCQP: Warning! Degenerate ellipsoid.\n');
237       fprintf('      Regularizing...\n');
238     end
239     Q = regularize(Q);
240   end
241   Q = ell_inv(Q);
242   Q = 0.5*(Q + Q');
243   
244   % YALMIP
245   x         = sdpvar(length(Q), 1);
246   objective = x'*Q*x + 2*(-Q*q)'*x + (q'*Q*q - 1);
247   F         = set([]);
248   
249   [m, n] = size(EA);
250   for i = 1:m
251     for j = 1:n
252       [q, Q] = parameters(EA(i, j));
253       if size(Q, 2) > rank(Q)
254         Q = regularize(Q);
255       end
256       Q = ell_inv(Q);
257       Q = 0.5*(Q + Q');
258 
259       % YALMIP
260       F = F + set(x'*Q*x + 2*(-Q*q)'*x + (q'*Q*q - 1) <= 0);
261     end
262   end
263 
264   % YALMIP 
265   status = solvesdp(F, objective, ellOptions.sdpsettings);
266   
267   if status.problem ~= 0
268     % problem is infeasible, or global minimum cannot be found
269     res = -1;
270     return;
271   end
272 
273   if double(objective) < ellOptions.abs_tol
274     res = 1;
275   else
276     res = 0;
277   end
278 
279   return;
280 
281 
282 
283 
284 
285 %%%%%%%%
286 
287 function [res, status] = lqcqp(EA, H)
288 %
289 % LQCQP - formulate quadratic programming problem with linear and quadratic constraints,
290 %         and invoke external solver.
291 %
292 
293   global ellOptions;
294 
295   [v, c] = parameters(H);
296   if c < 0
297     c = -c;
298     v = -v;
299   end
300 
301   % YALMIP
302   x         = sdpvar(size(v, 1), 1);
303   objective = v'*x - c;
304   F         = set([]);
305   
306   [m, n] = size(EA);
307   for i = 1:m
308     for j = 1:n
309       [q, Q] = parameters(EA(i, j));
310       if size(Q, 2) > rank(Q)
311         Q = regularize(Q);
312       end
313       Q  = ell_inv(Q);
314       Q  = 0.5*(Q + Q');
315 
316       % YALMIP
317       F = F + set(x'*Q*x - 2*q'*Q*x + (q'*Q*q - 1) <=