\datethis % print date on listing

\def\cite#1{[{\bf #1}]}

\centerline{\bf Generating 3D Discrete Morse Functions from Point Data}
\medskip
\centerline{\bf Henry C. King}
\bigskip

This code is a companion to the paper \cite{KKM} by King, Knudson, and Mramor
which explains the theory behind why it works. 
It takes a simplicial complex $K$,
which is a subcomplex of a compact 3 dimensional manifold\footnote{$^1$}
{Actually $M$ does not have to be a manifold as long as $M$
minus its vertices is a manifold and the link of each vertex is connected.
But probably in practice $M$ will be $S^3$ anyway.} $M$,
and an integer valued function $h$ on the vertices of $K$
and extends it to a discrete Morse function on $K$ in
the sense of Forman \cite{F}.
It also cancels pairs of critical $j$ and $j-1$ simplices
if they are joined by a single gradient path 
and the maximum values of $h$ on the two simplices
differ by less than a given persistence $p$.


We do not actually calculate a discrete Morse function on $K$,
but instead determine the crucial information given by a discrete Morse function.
This crucial information is the designation of some simplices as critical,
and a pairing of the remaining simplices so that each noncritical simplex is either
paired to one of its codimension one faces or is paired to a simplex
of which it is a codimension one face.

While \cite{KKM} assumes that $h$ is injective, this program does not.
The program gets around this assumption by, in effect, perturbing $h$.

We have tested this program with $M$ a triangulation of $B\times S^1$ where
$B$ is a fixed triangulation of the Klein bottle and we vary the triangulation
of the circle $S^1$ and take a product triangulation. Then we let $K$ be $M$ minus a disc
and put a random function on the vertices of $K$.
On a 933 MHz PowerPC G4 with 640 MB running Mac OS X with 11,960,000 simplices
it runs in about 30 seconds with persistence set at 5000 out of a total
range of 65536 for the value of $h$ on the vertices.
A average run took 19.44 seconds to generate a discrete Morse function,
.57 seconds to cancel 0 and 1 simplices, 1.71 seconds to cancel 2 and 3 simplices,
and 2.92 seconds to cancel 1 and 2 simplices.
An average run with half the number of simplices, 5,980,000 and the same persistence of 5000
took 9.21 seconds to generate a discrete Morse function,
.25 seconds to cancel 0 and 1 simplices, .46 seconds to cancel 2 and 3 simplices,
and 1.42 seconds to cancel 1 and 2 simplices.
All the above times are with the slower\_and\_finer option, which takes a bit longer 
to determine critical simplices around a given edge, but does it in such a way that
the critical triangles descend as steeply as possible from the edge.
For 11,960,000 simplices without the slower\_and\_finer option it's about .12 seconds
faster to generate a discrete Morse function but canceling 1 and 2 simplices
is about .39 seconds slower and canceling 2 and 3 simplices is .09 seconds slower.
For what it's worth, in the end about 4\% of the vertices were critical,
about 2.8\% of the edges were critical, about 1.7\% of the triangles were critical, and
about .7\% of the tetrahedra were critical.


In practice, the code looks linear in the number of simplices of $M$.
The major portion of the time is spent generating a discrete Morse function.
This is linear in the number of simplices if there is an a priori bound on the number
of simplices in the Star of a vertex.  Without such a bound I expect the time to 
be proportional to the number of vertices times the average square of the number of 
simplices in the Star of a vertex.
This square comes from local cancelation of critical simplices in the Star of a vertex,
since presumably the number of critical simplices in the Star of a vertex is on average
proportional to the number of all simplices in the Star of a vertex and the average path length
is likewise proportional.
I presume that cancelling critical points would tend to be quadratic in the number of
simplices since I would expect the number of cancelled critical points to be
proportional to the number of simplices and for each pair of cancelled critical
points the gradient path between them presumably has length proportional to the
number of simplices.  In experiments this is not evident except possibly when
cancelling 2 and 3 simplices.  I only say ``possibly" because there was a weird phenomenon
where the time to cancel 2 and 3 simplices for the same 11,960,000 simplex complex
was either about 1 second or about 1.8 seconds depending on when it was run.
Always the lower figure right after compilation and maybe a few more runs,
then the higher figure.  Go figure.




\bigskip
\centerline{Bibliography}
\medskip

