Copyright (c) [2019] Madeline Hamilton, Pavel Nadolsky (nadolsky@smu.edu)

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

1. The above copyright notice, citation policy, permission notice, and program description 
shall be included in all copies or substantial portions of the Software.

2. You shall cite the following research article on which this program is based in the publications 
that use any results obtained with this program:

Fast Ellipsoidal Fitting of Discrete Multi-Dimensional Data,  
by Rafey Anwar, Madeline Hamilton, Pavel Nadolsky
E-print: arXiv:1901.05511 [physics.data-an] (https://arxiv.org/abs/1901.05511)
SMU Journal of Undergraduate Research: Vol. 5 , Article 4 
https://doi.org/10.25172/jour.5.1.4


Permission notice: THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
 
Program description: Reconstruction of a 5-dimensional ellipsoid from quasi-Gaussian data
Version 2019-01-10

This Mathematica notebook is the fastest available version of the ellipsoid reconstruction algorithm developed by Anwar, Hamilton and Nadolsky. 
Its inputs are a dataset and a desired confidence level, and its outputs are graphs visualizing the ellipses fitted to each orthogonal projection of the data, 
matrices calculated using both the Covariance Matrix Method and the Convex Hull Method that characterize the reconstructed ellipsoids, 
and metrics that evaluate the quality of these ellipsoids.The notebook may be repurposed for any dataset and any confidence level, 
given that the data is entirely numerical takes the form of an array of equally sized arrays.
The notebook was developed for Mathematica 11.1 and is not guaranteed to work for earlier or later versions.

Visual and Convex Hull Setup

This section sets up the style of the notebook's visuals. This may be edited to taste. Matrices necessary for the Convex Hull Method are also initialized in this cell, so this cell must be executed before the fit.

(* Visual setup *) 
Switch[StringTake[$Version,{1,3}]//ToExpression//Round,
6,<<"BarCharts`";<<"Histograms`";<<"PieCharts`";<<"ErrorBarPlots`";SetOptions[{Plot,ListPlot,ErrorListPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
7,SetOptions[{Plot,ListPlot,ErrorListPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
8,SetOptions[{Plot,ListPlot,ListLogLinearPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
9,SetOptions[{Plot,ListPlot,ListLogLinearPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
10,SetOptions[{Plot,ListPlot,ListLogPlot,ListLogLinearPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
11,SetOptions[{Plot,ListPlot,ListLogPlot,ListLogLinearPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
_,<< Graphics`Graphics`;<<Graphics`FilledPlot`;
];
(* Matrices necessary for Convex Hull Method *)
CC=(({{0, 0, 2, 0, 0, 0}, {0, -1, 0, 0, 0, 0}, {2, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}}));
{U,L,V}=SingularValueDecomposition[CC];
U2={U[[1]],U[[2]],U[[3]],U[[6]],U[[5]],U[[4]]};
V2={V[[1]],V[[2]],V[[3]],V[[6]],V[[5]],V[[4]]};
CC-U.L.Transpose[V];
CC-U2.L.Transpose[V2];
U.Transpose[U];
V.Transpose[V];
Transpose[U].U;
Transpose[V].V;
cc1=CC[[1;;3,1;;3]];
u1=U[[1;;3,1;;3]];
v1=V[[1;;3,1;;3]];
l=L[[1;;3,1;;3]];
linv=Inverse[l];

Initialize Inputs

In this section, the algorithm's two main inputs, the dataset and the confidence level, are initialized. This is cell must be adapted if one wishes to repurpose this notebook for their own analyses.
The dataset "data" must be an array of arrays: all observations should be arrays with an equal number of elements, and these arrays should in turn be stored in an array. All data must be numerical. "nsigma" is a number that represents the confidence level outlined by the ellipsoid. Below are the proper "nsigma" values for the most commonly used confidence levels. A Z Table can be used to find the z-score of any other confidence level.

68% confidence: 1.
90% confidence: 1.645
95% confidence: 1.96
99% confidence: 2.575

(* plist initialization *)
file= OpenRead["data.txt"];
headers=Read[file,{String,String}];
data=ReadList[file,LF[Number,Number,Number,Number,Number]]; (* <- in this line, the number of "Number" entries must match the number of your dimensions *)
Close[file];

(* nsigma initialization *)
nsigma= 1;

Data Normalization

In this section, the data is normalized to ensure optimal results. This cell is not essential and may be skipped. However, if the scales of the input data's dimensions differ significantly, executing this cell is recommended. 
The data is normalized such that the mean value across each dimension is exactly 1.

MeanX=Table[data/.LF[a__]:>{a}[[i]]//Mean,{i,Length[data[[1]]]}];
reldata=Table[((data[[i]]/.LF:>List)/MeanX)/.List:>LF,{i,Length[data]}];
relMeanX=Table[1,{i,Length[MeanX]}];

Preparation for Reconstruction

Here, the final steps in preparing the algorithm's inputs are performed. "dsmear" is a variable related to the Convex Hull method; it is a measure to suppress the method's tendency to randomly fluctuate in the presence of outliers and noise. If "dsmear" is set to some positive value rather than 0, the method fits an ellipse not to the convex hull, but to the band of width (1±dsmear) around the boundary of the fitted region. "dsmear" should be set to 0 if the "default" Convex Hull method is desired. Otherwise, values between 0.04 - 0.1 usually produce good results. 

If the "Data Normalization" cell was executed, the rest of the cell does not need edited. If it was not executed, then the "mean" and "plist" variables will have to be edited before execution. "mean" must be supplied manually, as an array containing the mean of each of the data's columns. "plist" must be set to "data" rather than "reldata".

dsmear=0.07; 
mean=relMeanX; 
(* Example of "mean" for non-normalized 5-dimensional data: {1.65, .7, 10.6, -.4, 4.11} *)
plist=reldata;
(* plist = data *)

Projection Reconstruction

This section examines each unique two-dimensional projection of the data and fits ellipses to these projections with both the Convex Hull Method and the Covariance Matrix Method. This cell produces graphical output: if the data is d-dimensional, then d*(d-1)*.5 graphs of each two-dimensional projection and its two fitted ellipses are generated.

This cell is completely automatic; no edits are needed regardless of whether or not the data was normalized.

d = Length[plist[[1]]];

Print["Reconstruct the ellipsoid's projections:\naCM (correlation ellipse, solid red); aHull (convex hull, black dashed)"];
{aCM,aHull}={ConstantArray[0,d*(d-1)/2],ConstantArray[0,d*(d-1)/2]};
Module[{x1,y1,x0,y0,dx,dy,xmin,xmax,ymin,ymax,δx,δy,cosphi,phi,s02,xy,oneZScore,a,id},
id=0;
For[ix=1,ix≤d-1,ix++,
For[iy=ix+1,iy≤d,iy++,
Print["Plotting plane x"<>ToString[ix]<>"x"<>ToString[iy]];
id++;

(* Mean-subtract the data and calculate statistics necessary for both ellipse-fitting methods. *)
x1=plist/.LF[a__]:>{a}[[ix]]; 
y1=plist/.LF[a__]:>{a}[[iy]];  
{x0,y0}={mean[[ix]],mean[[iy]]};
Print["The center of the ellipse is at {"<>ToString[x0]<>","<>ToString[y0]<>"}"];
{dx, dy} = {x1-x0, y1 - y0};
{xmin,xmax}={Min[#],Max[#]}&@ x1;
{ymin,ymax}={Min[#],Max[#]}&@ y1;
{δx, δy} ={StandardDeviation[x1],StandardDeviation[y1]};

(* Fit the Covariance Matrix ellipse *)
cosphi=Mean[dx dy]/δx/δy;   
phi=ArcCos[cosphi];
{z1x,z1y} = {nsigma* δx, nsigma* δy};
s02 =((Sin[phi])^2) ;
aCM[[id]]= {{(1/(((δx)^2) s02)), -(cosphi/(δx δy s02))},{-(cosphi/(δx δy s02)),(1/(((δy)^2) s02))}}; 

(* Fit the Convex Hull ellipse *)
Print["-----------------------"];
xy = plist/.LF[a__]:>{{a}[[ix]],{a}[[iy]]}; 
zsmear=1.0+dsmear;
oneZScore= Select[xy,(((#[[1]]>(x0- zsmear z1x)) && (#[[1]]<(x0+zsmear z1x))) && ((#[[2]]>(y0-zsmear z1y)) && (#[[2]]<(y0+zsmear z1y)) ))&]; 
cH=ConvexHullMesh[oneZScore]; 
convexHullPoints=Cases[Show[cH],GraphicsComplex[pts_,___]:>pts][[1]]; 
If[Abs[dsmear]>(10^(-4)),
zsmear=1.0-dsmear;
oneZScore= Select[xy,(((#[[1]]>(x0- zsmear z1x)) && (#[[1]]<(x0+zsmear z1x))) && ((#[[2]]>(y0-zsmear z1y)) && (#[[2]]<(y0+zsmear z1y)) ))&]; 
cH=ConvexHullMesh[oneZScore]; (* finding convex hull of points *)
convexHullPoints=Join[convexHullPoints,Cases[Show[cH],GraphicsComplex[pts_,___]:>pts][[1]]]; 
];
pointsToFit= convexHullPoints/. {x_, y_} :> {x-x0,y-y0};

{d1,d2}={pointsToFit/.{x_,y_}:>{(x^2),x y, (y^2)},pointsToFit/.{x_,y_}:>{1}};
{s11,s12,s21,s22}={Transpose[d1].d1,Transpose[d1].d2,Transpose[d2].d1,Transpose[d2].d2};
{sv11,sv12,sv21,sv22}={linv.Transpose[u1].s11.v1,linv.Transpose[u1].s12,s21.v1,s22};

If[ Det[sv22]==0,
Print["Error: Singular sv22"],
ES=Sort[Eigensystem[sv11-sv12.Inverse[sv22].sv21]//Transpose,#1[[1]]>#2[[1]]&];
NumPositive=Select[Transpose[ES][[1]],#>0&]//Length;
If[NumPositive≠1,
Print["Error: number of positive eigenvalues is ",NumPositive],
av1=Transpose[{ES[[1]][[2]]}];
av2=-Inverse[sv22].sv21.av1;
{a1,a2}={v1.av1,av2};
{X0, Y0} = {0,0};
NormFac = 1;
]
];
a=Transpose[{Join[Transpose[a1][[1]],{0,0,a2[[1]][[1]]}]}]//Flatten;
constant2 = -(((a[[6]]))^(-1)); 
aHull[[id]] ={{constant2 a[[1]], (1/2)constant2 a[[2]]},{(1/2)constant2 a[[2]],constant2 a[[3]]}}; 
(* Visualize the projections and their fitted ellipses *)

parametricPlot1 = ParametricPlot[{x0+ z1x Cos[u],y0+ z1y Cos[u+phi]},{u,0,2π},PlotStyle->RGBColor[1,0,0]];
parametricPlotA=ContourPlot[1-((({{x-x0, y-y0}})).aHull[[id]].(({{x-x0}, {y-y0}})))[[1]][[1]]==0,{x,xmin,xmax},{y,ymin,ymax},ContourStyle->{RGBColor[0.0,0,0],Dashing[{0.06,0.02}],Thickness[0.006]},PlotPoints->100];
parametricPlotB=ContourPlot[(nsigma^2)+((({{x-x0, y-y0}})).aCM[[id]].(({{x-x0}, {y-y0}})))[[1]][[1]]==0,{x,xmin,xmax},{y,ymin,ymax},ContourStyle->{RGBColor[0.0,0.5,0.5],Dashing[{0.02,0.02}],Thickness[0.01]},PlotPoints->100];
pp2=ListPlot[plist/.LF[a__]:>{{a}[[ix]],{a}[[iy]]}, PlotMarkers -> {"●",4},PlotStyle->Black,AspectRatio->1,ImageSize->{500,Automatic},Frame->True,AxesOrigin->{0,0},Axes->True,FrameLabel->{"σ"<>ToString[ix]<>"/⟨σ"<>ToString[ix]<>"⟩","σ"<>ToString[iy]<>"/⟨σ"<>ToString[iy]<>"⟩"}];  
Print[Show[pp2,parametricPlot1,parametricPlotA, 
Graphics[{
RGBColor[0,0,1],Text["★",{x0,y0},{0,0}],
Dashing[{0.02,0.02}],
Line[{{x0-z1x,y0-z1y},{x0+z1x,y0-z1y},{x0+z1x,y0+z1y},{x0-z1x,y0+z1y},{x0-z1x,y0-z1y}}],
}], PlotRange->{{x0-1.2z1x,x0+1.2z1x},{y0-1.2 z1y,y0+1.2z1y}} ,AxesOrigin->{x0,y0},PlotLabel->"Approximate "<>ToString[nsigma]<>"σ countours"]];
](* For iy -> *)
](* For ix -> *)
](* Module -> *)

Ellipsoid Reconstruction and Evaluation

This cell reconstructions two matrices, one per ellipse-fitting method, that each characterize the quadratic form of the reconstructed confidence ellipsoid. It also calculates a “consistency coefficient” for each matrix, a measure of the spread of the coefficients along the diagonal of the ellipses’ characterizing matrices. Finally, the average relative difference between the two matrices are calculated to measure the divergence of the two ellipse-fitting methods. 

This cell is completely automatic; no edits are needed.

(* Initialize necessary variables. *)
i=1;
{aOutCM,aOutHull}={ConstantArray[0,{d,d}],ConstantArray[0,{d,d}]};
{adiagCM,adiagHull}={ConstantArray[{},d],ConstantArray[{},d]};
{daOutCM,daOutHull}={ConstantArray[0,d],ConstantArray[0,d]};
{aconCM,aconHull}={1,1};

(* Calculate the non-diagonal elements of the final matrices. *)
For[ix=1,ix≤d-1,ix++,
For[iy=ix+1,iy≤d,iy++,
{aInHull,aInCM}={-Inverse[aHull[[i]]],-(nsigma^2) Inverse[aCM[[i]]]};
aOutCM[[ix,iy]]=aOutCM[[iy,ix]]=aInCM[[1,2]];
aOutHull[[ix,iy]]=aOutHull[[iy,ix]]=aInHull[[1,2]];
adiagCM[[ix]]=Append[adiagCM[[ix]],aInCM[[1,1]]];
adiagCM[[iy]]=Append[adiagCM[[iy]],aInCM[[2,2]]];
adiagHull[[ix]]=Append[adiagHull[[ix]],aInHull[[1,1]]];
adiagHull[[iy]]=Append[adiagHull[[iy]],aInHull[[2,2]]];
i++
]
]

(* Calculate the diagonal elements of the final matrices, as well as statistics for the consistency coefficient. *)
For[ix=1,ix≤d,ix++,
{aOutCM[[ix,ix]],aOutHull[[ix,ix]]}={Mean[adiagCM[[ix]]],Mean[adiagHull[[ix]]]};
{daOutCM[[ix]],daOutHull[[ix]]}={StandardDeviation[adiagCM[[ix]]],StandardDeviation[adiagHull[[ix]]]};
aconCM*=daOutCM[[ix]]/aOutCM[[ix]][[ix]];
aconHull*=daOutHull[[ix]]/aOutHull[[ix]][[ix]];
]
{aconCM,aconHull}={(((Abs[aconCM]))^(1/d)),(((Abs[aconHull]))^(1/d))}//Chop;

(* Calculate the final matrices. *)
aFinalCM=Inverse[aOutCM];
aFinalHull=Inverse[aOutHull];

(* Compute the average relative difference between the matrices produced by the two methods. *)
delta = 0;
d = Length[plist[[1]]];
Table[
num1 = aFinalCM[[ix]][[iy]];
num2 = aFinalHull[[ix]][[iy]];
difference = Abs[num1-num2];
sum = Abs[(num1+num2)];
If[sum==0, Print["Error: Cannot calculate the relative difference because the sum of these two elements is zero."];,
relative = difference / sum;
delta+= relative;
], {ix, 1, d},{iy,1, ix}];
delta  = delta(2/d/(d+1));

Print["The reconstructed matrices are (Subscript[A, Hull])=",aFinalHull," and  (Subscript[A, cov . mat . ])=",aFinalCM];
Print["The consistency coefficients (average differences in diagonal elements) are (Subscript[acon, Hull])=",aconHull," and (Subscript[acon, CM])=",aconCM];
Print[" The average relative difference between the elements in (Subscript[A, Hull]) and their associated elements in (Subscript[A, cov . mat . ]) is avrd=",delta];