f8-6

8.6DeterminationofEquivalenceClasses 337Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).REAL arr(n),heap(m) C USES sort Returnsin heap(1:m)thelargest melementsofthearray arr(1:n),with heap(1)guar- anteedtobethethe mthlargestelement. Thearray arrisnotaltered. Forefficiency,this routineshouldbeusedonlywhen m/lessmuchn. INTEGER i,j,k REAL swapif (m.gt.n/2.or.m.lt.1) pause ’probable misuse of hpsel’do 11i=1,m heap(i)=arr(i) enddo 11 call sort(m,heap) Createinitialheapbyoverkill! Weassume m/lessmuchn. do12i=m+1,n Foreachremainingelement... if(arr(i).gt.heap(1))then Putitontheheap? heap(1)=arr(i)j=1 1 continue Siftdown. k=2*jif(k.gt.m)goto 2if(k.ne.m)then if(heap(k).gt.heap(k+1))k=k+1 endifif(heap(j).le.heap(k))goto 2 swap=heap(k) heap(k)=heap(j)heap(j)=swapj=k goto 1 2 continue endif enddo 12 return end CITED REFERENCES AND FURTHER READING: Sedgewick, R. 1988, Algorithms , 2nd ed. (Reading, MA: Addison-Wesley), pp. 126ff. [1] Knuth,D.E. 1973, SortingandSearching ,v ol.3of TheArtofComputerProgramming (Reading, MA: Addison-Wesley). 8.6 Determination of Equivalence Classes A number of techniques for sorting and searching relate to data structures whose details are beyond the scope of this book, for example, trees, linked lists, etc. These structures andtheir manipulations are the bread and butter of computer science, as distinct from numericalanalysis, and there is no shortage of books on the subject. Inworkingwithexperimentaldata,wehavefoundthatoneparticularsuchmanipulation, namely the determination of equivalence classes, arises sufficiently often to justify inclusionhere. The problem is this: There are N“elements” (or “data points” or whatever), numbered 1,...,N. You are given pairwise information about whether elements are in the same equivalenceclass of“sameness,”bywhatevercriterionhappenstobeofinterest. Forexample, you may have a list of facts like: “Element 3 and element 7 are in the same class; element19 and element 4 are in the same class; element 7 and element 12 are in the same class, ....” Alternatively, you may have a procedure, given the numbers of two elements jand k, for 338 Chapter8. SortingSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).deciding whether they are in the same class or different classes. (Recall that an equivalence relation can be anything satisfying the RST properties : reflexive, symmetric, transitive. This is compatible with any intuitive definition of “sameness.”) The desired output is an assignment to each of the Nelements of an equivalence class number, such that two elements are in the same class if and only if they are assigned thesame class number. Efficientalgorithmsworklikethis: Let F(j)betheclassor“family”numberofelement j. Start off with each element in its own family, so that F(j)= j. The array F(j)can be interpretedasatreestructure,where F(j)denotestheparentof j. Ifwearrangeforeachfamily to be its own tree, disjoint from all the other “family trees,” then we can label each family(equivalence class) by its most senior great-great- ...grandparent. The detailed topology of the tree doesn’t matter at all, as long as we graft each related element onto it somewhere . Therefore, we process each elemental datum “ jis equivalent to k” by (i) tracking j up to its highest ancestor, (ii) tracking kup to its highest ancestor, (iii) giving jtokas a new parent, or vice versa (it makes no difference). After processing all the relations, we go through all the elements jand reset their F(j)’s to their highest possible ancestors, which then label the equivalence classes. The following routine, based on Knuth [1], assumes that there are melemental pieces of information, stored in two arrays of length m,lista,listb , the interpretation being thatlista(j) andlistb(j) ,j=1...m, are the numbers of two elements which (we are thus told) are related. SUBROUTINE eclass(nf,n,lista,listb,m) INTEGER m,n,lista(m),listb(m),nf(n) Given mequivalencesbetweenpairsof nindividualelementsintheformoftheinputarrays lista(1:m)andlistb(1:m),thisroutinereturnsin nf(1:n)thenumberoftheequiva- lenceclassofeachofthe nelements,integersbetween 1andn(notallsuchintegersused). INTEGER j,k,l do11k=1,n Initializeeachelementitsownclass. nf(k)=k enddo 11 do12l=1,m Foreachpieceofinputinformation... j=lista(l) 1 if(nf(j).ne.j)then Trackfirstelementuptoitsancestor. j=nf(j) goto 1endifk=listb(l) 2 if(nf(k).ne.k)then Tracksecondelementuptoitsancestor. k=nf(k) goto 2 endif if(j.ne.k)nf(j)=k Iftheyarenotalreadyrelated,makethemso. enddo 12 do13j=1,n Finalsweepuptohighestancestors. 3 if(nf(j).ne.nf(nf(j)))then nf(j)=nf(nf(j)) goto 3 endif enddo 13 returnEND Alternatively, wemay be able toconstruct aprocedure equiv(j,k) thatreturns a value .true.if elements jandkare related, or .false. if they are not. Then we want to loop over all pairs of elements to get the complete picture. D. Eardley has devised a clever way ofdoing thiswhilesimultaneously sweeping thetreeuptohigh ancestors inamanner thatkeepsit current and obviates most of the final sweep phase: 8.6Determinationof EquivalenceClasses 339Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE eclazz(nf,n,equiv) INTEGER n,nf(n) LOGICAL equiv EXTERNAL equiv Givenauser-suppliedlogicalfunction equivwhichtellswhetherapairofelements,each intherange 1...n,arerelated,returnin nfequivalenceclassnumbersforeachelement. INTEGER jj,kknf(1)=1do 12jj=2,n Loopoverfirstelementofallpairs. nf(jj)=jj do11kk=1,jj-1 Loopoversecondelementofallpairs. nf(kk)=nf(nf(kk)) Sweepitupthismuch. if (equiv(jj,kk)) nf(nf(nf(kk)))=jj Goodexerciseforthereadertofigure out why this much ancestry is necessary!enddo 11 enddo 12 do13jj=1,n Onlythismuchsweepingisneededfinally. nf(jj)=nf(nf(jj)) enddo 13 returnEND CITED REFERENCES AND FURTHER READING: Knuth,D.E.1968, FundamentalAlgorithms ,vol.1ofTheArtofComputerProgramming (Reading, MA: Addison-Wesley),§2.3.3. [1] Sedgewick, R. 1988, Algorithms , 2nd ed. (Reading, MA: Addison-Wesley), Chapter 30.