diff --git a/build/FUSE_SRC/driver/functn.f90 b/build/FUSE_SRC/driver/functn.f90 index 9fc39a6..df4464f 100644 --- a/build/FUSE_SRC/driver/functn.f90 +++ b/build/FUSE_SRC/driver/functn.f90 @@ -4,6 +4,7 @@ FUNCTION FUNCTN(NOPT,A) ! -------- ! Martyn Clark, 2009 ! Modified by Cyril Thébault to allow different metrics as objective function, 2024 +! Modified by Cyril Thébault to allow parameter transformations, 7/2026 ! --------------------------------------------------------------------------------------- ! Purpose: ! -------- @@ -11,6 +12,7 @@ FUNCTION FUNCTN(NOPT,A) ! --------------------------------------------------------------------------------------- USE nrtype ! variable types, etc. USE sce_callback_context, only: ctx ! access FUSE data structures +USE parameter_transform_module, only: vector_to_physical_space USE fuse_evaluate_module, only: fuse_evaluate ! run model and compute the metric chosen as objective function USE multiforce, only: ncid_forc ! NetCDF forcing file ID USE fuse_fileManager,only:METRIC, TRANSFO ! metric and transformation requested in the filemanager @@ -22,10 +24,10 @@ FUNCTION FUNCTN(NOPT,A) REAL(MSP), DIMENSION(100), INTENT(IN) :: A ! model parameter set - can be bumped up to 100 elements ! internal -REAL(WP), DIMENSION(NOPT) :: SCE_PAR ! sce parameter set -INTEGER(I4B) :: IERR ! error code for allocate/deallocate -INTEGER(I4B) :: ERR ! error code for fuse_metric -CHARACTER(LEN=256) :: MESSAGE ! error message for fuse_metric +REAL(MSP), DIMENSION(NOPT) :: SCE_PAR_MSP ! physical parameters in SCE precision +REAL(WP), DIMENSION(NOPT) :: SCE_PAR ! physical parameters in FUSE precision +INTEGER(I4B) :: IERR ! transformation error code +CHARACTER(LEN=256) :: MESSAGE ! parameter transformation error message LOGICAL(LGT) :: OUTPUT_FLAG ! .TRUE. = write model time series REAL(WP) :: METRIC_VAL ! value of the metric chosen as objective function @@ -36,9 +38,31 @@ FUNCTION FUNCTN(NOPT,A) nFUSE_eval = nFUSE_eval + 1 -! get SCE parameter set -SCE_PAR(1:NOPT) = A(1:NOPT) ! convert from MSP used in SCE to WP used in FUSE -OUTPUT_FLAG=.FALSE. ! do not produce *runs.nc files only, param.nc files +! Convert the optimizer vector from search space back to physical space. +IF (.NOT. ALLOCATED(ctx%transform_codes)) THEN + STOP 'SCE parameter transformations are not initialized' +END IF + +IF (SIZE(ctx%transform_codes) /= NOPT) THEN + STOP 'Incorrect number of SCE parameter transformation codes' +END IF + +CALL vector_to_physical_space( & + A(1:NOPT), & + ctx%transform_codes, & + SCE_PAR_MSP, & + IERR, & + MESSAGE) + +IF (IERR /= 0) THEN + WRITE(*,'(A)') TRIM(MESSAGE) + STOP 'Unable to transform SCE parameters to physical space' +END IF + +! Convert from MSP used by SCE to WP used by FUSE. +SCE_PAR = SCE_PAR_MSP + +OUTPUT_FLAG = .FALSE. ! do not produce runs.nc files during calibration CALL FUSE_evaluate(SCE_PAR, ctx%info, ctx%work, ctx%domain, OUTPUT_FLAG, METRIC_VAL) diff --git a/build/FUSE_SRC/driver/sce_callback_context.f90 b/build/FUSE_SRC/driver/sce_callback_context.f90 index 1a74af6..e99efed 100644 --- a/build/FUSE_SRC/driver/sce_callback_context.f90 +++ b/build/FUSE_SRC/driver/sce_callback_context.f90 @@ -1,4 +1,5 @@ module sce_callback_context + use nrtype, only: I4B use info_types, only: fuse_info use data_types, only: domain_data use work_types, only: fuse_work @@ -11,23 +12,34 @@ module sce_callback_context type(fuse_info), pointer :: info => null() type(fuse_work), pointer :: work => null() type(domain_data), pointer :: domain => null() + integer(I4B), allocatable :: transform_codes(:) end type sce_context type(sce_context), save :: ctx contains - subroutine set_sce_context(info, work, domain) + subroutine set_sce_context(info, work, domain, transform_codes) type(fuse_info), target, intent(inout) :: info type(fuse_work), target, intent(inout) :: work type(domain_data), target, intent(inout) :: domain + integer(I4B), intent(in) :: transform_codes(:) ctx%info => info ctx%work => work ctx%domain => domain + + if (allocated(ctx%transform_codes)) deallocate(ctx%transform_codes) + + allocate(ctx%transform_codes(size(transform_codes))) + ctx%transform_codes = transform_codes + end subroutine subroutine clear_sce_context() nullify(ctx%info, ctx%work, ctx%domain) + if (allocated(ctx%transform_codes)) then + deallocate(ctx%transform_codes) + end if end subroutine end module sce_callback_context diff --git a/build/FUSE_SRC/driver/sce_driver.f90 b/build/FUSE_SRC/driver/sce_driver.f90 index 03898ca..d2d5d19 100644 --- a/build/FUSE_SRC/driver/sce_driver.f90 +++ b/build/FUSE_SRC/driver/sce_driver.f90 @@ -5,7 +5,7 @@ module sce_driver_MODULE use work_types, only: fuse_work use data_types, only: domain_data - use sce_callback_context, only: set_sce_context, clear_sce_context + use sce_callback_context, only: set_sce_context, clear_sce_context implicit none @@ -15,12 +15,19 @@ module sce_driver_MODULE contains subroutine sce_driver(info, work, domain, APAR, BL, BU) - USE multiparam, only: MAXN ! maximum number of trials before optimization is terminated + USE multiparam, only: MAXN ! maximum number of trials before optimization is terminated USE multiparam, only: KSTOP ! number of shuffling loops the value must change by PCENTO USE multiparam, only: PCENTO ! the percentage USE multiparam, only: NUMPAR ! # parameters + USE multiparam, only: LPARAM, PARATT + + USE GETPAR_STR_MODULE, only: GETPAR_STR + + USE parameter_transform_module, only: validate_transform + USE parameter_transform_module, only: vector_to_search_space + USE fuse_globaldata, only: isPrint ! used to turn of printing for calibration runs - USE fuse_globaldata, only: nFUSE_eval ! # FUSE evaluations + USE fuse_globaldata, only: nFUSE_eval ! # FUSE evaluations USE model_defn, only: FNAME_TEMPRY, FNAME_ASCII implicit none ! input/output @@ -35,7 +42,18 @@ subroutine sce_driver(info, work, domain, APAR, BL, BU) REAL(MSP), DIMENSION(:), ALLOCATABLE :: APAR_MSP ! ! lower bound of model parameters REAL(MSP), DIMENSION(:), ALLOCATABLE :: BL_MSP ! ! lower bound of model parameters REAL(MSP), DIMENSION(:), ALLOCATABLE :: BU_MSP ! ! upper bound of model parameters - REAL(MSP), DIMENSION(:), ALLOCATABLE :: URAND_MSP ! vector of quasi-random numbers U[0,1] + + REAL(MSP), DIMENSION(:), ALLOCATABLE :: APAR_PHYS_MSP + REAL(MSP), DIMENSION(:), ALLOCATABLE :: BL_PHYS_MSP + REAL(MSP), DIMENSION(:), ALLOCATABLE :: BU_PHYS_MSP + + INTEGER(I4B), DIMENSION(:), ALLOCATABLE :: TRANSFORM_CODES + + TYPE(PARATT) :: PARAM_META + INTEGER(I4B) :: IPAR + INTEGER(I4B) :: IERR + CHARACTER(LEN=256) :: MESSAGE + INTEGER(I4B) :: NOPT ! number of parameters to be optimized INTEGER(I4B) :: NGS ! # complexes in the initial population INTEGER(I4B) :: NPG ! # points in each complex @@ -45,8 +63,6 @@ subroutine sce_driver(info, work, domain, APAR, BL, BU) INTEGER(I4B) :: INIFLG ! 1 = include initial point in the population INTEGER(I4B) :: IPRINT ! 0 = supress printing INTEGER(I4B) :: ISCE ! unit number for SCE write - integer(i4b) :: NUMPSET ! number of parameter sets - REAL(MSP) :: FUNCTN ! function name for the model run INTEGER(KIND=4) :: ISEED ! seed for the random sequence NOPT = NUMPAR ! number of parameters to be optimized (NUMPAR in module multiparam) @@ -58,16 +74,86 @@ subroutine sce_driver(info, work, domain, APAR, BL, BU) INIFLG = 1 ! 1 = include initial point in the population IPRINT = 1 ! 0 = supress printing - NUMPSET=1.2*MAXN ! will be used to define the parameter set dimension of the NetCDF files - ! using 1.2MAXN since the final number of parameter sets produced by SCE is unknown - - ! convert from WP used in FUSE to MSP used in SCE - ALLOCATE(APAR_MSP(NUMPAR), BL_MSP(NUMPAR), BU_MSP(NUMPAR)) - APAR_MSP=APAR; BL_MSP=BL; BU_MSP=BU + ! Store physical-space values using the precision required by SCE. + ALLOCATE(APAR_PHYS_MSP(NUMPAR)) + ALLOCATE(BL_PHYS_MSP(NUMPAR)) + ALLOCATE(BU_PHYS_MSP(NUMPAR)) + + ALLOCATE(APAR_MSP(NUMPAR)) + ALLOCATE(BL_MSP(NUMPAR)) + ALLOCATE(BU_MSP(NUMPAR)) + + ALLOCATE(TRANSFORM_CODES(NUMPAR)) + + APAR_PHYS_MSP = APAR + BL_PHYS_MSP = BL + BU_PHYS_MSP = BU + + ! Retrieve and validate the transformation code for each model parameter. + DO IPAR = 1, NUMPAR + + CALL GETPAR_STR(LPARAM(IPAR)%PARNAME, PARAM_META) + + TRANSFORM_CODES(IPAR) = PARAM_META%PARVTN + + CALL validate_transform( & + LPARAM(IPAR)%PARNAME, & + BL_PHYS_MSP(IPAR), & + APAR_PHYS_MSP(IPAR), & + BU_PHYS_MSP(IPAR), & + TRANSFORM_CODES(IPAR), & + IERR, & + MESSAGE) + + IF (IERR /= 0) THEN + WRITE(*,'(A)') TRIM(MESSAGE) + STOP 'Invalid parameter transformation' + END IF + + END DO + + ! Transform the initial parameter set into optimizer search space. + CALL vector_to_search_space( & + APAR_PHYS_MSP, & + TRANSFORM_CODES, & + APAR_MSP, & + IERR, & + MESSAGE) + + IF (IERR /= 0) THEN + WRITE(*,'(A)') TRIM(MESSAGE) + STOP 'Unable to transform initial parameter set' + END IF + + ! Transform the lower parameter bounds into optimizer search space. + CALL vector_to_search_space( & + BL_PHYS_MSP, & + TRANSFORM_CODES, & + BL_MSP, & + IERR, & + MESSAGE) + + IF (IERR /= 0) THEN + WRITE(*,'(A)') TRIM(MESSAGE) + STOP 'Unable to transform lower parameter bounds' + END IF + + ! Transform the upper parameter bounds into optimizer search space. + CALL vector_to_search_space( & + BU_PHYS_MSP, & + TRANSFORM_CODES, & + BU_MSP, & + IERR, & + MESSAGE) + + IF (IERR /= 0) THEN + WRITE(*,'(A)') TRIM(MESSAGE) + STOP 'Unable to transform upper parameter bounds' + END IF ! pass the FUSE structures to the context setter ! NOTE: in sce_context_set, info/work/domain have the target attribute so can point to them - call set_sce_context(info, work, domain) + call set_sce_context(info, work, domain, TRANSFORM_CODES) ! open up ASCII output file ISCE = 96 ! (file unit) @@ -93,9 +179,11 @@ subroutine sce_driver(info, work, domain, APAR, BL, BU) ! nullify pointers in the context setter call clear_sce_context() - ! deallocate space for real32 vectors + ! deallocate SCE and transformation arrays DEALLOCATE(APAR_MSP, BL_MSP, BU_MSP) - + DEALLOCATE(APAR_PHYS_MSP, BL_PHYS_MSP, BU_PHYS_MSP) + DEALLOCATE(TRANSFORM_CODES) + end subroutine sce_driver end module sce_driver_MODULE diff --git a/build/FUSE_SRC/util/parameter_transform.f90 b/build/FUSE_SRC/util/parameter_transform.f90 new file mode 100644 index 0000000..c3150fd --- /dev/null +++ b/build/FUSE_SRC/util/parameter_transform.f90 @@ -0,0 +1,260 @@ +! --------------------------------------------------------------------------------------- +! Creator: +! -------- +! Cyril Thébault, 2026 +! --------------------------------------------------------------------------------------- +! Purpose: +! -------- +! Apply parameter transformations for optimization. +! +! Parameters remain in physical space inside FUSE. Optimizers may operate in a +! transformed search space. This module provides conversions between both spaces. +! --------------------------------------------------------------------------------------- + +module parameter_transform_module + + use nrtype, only : MSP, I4B + + implicit none + private + + ! Transformation codes used in the PARVTN column of zConstraints. + integer(I4B), parameter, public :: TRANS_NONE = 0 + integer(I4B), parameter, public :: TRANS_LOG10 = 1 + integer(I4B), parameter, public :: TRANS_LN = 2 + + public :: to_search_space + public :: to_physical_space + public :: vector_to_search_space + public :: vector_to_physical_space + public :: validate_transform + +contains + + ! ------------------------------------------------------------------------------------- + ! Transform one physical parameter to the optimizer search space. + ! ------------------------------------------------------------------------------------- + pure function to_search_space(x, transform_code) result(z) + + real(MSP), intent(in) :: x + integer(I4B), intent(in) :: transform_code + real(MSP) :: z + + select case (transform_code) + + case (TRANS_NONE) + z = x + + case (TRANS_LOG10) + z = log10(x) + + case (TRANS_LN) + z = log(x) + + case default + ! Unknown codes must be rejected by validate_transform before this function is used. + z = x + + end select + + end function to_search_space + + + ! ------------------------------------------------------------------------------------- + ! Transform one optimizer value back to the physical parameter space. + ! ------------------------------------------------------------------------------------- + pure function to_physical_space(z, transform_code) result(x) + + real(MSP), intent(in) :: z + integer(I4B), intent(in) :: transform_code + real(MSP) :: x + + select case (transform_code) + + case (TRANS_NONE) + x = z + + case (TRANS_LOG10) + x = 10.0_MSP**z + + case (TRANS_LN) + x = exp(z) + + case default + ! Unknown codes must be rejected by validate_transform before this function is used. + x = z + + end select + + end function to_physical_space + + + ! ------------------------------------------------------------------------------------- + ! Transform a complete parameter vector to the optimizer search space. + ! ------------------------------------------------------------------------------------- + subroutine vector_to_search_space(physical_values, transform_codes, search_values, & + ierr, message) + + real(MSP), intent(in) :: physical_values(:) + integer(I4B), intent(in) :: transform_codes(:) + real(MSP), intent(out) :: search_values(:) + integer(I4B), intent(out) :: ierr + character(len=*), intent(out) :: message + + integer(I4B) :: i + integer(I4B) :: npar + + ierr = 0 + message = '' + + npar = size(physical_values) + + if (size(transform_codes) /= npar .or. size(search_values) /= npar) then + ierr = 10 + message = 'Inconsistent array sizes in vector_to_search_space' + return + end if + + do i = 1, npar + + select case (transform_codes(i)) + + case (TRANS_NONE) + continue + + case (TRANS_LOG10, TRANS_LN) + + if (physical_values(i) <= 0.0_MSP) then + ierr = 11 + write(message,'(a,i0,a)') & + 'Cannot apply logarithmic transformation to parameter ', i, & + ': value must be strictly positive.' + return + end if + + case default + ierr = 12 + write(message,'(a,i0,a,i0)') & + 'Unknown transformation code ', transform_codes(i), & + ' for parameter ', i + return + + end select + + search_values(i) = to_search_space(physical_values(i), transform_codes(i)) + + end do + + end subroutine vector_to_search_space + + + ! ------------------------------------------------------------------------------------- + ! Transform a complete optimizer vector back to physical parameter space. + ! ------------------------------------------------------------------------------------- + subroutine vector_to_physical_space(search_values, transform_codes, physical_values, & + ierr, message) + + real(MSP), intent(in) :: search_values(:) + integer(I4B), intent(in) :: transform_codes(:) + real(MSP), intent(out) :: physical_values(:) + integer(I4B), intent(out) :: ierr + character(len=*), intent(out) :: message + + integer(I4B) :: i + integer(I4B) :: npar + + ierr = 0 + message = '' + + npar = size(search_values) + + if (size(transform_codes) /= npar .or. size(physical_values) /= npar) then + ierr = 20 + message = 'Inconsistent array sizes in vector_to_physical_space' + return + end if + + do i = 1, npar + + select case (transform_codes(i)) + + case (TRANS_NONE, TRANS_LOG10, TRANS_LN) + physical_values(i) = & + to_physical_space(search_values(i), transform_codes(i)) + + case default + ierr = 21 + write(message,'(a,i0,a,i0)') & + 'Unknown transformation code ', transform_codes(i), & + ' for parameter ', i + return + + end select + + end do + + end subroutine vector_to_physical_space + + + ! ------------------------------------------------------------------------------------- + ! Validate a transformation and its physical bounds. + ! ------------------------------------------------------------------------------------- + subroutine validate_transform(parname, lower, default_value, upper, & + transform_code, ierr, message) + + character(len=*), intent(in) :: parname + real(MSP), intent(in) :: lower + real(MSP), intent(in) :: default_value + real(MSP), intent(in) :: upper + integer(I4B), intent(in) :: transform_code + integer(I4B), intent(out) :: ierr + character(len=*), intent(out) :: message + + ierr = 0 + message = '' + + if (lower > upper) then + ierr = 1 + write(message,'(a,a)') & + 'Invalid parameter bounds for ', trim(parname) + return + end if + + if (default_value < lower .or. default_value > upper) then + ierr = 2 + write(message,'(a,a)') & + 'Default value outside bounds for ', trim(parname) + return + end if + + select case (transform_code) + + case (TRANS_NONE) + continue + + case (TRANS_LOG10, TRANS_LN) + + if (lower <= 0.0_MSP .or. default_value <= 0.0_MSP .or. & + upper <= 0.0_MSP) then + + ierr = 3 + write(message,'(a,a,a)') & + 'Logarithmic transformation requires strictly positive values for ', & + trim(parname), '.' + return + + end if + + case default + + ierr = 4 + write(message,'(a,i0,a,a)') & + 'Unknown transformation code ', transform_code, & + ' for parameter ', trim(parname) + return + + end select + + end subroutine validate_transform + +end module parameter_transform_module diff --git a/build/Makefile b/build/Makefile index 37b295e..8ba1ee6 100755 --- a/build/Makefile +++ b/build/Makefile @@ -242,6 +242,7 @@ FUSE_UTILMS += fuse_fileManager.f90 FUSE_UTILMS += alloc_domain.f90 FUSE_UTILMS += alloc_scratch.f90 FUSE_UTILMS += time_utils.f90 +FUSE_UTILMS += parameter_transform.f90 FUSE_UTILMS += metaoutput.f90 FUSE_UTILMS += metaparams.f90 FUSE_UTILMS += meta_stats.f90