@i common.w @* 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). @c #include "Morse.h" #include "globals.h" @@; @@; @@; @*1 Local Cancellation. |LocalCancel| cancels simplices in the lower Star of a vertex |v|. 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. The input parameter |n| is the lowest list index of all the critical 2 simplices in the lower star of |v|. As it happens, when this is called all succeeding critical 2 simplices will also be in the lower star of |v|. @+= void LocalCancel(simplex_id v,UInt32 n) {@/ simplex_id f,s; // triangles simplex_id e[2],ep[2]; // edges simplex_id tid[2],tidp[2]; // tetrahedra int j,i; Boolean in;@/ list_read_init_at(crit[2],n); while (!id_is_null(f=id_list_read(crit[2]))) // run through the list of new critical 2 simplices { @; @; if(!edges_equal(ep[0],ep[1])) // can we cancel? { @; list_read_delete(crit[2]); } else { // We can't cancel with an edge, so see if we can cancel with a 3 simplex for(i=0;i<2;i++) { @; } if(tid[0]!=tid[1]) { @; list_read_delete(crit[2]); } } } } @ @= e[0]=get_edge(f,1); e[1]=get_edge(f,2); @ @= for(i=0;i<2;i++) { ep[i]=e[i]; while(!edge_is_critical(ep+i)) { if(edge_is_paired_down(ep+i)) { ep[i]=NULL_ID; // runs into $B_1$ so no cancelation possible break; } s = r12(ep+i); // get paired triangle ep[i] = other_edge(v,ep[i],s); } } @ @= tidp[i] = tid[i] = coface(f,i); s = f; in = false; while (tetrahedron_is_in_K(tid[i])) // while we are in K { in = is_in_lower_Star(tid[i],v); if(!in) break; // see if exit lower Star if(tetrahedron_is_critical(tid[i])) break; // see if reached critical s = r32(tid[i]); tid[i] = other_coface(s,tid[i]); } if(!tetrahedron_is_in_K(tid[i]) || !in ) tid[i]=NULL_ID; @ @= if(id_is_null(ep[0])) LocalCancel12(v,ep[1],e[1],f); else if(id_is_null(ep[1])) LocalCancel12(v,ep[0],e[0],f); else if(min_value(ep[0],1) > min_value(ep[1],1)) LocalCancel12(v,ep[0],e[0],f); else LocalCancel12(v,ep[1],e[1],f); @ @= if (id_is_null(tid[0])) Cancel23(f,tidp[1],tid[1]); else if (id_is_null(tid[1])) Cancel23(f,tidp[0],tid[0]); else if(min_value(tid[1],3) > min_value(tid[0],3)) Cancel23(f,tidp[0],tid[0]); else Cancel23(f,tidp[1],tid[1]); @*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$. @+= void LocalCancel12(simplex_id v,simplex_id sigma,simplex_id sigmap,simplex_id tau) { simplex_id u; // edge simplex_id w, wp; // triangles u=sigmap; w=tau; unmake_triangle_critical(&w); while(!edge_is_critical(&u)) { wp = r12(&u); pair12(&u,w,true); w=wp; u=other_edge(v,u,w); } if(!edges_equal(u,sigma)) abort_message("localcancel12 internal error"); unmake_edge_critical(&u); pair12(&u,w,true); } @ The following routine cancels a critical 2 simplex $\sigma$ with a critical 3 simplex $\tau$. The 3 simplex |taup| contains $\sigma$ and is such that the gradient path from $\tau$ to $\sigma$ goes through |taup|. @+= void Cancel23(simplex_id sigma,simplex_id taup,simplex_id tau) { simplex_id s,sp; simplex_id t; s=sigma; t=taup; unmake_triangle_critical(&s); while(!tetrahedron_is_critical(t)) { sp=r32(t); pair23(&s,t,true); s=sp; t=other_coface(s,t); } if(t!=tau) abort_message("cancel23 internal error"); unmake_tetrahedron_critical(t); pair23(&s,t,true); } @ The following routine cancels a critical vertex |v| with a critical edge $\kappa$. The n is such that the gradient path from |v| to $\kappa$ goes through |get_vertex(kappa,n)|. @+= void Cancel01(simplex_id v,int n,simplex_id kappa) { simplex_id u; // vertex simplex_id e,ep; // edges u = get_vertex(kappa,n); e = kappa; unmake_edge_critical(&e); while(!vertex_is_critical(u)) // while non critical { ep = r01(u); pair01(u,e,true); e = ep; u = other_vertex_in_edge(u,e); } if(u!=v) abort_message("cancel01 internal error"); unmake_vertex_critical(v); pair01(v,e,true); } @ Cancel a critical 1 simplex $\sigma$ with a critical 2 simplex $\kappa$. It assumes that all the edges in the gradient path are in the list |grad_path|, in order. It pushes any 2 simplices with changed pairings to the stack |changed|. @+= void Cancel12(simplex_id sigma,simplex_id kappa,list *grad_path,list *changed) { simplex_id t,nextt; // triangles simplex_id e; // edge t = kappa; unmake_triangle_critical(&t); do { e = id_list_pop(grad_path); if(edge_is_critical(&e)) break; nextt = r12(&e); pair12(&e,t,true); id_list_push(changed,t); t = nextt; } while(true); if(!edges_equal(e,sigma)) abort_message("cancel12 internal error"); unmake_edge_critical(&e); pair12(&e,t,true); id_list_push(changed,t); } @ @(cancel.h@>= void LocalCancel(simplex_id v,UInt32 n); void LocalCancel12(simplex_id v,simplex_id sigma,simplex_id sigmap,simplex_id tau); void Cancel01(simplex_id v,int n,simplex_id kappa); void Cancel12(simplex_id sigma,simplex_id kappa,list *grad_path,list *changed); void Cancel23(simplex_id sigma,simplex_id taup,simplex_id tau); void ExtractCancel1(SInt32 p); void ExtractCancel2(SInt32 p); void ExtractCancel3(SInt32 p); simplex_id FindGrad01(simplex_id u,SInt32 m); simplex_id FindGrad23(simplex_id tau,SInt32 m); simplex_id FindGrad23orientation(simplex_id sigma,int *orientation); simplex_id FindGradPaths12(simplex_id sigma, SInt32 p, list *grad_path, int flags); Boolean splitrejoin(simplex_id t,list *changed); SInt32 find_all_grad12_paths(simplex_id sigma,hlist *edges,int options); SInt32 find_all_backward_grad12_paths(simplex_id tau,hlist *triangles,list *crits,int options); @*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. @= void ExtractCancel1(SInt32 p) {@/ hlist *c0; // list of vertices to cancel with hlist *c1; // list of edges which can be cancelled with vertices struct ccrit1 { simplex_id s; SInt32 c[2]; } s; int i; simplex_id v[2]; SInt32 thisk; SInt32 evalue; UInt8 persist_scaling; @@; hlist_initialize(&c0,sizeof(SInt32),113,1,0); hlist_initialize(&c1,sizeof(struct ccrit1),1+(p>>persist_scaling),1,persist_scaling); @; while (!hlist_is_empty(c1)) {@t\1@>@/ hlist_min(c1,&thisk,&s); // put best candidate in |s|, persistence in |thisk| @@; i = @; evalue = value(s.s); if (vertex_is_critical(v[i]) && thisk == evalue-vertex_value(v[i])) @@; else @@;@t\2@>@/ }@/ hlist_abandon(&c0); hlist_abandon(&c1);@t\2@>@/ } @ |hlist_min| will only find the actual minimum if its hash function is one to one (and order prerserving). If |p| is large, this would lead to very large hash tables, so instead we scale back the persistence precision so the hash table size is $<8096$. @= { UInt32 mask = 0xffffe000; if(p>65535) p = 65535; else if(p<0) p=0; for(persist_scaling=0;persist_scaling<4;persist_scaling++,mask<<=1) if((p&mask)==0) break; } @ @= { simplex_id e; list_read_init(crit[1]); while ((e=id_list_read(crit[1]))!=-1) // run through all critical edges { if(!edge_is_critical(&e)) list_read_delete(crit[1]); else { evalue = value(e); v[0] = FindGrad01(get_vertex(e,0),evalue-p); // find the two gradient paths v[1] = FindGrad01(get_vertex(e,1),evalue-p); // descending from |e| if(v[0]!=v[1]) { @; } } } } @ @= {@t\1@>@/ SInt32 link; SInt32 persist; link = -1; s.s=e; s.c[0] = (id_is_null(v[0]))? -1: hlist_find_add(c0,v,&link,NULL); s.c[1] = (id_is_null(v[1]))? -1: hlist_find_add(c0,v+1,&link,NULL);@/ i = @; persist = evalue - vertex_value(v[i]); hlist_add(c1,&persist ,&s);@t\2@>@/ } @ @= v[0] = (s.c[0]<0)? NULL_ID : hlist_get_key(c0,s.c[0],0); v[1] = (s.c[1]<0)? NULL_ID : hlist_get_key(c0,s.c[1],0); @ @= { Cancel01(v[i],i,s.s); *((SInt32 *)hlist_entry(c0,s.c[i])) = s.c[1-i]; // now any gradient path to |v[i]| goes to other vertex |v[1-i]| } @ @= {@t\1@>@/ SInt32 m; evalue = value(s.s); for(i=0;i<2;i++) { if(s.c[i]<0) continue; @@; } if( s.c[1]!=s.c[0]) {@t\1@>@/ i = @; m = evalue - vertex_value(v[i]); if( m < p) hlist_add(c1,&m,&s);@t\2@>@/ } @t\2@>@/ } @ @= { SInt32 a; for(a = s.c[i]; a >= 0; ) { v[i] = hlist_get_key(c0,a,0); if( vertex_is_critical(v[i])) break; a = *((SInt32 *)hlist_entry(c0,a)); } if(a!=s.c[i]) { *((SInt32 *)hlist_entry(c0,s.c[i])) = a; } if (vertex_value(v[i]) <= (evalue-p)) { v[i] = NULL_ID; s.c[i] = -1; } else s.c[i] = a; } @ @= (id_is_null(v[0]))? 1:@t}$\par${@> (id_is_null(v[1]))? 0:@t}$\par${@> (vertex_value(v[0]) > vertex_value(v[1]))? 0: @t}$\par${@> (vertex_value(v[0]) < vertex_value(v[1]))? 1: @t}$\par${@> vertex_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. @= void ExtractCancel3(SInt32 p) {@/ hlist *c3; // list of tetrahedra to cancel hlist *c2; // list of triangles which cancel with tetrahedra struct ccrit2 { simplex_id s; SInt32 c[2]; } r; int i; simplex_id t[2]; SInt32 thisk; SInt32 fvalue,tvalue; UInt8 persist_scaling; @@; hlist_initialize(&c3,sizeof(SInt32),113,1,0); hlist_initialize(&c2,sizeof(struct ccrit2),1+(p>>persist_scaling),1,persist_scaling); @; while (!hlist_is_empty(c2)) {@t\1@>@/ hlist_min(c2,&thisk,&r); // put best candidate in |r|, persistence in |thisk| @@; i = @; fvalue = value(r.s); tvalue = value(t[i]); if (tetrahedron_is_critical(t[i]) && thisk == tvalue-fvalue) @@; else @@;@t\2@>@/ }@/ hlist_abandon(&c2); hlist_abandon(&c3);@t\2@>@/ } @ @= {@t\1@>@/ simplex_id s; list_read_init(crit[2]); while ((s = id_list_read(crit[2]))!=-1) // go through all critical triangles {@/ if(!triangle_is_critical(&s)) list_read_delete(crit[2]); else { fvalue = value(s); // find the beginnings of the two gradient paths ending at |s| t[0] = FindGrad23(coface0(s),fvalue+p); t[1] = FindGrad23(coface1(s),fvalue+p);@/ if(t[0]!=t[1]) @@; } }@t\2@>@/ } @ @= {@t\1@>@/ SInt32 link; SInt32 persist; link = -1; r.s=s; r.c[0] = (id_is_null(t[0]))? -1: hlist_find_add(c3,t,&link,NULL); r.c[1] = (id_is_null(t[1]))? -1: hlist_find_add(c3,t+1,&link,NULL);@/ i = @; persist = value(t[i]) - fvalue; hlist_add(c2, &persist, &r);@t\2@>@/ } @ @= t[0] = (r.c[0]<0)? NULL_ID : hlist_get_key(c3,r.c[0],0); t[1] = (r.c[1]<0)? NULL_ID : hlist_get_key(c3,r.c[1],0); @ @= { Cancel23(r.s,coface(r.s,i),t[i]); *((SInt32 *)hlist_entry(c3,r.c[i])) = r.c[1-i]; } @ @= {@t\1@>@/ SInt32 m; fvalue = value(r.s); for(i=0;i<2;i++) { if(r.c[i]<0) continue; @@; } if(r.c[1]!=r.c[0]) {@t\1@>@/ i = @; m = value(t[i]) - fvalue; if( m < p) hlist_add(c2,&m,&r);@t\2@>@/ }@t\2@>@/ } @ @= { SInt32 a; for(a = r.c[i]; a >= 0; ) { t[i] = hlist_get_key(c3,a,0); if( tetrahedron_is_critical(t[i])) break; a = *((SInt32 *)hlist_entry(c3,a)); } *((SInt32 *)hlist_entry(c3,r.c[i])) = a; if(value(t[i])>=fvalue+p) { t[i] = NULL_ID; r.c[i] = -1; } else r.c[i] = a; } @ @= (id_is_null(t[0]))? 1:@t}$\par${@> (id_is_null(t[1]))? 0:@t}$\par${@> (value(t[0]) > value(t[1]))? 1: @t}$\par${@> (value(t[0]) < value(t[1]))? 0: @t}$\par${@> tetrahedron_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. @= void ExtractCancel2(SInt32 p) { list *changed; // triangles whose pairing was changed by Cancel12 list *grad_path; // edges in gradient path of a canceling pair hlist *goodpairs; // cancelable critical triangles UInt32 lastp; UInt8 persist_scaling; @@; list_initialize(&changed,sizeof(simplex_id)); list_initialize(&grad_path,sizeof(simplex_id)); hlist_initialize(&goodpairs,sizeof(simplex_id), 1+(p>>persist_scaling),1,persist_scaling); lastp = 0; do { @@; if(hlist_is_empty(goodpairs)) break; @@; } while(true); list_abandon(&changed); list_abandon(&grad_path); hlist_abandon(&goodpairs); } @ @= { simplex_id f; simplex_id e; UInt32 persist; list_read_init(crit[2]); while (!id_is_null(f = id_list_read(crit[2]))) // go through all critical triangles { if(!triangle_is_critical(&f)) list_read_delete(crit[2]); else { e = FindGradPaths12(f,p,NULL,0); if(!id_is_null(e)) // if there is a gradient path from |f| to |e| { persist = value(f)-value(e); hlist_add(goodpairs,&persist,&f); } } } } @ @= { simplex_id f; simplex_id e; UInt32 bp,thisp; while(!hlist_is_empty(goodpairs)) { hlist_min(goodpairs,&bp,&f); e = FindGradPaths12(f,p,grad_path,0); if(id_is_null(e)) continue; // no gradient paths from |t| found thisp = value(f)-value(e); if(thisp<=bp) { list_clear(changed); Cancel12(e,f,grad_path,changed); // cancel |e| and |t| @; if(thisp %d\n",lastp,thisp); lastp = thisp; } else hlist_add(goodpairs,&thisp,&f); } } @ 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. @= { int i,j,k; simplex_id t; k = 0; while (!list_is_empty(changed) && k<10000) { f = id_list_pop(changed); if(id_is_null(f)) continue; for (i=0;i<2;i++) { t = coface(f,i); if( tetrahedron_is_in_K(t) && !tetrahedron_is_critical(t) ) { if(splitrejoin(t,changed)) k++; } } } if (!list_is_empty(changed)) printf("May be infinite loop in split-rejoin"); } @ 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. @d debug9 1 @= Boolean splitrejoin(simplex_id t,list *changed) { int i,j; simplex_id e,ep; simplex_id s,sp; for(i=0;i<4;i++) { s = get_face(t,i); if(triangle_is_paired_down(&s)) { e = r21(&s); for(j=0;j<3;j++) { ep = get_edge(s,j); if(!edges_equal(ep,e) &&(!edge_is_paired_up(&ep)||triangle_in_tetrahedron(r12(&ep),t)<0)) break; } if(j==3) { sp = r32(t); @q+9 if(!triangles_equal(r32(t),sp)) abort_message("Split-rejoin error"); @> pair23(&s,t,true); pair12(&e,sp,true); id_list_push(changed,sp); return true; } } } return false; } @*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_ID| if $h(v)\le m$. @= simplex_id FindGrad01(simplex_id u,SInt32 m) { simplex_id v; simplex_id e; v=u; while(!vertex_is_critical(v) && vertex_value(v) > m) { e = r01(v); v = other_vertex_in_edge(v,e); } if (vertex_value(v) <= m) return NULL_ID; 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_ID| 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. @= simplex_id FindGrad23(simplex_id tau,SInt32 m) { simplex_id t; simplex_id s; if (id_is_null(tau)) return NULL_ID; t=tau; while( tetrahedron_is_in_K(t) && value(t) < m ) { if (tetrahedron_is_critical(t)) return t; // reached critical t s = r32(t); if (triangle_is_deadend(&s)) { make_tetrahedron_deadend(t); return NULL_ID; } t = other_coface(s,t); if (tetrahedron_is_deadend(t)) { make_triangle_deadend(&s); return NULL_ID; } } if(!tetrahedron_is_in_K(t) && t!=tau) make_triangle_deadend(&s); return NULL_ID; } @ Find the starting critical tetrahedron of a gradient path ending at the critical triangle $\sigma$ and going through the tetrahedron |id2tetra_id(sigma)|. Return $\pm 1$ in |orientation| depending on the change in orientation along the path. @= simplex_id FindGrad23orientation(simplex_id sigma,int *orientation) { simplex_id t; simplex_id s; int orient; t=id2tetra_id(sigma); orient = face_index(sigma); while( tetrahedron_is_in_K(t) ) { if (tetrahedron_is_critical(t)) { if(orient&1) *orientation=-1; else *orientation=1; return t; // reached critical t } s = r32(t); orient += face_index(s); if (triangle_is_deadend(&s)) { make_tetrahedron_deadend(t); return NULL_ID; } s=other_face_id(s); t = id2tetra_id(s); if (tetrahedron_is_deadend(t)) { make_triangle_deadend(&s); return NULL_ID; } orient++; orient += face_index(s); } return NULL_ID; } @ 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. 8/30/2005 Possible problem: perhaps there is a second path which dips down below sigvalue-p and thus is not detected? @= simplex_id FindGradPaths12(simplex_id sigma, SInt32 p, list *grad_path, int flags) { static hlist *graph; static list *to_do; static list *crits; static Boolean first = true; struct edge_graph { SInt32 up; // first uplink SInt32 count; // number of uplinks simplex_id e; } r,*q; simplex_id e; simplex_id f; simplex_id vlist[2]; SInt32 m; SInt32 sigvalue = value(sigma); @@; m = -1; // indicate |NULL| uplink e = NULL_ID; f = sigma; @@; while (!list_is_empty(to_do)) { list_pop(to_do,&m); q = (struct edge_graph *)hlist_entry(graph,m); e = q->e; f = r12(&e); // get the |f| which is paired with |e| @@; } if(flags&1) @@; else @@; return e; } @ @= if(first) { hlist_initialize(&graph,sizeof(struct edge_graph),703,2,0); list_initialize(&crits,sizeof(SInt32)); list_initialize(&to_do,sizeof(SInt32)); first=false; } else { hlist_clear(graph); list_clear(crits); list_clear(to_do); } @ @= { int i; SInt32 n; Boolean flag; simplex_id ep; int livecount; livecount=0; for(i=0;i<3;i++) { ep = get_edge(f,i); if (!edges_equal(ep,e) && !edge_is_paired_downx(ep) && !edge_is_deadendx(&ep)) { livecount++; // The following statement must be eliminated // otherwise one might think a pair is cancellable // when in fact there is a second path between them // which dips down below the limit and later comes back up. // Eliminating gives a big performancs hit though. // Perhaps there is a way to determine a valid lower bound // to replace |sigvalue - p| and thus regain most of the performance. // As it is, we end up determining the whole descending disc. // |if (value(ep) <= sigvalue - p) continue;| r.up = m; r.count = 0; r.e = ep; get_edge_vertices(ep,vlist); n = hlist_find_add(graph,vlist,&r,&flag); q = (struct edge_graph*)hlist_entry(graph,n); q->count++; if(!flag) { // |ep| is newly added to the list |graph| if(edge_is_critical(&ep)) list_push(crits,&n); else list_push(to_do,&n); // The next line might change q->e to the optimal descriptor of |ep| q->e = ep; } } } if(livecount == 0 && !id_is_null(e)) make_edge_deadend(&e); } @ @= { SInt32 best_value,val,bestm; bestm=-1; while(!list_is_empty(crits)) { list_pop(crits,&m); q = (struct edge_graph*)hlist_entry(graph,m); val = value(q->e); if (val <= sigvalue - p) continue; if(q->count==1 && (bestm<0 || val > best_value)) @e|@>@; } if(bestm>=0) { q = (struct edge_graph*)hlist_entry(graph,bestm); e = q->e; } else e = NULL_ID; if(!id_is_null(e) && grad_path!=NULL) @@; } @ @e|@>= { while(q->up >= 0) { q = (struct edge_graph*)hlist_entry(graph,q->up); if(q->count > 1) break; } if(q->count==1) { best_value = val; bestm = m; } } @ @= { list_clear(grad_path); m = bestm; while(m >= 0) { q = (struct edge_graph*)hlist_entry(graph,m); id_list_push(grad_path,q->e); m = q->up; } } @ Find all gradient paths from the triangle |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 2 keys, the ids of the edges vertices. Each edge |e| is either critical or paired with a triangle |t|. The fields |links| have the following meaning for an edge |e|: One of the links we call the uplink, and the remaining links are side links. If the edge is critical, the uplink is |link[0]|, otherwise it is |links[i]| if $i$ is the index of |e| in |t|. This uplink is the list index of an edge |ep| so that |e| is contained in the triangle paired with |ep|, in fact, the start of a linked list of such edges. The links of this list are in the side links. So for example, if the uplink of |e| points to |ep| and |e| has index |j| in the triangle paired with |ep|, then |links[j]| (of |ep|?) points to another edge |epp| paired with a triangle containing |e|, or is -1 if this is the last item on the list. In case |e| is critical, then the side link |links[1]| is the list index of another critical edge. We could also look at |edges| as the vertices of a directed graph $G$ without cycles. There is a directed path in $G$ from $ep$ to $e$ if $e$ is in the triangle paired with $ep$ and $e\ne ep$. Then the uplink of $e$ is an $ep$ pointing to $e$, and the other two are $epp$ 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 of |find_all_grad12_paths| 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. @(cancel.h@>= struct grad12_struct { SInt32 links[3]; SInt32 fake; SInt32 count; SInt32 count2; simplex_id e; int flags; }; @ @= void check_grad12_struct(hlist *edges,int options) { int i,k,kk,id; struct grad12_struct *p,*q; simplex_id e,ep,f,vlist[2]; for(i=0;ie; kk = p->flags&3; if(kk!=3) { if(!edge_is_paired_up(&e)) abort_message("check grad12 err 1"); f = r12(&e); for(k=0;k<3;k++) { if(k==kk) continue; ep = get_edge(f,k); if((options&3)!=0 && edge_is_deadend(&ep)) continue; @@; } } } } @ @= { simplex_id epp; int j; get_edge_vertices(ep,vlist); id = hlist_find(edges,vlist); if(id<0) { if(!edge_is_deadend(&ep) && !edge_is_paired_down(&ep) ) abort_message("check grad12 err 2"); } else { q = (struct grad12_struct *)hlist_entry(edges,id); id = q->links[(q->flags&3)%3]; while(id>=0) { if(id==i) break; q = (struct grad12_struct *)hlist_entry(edges,id); epp = q->e; j = edge_verts_in_triangle(vlist,r12(&epp)); if(j<0) abort_message("check grad12 err 3"); id = q->links[j]; } if(id<0) abort_message("check grad12 err 4"); } } @ @= SInt32 find_all_grad12_paths(simplex_id sigma,hlist *edges,int options) { int kk = -1; simplex_id ep; SInt32 todo = -2; SInt32 this; SInt32 critid = -1; struct grad12_struct *p; simplex_id f; hlist_clear(edges); do { @@; @@; } while(todo>=0); check_grad12_struct(edges,options); if((options&3)==2) @@; if(options&4) @@; if(options&8) @@; return critid; } @ If nonnegative, |todo| is the list index of an edge paired to an unprocessed triangle. All these unprocessed triangles are in a list, linked by either |links[0]| or |links[1]|. @= { this = todo; if(this>=0) { p = (struct grad12_struct *)hlist_entry(edges,this); ep = p->e; 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; } @ @= { int k; simplex_id e; int livecount = 0; for(k=0;k<3;k++) { if(k==kk) continue; e = get_edge(f,k); @@; } if(livecount == 0 && (options&3)!=0 && this>=0) { make_triangle_deadend(&f); make_edge_deadend(&ep); } } @ If the edge |e| is paired up or critical we increment |livecount|, we point the uplink of |e| to |ep| (which has list index |this|) and we set the sidelink of |ep| to the former uplink of |e|, thus linking up the list. If |e| is paired up but not on the list already, we link |e| into the todo list. On the other hand, if |e| is critical but not on the list already we put it at the start of the linked list of critical edges. @= { struct grad12_struct r,*q; simplex_id s,vlist[2]; Boolean flag; int j; SInt32 id; // For efficiency, first do the following tests // which don't require an optimal edge id. if( ( (options&3)==0 || !edge_is_deadend(&e)) && !edge_is_paired_down(&e) ) { livecount++; get_edge_vertices(e,vlist); r.fake=r.count=r.count2=0; id = hlist_find_add(edges,vlist,&r,&flag); q = (struct grad12_struct *)hlist_entry(edges,id); if(flag) // already on list { p = (struct grad12_struct *)hlist_entry(edges,this); p->links[k] = q->links[(q->flags&3)%3]; q->links[(q->flags&3)%3] = this; } else // newly added to list { // note the following test will change |e| to the optimal id of the edge // as a side effect of the test. if(edge_is_paired_up(&e)) { q->links[1] = q->links[2] = -1; s = r12(&e); q->e = e; q->flags = edge_in_triangle(e,s); if(q->flags<0) abort_message("help"); q->links[q->flags] = this; if(q->flags) q->links[0] = todo; else q->links[1] = todo; todo = id; } else if(edge_is_critical(&e)) { q->e = e; q->links[0] = this; q->links[1] = critid; q->flags = 3; critid = id; } } } } @ @= { SInt32 j; simplex_id e,s; list *todo_list; list_initialize(&todo_list,sizeof(SInt32)); todo = critid; while(todo>=0) { @@; p = (struct grad12_struct *)hlist_entry(edges,todo); todo = p->links[1]; } list_abandon(&todo_list); for(j=0;jflags & 4) continue; e = p->e; make_edge_deadend(&e); s = r12(&e); make_triangle_deadend(&s); } } @ @= { SInt32 next,id; struct grad12_struct *q; simplex_id vl[2]; 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 *)hlist_entry(edges,next); if(p->flags&4) continue; p->flags|=4; i = p->flags&3; id = p->links[i%3]; e = p->e; vl[0] = hlist_get_key(edges,next,0); vl[1] = hlist_get_key(edges,next,1); while(id>=0) { q = (struct grad12_struct *)hlist_entry(edges,id); ep = q->e; if((q->flags&4) == 0) { list_push(todo_list,&id); } i = edge_verts_in_triangle(vl,r12(&ep)); if(i<0) abort_message("error"); id = q->links[i]; } } } @ @= { list *todo_list; SInt32 id = -1,*idp; struct grad12_struct *q; int error_check; list_initialize(&todo_list,sizeof(SInt32)); list_push(todo_list,&id); while(!list_is_empty(todo_list)) { list_read_init(todo_list); error_check = 1; while((idp = (SInt32 *)list_read(todo_list))!=NULL) { if(*idp<0) { list_read_delete(todo_list); kk = -1; f = sigma; } else { p = (struct grad12_struct *)hlist_entry(edges,*idp); if(p->flags&8) {list_read_delete(todo_list); continue;} //if(p->flags&24) continue; @@; } if(!id_is_null(f)) { error_check = 0; @@; } } if(error_check && !list_is_empty(todo_list)) { // debugging code follows @@; // I don't understand it, sometimes the todo list // is not empty but will never be shortened. // here we check that nevertheless, all critical edges // have a valid count. // If so, there must be some deadend edges unmarked I think. todo = critid; while(todo>=0) { p = (struct grad12_struct *)hlist_entry(edges,todo); if(p->flags&8==0) break; todo = p->links[1]; } if(todo>=0) { // If we reach here not all critical edges have a valid count. list_read_init(todo_list); while((idp = (SInt32 *)list_read(todo_list))!=NULL) { simplex_id vl[2]; simplex_id e; int i; p = (struct grad12_struct *)hlist_entry(edges,*idp); ep = p->e; vl[0] = hlist_get_key(edges,*idp,0); vl[1] = hlist_get_key(edges,*idp,1); printf("%d %d ",*idp,p->flags); id = p->links[(p->flags&3)%3]; while(id>=0) { q = (struct grad12_struct *)hlist_entry(edges,id); printf("%d-",q->flags); e = q->e; i = edge_verts_in_triangle(vl,r12(&e)); id = q->links[i]; } printf("\n"); } abort_message("count error"); } else {printf("extra todo_list\n"); break;} } } list_abandon(&todo_list); } @ @d edge_orient(i) (((i)==1)?-1:1) @= { simplex_id e; int i; simplex_id vl[2]; SInt32 ct=0; ep = p->e; vl[0] = hlist_get_key(edges,*idp,0); vl[1] = hlist_get_key(edges,*idp,1); kk = p->flags&3; id = p->links[kk%3]; while(id>=0) { q = (struct grad12_struct *)hlist_entry(edges,id); if((q->flags&8)==0) break; e = q->e; i = edge_verts_in_triangle(vl,r12(&e)); @q+9 if(i<0) abort_message("count error2");@> ct += edge_orient(i)*q->count; id = q->links[i]; } if(id<0) @@; else p->flags |= 16; if(id>=0 || edge_is_critical(&ep)) f = NULL_ID; else f = r12(&ep); } @ @= { int k; simplex_id e; for(k=0;k<3;k++) { if(k==kk) continue; e = get_edge(f,k); @@; } } @ @= { int flag; SInt32 qid; simplex_id vll[2]; if(!edge_is_paired_down(&e) && ( (options&3)==0 || !edge_is_deadend(&e) )) { get_edge_vertices(e,vll); qid = hlist_find(edges,vll); @q+9 if(qid<0) abort_message("count error 2");@> q = (struct grad12_struct *)hlist_entry(edges,qid); if((q->flags & 32)==0) { list_read_insert(todo_list,&qid); q->flags |= 32; } q->flags &= 0xffffffef; // indicate count field perhaps can now be filled, // since another incoming triangle has been counted. } } @ @= { p->flags &= 0xffffffef; p->flags |= 8; list_read_delete(todo_list); if (id==-2) { // |sigma| is a coface of |ep| so count it ct += edge_orient(edge_in_triangle(ep,sigma)); } if(edge_is_critical(&ep)) p->count = ct; else p->count = -ct * edge_orient(p->flags&3); if(p->count>1000 || p->count<-1000) abort_message("large count"); } @ @= { simplex_id vl[2]; simplex_id e; int i; printf("\nbacktrack: "); id = *(SInt32 *)list_entry(todo_list,0); while(id>=0) { p = (struct grad12_struct *)hlist_entry(edges,id); ep = p->e; vl[0] = hlist_get_key(edges,id,0); vl[1] = hlist_get_key(edges,id,1); printf(" %d %d,",id,p->flags); if(p->flags & 4096) break; p->flags |= 4096; id = p->links[(p->flags&3)%3]; while(id>=0) { p = (struct grad12_struct *)hlist_entry(edges,id); if((p->flags&8)==0) break; e = p->e; i = edge_verts_in_triangle(vl,r12(&e)); id = p->links[i]; } } printf("\n"); if(p->flags & 4096) abort_message("cyclic gradient field"); } @ @= { simplex_id vl[2]; simplex_id e; SInt32 id; printf("\nbacktrack2: "); id = *(SInt32 *)list_entry(todo_list,0); while(id>=0) { p = (struct grad12_struct *)hlist_entry(edges,id); ep = p->e; vl[0] = hlist_get_key(edges,id,0); vl[1] = hlist_get_key(edges,id,1); printf(" %d %d,",id,p->flags); id = p->links[(p->flags&3)%3]; while(id>=0) { p = (struct grad12_struct *)hlist_entry(edges,id); if((p->flags&1024)==0) break; e = p->e; i = edge_verts_in_triangle(vl,r12(&e)); id = p->links[i]; } } printf("\n"); if(p->flags & 4096) abort_message("cyclic gradient field"); } @ An alternative version since something is screwy here. set |flags&256| if an edge is added to todo list set |flags&1024| if an edge has a valid count2 field set |flags&2048| if an edge cannot yet be given a valid count2 field @= { list *todo_list; SInt32 *idp; int error_check; int i; list_initialize(&todo_list,sizeof(SInt32)); // first push all critical edges to todo_list todo = critid; while(todo>=0) { p = (struct grad12_struct *)hlist_entry(edges,todo); p->flags |= 256; list_push(todo_list,&todo); todo = p->links[1]; } // now push all ancestors of critical edges to todo_list for(i=0;i@; } // repeatedly read through, counting those we can, // and eliminating them from |todo_list| while(!list_is_empty(todo_list)) { list_read_init(todo_list); error_check = 1; while((idp = (SInt32 *)list_read(todo_list))!=NULL) { p = (struct grad12_struct *)hlist_entry(edges,*idp); if(p->flags&1024) {list_read_delete(todo_list); continue;} //if(p->flags&2048) continue; @@; } if(error_check && !list_is_empty(todo_list)) { @@; abort_message("count2 error"); } } list_abandon(&todo_list); } @ @= { SInt32 id; SInt32 ct=0; struct grad12_struct *q; simplex_id e; simplex_id vl[2]; ep = p->e; vl[0] = hlist_get_key(edges,*idp,0); vl[1] = hlist_get_key(edges,*idp,1); kk = p->flags&3; id = p->links[kk%3]; while(id>=0) { q = (struct grad12_struct *)hlist_entry(edges,id); if((q->flags&1024)==0) break; e = q->e; i = edge_verts_in_triangle(vl,r12(&e)); @q+9 if(i<0) abort_message("count error2");@> ct += edge_orient(i)*q->count2; id = q->links[i]; } if(id<0) @@; else p->flags |= 2048; if(id<0 && !edge_is_critical(&ep)) @@; } @ @= { p->flags &= 0xfffff7ff; p->flags |= 1024; list_read_delete(todo_list); if (id==-2) { // one incoming face is |sigma|, so count it too ct += edge_orient(edge_in_triangle(ep,sigma)); } if(edge_is_critical(&ep)) p->count2 = ct; else p->count2 = -ct * edge_orient(p->flags&3); error_check = 0; } @ @= { int k; f = r12(&ep); for(k=0;k<3;k++) { if(k==kk) continue; e = get_edge(f,k); @@; } } @ @= { SInt32 qid; simplex_id vll[2]; if(!edge_is_paired_down(&e)) { get_edge_vertices(e,vll); qid = hlist_find(edges,vll); if(qid>=0) { q = (struct grad12_struct *)hlist_entry(edges,qid); q->flags &= 0xfffff7ff; // indicate count2 field perhaps can now be filled, // since another incoming triangle has been count2ed. } else if( (options&3)==0 || !edge_is_deadend(&e) ) abort_message("count2 error 37"); } } @ @= { SInt32 id; struct grad12_struct *q; simplex_id vl[2],e; kk = p->flags&3; id = p->links[kk%3]; e = p->e; vl[0] = hlist_get_key(edges,todo,0); vl[1] = hlist_get_key(edges,todo,1); while(id>=0) { q = (struct grad12_struct *)hlist_entry(edges,id); ep = q->e; if((q->flags&256) == 0) { list_push(todo_list,&id); q->flags |= 256; } kk = edge_verts_in_triangle(vl,r12(&ep)); if(kk<0) abort_message("error"); id = q->links[kk]; } } @ Find all gradient paths ending at the edge tau. Now the hlist key is a triangle id but the entries are still grad12_structs. The list |crits| is a list of indices of items in |triangles| which are critical triangles. @= SInt32 find_all_backward_grad12_paths(simplex_id tau,hlist *triangles,list *crits,int options) { int kk = -1; simplex_id f; simplex_id ep; SInt32 todo = -1; SInt32 this; struct grad12_struct *p; SInt32 start = -1; list *todo_list; int tab[6] = {3,3,2,3,2,1}; hlist_clear(triangles); if(crits!=NULL) list_clear(crits); if((options&3)==2) list_initialize(&todo_list,sizeof(SInt32)); do { this = todo; if(this>=0) { p = (struct grad12_struct *)hlist_entry(triangles,this); f = hlist_get_key(triangles,this,0); ep = r21(&f); kk = p->flags&3; todo = p->links[kk]; p->links[kk] = -1; } else ep = tau; @@; } while(todo>=0); if((options&3)==2) @@; return start; } @ @= { simplex_id t,t1; simplex_id vlist[3]; int link_count = 0; t = t1 = id2tetra_id(ep); get_edge_vertices(ep,vlist); vlist[2] = id2tetrahedron(t)->vertex_ids[tab[edge_index(ep)]]; f = vertices2trianglex(vlist,t); do { if(id_is_null(f)) abort_message("backgrad error 3"); if(triangle_is_in_K(f) && (triangle_is_critical(&f) || ((!triangle_is_xdeadend(&f) || (options&3)==0 ) && triangle_is_paired_down(&f) && !edges_equal(ep,r21(&f))))) { link_count++; @@; } f = next_face(vlist,f); } while(t!=t1); if(link_count==0 && (options&3) && this>=0) { make_edge_xdeadend(&ep); f = r12(&ep); make_triangle_xdeadend(&f); } } @ @= { struct grad12_struct r,*q; SInt32 id; Boolean flag; int k; f = best_face_id(f); id = hlist_find_add(triangles,&f,&r,&flag); q = (struct grad12_struct *)hlist_entry(triangles,id); k = edge_in_triangle(ep,f); if(this>=0) p = (struct grad12_struct *)hlist_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(triangle_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; } } } @ @= { SInt32 id; int k; simplex_id e; struct grad12_struct r,*q; simplex_id f1; Boolean flag; while(!list_is_empty(todo_list)) { list_pop(todo_list,&id); p = (struct grad12_struct *)hlist_entry(triangles,id); if(p->flags&4) continue; p->flags|=4; kk = p->flags&3; f = hlist_get_key(triangles,id,0); for(k=0;k<3;k++) { if(k==kk || p->links[k] < -1) continue; e = get_edge(f,k); if(edge_is_paired_up(&e) && !edge_is_deadend(&e) ) { if(edge_is_xdeadend(&e)) abort_message("back grad error 3"); f1 = r12(&e); id = hlist_find_add(triangles,&f1,&r,&flag); q = (struct grad12_struct *)hlist_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;idflags & 4) continue; f = hlist_get_key(triangles,id,0); make_triangle_xdeadend(&f); e = r21(&f); make_edge_xdeadend(&e); } } @ @= { simplex_id ee; e = NULL_ID; while(!list_is_empty(crits)) { list_pop(crits,&m); q = (struct edge_graph*)hlist_entry(graph,m); ee = q->e; if(false) // |if(edge_is_persistent(&ee))| disable this for now { e = ee; break; } } }