(* ::Package:: *)

(***
    this document contains some routines to compute the persistent rational
    cohomology of point clouds.  it's not super fast, but it seems to do kind of OK.

    the main goal is to produce minimal free cDGA models for persistent rational
    cohomology (\[AGrave] la Sullivan); everything else seen here is pretty standard.

    by: Matt Pancia, Eric Peterson
***)

(*BeginPackage["PersistentRationalHomotopy`"];*)

(* documentation goes here *)
BlockSum::usage =
	"BlockSum[A,B,C,D] takes the block sum {{A, B}, {C, D}}. Allows for null matrices.";

(* set up life *)
Unprotect[DGAMonomial, DGADiff, DGAMap];
Clear[Evaluate[Context[] <> "*"]]
$RecursionLimit=1024;

(* eventually we will want to mark most everything private. *)
(*Begin["Private`"];
End[];*)


(* some utilities *)

DPrint[expr_] := (Print[expr]; expr)
TopMinor[matrix_] :=
	matrix[[1 ;; Length[matrix]-1, 1 ;; Length[matrix]-1]]
DeltaColumn[pos_,exp_,len_] :=
	Join[Table[0, {i, 1, pos-1}], {(-1)^exp}, Table[0, {i, pos+1, len}]]
NullMatrix[m_] := m===0 || m==={}
BlockSum[A_,B_,C_,D_] := 
	Which[
		NullMatrix[A]&&NullMatrix[B]&&NullMatrix[C], D,
		NullMatrix[A]&&NullMatrix[B]&&NullMatrix[D], C,
		NullMatrix[A]&&NullMatrix[C]&&NullMatrix[D], B,
		NullMatrix[B]&&NullMatrix[C]&&NullMatrix[D], A,
		NullMatrix[A]&&NullMatrix[B], ArrayFlatten[{{C,D}}],
		NullMatrix[C]&&NullMatrix[D], ArrayFlatten[{{A,B}}],
		NullMatrix[A]&&NullMatrix[C], ArrayFlatten[{{B},{D}}],
		NullMatrix[B]&&NullMatrix[D], ArrayFlatten[{{A},{C}}],
		True, ArrayFlatten[{{A,B},{C,D}}]]

