CechOpsStable.wl

(* ::Package:: *)

BeginPackage["CechOpsStable`"];


(* ::Title:: *)
(*Descriptions of Public (Unhidden) Functions*)


(* ::Section:: *)
(*Cech Nerve Operations and Indexing*)


partitionToCover::usage="partitionToCover[part] Takes in a collection of disjoint subsets part={lambda_1,...lambda_n}, where lambda_1 are sets themselves
denoted by lists, (a partition of X = Union[{lambda_1,...lambda_n}]) and outputs the associated 'complementary' cover {X-lambda_1,...X-lambda_n}.";

complementaryCover::usage="complementaryCover[n] takes in an integer n and outputs the cover of the set {1,...,n} given by taking complements of each element.";

numPrim::usage="Takes in a cover {U_1, ... U_n} and outputs the length n: the number of 'primitive' elements of the cover.";

indexToSimplex::usage="indexToSimplex[simpnum,deg,cover] takes in an index I=simpnum and outputs the I^th simplex of the cover (according to lexicographical ordering).";

simplexToIndex::usage="simplexToIndex[simp,cover] takes in a simplex simp and outputs its associated index with respect to cover (using lexicographical ordering).";

simpFace::usage="simpFace[simp,k] takes the kth face of simp.  Here k must be an element of {0,...,size(simp)-1}.";

support::usage="support[simp,cover] outputs the support of the simplex simp.  E.g. if cover = {{1,2},{2,3}} and simp = {1,2}, then support[simp,cover] = {2}.";

numSimp::usage="numSimp[deg,cover] outputs the number of simplicies of degree deg with respect to the cover.";

simpList::usage="simpList[deg,cover] outputs a list of the simplices of degree 'deg' with respect to the cover 'cover'.";

simpSuppList::usage="simpSuppList[deg,cover] outputs a list of the supports of simplices of degree deg with respect to the cover.";


(* ::Section:: *)
(*Face and Alternating Boundary Operator*)


face::usage="face[chain,k,deg,funMor,cover] is the kth face map on the simplicial object formed by the Cech complex 
generated by the functor {funObj,funMor} and the cover 'cover', it is evaluated on a chain 'chain' of degree 'deg-1'; the output is a chain 
(thought of as a function on simplices with respect to the cover 'cover'): given a simplex 'simpeval' of degree
deg-1, face[chain,k,deg,funMor,cover][simpeval] is an element of funObj(support(simpeval)).";


bdry::usage="bdry[chain,deg,funMor,cover] is the boundary operator, evaluated on the chain 'chain' of degree 'deg',
 on the alternating face-map complex associated to the Cech complex generated by the functor {funObj,funMor} and the cover 'cover'.
Its output is a chain of degree deg-1.";


(* ::Section::Closed:: *)
(*Natural Transformation to Free Module valued functor*)


funObjVects::usage="funObjVects[sourceobj,funObj,inprod] can be thought of in two ways: the Cartan matrix associated
to an inner product on the linear span of the list funObj[sourceobj], or as a morphism from the linear span of funObj[sourceob]
into a standard vector space equipped with the standard inner product (each element of funObjVects is the image of an element of
funObj in the standard vector space).";

eltToVect::usage="eltToVect[elt,sourceobj,funObj,inprod] takes the element 'elt' of funObj[sourceob] to its image
in the 'standard vector space equipped with the standard inner product' using the inner product 'inprod'.  Inverse to vectToElt";

vectToElt::usage="vectToElt[vect,sourceobj,funObj,inprod] takes a vector in the standard vector space with standard inner product
 to an element of funObj[sourceobj] equipped with inner product 'inprod'.  Inverse to eltToVect.";

morToMat::usage="morToMat[src,tgt,mor,funObj,inprod] takes a morphism: funObj[src] -> funObj[tgt] to a matrix after identifying
funObj with a module/vector space valued function using 'inprod'.";

funMorMat::usage="funMorMat[funObj_,funMor_,inprod_] is application of morToMat to funMor: funMorMat[funObj_,funMor_,inprod_][src,tgt]
outputs the matrix associated to funMor[src,tgt].";


(* ::Section:: *)
(*Face and Alternating Boundary Matrices*)


chainToVect::usage="chainToVect[chain,deg,funObj,cover,inprod] takes the 'chain' (an element of the associated Cech complex given the data 
'funObj' and 'cover') to its image in a 'standard vector space equipped with the standard inner product' using the inner product 'inprod'.
  Inverse to vectToChain.";

vectToChain::usage="";

morOfChainToMat::usage="";

rkChain::usage="";

faceMatSlow::usage="faceMatSlow[deg,k,funMor,funObj,cover,inprod] computes the matrix form of the 
'k'th face map at degree 'deg' by applying the function morOfChainToMat to the function 'face'.  As the name suggests, this function is slower than the more
direct (but less elegantly implemented) function 'faceMat': it uses the function morOfChainToMat.";

bdryMat::usage="";

faceMat::usage="faceMat[deg,k,funMor,funObj,cover,inprod] computes the matrix form of the 
'k'th face map at degree 'deg'.  For a more elegant (but slower) implementation of this function
see the function 'faceMatSlow'.";

faceMatSparse::usage="";


(* ::Section:: *)
(*Semi-Simplicial Alternating Complex Homology*)


homologyRank::usage="homologyRank[deg, funObj, funMor, cover, inprod]";


homologyVect::usage="homologyVect[deg,funObj,funMor,cover,inprod]";


homologyObj::usage="homologyObj[deg,funObj,funMor,cover,inprod]";


(* ::Title:: *)
(*Function Definitions*)


Begin["`Private`"];


(* ::Section:: *)
(*Miscellaneous Operations*)


(* ::Subsection:: *)
(*Linear Algebraic Operations*)


stdBasis[dim_]:=If[dim==0, {{0}} , IdentityMatrix[dim]];


zeroMat[n_,m_]:=ConstantArray[0,{n,m}];


blockMatrix[indices_,values_,dims_]:=ArrayFlatten@ReleaseHold@Normal@SparseArray[indices->Map[Hold,values],dims];


inprodToMat[inprod_,basis_]:=Outer[inprod,basis,basis,1];


genAdjoint[M_,inprodsrc_,inprodtgt_]:= Module[{ipMatSrc, ipMatTgt},
ipMatSrc=inprodToMat[inprodsrc,IdentityMatrix[Dimensions[M][[1]]]];
ipMatTgt=inprodToMat[inprodtgt,IdentityMatrix[Dimensions[M][[2]]]];
PseudoInverse[ipMatSrc] . Transpose[M] . ipMatTgt
];


orthoComp[vspace_,ip_]:=Module[{ipMat=inprodToMat[ip,IdentityMatrix[Dimensions[vspace][[2]]]]}, NullSpace[vspace . ipMat]];


rowSpace[M_]:=RowReduce[M][[1;;MatrixRank[M]]];


(*These work over the reals -- to work over the integers use HermiteDecomposition**)


imageSpace[M_]:=rowSpace@Transpose[M];


kernelSpace[M_]:=If[NullSpace[M]==={}, {ConstantArray[0,Dimensions[M][[2]]]} , NullSpace[M]];


