Commit 74a15fa7 authored by Stephan Thober's avatar Stephan Thober
Browse files

parsed parameters from NelMin optimization back to program-level

parent c0870758
......@@ -40,6 +40,7 @@ program ED_Kriging
integer(i4) :: year, month, day ! current date
integer(i4) :: netcdfid ! id of netcdf files
real(dp), dimension(:,:), allocatable :: tmp_array ! temporal array for output
real(dp) :: param(3) ! variogram parameters
call Timer
call ReadDataMain
......@@ -56,7 +57,8 @@ program ED_Kriging
call ReadDataMeteo
print*, 'finished reading of meteorological data'
! estimate variogram
call setVario
call setVario(param)
stop 'TESTING OPTIMIZATION'
! write variogram
if (flagVario) call WriteDataMeteo(0,0,2)
!
......
......@@ -50,12 +50,14 @@ subroutine ReadDataMain
ios = 0
allocate (beta(nParam))
filename = trim(fNameVARIO)
open (unit=20, file=filename, STATUS='OLD', ACTION='READ')
if ( .not. flagVario ) then
filename = trim(fNameVARIO)
open (unit=20, file=filename, STATUS='OLD', ACTION='READ')
read (20,1) dummy
read (20, *,iostat=ios) (beta(i), i=1,nParam), i
close (20)
close (20)
end if
! if vario is read from file take read in vario type
if ( .not. flagVario ) then
......@@ -288,6 +290,7 @@ subroutine ReadDataMeteo
! read yearly data file
write (dummy, 2) MetSta(i)%Id
fileName = trim(dataPathIn)//trim(dummy)
print *, 'read file: '//trim(fileName)
open (60, file=fileName, status='old', action='read', iostat=ios)
read (60, *) dummy
!
......
!*****************************************************************************
! Main_opti.for
! AUTHOR
! Luis E. Samaniego-Eguiguren, IREUS, 28.05.1999
! DESCRIPTION
! Initialization of the Nonlinear Optimization Subroutine GRG2
! The function to be optimized is suplied in subroutine GCOMP
! DEFINITION OF INPUT/OUTPUT-CHANNELS
! SCREEN: *
!
! SCREEN OUTPUT: *
! OUTPUT FILE: 6
!******************************************************************************
subroutine OPTI
use VarFit
use mo_kind, only: i4, dp
use mo_obj_func, only: obj_func
use mo_nelmin, only: nelmin
! parameters for Nelder-Mead algorithm
real(dp) :: pstart(3) ! Starting point for the iteration.
real(dp) :: pmin(3) ! optimized parameter set
real(dp) :: prange(3, 2) ! Range of parameters (upper and lower bound).
real(dp) :: varmin ! the terminating limit for the variance of the function values. varmin>0 is required
real(dp) :: step(3) ! determines the size and shape of the initial simplex. The relative magnitudes of its elements should reflect the units of the variables. size(step)=size(start)
integer(i4) :: konvge ! the convergence check is carried out every konvge iterations
integer(i4) :: maxeval ! the maximum number of function evaluations. default: 1000
real(dp) :: funcmin
integer(i4) :: neval ! the number of function evaluations used.
integer(i4) :: numrestart ! the number of restarts.
integer(i4) :: ierror ! error indicator.
! 0: no errors detected.
! 1: varmin or konvge have an illegal value.
! 2: iteration terminated because maxeval was exceeded without convergence.
real(dp), allocatable :: history(:)
! ! inputs for GRG2
! IMPLICIT DOUBLE PRECISION(A-H,O-Z), INTEGER(I,J,L,M,N)
! INTEGER*4 NCORE,NNVARS,NFUN,MAXBAS,MAXHES,LASTZ
! INTEGER*4 M,N,MP1,NPMP1,NBMAX,NNBMAX,NPNBMX,MPNBMX,NRTOT
! LOGICAL MAXIM,INPRNT,OTPRNT
! DIMENSION BLVAR(100),BUVAR(100),BLCON(10),BUCON(10)
! DIMENSION RAMCON(10),RAMVAR(10),INBIND(10),Z(5000)
! DIMENSION NONBAS(10),REDGR(10),DEFAUL(20),TITLE(19)
! DIMENSION XX(100),FCNS(10),RMULTS(10)
! DATA BIG/1.D31/
! Initialization of Nelder-Mead
pstart = (/0.05, 100., 2./) ! Starting point for the iteration.
prange(:, 1) = (/0., 0., 0./) ! Range of parameters (lower bound).
prange(:, 2) = (/2., 300., 5./) ! Range of parameters (upper bound).
varmin = 0.01 ! the terminating limit for the variance of the function values. varmin>0 is required
step = (/0.001, 0.05, 0.01/) ! determines the size and shape of the initial simplex. The relative magnitudes of its elements should reflect the units of the variables. size(step)=size(start)
konvge = 10 ! the convergence check is carried out every konvge iterations
maxeval = 2000 ! the maximum number of function evaluations. default: 1000
! Call Nelder-Mead optimizer to reduce GCOMP
pmin = nelmin(obj_func, pstart, varmin, step, konvge, maxeval, &
funcmin, neval, numrestart, ierror, history)
! pmin = nelminrange(obj_func, pstart, prange, varmin, step, konvge, maxeval, &
! funcmin, neval, numrestart, ierror)
! scale up distance h
where (gamma(:,1) > 0._dp) gamma(:,1) = gamma(:,1) * gmax(1)
! scale back parameters (only for range a)
pmin(2)=pmin(2)*gmax(1)
print *, neval, ' of ', maxeval
print *, "funcmin: ", funcmin
print *, "p_obj: ", pmin
print *, 'error: ', ierror
print *, 'varmin: ', varmin
print *, 'history: ', history(1), history(size(history))
! stop 'Testing Nelder Mead algorithm'
! ! estimate statistics
! call stats
! !
! ! scale up distance h
! where (gamma(:,1) > 0._dp) gamma(:,1) = gamma(:,1) * gmax(1)
! ! scale back parameters (only for range a)
! beta(3)=beta(3)*gmax(1)
! !print*, 'opt co =', beta(1)
! !print*, 'opt c =', beta(2)
! !print*, 'opt a =', beta(3)
end subroutine OPTI
!*****************************************************************************
! Main_opti.for
! AUTHOR
! Luis E. Samaniego-Eguiguren, IREUS, 28.05.1999
! DESCRIPTION
! Initialization of the Nonlinear Optimization Subroutine GRG2
! The function to be optimized is suplied in subroutine GCOMP
! DEFINITION OF INPUT/OUTPUT-CHANNELS
! SCREEN: *
!
! SCREEN OUTPUT: *
! OUTPUT FILE: 6
!******************************************************************************
subroutine OPTI(pmin)
use VarFit
use mo_kind, only: i4, dp
use mo_obj_func, only: obj_func
use mo_nelmin, only: nelmin
! parameters for Nelder-Mead algorithm
real(dp) :: pstart(3) ! Starting point for the iteration.
real(dp), intent(out) :: pmin(3) ! optimized parameter set
real(dp) :: prange(3, 2) ! Range of parameters (upper and lower bound).
real(dp) :: varmin ! the terminating limit for the variance of the function values. varmin>0 is required
real(dp) :: step(3) ! determines the size and shape of the initial simplex. The relative magnitudes of its elements should reflect the units of the variables. size(step)=size(start)
integer(i4) :: konvge ! the convergence check is carried out every konvge iterations
integer(i4) :: maxeval ! the maximum number of function evaluations. default: 1000
real(dp) :: funcmin
integer(i4) :: neval ! the number of function evaluations used.
integer(i4) :: numrestart ! the number of restarts.
integer(i4) :: ierror ! error indicator.
! 0: no errors detected.
! 1: varmin or konvge have an illegal value.
! 2: iteration terminated because maxeval was exceeded without convergence.
real(dp), allocatable :: history(:)
! ! inputs for GRG2
! IMPLICIT DOUBLE PRECISION(A-H,O-Z), INTEGER(I,J,L,M,N)
! INTEGER*4 NCORE,NNVARS,NFUN,MAXBAS,MAXHES,LASTZ
! INTEGER*4 M,N,MP1,NPMP1,NBMAX,NNBMAX,NPNBMX,MPNBMX,NRTOT
! LOGICAL MAXIM,INPRNT,OTPRNT
! DIMENSION BLVAR(100),BUVAR(100),BLCON(10),BUCON(10)
! DIMENSION RAMCON(10),RAMVAR(10),INBIND(10),Z(5000)
! DIMENSION NONBAS(10),REDGR(10),DEFAUL(20),TITLE(19)
! DIMENSION XX(100),FCNS(10),RMULTS(10)
! DATA BIG/1.D31/
! Initialization of Nelder-Mead
pstart = (/0.05, 10., 0.5/) ! Starting point for the iteration.
prange(:, 1) = (/0., 0., 0./) ! Range of parameters (lower bound).
prange(:, 2) = (/2., 10., 2./) ! Range of parameters (upper bound).
varmin = 0.001 ! the terminating limit for the variance of the function values. varmin>0 is required
step = (/0.001, 1., 0.01/) ! determines the size and shape of the initial simplex. The relative magnitudes of its elements should reflect the units of the variables. size(step)=size(start)
konvge = 100 ! the convergence check is carried out every konvge iterations
maxeval = 2000 ! the maximum number of function evaluations. default: 1000
! Call Nelder-Mead optimizer to reduce GCOMP
pmin = nelmin(obj_func, pstart, varmin, step, konvge, maxeval, &
funcmin, neval, numrestart, ierror, history)
! pmin = nelminrange(obj_func, pstart, prange, varmin, step, konvge, maxeval, &
! funcmin, neval, numrestart, ierror)
! scale up distance h
where (gamma(:,1) > 0._dp) gamma(:,1) = gamma(:,1) * gmax(1)
! scale back parameters (only for range a)
pmin(3)=pmin(3)*gmax(1)
beta = pmin
print *, '=============================='
print *, ' Results of Nelder-Mead optimization '
print *, neval, ' of ', maxeval
print *, "funcmin: ", funcmin
print *, "p_obj: ", pmin
print *, 'error: ', ierror
print *, 'varmin: ', varmin
print *, 'history: ', history(1), history(size(history))
print *, 'gmax: ', gmax(1)
! estimate statistics
call stats
!
! scale up distance h
where (gamma(:,1) > 0._dp) gamma(:,1) = gamma(:,1) * gmax(1)
end subroutine OPTI
......@@ -10,13 +10,14 @@
! Created Sa 19.02.2004 main structure
! Last update 12.04.2006
!**********************************************************************************
subroutine setVario
subroutine setVario(param)
use runControl
use VarFit
use mainVar
use mo_kind, only: i4, dp
implicit none
integer(i4) :: jd, y, y0
real(dp), intent(out) :: param(3)
integer(i4) :: jd, y, y0
!
! Estimation
if (flagVario) then
......@@ -37,12 +38,10 @@ subroutine setVario
! variance
v0 = v0 / real(nobs, dp)
! optimize parameters
open(UNIT = 6,FILE = 'Report_OPT.sol', STATUS='REPLACE')
! open(UNIT = 6,FILE = 'Report_OPT.sol', STATUS='REPLACE')
print *, "ST: replace old GRG2 opti with something better"
! call opti
! fit a theoretical variogram
call opti(param)
end if
end subroutine setVario
! **************************************************************************
!
! SUBROUTINE statistics
!
! **************************************************************************
subroutine stats
use varFit, only : E, beta, gamma, nbins
use mo_kind, only : i4, dp
implicit none
integer(i4), parameter :: incx = 1
integer(i4) :: k
integer(i4) :: ne
real(dp) :: SSE
real(dp) :: zObsMean, zCalMean
real(dp) :: zObsVar, zCalVar, sumP, NSE_denom
real(dp), dimension(:), allocatable :: error, denom
real(dp), dimension(:), allocatable :: zCal, zObs
real(dp), parameter :: small = -9.999d3
real(dp) :: tVar
!
!Initialize
ne = count(gamma(:,1) > 0.0_dp)
allocate (error(ne), denom(ne), zCal(ne), zObs(ne))
zObs = gamma(1:ne,2)
zCal = 0.0_dp
do k=1,ne
if (gamma(k,1) > 0._dp) zCal(k) = tvar(gamma(k,1),beta(1),beta(2),beta(3))
end do
!
error = zObs-zCal
if ( ne > 1 ) then
zObsMean = sum(zObs)/real(ne, dp)
zCalMean = sum(zCal)/real(ne, dp)
sumP = dot_product(zObs,zCal)
zObsVar = dot_product(zObs,zObs) - real(ne, dp) * zObsMean * zObsMean
zCalVar = dot_product(zCal,zCal) - real(ne, dp) * zCalMean * zCalMean
SSE = dot_product(error,error)
denom = zObs - zObsMean
NSE_denom = dot_product(denom,denom)
else
zObsMean = small
zCalMean = small
zObsVar = small
zCalVar = small
end if
! ****************
! Quality measures
! ****************
if ( ne > 0 ) then
!
! BIAS
E(1) = zCalMean - zObsMean
!
! MSE
E(2) = SSE/real(ne, dp)
!
! RMSE
if ( E(2) > 0.0_dp ) then
E(3) = dsqrt(E(2))
else
E(3) = small
end if
!
! RRMSE
if ( E(3) > 0.0_dp ) then
E(4)=E(3)/zObsMean
else
E(4)= small
end if
!
! MAE
E(5)= sum(abs(error))
E(5)= E(5)/real(ne, dp)
!
! RMAE
E(6)=E(5)/zObsMean
!
! r
E(7)= (sumP-real(ne, dp) * zCalMean * zObsMean) / dsqrt(zCalVar * zObsVar)
!
! NSE
E(8)= 1.0_dp - (SSE/NSE_denom)
else
E = small
end if
deallocate (error, denom, zCal, zObs)
end subroutine stats
! **************************************************************************
!
! SUBROUTINE statistics
!
! **************************************************************************
subroutine stats
use varFit, only : E, beta, gamma, nbins
use mo_kind, only : i4, dp
implicit none
integer(i4), parameter :: incx = 1
integer(i4) :: k
integer(i4) :: ne
real(dp) :: SSE
real(dp) :: zObsMean, zCalMean
real(dp) :: zObsVar, zCalVar, sumP, NSE_denom
real(dp), dimension(:), allocatable :: error, denom
real(dp), dimension(:), allocatable :: zCal, zObs
real(dp), parameter :: small = -9.999d3
real(dp) :: tVar
!
!Initialize
ne = count(gamma(:,1) > 0.0_dp)
allocate (error(ne), denom(ne), zCal(ne), zObs(ne))
zObs = gamma(1:ne,2)
zCal = 0.0_dp
do k=1,ne
if (gamma(k,1) > 0._dp) zCal(k) = tvar(gamma(k,1),beta(1),beta(2),beta(3))
end do
!
error = zObs-zCal
if ( ne > 1 ) then
zObsMean = sum(zObs)/real(ne, dp)
zCalMean = sum(zCal)/real(ne, dp)
sumP = dot_product(zObs,zCal)
zObsVar = dot_product(zObs,zObs) - real(ne, dp) * zObsMean * zObsMean
zCalVar = dot_product(zCal,zCal) - real(ne, dp) * zCalMean * zCalMean
SSE = dot_product(error,error)
denom = zObs - zObsMean
NSE_denom = dot_product(denom,denom)
else
zObsMean = small
zCalMean = small
zObsVar = small
zCalVar = small
end if
! ****************
! Quality measures
! ****************
if ( ne > 0 ) then
!
! BIAS
E(1) = zCalMean - zObsMean
print *, 'BIAS: ', E(1)
!
! MSE
E(2) = SSE/real(ne, dp)
print *, 'MSE: ', E(2)
!
! RMSE
if ( E(2) > 0.0_dp ) then
E(3) = dsqrt(E(2))
else
E(3) = small
end if
print *, 'RMSE: ', E(3)
!
! RRMSE
if ( E(3) > 0.0_dp ) then
E(4)=E(3)/zObsMean
else
E(4)= small
end if
print *, 'PRMSE:', E(4)
!
! MAE
E(5)= sum(abs(error))
E(5)= E(5)/real(ne, dp)
print *, 'MAE: ', E(5)
!
! RMAE
E(6)=E(5)/zObsMean
print *, 'RMAE: ', E(6)
!
! r
E(7)= (sumP-real(ne, dp) * zCalMean * zObsMean) / dsqrt(zCalVar * zObsVar)
print *, 'r: ', E(7)
!
! NSE
E(8)= 1.0_dp - (SSE/NSE_denom)
print *, 'NSE: ', E(8)
else
E = small
end if
deallocate (error, denom, zCal, zObs)
end subroutine stats
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment