f8-5

8.5SelectingtheMthLargest 333Sample 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 rank(n,indx,irank) INTEGER n,indx(n),irank(n) Given indx(1:n) asoutput fromtheroutine indexx,thisroutinereturns anarray irank(1:n) , the corresponding table of ranks. INTEGER j do11j=1,n irank(indx(j))=j enddo 11 return END Figure 8.4.1 summarizes the concepts discussed in this section. 8.5 Selecting the Mth Largest Selectionissorting’sausteresister. (Say thatfivetimesquickly!) Wheresorting demandstherearrangementofanentiredataarray,selectionpolitelyasksforasinglereturnedvalue: Whatisthe kthsmallest(or,equivalently,the m=N+1−kthlargest) element out of Nelements? The fastest methods for selection do, unfortunately, rearrangethearrayfortheirowncomputationalpurposes,typicallyputtingallsmaller elements to the left of the kth, all larger elements to the right, and scrambling the order within each subset. This side effect is at best innocuous, at worst downrightinconvenient. Whenthearrayisverylong,sothatmakingascratchcopyofitistaxing on memory, or when the computational burden of the selection is a negligible part of a larger calculation, one turns to selection algorithms without side effects, whichleavetheoriginalarrayundisturbed. Such inplaceselectionisslowerthanthefaster selection methodsbya factor ofabout10. We giveroutinesofboth types,below. The most common use of selection is in the statistical characterization of a set of data. One often wants to know the median element in an array, or the top and bottom quartile elements. When Nis odd, the median is the kth element, with k=( N+1 ) /2. When Niseven,statisticsbooksdefinethemedianasthearithmetic mean of the elements k=N/2and k=N/2+1(that is, N/2from the bottom and N/2fromthetop). Ifyouacceptsuchpedantry,youmust performtwo separate selections to find these elements. For N> 100we usually define k=N/2to be the median element, pedants be damned. The fastest general method for selection, allowing rearrangement,is partition- ing, exactly as was done in the Quicksort algorithm ( §8.2). Selecting a “random” partition element, one marches through the array, forcing smaller elements to theleft, larger elements to the right. As in Quicksort, it is important to optimize the inner loop, using “sentinels” ( §8.2) to minimize the number of comparisons. For sorting, one would then proceed to further partition both subsets. For selection,we can ignore one subset and attend only to the one that contains our desired kth element. Selectionby partitioningthusdoes notneeda stackofpendingoperations, and its operations count scales as Nrather than as Nlog N(see [1]). Comparison with sortin§8.2 should make the following routine obvious: 334 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).FUNCTION select(k,n,arr) INTEGER k,n REAL select,arr(n) Returns the kth smallest value in the array arr(1:n) . The input array will be rearranged to have this value in location arr(k), with all smaller elements moved to arr(1:k-1) (in arbitrary order) and all larger elements in arr[k+1..n] (also in arbitrary order). INTEGER i,ir,j,l,midREAL a,templ=1 ir=n 1 if(ir-l.le.1)then Active partition contains 1 or 2elements. if(ir-l.eq.1)then Active partition contains 2elements. if(arr(ir).lt.arr(l))then temp=arr(l) arr(l)=arr(ir)arr(ir)=temp endif endifselect=arr(k)return else mid=(l+ir)/2 Choose median of left, center, and right elements as par- titioning element a. Also rearrange so that arr(l) ≤ arr(l+1) ,arr(ir) ≥arr(l+1) .temp=arr(mid) arr(mid)=arr(l+1) arr(l+1)=tempif(arr(l).gt.arr(ir))then temp=arr(l) arr(l)=arr(ir) arr(ir)=temp endif if(arr(l+1).gt.arr(ir))then temp=arr(l+1)arr(l+1)=arr(ir)arr(ir)=temp endif if(arr(l).gt.arr(l+1))then temp=arr(l)arr(l)=arr(l+1) arr(l+1)=temp endifi=l+1 Initialize pointers for partitioning. j=ir a=arr(l+1) Partitioning element. 3 continue Beginning of innermost loop. i=i+1 Scan up to find element >a. if(arr(i).lt.a)goto 3 4 continue j=j-1 Scan down to find element <a. if(arr(j).gt.a)goto 4 if(j.lt.i)goto 5 Pointers crossed. Exit with partitioning complete. temp=arr(i) Exchange elements. arr(i)=arr(j) arr(j)=temp goto 3 End of innermost loop. 5 arr(l+1)=arr(j) Insert partitioning element. arr(j)=a if(j.ge.k)ir=j-1 Keep active the partition that contains the kth element. if(j.le.k)l=i endif goto 1 END 8.5SelectingtheMthLargest 335Sample 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).In-place, nondestructive, selection is conceptually simple, but it requires a lot of bookkeeping,and it is correspondinglyslower. The general idea is to pick somenumber Mof elements at random, to sort them, and then to make a pass through the array counting how many elements fall in each of the M+1intervals defined by these elements. The kth largest will fall in one such interval — call it the “live” interval. Onethendoesasecondround,first picking Mrandomelementsinthelive interval,andthendeterminingwhichof thenew,finer, M+1intervalsall presently live elements fall into. And so on, until the kth element is finally localized within a single array of size M, at which point direct selection is possible. How shall we pick M? The number of rounds, log MN=l o g2N/ log2M, will be smaller if Mis larger; but the work to locate each element among M+1 subintervals will be larger, scaling as log2Mfor bisection, say. Each round requires looking at all Nelements, if only to find those that are still alive, while the bisections are dominated by the Nthat occur in the first round. Minimizing O(NlogMN)+ O(Nlog2M)thus yields the result M∼2√ log2N(8.5.1 ) Thesquarerootofthelogarithmissoslowlyvaryingthatsecondaryconsiderationsof machinetimingbecomeimportant. We use M=6 4as aconvenientconstantvalue. Twominoradditionaltricksinthefollowingroutine, selip,are(i)augmenting the set of Mrandom values by an M+1st, the arithmetic mean, and (ii) choosing the Mrandomvalues“onthefly”inapassthroughthedata,byamethodthatmakes later values no less likely to be chosen than earlier ones. (The underlyingidea is to giveelement m>ManM/mchanceof beingbroughtintothe set. Youcan prove by induction that this yields the desired result.) FUNCTION selip(k,n,arr) INTEGER k,n,M REAL selip,arr(n),BIG PARAMETER (M=64,BIG=1.E30) Returns the kth smallest value in the array arr(1:n) . The input array is not altered. C USES shell INTEGER i,j,jl,jm,ju,kk,mm,nlo,nxtmm,isel(M+2)REAL ahi,alo,sum,sel(M+2)if(k.lt.1.or.k.gt.n.or.n.le.0) pause ’bad input to selip’kk=k ahi=BIG alo=-BIG 1 continue Mainiterationloop, untildesiredelement isisolated. mm=0 nlo=0sum=0.nxtmm=M+1 do 11i=1,n Make a pass through the whole array. if(arr(i).ge.alo.and.arr(i).le.ahi)then Consideronlyelements inthecur- rent brackets. mm=mm+1 if(arr(i).eq.alo) nlo=nlo+1 In case of ties for low bracket. if(mm.le.M)then Statistical procedure forselecting min-range elements with equal probability, even without knowing inadvance how many there are!sel(mm)=arr(i) else if(mm.eq.nxtmm)then nxtmm=mm+mm/M sel(1+mod(i+mm+kk,M))=arr(i) Themodfunction providesasome- what random number. endif sum=sum+arr(i) 336 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).endif enddo 11 if(kk.le.nlo)then Desired element is tied for lower bound; return it. selip=aloreturn else if(mm.le.M)then All in-range elements were kept. So return answer by direct method. call shell(mm,sel) selip=sel(kk)return endif Augment selected set by mean value (fixes degenera- cies), and sort it. sel(M+1)=sum/mm call shell(M+1,sel)sel(M+2)=ahi do 12j=1,M+2 Zero the count array. isel(j)=0 enddo 12 do13i=1,n Make another pass through the whole array. if(arr(i).ge.alo.and.arr(i).le.ahi)then For each in-range element.. jl=0ju=M+2 2 if(ju-jl.gt.1)then ...find its position among the select by bisection... jm=(ju+jl)/2if(arr(i).ge.sel(jm))then jl=jm else ju=jm endif goto 2 endifisel(ju)=isel(ju)+1 ...and increment the counter. endif enddo 13 j=1 Now we can narrow the bounds to just one bin, that is, by a factor of order m. 3 if(kk.gt.isel(j))then alo=sel(j) kk=kk-isel(j)j=j+1 goto 3 endif ahi=sel(j) goto 1 END Approximate timings: selipis about 10 times slower than select. Indeed, for Nin the range of ∼105,selipis about 1.5 times slower than a full sort with sort, while selectis about 6 times faster than sort. You should weigh time against memory and convenience carefully. Of course neither of the above routines should be used for the trivial cases of finding the largest, or smallest, element in an array. Those cases, you code by handas simple doloops. Thereare also goodways to codethe case where kis modestin comparison to N, so that extra memory of order kis not burdensome. An example is to use the method of Heapsort ( §8.3) to make a single pass through an array of length Nwhile saving the mlargestelements. The advantage of the heap structure isthatonly log m,ratherthan m,comparisonsarerequiredeverytimeanewelement is added to the candidate list. This becomes a real savings when m>O (√ N),b u t it neverhurts otherwiseandis easy tocode. Thefollowingprogramgivesthe idea. SUBROUTINE hpsel(m,n,arr,heap) INTEGER m,n 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 Returns in heap(1:m) the largest melements of the array arr(1:n) ,w i t h heap(1) guar- anteed to be the the mth largest element. The array arris not altered. For efficiency, this routine should be used only when 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) Create initial heap by overkill! We assume m/lessmuchn. do12i=m+1,n For each remaining element... if(arr(i).gt.heap(1))then Put it on the heap? heap(1)=arr(i)j=1 1 continue Sift down. 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