(*Number of points used in the computation. Use max 16 for trial purposes and max 22 for the computation not to last more than 1 day on average PC.*)
varNOfPoints=12;

(*The following option is needed for experimental states 7 and 8, the product of the numbers should be equal to varNOfPoints*)
varLatticeDimensions={3,4};

(*Important note: varDimensions fields correspond to the subsystems from left to right (because also )*)
varDimensions=Table[2,{i,varNOfPoints}];

(*Check, if length of the list Dimensions is equal to the number of points.*)
"Number of dimensions is equal to the number of points: "<>ToString[varNOfPoints==Length[varDimensions]]

(*Check, if the dimensions are integer*)
"Dimensions are integer: "<>ToString[ArrayQ[varDimensions,_,IntegerQ]]

(*Number of basis states, characterizes well the overall complexity of the computation*)
varComplexity=Times@@varDimensions;

(*Converts integer representing an ordinal number of a basis vector to that basis vector*)
funcIntegerToList[i_Integer]:=Module[{locvarInt=i,locvarOutput={},j=1},
While[
j≤varNOfPoints,
PrependTo[locvarOutput,Mod[locvarInt,varDimensions[[-j]]]];
locvarInt=Quotient[locvarInt,varDimensions[[-j]]];
j++
];
locvarOutput
];

(*Generates all basis vectors*)
varBasis=funcIntegerToList[#]&/@(Range[0,varComplexity-1]);

(*Converts a basis vector to its ordinal number*)
funcListToInteger[l_List]:=Total[Table[l [[-k]](varDimensions[[-k]]^(k-1)),{k,1,Length[l]}]];

(*The following function computes the reduced matrix based on the full density matrix corresponding to the state 'varState'. It computes the trace over all of the basis subsystems but those specified by the integers in the input list 'points'*)

funcReduceMatrixFromState[points_List,state_List]:=Module[{locvarOrderedList=Sort[points,Greater],locvarStartIterators={},locvarIterators={},locvarIndexOfTraced=Complement[Range[varNOfPoints],points],locvarStart,locvarStartJumps,locvarVectorJumps,locvarVector,locvarOutput},

If[locvarIndexOfTraced=={},Outer[Times,state,Conjugate[state]],

locvarStartIterators=Table[{j[i],0,varDimensions[[i]]-1},{i,locvarIndexOfTraced}];

locvarStartJumps=Total@((j[#]Times@@varDimensions[[1;;#-1]])&/@locvarIndexOfTraced);

locvarStart=Flatten[Table[locvarStartJumps,Evaluate[Sequence@@locvarStartIterators]]];

locvarIterators=Table[{j[i],0,varDimensions[[i]]-1},{i,points}];

locvarVectorJumps=Total@((j[#]Times@@varDimensions[[1;;#-1]])&/@points);

locvarVector=Flatten[Table[#+locvarVectorJumps,Evaluate[Sequence@@locvarIterators]]]&/@locvarStart;

locvarOutput = Total[Outer[Times,#,Conjugate[#]]&/@(state[[#]]&/@(locvarVector+1))];

locvarOutput
]
];

(*Computes the distance based on information. There is a hard set maximal distance, in order to avoid infinites in the cases of degenerate results. Such results are not interesting in any case, but in this way one avoid clutter with infinites and indeterminate expressions*)
funcDistanceFunction[info_,order_]:=If[Chop[info[[order,1]]]==0,30,Chop[-Log[info[[order,1]]/(2Log[2])]]]

(*Measurment of the time the computation of the reduced matrices starts here (there are some commands before the actual tracing starts, but these take little time to influence the result).*)
varTime=AbsoluteTime[];

(*Library of states*)

(*(1) - State with fully random amplitude and phase for each basis vector.*)
varStateAllRandom=Normalize[Table[(#1 Exp[2 Pi I #2])&@@RandomReal[{0,1},2],{i,varComplexity}]];

(*(2) - State with uniform amplitude and random phase*)
varStateRandomPhase=Normalize[Table[Exp[2 Pi I #]&@@RandomReal[{0,1}],{i,varComplexity}]];

(*(3) - State with linearly rising amplitude and constant phase*)
varStateLinearAmp=Normalize[Table[i,{i,varComplexity}]];

(*(4)- State with gaussian distribution of amplitudes and constant phase*)
varStateGaussian=Normalize[Table[N[PDF[NormalDistribution[Quotient[varComplexity,2],Quotient[varComplexity,6]]][i]],{i,varComplexity}]];

(*Excluded (5)- State with constant amplitudes but opposite phases between each two neighbours.
Reason for excluding - gives zero entropy*)
varStateOppositePhase=Normalize[Table[(-1)^i,{i,varComplexity}]];

(*Excluded (6)- State with constant phase and amplitude
Reason for excluding - gives zero entropy*)
varStateConstant=Normalize[Table[1,{i,varComplexity}]];

(*(7) - Attempting to find a state for lattice-like positioning of the points with binary Hilbert space at each point.*)
varLatticeBasis=Module[{locvarInstructions},

(*This contains the instructions, which "spins" to flip for each vector creating the total vector used for further computation*)
locvarInstructions=Flatten[Table[

varTMP1=i+ varLatticeDimensions[[1]]j;

{1+varTMP1->1, (*choosen point*)
1+varLatticeDimensions[[1]]j+Mod[varTMP1+1,varLatticeDimensions[[1]]]->1,(*the point to the right*)
1+varLatticeDimensions[[1]]j+Mod[varTMP1-1,varLatticeDimensions[[1]]]->1, (*the point to the left*)1+(j-Mod[j,varLatticeDimensions[[2]]]+Mod[j+1,varLatticeDimensions[[2]]])varLatticeDimensions[[1]]+i->1,(*the point above*)
1+(j-Mod[j,varLatticeDimensions[[2]]]+Mod[j-1,varLatticeDimensions[[2]]])varLatticeDimensions[[1]]+i->1(*the point below*)},

{i,0,varLatticeDimensions[[1]]-1},{j,0,varLatticeDimensions[[2]]-1}]

,1];

Reverse[ReplacePart[Table[0,{i,Times@@varLatticeDimensions}],#]]&/@locvarInstructions

];



varExperimentalState1=Normalize[ReplacePart[Table[0,{i,varComplexity}],#->1&/@(funcListToInteger/@varLatticeBasis)]];

(*(8) - Attempting to find a state for lattice-like positioning of the points. (Alternate approach - the central point is omitted)*)
varLatticeBasis=Reverse[ReplacePart[Table[0,{i,Times@@varLatticeDimensions}],#]]&/@Flatten[Table[

varTMP1=i+ varLatticeDimensions[[1]]j;

{1+varLatticeDimensions[[1]]j+Mod[varTMP1+1,varLatticeDimensions[[1]]]->1,(*the point to the right*)
1+varLatticeDimensions[[1]]j+Mod[varTMP1-1,varLatticeDimensions[[1]]]->1, (*the point to the left*)1+(j-Mod[j,varLatticeDimensions[[2]]]+Mod[j+1,varLatticeDimensions[[2]]])varLatticeDimensions[[1]]+i->1,(*the point above*)
1+(j-Mod[j,varLatticeDimensions[[2]]]+Mod[j-1,varLatticeDimensions[[2]]])varLatticeDimensions[[1]]+i->1(*the point below*)},

{i,0,varLatticeDimensions[[1]]-1},{j,0,varLatticeDimensions[[2]]-1}],

1];


varExperimentalState2=Normalize[ReplacePart[Table[0,{i,varComplexity}],#->1&/@(funcListToInteger/@varLatticeBasis)]];

(*(9) - Weighted state: Each basis vector is given the weight corresponding to the number of 1's in its binary representation.*)
varWeightedState=Normalize[Table[Total[varBasis[[i]]],{i,Length[varBasis]}]];

(*(10) - Modified weighted state: Chosen basis vectors are attributed weight corresponding to their binary representation.*)
varWeightedState2=Table[0,{i,varComplexity}];
(varWeightedState2[[#+1]]=1/#)&/@Table[funcListToInteger@Table[If[j≥i,1,0],{j,varNOfPoints}],{i,varNOfPoints}];
varWeightedState2=Normalize[varWeightedState2];

(*(11) - Modified weighted state with random phase: Chosen basis vectors are attributed weight corresponding to their binary representation and random phase.*)
varWeightedState3=Table[0,{i,varComplexity}];
(varWeightedState3[[#+1]]=(1/# )Exp[I RandomReal[{0,2Pi}]])&/@Table[funcListToInteger@Table[If[j≥i,1,0],{j,varNOfPoints}],{i,varNOfPoints}];
varWeightedState3=Normalize[varWeightedState3];

(*In the following variable all the states are stored in denoted order, so the state actually chosen for the computation can be easily reffered to in the following.*)
varStateLibrary={varStateAllRandom,varStateRandomPhase,varStateLinearAmp,varStateGaussian,varStateOppositePhase,varStateConstant,varExperimentalState1,varExperimentalState2,varWeightedState,varWeightedState2,varWeightedState3};

varChosenState=1;


(*Choosing the state for the following computation*)
varState=varStateLibrary[[varChosenState]];

(*The two following calls compute the entropy of a single point and between two points. For each computation the result is a two-field array, first field containing the result itself and the other indices of the subsystems the entropy corresponds to.*)
varPairwiseEntropy=Chop[{Total[N[-Log[#]#]&/@DeleteCases[Eigenvalues[funcReduceMatrixFromState[#,varState]],N[0]|0]],#}&/@Subsets[Range[varNOfPoints],{2}]];
varPointwiseEntropy =Chop[{Total[N[-Log[#]#]&/@DeleteCases[Eigenvalues[funcReduceMatrixFromState[#,varState]],N[0]|0]],#}&/@Subsets[Range[varNOfPoints],{1}]];

(*The following function computes the pairwise information for all pairs while preserving also the information of which pair does the result correspond to.*)
funcComputeInformation[pairwise_List,pointwise_List]:=Module[{locvarOutput={}},
For[k=1,k≤Length[pairwise],k++,
AppendTo[locvarOutput,{-pairwise[[k,1]]+pointwise[[pairwise[[k,2,1]],1]]+pointwise[[pairwise[[k,2,2]],1]],pairwise[[k,2]]}];
];
locvarOutput
];

varInformationOfPairs= funcComputeInformation[varPairwiseEntropy,varPointwiseEntropy];

(*This command prepeares the table of the distances to be filled in in the next step.*)
varDistances=Table[0,{i,varNOfPoints},{j,varNOfPoints}];

(*This fills in the table of varDistances symmetrically. Currently the formula is distance=-Ln[mutual information].*)
For[k=1,k<=Length[varInformationOfPairs],k++,
varDistances[[varInformationOfPairs[[k,2,1]],varInformationOfPairs[[k,2,2]]]]=varDistances[[varInformationOfPairs[[k,2,2]],varInformationOfPairs[[k,2,1]]]]=funcDistanceFunction[varInformationOfPairs,k];
];

varDistances//MatrixForm;

"The calculation took "<>ToString[AbsoluteTime[]-varTime]<>" seconds"

(*Choose, which point will be the origin (can be either set or randomly chosen). The number identifies the point by its column/row in the table of distances.
Choosing the point in such a way can lead to translation of the final set of points after the computation, but should not affect the result otherwise.*)
varOriginIndex=1;

(*Definition of the triple product function accoring to the article.*)
funcTP[dist_,i_Integer,j_Integer,k_Integer]:=1/2 (dist[[i,k]]^2+dist[[j,k]]^2-dist[[i,j]]^2);

(*A boole to check later, if the dimension of the space into which the embedding will be done, will have been computed during the calculation or will be preset*)
varMaxDimComputed=True;

(*Determination of volumes of all simplices up to the order of the number of the vertices. From the end:
The last Range creates the set of all available orders of the simplices.
Then the Subsets command takes appropriate n-tuples of vertices while making sure to not pick origin,or otherwise the volume would need to be 0.
The the main function produces the Table of triple products between all pairs from given n-tuple and the origin and appropriate Determinant, square Root and Normalization is done
The results are rounded (Chopped) for very small values, for the next steps to be sensible.*)
varNOfSamplePoints=varNOfPoints-1;
"Number of sample points ≤ total number of points: "<>ToString[varNOfSamplePoints≤ varNOfPoints]

(*Sample points indices*)
varSamplePoints=RandomSample[DeleteCases[Range[varNOfPoints],varOriginIndex],varNOfSamplePoints];

(*Simplices volumes*)
varSimplexVolume=Chop[(1/#1!(Sqrt[Det[Table[funcTP[varDistances,i,j,varOriginIndex],{i,#1},{j,#1}]]]&/@Subsets[varSamplePoints,{#1}])&/@Range[varNOfSamplePoints]),10^-7];

varCounter=0;

"Points not satisfying the triangle inequality (if any):"
varNote=Table[
If[varDistances[[i,j]]+varDistances[[j,k]]≥ varDistances[[i,k]],True,Print[{i,j,k,varCounter}];varCounter++;False],
{i,varNOfPoints},{j,varNOfPoints},{k,varNOfPoints}];

(*Check if triangle inequality is satisfied.*)
boolTriangleInequality=And@@Flatten[varNote];
"Triangle equality is satisfied: "<>ToString[boolTriangleInequality]

(*Check if all volumes are real.*)
boolZeroImaginaryPart=Map[Im[#1]==0&,varSimplexVolume,{2}];
boolFlatGraph=And@@@boolZeroImaginaryPart;
varHighestFlatDimension=Module[{locvarOutput=0},
While[ locvarOutput<Length[boolFlatGraph],
If[boolFlatGraph[[locvarOutput+1]]==True,locvarOutput++,Break[]] ;
];
locvarOutput
];

"The highest flat dimension of the graph is: "<>ToString[varHighestFlatDimension]

(*Determination of the highest order of simplices for which there exists at least one with nonzero volume. It creates a List and into each field it puts either the ordinal number of the field, if there exists a simplex of nonzero volume for given order, or 0, if all simplices of given order have 0 volume. Then it takes the maximum as the dimension of the space, into which the points are to be embedded. It also writes this dimension in the end.*)
varNonzeroSimplices=Table[If[varSimplexVolume[[i]]≠Array[0&,Length[varSimplexVolume[[i]]]],i,0],{i,varNOfSamplePoints}];
varMaxDim=Module[{locvarOutput={},k=0},
While[varNonzeroSimplices[[k+1]]≠ 0 &&boolFlatGraph[[k+1]],
k++;
];
k
];
varMaxDimComputed=True;

(*For reason of caluclational difficulty the acquisition of the dimension from the simplices volume can be skipped, but then the dimension of the space into which we want to embedd the points needs to be specified manually*)
varPresetDim=False;
varMaxDim=If[varMaxDimComputed,varMaxDim,varPresetDim];

"Maximal dimension is "<>ToString[varMaxDim]

(*The basis simplex is chosen to be the one, whose points are in distance table specified by rows/columns {1,2,3,...} with possibly skipping the origin, if needed.
The function then assings to each point its coordinates in n-dimensional space within the basis given by the above chosen simplex.*)
varBasisSimplex=DeleteCases[Range[varNOfPoints],varOriginIndex][[1;;varMaxDim]];

funcF[x_Integer]:=Table[funcTP[varDistances,x,varBasisSimplex[[i]],varOriginIndex],{i,varMaxDim}];
varCoordinates1=Table[funcF[i],{i,varNOfPoints}]\[Transpose];

(*Turning the coordinates back to the cartesian coordinates
1st step: the metric of the current space is g=<x_i,x_j,x_0>^-1, i, j goes through the vertices of the reference simplex.*)

varTableOfTP=varCoordinates1\[Transpose][[varBasisSimplex]];

(*step 2: finding the transformation matrix to cartesian space (basically, (u,v)=u.g.v, where u,v are vectors with coordinates with respect to the basis simplex' edges. What we do is, that we find R such, that g=RR\[Transpose] and so (u,v)=u.R.R\[Transpose].v and so R\[Transpose].u and R\[Transpose].v are coordinates of u and v with respect to cartesian unit basis. This is called the Cholesky Decomposition*)
varTransformationMatrix=CholeskyDecomposition[(#+Transpose[#])/2&@Inverse[varTableOfTP]];

(*step 3: applying the transformation to acquire the coordinates of the points in euclidean space with unit metric. The transpose is there in order to have the coordinates in rows for easier computation.*)
varCoordinates2=(varTransformationMatrix.varCoordinates1)\[Transpose];

(*Move the points to the center, so that the average of each of their coordinates is 0*)

varCoordinatesTranslation=Table[-Total[varCoordinates2[[All,i]]]/varNOfPoints,{i,varMaxDim}];
varCoordinates2=Map[#+varCoordinatesTranslation&,varCoordinates2,1];

varGraphOfEmbedding3D=If[varMaxDim==3,ListPointPlot3D[varCoordinates2,AspectRatio->1,ImageSize->Large]]
varGraphOfEmbedding2D=If[varMaxDim==2,ListPlot[varCoordinates2,AspectRatio->1,ImageSize->Large]]

(*Implementation of the multidimensional scaling*)
varB=Table[-1/2(varDistances[[i,j]]^2-Total[#^2&/@varDistances[[i]]]/varNOfPoints-Total[#^2&/@varDistances[[All,j]]]/varNOfPoints+Sum[varDistances[[k,l]]^2,{k,varNOfPoints},{l,varNOfPoints}]/varNOfPoints^2),{i,varNOfPoints},{j,varNOfPoints}];

varESystem=Eigensystem[varB];
(*Note: in transpose form the Eigensystem produces a list of pairs, 1st member of the pair being an eigenvalue and 2nd the corresponding eigenvector. So the zero eigenvalue and their eigenvectors can be easily deleted.*)
varESystem=DeleteCases[Chop[varESystem\[Transpose]],{0,_}]\[Transpose];

varCoordinatesFromMDS=Table[Sqrt[varESystem[[1,i]]]Normalize[varESystem[[2,i]]],{i,Length[varESystem[[1]]]}]\[Transpose];

(*Comparison of results of the coordinates. Each ROW! corresponds to the coordinates of a point.*)
Chop[varCoordinates2]//MatrixForm
Chop[varCoordinatesFromMDS]//MatrixForm

(*The two following lines show that the embeddings are centered at origin*)
Table[Chop[Total[varCoordinates2[[All,i]]]],{i,Length[varCoordinates2[[1]]]}]
Table[Chop[Total[varCoordinatesFromMDS[[All,i]]]],{i,Length[varCoordinatesFromMDS[[1]]]}]

(*Computing the distances from origin for all available embedding methods.*)
varDistancesFromOrigin1=Chop[Total[#^2&/@varCoordinates2,{2}]]
varDistancesFromOrigin2=Chop[Total[#^2&/@varCoordinatesFromMDS,{2}]]

(*Comparison of distances of the points from the origin for different embeddings.*)
varComparisonOfDistances=Table[varDistancesFromOrigin1[[i]]/varDistancesFromOrigin2[[i]],{i,Length[varCoordinatesFromMDS[[1]]]}]

(*Computing the distorsion. If n-th number from left is sufficiently small, n is the number of dimensions of the space, into which the points are well embedded*)
varDistortionParameter=Table[1-Total[varESystem[[1,1;;i]]]/Total[varESystem[[1]]],{i,Length[varESystem[[1]]]}]