! ************************************** ! * Module phil2 * ! * R. J. Purser NOAA/NCEP/EMC 2017 * ! * jim.purser@noaa.gov * ! ************************************** ! ! Module procedures pertaining to the project, sorting, and density ! estimation by B-spline smoothing the index delta functions of projected ! locations. The particular application that makes this approach potentially ! advantageous compared to conventional alternatives is the estimation ! of the local spatial density of data in those cases where the data tend to be ! inhomogeneously clumped, such as aircraft data, which are typically ! clustered around densely populated areas and along the heavily-trafficked ! flight tracks. ! ! Direct dependencies ! Libraries: psort,pmat ! Modules: phil,psort,peuc,kinds, pietc ! !============================================================================= module phil2 !============================================================================= ! B-spline smoothers of degree 0 and 1 are provided by module procedures ! whose interface names are bsmoo0 and bsmoo1 respectively. The versions ! bsmoo0s and bsmoo1s are special forms that assume that the data locations ! have already been presorted and arranged consecutively; otherwise, the ! versions bsmoo0 and bsmoo1 include an integer index array argument, rank, that ! allows the ordered data to be accessed by indirect addressing. ! The routine, denest, performs such volumetric density estimations in a thin ! spherical shell nrand times and avergaes the results, where nrand can be ! any integerup to 30, or 100, or 300, and uses precomputed pseudo-random ! orientations of the cubic framework within which the hilbert curve is ! constructed over the sphere. The vertical randomization is performed by ! steps that are regular powers of (1/3) times one quarter of the active ! thickness (containing all the data) of the shell. ! The precomputed random rotations are contained in the ascii file, qsets.asc, ! encoded as quaternions. !============================================================================= use kinds, only: dp,i_kind implicit none private public:: bsmoo0,bsmoo1,denest,denest2d,& getqset5,getqset7,getqset8,getqset13 ! <- temporarily public interface bsmoo0; module procedure bsmoo0,bsmoo0s; end interface interface bsmoo1; module procedure bsmoo1,bsmoo1s; end interface interface denest; module procedure denest,denestx; end interface interface denest2d; module procedure denest2d,denest2dx; end interface interface getqset5; module procedure getqset5; end interface interface getqset7; module procedure getqset7; end interface interface getqset8; module procedure getqset8; end interface interface getqset13;module procedure getqset13; end interface contains !============================================================================= subroutine bsmoo0(nob,span,sob,rank,dob)! [bsmoo0] !============================================================================= ! Perform a smoothing of an irregular one-dimensional distribution of nob unit ! impulses at parameter locations, sob, using a boxcar function of total width, ! span. The resulting 'density' at each datum location is dob. Note that this ! is the trivial instance of B-spline smoothing where the degree of the spline ! is zero. !============================================================================= implicit none integer(i_kind), intent(in ):: nob real(dp), intent(in ):: span real(dp), dimension(nob),intent(in ):: sob integer(i_kind) ,dimension(nob),intent(in ):: rank real(dp), dimension(nob),intent(inout):: dob !----------------------------------------------------------------------------- real(dp):: spanh,st,s1,s2 integer(i_kind) :: i,j,i1,i2,L1,L2 !============================================================================= spanh=span/2.0_dp i1=1 i2=1 do i=1,nob j=rank(i) st=sob(j); s1=st-spanh; s2=st+spanh L1=i1; L2=i2 do i1=l1,nob; if(sob(rank(i1))>s1)exit; enddo do i2=l2,nob; if(sob(rank(i2))>s2)exit; enddo dob(j)=dob(j)+(i2-i1)/span enddo end subroutine bsmoo0 !============================================================================= subroutine bsmoo0s(nob,span,sob,dob)! [bsmoo0] !============================================================================= implicit none integer(i_kind), intent(in ):: nob real(dp), intent(in ):: span real(dp),dimension(nob),intent(in ):: sob real(dp),dimension(nob),intent(inout):: dob !----------------------------------------------------------------------------- real(dp):: spanh,st,s1,s2 integer(i_kind) :: i,i1,i2,L1,L2 !============================================================================= spanh=span/2.0_dp i1=1 i2=1 do i=1,nob st=sob(i); s1=st-spanh; s2=st+spanh L1=i1; L2=i2 do i1=l1,nob; if(sob(i1)>s1)exit; enddo do i2=l2,nob; if(sob(i2)>s2)exit; enddo dob(i)=dob(i)+(i2-i1)/span enddo end subroutine bsmoo0s !============================================================================= subroutine bsmoo1(nob,span,sob,rank,dob)! [bsmoo1] !============================================================================= ! Perform a smoothing of an irregular one-dimensional distribution of nob unit ! impulses at parameter locations, sob, using a hat function of half-width, ! span. The resulting 'density' at each datum location is dob. Note that this ! is the instance of a B-spline smoothing where the degree of the spline ! is one. ! ! The general idea of this algorithm is that a piecewise linear spline p(s) ! continuously links nodes at each sb=sob(rank(ib)), ! where the slope changes by +1. ! The distribution of delta functions smoothed by the B-spline of half-width ! span leads to a smooth estimate, dob, at each location, s=sb, given by ! dob(ib) = [p(sb-span) -2*p(sb) +p(sb+span)]/span**2 ! ! However, because this formula, as stated, would tend to involve a "small" ! result coming from differences between approximately equally "large" numbers ! when the count, nob, of data is very large, it is preferable to ! transform the algorithm into a better-conditioned equivalent one in order ! to keep the round-off errors small. We do this by exploiting the fact that ! the result is formally unchanged by adding ANY linear polynomial uniformly ! to the whole spline function, p(s). We choose this polynomial, at each ! stage, to maintain the adjusted spline function p just left of each targetted ! datum, ib, the zero function, and therefore the spline function just right ! of it equal to p(s) = (s-sb), where sb denotes the location of datum ib. ! This keeps the pairs of numerical polynomial coefficients at sa=sb-span and at ! sc=sb+span relatively small. The "slope" coefficient of each linear ! polynomial is always integral, and is therefore not subject to any round-off ! error, but the floating-point constant coefficients, pa0 and pc0, are ! in principle subject to cumulative round-off in cases where the number, nob, ! of data is very large. As a safeguard, we shrink the new pa0 and pc0 values ! towards zero by a very tiny amount at each new datum, but by a degree that ! should be enough to counteract any tendency for cumulative round-off to cause ! their values to spuriously diverge. !============================================================================= use pietc, only: u1 use kinds, only: dp,i_kind implicit none integer(i_kind), intent(in ):: nob real(dp), intent(in ):: span real(dp), dimension(nob),intent(in ):: sob integer(i_kind), dimension(nob),intent(in ):: rank real(dp), dimension(nob),intent(inout):: dob !----------------------------------------------------------------------------- real(dp),parameter:: eps=1.0e-13_dp,shrink=1.0_dp-eps real(dp) :: sa,sb,sc,pa0,pc0,sbold,ds,soba,sobc,spansi integer(i_kind) :: ia,ib,ic,jb,La,Lc,pa1,pb1,pc1 !============================================================================= spansi=u1/span**2 ia=1 ic=1 sbold=0 ! Initially, all the polynomial coefficients vanish: pa0=0.0_dp pa1=0 pb1=0 pc0=0.0_dp pc1=0 do ib=1,nob jb=rank(ib) sb=sob(jb); sa=sb-span; sc=sb+span ds=sb-sbold; sbold=sb ! Move the new coordinate origin to location, sb=sob(rank(ib)) and ! subtract from polynomials A and C the (old) left-polynomial at b: pa1=pa1-pb1 pc1=pc1-pb1 pa0=(pa0+ds*pa1)*shrink ! Shrink by imperceptible amount to stabilize pc0=(pc0+ds*pc1)*shrink ! huge computations against cumulative round-off ! New pb0 is always taken implicitly to be zero and new pb1 becomes simply: pb1=1! < (new) rt-polynomial p(s) at B is trivially p=(s-sob(rank(ib))) ! Update the new locations of A and C (straddling B by distance, span) and ! the corresponding polynomial coefficients, {pa0,pa1} and {pc0,pc1} valid ! there: La=ia; Lc=ic do ia=La,nob soba=sob(rank(ia)); if(soba>sa)exit; pa0=pa0+sb-soba; pa1=pa1+1 enddo do ic=Lc,nob sobc=sob(rank(ic)); if(sobc>sc)exit; pc0=pc0+sb-sobc; pc1=pc1+1 enddo ! The formula, [p(sb-span) -2*p(sb) +p(sb+span)], simplies to: dob(jb)=dob(jb)+(pa0+pc0+span*(pc1-pa1))*spansi enddo end subroutine bsmoo1 !============================================================================= subroutine bsmoo1s(nob,span,sob,dob)! [bsmoo1] !============================================================================= use pietc, only: u1 use kinds, only: dp,i_kind implicit none integer(i_kind), intent(in ):: nob real(dp), intent(in ):: span real(dp),dimension(nob),intent(in ):: sob real(dp),dimension(nob),intent(inout):: dob !----------------------------------------------------------------------------- real(dp),parameter:: eps=1.e-13_dp,shrink=1.0_dp-eps real(dp) :: sa,sb,sc,pa0,pc0,sbold,ds,spansi integer(i_kind) :: ia,ib,ic,La,Lc,pa1,pb1,pc1 !============================================================================= spansi=u1/span**2 ia=1 ic=1 sbold=0.0_dp ! Initially, all the polynomial coefficients vanish: pa0=0.0_dp pa1=0 pb1=0 pc0=0.0_dp pc1=0 do ib=1,nob sb=sob(ib); sa=sb-span; sc=sb+span ds=sb-sbold; sbold=sb ! Move the new coordinate origin to location, sb=sob(ib) and ! subtract from polynomials A and C the (old) left-polynomial at b: pa1=pa1-pb1 pc1=pc1-pb1 pa0=(pa0+ds*pa1)*shrink ! Shrink by imperceptible amount to stabilize pc0=(pc0+ds*pc1)*shrink ! huge computations against cumulative round-off ! New pb0 is always taken implicitly to be zero and new pb1 becomes simply: pb1=1! < (new) rt-polynomial p(s) at B is trivially p=(s-sob(ib)) ! Update the new locations of A and C (straddling B by distance, span) and ! the corresponding polynomial coefficients, {pa0,pa1} and {pc0,pc1} valid ! there: La=ia; Lc=ic do ia=La,nob; if(sob(ia)>sa)exit; pa0=pa0+sb-sob(ia); pa1=pa1+1; enddo do ic=Lc,nob; if(sob(ic)>sc)exit; pc0=pc0+sb-sob(ic); pc1=pc1+1; enddo ! The formula, [p(sb-span) -2*p(sb) +p(sb+span)], simplies to: dob(ib)=dob(ib)+(pa0+pc0+span*(pc1-pa1))*spansi enddo end subroutine bsmoo1s !============================================================================= subroutine denest(nob,nrand,nor, &! [denest] re,dentrip,hscale,vscale,vmin,vmax,& latob,lonob,altob,wtob) !============================================================================= ! Use nrand randomized Hilbert curves over a sufficiently thick spherical ! shell to estimate the density of data relative to the "unit" volume ! implied by the vertical smoothing scale time the square of the horizontal ! smoothing scale. ! ! nob: the number of data ! nrand: the number of distinct pseudo-random Hilbert curves ! nor: the order of the B-spline smoother (either 0 or 1) ! re: the nominal radius of the earth ! dentrip:density criterion tripping the downward adjustment to the wt factor ! hscale: (in the same units) the characteristic horizontal scale of averaging ! vscale: (in the same units) the characteristic vertical scale of averaging ! vmin: (in same units) lowest altitude bound for valid data ! vmax: (in same units) greatest altitude bound for valid data ! latob: array of observation latitudes (degrees north) ! lonob: array of observation longitudes (degrees east) ! altob: array of observation altitudes (same units as other distance scales) ! troflg: logic parameter to determin wether dentrip in tropic different from ! other region ! wtob: Weight factor implied by the estimated data density and dentrip ! ! Internally, it is convenient to relate all distance units to "Hilbert ! parameter units" where 24 of these units fill the 24 fundamental ! quadrilateral panels that form the horizontal footprint ! of the spherical shell. This means that one Hilbert unit corresponds to ! a horizontal distance of sqrt(pi/6)*re, or about Uh=4610 km. The thickness, ! in Hilbert distance, of the valid vertical range is ! (vmax-vmin)*hscale/(uh*vscale), which must be less than one for the ! present "thin shell" Hilbert curve construction. We inflate this thickness ! by a small margin to leave room for vertical location randomization of ! the Hilbert curves, and round up to the next integer(i_kind) power, ngen4, of 1/2. ! Ngen4 effectively determines the generation at which the Hilbert curve ! transitions from being 2-dimensional (at coarse scales) to being ! 3-dimensional (at and below the scale of the effective thickness of the ! shell that the curve fills). For generations, 1--ngen4, the expansion is ! at base-4; at ngen4--ngen the expansion is base-8. !============================================================================= use pietc, only: u0,u1,u2,o2,pi,dtor use kinds, only: dp,i_kind use peuc, only: qtorot,mulqq use phil0, only: xs_to_hil48,hil8_to_r,hil4_to_r use psort, only: sort, invertperm use qcmod, only:troflg,lat_c implicit none integer(i_kind), intent(in ):: nob,nrand,nor real(dp), intent(in ):: re,dentrip,hscale,vscale,vmin,vmax real(dp),dimension(nob),intent(in ):: latob,lonob,altob real(dp),dimension(nob),intent(out):: wtob !----------------------------------------------------------------------------- integer(i_kind),parameter :: ngen=14 real(dp), dimension(3,nob) :: xob real(dp), dimension(nob) :: rob,sob,latobc real(dp), dimension(0:3) :: qset real(dp), dimension(0:3,3) :: qset3 real(dp), dimension(0:3,5) :: qset5 real(dp), dimension(0:3,7) :: qset7 real(dp), dimension(0:3,13) :: qset13 real(dp), dimension(3,3) :: rot,rotnew,rotold real(dp) :: span,uv,uh,vh,vhq,vhp,v,vrand,& rlat,clat,slat,rlon,clon,slon integer(i_kind),dimension(nob) :: rank integer(i_kind),dimension(0:ngen) :: hilr integer(i_kind) :: i,irand,j,k,L,ngen4,ntri,dentripc !============================================================================= uh=sqrt(pi/6.0_dp)*re uv=(uh*vscale)/hscale vh=(vmax-vmin)/uv vhq=vh/4.0_dp vhp=vh+vhq if(vhp>=u1)then print& '("In denest; vmax-vmin too large for thin-shell assumption to be valid")' print'("Make this vertical range smaller, or vscale/hscale larger")' stop endif v=o2 do ngen4=0,ngen-1 if(v=nrand)exit enddo vrand=vhq/i !vrand=3*vhq/i ! Convert lat and lon to more convenient unit cartesian vectors: do L=1,nob rlat=latob(L)*dtor; clat=cos(rlat); slat=sin(rlat) if(troflg) latobc(L)=latob(L) rlon=lonob(L)*dtor; clon=cos(rlon); slon=sin(rlon) xob(:,L)=(/clat*clon,clat*slon,slon/) rob(L)=(altob(L)-vmin)/uv-vrand! <- altitudes in hilbert vertical units enddo if(nrand<1 .or. nrand>273)stop 'nrand is invalid' if(nrand>5)then; call getqset7( qset7); if(nrand>7)call getqset13(qset13) else; call getqset5(nrand,qset5) endif if(nrand>91) call getqset5(3,qset3(:,:)) ! Project the data onto nrand differently-oriented Hilbert curves and sum ! the different estimated density estimates at the ob points in array wtob: rotnew=u0; do i=1,3; rotnew(i,i)=u1; enddo ! <- Identity matrix. wtob=0.0_dp do irand=1,nrand rotold=rotnew select case(nrand) case(1:5) qset=qset5(:,irand) case(6:7) qset=qset7(:,irand) case(8:13) qset=qset13(:,irand) case(14:91) i=1+mod(irand-1,13) j=1+(irand-1)/13 qset=mulqq(qset7(:,j),qset13(:,i)) case(92:273) i=1+mod(irand-1,13) j=1+(irand-1)/13 k=1+(j-1)/7 j=1+mod(j-1,7) qset=mulqq(mulqq(qset13(:,i),qset3(:,k)),qset7(:,j)) end select call qtorot(qset,rotnew) ! convert to an orthogonal matrix ! Form the relative rotation, rot (relative to previous irand iteration): rot=matmul(rotnew,transpose(rotold)) ! Get a new fraction of vh/4, based on ternary subdivisions such that the ! binary expansions (in units of vh) are non-terminating. This helps ensure ! that the resulting Hilbert curve projections are not accidentally similar. do L=1,nob xob(:,L)=matmul(rot,xob(:,L)) rob(L) =rob(L)+vrand ! Randomize by new vertical displacement enddo do L=1,nob sob(L)=0_dp call xs_to_hil48(ngen4,ngen,xob(:,L),rob(L), hilr) call hil8_to_r(ngen4+1,ngen,hilr(ngen4+1:ngen),sob(L)) call hil4_to_r(1,ngen4,hilr(1:ngen4),sob(L)) sob(L)=sob(L)+hilr(0) enddo! L ! Sort the data implicitly by assigning each a "rank": call sort(sob,rank); call invertperm(rank) select case(nor) case(0) call bsmoo0(nob,span,sob,rank,wtob) case(1) call bsmoo1(nob,span,sob,rank,wtob) case default stop 'In denest; this value of B-spline order, nor, is not supported' end select enddo! irand ! Convert the sum of Hilbert-parameter-relative densities to an average, ! and convert the units of density to number-per-span: wtob=wtob*span/nrand do L=1,nob dentripc=dentrip if(troflg) then if(abs(latobc(L)) =u1)then print& '("In denest; vmax-vmin too large for thin-shell assumption to be valid")' print'("Make this vertical range smaller, or vscale/hscale larger")' stop endif v=o2 do ngen4=0,ngen-1 if(v=nrand)exit enddo vrand=3.0_dp*vhq/i ! Convert lat and lon to more convenient unit cartesian vectors: do L=1,nob if(troflg) latobc(L)=latob(L) rlat=latob(L)*dtor; clat=cos(rlat); slat=sin(rlat) rlon=lonob(L)*dtor; clon=cos(rlon); slon=sin(rlon) xob(:,L)=(/clat*clon,clat*slon,slon/) rob(L)=(altob(L)-vmin)/uv-vrand! <- altitudes in hilbert vertical units enddo if(nrand<1 .or. nrand>104)stop 'nrand is invalid' if(nrand>5)then; call getqset8( qset8); if(nrand>8)call getqset13(qset13) else; call getqset5(nrand,qset5) endif ! Project the data onto nrand differently-oriented Hilbert curves and sum ! the different estimated density estimates at the ob points in array denob: rotnew=u0; do i=1,3; rotnew(i,i)=u1; enddo ! <- Identity matrix. denob=0_dp do irand=1,nrand rotold=rotnew select case(nrand) case(1:5) qset=qset5(:,irand) case(6:8) qset=qset8(:,irand) case(9:13) qset=qset13(:,irand) case(14:104) i=1+mod(irand-1,13) j=1+(irand-1)/13 qset=mulqq(qset8(:,j),qset13(:,i)) end select call qtorot(qset,rotnew) ! convert to an orthogonal matrix ! Form the relative rotation, rot (relative to previous irand iteration): rot=matmul(rotnew,transpose(rotold)) ! Get a new fraction of vh/4, based on ternary subdivisions such that the ! binary expansions (in units of vh) are non-terminating. This helps ensure ! that the resulting Hilbert curve projections are not accidentally similar. do L=1,nob xob(:,L)=matmul(rot,xob(:,L)) rob(L) =rob(L)+vrand ! Randomize by new vertical displacement enddo do L=1,nob sob(L)=0_dp call xs_to_hil48(ngen4,ngen,xob(:,L),rob(L), hilr) call hil8_to_r(ngen4+1,ngen,hilr(ngen4+1:ngen),sob(L)) call hil4_to_r(1,ngen4,hilr(1:ngen4),sob(L)) sob(L)=sob(L)+hilr(0) enddo! L ! Sort the data implicitly by assigning each a "rank": call sort(sob,rank); call invertperm(rank) select case(nor) case(0) call bsmoo0(nob,span,sob,rank,denob) case(1) call bsmoo1(nob,span,sob,rank,denob) case default stop 'In denest; this value of B-spline order, nor, is not supported' end select enddo! irand ! Convert the sum of Hilbert-parameter-relative densities to an average, ! and convert the units of density to number-per-span: denob=denob*span/nrand do L=1,nob dentripc=dentrip if(troflg) then if(abs(latobc(L))