f8-4

8.4 IndexingandRanking 329Sample 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 hpsort(n,ra) INTEGER n REAL ra(n) Sorts an array ra(1:n) into ascending numerical order using the Heapsort algorithm. nis input; rais replaced on output by its sorted rearrangement. INTEGER i,ir,j,l REAL rraif (n.lt.2) return The index lwill be decremented from its initial value down to 1 during the “hiring” (heap creation) phase. Once it reaches 1, the index irwill be decremented from its initial value down to 1 during the “retirement-and-promotion” (heap selection) phase. l=n/2+1ir=n 10 continue if(l.gt.1)then Still in hiring phase. l=l-1 rra=ra(l) else In retirement-and-promotion phase. rra=ra(ir) Clear a space at end of array. ra(ir)=ra(1) Retire the top of the heap into it. ir=ir-1 Decrease the size of the corporation. if(ir.eq.1)then Done with the last promotion. ra(1)=rra The least competent worker of all! return endif endifi=l Whether in the hiring phase or promotion phase, we here set up to sift down element rra to its proper level. j=l+l 20 if(j.le.ir)then “Do while j.le.ir :” if(j.lt.ir)then if(ra(j).lt.ra(j+1))j=j+1 Compare to the better underling. endifif(rra.lt.ra(j))then Demote rra. ra(i)=ra(j) i=j j=j+j else This is rra’s level. Set jto terminate the sift-down. j=ir+1 endif goto 20endif ra(i)=rra Put rra into its slot. goto 10END CITED REFERENCES AND FURTHER READING: Knuth,D.E. 1973, SortingandSearching ,v ol.3of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), §5.2.3. [1] Sedgewick, R. 1988, Algorithms , 2nd ed. (Reading, MA: Addison-Wesley), Chapter 11. [2] 8.4 Indexing and Ranking The conceptof keysplays a prominentrole in the managementof data files. A datarecordin sucha file maycontainseveralitems, or fields. Forexample,a record in a file of weather observations may have fields recording time, temperature, and 330 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).15 63 57 432 38 2114 3 66 51 42 34 215 5 61 52 46 33 214 32 615 514 48 37 213original arrayindex tablerank tablesorted array (a) (b) (c) (d) Figure 8.4.1. (a) An unsorted array of six numbers. (b) Index table, whose entries are pointers to the elements of (a) in ascending order. (c) Rank table, whose entries are the ranks of the corresponding elements of (a). (d) Sorted array of the elements in (a). wind velocity. When we sort the records, we must decide which of these fields we want to be brought into sorted order. The other fields in a record just come along for the ride, and will not, in general, end up in any particular order. The field on which the sort is performed is called the keyfield. For a data file with many records and many fields, the actual movement of N recordsintothesortedorderoftheirkeys K i,i=1 ,...,N,canbea dauntingtask. Instead, one can construct an index table Ij,j =1 ,...,N, such that the smallest Kihas i=I1, the second smallest has i=I2, and so on up to the largest Kiwith i=IN. In other words, the array KIj j=1 ,2,...,N (8.4.1 ) is insortedorderwhenindexedby j. Whenanindextableisavailable,oneneednot move records from their original order. Further, different index tables can be made from the same set of records, indexing them to different keys. The algorithm for constructing an index table is straightforward: Initialize the index array with the integers from 1 to N, then perform the Quicksort algorithm, movingtheelementsaround asifoneweresortingthekeys. Theintegerthatinitially numberedthe smallest key thus ends up in the numberone position, and so on. SUBROUTINE indexx(n,arr,indx) INTEGER n,indx(n),M,NSTACK REAL arr(n) PARAMETER (M=7,NSTACK=50) Indexes an array arr(1:n) , i.e., outputs the array indx(1:n) such that arr(indx(j)) is in ascending order for j=1 ,2,...,N . The input quantities nandarr are not changed. 8.4 IndexingandRanking 331Sample 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).INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK) REAL a do11j=1,n indx(j)=j enddo 11 jstack=0l=1ir=n 1 if(ir-l.lt.M)then do 13j=l+1,ir indxt=indx(j)a=arr(indxt)do 12i=j-1,l,-1 if(arr(indx(i)).le.a)goto 2 indx(i+1)=indx(i) enddo 12 i=l-1 2 indx(i+1)=indxt enddo 13 if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1)jstack=jstack-2 else k=(l+ir)/2itemp=indx(k)indx(k)=indx(l+1) indx(l+1)=itemp if(arr(indx(l)).gt.arr(indx(ir)))then itemp=indx(l) indx(l)=indx(ir) indx(ir)=itemp endifif(arr(indx(l+1)).gt.arr(indx(ir)))then itemp=indx(l+1) indx(l+1)=indx(ir)indx(ir)=itemp endif if(arr(indx(l)).gt.arr(indx(l+1)))then itemp=indx(l)indx(l)=indx(l+1) indx(l+1)=itemp endifi=l+1j=ir indxt=indx(l+1) a=arr(indxt) 3 continue i=i+1 if(arr(indx(i)).lt.a)goto 3 4 continue j=j-1 if(arr(indx(j)).gt.a)goto 4 if(j.lt.i)goto 5itemp=indx(i) indx(i)=indx(j) indx(j)=itempgoto 3 5 indx(l+1)=indx(j) indx(j)=indxt jstack=jstack+2if(jstack.gt.NSTACK)pause ’NSTACK too small in indexx’if(ir-i+1.ge.j-l)then istack(jstack)=ir 332 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).istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1istack(jstack-1)=l l=i endif endifgoto 1 END If you want to sort an array while making the correspondingrearrangementof several or many other arrays, you should first make an index table, then use it to rearrange each array in turn. This requires two arrays of working space: one to hold the index, and another into which an array is temporarily moved, and fromwhich it is redeposited back on itself in the rearranged order. For 3 arrays, the procedure looks like this: SUBROUTINE sort3(n,ra,rb,rc,wksp,iwksp) INTEGER n,iwksp(n)REAL ra(n),rb(n),rc(n),wksp(n) C USES indexx Sorts an array ra(1:n) into ascending numerical order while making the corresponding rearrangements of the arrays rb(1:n) andrc(1:n) . An index table is constructed via the routine indexx . INTEGER jcall indexx(n,ra,iwksp) Make the index table. do 11j=1,n Save the array ra. wksp(j)=ra(j) enddo 11 do12j=1,n Copy it back in the rearranged order. ra(j)=wksp(iwksp(j)) enddo 12 do13j=1,n Ditto rb. wksp(j)=rb(j) enddo 13 do14j=1,n rb(j)=wksp(iwksp(j)) enddo 14 do15j=1,n Ditto rc. wksp(j)=rc(j) enddo 15 do16j=1,n rc(j)=wksp(iwksp(j)) enddo 16 return END The generalizationto any other numberof arrays is obviouslystraightforward. Aranktable is differentfroman indextable. A ranktable ’sjth entrygivesthe rankof the jth element ofthe originalarrayof keys,rangingfrom1 (if that element was the smallest) to N(if that element was the largest). One can easily construct a rank table from an index table, however: 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) as output from the routine indexx , this routine returns an array 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,statisticsbooksde finethemedianasthearithmetic 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 de finek=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: