third_party/cddlib/examples-ml/UnfoldPolytope.m - RealtimeRoboticsGroup/test - Gitiles

 (*

         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[]
	(*

	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[]