#!/usr/bin/python # Simple 2D univariate analysis - F.Bouttier Dec 2015 import random,math,sys,os import numpy as np import subprocess from scipy.optimize import minimize import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt ### technical params verb=1 # 1 for verbose output epsilon=1.e-10 # to avoid div-by-zero osname=subprocess.Popen('uname',stdout=subprocess.PIPE).stdout.read() ######## user parameters ############ ### define 1D model state ni=1 ; nj=51 # model grid size imin=0. ; imax=0. ; jmin=0. ; jmax=10. # geographical model boundaries ### define 2D model state #ni=23 ; nj=23 # model grid size #imin=0. ; imax=10. ; jmin=0. ; jmax=10. # geographical model boundaries ### define obs network #obsinit='onemiddle' # obs initialization method obsinit='twomiddle' # obs initialization method errobs=1. # R obs standard error ### define background bgtype='zero' # background type: zeroes everywhere ### background error model sigmab=1. # B uniform background standard error bdist=0.5 # characteristic distance for bg error correlations bcortype='gaussian' # gaussian error correlation #bcortype='triangle' # triangular error correlation ### define analysis method #anamethod='cressman' ; minweight=0.1 # Cressman method anamethod='3dvar' ; niter=10 # variational ### plotting params PLOTFILE='anaplot.png' PLOT=1 # 1 to generate plot VIEW=1 # 1 to view plotfile on screen VIEWER='gqview 2>/dev/null' # linux shell command to display png file on screen if osname=='Darwin\n' : VIEWER='open' # Mac OSX compatibility print '\n--- welcome to ana2d objective analysis ---\n' # precompute model resolution di=0. ; dj=0. if ni>1 : di=(imax-imin)/float(ni-1) if nj>1 : dj=(jmax-jmin)/float(nj-1) def idx2real( ij ): # convert gridpoint indices to physical coordinates (i,j)=ij zi=imin+i*di zj=jmin+j*dj return(zi,zj) ######################## setup obs #################### if obsinit=='onemiddle': # obs = 1 in middle of domain nobs=1 y =[ 1. ] yi=[ (imin+imax)/2. ] yj=[ (jmin+jmax)/2. ] if obsinit=='twomiddle': # obs = 1 in middle of domain nobs=2 y =[ 0.6 , 1. ] yi=[ 0.6*imin + 0.4*imax , 0.4*imin + 0.6*imax ] yj=[ 0.6*jmin + 0.4*jmax , 0.4*jmin + 0.6*jmax ] ### print obs ### if verb : print '---observations:---' print 'obsinit=',obsinit,' Nobs=',nobs for i in range(nobs): print '%10.3f at ( %10.3f , %10.3f )' % ( y[i],yi[i],yj[i] ) print '-------------------\n' ######################## setup background ##################### ### setup background state ##### print '---background state:--- (ni,nj)=',ni,nj if bgtype=='zero': print 'xb=zero everywhere' xb = np.zeros( (ni,nj) ) print '-----------------------\n' ### setup background error covs ### def bcov(ia,ja,ib,jb): # compute bg error cov between points (ia,ja) and (ib,jb) d=math.sqrt( (ia-ib)**2 + (ja-jb)**2 ) if bcortype=='triangle': c = max( 0. , 1. - d / (2.*bdist) ) # linear triangular correl elif bcortype=='gaussian': c = math.exp( - d**2 / ( 2. * bdist**2 ) ) # gaussian correl else: print 'error, unknown bcortype=',bcortype sys.exit(1) c = c * sigmab * sigmab # convert correlation to covariance return(c) ### print background covariance ### strucfun=None if verb: print '--- B cov in middle of domain:--- (j cross-section)' # define reference point zia=(imax-imin)/2. zja=(jmax-jmin)/2. # travel the cross-section print ' idx x cov(x)' strucfun=[None]*nj for jb in range(nj): (zib,zjb)=idx2real( (ni/2,jb) ) strucfun[jb]=bcov( zia,zja,zib,zjb ) print '%6.6d %10.3f %10.5f' % ( jb,zjb,strucfun[jb] ) print '--------------------------------\n' ###################### define H the observation operator ########## def nearestgp(i,j): # compute nearest gridpoint from real coords (i,j) kio=int( float(ni)*(i-imin)/(imax-imin+epsilon) + 0.5 ) kjo=int( float(nj)*(j-jmin)/(jmax-jmin+epsilon) + 0.5 ) assert( (kio>=0) and (kio<=ni) ),'out-of-domain obs i-coord kio='+repr(kio) assert( (kjo>=0) and (kjo<=nj) ),'out-of-domain obs j-coord kjo='+repr(kjo) return(kio,kjo) def hobsop(x,i,j): # interpolate state x at real coords (i,j) if verb: print ' obsop: i,imin,imax=',i,imin,imax if verb: print ' obsop: j,jmin,jmax=',j,jmin,jmax (kio,kjo)=nearestgp(i,j) if verb: print ' obsop: nearest gridpoint (i,j)=',kio,kjo hx=x[kio,kjo] if verb: print ' obsop: kio,kjo,h(x)=',kio,kjo,hx return( hx ) # interpolate using nearest-neighbour method ####################### run analysis ##################### ### compute innovations ### if verb: print '---compute innovations---' dep=[None]*nobs for k in range(nobs): hxb=hobsop( xb , yi[k] , yj[k] ) dep[k] = y[k] - hxb if verb : print 'idx,y,h(xb),y-h(xb)= %5.5d %12.3f %12.3f %12.3f' % (k,y[k],hxb,dep[k]) if verb: print '-------------------------\n' ########################################### ### analysis solver: Cressman increment ### ########################################### if anamethod=='cressman': inc = np.zeros( (ni,nj) ) for i in range(ni): for j in range(nj): zi=imin+i*di # real coords of gridpoint(i,j) zj=jmin+j*dj ## do the following for all model gridpoints ## zinc=0. ; zsum=0. for k in range(nobs): if verb: print ' obs idx=',k zw=bcov( zi,zj , yi[k],yj[k] ) # weight between model & obs points zinc = zinc + zw * dep[k] zsum = zsum + zw if zsum>minweight : inc[i,j] = zinc / zsum if verb: print 'Cressman ana at i,j=',i,j,' zi,zj=',zi,zj,' zsum,inc=',zsum,inc[i,j] xa = xb + inc ########################################### ### analysis solver: variational ########## ########################################### if anamethod=='3dvar': ### define mapping between 2D (i,j) fields and control variables ### ncv=ni*nj print '3DVar: size of control variable=',ncv def idx_ij2cv(i,j): # convert indices: 2D -> CV return( i*nj + j ) def idx_cv2ij(k): # convert indices: CV -> 2D i = int( float(k)/float(nj) ) j = k - i*nj return(i,j) def vec_cv2ij(xcv): # convert vector: CV -> 2D x2d = np.zeros( (ni,nj) ) for k in range(ncv): (i,j) = idx_cv2ij(k) x2d[i,j] = xcv[k] return(x2d) def vec_ij2cv(x2d): # convert vector: 2D -> CV xcv=np.zeros(ncv) for k in range(ncv): (i,j) = idx_cv2ij(k) xcv[k] = x2d[i,j] return(xcv) #### build explicit H matrix ### print '3DVar: building H matrix' hmatrix = np.zeros( (nobs,ncv) ) for k in range(nobs): (io,jo)=nearestgp( yi[k],yj[k] ) hmatrix[ k , idx_ij2cv(io,jo) ] = 1. # interpolation coeff for obs k if verb : print hmatrix ### build explicit B matrix ### print '3DVar: building B matrix' bmatrix = np.zeros( (ncv,ncv) ) for k1 in range(ncv): for k2 in range(ncv): (i1,j1)=idx2real( idx_cv2ij(k1) ) (i2,j2)=idx2real( idx_cv2ij(k2) ) bmatrix[k1,k2] = bcov( i1,j1,i2,j2 ) if verb : print bmatrix # define inverse of R (diagonal) ri = 1./(errobs*errobs) # departure computation from state vector (not preconditioned) def calcdep(xcv): x2d=vec_cv2ij(xcv) dep=np.zeros( (nobs) ) for k in range(nobs): hx=hobsop( x2d , yi[k] , yj[k] ) dep[k] = y[k] - hx return(dep) ### define cost-function and its gradient ### precon='B' # preconditioner B: background is B, slow convergence #precon='sqrtB' # preconditioner L, LLT=B: background is BB, fast convergence def jcost(u): # u is a B-preconditioned vector in control variable space if precon=='B': bu = np.dot( bmatrix , u ) jb = np.dot( np.transpose(u) , bu ) # Jb = u^T B u gradjb = 2. * bu # Jb' = 2 B u else: jb = np.dot( np.transpose(u) , u ) gradjb = 2. * u x = vec_ij2cv(xb) + np.dot( bmatrix , u ) # x = xb + B u dep = calcdep(x) # d = y - Hx jo = np.dot( np.transpose(dep) , ri*dep) # Jo = (y-Hx)^T R-1 (y-Hx) gradjo = -2. * np.dot( np.transpose( hmatrix ) , ri * dep ) #dJo/dx gradjo = np.dot( bmatrix , gradjo ) # dJo/du = dJo/dx.dx/du j = jb + jo gradj = gradjb + gradjo return( j , gradj ) def jvalue(u) : return( jcost(u)[0] ) # return only J(u) def jgrad(u) : return( jcost(u)[1] ) # return only gradJ(u) def printj(u) : # print J and ||J'|| (ju,gradju) = jcost( u ) normgradju=np.linalg.norm(gradju) print ' j=',ju,' normgradJ=',normgradju return(ju,gradju,normgradju) # test J on xb ub=np.zeros(ncv) # u = B(x-xb) print 'xb ->', (jub,gradjub,normgradjub) = printj(ub) # test correctness of the gradient (Taylor formula) print 'gradient test: a deltaJ gradJ.h ratio' verb=0 for k in range(-1,-8,-1): a = 10.**k h = a*gradjub # perturbation h (ju,gradju) = jcost( ub+h ) # J( ub + h ) varlin= np.dot( gradjub , h ) # linear variation of J predicted by gradjub print ' %10.3g %12.3g %12.3g %15.8f' % (a, ju-jub, varlin , varlin/( ju-jub+epsilon) ) # minimization inc=np.zeros(ncv) niter=30 print '#### start BFGS minim ####' result=minimize(jvalue, ub, method='BFGS', jac=jgrad, callback=printj, \ options={'disp':True,'maxiter':niter,'gtol':normgradjub/100.} ) print 'success:',result.success print 'message:',result.message print '#### end BFGS minim ####' ua=result.x # this is the analysis print 'xa ->', (jua,gradjua,normgradjua) = printj(ua) print " ratio of gradient norms: J'a/J'b=",normgradjua/(epsilon+normgradjub) inc = np.dot( bmatrix , result.x ) xa = xb + vec_cv2ij(inc) ### plot the analysis ### if PLOT : fig,ax = plt.subplots(1) ax.grid(True) if ni==1 : # matplotlib x=f(j) 1D plot ax.set_ylim([-0.1,1.1] ) ax.set_title('1D analysis') ax.set_xlabel('space coord') ax.set_ylabel('field value') jcoords=np.linspace(jmin,jmax,nj) # list of i-coordinates in x # plot xa(j) as curve plt.plot( jcoords,xa[0,:] ,'o-',linewidth=5,color='cyan') # plot observations plt.plot( yj,y ,marker=r'$\asterisk$',linestyle='', \ color='green',markeredgecolor=None,markersize=40) # plot B structure function if available if strucfun: plt.plot( jcoords,strucfun[:] ,'-',linewidth=3,color='blue') else : # matplotlib x=f(i,j) 2D plot ax.set_xlabel('i') ax.set_xlabel('j') ax.set_xlim([0,ni-1]) ax.set_ylim([0,nj-1]) plt.colorbar( ax.pcolor(xa,cmap=plt.cm.Blues) ) # plot backend: make PNG file fig.set_size_inches(20, 10) fig.savefig(PLOTFILE,bbox_inches='tight') os.system('ls -l '+PLOTFILE) if VIEW: os.system(VIEWER+' '+PLOTFILE) ### terminate script ### print 'finished ',sys.argv[0],' ok.\n' sys.exit(0)