vectSpInt[vsp1_,vsp2_]:=Module[{ker,coeffs,intBasis},
ker=NullSpace@Transpose@Join[vsp1,vsp2];
coeffs=ker[[All,1;;Length@vsp1]];
 intBasis=rowSpace@Map[# . vsp1&,coeffs];
If[intBasis==={}||coeffs==={},{ConstantArray[0,Length@vsp1[[1]]]},intBasis]
];


vectSpIntMult[vspaceList_]:=Fold[vectSpInt,vspaceList];


(* ::Subsection:: *)
(*Partition List via list of lengths (see http://forums.wolfram.com/mathgroup/archive/1992/Nov/msg00096.html)*)


splitFast[list_,parts_]:=Block[{accumulation=FoldList[Plus,0,parts]},Inner[Take[list,{#1,#2}]&,Drop[accumulation,-1]+1,Rest[accumulation],List]];


(* ::Section:: *)
(*Cech Nerve Operations and Indexing*)


partitionToCover[part_]:=If[
 Length@part==1, (*trivial partition becomes the trivial cover*)
 part,
 Map[Complement[Union@@part,#]&,part]
];


complementaryCover[n_]:=partitionToCover[Table[{i},{i,1,n}]];


numPrim[cover_]:=Length[cover];


numSimp[deg_,cover_]:=Binomial[numPrim[cover],deg+1];


simpList[deg_,cover_]:=Subsets[Range@numPrim[cover],{deg+1}];


indexToSimplex[simpnum_,deg_,cover_]:=simpList[deg,cover][[simpnum]];


simplexToIndex[simp_,cover_]:=If[simp==={},1,FirstPosition[simpList[Length[simp]-1,cover],simp][[1]]];


simpFace[simp_,k_]:=Delete[simp,k+1];


support[simp_,cover_]:=If[simp==={}, Union@@cover, Intersection@@cover[[simp]] ];


simpSuppList[deg_,cover_]:= Map[support[#,cover]&, simpList[deg,cover]];


simpMapToList[map_,deg_,cover_]:=Map[map@support[#,cover]&, simpList[deg,cover]];


(* ::Section:: *)
(*Face and Alternating Boundary Operator*)


face[chain_,k_,deg_,funMor_,cover_][simpeval_] := If[deg<0, 0, 
(*else*) Module[{simpFacePreIm, chainExtend},
simpFacePreIm=Flatten@Position[Map[simpFace[#,k]&,simpList[deg, cover]],simpeval];
chainExtend=funMor[support[#,cover],support[simpFace[#,k],cover]][chain[#]]&;
Total@Map[chainExtend@indexToSimplex[#,deg,cover]&,simpFacePreIm]
]
];


bdry[chain_,deg_,funMor_,cover_]:=If[deg<0, 0&,
(*else*) Sum[(-1)^(k) face[chain,k,deg,funMor,cover][#],{k,0,deg}]&
];


(* ::Section::Closed:: *)
(*Natural Transformation to Free Module valued functor*)


funObjVects[sourceobj_,funObj_,inprod_]:=Outer[inprod[sourceobj],funObj[sourceobj],funObj[sourceobj],1];


vectInprod[sourceobj_,funObj_,inprod_]:=#1 . funObjVects[sourceobj,funObj,inprod] . #2&;


eltToVect[elt_,sourceobj_,funObj_,inprod_]:=If[elt===0,ConstantArray[0,Length@funObj[sourceobj]], 
Map[inprod[sourceobj][#,elt]&,funObj[sourceobj]]
];


vectToElt[vect_,sourceobj_,funObj_,inprod_]:=(PseudoInverse[funObjVects[sourceobj,funObj,inprod]] . vect) . funObj[sourceobj];


morToMat[src_,tgt_,mor_,funObj_,inprod_]:=Module[{morVect,dimVectSrc},
morVect=Composition[eltToVect[#,tgt,funObj,inprod]&,mor,vectToElt[#,src,funObj,inprod]&];
dimVectSrc=Length[funObj[src]];
Transpose@Map[morVect, stdBasis[dimVectSrc]]
];

funMorMat[funObj_,funMor_,inprod_]:=morToMat[#1,#2,funMor[#1,#2],funObj,inprod]&;


(* ::Section:: *)
(*Face Maps as Matrices*)


(* ::Subsection:: *)
(*Ranks of Simplicial Objects/Alternating Complex Components*)


rkChain[deg_,funObj_,cover_]:=If[deg<-1,0,Total@simpMapToList[Length@*funObj,deg,cover]];


(* ::Subsection:: *)
(*Elegant Implementation of Face Map Matrix*)


chainToVect[chain_,deg_,funObj_,cover_,inprod_]:= Join@@Map[eltToVect[chain[#],support[#,cover],funObj,inprod]&, simpList[deg,cover]];


vectToChain[vect_,deg_,funObj_,cover_,inprod_]:=Module[{dimsOfFactors,vectPartitioned},
dimsOfFactors=simpMapToList[Length@*funObj,deg,cover];
vectPartitioned=splitFast[vect,dimsOfFactors];
vectToElt@@{Part[vectPartitioned, simplexToIndex[#,cover]],support[#,cover],funObj,inprod}&
];


chainVectInprod[deg_,funObj_,cover_,inprod_]:=Module[{Factors,srcObjFactors,dimsOfFactors,vectPartitioned,inprodShuffled},
srcObjFactors=simpSuppList[deg,cover];
Factors=simpMapToList[funObj,deg,cover];
dimsOfFactors=Map[Length,Factors];
vectPartitioned=splitFast[#,dimsOfFactors]&;
inprodShuffled=vectInprod[srcObjFactors[[First@#2]],funObj,inprod][First@#1,Last@#1]&;
Total@MapIndexed[inprodShuffled,Transpose[{vectPartitioned[#1],vectPartitioned[#2]}],1]&
];


morOfChainToMat[srcdeg_,tgtdeg_,mor_,funObj_,cover_,inprod_]:=
Piecewise[{{zeroMat[1,rkChain[-1,funObj,cover]],srcdeg==-1&&tgtdeg<-1},
{zeroMat[rkChain[-1,funObj,cover]],srcdeg<-1&&tgtdeg==-1},
{zeroMat[1,1],srcdeg<-1&&tgtdeg<-1}},
(*else*)
Module[{chainMorVect,dimVectSrc},
chainMorVect=Composition[chainToVect[#,tgtdeg,funObj,cover,inprod]&,mor,vectToChain[#,srcdeg,funObj,cover,inprod]&];
dimVectSrc=rkChain[srcdeg,funObj,cover];
Transpose@Map[chainMorVect, stdBasis[dimVectSrc]]
]
];


faceMatSlow[deg_,k_,funMor_,funObj_,cover_,inprod_]:=morOfChainToMat[deg,deg-1,face[#,k,deg,funMor,cover]&,funObj,cover,inprod];


(* ::Subsection::Closed:: *)
(*Faster Implementation of Face Map Matrix*)


faceMat[deg_,k_,funMor_,funObj_,cover_,inprod_]:=
Piecewise[{{zeroMat[1,rkChain[-1,funObj,cover]],deg==-1&&k==0},
{zeroMat[1,1],deg<0&&k==0}},
(*else*)
Module[{simpsSrc,simpsTgt,simpsSrcInds,simpsTgtInds,simpsSrcSupp,simpsSrcFaces,
simpsSrcFacesSupps,simpToIndLocal,simpsSrcFacesInds,nonZeroIndices,simpsTgtNoFaceInds,
zeroIndices,zeroSimplices,zeroPart,matSparseIndices,matSparseValues},
    simpsSrc=simpList[deg,cover];
    simpsTgt=simpList[deg-1,cover];
    simpsSrcInds=Range@numSimp[deg,cover];
    simpsTgtInds=Range@numSimp[deg-1,cover];    
    simpsSrcSupp=simpSuppList[deg,cover];
    simpsSrcFaces=Map[simpFace[#,k]&, simpsSrc];
    simpsSrcFacesSupps=Map[support[#,cover]&,simpsSrcFaces];
    simpToIndLocal[simpyMcsimp_]:= If[simpyMcsimp==={}, 1, FirstPosition[simpsTgt,simpyMcsimp][[1]]];
    simpsSrcFacesInds=Map[simpToIndLocal,simpsSrcFaces];
    nonZeroIndices=Transpose@{simpsSrcFacesInds,simpsSrcInds};
    zeroIndices=Complement[Flatten[Outer[List,simpsTgtInds,simpsSrcInds],1], nonZeroIndices];
    zeroSimplices=Map[{indexToSimplex[First[#],deg-1,cover],indexToSimplex[Last[#],deg,cover]}&,zeroIndices];
    matSparseIndices=Join[nonZeroIndices,zeroIndices];
    zeroPart= Map[ConstantArray[0,#]&,Map[Length@funObj@support[#,cover]&,zeroSimplices,{2}]];
    matSparseValues=Join[Map[Apply@funMorMat[funObj,funMor,inprod],Transpose@{simpsSrcSupp,simpsSrcFacesSupps}],zeroPart];
    blockMatrix[matSparseIndices,matSparseValues,{numSimp[deg-1,cover],numSimp[deg,cover]}]
]
];


faceMatSparse[deg_,k_,funMor_,funObj_,cover_,inprod_]:=
Piecewise[{{zeroMat[1,rkChain[-1,funObj,cover]],deg==-1&&k==0},
{zeroMat[1,1],deg<0&&k==0}},
(*else*)
Module[{simpsSrc,simpsTgt,simpsSrcInds,simpsTgtInds,simpsSrcSupp,simpsSrcFaces,
simpsSrcFacesSupps,simpToIndLocal,simpsSrcFacesInds,nonZeroIndices,simpsTgtNoFaceInds,
zeroIndices,zeroSimplices,zeroPart,matSparseIndices,matSparseValues},
    simpsSrc=simpList[deg,cover];
    simpsTgt=simpList[deg-1,cover];
    simpsSrcInds=Range@numSimp[deg,cover];
    simpsTgtInds=Range@numSimp[deg-1,cover];    
    simpsSrcSupp=simpSuppList[deg,cover];
    simpsSrcFaces=Map[simpFace[#,k]&, simpsSrc];
    simpsSrcFacesSupps=Map[support[#,cover]&,simpsSrcFaces];
    simpToIndLocal[simpyMcsimp_]:= If[simpyMcsimp==={}, 1, FirstPosition[simpsTgt,simpyMcsimp][[1]]];
    simpsSrcFacesInds=Map[simpToIndLocal,simpsSrcFaces];
    nonZeroIndices=Transpose@{simpsSrcFacesInds,simpsSrcInds};
    zeroIndices=Complement[Flatten[Outer[List,simpsTgtInds,simpsSrcInds],1], nonZeroIndices];
    zeroSimplices=Map[{indexToSimplex[First[#],deg-1,cover],indexToSimplex[Last[#],deg,cover]}&,zeroIndices];
    matSparseIndices=Join[nonZeroIndices,zeroIndices];
    zeroPart= Map[ConstantArray[0,#]&,Map[Length@funObj@support[#,cover]&,zeroSimplices,{2}]];
    matSparseValues=Join[Map[Apply@funMorMat[funObj,funMor,inprod],Transpose@{simpsSrcSupp,simpsSrcFacesSupps}],zeroPart];
   SparseArray`SparseBlockMatrix[Thread[matSparseIndices->matSparseValues],{rkChain[deg-1,funObj,cover],rkChain[deg,funObj,cover]}]
]
];


(* ::Section::Closed:: *)
(*Alternating Complex*)


altComplexVect[deg_,funObj_,funMor_,cover_,inprod_]:=Module[{dimVectSrc},
If[deg<-1, {{0}},
(*else*)
dimVectSrc=Total@simpMapToList[Length@*funObj,deg,cover];
stdBasis[dimVectSrc]
]
];


altComplexOb[deg_,funObj_,funMor_,cover_,inprod_]:=Map[vectToChain[#,deg,funObj,cover,inprod]&,altComplexVect[deg,funObj,funMor,cover,inprod]];


(* ::Subsection:: *)
(*Alternating Complex Boundary Matrix*)


bdryMat[deg_,funMor_,funObj_,cover_,inprod_]:=If[deg<=-1,faceMat[deg,0,funMor,funObj,cover,inprod],
Sum[(-1)^(k) faceMat[deg,k,funMor,funObj,cover,inprod],{k,0,deg}] 
];


(* ::Section::Closed:: *)
(*Semi-Simplicial Alternating Complex Homology*)


dimKer::usage="dimKer[M] calculates the dimension of the kernel of a Matrix M using the rank-nullity theorem; should equal Length[NullSpace[M]].";

dimKer[M_]:=Dimensions[M][[2]]-MatrixRank[M];


homologyRank[deg_, funObj_, funMor_, cover_, inprod_] := Module[{bdryprev, bdry},
   bdry = bdryMat[deg, funMor, funObj, cover, inprod];
   bdryprev = bdryMat[deg + 1, funMor, funObj, cover, inprod];
   dimKer[bdry] - MatrixRank[bdryprev] ];


homologyVect[deg_,funObj_,funMor_,cover_,inprod_]:=Module[{bdryprev,bdry,rowSpOrthoComp},
   bdry = bdryMat[deg, funMor, funObj, cover, inprod];
   bdryprev = bdryMat[deg + 1, funMor, funObj, cover, inprod];
   rowSpOrthoComp=orthoComp[imageSpace[bdryprev],chainVectInprod[deg,funObj,cover,inprod]];
   vectSpInt[kernelSpace[bdry],rowSpOrthoComp] 
];


homologyObj[deg_,funObj_,funMor_,cover_,inprod_]:=Map[vectToChain[#,deg,funObj,cover,inprod]&, homologyVect[deg,funObj,funMor,cover,inprod]];


(* ::Title:: *)
(*End Matter*)


End[];


EndPackage[]