var/da/da_test/da_check_gradient.inc

   1 subroutine da_check_gradient(grid, config_flags, cv_size, xhat, cv, pdx, itertest, xbx, be, iv, y, re, j_cost)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Gradient test
   5    ! Adopted from grtest.f90 of GSI by Xin Zhang , September, 2011
   6    !-----------------------------------------------------------------------
   7
   8    implicit none
   9
  10    type(domain),               intent(inout) :: grid
  11    type(grid_config_rec_type), intent(inout) :: config_flags
  12
  13    real,    intent(in)    ::   pdx
  14    integer, intent(in   ) :: itertest
  15    integer, intent(in)           :: cv_size ! Size of cv array.
  16    real,    intent(inout)        :: xhat(1:cv_size)  ! control variable (local).
  17    real,    intent(inout)           :: cv(1:cv_size)    ! control variable (local).
  18    type (xbx_type),   intent(inout) :: xbx   ! Header & non-gridded vars.
  19    type (be_type),    intent(in) :: be    ! background error structure.
  20    type (iv_type),    intent(inout) :: iv    ! ob. increment vector.
  21    type (y_type),     intent(inout) :: y             ! y = h (xa)
  22    type (y_type), intent(inout)  :: re    ! residual (o-a) structure.
  23    type (j_type), intent(inout)  :: j_cost                 ! cost function
  24
  25    real    :: xdir(1:cv_size) , grad(1:cv_size), yhat(1:cv_size)
  26    real    :: zfy,zf0,zdf0,za,zfa,zdfa
  27    real    :: zabuf(itertest),zfabuf(itertest),ztf2buf(itertest)
  28    real    :: ZT1,ZB,ZFB,ztco,ZTC1,ZT2,ZTC1A,ZTC2,ZTF2
  29    real    :: ZAL,ZFAL,ZBL,ZFBL,ZTF2L
  30    real    :: ZTC00,ZTC02,ZTC10,ZTC12
  31    real    :: ZERMIN,ZT1TST,ZREF
  32    integer                           :: jp_start, jp_end       ! Start/end indices of Jp.
  33    integer :: ibest,idig
  34    integer :: ii, jj
  35
  36    call da_trace_entry("da_check_gradient")
  37
  38    jp_start   = be % cv % size_jb + be % cv % size_je + 1
  39    jp_end     = be % cv % size_jb + be % cv % size_je + be % cv % size_jp
  40
  41    call da_message((/' gradient test starting'/))
  42
  43    if (pdx<=EPSILON(pdx)) then
  44       if (rootproc) write(6,*)'grtest, pdx=',pdx
  45       write(6,*)'grtest: pdx too small',pdx
  46       call da_trace_exit("da_check_gradient")
  47       return
  48    endif
  49
  50    !----------------------------------------------------------------------------
  51    ! [1] Initial point
  52    !----------------------------------------------------------------------------
  53
  54    call da_initialize_cv (cv_size, cv)
  55    call da_initialize_cv (cv_size, xhat)
  56
  57    ! Initialize cv values with random data:
  58    call random_number(xhat(:))
  59    xhat(:) = xhat(:) - 0.5
  60    if (rootproc) write(6,*)'grtest: use random_number(xhat)'
  61
  62    yhat=xhat
  63
  64    call da_calculate_j(1, 1, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, &
  65                        be%cv%size_jl, be%cv%size_jt, xbx, be, iv, yhat, cv, re, y, j_cost, grid, config_flags)
  66
  67    call da_calculate_gradj(1,1,cv_size,be%cv%size_jb,be%cv%size_je,be%cv%size_jp, &
  68                            be%cv%size_jl, be%cv%size_jt, xbx, be, iv, yhat+cv, y, grad, grid, config_flags, re)
  69
  70    zfy = j_cost%total
  71
  72    !----------------------------------------------------------------------------
  73    ! [1.1] Define perturbation direction ZH
  74    !----------------------------------------------------------------------------
  75
  76    call da_message((/' The test direction is the opposite of the gradient '/))
  77
  78    xdir  = -1.0  * grad
  79    !----------------------------------------------------------------------------
  80    ! [1.2] Set function f value and derivative at origin
  81    !----------------------------------------------------------------------------
  82
  83    zf0=zfy
  84    zdf0=da_dot_cv(cv_size,grad,xdir,grid,be%cv_mz,be%ncv_mz &
  85 #if (WRF_CHEM == 1)
  86                ,be%cv_mz_chemic &
  87 #endif
  88                ,jp_start,jp_end)
  89    if (rootproc) write(6,*)'grtest: F(0)=',zf0,' DF(0)=',zdf0
  90
  91    IF (ZDF0>0.0) write(6,*) 'GRTEST Warning, DF should be negative'
  92    IF (ABS(ZDF0) < SQRT(EPSILON(ZDF0))) THEN
  93       if (rootproc) write(6,*) 'GRTEST WARNING, DERIVATIVE IS TOO SMALL'
  94    ENDIF
  95
  96    !----------------------------------------------------------------------------
  97    ! [2] Loop on test point
  98    !----------------------------------------------------------------------------
  99
 100    ztf2buf(1)=0.0
 101
 102    DO jj=1,itertest
 103
 104       za=pdx*(10.0**(jj-1))
 105
 106       if (rootproc) write(6,*)'grtest iter=',jj,' alpha=',za
 107
 108       !----------------------------------------------------------------------------
 109       ! [2.1] Compute f and df at new point y=x+a.h
 110       !----------------------------------------------------------------------------
 111
 112       do ii=1,cv_size
 113          yhat(ii) = xhat(ii) + za * xdir(ii)
 114       end do
 115
 116       call da_calculate_j(1, 1, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, &
 117                           be%cv%size_jl, be%cv%size_jt, xbx, be, iv, yhat, cv, re, y, j_cost, grid, config_flags)
 118
 119       call da_calculate_gradj(1,1,cv_size,be%cv%size_jb,be%cv%size_je,be%cv%size_jp, &
 120                               be%cv%size_jl, be%cv%size_jt, xbx, be, iv, yhat+cv, y, grad, grid, config_flags, re)
 121
 122       zfy = j_cost%total
 123
 124       zfa=zfy
 125       zdfa=da_dot_cv(cv_size,grad,xdir,grid,be%cv_mz,be%ncv_mz &
 126 #if (WRF_CHEM == 1)
 127                ,be%cv_mz_chemic &
 128 #endif
 129                ,jp_start,jp_end)
 130
 131       if (rootproc) write(6,*)'grtest: alpha=',za,' F(a)=',zfa,' DF(a)=',zdfa
 132
 133       zabuf(jj)=za
 134       zfabuf(jj)=zfa
 135
 136       !----------------------------------------------------------------------------
 137       ! [2.2] Quantity TC0=f(a)/f(0)-1
 138       !----------------------------------------------------------------------------
 139
 140 !         if f is continuous then TC0->1 at origin,
 141 !         at least linearly with a.
 142
 143       IF (ABS(zf0)<=TINY(zf0)) THEN
 144 !           do not compute T1 in this unlikely case
 145          if (rootproc) write(6,*) 'grtest: Warning: zf0 is suspiciously small.'
 146          if (rootproc) write(6,*) 'grtest: F(a)-F(0)=',zfa-zf0
 147       ELSE
 148          ztco=zfa/zf0-1.0
 149          if (rootproc) write(6,*)'grtest: continuity TC0=',ztco
 150       ENDIF
 151
 152       !----------------------------------------------------------------------------
 153       !                     f(a)-f(0)
 154       ! [2.3] Quantity T1=-----------
 155       !                      a.df(0)
 156       !----------------------------------------------------------------------------
 157
 158 !         if df is the gradient then T1->1 at origin,
 159 !         linearly with a. T1 is undefined if df(0)=0.
 160
 161       IF (ABS(za*zdf0)<=SQRT(TINY(zf0))) THEN
 162          if (rootproc) write(6,*)'grtest: Warning: could not compute ',&
 163           & 'gradient test T1, a.df(0)=',za*zdf0
 164       ELSE
 165          zt1=(zfa-zf0)/(za*zdf0)
 166          if (rootproc) write(6,*)'grtest: gradient T1=',zt1
 167       ENDIF
 168
 169       !----------------------------------------------------------------------------
 170       ! [2.4] Quantity TC1=( f(a)-f(0)-a.df(0) )/a
 171       !----------------------------------------------------------------------------
 172
 173 !         if df is the gradient and df is continuous,
 174 !         then TC1->0 linearly with a.
 175       ZTC1=(ZFA-ZF0-ZA*ZDF0)/ZA
 176       if (rootproc) write(6,*)'grtest: grad continuity TC1=',ZTC1
 177
 178       !----------------------------------------------------------------------------
 179       ! [2.5] Quantity T2=( f(a)-f(0)-a.df(0) )*2/a**2
 180       !----------------------------------------------------------------------------
 181
 182 !         if d2f exists then T2 -> d2f(0) linearly with a.
 183       ZT2=(ZFA-ZF0-ZA*ZDF0)*2.0/(ZA**2)
 184       if (rootproc) write(6,*)'grtest: second derivative T2=',ZT2
 185
 186       !----------------------------------------------------------------------------
 187       ! [2.6] Quantity TC1A=df(a)-df(0)
 188       !----------------------------------------------------------------------------
 189
 190 !         if df is the gradient in a and df is continuous,
 191 !         then TC1A->0 linearly with a.
 192       ZTC1A=ZDFA-ZDF0
 193       if (rootproc) write(6,*)'grtest: a-grad continuity TC1A=',ZTC1A
 194
 195       !----------------------------------------------------------------------------
 196       ! [2.7] Quantity TC2=( 2(f(0)-f(a))+ a(df(0)+df(a))/a**2
 197       !----------------------------------------------------------------------------
 198
 199 !         if f is exactly quadratic, then TC2=0, always: numerically
 200 !         it has to -> 0 when a is BIG. Otherwise TC2->0 linearly for
 201 !         small a is trivially implied by TC1A and T2.
 202       ZTC2=(2.0*(ZF0-ZFA)+ZA*(ZDF0+ZDFA))/(ZA**2)
 203       if (rootproc) write(6,*)'grtest: quadraticity TC2=',ZTC2
 204
 205       !----------------------------------------------------------------------------
 206       !                     2   f(0)-f(b)   f(a)-f(b)
 207       ! [2.8] Quantity TF2=---( --------- + --------- )
 208       !                     a       b          a-b
 209       !----------------------------------------------------------------------------
 210
 211 !         if 0, a and b are distinct and f is quadratic then
 212 !         TF2=d2f, always. The estimate is most stable when a,b are big.
 213 !         This test works only after two loops, but it is immune against
 214 !         gradient bugs.
 215
 216       IF (jj>=2) THEN
 217          ZB =ZABUF (jj-1)
 218          ZFB=ZFABUF(jj-1)
 219          ZTF2=2.0/ZA*((ZF0-ZFB)/ZB+(ZFA-ZFB)/(ZA-ZB))
 220          if (rootproc) write(6,*)'grtest: convexity ZTF2=',ZTF2
 221          ztf2buf(jj)=ztf2
 222       ENDIF
 223
 224 ! End loop
 225    ENDDO
 226
 227    !----------------------------------------------------------------------------
 228    ! [3] Comment on the results
 229    !----------------------------------------------------------------------------
 230
 231 !       TC0(0)/TC0(2)<.011 -> df looks continuous
 232 !       item with (T1<1 and 1-T1 is min) = best grad test item
 233 !       reldif(TF2(last),TF2(last-1)) = precision on quadraticity
 234
 235    !----------------------------------------------------------------------------
 236    !       3.1 Fundamental checks
 237    !----------------------------------------------------------------------------
 238
 239    if (rootproc) then
 240       write(6,*) 'GRTEST: TENTATIVE CONCLUSIONS :'
 241
 242       ZTC00=ABS(zfabuf(1)-zf0)
 243       ZTC02=ABS(zfabuf(3)-zf0)
 244       IF( ZTC00/zabuf(1)  <=  1.5*(ZTC02/zabuf(3)) )THEN
 245          write(6,*) 'GRTEST: function f looks continous.'
 246       ELSE
 247          write(6,*) 'GRTEST: WARNING f does not look continuous',&
 248           & ' (perhaps truncation problem)'
 249       ENDIF
 250
 251    !----------------------------------------------------------------------------
 252    !       3.2 Gradient quality
 253    !----------------------------------------------------------------------------
 254
 255       IF (ABS(zdf0)<=SQRT(TINY(zf0))) THEN
 256          write(6,*) 'GRTEST: The gradient is 0, which is unusual !'
 257          ZTC10=ABS(zfabuf(1)-zf0)
 258          ZTC12=ABS(zfabuf(3)-zf0)
 259          IF( ZTC10/zabuf(1)**2  <=  1.1*ZTC12/zabuf(3)**2)THEN
 260             write(6,*)'GRTEST: The gradient looks good anyway.'
 261          ENDIF
 262       ELSE
 263 !        Find best gradient test index
 264          ZERMIN=HUGE(0.0)
 265          ibest=-1
 266          DO jj=1,itertest
 267             ZT1TST=(zfabuf(jj)-zf0)/(zabuf(jj)*zdf0)
 268             ZT1TST=ABS(ZT1TST-1.0)
 269             IF (ZT1TST<ZERMIN) THEN
 270                ibest=jj
 271                ZERMIN=ZT1TST
 272             ENDIF
 273          ENDDO
 274          IF(ibest == -1)THEN
 275             write(6,*)'GRTEST: gradient test problem : bad ',&
 276              & 'gradient, non-convex cost, or truncation errors ?'
 277          ELSE
 278             idig=INT(-LOG(ZERMIN+TINY(ZERMIN))/LOG(10.0))
 279             write(6,*)'GRTEST: the best gradient test found has ',&
 280              & idig,' satisfactory digits.'
 281             IF(idig <= 1)THEN
 282                write(6,*)'GRTEST: SAYS: THE GRADIENT IS VERY BAD.'
 283             ELSEIF(idig <= 3)THEN
 284                write(6,*)'GRTEST: SAYS: THE GRADIENT IS SUSPICIOUS.'
 285             ELSEIF(idig <= 5)THEN
 286                write(6,*)'GRTEST: SAYS: THE GRADIENT IS ACCEPTABLE.'
 287             ELSE
 288                write(6,*)'GRTEST: SAYS: THE GRADIENT IS EXCELLENT.'
 289             ENDIF
 290
 291             IF (ibest<=itertest-2) THEN
 292                ZTC10=ABS(zfabuf(ibest         )-zf0-zabuf(ibest         )*zdf0)/zabuf(ibest         )
 293                ZTC12=ABS(zfabuf(ibest+2)-zf0-zabuf(ibest+2)*zdf0)/zabuf(ibest+2)
 294                IF(ZTC10/zabuf(ibest) <=  1.1*ZTC12/zabuf(ibest+2) )THEN
 295                   write(6,*)'GRTEST: Gradient convergence looks good.'
 296                ELSE
 297                   write(6,*)'GRTEST: Gradient convergence is suspicious.'
 298                ENDIF
 299             ELSE
 300                write(6,*)'GRTEST: could not check grad convergence.'
 301             ENDIF
 302          ENDIF
 303       ENDIF
 304
 305    !----------------------------------------------------------------------------
 306    !         3.3 Quadraticity
 307    !      finite difference quadraticity test (gradient-free)
 308    !----------------------------------------------------------------------------
 309
 310       ZTF2=ztf2buf(itertest)
 311       ZTF2L=ztf2buf(itertest-1)
 312       write(6,*) 'GRTEST: finite diff. d2f estimate no1:',ZTF2
 313       write(6,*) 'GRTEST: finite diff. d2f estimate no2:',ZTF2L
 314       ZREF=(ABS(ZTF2L)+ABS(ZTF2))/2.0
 315       IF (ZREF<=TINY(ZREF)) THEN
 316          write(6,*) 'GRTEST: they are too small to decide whether ',&
 317           & 'they agree or not.'
 318       ELSE
 319          idig=INT(-LOG(ABS(ZTF2L-ZTF2)/ZREF+TINY(ZTF2))/LOG(10.0))
 320          write(6,*) 'GRTEST: the fin.dif. estimates of d2f ',&
 321           & 'have ',idig,' satisfactory digits.'
 322          IF(idig <= 1)THEN
 323             write(6,*) 'GRTEST: THE FD-QUADRATICITY IS BAD.'
 324          ELSEIF(idig <= 3)THEN
 325             write(6,*) 'GRTEST:: THE FD-QUADRATICITY IS SUSPICIOUS.'
 326          ELSEIF(idig <= 5)THEN
 327             write(6,*) 'GRTEST: THE FD-QUADRATICITY IS ACCEPTABLE.'
 328          ELSE
 329             write(6,*) 'GRTEST: THE FD-QUADRATICITY IS EXCELLENT.'
 330          ENDIF
 331       ENDIF
 332
 333       write(6,*) 'grtest: Goodbye.'
 334    endif
 335
 336    call da_trace_exit("da_check_gradient")
 337
 338 end subroutine da_check_gradient
 339
 340