\item{\cite{F}} R.~Forman, {\it A User's Guide to Discrete Morse Theory}\/.

\smallskip
\item{\cite{KKM}} H.~King, K.~Knudson, and N.~Mramor, 
{\it Generating Discrete Morse Functions from Point Data}\/, preprint.


@* Extract. 
Portability issue: This code requires |sizeof(int) = sizeof(void *)|.
This could be fixed by duplicating the ordered list routines so that there is a version using
pointer keys rather than integer keys, or else always converting any address keys
(which are always pointers to simplexes) to the index of the simplex.

This documentation is produced by cweave which uses the following symbols:
\item{}|NULL| is a NULL pointer
\item{}|&&| is logical and (\AM\AM)
\item{} |||| is logical or ($\mid\mid$)
\item{} |!| is logical not (!)

The program has the following structure:

@c
@<Header files to include@>@/
@<Subroutine prototypes@>@/
@<Global variables@>@/
@<List functions@>@/
@<Simplex functions@>@/
@<Debugging output@>@/
@<Subroutines for constructing complexes@>@/
@<Subroutines finding gradient paths@>@/
@<Subroutines canceling a single pair of critical simplices@>@/
@<Subroutines canceling many pairs of critical simplices @>@/
@<The main routine |Extract|@>@/

@ @<Global variables@>=
list *crit[4];  // lists of critical simplices of dimensions 0,1,2, and 3
long num_simp[4];  // number of simplices
vertex *vertexlist;
edge *elist;
triangle *flist;
tetrahedron *tlist;


@ And now the main routine |Extract|. 
The parameter |vfirst| is the first vertex of the ambient manifold $M$ containing $K$,
and |p| is the persistence.
The complex $M$ will have information in it which allows us to determine $K$ 
and the function $h$ on the vertices of $K$.
We cancel 1 and 2 simplices last because this cancelation takes longer,
so it is best to reduce the number of 1 and 2 simplices as much as possible beforehand.

@<The main...@>=
void Extract(int p)
{@/
	long i;
	long t1,t2,t3,t4,t5,t6;
	
	t1 = clock();
	for(i=0;i<4;i++) list_initialize(crit+i,sizeof(void *));@/

	@<Find a nonoptimized discrete Morse function on $K$@>;
	t2 = clock();


	ExtractCancel1(p); //Cancel 0 and 1 simplices with persistence less than |p|
	t3 = clock();


	ExtractCancel3(p); //Cancel 2 and 3 simplices with persistence less than |p|
	t4 = clock();


	ExtractCancel2(p); //Cancel 1 and 2 simplices with persistence less than |p|
	t5 = clock();
	
#ifdef verbose
	
	printf(
	"\nExtract Raw = %d, Cancel1 = %d, Cancel2 = %d, Cancel3 = %d\n",  
	t2-t1,t3-t2,t5-t4,t4-t3);
#endif
}



@ We called this section of code |ExtractRaw| in \cite{KKM}.
For each vertex |v| of $K$,
it finds an optimal discrete Morse function on the lower Star of |v|.
(The lower Star of |v| is all simplices for which |v|
is the vertex of maximal value.)
It then takes one critical edge |beste| in the lower Star of |v|,
makes it noncritical, and pairs it with |v|.
It cancels all pairs of critical simplices it can
in the lower Star so that in the end you get as few of them
as possible.  There are other algorithms which find optimal
discrete Morse functions and are presumably faster if the
lower Star has a great many simplices, but this is simple
and works fine for reasonably sized links.
It also chooses critical simplices based on the given 
function values on the vertices
(especially if |slower_and_finer| is set), which other
algorithms we know of do not.
Note that if the lower Star of |v| is empty, you make the 
vertex |v| critical, since |v| is a local minimum.

When we say ``steepest" below we are assuming all edges
have equal length.
If we wish to account for edges of unequal length
and pick steepest edges, triangles, etc., it would be necessary
to modify the algorithm.


@<Find a nonoptimized discrete Morse function on $K$@>=
{
	edge *beste; // this will be the steepest edge from |v|
	vertex *v;
	int is_lower_Star_empty; 
	long vid;
	list *lcrit2;  // list of critical 2 simplices in lower Star of |v|
	list *edges_todo;  // list of unprocessed edges containing |v|

	list_initialize(&lcrit2,sizeof(triangle *)); 
	list_initialize(&edges_todo,sizeof(edge *));
	
	for(vid=0; vid < num_simp[0]; vid++)
	{
		v = id2vertex(vid);
		if(!is_in_K(v)) continue;
		beste = NULL;  

		@<extract the lower Star of |v|@>;
		if(is_lower_Star_empty) make_critical(v);
		else
		{
			pair01(v,beste);  // pair |v| with the steepest edge down
			LocalCancel(v,lcrit2);  // cancel critical simplices in the lower Star of |v|
		}
	}
	list_abandon(&lcrit2);
	list_abandon(&edges_todo);
}


@ Now we determine the lower Star of |v| and 
find a discrete Morse function on it.
We do this by taking each edge |e| in the lower Star of |v|
and putting a discrete Morse function on its lower Star.
The lower Star of |e| is all simplices $\sigma$ of $K$ containing 
|e| so that the maximal and next--to--maximal vertices
of $\sigma$ are in |e|.
@d mark_edge_done(v,e) ((e)->type |= (((v)==get_vertex(e,0))? 0x200 : 0x400))
@d edge_marked_done(v,e) ((e)->type & (((v)==get_vertex(e,0))? 0x200 : 0x400))


@<extract the lower Star of |v|@>=
{
	edge *e;
	vertex *w;
	int val,bestval,flag; 
		
	val = value(v);
	
	list_clear(edges_todo);

	plist_push(edges_todo,v->links[0]);  // push an edge containing |v|
	
	is_lower_Star_empty=1;    // flag to test if lower Star(v) is empty
	while(!list_is_empty(edges_todo))
	{
		@<pop the next item |e| from |edges_todo| which we haven't processed yet@>@;
		@<set |flag| and reset |is_lower_Star_empty| if |e| is in lower Star(v)@>@;
		
		@<add adjacent edges of |e| to |edges_todo| and Extract lower Star(e) if |flag| is set@>@;
	}
}


@ @<pop the next item |e| from |edges_todo| which we haven't processed yet@>=
e = plist_pop(edges_todo);  
if(edge_marked_done(v,e)) continue;  
mark_edge_done(v,e);


@ @<set |flag| and reset |is_lower_Star_empty| if |e| is in lower Star(v)@>=
flag = 0; // default
if(is_in_K(e))
{
	w = smallest_vertex_in_edge(e);
	
	if(w!=v) // is |e| in lower Star(v)?
	{
		set_value(e,val);  // fill in max value on e
		is_lower_Star_empty=0;   // indicate that lower Star(v) is not empty
		flag = 1;
	}
}


@ This section of code puts a discrete Morse function on the lower Star of the edge |e|.
Since this essentially means putting a discrete Morse function on a
subcomplex of a circle, we can use a quick algorithm to do this.
The drawback is that the choice of critical simplices has little to
do with the values of the given function on the vertices, beyond
the fact that this is used to determine the lower Star of |e|.
A slightly slower algorithm is used if you define |slower_and_finer|.
It chooses critical triangles 
so that the value at the vertex opposite |e| is a minimum in its connected
component of lower Star(e).
It takes a bit longer here, but experiments indicate it tends to speed up
canceling 1 and 2 simplices, presumably by making descending discs descend faster.

 The |state| variable has the following significance:
\item{0} means a triangle in lower Star(e) has not yet been encountered.
\item{1} means all triangles and tetrahedra so far are in lower Star(e)
\item{2} means the last triangle was not in lower Star(e), but some previous triangle was.
\item{3} means the last triangle was in lower Star(e), 
			but some previous triangles or tetrahedra weren't.
			
@d slower_and_finer 1

@<add adjacent edges of |e| to |edges_todo| and Extract lower Star(e) if |flag| is set@>=
{
	triangle *s,*bests,*sp;
	tetrahedron *t;
	edge *ep;
	int state=0;
	int bestsval;
	int first_in_lower_Star;
	vertex *min_vertex;  //  the vertex of |s| with minimum value

#ifdef slower_and_finer	
	triangle *bests_in_run,*first_bests_in_run; 
	tetrahedron *bestt_in_run,*first_bestt_in_run;
	int bestsval_in_run,first_bestsval;
#else
	vertex *first_min_vertex; //  the vertex of |sp| with minimum value
#endif	

	s=sp=e->links[2];  // a triangle containing |e|
	t = coface(s,0);  // a tetrahedron containing |s|
	
	do
	{
		ep=other_edge(v,e,s);
		if(!edge_marked_done(v,ep)) plist_push(edges_todo,ep);
#ifdef slower_and_finer	
		if(flag) @<xextract s and t@>@;
#else
		if(flag) @<extract s and t@>@;
#endif	
		t = other_coface(s,t);
		s = other_face(e,s,t);  
	} while (s!=sp);
#ifdef slower_and_finer	
	if(flag) @<xlast extract of s and t@>@;
#else
	if(flag) @<last extract of s and t@>@;
#endif	
}

@ @<extract s and t@>=
{
	int in_lower_Star_of_e;
	
	@<see if |s| is in lower Star(e) and set |min_vertex| and |in_lower_Star_of_e| accordingly@>;
	
	switch(state)
	{
		case 0:
			if(in_lower_Star_of_e) 
			{
				first_in_lower_Star = (s==sp && is_in_K(t));
				
				if(first_in_lower_Star) 
				{
					state = 1; 
					bests = NULL;
					first_min_vertex = min_vertex;
				}
				else 
				{
					state = 3;
					bests = s;
					bestsval = value(min_vertex);
				}
			}
			break;
		case 1:
		case 3:
			if(!in_lower_Star_of_e) 
			{
				state = 2;
				break;
			} 
			else if(is_in_K(t))
			{
				set_value(t,val);  // max value on |t| must be at vertex |v|
				pair23(s,t);
				break;
			}
			// else pass through to case 2
		case 2:
			if(in_lower_Star_of_e) 
			{
				@<make s critical xor make |bests| critical and replace |bests|@>@;
				state = 3;
			}
			break;
	}
}

@ @<see if |s| is in lower Star(e) and set |min_vertex| and |in_lower_Star_of_e| accordingly@>=
in_lower_Star_of_e=0;  // default
if(is_in_K(s))
{
	vertex *vlist[3];

	get_triangle_vertices(s,vlist);
	min_vertex = vlist[smallest_vertex(vlist,3)];
	in_lower_Star_of_e = (vertex_in_edge(min_vertex,e)<0);
	if(in_lower_Star_of_e) set_value(s,val);
}

@ @<make s critical xor make |bests| critical and replace |bests|@>=
if(bests==NULL)
{
	bests=s;
	bestsval = value(min_vertex);
}
else if(value(min_vertex)<bestval)
{
	make_critical(bests);
	plist_push(lcrit2,bests);
	bests=s;
	bestsval = value(min_vertex);
}
else
{
	make_critical(s);
	plist_push(lcrit2,s); // remember all critical triangles for |LocalCancel|
}


@ @<last extract of s and t@>=
switch(state)
{
	case 0:
		@<make |e| critical unless it is the best candidate so far to pair with |v|@>@;
		break;
	case 1:@/
		make_critical(t);
		set_value(t,val);
		pair12(e,s);
		break;
	case 3:
		if(first_in_lower_Star)
		{
			pair23(s,t);
			set_value(t,val);
		}
		pair12(e,bests);
		break;
	case 2:
		if(first_in_lower_Star)
		{
			min_vertex = first_min_vertex;
			@<make s critical xor make |bests| critical and replace |bests|@>@;
		}
		pair12(e,bests);
		break;
}


		


@ @<make |e| critical unless it is the best candidate so far to pair with |v|@>=
if(beste==NULL) // is it the first?
{
	beste = e;
	bestval = value(w);
}
else if(value(w)<bestval)
{
	make_critical(beste);
	beste = e;
	bestval = value(w);
}
else make_critical(e);


@ @<xextract s and t@>=
{
	int in_lower_Star_of_e;
	
	@<see if |s| is in lower Star(e) and set |min_vertex| and |in_lower_Star_of_e| accordingly@>;
	
	switch(state)
	{
		case 0:
			if(in_lower_Star_of_e) 
			{
				first_in_lower_Star = (s==sp && is_in_K(t));
				
				if(first_in_lower_Star) state = 1; 
				else	state = 3;

				bests = NULL;
				bests_in_run = s;
				bestt_in_run = t;
				bestsval_in_run = value(min_vertex);
			}
			break;
		case 1:
		case 3:
			if(!in_lower_Star_of_e || !is_in_K(t)) 
			{
				@<make |bests_in_run| or |bests| critical and update |bests|@>@;
			}
			if(!in_lower_Star_of_e) 
			{
				state = 2;
				break;
			} 
			else if(is_in_K(t))
			{
				set_value(t,val);  // max value on |t| must be at vertex |v|
				
				if(value(min_vertex) < bestsval_in_run)
					@<replace |bests_in_run| by |s|@>@;
				else pair23(s,t);
				break;
			}
			// else pass through to case 2
		case 2:
			if(in_lower_Star_of_e) 
			{
				bests_in_run = s;
				bestt_in_run = t;
				bestsval_in_run = value(min_vertex);
				state = 3;
			}
			break;
	}
}

@ @<make |bests_in_run| or |bests| critical and update |bests|@>=
if(bests==NULL)
{
	if(state==1) 
	{
		first_bests_in_run = bests_in_run;
		first_bestt_in_run = bestt_in_run;
		first_bestsval = bestsval_in_run;
	}
	else
	{
		bests = bests_in_run;
		bestsval = bestsval_in_run;
	}
}
else if(bestsval_in_run < bestval)
{
	make_critical(bests);
	plist_push(lcrit2,bests);
	bests = bests_in_run;
	bestsval = bestsval_in_run;
}
else
{
	make_critical(bests_in_run);
	plist_push(lcrit2,bests_in_run); // remember all critical triangles for |LocalCancel|
}

@ @<replace |bests_in_run| by |s|@>=
{
	triangle *ss;
	tetrahedron *tt;
	
	do 
	{
		tt = other_coface(bests_in_run,bestt_in_run);
		ss = r32(tt);
		pair23(bests_in_run,tt);
		bests_in_run = ss;
		bestt_in_run = tt;
	} while(bestt_in_run != t);
	bests_in_run = s;
	bestsval_in_run = value(min_vertex);
}

@ @<xlast extract of s and t@>=
{
	triangle *sss;
	
	switch(state)
	{
		case 0:
			@<make |e| critical unless it is the best candidate so far to pair with |v|@>@;
			break;
		case 1:@/
			make_critical(t);
			set_value(t,val);
			pair12(e,bests_in_run);
			break;
		case 3:
			if(first_in_lower_Star)
			{
				set_value(t,val);
				if(first_bestsval > bestsval_in_run)
				{
					sss = bests_in_run;
					bestt_in_run = other_coface(first_bests_in_run,first_bestt_in_run);
					bests_in_run = first_bests_in_run;
				}
				else sss = first_bests_in_run;
				@<replace |bests_in_run| by |s|@>@;
				bests_in_run = sss;
			}
			@<make |bests_in_run| or |bests| critical and update |bests|@>@;
			pair12(e,bests);
			break;
		case 2:
			if(first_in_lower_Star)
			{
				bests_in_run = first_bests_in_run;
				bestsval_in_run = first_bestsval;
				@<make |bests_in_run| or |bests| critical and update |bests|@>@;
			}
			pair12(e,bests);
			break;
	}
}





@ A few helpful functions.

@d min(x,y) (((x)>(y))?y:x)
@d max(x,y) (((x)>(y))?x:y)
@d in_MorseExtract 1

@<Header files to include@>=
#include <stdio.h>
#include "MorseExtract.h"

@ @<Subroutine prototypes@>=
void abort_message(char *s);



@* Cancel Routines.  
Sometimes a critical $m-1$ simplex and a critical $m$ simplex can be cancelled.
You can do this if there is exactly one gradient path between them.
A gradient path is a sequence of simplices 
$\sigma_i,\sigma_{i+1},\sigma_{i+2},\ldots ,\sigma_k$ so that
for $j$ even, $\sigma_j$ is an $m$ simplex,
and $\sigma_{j-1}$ is paired with $\sigma_j$, and $\sigma_{j+1}$
is a codimension one face of $\sigma_j$.
To cancel a pair of critical simplices $\sigma$ and $\tau$, 
you take the unique gradient path
$\sigma = \sigma_0,\ldots ,\sigma_{2k-1}=\tau$,
and you pair each $\sigma_{2i}$ with $\sigma_{2i+1}$
(rather than with $\sigma_{2i-1}$ as was done formerly).

@*2 Local Cancellation.
 |LocalCancel| cancels simplices in the lower Star of a vertex.
After doing it, you can show that you are left with a
minimal number of critical simplices.
For example, at a local maximum you would end up with only one
critical simplex, a 3 simplex.
Likewise at PL approximations of traditional smooth
Morse saddles of index $i$ you would end up with a single critical
simplex, an $i$ simplex.
@<Subroutines canceling many pairs of critical simplices @>=
void
LocalCancel(vertex *v,list *lcrit2)
{@/
	triangle *t,*s;
	tetrahedron *te[2];
	edge *e[2],*ep[2];
	int j,i,in;@/
	
	
	while (!list_is_empty(lcrit2))  // run through the list of new critical 2 simplices
	{
		t = plist_pop(lcrit2);
		@<find the edges |e[0]| and |e[1]| of |t| containing |v|@>;
		@<find the ends |ep[i]| of the gradient paths starting at |e[i]|@>;
				
		if(ep[0] != ep[1])  // can we cancel?
		{
			@<Cancel the |t| with the best of the two edges@>; 
		}
		else
		{
			// We can't cancel with an edge, so see if we can cancel with a 3 simplex
			for(i=0;i<2;i++)
			{
				@<find the 3 simplex |te[i]| connected to |coface(t,i)| by a gradient path@>;
			}
			if(te[0]!=te[1])
			{
				@<Cancel the |t| with the best of the two tetrahedra@>; 
			}
		}
	}
}


@ @<find the edges |e[0]| and |e[1]| of |t| containing |v|@>=

for(i=0,j=0;j<3;j++)
{
	if(vertex_in_edge(v,get_edge(t,j))>=0)
	{
		e[i++]=get_edge(t,j);
	}
}

@ @<find the ends |ep[i]| of the gradient paths starting at |e[i]|@>=
		
for(i=0;i<2;i++)	
{
	ep[i]=e[i];
	while(!is_critical(ep[i]))  
	{
		if(is_paired_down(ep[i]))
		{
			ep[i]=NULL;  // runs into $B_1$ so no cancelation possible
			break;
		}
		s = r12(ep[i]);  // get paired triangle
		ep[i] = other_edge(v,ep[i],s);
	}
}		

@ @<find the 3 simplex |te[i]| connected to |coface(t,i)| by a gradient path@>=
te[i] = coface(t,i);
s = t;
while (is_in_K(te[i]))  // while we are in K
{
	in = is_in_lower_Star(te[i],v);
	if(!in)  break;  // see if exit lower Star
	if(is_critical(te[i])) break; // see if reached critical
	s = r32(te[i]);
	te[i] = other_coface(s,te[i]);
}
if(!is_in_K(te[i]) || !in )		
	  te[i]=NULL;   

@ @<Cancel the |t| with the best of the two edges@>= 
if(ep[0]==NULL) LocalCancel12(v,ep[1],e[1],t);
else if(ep[1]==NULL) LocalCancel12(v,ep[0],e[0],t);
else if(min_value(ep[0],1) > min_value(ep[1],1)) 
	LocalCancel12(v,ep[0],e[0],t);
else LocalCancel12(v,ep[1],e[1],t);

@ @<Cancel the |t| with the best of the two tetrahedra@>= 
if (te[0]==NULL) Cancel23(t,1,te[1]);
else if (te[1]==NULL) Cancel23(t,0,te[0]);
else if(min_value(te[1],3) > min_value(te[0],3)) 
	Cancel23(t,0,te[0]);
else Cancel23(t,1,te[1]);  







@*2 Canceling vertices and edges. 
Find all critical edges and vertices connected by a single gradient path.
Cancel the pair with smallest difference in value, as long as it's less than |p|.
Keep on doing this until you can't do any more.


@<Subroutines canceling many pairs of critical simplices @>=

void ExtractCancel1(int p)
{@/ 
	olist *c0;  // list of vertices to cancel with
	olist *c1;  // list of edges which can be cancelled with vertices
	struct ccrit1
	{
		edge *s;
		long c[2];
	} s;
	int i;
	vertex *v[2];
	int thisk;

	olist_initialize(&c0,sizeof(long));
	olist_initialize(&c1,sizeof(struct ccrit1));

	@<Find all cancelable pairs and put them on lists |c0| and |c1|@>;
	
	while (!olist_is_empty(c1))
	{@t\1@>@/ 
		thisk = olist_min(c1,&s);  // put best candidate in |s|
		@<let |v[i]| be the two vertices with which |s| can cancel@>@;
		i = @<index of best vertex of |s.s| to cancel@>;
		if (is_critical(v[i]) && thisk == value(s.s)-value(v[i])) @<cancel with the vertex@>@;
		else @<see if we can still cancel with some other vertex and replace |s| on |c1|@>@;@t\2@>@/ 
	}@/
	olist_abandon(&c0);
	olist_abandon(&c1);@t\2@>@/
}

@ @<Find all cancelable pairs and put them on lists |c0| and |c1|@>=
{
	edge *e;

	list_read_init(crit[1]);
	while ((e=(edge *)plist_read(crit[1]))!=NULL)  // run through all critical edges
	{
		if(!is_critical(e)) list_read_delete(crit[1]);
		else
		{
			// find the two gradient paths descending from |e|
			v[0] = FindGrad01(get_vertex(e,0),value(e)-p);
			v[1] = FindGrad01(get_vertex(e,1),value(e)-p);
			if(v[0]!=v[1])
			{
				@<put |e| on the list |c1| and put |v[i]| on the list |c0|@>;
			}
		}
	}
}

@ @<put |e| on the list |c1| and put |v[i]| on the list |c0|@>=
{@t\1@>@/ 
	long link;

	link = -1;
	s.s=e;
	s.c[0] = (v[0]==NULL)? -1: olist_find_add(c0,(long)v[0],&link,NULL);
	s.c[1] = (v[1]==NULL)? -1: olist_find_add(c0,(long)v[1],&link,NULL);@/
	i = @<index of best vertex of |s.s| to cancel@>;
	olist_add(c1, value(e) - value(v[i]),&s);@t\2@>@/ 
}

@ @<let |v[i]| be the two vertices with which |s| can cancel@>=
v[0] = (s.c[0]<0)? NULL : (vertex *)olist_get_key(c0,s.c[0]);
v[1] = (s.c[1]<0)? NULL : (vertex *)olist_get_key(c0,s.c[1]);


@ @<cancel with the vertex@>=
		
{
	Cancel01(v[i],i,s.s);
	*((long *)olist_entry(c0,s.c[i])) = s.c[1-i];  
	// now any gradient path to |v[i]| goes to other vertex |v[1-i]|
}
		
@ @<see if we can still cancel with some other vertex and replace |s| on |c1|@>= 
{@t\1@>@/ 
	int m;
	
	for(i=0;i<2;i++)
	{
		if(s.c[i]<0) continue;
		@<update |s.c[i]| to point to a critical vertex@>@;
	}
	
	
	if( s.c[1]!=s.c[0])
	{@t\1@>@/
		i = @<index of best vertex of |s.s| to cancel@>;
		m = value(s.s) - value(v[i]);

		if( m < p)	olist_add(c1,m,&s);@t\2@>@/
	} @t\2@>@/
}

@ @<update |s.c[i]| to point to a critical vertex@>=
{
	long a;

	for(a = s.c[i]; a >= 0; )
	{
		v[i] = (vertex *)olist_get_key(c0,a);
		if( is_critical(v[i])) break;
		a = *((long *)olist_entry(c0,a));
	}
	*((long *)olist_entry(c0,s.c[i])) = a;
	if (value(v[i])<=value(s.s)-p)
	{
		v[i] = NULL;
		s.c[i] = -1;
	}
	else s.c[i] = a;
}

@ @<index of best vertex of |s.s| to cancel@>=
(v[0]==NULL)? 1:@t}$\par${@>
	(v[1]==NULL)? 0:@t}$\par${@>
		(value(v[0]) > value(v[1]))? 0: @t}$\par${@>
			(value(v[0]) < value(v[1]))? 1: @t}$\par${@>
				is_critical(v[0])?  1:0




@*2 Canceling triangles and tetrahedra.
 Find all critical triangles and tetrahedra connected by a single gradient path.
Cancel the pair with smallest difference in value, as long as it's less than |p|.
Keep on doing this until you can't do any more.
The code is essentially the same as that used to cancel edges and vertices.

@<Subroutines canceling many pairs of critical simplices @>=

void ExtractCancel3(int p)
{@/
	olist *c3;  // list of tetrahedra to cancel
	olist *c2;  // list of triangles which cancel with tetrahedra
	struct ccrit2
	{
		triangle *s;
		long c[2];
	} r;
	int i;
	tetrahedron *t[2];
	int thisk;

	olist_initialize(&c2,sizeof(struct ccrit2));
	olist_initialize(&c3,sizeof(long));

	@<Find all cancelable pairs and put them on lists |c2| and |c3|@>;
	while (!olist_is_empty(c2))
	{@t\1@>@/
		thisk = olist_min(c2,&r);
		@<let |t[i]| be the two tetrahedra with which we can cancel@>@;
		i = @<index of best coface of |r.s| to cancel@>;
		if (is_critical(t[i]) && thisk == value(t[i])-value(r.s))  @<cancel with the tetrahedron@>@;
		else @<see if we can still cancel with some other tetrahedron and replace on |c2|@>@;@t\2@>@/ 
	}@/
	olist_abandon(&c2);
	olist_abandon(&c3);@t\2@>@/
}

@ @<Find all cancelable pairs and put them on lists |c2| and |c3|@>=
{@t\1@>@/
	triangle *s;

	list_read_init(crit[2]);
	while ((s = (triangle *)plist_read(crit[2]))!=NULL) // go through all critical triangles
	{@/
		if(!is_critical(s)) list_read_delete(crit[2]);
		else
		{
			// find the beginnings of the two gradient paths ending at |s|
			t[0] = FindGrad23(coface(s,0),value(s)+p);
			t[1] = FindGrad23(coface(s,1),value(s)+p);@/
			if(t[0]!=t[1])
				@<put |s| on |c2| and |t[0]| and |t[1]| on |c3|@>@;
		}
	}@t\2@>@/
}

@ @<put |s| on |c2| and |t[0]| and |t[1]| on |c3|@>=
{@t\1@>@/
	long link;
	
	link = -1;
	r.s=s;
	r.c[0] = (t[0]==NULL)? -1: olist_find_add(c3,(long)t[0],&link,NULL);
	r.c[1] = (t[1]==NULL)? -1: olist_find_add(c3,(long)t[1],&link,NULL);@/
	
	i = @<index of best coface of |r.s| to cancel@>;
	olist_add(c2, value(t[i]) - value(s), &r);@t\2@>@/
}

@ @<let |t[i]| be the two tetrahedra with which we can cancel@>=
t[0] = (r.c[0]<0)? NULL : (tetrahedron *)olist_get_key(c3,r.c[0]);
t[1] = (r.c[1]<0)? NULL : (tetrahedron *)olist_get_key(c3,r.c[1]);
	
@ @<cancel with the tetrahedron@>=
{
	Cancel23(r.s,i,t[i]);
	*((long *)olist_entry(c3,r.c[i])) = r.c[1-i]; 
}

@ @<see if we can still cancel with some other tetrahedron and replace on |c2|@>= 		
{@t\1@>@/
	int m;
	
	for(i=0;i<2;i++)
	{
		if(r.c[i]<0) continue;
		@<update |r.c[i]| to point to a critical tetrahedron@>@;
	}
	
	
	if(r.c[1]!=r.c[0])
	{@t\1@>@/
		i = @<index of best coface of |r.s| to cancel@>;
		
		m = value(t[i]) - value(r.s);
		
		if( m < p)	olist_add(c2,m,&r);@t\2@>@/
	}@t\2@>@/
}

@ @<update |r.c[i]| to point to a critical tetrahedron@>=
{
	long a;

	for(a = r.c[i]; a >= 0; )
	{
		t[i] = (tetrahedron *)olist_get_key(c3,a);
		if( is_critical(t[i])) break;
		a = *((long *)olist_entry(c3,a));
	}
	*((long *)olist_entry(c3,r.c[i])) = a;
	if(value(t[i])>=value(r.s)+p)
	{
		t[i] = NULL;
		r.c[i] = -1;
	}
	else r.c[i] = a;
}
	
@ @<index of best coface of |r.s| to cancel@>=
(t[0]==NULL)? 1:@t}$\par${@>
	(t[1]==NULL)? 0:@t}$\par${@>
		(value(t[0]) > value(t[1]))? 1: @t}$\par${@>
			(value(t[0]) < value(t[1]))? 0: @t}$\par${@>
				is_critical(t[0])? 1: 0





@*2 Canceling edges and triangles.
 Find and cancel pairs of 1 and 2 simplices whose values differ by less than p.
This version speeds things up by not always canceling the pair with least persistence.
In particular, after one cancelation it may become possible for
some other pair of critical simplices to cancel which could not have been cancelled
before.  If this happens, this routine could be unaware until it has
cancelled all the pairs it knows about already.
Checking for this situation would slow down this routine dramatically.
In fact such occurrences appear to be quite rare.
For random functions on |K| with millions of simplices and persistence 1/13 of
the range of the values of vertices, it tends to happen at most once or twice.

@<Subroutines canceling many pairs of critical simplices @>=
void ExtractCancel2(int p)
{
	list *changed;  // triangles whose pairing was changed by Cancel12
	list *grad_path;  // edges in gradient path of a canceling pair
	olist *goodpairs; // cancelable critical triangles
	int lastp;

	list_initialize(&changed,sizeof(triangle *));
	list_initialize(&grad_path,sizeof(edge *));
	olist_initialize(&goodpairs,sizeof(triangle *));
	lastp = 0;
	do
	{
		@<fill list |goodpairs| with possible cancels@>@;
		if(olist_is_empty(goodpairs)) break;
		@<cancel all pairs on |goodpairs|@>@;
	} while(1);
	list_abandon(&changed);
	list_abandon(&grad_path);
	olist_abandon(&goodpairs);
}

@ @<fill list |goodpairs| with possible cancels@>=
{
	triangle *t;
	edge *e;

	list_read_init(crit[2]);
	while ((t = (triangle *)plist_read(crit[2]))!=NULL) // go through all critical triangles
	{
		if(!is_critical(t)) list_read_delete(crit[2]);
		else
		{
			e = FindGradPaths12(t,p,NULL,0);
			if(e!=NULL) // if there is a gradient path from |t| to |e|
				olist_add(goodpairs,value(t)-value(e),&t);
		}
	}
}




@ @<cancel all pairs on |goodpairs|@>=
{
	triangle *t;
	edge *e;
	int bp,thisp;

	while(!olist_is_empty(goodpairs))
	{
		bp = olist_min(goodpairs,&t);
		e = FindGradPaths12(t,p,grad_path,0);
		if(e==NULL)  continue;  // no gradient paths from |t| found
		thisp = value(t)-value(e);
		if(thisp<=bp) 
		{
			list_clear(changed);
			Cancel12(e,t,grad_path,changed);  // cancel |e| and |t|
			@<fix up split-rejoin paths@>;
#ifdef verbose
			if(thisp<lastp) printf("\n persistence out of order: %d > %d\n",lastp,thisp);
#endif
			lastp = thisp;
		}
		else olist_add(goodpairs,thisp,&t);
	}
}



	


@ After canceling a pair of one and two simplices,
there may be new pairs of gradient paths which differ only on the
boundary of a single tetrahedron.
This code modifies the simplex pairings to eliminate this happening.
It is not clear whether it is worth doing this.
Experiments show that this tends to allow cancellation of a few more pairs
of critical one and two simplices, but not many.
Experiments also show that very little time is spent in this code.
So I have left it in.
We have not yet attempted to prove that this fixing up process always
terminates, so to be careful I have put a limit of 10000 iterations which has so far never
been exceeded.

@<fix up split-rejoin paths@>=
{
	int i,j,k;
	tetrahedron *te;

	k = 0;
	while (!list_is_empty(changed) && k<10000)
	{
		t = plist_pop(changed);
		if(t==NULL) continue;
		for (i=0;i<2;i++)
		{
			te = coface(t,i);
			if( is_in_K(te) && !is_critical(te) )
			{
				j = splitrejoin(te);
				if(j>=0)
				{
					plist_push(changed,unsplitrejoin(te,j));
					k++;
				}
			}
		}
	}
	if (!list_is_empty(changed)) printf("May be infinite loop in split-rejoin");
}




@*1 Canceling pairs of Critical Simplices.  
The following routine cancels a critical 1 simplex $\sigma$ containing |v|
 with a critical 2 simplex $\tau$ containing |v|.
 It speeds up the ordinary Cancel12 by assuming that 
 the unique gradient path from $\tau$ to $\sigma$
 has all of its simplices containing |v|.
 The gradient path from $\tau$ to $\sigma$ goes through the face |sigmap| of $\tau$.

@<Subroutines canceling a single pair of critical simplices@>=
void LocalCancel12(vertex *v,edge *sigma,edge *sigmap,triangle *tau)
{
	edge *u;
	triangle *w,*wp;
	
	u=sigmap;
	w=tau;
	unmake_critical(tau);
	while(!is_critical(u)) // while u is noncritical
	{
		wp = r12(u);
		pair12(u,w);
		w=wp;
		u=other_edge(v,u,w);
	}
	unmake_critical(sigma);
	pair12(sigma,w);
}

@ The following routine cancels a critical 2 simplex $\sigma$ 
with a critical 3 simplex $\tau$. 
The n is such that  the gradient path from $\tau$ to $\sigma$
goes through |coface(sigma,n)|.
It does not delete $\sigma$ and $\tau$ 
from |crit[2]| and |crit[3]|.

@<Subroutines canceling a single pair of critical simplices@>=
void Cancel23(triangle *sigma,int n,tetrahedron *tau)
{
	triangle *s,*sp;
	tetrahedron *t;

	s=sigma;
	t=coface(s,n);
	unmake_critical(sigma);
	
	while(!is_critical(t))  // while t is not critical
	{
		sp=r32(t);
		pair23(s,t);
		s=sp;
		t=other_coface(s,t);
	}
	unmake_critical(tau);
	pair23(s,tau);
}

@ The following routine cancels a critical vertex |v| 
with a critical edge $\kappa$. It does not delete |v| and $\kappa$ 
from |crit[0]| and |crit[1]| however.
The n is such that  the gradient path from |v| to $\kappa$
goes through |get_vertex(kappa,n)|.

@<Subroutines canceling a single pair of critical simplices@>=


void Cancel01(vertex *v,int n,edge *kappa)
{
	vertex *u;
	edge *e,*ep;
	
	u = get_vertex(kappa,n);
	
	e = kappa;
	unmake_critical(kappa);
	
	while(!is_critical(u)) // while non critical
	{
		ep = r01(u);
		pair01(u,e);
		e = ep;
		u = other_vertex_in_edge(u,e);
	}
	unmake_critical(v);
	pair01(v,e);
}


@ Cancel a critical 1 simplex $\sigma$ with a critical 2 simplex $\kappa$,
using the list of edges in the gradient path |grad_path|.
Push any 2 simplices with changed pairings to the stack |changed|.
Note $\sigma$ and $\kappa$ are not deleted from |crit[1]| and |crit[2]|. 

@<Subroutines canceling a single pair of critical simplices@>=
void Cancel12(edge *sigma,triangle *kappa,list *grad_path,list *changed)
{
	triangle *t,*nextt;
	edge *e;

	unmake_critical(kappa);
	unmake_critical(sigma);

	t = kappa;
	
	do
	{
		e = plist_pop(grad_path);
		if(e==sigma) break;
		nextt = r12(e);
		pair12(e,t);
		plist_push(changed,t);
		t = nextt;
	} while(1);
	pair12(e,t);
	plist_push(changed,t);
}




@*1 Finding Gradient Paths. Find the end of a gradient path starting at $u\in K_0$.
Return the critical vertex $v$ at the end of the path,
or return |NULL| if $h(v)\le m$.

@<Subroutines finding gradient paths@>=
vertex *FindGrad01(vertex *u,int m)
{
	vertex *v;
	edge *e;

	v=u;
	while(!is_critical(v) && value(v) > m)
	{
		e = r01(v);
		v = other_vertex_in_edge(v,e);
	}
	if (value(v) <= m) return NULL;
	return v;
}

@ Find the start of a gradient path passing through a
tetrahedron $\tau$.
Return the critical 3 simplex |t| at the start of the path,
or return |NULL| if the value of |t| is  $\ge m$ or if the path starts
at a boundary 2 simplex.  
Mark simplices deadend if we know the only gradient paths going to them
start at a boundary 2 simplex.
This is a lazy way of shortening future gradient path searches.


@<Subroutines finding gradient paths@>=
tetrahedron *FindGrad23(tetrahedron *tau,int m)

{
	tetrahedron *t;
	triangle *s;
	
	if (tau==NULL) return NULL;  // obsolete?
	t=tau;
	while( is_in_K(t) && value(t) < m )
	{
		if (is_critical(t)) return t;  // reached critical t
		s = r32(t);
		if (is_deadend(s)) 
		{
			make_deadend(t);
			return NULL;
		}
		t = other_coface(s,t);
		if (is_deadend(t)) 
		{
			make_deadend(s);
			return NULL;
		}
	}
	if(!is_in_K(t) && t!=tau) make_deadend(s);
	return NULL;
}

@ Find all of the gradient paths starting at $\sigma$ which descend less than $p$,
and return the best one.
If |grad_path| is not |NULL| return a list of the edges in the best gradient path. 
The |flags| parameter, if odd, is experimental code which
changes the behavior by returning if possible a critical edge 
which is connected to $\sigma$  by at least one gradient path 
and could possibly generate persistent homology.

For performance reasons, the lists used are declared static
so they do not have to be initialized each time the routine is called.
I found that without doing this this routine spent
much of its time calling |malloc|.
As a consequence this routine is not reentrant.
If this is ever a problem, just delete the words |static| below.

Note that here and elsewhere we are essentially just searching for paths in a graph,
in this case for nodes connected by exactly one path.
No doubt there are sophisticated and well-known algorithms for doing this
which would improve performance.
As an excercise the reader can improve on these naive implementations.

@<Subroutines finding gradient paths@>=
edge *FindGradPaths12(triangle *sigma, int p, list *grad_path, int flags)
{
	static olist *graph;
	static list *to_do;
	static list *crits;
	static int first = 0;
	struct edge_graph
	{
		long up;  // first uplink
		long count; // number of uplinks
	} r,*q;
	edge *e;
	triangle *t;
	long m;
	
	@<initialize lists |graph|, |crits|, and |to_do|@>@;
	
	m = -1;  // indicate |NULL| uplink
	e = NULL;
	t = sigma;
	@<put eligible edges of |t| on |graph|, |crits|, and |to_do|@>@;
	
	while (!list_is_empty(to_do))
	{
		list_pop(to_do,&m);
		e = (edge *)olist_get_key(graph,m);
		t = r12(e); // get the |t| which is paired with |e|
		@<put eligible edges of |t| on |graph|, |crits|, and |to_do|@>@;
	}
	
	if(flags&1) @<find a persistent critical edge |e| in |crits|@>@;
	else @<find the critical edge |e| with only one grad path so |value(e)| is maximized@>@;
	
	return e;
}

@ @<initialize lists |graph|, |crits|, and |to_do|@>=
if(first==0)
{
	olist_initialize(&graph,sizeof(struct edge_graph));
	list_initialize(&crits,sizeof(long));
	list_initialize(&to_do,sizeof(long));
	first=1;
}
else
{
	olist_clear(graph);
	list_clear(crits);
	list_clear(to_do);
}


@ @<put eligible edges of |t| on |graph|, |crits|, and |to_do|@>=
{
	int i;
	long n;
	int flag;
	edge *ep;
	int livecount;
	
	livecount=0;
	
	for(i=0;i<3;i++)
	{
		ep = get_edge(t,i);
		if (e!=ep && ((is_paired_up(ep)  && !is_deadend(ep)) || is_critical(ep))) 
		{
			livecount++;
						
			if (value(ep) <= value(sigma) - p)  continue;
			
			r.up = m;
			r.count = 0;
			
			n = olist_find_add(graph,(long)ep,&r,&flag);
			q = (struct edge_graph*)olist_entry(graph,n);
			q->count++; 
			
			if(!flag)
			{
				if(is_critical(ep)) list_push(crits,&n);
				else list_push(to_do,&n);
			}
		}
	}
	if(livecount == 0 && e!=NULL) make_deadend(e);
}

@ @<find the critical edge |e| with only one grad path so |value(e)| is maximized@>=
{
	int best_value,val,bestm;
	
	bestm=-1;
	
	while(!list_is_empty(crits))
	{
		list_pop(crits,&m);
		q = (struct edge_graph*)olist_entry(graph,m);
		val = value((edge *)olist_get_key(graph,m));
		if(q->count==1 && (bestm<0 || val > best_value))
			@<replace |bestm| by |m| if there's just one path to it@>@;
	}
	if(bestm>=0) e = (edge *)olist_get_key(graph,bestm);
	else e = NULL;
	if(e!=NULL && grad_path!=NULL) @<put gradient path to |e| on |grad_path|@>@;
	
}

@ @<replace |bestm| by |m| if there's just one path to it@>=
{
	
	while(q->up >= 0)
	{
		q = (struct edge_graph*)olist_entry(graph,q->up);
		if(q->count > 1) break;
	}
	if(q->count==1)
	{
		best_value = val;
		bestm = m;
	}
}

@ @<put gradient path to |e| on |grad_path|@>=
{
	list_clear(grad_path);
	
	m = bestm;
	
	while(m >= 0)
	{
		plist_push(grad_path,(edge *)olist_get_key(graph,m));
		m = ((struct edge_graph*)olist_entry(graph,m))->up;
	}
}


@ Check if |t| is a bad tetrahedron where a gradient path could split and rejoin.
It returns -1 if not, otherwise it returns |i| so that |get_face(t,i)| is bad.

@<Subroutines finding gradient paths@>=
int
splitrejoin(tetrahedron *t)
{
	int i,j;
	edge *e,*ep;
	triangle *s;
	
	for(i=0;i<4;i++)
	{
		s = get_face(t,i);
		if(is_paired_down(s))
		{
			e = r21(s);
			for(j=0;j<3;j++)
			{
				ep = get_edge(s,j);
				if(ep != e
					&&(!is_paired_up(ep)||triangle_in_tetrahedron(r12(ep),t)<0)) break;
			}
			if(j==3) return i;
		}
	}
	return -1;
}

@ Fix up a bad split rejoin tetrahedron |t| with bad face |get_face(t,n)|.
Return the new triangle to which the bad edge is newly paired.

@<Subroutines finding gradient paths@>=
triangle *unsplitrejoin(tetrahedron *t,int n)
{
	triangle *s,*sp;
	edge *e;
	
	s = get_face(t,n);
	e = r21(s);
	sp = other_face(e,s,t);
	
	if(!is_in_K(t)) return NULL;  // can't fix if not in K
	if(is_critical(sp)) return NULL;  // can't fix if critical
	if(is_paired_down(sp)) 
		abort_message("tetrahedron surrounded by 2->1 triangles");  // I think this is impossible
	pair23(s,t);
	pair12(e,sp);
	return sp;
}

@ Find all gradient paths from |sigma| and return them encoded in |edges|.
If |options&3 == 0|, then no pruning is done,
if |options&3 == 1|, then lazy pruning is done,
if |options&3 == 2|, then full pruning is done.
Pruning deletes edges which are not connected to a critical edge
by some gradient path.
If |options&4| is set, then the count field in the edges data structure
will be filled in.
The ordered list |edges| is a list of structures |grad12_struct| below,
one for each edge connected to |sigma|, with key the list index of the edge.
The fields |links| have the following meaning for an edge |e|:
First suppose that |e| is paired with a triangle |t|.
Let $i=0,1,2$ be the index of |e| in |t|.
Then |links[i]| is the list index of an edge |ep| so that
|e| is contained in the triangle paired with |ep|.
For |j!=i|, |links[j]| is the list index of another edge |epp|
so that the |j|-th edge of |t| is contained in the triangle paired with |epp|.
In case |e| is critical, then |links[0]| is the list index
of an edge |ep| so that |e| is in the triangle paired with |ep|
and |links[1]| is the list index of another critical edge.
In other words we have a directed graph $G$ without cycles
so that the vertices of $G$ are edges of our complex which are paired to
triangles, or are critical.
There is a directed path in $G$ from $e'$ to $e$ if $e$ is in the
triangle paired with $e'$ and $e\ne e'$.
Then one of the links of $e$ is an $e'$ pointing to $e$,
and the other two are $e''$ pointing to an edge also pointed to by $e$.
It still sounds confusing but it works.
|flags&3| is the index of the edge in its paired triangle, or 3 if critical;
If complete pruning is requested (|options&3 == 2|) and
|flags&4| is set, then the edge is connected to a critical
edge by a gradient path, i.e., it cannot be pruned.
The return value is the list index of a critical edge,
so that all critical edges are obtained by following the |links[1]|.
The |count| field gives the signed number of gradient paths
going through the edge (counting orientation).
If |flags&|8 is set, then the count field in the edges data structure
is filled in.
If |flags&16| is set, then the count field in the edges data structure
cannot yet be filled in.

@(MorseExtract.h@>=
struct grad12_struct
{
	long links[3];
	long count;
	int flags;
};


@ @<Subroutines finding gradient paths@>=

long
find_all_grad12_paths(triangle *sigma,olist *edges,int options)
{
	int kk = -1;
	edge *ep;
	long todo = -2;
	long this;
	long crit = -1;
	struct grad12_struct *p;
	triangle *f;
	
	olist_clear(edges);
	
	do
	{
		this = todo;
		if(this>=0)
		{
			p = (struct grad12_struct *)olist_entry(edges,this);
			ep = id2edge(olist_get_key(edges,this));
			kk = p->flags&3;
			if(kk)
			{
				todo = p->links[0];
				p->links[0] = -1;
			}
			else
			{
				todo = p->links[1];
				p->links[1] = -1;
			}
			f = r12(ep);
		}
		else	f = sigma;
		@<find-add edges of |f| to |edges| if needed@>@;
	} while(todo>=0);
	
	if((options&3)==2)
		@<prune |edges|@>@;
		
	if(options&4)
		@<fill in count field@>@;
	
	return crit;
}

@ @<find-add edges of |f| to |edges| if needed@>=
{
	int k;
	edge *e;
	int livecount = 0;
	
	for(k=0;k<3;k++)
	{
		if(k==kk) continue;
		e = get_edge(f,k);
		@<find-add |e| to |edges| if needed@>@;
	}	
	if(livecount == 0 && (options&3) && this>=0) 
	{
		make_deadend(f);
		make_deadend(ep);
	}
}


@ @<find-add |e| to |edges| if needed@>=
{
	struct grad12_struct r,*q;
	triangle *t;
	int flag;
	int j;
	long id;
	
	if(is_paired_up(e) && (!is_deadend(e) || (options&3)==0 ))
	{
		livecount++;
		id = olist_find_add(edges,edge_id(e),&r,&flag);
		q = (struct grad12_struct *)olist_entry(edges,id);
		if(flag) // already on list
		{
			p->links[k] = q->links[q->flags&3];
			q->links[q->flags&3] = this;
		}
		else // newly added to list
		{
			q->links[1] = q->links[2] = -1;
			t = r12(e);
			q->flags = edge_in_triangle(e,t);
			q->links[q->flags] = this;
			if(q->flags) q->links[0] = todo;
			else q->links[1] = todo;
			todo = id;
		}
	}
	else if(is_critical(e))
	{
		livecount++;
		id = olist_find_add(edges,edge_id(e),&r,&flag);
		q = (struct grad12_struct *)olist_entry(edges,id);
		if(flag)
		{
			p->links[k] = q->links[0];
			q->links[0] = this;
		}
		else
		{
			q->links[0] = this;
			q->links[1] = crit;
			q->flags = 3;
			crit = id;
		}
	}
}

@ @<prune |edges|@>=
{
	long j;
	struct grad12_struct *p;
	edge *e;
	list *todo_list;
	
	list_initialize(&todo_list,sizeof(long));
	
	todo = crit;
	
	while(todo>=0)
	{
		@<mark all edges on critical paths to this critical edge@>@;
		p = (struct grad12_struct *)olist_entry(edges,todo);
		todo = p->links[1];
	}
	
	list_abandon(&todo_list);

	for(j=0;j<olist_count(edges);j++)
	{
		p = (struct grad12_struct *)olist_entry(edges,j);
		if(p->flags & 4) continue;
		e = id2edge(olist_get_key(edges,j));
		make_deadend(e);
		make_deadend(r12(e));
	}
}

@ @<mark all edges on critical paths to this critical edge@>=
{
	long next,id;
	struct grad12_struct *q;
	int i;
	
	list_clear(todo_list);
	list_push(todo_list,&todo);
	
	while(!list_is_empty(todo_list))
	{
		list_pop(todo_list,&next);
		p = (struct grad12_struct *)olist_entry(edges,next);
		
		if(p->flags&4) continue;
		p->flags|=4;
		
		i = p->flags&3;
		if(i==3) id = p->links[0];
		else id = p->links[i];
		
		e = id2edge(olist_get_key(edges,next));
		
		while(id>=0)
		{
			q = (struct grad12_struct *)olist_entry(edges,id);
			
			ep = id2edge(olist_get_key(edges,id));
			
			if((q->flags&4) == 0)
			{
				list_push(todo_list,&id);
			}
			i = edge_in_triangle(e,r12(ep));
			if(i<0) abort_message("error");
			id = q->links[i];
		}
	}
}


@ @<fill in count field@>=
{
	list *todo_list;
	long id = -1,*idp;
	struct grad12_struct *q;
	int error_check;
	
	list_initialize(&todo_list,sizeof(long));
	list_push(todo_list,&id);
	
	while(!list_is_empty(todo_list))
	{
		list_read_init(todo_list);
		error_check = 1;
		while((idp = (long *)list_read(todo_list))!=NULL)
		{
			if(*idp<0)
			{
				list_read_delete(todo_list);
				kk = -1;
				f = sigma;
			}
			else
			{
				p = (struct grad12_struct *)olist_entry(edges,*idp);
				if(p->flags&24) continue;

				@<see if all incoming faces are counted@>@;
			}
				
			if(f!=NULL)
			{
				error_check = 0;
				@<add edges of |f| to |todo_list|@>@;
			}
		}
		if(error_check && !list_is_empty(todo_list)) abort_message("count error");
	} 
	
	list_abandon(&todo_list);
}

@ @<see if all incoming faces are counted@>=
{	
	edge *e;
	int i;
	
	ep = id2edge(olist_get_key(edges,*idp));
	kk = p->flags&3;
	if(kk==3) id = p->links[0];
	else id = p->links[kk];
	
	while(id>=0)
	{
		q = (struct grad12_struct *)olist_entry(edges,id);
		if((q->flags&8)==0) break;
		e = id2edge(olist_get_key(edges,id));
		i = edge_in_triangle(ep,r12(e));
		id = q->links[i];
	}
	
	if(id<0)
		@<count incoming faces to |ep|@>@;
	else p->flags |= 16;
	
	if(id>=0 || is_critical(ep)) f = NULL;
	else	f = r12(ep);
}

@ @<add edges of |f| to |todo_list|@>=
{
	int k;
	edge *e;
	
	for(k=0;k<3;k++)
	{
		if(k==kk) continue;
		e = get_edge(f,k);
		@<add |e| to |todo_list|@>@;
	}	
}

@ @<add |e| to |todo_list|@>=
{
	int flag;
	long qid;
	
	if(is_critical(e) || (is_paired_up(e) && (!is_deadend(e) || (options&3)==0 )))
	{
		qid = olist_find_add(edges,edge_id(e),NULL,&flag);
		if(!flag) abort_message("count error 2");
		q = (struct grad12_struct *)olist_entry(edges,qid);
		if((q->flags & 32)==0)
		{
			list_read_insert(todo_list,&qid);
			q->flags |= 32;
		}
		else q->flags &= 0xffffffef;
	}
}

@ @<count incoming faces to |ep|@>=
{
	edge *e;
	int i;
	long ct=0;
	long orient;
	triangle *t;
	
	p->flags &= 0xffffffef;
	p->flags |= 8;
	list_read_delete(todo_list);

	kk = p->flags&3;
	if(kk==3) id = p->links[0];
	else id = p->links[kk];
	
	while(id>=0)
	{
		q = (struct grad12_struct *)olist_entry(edges,id);
		e = id2edge(olist_get_key(edges,id));
		t = r12(e);
		i = edge_in_triangle(ep,t);
		id = q->links[i];
		ct += edge_orient(t,i)*q->count;
	}
	if (id==-2)
	{
		ct += edge_orient(sigma,edge_in_triangle(ep,sigma));
	}
	if(is_critical(ep)) p->count = ct;
	else p->count = -ct * edge_orient(r12(ep),p->flags&3);
}

@ Find all gradient paths ending at tau.
Only now the key is a triangle id, rather than an edge id.
The list |crits| is a list of indices of items in |triangles|
which are critical triangles.
@d is_xdeadend(s) ((s)->type & 0x1000)
@d make_xdeadend(s) ((s)->type |= 0x1000)

@<Subroutines finding gradient paths@>=
long
find_all_backward_grad12_paths(edge *tau,olist *triangles,list *crits,int options)
{
	int kk = -1;
	triangle *f;
	edge *ep;
	long todo = -1;
	long this;
	struct grad12_struct *p;
	long start = -1;
	list *todo_list;

	olist_clear(triangles);
	if(crits!=NULL) list_clear(crits);
	if((options&3)==2) list_initialize(&todo_list,sizeof(long));

	do
	{
		this = todo;
		if(this>=0)
		{
			p = (struct grad12_struct *)olist_entry(triangles,this);
			f = id2triangle(olist_get_key(triangles,this));
			ep = r21(f);
			kk = p->flags&3;
			todo = p->links[kk];
			p->links[kk] = -1;
		}
		else ep = tau;
		@<find-add cofaces of |ep| to |triangles| if needed@>@;
	} while(todo>=0);
	if((options&3)==2)
		@<prune |triangles|@>@;
	return start;
}

@ @<find-add cofaces of |ep| to |triangles| if needed@>=
{
	tetrahedron *t;
	int link_count = 0;
	
	f = ep->links[2];
	t = coface(f,0);
	do
	{
		if(is_in_K(f) && (is_critical(f) ||
			((!is_xdeadend(f) || (options&3)==0 ) &&
				is_paired_down(f) && ep!=r21(f))))
		{
			link_count++;
			@<find-add |f| to |triangles|@>@;
		}
		t = other_coface(f,t);
		f = other_face(ep,f,t);  
	} while (f!=ep->links[2]);
	if(link_count==0 && (options&3) && this>=0)
	{
		make_xdeadend(ep);
		make_xdeadend(r12(ep));
	}
}

@ @<find-add |f| to |triangles|@>=
{
	struct grad12_struct r,*q;
	long id;
	int flag;
	int k;
	
	id = olist_find_add(triangles,triangle_id(f),&r,&flag);
	q = (struct grad12_struct *)olist_entry(triangles,id);
	k = edge_in_triangle(ep,f);
	if(this>=0)
		p = (struct grad12_struct *)olist_entry(triangles,this);
	if(flag) // already on list
	{
		if(q->links[k]>=0) abort_message("backgrad error 1");
		if(this<0) abort_message("backgrad error 2");
		q->links[k] = p->links[kk];
		p->links[kk] = id;
	}
	else // newly added to list
	{
		q->links[0] = q->links[1] = q->links[2] = -2;
		if(this>=0) 
		{
			q->links[k] = p->links[kk];
			p->links[kk] = id;
		}
		else
		{
			q->links[k] = start;
			start = id;
		}
		if(is_critical(f))
		{
			q->flags = 3;
			if(crits!=NULL) list_push(crits,&id);
			if((options&3)==2) list_push(todo_list,&id);
		}
		else
		{
			q->flags = edge_in_triangle(r21(f),f);
			q->links[q->flags] = todo;
			todo = id;
		}
	}
}

@ @<prune |triangles|@>=
{
	long id;
	int k;
	edge *e;
	struct grad12_struct r,*q;
	triangle *f;
	int flag;
	
	while(!list_is_empty(todo_list))
	{
		list_pop(todo_list,&id);
		p = (struct grad12_struct *)olist_entry(triangles,id);
		if(p->flags&4) continue;
		p->flags|=4;
		kk = p->flags&3;
		f = id2triangle(olist_get_key(triangles,id));
		for(k=0;k<3;k++)
		{
			if(k==kk || p->links[k] < -1) continue;
			e = get_edge(f,k);
			if(is_paired_up(e) && !is_deadend(e) )
			{
				if(is_xdeadend(e)) abort_message("back grad error 3");
				id = olist_find_add(triangles,triangle_id(r12(e)),&r,&flag);
				q = (struct grad12_struct *)olist_entry(triangles,id);
				if(flag==0) 
				{
					q->flags = edge_in_triangle(e,r12(e));
					printf("back grad puzzle\n");
				}
				if((q->flags&4)==0) list_push(todo_list,&id);
			}
		}	
	}
	list_abandon(&todo_list);

	for(id=0;id<olist_count(triangles);id++)
	{
		p = (struct grad12_struct *)olist_entry(triangles,id);
		if(p->flags & 4) continue;
		f = id2triangle(olist_get_key(triangles,id));
		make_deadend(f);
		make_deadend(r21(f));
	}
}

@* Persistent critical simplices.  This is experimental code, trying to identify
critical simplices which might generate persistent homology.
So far in random examples all but a few critical simplices end up being persistent.

@<find persistent critical simplices@>=

@<make all critical vertices persistent@>@;
@<find persistent critical edges@>@;
@<find persistent critical triangles@>@;
@<find persistent critical tetrahedra@>@;



@ @<make all critical vertices persistent@>=
{
	vertex *v;

	list_read_init(crit[0]);
	while ((v = (vertex *)plist_read(crit[0]))!=NULL) 
	{
		if(!is_critical(v)) list_read_delete(crit[0]);
		else make_persistent(v);
	}
}

@ This is called after canceling is done, so if |v| is not |NULL|
then we know that the two gradient paths from |e| end at the same
vertex |v| and |value(v)| $\ge$ |value(e)-p|.
@<find persistent critical edges@>=
{
	edge *e;
	vertex *v;

	list_read_init(crit[1]);
	while ((e = (edge *)plist_read(crit[1]))!=NULL) 
	{
		if(!is_critical(e)) list_read_delete(crit[1]);
		else 
		{
			v = FindGrad01(get_vertex(e,0),value(e)-p);
			if(v==NULL) make_persistent(e);
		}
	}
}

@ @<find persistent critical tetrahedra@>=
@<make all critical tetrahedra persistent@>@;
{
	triangle *t;
	tetrahedron *te;

	list_read_init(crit[2]);
	while ((t = (triangle *)plist_read(crit[2]))!=NULL) 
	{
		if(is_persistent(t)) 
		{
			te = FindGrad23(coface(t,0),value(t)+p);
			if(te!=NULL) unmake_persistent(te);
		}
	}
}

@ @<make all critical tetrahedra persistent@>=
{
	tetrahedron *t;

	list_read_init(crit[3]);
	while ((t = (tetrahedron *)plist_read(crit[3]))!=NULL) 
	{
		if(!is_critical(t)) list_read_delete(crit[3]);
		else make_persistent(t);
	}
}

@ @<find persistent critical triangles@>=
{
	triangle *t;
	void *q;

	list_read_init(crit[2]);
	while ((t = (triangle *)plist_read(crit[2]))!=NULL) 
	{
		if(!is_critical(t)) list_read_delete(crit[2]);
		else 
		{
			q = FindGradPaths12(t,p,NULL,1);
			if(q==NULL) make_persistent(t);
		}
	}
}

@ @<find a persistent critical edge |e| in |crits|@>=
{
	edge *ee;
	
	e = NULL;
	while(!list_is_empty(crits))
	{
		list_pop(crits,&m);
		ee = (edge*)olist_get_key(graph,m);
		if(is_persistent(ee))
		{
			e = ee;
			break;
		}
	}
}

@* Unordered Lists.
Unordered lists work as LIFO stacks.
We can push to and pop from a |list|.
We can find the $n$-th entry on a |list|.
We can read the entries of a |list| one by one, deleting those we no longer want.

@f list int

@(MorseExtract.h@>=
typedef struct 
{
	long *body;
	unsigned int size;  // size of each entry, not including padding
	unsigned int padded_size;  // size of each entry, in long words
	unsigned long body_length;  // capacity of body
	unsigned long length;	// number of entries on list
	unsigned long read_index;
	unsigned long read_delete_index;
} list;
#define list_count(l) ((l)->length)
#define list_is_empty(l) ((l)->length == 0)
#define list_clear(l)  ((l)->length = 0)
#define list_entry(l,n) ((l)->body + (n)*((l)->padded_size))

@ Initialize a list, setting up storage
@<List functions@>=
void list_initialize(list **l,unsigned int sz)
{
	*l = (list *)malloc(sizeof(list));
	if((*l)==NULL) abort_message("Out of memory");
	(*l) -> size = sz;
	(*l) -> padded_size = (sz+sizeof(long)-1)/sizeof(long);
	(*l) -> length = 0;
	(*l) -> body = NULL;
	(*l) -> body_length = 0;
}

  
@ Abandon a list, freeing storage used.
@<List functions@>=

void list_abandon(list **l)
{
	if ((*l)->body_length>0) free((*l)->body);
	free(*l);
	*l = NULL;
}



@ |list_push| pushes a new entry pointed to by |q| to the top of the list.
The variant |plist_push| pushes the pointer |q| directly to the list
(as opposed to the contents pointed to by |q|).

@<List functions@>=

void *list_push(list *l,void *q)
{
	long *p;
	if(l->body_length <= l->length)
	{
		l->body_length  = l->length + 1000;
		l->body = (long *)realloc(l->body,l->body_length*l->padded_size*sizeof(long));
		if(l->body==NULL) abort_message("out of memory");
	}
	
	p = list_entry(l, l->length);

	memcpy(p,q,l->size);
	l->length++;
	return p;
}

void plist_push(list *l, void *q)
{
	void *p = q;
	list_push(l,&p);
}


@ |list_pop| pops the top entry into the region pointed to by |q|,
and deletes it from the list.
It copies the top entry to |q|.
The variant |plist_pop| returns the entry from a pointer list directly.

@<List functions@>=

void list_pop(list *l,void *q)
{
	long *p;
	
	if (l->length==0) abort_message("pop from empty list");
	
	l->length--;
	p = list_entry(l, l->length);
	memcpy(q,p,l->size);
}

void *plist_pop(list *l)
{
	void *p;
	
	list_pop(l,&p);
	return p;
}


@ Routine to read the entries of a list one by one,
and deleting those you don't wish to retain on the list.
To use it, first call |list_read_init|.
Then repeated calls to |list_read| will return the entries on the list,
until a |NULL| is returned signifying the end.
If you wish to delete the entry you last read, just call |list_read_delete|.
If your list is a list of pointers, then 
the companion pointer list routine |plist_read| returns the pointer,
rather than its address.
Since the end of the list is signaled by returning a |NULL| pointer
it will only work well if the list of pointers
contains no |NULL| pointers, otherwise you might stop too early.
So in this case you would need some alternate method to
recognize the end of the list.

{\bf Warnings:} if you use |list_read_delete| then you must read through
to the end of the list.  
@(MorseExtract.h@>=
#define list_read_init(l) ((l)->read_index = (l)->read_delete_index = 0)
#define list_read_delete(l) ((l)->read_delete_index--)

@ @<List functions@>=
void *list_read(list *l)
{
	void *p,*pp;
	
	if(l->read_index >= l->length)
	{
		// we are at the end of the list
		l->length = l->read_index = l->read_delete_index;
		return(NULL);
	}
	
	pp = list_entry(l,l->read_delete_index);
	

	if(l->read_index != l->read_delete_index) 
	{
		p = list_entry(l,l->read_index);
		memcpy(pp, p, l->size);
	}
	
	l->read_delete_index++;
	l->read_index++;
	return pp;
}

void *plist_read(list *l)
{
	void **r;
	r = list_read(l);
	if(r!=NULL) return *r;
	else return NULL; 
}

@ The following |list_read_insert| will insert an item in a list
you are reading, placing it just before the most recently read item
if there is room, or if not, placing it at the end.
@<List functions@>=
void *list_read_insert(list *l,void *p)
{
	if(l->read_index == l->read_delete_index) 
		list_push(l,p);
	else
	{
		memcpy(list_entry(l,l->read_delete_index), p, l->size);
		l->read_delete_index++;
	}
}


@* Ordered Lists. 
Ordered lists are implemented naively and could no doubt be made more efficient.
Items on an |olist| are ordered by an integer key.

Portability issue: these are sometimes used with a simplex pointer
as the key, so that requires |sizeof(long) = sizeof(simplex *)|.
When this is done, of course we don't really care about the resulting order,
we are just using ordered lists to find entries quickly.
One could also implement such lists by adding auxiliary fields to
each simplex, but then memory usage is increased.
@f olist int

@ @(MorseExtract.h@>=
struct olist_key
{
	long high;
	long low;
	long eq;
	long k;  // the list is ordered by this key |k|
};

typedef struct 
{
	list *keys;
	list *entries;
	long top;
	long free;
	int sz;
} olist;
#define olist_is_empty(l) ((l)->top < 0)
#define olist_entry(l,n) list_entry((l)->entries,n)
#define olist_get_key(l,n)  (((struct olist_key *)list_entry((l)->keys,n))->k)
#define olist_count(l) list_count((l)->entries)

@ @<List functions@>=
void olist_initialize(olist **l,int sz)
{
	*l = (olist *)malloc(sizeof(olist));
	if((*l)==NULL) abort_message("Out of memory");
	list_initialize(&((*l)->keys),sizeof(struct olist_key));
	list_initialize(&((*l)->entries),sz);
	(*l)->top = -1;
	(*l)->free = -1;
	(*l)->sz = sz;
}

@ @<List functions@>=
void olist_abandon(olist **l)
{
	list_abandon(&((*l)->keys));
	list_abandon(&((*l)->entries));
	free(*l);
	*l = NULL;
}

@ @<List functions@>=
void olist_clear(olist *l)
{
	list_clear(l->keys);
	list_clear(l->entries);
	l->top = -1;
	l->free = -1;
}


@ Pop the smallest entry from the list, return its key, and put the list entry in |p|.
@<List functions@>=
long olist_min(olist *l,void *p)
{
	struct olist_key *q,*r;
	long min_index;

	
	min_index = l->top;
	if(min_index<0) abort_message("min of empty list");
	q = (struct olist_key *)list_entry(l->keys,min_index);
	if(q->low >= 0)
	{
		while(q->low >= 0)
		{
			r = q;
			min_index = q->low;
			q = (struct olist_key *)list_entry(l->keys,min_index);
		}
		if(q->eq>=0)
		{
			r=q;
			min_index = q->eq;
			q = (struct olist_key *)list_entry(l->keys,min_index);
			r->eq = q->eq;
		}
		else r->low = q->high;
	}
	else if(q->eq>=0)
	{
		r=q;
		min_index = q->eq;
		q = (struct olist_key *)list_entry(l->keys,min_index);
		r->eq = q->eq;
	}
	else l->top = q->high;

	memcpy(p,list_entry(l->entries,min_index),l->sz);
	
	q->low = l->free;
	l->free = min_index;
	
	return q->k;
}


@ Add an entry |p| with key |m| to the ordered list.
Return the index of the entry.
@<List functions@>=

long olist_add(olist *l,long m,void *p)
{
	struct olist_key *q,*r;
	long b,*last;
	

	b = olist_next_free(l);
	@<point |r| the |b|-th entry and copy |p| to it@>@;
	
	last = &(l->top);
	
	while((*last)>=0)
	{
		q = (struct olist_key *)list_entry(l->keys,*last);
		if (m < q->k) last = &(q->low);
		else if(m > q->k) last = &(q->high);
		else 
		{
			r->eq = *last;
			r->low = q->low;
			r->high = q->high;
			q->low = -1;
			q->high = -1;
			break;
		}
	}
	*last = b;
	return b;
}

@ This searches the ordered list |l| for an entry with key |m|.
If it finds one, it returns its index and sets |flag|.
If it doesn't find one, and |p| is not NULL,
it creates an entry with key |m|
and moves |*p| into that entry and resets |flag|.
@d olist_next_free(l) ((l)->free>=0)? (l)->free: list_count((l)->entries)
@<List functions@>=

long olist_find_add(olist *l,long m,void *p,int *flag)
{
	struct olist_key  *q,*r;
	long b,*last;

	last = &(l->top);
	
	while((*last)>=0)
	{
		q = (struct olist_key *)list_entry(l->keys,*last);
		if (m < q->k) last = &(q->low);
		else if(m > q->k) last = &(q->high);
		else // if duplicate key, return entry
		{
			if(flag!=NULL) *flag =1;
			return *last;
		}
	}
	if(p!=NULL)
	{
		b = *last = olist_next_free(l);
		@<point |r| the |b|-th entry and copy |p| to it@>@;
	}
	else b = -1;
	
	if(flag!=NULL) *flag = 0;
	return b;
}


@ @<point |r| the |b|-th entry and copy |p| to it@>=
{
	struct olist_key rr;
	
	rr.high = rr.low = rr.eq = -1;
	rr.k = m;
	
	if(l->free>=0)
	{
		memcpy(list_entry(l->entries,b),p,l->sz);
		r = (struct olist_key *)list_entry(l->keys,b);
		l->free = r->low;
		memcpy(r,&rr,sizeof(struct olist_key));
	}
	else
	{
		list_push(l->entries,p);
		r = (struct olist_key *)list_push(l->keys,&rr);
	}
}

@* Navigating Simplicial Complexes.
Simplices of all dimensions have a common initial part, 
there is a word with all sorts of bit flags.
Yes I know I should implement this with all these bit fields as different
entries in the struct.  Maybe later.
The meaning of the bits in the type field is:
\item{$\cdot$} Bits 0 and 1 give the dimension.
\item{$\cdot$} Bit 2 is set if the simplex is in $K$.
\item{$\cdot$} Bit 3 is set if the simplex is critical.
\item{$\cdot$} Bit 4 is set if the simplex is paired with a codimension one face.
\item{$\cdot$} Bits 5 and 6 give the index of that face if bit 4 is set. 
	If bit 3 is set then bit 5 is set the simplex is persistent critical.
	If bits 3 and 4 are not set and the simplex is a triangle
	then bit 5 gives the index of the coface to which the simplex is paired.
\item{$\cdot$} Bit 7 is set if the |h| field is valid.
\item{$\cdot$} Bit 8 is set in an edge if the edge  is paired up and it is known 
	that all gradient paths from the edge do not lead to critical edges.
\item{$\cdot$} Bit 9 is set in an edge |e| when that edge is processed by
   |@<extract the lower Star of ...@>| with |v= get_vertex(e,0)|.
\item{$\cdot$} Bit 10 is set in an edge |e| when that edge is processed by
   |@<extract the lower Star of ...@>| with |v= get_vertex(e,1)|.
\item{$\cdot$} Bits 11-15 are available for other uses.  For example, if
this program were adapted to use
the CGAL representation of a complex, edges and triangles are not represented
explicitly but the unused bits could indicate, say, which face of a tetrahedron
the triangle is, or which of its two vertices an edge is.

The field |h| is the maximum of the values of the vertices of the simplex,
if bits 2 and 7 of |type| are set.  Not essential, but there was this unused
space, so what the heck.
I expect 16 bit resolution of the function is good enough.

The remaining  field |links| is a variable sized array of pointers to other simplices
 and is used to navigate around the complex.
The entries in |links| have the following meaning:
\item{} For a vertex:
\itemitem{$\cdot$}{|links[0]|} is an edge containing the vertex, or |NULL| if there is none.
\item{} For an edge:
\itemitem{$\cdot$}{|links[0-1]|} are the vertices in the edge.
\itemitem{$\cdot$}{|links[2]|} is a triangle containing the edge.
\item{} For a triangle:
\itemitem{$\cdot$}{|links[0-2]|} are the edges contained in the triangle.
\itemitem{$\cdot$}{|links[3-4]|} are the two tetrahedra which contain the triangle.
\item{} For a tetrahedron:
\itemitem{$\cdot$}{|links[0-3]|} are the triangles contained in the tetrahedron.

@f vertex int
@f edge int
@f triangle int
@f tetrahedron int

@ @(MorseExtract.h@>=

typedef struct vertexstruct
{
	short type;  
	short h;
	
	void *links[1];
} vertex;

typedef struct edgestruct
{
	short type;  
	short h;
	
	void *links[3];
} edge;

typedef struct trianglestruct
{
	short type;  
	short h;
	
	void *links[5];
} triangle;

typedef struct tetrahedronstruct
{
	short type;  
	short h;
	
	void *links[4];
} tetrahedron;

#define get_vertex(e,i) (vertex *)((e)->links[i])
#define get_edge_vertices(e,vl) ((vl)[0]=get_vertex(e,0), (vl)[1]=get_vertex(e,1))
#define dimension(s) ((s)->type&3)
#define is_critical(t) ((t)->type&8)
#define other_vertex_in_edge(u,e) ((get_vertex(e,0) == u)?get_vertex(e,1):get_vertex(e,0))
#define coface(t,i) (tetrahedron *)((t)->links[3+i])
#define get_edge(t,i) (edge *)((t)->links[i])
#define get_face(t,i) (triangle *)((t)->links[i])
#define is_in_K(v) ((v)->type&4)
#define value(t) (int)((((t)->type&0x80)!=0)?(t)->h:max_value(t,dimension(t)))
#define other_coface(s,t) 	(tetrahedron *)(((t) == (s)->links[3])? (s)->links[4]: (s)->links[3])
#define is_deadend(s) ((s)->type & 0x100)
#define make_deadend(s) ((s)->type |= 0x100)



@ Now some functions to manipulate and navigate these simplices

@d set_value(t,f) ((t)->type|=0x80, (t)->h= f)
@<Simplex functions@>=

void get_triangle_vertices(triangle *t,vertex *vlist[3])
// return a list of the vertices of |t|
{
	vertex *v;
	edge *e;
	
	e = get_edge(t,0);

	get_edge_vertices(e,vlist);
	
	e = get_edge(t,1);
	
	v= get_vertex(e,1);
	if(v!=vlist[0] && v!=vlist[1]) vlist[2]=v;
	else vlist[2]= get_vertex(e,0);
}

void get_tetrahedron_vertices(tetrahedron *t,vertex *vlist[4])
// return a list of the vertices of |t|
{
	vertex *vlistp[3];
	int i;
	
	
	get_triangle_vertices(get_face(t,0),vlist);
	get_triangle_vertices(get_face(t,1),vlistp);
	
	for(i=0;i<3;i++)
	{
		if(vlistp[i]!=vlist[0] && vlistp[i]!=vlist[1] && vlistp[i]!=vlist[2])
		{
			vlist[3]=vlistp[i];
			break;
		}
	}
}

@ Return $\pm1$ depending on the orientation of the i-th edge of |t|.
@<Simplex functions@>=
int edge_orient(triangle *t,int i)
{
	vertex *vl[2],*vlp[2];

	if(i==0) return 1;
	get_edge_vertices(get_edge(t,0),vl);
	get_edge_vertices(get_edge(t,i),vlp);
	if(vlp[0]==vl[1] || vlp[1]==vl[0]) return 1;
	return -1;
}


@ a test for when the value of one vertex is greater
than another. The tiebreaker is just the address of the vertex.
@<Simplex functions@>=
long vertex_compare(vertex *v,vertex *w)
{
	if(value(v) < value(w)) return -1;
	else if (value(v) > value(w)) return 1;
	else return (long)v - (long)w;
}

@ Find the smallest or largest valued vertex in a list of vertices of length n
@d smallest_vertex_in_edge(e) (vertex_compare(get_vertex(e,0),get_vertex(e,1))<0 ? 
	get_vertex(e,0) : get_vertex(e,1))
@d largest_vertex_in_edge(e) (vertex_compare(get_vertex(e,0),get_vertex(e,1))<0 ? 
	get_vertex(e,1) : get_vertex(e,0))
@<Simplex functions@>=

int smallest_vertex(vertex *v[2],int n)
{
	int i,j;
	for(i=0,j=1;j<n;j++)
	{
		if(vertex_compare(v[j],v[i])<0) i=j;
	}
	return i;
}

int largest_vertex(vertex *v[2],int n)
{
	int i,j;
	for(i=0,j=1;j<n;j++)
	{
		if(vertex_compare(v[j],v[i])>0) i=j;
	}
	return i;
}

@ This routine tests whether the simplex |t| is in the lower Star of the vertex |v|.
It does tiebreaking if two vertices in |t| have the same value.

@<Simplex functions@>=

int is_in_lower_Star(tetrahedron *t,vertex *v)
{
	vertex *vlist[4];
	
	get_tetrahedron_vertices(t,vlist);
	
	return v == vlist[largest_vertex(vlist,4)];
}




@ These routines test whether a given simplex is a codimension one face of another,
and if it is, returns the index of that face.
If it is not, then -1 is returned.
@(MorseExtract.h@>=
#define vertex_in_edge(v,e) (((v) == get_vertex(e,0))?0:(((v) == get_vertex(e,1))?1:-1))
#define edge_in_triangle(e,t) (((e) == get_edge(t,0))?0:(((e) == get_edge(t,1))?1: \
(((e) == get_edge(t,2))?2:-1)))


@ @<Simplex functions@>=
int triangle_in_tetrahedron(triangle *s,tetrahedron *t)
{
	int i;
	
	for(i=0;i<4;i++)
	{
		if(s == get_face(t,i)) return i;
	}
	return -1;
}


@ These routines take a codimension one face of a simplex
and find the unique other codimension one face which also contains
the given codimension 2 face.


@<Simplex functions@>=
edge *other_edge(vertex *v,edge *f,triangle *t)
// return the edge of t which contains v but is not f
{
	int i;
	edge *e;
	
	
	for(i=0;i<3;i++)
	{
		e = get_edge(t,i);
		if(e == f) continue;
		if(get_vertex(e,0)==v || get_vertex(e,1)==v) 	return e;
	}
	
	return NULL;  // error, v not in t
}

triangle *
other_face(edge *e,triangle *f,tetrahedron *t)
// Find the face of t which contains e but is not f
{
	int i;
	triangle *s;
	
	for(i=0;i<4;i++)
	{
		s = get_face(t,i);
		if(f!=s && edge_in_triangle(e,s)>=0) return s;
	}
	return NULL;  // error, e not in f or f not in t
}



@ Here are various routines for critical simplices.
In |unmake_critical|, the 
formerly critical simplex is not actually removed from
the list of critical simplices.
So when running through |crit| you may need to 
check if the entries are in fact critical.
The routine |clean_crit| will delete all the noncritical entries from |crit|,

@d unmake_critical(s) (s)->type &= 0xff87
@d make_critical(s) ((s)->type &= 0xff87, (s)->type |= 8, plist_push(crit[dimension(s)],s))
@d is_persistent(t) ((t)->type&0x20)
@d make_persistent(t) ((t)->type|=0x20)
@d unmake_persistent(t) ((t)->type&=0xffdf)
@<Simplex functions@>=

void clean_crit(void)
{
	tetrahedron *t;
	triangle *f;
	edge *e;
	vertex *v;
	list *l;

	l = crit[0];
	list_read_init(l);
	while ((v = (vertex *)plist_read(l))!=NULL)	
	{
		if(!is_critical(v)) list_read_delete(l);
	}
	l = crit[1];
	list_read_init(l);
	while ((e = (edge *)plist_read(l))!=NULL)	
	{
		if(!is_critical(e)) list_read_delete(l);
	}
	l = crit[2];
	list_read_init(l);
	while ((f = (triangle *)plist_read(l))!=NULL)	
	{
		if(!is_critical(f)) list_read_delete(l);
	}
	l = crit[3];
	list_read_init(l);
	while ((t = (tetrahedron *)plist_read(l))!=NULL)	
	{
		if(!is_critical(t)) list_read_delete(l);
	}

	
}


@  Routines to pair simplices.  
The routine |is_paired_up| finds whether the simplex |t| is paired with a
simplex of one higher dimension.
Likewise |is_paired_down| finds whether the simplex |t| is paired with a
simplex of one lower dimension.
The routines |rij| return the dimension |j| simplex paired with the
dimension |i| simplex |t|.
The routines |pairij| pair up the i simplex |s0| and the j simplex |s1|.

@(MorseExtract.h@>=
#define is_paired_up(t) (((t)->type&0x1C)==4)
#define is_paired_down(t) (((t)->type&0x1C)==0x14)
#define r32(t) get_face(t,((t)->type&0x60)>>5)
#define r21(t) get_edge(t,((t)->type&0x60)>>5)
#define r10(t) get_vertex(t,((t)->type&0x20)>>5)
#define r23(t) coface(t,((t)->type&0x20)>>5)
#define r12(t) (triangle *)((t)->links[2])
#define r01(t) (edge *)((t)->links[0])

@ @<Simplex functions@>=
void pair01(vertex *s0,edge *s1)
{
	int n;
	
	s0->type&=0xFF87;
	s1->type&=0xFF87;
	
	s0->links[0] = s1;
	n = vertex_in_edge(s0,s1);
	s1->type|=(16+(n<<5));
} 

void pair12(edge *s0,triangle *s1)
{
	int n;
	
	s0->type&=0xFF87;
	s1->type&=0xFF87;
	
	s0->links[2] = s1;
	n = edge_in_triangle(s0,s1);
	s1->type|=(16+(n<<5));
} 

void pair23(triangle *s0,tetrahedron *s1)
{
	int n;
	
	s0->type&=0xFF87;
	s1->type&=0xFF87;
	
	if(s0->links[4] == s1)  s0->type|=0x20;
	n = triangle_in_tetrahedron(s0,s1);
	s1->type|=(16+(n<<5));
} 

@ find the maximum value of the vertices of the simplex |s|
@<Simplex functions@>=
int max_value(void *s,int d)
{
	tetrahedron *te;
	triangle *t;
	edge *e;
	int m0,m1;
	
	if(s==NULL) return -0x8000;  // smallest value
	
	switch(d)
	{
		case 3:
			te = (tetrahedron *)s;
			m0 = value(get_face(te,0));
			m1 = value(get_face(te,1));
			break;
		case 2:
			t = (triangle *)s;
			m0 = value(get_edge(t,0));
			m1 = value(get_edge(t,1));
			break;
		case 1:
			e = (edge *)s;
			m0 = value(get_vertex(e,0));
			m1 = value(get_vertex(e,1));
			break;
		default:
			abort_message("max value of simplex of unknown dimension");
			break;
	}
	return max(m0,m1);
}

@ find the minimum value of the vertices of the simplex |s|
@<Simplex functions@>=
int min_value(void *s,int d)
{
	tetrahedron *te;
	triangle *t;
	edge *e;
	int m0,m1;
	
	if(s==NULL) return 0x7fff;  // largest value
	switch(d)
	{
		case 3:
			te = (tetrahedron *)s;
			m0 = min_value(get_face(te,0),2);
			m1 = min_value(get_face(te,1),2);
			break;
		case 2:
			t = (triangle *)s;
			m0 = min_value(get_edge(t,0),1);
			m1 = min_value(get_edge(t,1),1);
			break;
		case 1:
			e = (edge *)s;
			m0 = value(get_vertex(e,0));
			m1 = value(get_vertex(e,1));
			break;
		default:
			abort_message("min value of simplex of unknown dimension");
			break;
	}
	return min(m0,m1);
}


@* Debugging routines.
@<Debugging output@>=
void abort_message(char *s)
{
	printf("Error: %s",s);
	exit(1);
}





@* Constructing Complexes. 
The following routines allocate a new simplex.
The last parameter |flag| is $>0$ if the simplex is in |K|,
is $=0$ if the simplex is not in |K|, and is $<0$
if you want it to be in |K| if and only if all its
faces are in |K|.

@(MorseExtract.h@>=
#define vertex_id(v) ((v)-vertexlist)
#define edge_id(e) ((e)-elist)
#define triangle_id(f) ((f)-flist)
#define tetrahedron_id(t) ((t)-tlist)
#define id2vertex(id) (vertexlist+(id))
#define id2edge(id) (elist +(id))
#define id2triangle(id) (flist +(id))
#define id2tetrahedron(id) (tlist +(id))

@ @<Subroutines for constructing complexes@>=
vertex *new_vertex(long id,short val,short flag)
{
	vertex *v;
	
	v = id2vertex(id);
	
	v->type=(flag)?0x84:0x80; // indicate in $K$ or not
	v->links[0]=NULL;
	v->h = val;
	
	return v;
}

@ Construct a new edge
@<Subroutines for constructing complexes@>=
edge *new_edge(long id,long vids[2],short flag)
{
	edge *e;
	short flagp;
	vertex *v0,*v1;
	
	e = id2edge(id);
	v0 = id2vertex(vids[0]);
	v1 = id2vertex(vids[1]);
	e->links[0]=v0;
	e->links[1]=v1;
	e->links[2]=NULL;
	
	
	e->h = max(v0->h,v1->h);
	
	flagp = (flag && is_in_K(v0) && is_in_K(v1));
	
	
	e->type=(flagp)?5:1;
	
	if(flag>0 && flagp==0) 
		abort_message("vertices of edge in K must be in K");
		
	if(v0->links[0]==NULL) v0->links[0]=e;
	if(v1->links[0]==NULL) v1->links[0]=e;
	
	return e;
}

@ @<Subroutines for constructing complexes@>=
triangle *new_triangle(long id,long eids[3],short flag)
{
	triangle *t;
	vertex *vl[3],*v[2];
	edge *e0,*e1,*e2;
	short flagp;
	
	e0 = id2edge(eids[0]);
	e1 = id2edge(eids[1]);
	e2 = id2edge(eids[2]);
	
	get_edge_vertices(e0,vl);
	get_edge_vertices(e1,v);
	if(vertex_in_edge(v[0], e0)<0) vl[2]=v[0];
	else if(vertex_in_edge(v[1], e0)<0) vl[2]=v[1];
	else abort_message("new triangle has more than three vertices");
	
	get_edge_vertices(e2,v);
	if(v[0]!=vl[0] && v[0]!=vl[1] &&v[0]!=vl[2])
		abort_message("new triangle has more than three vertices-1");
	if(v[1]!=vl[0] && v[1]!=vl[1] && v[1]!=vl[2]) 
		abort_message("new triangle has more than three vertices-2");
	
	
	t = id2triangle(id);
	t->links[0]=e0;
	t->links[1]=e1;
	t->links[2]=e2;
	
	t->h = max(e0->h,e1->h);

	flagp = (flag && is_in_K(e0) && is_in_K(e1) && is_in_K(e2));

	t->type = (flagp)?6:2;
	t->links[3]=NULL;
	t->links[4]=NULL;
	
	if(e0->links[2]==NULL) e0->links[2]=t;
	if(e1->links[2]==NULL) e1->links[2]=t;
	if(e2->links[2]==NULL) e2->links[2]=t;

	if(flag >0 && flagp==0) 
		abort_message("edges of triangle in K must be in K");
	return t;
}

@ @<Subroutines for constructing complexes@>=
tetrahedron *new_tetrahedron(long id,long fids[4],short flag)
{
	tetrahedron *t;
	vertex *vl[4],*vlp[3];
	int i;
	short flagp;
	triangle *f0,*f1,*f2,*f3;
	
	f0 = id2triangle(fids[0]);
	f1 = id2triangle(fids[1]);
	f2 = id2triangle(fids[2]);
	f3 = id2triangle(fids[3]);
	
	get_triangle_vertices(f0,vl);
	get_triangle_vertices(f1,vlp); 
	for(i=0;i<3;i++)
	{
		if(vlp[i]!=vl[0] && vlp[i]!=vl[1] && vlp[i]!=vl[2]) 
		{
			vl[3] = vlp[i];
			break;
		}
	}
	@<make sure |vlp| entries are all in |vl|@>@;
	get_triangle_vertices(f2,vlp); 
	@<make sure |vlp| entries are all in |vl|@>@;
	get_triangle_vertices(f3,vlp);
	@<make sure |vlp| entries are all in |vl|@>@;

	
	t = id2tetrahedron(id);


	t->links[0]=f0;
	t->links[1]=f1;
	t->links[2]=f2;
	t->links[3]=f3;
	
	t->h = max(f0->h,f1->h);

	flagp = (flag && is_in_K(f0) && is_in_K(f1) && is_in_K(f2) && is_in_K(f3));

	t->type=(flagp)?7:3; 
	
	if(f0->links[3]==NULL) f0->links[3]=t;
	else if(f0->links[4]==NULL) f0->links[4]=t;
	else abort_message("a triangle can be in only two tetrahedra");
	if(f1->links[3]==NULL) f1->links[3]=t;
	else if(f1->links[4]==NULL) f1->links[4]=t;
	else abort_message("a triangle can be in only two tetrahedra-1");
	if(f2->links[3]==NULL) f2->links[3]=t;
	else if(f2->links[4]==NULL) f2->links[4]=t;
	else abort_message("a triangle can be in only two tetrahedra-2");
	if(f3->links[3]==NULL) f3->links[3]=t;
	else if(f3->links[4]==NULL) f3->links[4]=t;
	else abort_message("a triangle can be in only two tetrahedra-3");
	
	if(flag>0 && flagp==0) 
		abort_message("faces of tetrahedron in K must be in K");
	return t;
}

@ @<make sure |vlp| entries are all in |vl|@>=
for(i=0;i<3;i++)
{
	if(vlp[i]!=vl[0] && vlp[i]!=vl[1] && vlp[i]!=vl[2] && vlp[i]!=vl[3]) 
		abort_message("new tetrahedron has more than four vertices");
}






@ @<output simplex pairings@>=
{
	FILE *df;
	vertex *v;
	edge *e;
	triangle *f;
	tetrahedron *t;
	
	df = fopen(argv[2],"w");
	if(df==NULL) abort_message("Cannot open output file");

	for(i=0;i<4;i++) // output number if i simplices and critical i simplices
	{
		fprintf(df,"%d %d\n",num_simp[i],n[i]);
	}
	
	for(i=0,v=vertexlist;i<num_simp[0];i++,v++)
	{
		if(!is_in_K(v)) fprintf(df,"x%d\n",i);
	}
	fprintf(df,"z\n");
	for(i=0,e=elist;i<num_simp[1];i++,e++)
	{
		if(!is_in_K(e)) fprintf(df,"x\n");
		else if(is_critical(e)) fprintf(df,"c\n");
		else if(is_paired_up(e)) fprintf(df,"a%d\n",triangle_id(r12(e)));
		else fprintf(df,"b%d\n",vertex_id(r10(e)));
	}
	fprintf(df,"z\n");
	for(i=0,f=flist;i<num_simp[2];i++,f++)
	{
		if(!is_in_K(f)) 
		{
			if(is_in_K(get_edge(f,0))&&is_in_K(get_edge(f,1))&&is_in_K(get_edge(f,2)))
				fprintf(df,"x%d\n",i);
		}
	}
	fprintf(df,"z\n");
	for(i=0,t=tlist;i<num_simp[3];i++,t++)
	{
		if(!is_in_K(t)) fprintf(df,"x\n");
		else if(is_critical(t)) fprintf(df,"c\n");
		else fprintf(df,"b%d\n",triangle_id(r32(t)));
	}
	fclose(df);
}

@ @<Subroutines for constructing complexes@>=
void
read_in_complex(FILE *df)
{
	long i,n[4];
	int res;

	@<read in vertices@>@;
	@<read in edges@>@;
	@<read in triangles@>@;
	@<read in tetrahedra@>@;
}	

@ @<read in vertices@>=
{
	int val;
	int vertex_in_K,c;

	res = fscanf(df,"%ld\n",n);
	if(res!=1) abort_message("bad number of vertices");
	
	
	num_simp[0]=n[0];
	vertexlist = (vertex *)malloc(n[0]*sizeof(vertex));
	if(vertexlist==NULL) abort_message("Out of memory");
	
	for(i=0;i<n[0];i++)
	{
		c = getc(df);
		if(c=='x') vertex_in_K=0;
		else
		{
			vertex_in_K = 1;
			ungetc(c,df);
		}
		res = fscanf(df,"%ld\n",&val);
		if(res!=1) abort_message("bad vertex");
		new_vertex(i,val,vertex_in_K);
	}
}

@ @<read in edges@>=
{
	long ids[2];
	int edge_in_K,c;
	
	res = fscanf(df,"%ld\n",n+1);
	if(res!=1) abort_message("bad number of edges");
	num_simp[1]=n[1];
	
	elist = (edge *)malloc(n[1]*sizeof(edge));
	if(elist==NULL) abort_message("Out of memory");
	
	for(i=0;i<n[1];i++)
	{
		c = getc(df);
		if(c=='x') edge_in_K=0;
		else
		{
			edge_in_K = -1;
			ungetc(c,df);
		}
		res = fscanf(df,"%ld %ld\n",ids,ids+1);
		if(res!=2) abort_message("bad edge");
		new_edge(i,ids,edge_in_K);
	}
}

@ @<read in triangles@>=
{
	long ids[3];
	int face_in_K,c;
	
	res = fscanf(df,"%ld\n",n+2);
	if(res!=1) abort_message("bad number of triangles");
	num_simp[2]=n[2];
	
	flist = (triangle *)malloc(n[2]*sizeof(triangle));
	if(flist==NULL) abort_message("Out of memory");
	
	for(i=0;i<n[2];i++)
	{
		c = getc(df);
		if(c=='x') face_in_K=0;
		else
		{
			face_in_K = -1;
			ungetc(c,df);
		}
		res = fscanf(df,"%ld %ld %ld\n",ids,ids+1,ids+2);
		if(res!=3) abort_message("bad triangle");
		new_triangle(i,ids,face_in_K);
	}
}

@ @<read in tetrahedra@>=
{
	long ids[4];
	int tetrahedron_in_K,c;
	
	res = fscanf(df,"%ld\n",n+3);
	if(res!=1) abort_message("bad number of tetrahedra");
	num_simp[3]=n[3];
	
	tlist = (tetrahedron *)malloc(n[3]*sizeof(tetrahedron));
	if(tlist==NULL) abort_message("Out of memory");
	
	for(i=0;i<n[3];i++)
	{
		c = getc(df);
		if(c=='x') tetrahedron_in_K=0;
		else
		{
			tetrahedron_in_K = -1;
			ungetc(c,df);
		}
		res = fscanf(df,"%ld %ld %ld %ld\n",ids,ids+1,ids+2,ids+3);
		if(res!=4) abort_message("bad tetrahedron");
		new_tetrahedron(i,ids,tetrahedron_in_K);
	}
}

@ @(MorseExtract.h@>=
#ifndef in_MorseExtract
extern list *crit[4];
extern long num_simp[4];  // number of simplices
extern vertex *vertexlist;
extern edge *elist;
extern triangle *flist;
extern tetrahedron *tlist;
#endif


void Extract(int p);
void ExtractCancel1(int p);
void ExtractCancel2(int p);
void ExtractCancel3(int p);
void read_in_complex(FILE *df);
void *plist_read(list *l);
void get_triangle_vertices(triangle *t,vertex *vlist[3]);
void get_tetrahedron_vertices(tetrahedron *t,vertex *vlist[4]);
vertex *FindGrad01(vertex *u,int m);
tetrahedron *FindGrad23(tetrahedron *tau,int m);
long find_all_grad12_paths(triangle *sigma,olist *edges,int flags);
void list_initialize(list **l,unsigned int sz);
void list_abandon(list **l);
void *list_push(list *l,void *q);
void list_pop(list *l,void *q);
void *list_read(list *l);
void plist_push(list *l, void *q);
void *plist_pop(list *l);
void *plist_read(list *l);
void olist_initialize(olist **l,int sz);
void olist_abandon(olist **l);
void olist_clear(olist *l);
long olist_min(olist *l,void *p);
long olist_add(olist *l,long m,void *p);
long olist_find_add(olist *l,long m,void *p,int *flag);
void abort_message(char *s);

@t}\vfil\eject{@>