diff --git a/.gitignore b/.gitignore index b568b28..fb31806 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,7 @@ *.o *.mod +Applications +example_* +.gitignore +temp diff --git a/example_Voronoi_clean/DSurfTomo.in b/example_Voronoi_clean/DSurfTomo.in index 7f28655..16a444a 100644 --- a/example_Voronoi_clean/DSurfTomo.in +++ b/example_Voronoi_clean/DSurfTomo.in @@ -18,5 +18,5 @@ surfdataTB.dat c: data file 0 c: kmaxLg 0 c: synthetic flag(0:real data,1:synthetic) 0.02 c: noiselevel -0.05 c: threshold -1 100 50 c: vorotomo,ncells,nrelizations +1.5 c: outlier factor (>=1.5) +1 200 50 c: vorotomo,ncells,nrelizations diff --git a/src/Makefile b/src/Makefile index f295801..00cc810 100644 --- a/src/Makefile +++ b/src/Makefile @@ -4,7 +4,7 @@ FFLAGS = -O -ffixed-line-length-none -ffloat-store\ -fbounds-check -m64 -mcmodel=medium F90SRCS = lsmrDataModule.f90 lsmrblasInterface.f90\ lsmrblas.f90 lsmrModule.f90 delsph.f90\ - aprod.f90 gaussian.f90 voronoiproj.f90 + aprod.f90 gaussian.f90 voronoiproj.f90 getpercentile.f90 FSRCS = surfdisp96.f slarnv.f slaruv.f OBJS = $(F90SRCS:%.f90=%.o) $(FSRCS:%.f=%.o) CalSurfG.o main.o all:$(CMD) diff --git a/src/getpercentile.f90 b/src/getpercentile.f90 new file mode 100644 index 0000000..c0c8555 --- /dev/null +++ b/src/getpercentile.f90 @@ -0,0 +1,51 @@ +SUBROUTINE getpercentile(N,array,q25,q75) + real q25 + real q75 + integer N + real,intent(in) :: array(*) + integer idx + + real RA(N) + RA(1:N) = array(1:N) + L=N/2+1 + IR=N + !The index L will be decremented from its initial value during the + !"hiring" (heap creation) phase. Once it reaches 1, the index IR + !will be decremented from its initial value down to 1 during the + !"retirement-and-promotion" (heap selection) phase. +10 continue + if(L > 1)then + L=L-1 + RRA=RA(L) + else + RRA=RA(IR) + RA(IR)=RA(1) + IR=IR-1 + if(IR.eq.1)then + RA(1)=RRA + idx = int(0.25*N) + q25 = RA(idx) + idx = int(0.75*N) + q75 = RA(idx) + return + end if + end if + I=L + J=L+L +20 if(J.le.IR)then + if(J < IR)then + if(RA(J) < RA(J+1)) J=J+1 + end if + if(RRA < RA(J))then + RA(I)=RA(J) + I=J; J=J+J + else + J=IR+1 + end if + goto 20 + end if + RA(I)=RRA + goto 10 + +END + diff --git a/src/main.f90 b/src/main.f90 index 51a0807..e9a4334 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -100,15 +100,15 @@ program SurfTomo integer ifsyn real averdws real maxnorm - real threshold0 + real threshold0,q25,q75 !For Poisson Voronoi inverison integer iproj,vorotomo,ncells,nrealizations,idx real hvratio ! OPEN FILES FIRST TO OUTPUT THE PROCESS - nout=36 - open(nout,file='lsmr.txt') + !nout=36 + !open(nout,file='lsmr.txt') ! OUTPUT PROGRAM INFOMATION write(*,*) @@ -366,10 +366,18 @@ program SurfTomo cbst(i) = obst(i) - dsyn(i) enddo + call getpercentile(dall,cbst,q25,q75) + datweight = 1.0 do i = 1,dall - datweight(i) = 0.01+1.0/(1+0.05*exp(cbst(i)**2*threshold0)) - cbst(i) = cbst(i)*datweight(i) + if (cbst(i)q75*threshold0) then + datweight(i) = 0.0 + cbst(i) = 0 + endif enddo + ! do i = 1,dall + ! datweight(i) = 0.01+1.0/(1+0.05*exp(cbst(i)**2*threshold0)) + ! cbst(i) = cbst(i)*datweight(i) + ! enddo do i = 1,nar rw(i) = rw(i)*datweight(iw(1+i)) @@ -531,14 +539,14 @@ program SurfTomo residual after weighting: ',mean*1000,'ms ',1000*std_devs,'ms ',& dnrm2(dall,cbst,1)/sqrt(real(dall)) - do i =1,dall - cbst(i)=cbst(i)/datweight(i) - enddo + !do i =1,dall + !cbst(i)=cbst(i)/datweight(i) + !enddo - mean = sum(cbst(1:dall))/dall - std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2) - write(*,'(a,f8.1,a,f8.2,a,f8.3)')' residual before weighting: ',mean*1000,'ms ',1000*std_devs,'ms ',& - dnrm2(dall,cbst,1)/sqrt(real(dall)) + !mean = sum(cbst(1:dall))/dall + !std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2) + !write(*,'(a,f8.1,a,f8.2,a,f8.3)')' residual before weighting: ',mean*1000,'ms ',1000*std_devs,'ms ',& + ! dnrm2(dall,cbst,1)/sqrt(real(dall)) write(66,'(i2,a)')iter,'th iteration...' ! write(66,'(a,f7.3)')'weight is:',weight write(66,'(a,f8.1,a,f8.2,a,f8.3)')'mean,std_devs and rms of & @@ -639,8 +647,8 @@ program SurfTomo endif - close(40) - close(nout) !close lsmr.txt + !close(40) + !close(nout) !close lsmr.txt close(66) !close surf_tomo.log deallocate(obst) diff --git a/src/voronoiproj.f90 b/src/voronoiproj.f90 index b57f347..8eb18c6 100644 --- a/src/voronoiproj.f90 +++ b/src/voronoiproj.f90 @@ -6,10 +6,10 @@ subroutine voronoiproj(leniw,lenrw,colg,nrow,rw,dres,goxd,dvxd,gozd,dvzd,depz,& integer leniw,lenrw integer nx,ny,nz ! integer iw(leniw) - integer colg(leniw),nrow(nd) + integer colg(lenrw),nrow(nd) real depz(nz) real rw(lenrw) - integer ncells + integer ncells,acells real dv(*),dres(*) real goxd,gozd,dvxd,dvzd real damp @@ -27,6 +27,7 @@ subroutine voronoiproj(leniw,lenrw,colg,nrow,rw,dres,goxd,dvxd,gozd,dvzd,depz,& integer maxnar,nzid integer iseed(4) real xs,ys,zs + real rx real atol,btol real conlim @@ -58,10 +59,10 @@ subroutine voronoiproj(leniw,lenrw,colg,nrow,rw,dres,goxd,dvxd,gozd,dvzd,depz,& !rad(ii) = cmb+depz(ii)*hvratio enddo - allocate(theta(ncells),phi(ncells),rrad(ncells),dws(ncells),norm(ncells)) + allocate(theta(ncells),phi(ncells),rrad(ncells),norm(ncells)) allocate(xpts(ncells),ypts(ncells),zpts(ncells),dis(ncells),xunknown(ncells)) allocate(rw_p(ndim)) - allocate(iw_p(2*ndim+1),row(ndim),col(ndim)) + allocate(iw_p(2*ndim+1),row(ndim),col(ndim),dws(ndim)) iseed(1:3) = (/38,62,346/) iseed(4) = 2*iproj+1 @@ -72,6 +73,23 @@ subroutine voronoiproj(leniw,lenrw,colg,nrow,rw,dres,goxd,dvxd,gozd,dvzd,depz,& call slarnv(1,iseed,ncells,rrad) rrad = radius-rrad*depz(nz-1)*hvratio + ! adaptive cells based on dws, assume 1/2 of all ncells are used + ! as adaptive cells + dws = 0 + do ii = 1,lenrw + dws(colg(ii)) = dws(colg(ii))+abs(rw(ii)) + enddo + acells = int(ncells/2.0) + do ii = ncells-acells,ncells + call random_index(idx,dws) + ix = mod(idx,nx-2) + iy = idx/(nx-2) + iz = idx/((nx-2)*(ny-2)) + theta(ii) = (gozd+(ix+1)*dvzd)*pi/180 + phi(ii) = pi/2-(goxd-(iy+1)*dvxd)*pi/180 + rrad(ii) = radius-depz(iz+1)*hvratio + enddo + xpts = rrad*sin(phi)*cos(theta) ypts = rrad*sin(phi)*sin(theta) zpts = rrad*cos(phi) @@ -127,10 +145,6 @@ subroutine voronoiproj(leniw,lenrw,colg,nrow,rw,dres,goxd,dvxd,gozd,dvzd,depz,& iwgp(1) = lenrwgp iwgp(nzid+2:nzid*2+1) = colgp(1:nzid) - dws = 0 - !do ii=1,nzid - !dws(iwgp(1+ii+nzid)) = dws(iwgp(1+ii+nzid))+abs(rwgp(ii)) - !enddo norm = 0 do ii=1,nzid norm(iwgp(1+ii+nzid)) = norm(iwgp(1+ii+nzid))+rwgp(ii)**2 @@ -169,6 +183,14 @@ subroutine voronoiproj(leniw,lenrw,colg,nrow,rw,dres,goxd,dvxd,gozd,dvzd,depz,& do ii = 1,ncells xunknown(ii) = xunknown(ii)/norm(ii) enddo + norm = (norm**2-0.01)*nd + do ii = 1,ncells + if (norm(ii)<0.01) then + call random_number(rx) + xunknown(ii) = xunknown(ii)+rx-0.5 + endif + enddo + dv(1:ndim) = 0 call aprod(2,ncells,ndim,dv,xunknown,leniw_p,lenrw_p,iw_p,rw_p) deallocate(grow,gcol,subrow) @@ -178,4 +200,23 @@ subroutine voronoiproj(leniw,lenrw,colg,nrow,rw,dres,goxd,dvxd,gozd,dvzd,depz,& deallocate(lat,lon,rad) deallocate(iwgp,colgp,rwgp) - end subroutine +contains +subroutine random_index( idx, weights ) + integer :: idx + real, intent(in) :: weights(:) + + real x, wsum, prob + + wsum = sum( weights ) + + call random_number( x ) + + prob = 0 + do idx = 1, size( weights ) + prob = prob + weights( idx ) / wsum !! 0 < prob < 1 + if ( x <= prob ) exit + enddo +end subroutine random_index + + +end subroutine