Added cddlib-094h from http://www.inf.ethz.ch/personal/fukudak/cdd_home/
Change-Id: I64519509269e434b1b9ea87c3fe0805e711c0ac9
diff --git a/third_party/cddlib/examples-ml/UnfoldPolytope.m b/third_party/cddlib/examples-ml/UnfoldPolytope.m
new file mode 100755
index 0000000..64dcb22
--- /dev/null
+++ b/third_party/cddlib/examples-ml/UnfoldPolytope.m
@@ -0,0 +1,423 @@
+(*
+
+ Unfolding Convex Polytopes
+
+ Version 1.0 Beta
+ February, 1992
+
+ Copyright (c) 1992 by
+ Makoto Namiki
+
+This package contains Mathematica implementations
+of unfolding 3-dimsional polytope.
+
+This package is copyright 1992 by Makoto Namiki.
+This package may be copied in its entirety for
+nonprofit purposes only. Sale, other than for
+the direct cost of the media, is prohibited.
+This copyright notice must accompany all copies.
+
+The authors make no representations, express or
+implied, with respond to this documentation, of
+the software it describes and contains, including
+without limitations, any implied warranties of
+mechantability or fitness for a particular purpose,
+all of which are expressly disclaimed. The authors
+shall in no event be liable for any indirect,
+incidental, or consequential damages.
+
+This beta release is designed to run under
+Version 1.2 & 2.0 of Mathematica. Any comments,
+bug reports, or requests to get on the
+UnfoldPolytope mailing list should be
+forwarded to:
+
+ Makoto Namiki
+ Department of Information Science,
+ Tokyo Institute of Technology, Tokyo
+ 2-12-1 Oh-okayama, Meguro-ku
+ Tokyo 145, Japan
+
+ +81-3-3726-1111 (Ex. 4131)
+ namiki@titisna.is.titech.ac.jp
+
+*)
+
+BeginPackage["UnfoldPolytope`"]
+
+UnfoldPolytope::usage="UnfoldPolytope[facets] gives graphics of unfoldings of a polytope which is given as a list of facets.";
+
+OrderFacets::usage="OrderFacets[facets] returns a correctly order lists of facets.";
+
+ProperSubsetQ::usage="ProperSubsetQ[t,s] returns True if t is a proper subset of s.";
+
+MakeEdgesFromFacets::usage="MakeEdgesFromFacets[faces] returns a cordinats of lower faces and a adjacency of faces.";
+
+MakeFacetsFromZerosets::usage="MakeFacetsFromZerosets[v,z] returns a list of facets consisting of cordinates of vertices, where v denotes a list of vertices and z denotes a list of zerosets.";
+
+MakeAdjTable::usage="MakeAdjTable[n,l] or MakeAdjTable[n,l,costs] return an adjacency table where l denotes a list of unordered pairs consist of {1,2,...,n} and a costs denotes an adjacency cost.";
+
+BreadthFirstSearch::usage="BreadthFirstSearch[adj_list,n] gives a incidence matrix representing a tree of a graph obtained by breadth first search from a vertex n, where this graph is given as a incidence matrix of adj_list.";
+
+DepthFirstSearch::usage="DepththFirstSearch[adj_list,n] gives a incidence matrix representing a tree of a graph obtained by depth first search from a vertex n, where this graph is given as a incidence matrix of adj_list.";
+
+IntersectionOfList::usage="IntersectionOfList[{l_1,l_2,...l_n}] returns a Intersection[l_1,l_2,...,l_n].";
+
+SelectLeaf::usage="SelectLeaf[adj_Mat] returns {i,j} such that node i is a child of node j for adjacent matrix adj_Mat.";
+
+ProjectVertex::usage="ProjectVetex[v1,v2,v3,v4] returns coordinates which is a projection of v4 onto affine space {v1,v2,v3} around the line {v1,v2}.";
+
+ProjectFaces::usage="ProjectFaces[child,parent,faces,vervect] argument child is a list of faces, parent is a face, and faces is a list of all faces. This procedure project faces of child onto a face of parent by using ProjectVertices.";
+
+IntersectOfFaces::usage="IntersectOfFaces[f1_List,f2_List] returns veticeswhich commonly belong to f1 and f2.";
+
+VectEQ::usage="VectEQ[v1,v2] returns True if v1 == v2 otherwise False";
+
+Unfold::usage="Unfold[sptree,tree,maxf] open one of a leaf of Tree of Facet.";
+
+Zonotope::usage="Zonotope[vlist] returns a list extreme points of zonotope which is defined by vlist, a list of vectors.";
+
+VerticalVector::usage="VerticalVector[maxf_List] returns a list of Vertical vector of each face.";
+
+OrderFace::usage="OrderFace[facet]: reorder a facet correctly.
+ where argument facet consits of {v_1,v_2,..,v_n}
+ ,v_i is a cordinat of a vetex.";
+
+Begin["`private`"]
+
+UnfoldPolytope[facets_List]:=
+ Block[{odfacets,edg,faAdj,vertices,veAdj,
+ t,tree,cotree,sptree,i,tr,vervec},
+ odfacets = facets;
+ {edg,faAdj} = MakeEdgesFromFacets[odfacets];
+ vertices = Union[Flatten[facets,1]];
+ veAdj = Flatten[{Position[vertices,#[[1]]],
+ Position[vertices,#[[2]]]}]& /@ edg;
+ t = BreadthFirstSearch[MakeAdjTable[Length[vertices],veAdj],1];
+ tree = Flatten[Position[veAdj,#]& /@ Position[t,1]];
+ cotree = Complement[Table[i,{i,1,Length[edg]}],tree];
+ sptree = MakeAdjTable[Length[odfacets],faAdj[[cotree]]];
+ tr = Table[{i},{i,Length[facets]}];
+ vervec = VerticalVector[odfacets];
+ Show[Graphics3D[Polygon /@ odfacets],Boxed->False];
+ Do[{sptree,tr,odfacets,vervec} = Unfold[sptree,tr,N[odfacets],vervec];
+ Show[Graphics3D[(Polygon /@ odfacets)],Boxed->False];
+ ,{i,1,Length[odfacets]-1}];
+ Show[Graphics3D[(Polygon /@ odfacets)],Boxed->False,ViewPoint->vervec[[1]][[1]]];
+ ];
+
+
+OrderFacets[facets_List]:=
+ Block[{i,edges,facetsAdj,f,out,cyc,cand,cand2cyc,ff},
+ {edges,facetsAdj} = MakeEdgesFromFacets[facets];
+ out = {};
+ For[i=1,i<=Length[facets],i++,
+ f = edges[[Transpose[Position[facetsAdj,i]][[1]]]];
+ cand = Drop[f,1];
+ cyc = {First[f]};
+ ff = {};
+ While[(cand2cyc = Select[cand,(Length[Intersection[Last[cyc],#]]==1)&]) != {},
+ AppendTo[ff,Intersection[Last[cyc],cand2cyc[[1]]][[1]]];
+ AppendTo[cyc,cand2cyc[[1]]];
+ cand = Complement[cand,{cand2cyc[[1]]}];
+ ];
+ AppendTo[ff,Intersection[Last[cyc],First[cyc]][[1]]];
+ AppendTo[out,ff];
+ ];
+ Return[out];
+ ];
+
+NonNegativeVectorQ[v_List]:=
+ Block[{i},
+ For[i=1, i<= Length[v], ++i,
+ If [v[[i]] < 0, Return[False]];
+ ];
+ Return[True];
+ ];
+
+
+MakeEdgesFromFacets[l_List] :=
+ Block[{i,j,lf={},adj={}},
+ For[i=1,i<Length[l],i++,
+ For[j=i+1,j<=Length[l],j++,
+ If[Length[Intersection[l[[i]],l[[j]]]] == 2,
+ AppendTo[lf,Intersection[l[[i]],l[[j]]]];
+ AppendTo[adj,{i,j}];
+ ];
+ ];
+ ];
+ Return[{lf,adj}];
+ ];
+
+MakeFacetsFromZerosets[v_List,z_List]:=
+ Block[{i,j,flist,maxf,flag,out},
+ flist = Transpose[Position[z,#]][[1]]& /@ Union[Flatten[z]];
+ maxf = {};
+ For[i=1,i<=Length[flist],i++,
+ flag = 0;
+ For[j=1,j<=Length[flist],j++,
+ If[ProperSubsetQ[flist[[i]],flist[[j]]],
+ flag = 1,
+ If[SubsetQ[flist[[i]],flist[[j]]] && i<j,
+ flag = 1;
+ ];
+ ];
+ ];
+ If[flag == 0,
+ AppendTo[maxf,flist[[i]]];
+ ];
+ ];
+ (* Return[maxf]; *)
+ out = v[[#]]& /@ (MakeOrder[#,z]& /@ maxf);
+ Return[out];
+ ];
+
+MakeOrder[ver_List,zerova_List]:=
+ Block[{candidates, cycle ,candidate2cycle},
+ candidates = Drop[ver,1];
+ cycle = { First[ver] };
+ (* when one is list, comparing needs other is list *)
+ While[{} != (candidate2cycle = Select[candidates,AdjQ[Last[cycle],#,zerova] &]),
+ AppendTo[cycle,candidate2cycle[[1]]];
+ candidates = Complement[candidates,candidate2cycle[[{1}]]];
+ ];
+ While[{} != (candidate2cycle = Select[candidates,AdjQ[First[cycle],#,zerova] &]),
+ PrependTo[cycle,candidate2cycle[[1]]];
+ candidates = Complement[candidates,candidate2cycle[[{1}]]];
+ ];
+ cycle
+ ];
+
+
+AdjQ[i_Integer,j_Integer,l_List]:=
+ Block[{subface},
+ subface=Intersection[l[[i]],l[[j]]];
+ Length[Select[l,SubsetQ[subface,#]&]]==2];
+
+ProperSubsetQ[s_List,t_List]:=
+ Length[s]==Length[Intersection[s,t]] && Length[s]<Length[t];
+
+SubsetQ[s_List,t_List]:=
+ Length[s]==Length[Intersection[s,t]];
+
+MakeAdjTable[n_Integer,l_List]:=
+ Block[{t=Table[Table[0,{n}],{n}],i},
+ For[i=1,i<=Length[l],i++,
+ t[[l[[i,1]],l[[i,2]]]] = 1;
+ t[[l[[i,2]],l[[i,1]]]] = 1;
+ ];
+ Return[t];
+ ];
+
+MakeAdjTable[n_Integer,l_List,c_List]:=
+ Block[{t=Table[Table[0,{n}],{n}],i},
+ For[i=1,i<=Length[l],i++,
+ t[[l[[i,1]],l[[i,2]]]] = c[[i]];
+ t[[l[[i,2]],l[[i,1]]]] = c[[i]];
+ ];
+ Return[t];
+ ];
+
+BreadthFirstSearch[adj_?MatrixQ,n_Integer]:=
+ Block[{out={},search,find={},output,
+ i,obj,cd,l=Length[adj]},
+ search = {n};
+ While[search != {},
+ obj = search[[1]];
+ search = Drop[search,1];
+ cd = {};
+ For[i=1, i<=l, i++,
+ If[adj[[obj,i]] != 0, AppendTo[cd,i]];
+ ];
+ cd = Complement[cd,search];
+ cd = Complement[cd,find];
+ For[i=1, i<=Length[cd],i++,
+ AppendTo[search,cd[[i]]];
+ AppendTo[out,{obj,cd[[i]]}];
+ AppendTo[out,{cd[[i]],obj}];
+ ];
+ AppendTo[find,obj];
+ ];
+ Return[MakeAdjTable[l,out]];
+ ];
+
+DepthFirstSearch[adj_?MatrixQ,n_Integer]:=
+ Block[{out={},search,find={},output,
+ i,obj,cd,l=Length[adj]},
+ search = {n};
+ While[search != {},
+ obj = search[[1]];
+ search = Drop[search,1];
+ cd = {};
+ For[i=1, i<=l, i++,
+ If[adj[[obj,i]] != 0, AppendTo[cd,i]];
+ ];
+ cd = Complement[cd,search];
+ cd = Complement[cd,find];
+ search = Flatten[Append[cd,search]];
+ For[i=1, i<=Length[cd],i++,
+ AppendTo[out,{obj,cd[[i]]}];
+ AppendTo[out,{cd[[i]],obj}];
+ ];
+ AppendTo[find,obj];
+ ];
+ Return[MakeAdjTable[l,out]];
+ ];
+
+
+IntersectionOfList[l_List]:=
+ Block[{},
+ If[Length[l]==1,Return[l[[1]]],
+ Return[Intersection[l[[1]],
+ IntersectionOfList[Drop[l,{1}]]]];
+ ];
+ ];
+
+Maximals[l_List]:=
+ Block[{i,maximals={}},
+ Do[If[MaximalQ[i,l],
+ AppendTo[maximals,l[[i]]]
+ ], {i,Length[l]}
+ ];
+ maximals]
+
+MaximalQ[i_Integer,l_List]:=
+ Block[{candidate},
+ candidate=l[[i]];
+ Length[Union[Select[l,ProperSubsetQ[candidate,#]&]]]==0]
+
+SubsetQ[s_List,t_List]:=
+ Length[s]==Length[Intersection[s,t]];
+
+Is[s_List]:=
+ Block[{},
+ If [Length[s] == 1,
+ Return[s[[1]]],
+ Return[Intersection[s[[1]],Is[Drop[s,1]]]]
+ ];
+ ];
+
+SelectLeaf[g_List]:=
+ Block[{size,i,j,deg,child,par,cp},
+ size = Length[g];
+ For[i=2,i<=size,i++,
+ deg = 0;
+ child = i;
+ For[j=1,j<=size,j++,
+ If[g[[i]][[j]] > 0,
+ deg = deg + 1;
+ par = j;
+ ];
+ ];
+ If[deg == 1,
+ Return[{child,par}];
+ ];
+ ];
+ ]
+
+ProjectFaces[child_List,par_Integer,face_List,vervec_List]:=
+ Block[{v1,v2,nch,i,fl},
+ {v1,v2} = IntersectOfFaces[face[[child[[1]]]],face[[par]]];
+ For[i=1, i<=Length[face[[par]]], ++i,
+ If[!VectEQ[face[[par]][[i]],v1] && !VectEQ[face[[par]][[i]],v2],
+ v3 = face[[par]][[i]];
+ Break[];
+ ];
+ ];
+ fl = {};
+ For[i=1, i<=Length[child], i++,
+ AppendTo[fl,ProjectVertex[v1,v2,v3,#,vervec]& /@ face[[child[[i]]]]];
+ ];
+ Return[fl];
+ ]
+
+IntersectOfFaces[f1_List,f2_List]:=
+ Block[{i,j,lvect},
+ lvect = {};
+ For[i=1,i<=Length[f1],i++,
+ For[j=1,j<=Length[f2],j++,
+ If[VectEQ[f1[[i]],f2[[j]]],AppendTo[lvect,f1[[i]]]];
+ ];
+ ];
+ Return[{lvect[[1]],lvect[[2]]}];
+ ]
+
+VectEQ[v1_List,v2_List]:=
+ Block[{i},
+ For[i=1,i<=Length[v1],i++,
+ If[v1[[i]] != v2[[i]],Return[False]];
+ ];
+ Return[True];
+ ]
+
+
+ProjectVertex[v1_List,v2_List,v3_List,v4_List,a_List]:=
+ Block[{x,y},
+ If [N[a[[1]].v4] <= a[[2]] + 10^(-10),
+ x = (v2-v1).(v2-v4)/(v2-v1).(v2-v1)*v1 + (v1-v2).(v1-v4)/(v1-v2).(v1-v2)*v2;
+ y = (v2-v1).(v2-v3)/(v2-v1).(v2-v1)*v1 + (v1-v2).(v1-v3)/(v1-v2).(v1-v2)*v2;
+ Return[x + Sqrt[(v4-x).(v4-x)/(y-v3).(y-v3)]*(y-v3)]
+ ];
+ If [N[a[[1]].v4] >= a[[2]] - 10^(-10),
+ x = (v2-v1).(v2-v4)/(v2-v1).(v2-v1)*v1 + (v1-v2).(v1-v4)/(v1-v2).(v1-v2)*v2;
+ y = (v2-v1).(v2-v3)/(v2-v1).(v2-v1)*v1 + (v1-v2).(v1-v3)/(v1-v2).(v1-v2)*v2;
+ Return[x + Sqrt[(v4-x).(v4-x)/(y-v3).(y-v3)]*(v3-y)]
+ ];
+ ]
+
+Unfold[sptree1_List,tree1_List,maxf1_List,vervec_List]:=
+ Block[{j,child,par,posch,ch,of,pospa,sptree=sptree1,tree=tree1,maxf=maxf1,vv=vervec},
+ {child, par} = SelectLeaf[sptree];
+ sptree[[child,par]] = -1;
+ sptree[[par,child]] = -1;
+
+ {posch} = Transpose[Position[tree,child]][[1]];
+ ch = tree[[posch]];
+ of = ProjectFaces[ch,par,maxf,vv[[par]]];
+ For[j=1, j<=Length[ch],j++,
+ maxf[[ch[[j]]]] = of[[j]];
+ ];
+ {pospa} = Transpose[Position[tree,par]][[1]];
+ tree[[pospa]] = (AppendTo[tree[[pospa]],#]& /@ ch)[[Length[ch]]];
+ tree = Drop[tree,{posch}];
+ Return[{sptree,tree,maxf,vv}];
+ ]
+
+Zonotope[vectors_List]:=
+ Block[{vlist},
+ vlist = Sum[(#)[[i]],{i,Length[#]}]& /@ ((vectors * (#)&) /@ Pow[Length[vectors]]);
+ Return[vlist];
+ ]
+
+Pow[n_Integer]:=
+ Block[{p,a,b},
+ If[n == 1, Return[{{1},{-1}}]
+ ];
+ p = Pow[n-1];
+ a = Append[#,1]& /@ p;
+ b = Append[#,-1]& /@ p;
+ Return[(AppendTo[a,#]& /@ b)[[Length[b]]]];
+ ]
+
+NonEmptyFaces[zerova_List]:=
+ Map[Function[Transpose[Position[zerova,#]][[1]]],Union[Flatten[zerova]]]
+
+VerticalVector[maxf_List]:=
+ Block[{i,vervec,a1,a2,x,y,z,b,v,vv={}},
+ For[i=1,i<=Length[maxf],i++,
+ a1 = maxf[[i]][[3]]-maxf[[i]][[2]];
+ a2 = maxf[[i]][[1]]-maxf[[i]][[3]];
+ Clear[x,y,z];
+ vervec = {x,y,z};
+ {vervec} = vervec /. Solve[{a1.vervec == 0, a2.vervec == 0}];
+ x = 1; y = 1; z = 1;
+ b = vervec.maxf[[i]][[1]];
+ v = Complement[Union[Flatten[maxf,1]],maxf[[i]]][[1]];
+ If[v.vervec > b, vervec = -vervec; b = -b];
+ AppendTo[vv,{vervec,b}];
+ ];
+ Return[vv];
+ ]
+
+End[]
+EndPackage[]