(* TODO: it would be nice if this could be rewritten to work with SparseArrays.  the naive thing of Fold[]ing through and looking for columns that increase matrix rank is not sufficient.  Fold thinks that a SparseArray is an atomic expression where it was expecting a list. *)
MatrixImage[mat_] :=
	If[mat==={} || mat===0, {},
		With[
			{step = Select[RowReduce[Transpose[mat]], Or @@ (# != 0 &/@ #)&]},
			If[step=={}, {}, Transpose[step]]]]

(* from Leonid Shifrin *)
ClearAll[debug];
SetAttributes[debug, HoldAll];
debug[code_] :=
 Internal`InheritedBlock[{Message},
   Module[{inMessage},
     Unprotect[Message];        
     Message[args___] /; ! MatchQ[First[Hold[args]], _$Off] :=
       Block[{inMessage = True},
         Print[{
            Shallow /@ Replace[#, HoldForm[f_[___]] :> HoldForm[f], 1],
            Style[Map[Short, Last[#], {2}], Red]
           } &@Drop[Drop[Stack[_], -7], 4]
         ];
         Message[args];
         Throw[$Failed, Message];
       ] /; ! TrueQ[inMessage];
    Protect[Message];
   ];
   Catch[StackComplete[code], Message]]


(* metric spaces: a symmetric matrix of distances *)
(* here are some routines which procedurally generate some common spaces *)
SphereSpace[dim_ /; dim==0] := {{0,2},{2,0}}
SphereSpace[dim_ /; dim>0] :=
	ArrayFlatten[{{SphereSpace[dim-1], Table[1, {i, 1, 2dim}, {j, 1, 2}]},
				  {Table[1, {i, 1, 2}, {j, 1, 2dim}], SphereSpace[0]}}]

TorusPoints[radius_, freq1_, thickness_, freq2_] :=
	Join@@Table[
		{Cos[u](radius+thickness Cos[v]),
		 Sin[u](radius+thickness Cos[v]),
		 thickness Sin[v]}//N,
			{u, 2Pi/freq1, 2Pi, 2Pi/freq1},
			{v, 2Pi/freq2, 2Pi, 2Pi/freq2}]
TorusSpace[radius_, freq1_, thickness_, freq2_] :=
	With[
		{points = TorusPoints[radius, freq1, thickness, freq2]},
		Table[With[
				{vector = points[[i]]-points[[j]]},
				Sqrt[Plus@@(vector*vector)]],
			{i, 1, Length[points]}, {j, 1, Length[points]}]]

WedgeSpace[m1_, m2_] :=
	With[
		{corner = m2[[1, 2 ;; Length[m2]]] + #&/@m1[[1, All]]},
		BlockSum[m1, corner,
				 Transpose[corner], m2[[2 ;; Length[m2], 2 ;; Length[m2]]]]]


(* graph: {number of vertices, list of edges} *)
(* these routines build graphs, in various presentations, from the data of a metric space *)
FormGraph[metricspace_,threshold_]:=
	Module[{i, j, edges={}},
		For[i=1, i<Length[metricspace], i++,
			Module[{workingedges={}},
				For[j=i+1, j<=Length[metricspace], j++,
					If[metricspace[[i,j]]<=threshold,
						workingedges=Join[workingedges,{{i,j}}]]];
				edges = Join[edges, workingedges]]];
		{Length[metricspace], edges}]

FormAdjacencyMatrix[metricspace_, threshold_] :=
	Map[If[# <= threshold, 1, 0]&, metricspace, {2}] - IdentityMatrix[Length[metricspace]]


(* cochain complex: list of coboundary maps, each described as a matrix *)
(* simplicial complex: list of lists of simplices, each described as an ordered list of vertices *)
(* returns: (simplicial complex, cochain complex, cochain ranks) *)
FormComplexData[adjacency_] :=
	If[Length[adjacency] == 1, {{{{1}}},{{}}},
	With[{recursion = FormComplexData[TopMinor[adjacency]],
		  newEdges=Last[adjacency]},
		Module[{dim,
				parentIndices = {},
				newSimplices = {List /@ Range[1, Length[adjacency]]},
				newDifferentials = {},
				oldBigSimplices, oldSmallSimplices, oldBigDiffl, oldSmallDiffl,
				j, A, B, C, D,
				workingSimplices, workingIndices,
				simplex},
			For[dim = 1, dim <= Length[recursion[[1]]], dim++,
				oldBigSimplices = If[dim != Length[recursion[[1]]],recursion[[1, dim+1]],{}];
				oldSmallSimplices = recursion[[1, dim]];
				oldBigDiffl = If[dim < Length[recursion[[1]]], recursion[[2, dim]], {}];
				oldSmallDiffl =
					If[dim == 1,
						SparseArray[_ -> -1, {Length[oldSmallSimplices], 1}],
						recursion[[2, dim-1]]];
				workingSimplices = {};
				workingIndices = {};
				For[j = 1, j <= Length[oldSmallSimplices], j++,
					simplex = oldSmallSimplices[[j]];
					If[And @@ (newEdges[[#]] == 1 & /@ simplex),
						AppendTo[workingSimplices, Join[simplex, {Length[adjacency]}]];
						AppendTo[workingIndices,j]]];
				A = oldBigDiffl;
				B = If[Length[oldBigSimplices] == 0 || (dim > 1 && Length[parentIndices[[dim-1]]] == 0),
						{},
						SparseArray[_ -> 0, {Length[A], If[dim == 1, 1, Length[parentIndices[[dim-1]]]]}]];
				C = If[Length[oldSmallSimplices] == 0 || Length[workingSimplices] == 0,
						{},
						With[{wI = workingIndices} (* we need a local copy to persist with the array *),
							SparseArray[{i_,j_} :> -Boole[wI[[i]] == j],
								{Length[wI], Length[oldSmallSimplices]}]]];
				D = If[Length[workingSimplices] == 0 || (dim > 1 && Length[newSimplices[[dim-1]]] == 0),
						{},
						-1*oldSmallDiffl[[workingIndices, If[dim != 1, parentIndices[[dim-1]], {1}]]]];
				AppendTo[newDifferentials, BlockSum[A, B, C, D]];
				AppendTo[newSimplices, Join[oldBigSimplices, workingSimplices]];
				AppendTo[parentIndices, workingIndices];
			];
		If[Mod[Length[adjacency], 10] == 0, Print[Length[adjacency]]];
		If[Last[newSimplices] == {},
			{Drop[newSimplices, -1],
			 Drop[newDifferentials, -1],
			 Length /@ Drop[newSimplices, -1]},
			{newSimplices,
			 newDifferentials,
			 Length /@ newSimplices}]]]]
FormSimplicialComplex[adjacency_] := FormComplexData[adjacency][[1]]
FormCochainComplex[adjacency_] := FormComplexData[adjacency][[2 ;; 3]]


FormCocycles[complex_,ranks_] := 
	MapThread[
		Which[
			#1=={} && #2==0, 0,
			#1=={}, IdentityMatrix[#2],
			True, With[{ns=NullSpace[#1]}, If[ns=={}, {}, Transpose[ns]]]]&,
		{Append[complex,{}],ranks}]

FormCoboundaries[complex_,ranks_] :=
	If[#=={}, {}, MatrixImage[#]]& /@ Prepend[complex,{}]

BettiNumbers[z_,b_,ranks_] :=
	MapThread[
		Which[
			#1===0, #3,
			#1==={}, 0,
			True, Dimensions[#1][[2]]] - If[#2!={}, Dimensions[#2][[2]], 0]&,
		{z,b,ranks}]


(* it's useful to have actual cocycle representations of cohomology classes, rather than just the Betti numbers.  these routines build such representatives. *)
ComplementSubspace[m_] :=
	Module[{work = m, index = 1},
		While[MatrixRank[work] < Length[m],
		With[{nextwork = ArrayFlatten[{{work, SparseArray[{i_, _} /; i == index :> 1, {Length[m], 1}]}}]},
			If[MatrixRank[nextwork] > MatrixRank[work], work = nextwork];
			index = index+1]];
		If[work == m, {},
			work[[1 ;; Length[m], Dimensions[m][[2]]+1 ;; Dimensions[work][[2]]]]]]
PullbackSpan[cospanLeft_, cospanRight_] :=
	If[cospanLeft === {} || cospanLeft === 0 || cospanRight === {} || cospanRight === 0, {{},{}},
		With[{kernel = Transpose[NullSpace[ArrayFlatten[{{cospanLeft, -cospanRight}}]]]},
			{kernel[[1 ;; Dimensions[cospanLeft][[2]], All]],
			 kernel[[Dimensions[cospanLeft][[2]]+1 ;; Length[kernel], All]]}]]
HomologyRepresentatives[z_, b_, dim_] :=
	With[{cocycles = z[[dim+1]], coboundaries = b[[dim+1]]},
	With[{pullback = PullbackSpan[cocycles, MatrixImage[coboundaries]]},
		If[pullback === {{}, {}}, MatrixImage[cocycles],
			With[{pz = pullback[[1]], pb = pullback[[2]]},
			With[{step = ComplementSubspace[pz.Inverse[pb]]},
				If[step === {} || step === 0, {}, cocycles.step]]]]]]


(* build the cochain cup product of two cochains belonging to a simplicial complex *)
CupCochains[simplicialcomplex_, c1_, dim1_, c2_, dim2_] :=
	If[dim1+dim2+1 > Length[simplicialcomplex], {},
		Transpose[{Plus @@ With[
			{cochains1 = simplicialcomplex[[dim1+1]], 
			 cochains2 = simplicialcomplex[[dim2+1]],
			 cochainsT = simplicialcomplex[[dim1+dim2+1]]},
			With[{cochain1 = #}, Plus @@ (With[{cochain2 = #},
				If[Last[cochain1] != First[cochain2], 0,
					With[{posT = Position[cochainsT, Join[cochain1, Rest[cochain2]]]},
						If[posT=={}, 0,
						With[{w1 = c1[[Position[cochains1, cochain1][[1, 1]], 1]],
							  w2 = c2[[Position[cochains2, cochain2][[1, 1]], 1]]},
							w1 * w2 * DeltaColumn[posT[[1, 1]], 1, Length[cochainsT]]]]]]]&
				/@ cochains2)]& /@ cochains1]}]]


(* routine to convert any cocycle into a linear combination of representative cohomology classes. *)
(* returns a list of coefficients to scale up the homology reps *)
ConvertCocycleToStdRep[cocycle_, coboundaries_, homologyreps_] :=
	If[cocycle === 0 || cocycle === {} || homologyreps === 0 || homologyreps === {}, 0,
		With[{dim1 = Dimensions[homologyreps], dim2 = Dimensions[coboundaries]},
			PullbackSpan[
				ArrayFlatten[{{homologyreps, coboundaries}}],
				cocycle][[1]][[1 ;; Dimensions[homologyreps][[2]], 1]]]]


AddSymbToMonomial[symb_ (*triple*), monomial_ (*list of triples*)] :=
	With[{existingposition = Position[monomial, {symb[[1]], symb[[2]], _}]},
		Which[
			existingposition === {},
			With[{ret = SortBy[Prepend[monomial, symb], Drop[#, -1]&]},
				{ret,
				 (-1)^(symb[[1]]*Plus@@(Times@@Drop[#, {2}]& /@
					 ret[[1 ;; Position[ret, symb][[1, 1]]-1, All]]))}],
			EvenQ[symb[[1]]],
			Module[{ret = monomial},
				AddTo[ret[[existingposition[[1, 1]]]], {0, 0, symb[[3]]}];
				{ret, 1}],
			True,
				{0, 1}]]


(* DGAMonomial is a Head for lists of triples (dimension, ID, power). these lists are assumed to be dictionary-ordered ascending. *)
DGAMonomial[{}] := 1;
Protect[DGAMonomial];
(* elements in a cDGA graded-commute with one another.  this behavior is not supported by the multiplication symbol * so we introduce a new binary symbol DGAMult[]. *)
(* distributivity*)
DGAMult[x_ + y_, z_] := DGAMult[x, z] + DGAMult[y, z];
DGAMult[x_, y_ + z_] := DGAMult[x, y] + DGAMult[x, z];
(* linearity against scaling factors *)
DGAMult[k_Integer x_, y_] := k DGAMult[x, y];
DGAMult[x_, k_Integer y_] := k DGAMult[x, y];
(* special multiplicative values *)
DGAMult[x_, 1] := x;
DGAMult[1, x_] := x;
DGAMult[x_, 0] := 0;
DGAMult[0, x_] := 0;
(* monomial collection *)
DGAMult[x_DGAMonomial, y_DGAMonomial] :=
	With[{list1 = Last[x], list2 = Last[y]},
	With[{dirtywork = AddSymbToMonomial[Last[list1],list2]},
		DGAMult[
			DGAMonomial[Drop[list1, -1]],
			dirtywork[[2]] * If[dirtywork[[1]] === 0, 0, DGAMonomial[dirtywork[[1]]]]]]];


(* sullivan dictionary: list of DGAMonomials, arranged by total dimension and then by dictionary order. *)
(* generally, a sullivan model will be presented as a sequence of classes, DGASymb[dim, id], together with formulas for their differentials.  to compare this with the cochain complexes we were computing before, we have to coordinatize these symbolic formulas into matrices.  a sullivan dictionary gives instructions for forming such a coordinatization, by giving an ordered listing of the monomials involved. *)
AddToSullivanDictionary[dict_, dim_, id_] :=
	Module[{working = dict, n},
		AppendTo[working[[dim]], DGAMonomial[{{dim, id, 1}}]];
		For[n = dim + 1, n <= Length[dict], n++,
			(* multiply everything through by the new symbol *)
			working[[n]] =
				Join[working[[n]],
				Select[
					DGAMonomial[AddSymbToMonomial[{dim, id, 1}, Last[#]][[1]]]& /@ working[[n-dim]],
					# =!= DGAMonomial[0]&]]];
		working]


(* a little rewrite grammar for computing differentials of product classes.  this is actually one of the cooler features of Mathematica.  NOTE--i basically reused the stuff from the Steenrod rewrite project from 2009. *)

DGADegree[DGAMonomial[list_]] := Plus @@ (Times @@ Drop[#, {2}]& /@ list);
DGADegree[x_ + y_] := Max[DGADegree[x], DGADegree[y]];
DGADegree[x_^n_] := n * DGADegree[x];
DGADegree[DGAMult[x_, y_]] := DGADegree[x] + DGADegree[y];

Protect[DGADiff];
LinearityRule = DGADiff[k_Integer x_] -> k DGADiff[x];
SumRule = DGADiff[x_ + y_] -> DGADiff[x] + DGADiff[y];
LeibnizRule = DGADiff[DGAMonomial[list_List]] /; Length[list] > 1 :> With[{newpiece = First[list]}, DGAMult[DGADiff[DGAMonomial[{newpiece}]], DGAMonomial[Rest[list]]] + (-1)^(newpiece[[1]]*newpiece[[3]])DGAMult[DGAMonomial[{newpiece}],DGADiff[DGAMonomial[Rest[list]]]]];
PowerRule = DGADiff[DGAMonomial[{{dim_, id_, power_}}]] /; power > 1 :> power DGAMult[DGAMonomial[{{dim,id,power-1}}],DGADiff[DGAMonomial[{{dim,id,1}}]]];
BasicRules = {LinearityRule, SumRule, LeibnizRule, PowerRule};
ApplyDifferential[rules_, class_] :=
	With[{step = Expand[class /. Join[rules, BasicRules]]},
		If[step===class, step, ApplyDifferential[rules, step]]]


(* sullivan model: list of list of rules indicating the action of the differential *)
FormDictionaryFromSullivanModel[diffs_, bound_] :=
	Fold[AddToSullivanDictionary[#1, #2[[1]], #2[[2]]]&,
		Table[{}, {i, 1, bound + 1}],
		Join @@ MapIndexed[(#2)&, diffs, {2}]]
FormCochainComplexFromSullivanModel[diffs_, bound_] :=
	Module[{ret = Table[{}, {i, 1, bound}], n},
	With[{dictionary = FormDictionaryFromSullivanModel[diffs, bound]},
		For[n = 1, n <= bound, n++,
			With[{results = ApplyDifferential[Join @@ diffs, DGADiff[#]]& /@ dictionary[[n]]},
			ret[[n]] = Table[Coefficient[results[[j]], dictionary[[n+1,i]]],
							{i, 1, Length[dictionary[[n + 1]]]},
							{j, 1, Length[results]}];
			If[ret[[n]] != {} && ret[[n,1]] == {}, ret[[n]] = {}]]];
		Prepend[ret,
			If[Length[dictionary[[1]]] == 0, {},
			   Table[0, {i, 1, Length[dictionary[[1]]]}, {j, 1, 1}]]]]]


(* forms a sullivan algebra presentation of the cohomology of a simplicial complex!
   s_ is really the only input; c_, z_, and b_ are readily computed using routines
   above, but as they're often already computed it seems prudent to pass them in
   and save the effort. *)

(* decompose a product monomial into two factors. guarantees deg(fact1) \[LessEqual] deg(fact2). *)
DecomposeClass[DGAMonomial[list_]] := 
	If[Length[list]==1,
		With[{tuple = list[[1]]},
			{DGAMonomial[{{tuple[[1]], tuple[[2]], 1}}],
				DGAMonomial[{tuple-{0, 0, 1}}]}],
			{DGAMonomial[{Join[First[list][[1 ;; 2]],{1}]}],
				DGAMonomial[If[list[[1, 3]] == 1, Rest[list],
					list - SparseArray[{1, 3} -> 1, {Length[list], 3}]]]}]

FormSullivanAlgebra[s_, c_, z_, b_, hin_, bound_] :=
	Module[{
	h = Join[hin, Table[{}, {i, Length[h]+1, bound+1}]],
	n = 1,
	diffs = Table[{}, {i, 1, bound+1}],
	horizontalmap = Table[{}, {i, 1, bound+1}],
	cohomologymap = Table[{}, {i, 1, bound+1}],
	DGAH = Table[{}, {i, 1, bound+1}]},
	For[n=1, n <= bound, n++,
		Print["n: " <> ToString[n]];
		Module[{dict, DGAc, DGAz, DGAb, j},
		dict = FormDictionaryFromSullivanModel[diffs, bound];
		DGAc = FormCochainComplexFromSullivanModel[diffs, bound];
		DGAz = FormCocycles[DGAc, Prepend[Length /@ dict, 1]];
		DGAb = FormCoboundaries[DGAc, Prepend[Length /@ dict, 1]];
		(* compute the homology of the cDGA in dimension n *)
		DGAH[[n]] = HomologyRepresentatives[DGAz, DGAb, n];

		(* compute the action of the map C^n DGA --> C^n Space by multiplicativity
		   NOTE: every class in degree n at this moment is decomposable. *)
		(* TODO: this iterated BlockSum is not very efficient. *)
		Print["computing chain map"];
		Module[{decomp, degrees, col1, col2},
		For[j=1, j<=Length[dict[[n]]], j++,
			decomp = DecomposeClass[dict[[n,j]]];
			degrees = DGADegree /@ decomp;
			col1 = Position[dict[[degrees[[1]]]], decomp[[1]]][[1,1]];
			col2 = Position[dict[[degrees[[2]]]], decomp[[2]]][[1,1]];
			horizontalmap[[n]] =
				BlockSum[
					horizontalmap[[n]],
					CupCochains[s,
								horizontalmap[[degrees[[1]]]][[All,col1;;col1]],
								degrees[[1]],
								horizontalmap[[degrees[[2]]]][[All,col2;;col2]],
								degrees[[2]]],
					{},
					{}]
		]];

		(* compute the action of the map H^n DGA --> H^n Space by restriction.
		   recall that CCTSR actually outputs such a matrix, not one valued in
		   chains! *)
		Print["computing cohomology map"]
		If[DGAH[[n]]=!={} && DGAH[[n]]=!=0 && h[[n]]=!={},
			With[{dimhmap = Dimensions/@horizontalmap, dimdgah = Dimensions/@DGAH},
			Module[{HtoCmap =
				If[horizontalmap[[n]]=!={} && horizontalmap[[n]]=!=0,
					horizontalmap[[n]].DGAH[[n]], 0]},
				For[j=1, j<=Dimensions[DGAH[[n]]][[2]], j++,
					cohomologymap[[n]] = BlockSum[cohomologymap[[n]],
						Transpose[{ConvertCocycleToStdRep[HtoCmap[[All,j;;j]],
														 b[[n+1]],h[[n]]]}],
						{},{}]
		]]]];

		(* KERNEL FUSSING *)
		(* select a basis for the kernel, and select DGA representatives for them *)
		Print["correcting kernel"];
		If[DGAH[[n]] =!= {},
			With[{kernel = If[cohomologymap[[n]]=!={} && cohomologymap[[n]]=!=0,
					NullSpace[cohomologymap[[n]]],
					IdentityMatrix[Dimensions[DGAH[[n]]][[2]]]]},
				If[kernel=!=0 && kernel=!={} && kernel=!={0},
					For[j=1, j<=Dimensions[kernel][[2]], j++,
						(* introduce elements in dimension (n-1) with
						   differential onto these reps *)
						(* XXX: it seems likely that kernel[[j]] needs to be
						   rotated or something. *)
						AppendTo[diffs[[n-1]], 
							DGADiff[DGAMonomial[{{n-1,Length[diffs[[n-1]]]+1,1}}]]->DGAH[[n]].kernel[[j]].dict[[n]]];
						(* let the horizontal map act by zero. *)
						horizontalmap[[n-1]] = 
							BlockSum[horizontalmap[[n-1]],
								Table[0, {i, 1, If[n>Length[s], 0, Length[s[[n]]]]},
										 {j, 1, Dimensions[kernel][[2]]}], {}, {}]
		]]]];

		(* COKERNEL FUSSING *)
		(* select a basis for the cokernel, and find space cochain representatives *)
		Print["correcting cokernel"];
		If[h[[n]] =!= {},
			With[{cokernel =
				If[cohomologymap[[n]] =!= {},
					With[{image = MatrixImage[cohomologymap[[n]]]},
						If[image=!={},
							ComplementSubspace[image],
							IdentityMatrix[Length[cohomologymap[[n]]]]]],
					IdentityMatrix[Dimensions[h[[n]]][[2]]]]},
				If[cokernel=!=0 && cokernel=!={} && cokernel=!={0},
					For[j=1, j<=Dimensions[cokernel][[2]], j++,
						(* introduce elements in dimension n with vanishing
						   differential *)
					    AppendTo[diffs[[n]],
								DGADiff[DGAMonomial[{{n,Length[diffs[[n]]]+1,1}}]]->0]];
					(* let the horizontal map act on these by landing on the
					   cokernel reps *)
					horizontalmap[[n]] =
						BlockSum[horizontalmap[[n]], h[[n]].cokernel, {}, {}];
					cohomologymap[[n]] =
						BlockSum[cohomologymap[[n]], cokernel, {}, {}]]]]
	]]; (* end of For loop on n *)
	(* return the stuff we built! *)
	{diffs,horizontalmap}]


InsertionPoints[smallList_, bigList_] :=
	Module[{i, j, ret = {}},
	For[i = 1; j = 1, i<=Length[smallList], j++,
		If[smallList[[i]] == bigList[[j]],
			AppendTo[ret, j];
			i++]];
	ret]

InsertionMatrix[smallList_, bigList_] :=
	With[{points = InsertionPoints[smallList,bigList]},
		SparseArray[{i_, j_} :> Boole[points[[j]] == i],
				    Length /@ {bigList, smallList}]]

ThresholdCochainMap[s1_, s2_] :=
	With[{pad1 = Join[s1, Table[{}, {i, 1, Length[s2] - Length[s1]}]],
		  pad2 = Join[s2, Table[{}, {i, 1, Length[s1] - Length[s2]}]]},
		MapThread[If[NullMatrix[#1] || NullMatrix[#2], {}, InsertionMatrix[#1, #2]]&,
				  {pad1,pad2}]]


(* we're going to build maps of Sullivan models in the next section and we'll need some way to store them.  this seems most sanely done in a manner wholly similar to how we store DGADiff[]s. *)
Protect[DGAMap];
MapLinearityRule = DGAMap[k_Integer x_] -> k DGAMap[x];
MapSumRule       = DGAMap[x_+y_] -> DGAMap[x]+DGAMap[y];
MapProductRule   = DGAMap[DGAMonomial[list_List]] /; Length[list]>1 :> DGAMult[DGAMap[{First[list]}], DGAMap[Rest[list]]];
MapPowerRule     = DGAMap[DGAMonomial[{{dim_,id_,power_}}]] /; power>1 :> DGAMult[DGAMap[DGAMonomial[{{dim,id,power-1}}]], DGAMap[DGAMonomial[{{dim,id,1}}]]];
BasicMapRules    = {MapLinearityRule, MapSumRule, MapProductRule, MapPowerRule};

ApplyDGAMap[rules_,class_] := 
	With[
		{step = Expand[class /. Join[rules, BasicMapRules]]},
		If[step===class, step, ApplyDGAMap[rules, step]]]


(* the construction of minimal sullivan models for cDGAs is not functorial as minimality throws it off. it can however be made functorial in the homotopy category. given minimal sullivan models for simplicial complexes C and D and a map f this routine computes a lift of f (up to homotopy) between the Sullivan models. *)
(* XXX: it's conceivable that these things in the return array "map" should be separated into dimensions. *)
LiftSullivanMap[diffs1_, m1_, c1_, diffs2_, m2_, c2_, bound_, f_] :=
	Module[{n, map={}, j},
	With[{dict1 = FormDictionaryFromSullivanModel[diffs1, bound],
		  dict2 = FormDictionaryFromSullivanModel[diffs2, bound],
          sc2 = FormCochainComplexFromSullivanModel[diffs2,bound]},
		(* iterate over the active dimensions *)
		For[n = 1, n <= bound, n++,
			(* in this dimension, indecomposable classes either
			   have zero or nonzero internal boundary. *)
			For[j = 1, j <= Length[diffs1[[n]]], j++,
				With[{termname = (diffs1[[n,j]]) /. ((DGADiff[term_] -> _) :> term),
					  target = (diffs1[[n,j]]) /. ((_ -> target_) :> target)},
				If[MatchQ[diffs1[[n,j]], _ -> 0],
					(* for those with zero boundary: *)
					(* apply f.m1 to them... *)
					With[{classInM2 = f[[n+1]].({#}& /@ m1[[n, All, Position[dict1[[n]], termname][[1, 1]]]])},
					(* lift that chain into cDGA2 by taking m2^-1 up to boundaries. *)
					With[{alpha =
						PullbackSpan[
							BlockSum[m2[[n]],
									 If[n > 1, c2[[n]], 0],
									 0,
									 0],
							classInM2][[1]]},
					AppendTo[map,
						DGAMap[termname] -> Plus @@ (alpha[[1 ;; Length[dict2[[n]]], 1]]*dict2[[n]])]]],
				(* } else { *)
					With[{s2class = ApplyDGAMap[map, DGAMap[target]]},
					(* this class must be a boundary in cDGA2. find a preimage. *)
					With[{alpha = PullbackSpan[sc2[[n]], {Coefficient[s2class, #]}& /@ dict2[[n+1]]][[1]]},
					AppendTo[map, DGAMap[termname] -> Plus @@ (alpha[[All, 1]]*dict2[[n]])]]]
		]]]];
		map]]


(*EndPackage[];*)


(*

(* DEMO: compute the inclusion of a lower-threshold sampled torus into a higher-threshold sampled torus. :) *)
radius=1;
freq1=16;
thickness=1/3;
freq2=8;
threshold=30/48; (* 7/12 is better *)
m=TorusSpace[radius,freq1,thickness,freq2];
Length[m];
FormGraph[m,threshold];
ListPointPlot3D[TorusPoints[radius,freq1,thickness,freq2]]
g=FormAdjacencyMatrix[m,threshold];
AdjacencyGraph[g,GraphLayout\[Rule]"SpringEmbedding"]
d=FormComplexData[g];
s=d[[1]];
c=d[[2]];
z=FormCocycles[c,d[[3]]];
Print["done with cocycles"]
b=FormCoboundaries[c,d[[3]]];
Print["done with coboundaries"]
BettiNumbers[z,b,d[[3]]]
h=HomologyRepresentatives[z,b,#]&/@Range[1,Length[s]-1];
Print["done with homology reps"]
a=FormSullivanAlgebra[s,c,z,b,h,4];
%[[1]]

*)
