From 7a04203c31554dbc0087f61a9e2769a3d4d66af3 Mon Sep 17 00:00:00 2001 From: CyrilThebault Date: Wed, 29 Jul 2026 14:38:41 +1200 Subject: [PATCH 1/3] Add parameter transformation utilities --- build/FUSE_SRC/util/parameter_transform.f90 | 260 ++++++++++++++++++++ build/Makefile | 1 + 2 files changed, 261 insertions(+) create mode 100644 build/FUSE_SRC/util/parameter_transform.f90 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 From 30166e841c655d2b0e99a41d9acb909718d29cdc Mon Sep 17 00:00:00 2001 From: CyrilThebault Date: Wed, 29 Jul 2026 15:38:04 +1200 Subject: [PATCH 2/3] Integrate parameter transformations into SCE calibration --- build/FUSE_SRC/driver/functn.f90 | 38 ++++-- .../FUSE_SRC/driver/sce_callback_context.f90 | 14 ++- build/FUSE_SRC/driver/sce_driver.f90 | 118 +++++++++++++++--- 3 files changed, 147 insertions(+), 23 deletions(-) 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 From d5d32432efd50a699daea65a42cd0912bccda9e4 Mon Sep 17 00:00:00 2001 From: CyrilThebault Date: Tue, 4 Aug 2026 09:42:03 +1200 Subject: [PATCH 3/3] Refactor SCE parameter transformations --- build/FUSE_SRC/driver/sce_driver.f90 | 189 +++++++++++--------- build/FUSE_SRC/util/parameter_transform.f90 | 53 ++---- 2 files changed, 113 insertions(+), 129 deletions(-) diff --git a/build/FUSE_SRC/driver/sce_driver.f90 b/build/FUSE_SRC/driver/sce_driver.f90 index d2d5d19..fb8702f 100644 --- a/build/FUSE_SRC/driver/sce_driver.f90 +++ b/build/FUSE_SRC/driver/sce_driver.f90 @@ -19,12 +19,6 @@ subroutine sce_driver(info, work, domain, APAR, BL, BU) 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 @@ -43,14 +37,8 @@ subroutine sce_driver(info, work, domain, APAR, BL, BU) 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 :: 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 @@ -74,82 +62,15 @@ subroutine sce_driver(info, work, domain, APAR, BL, BU) INIFLG = 1 ! 1 = include initial point in the population IPRINT = 1 ! 0 = supress printing - ! 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 + call setup_parameter_transforms( APAR, BL, BU, & + APAR_MSP, BL_MSP, BU_MSP, & + TRANSFORM_CODES, & + IERR, MESSAGE ) + + if (IERR /= 0) then + write(*,'(A)') trim(MESSAGE) + stop 'Unable to set up parameter transformations' + 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 @@ -181,9 +102,99 @@ subroutine sce_driver(info, work, domain, APAR, BL, BU) ! 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 + ! ------------------------------------------------------------------------ + ! Prepare parameter transformations for SCE optimization. + ! ------------------------------------------------------------------------ + subroutine setup_parameter_transforms( apar_phys, bl_phys, bu_phys, & + apar_search, bl_search, bu_search, & + transform_codes, ierr, message) + + use multiparam, only: NUMPAR, 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 + + implicit none + + real(wp), intent(in) :: apar_phys(:) + real(wp), intent(in) :: bl_phys(:) + real(wp), intent(in) :: bu_phys(:) + + real(MSP), allocatable, intent(out) :: apar_search(:) + real(MSP), allocatable, intent(out) :: bl_search(:) + real(MSP), allocatable, intent(out) :: bu_search(:) + + integer(I4B), allocatable, intent(out) :: transform_codes(:) + + integer(I4B), intent(out) :: ierr + character(len=*), intent(out) :: message + + real(MSP), allocatable :: apar_phys_msp(:) + real(MSP), allocatable :: bl_phys_msp(:) + real(MSP), allocatable :: bu_phys_msp(:) + + type(PARATT) :: param_meta + integer(I4B) :: ipar + + ierr = 0 + message = '' + + allocate(apar_phys_msp(NUMPAR)) + allocate(bl_phys_msp(NUMPAR)) + allocate(bu_phys_msp(NUMPAR)) + + allocate(apar_search(NUMPAR)) + allocate(bl_search(NUMPAR)) + allocate(bu_search(NUMPAR)) + + allocate(transform_codes(NUMPAR)) + + apar_phys_msp = apar_phys + bl_phys_msp = bl_phys + bu_phys_msp = bu_phys + + 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 + message = "Parameter validation: "//trim(message) + return + end if + + end do + + call vector_to_search_space(apar_phys_msp, transform_codes, apar_search, ierr, message) + + if (ierr /= 0) then + message = 'Initial parameter set: '//trim(message) + return + end if + + call vector_to_search_space(bl_phys_msp, transform_codes, bl_search, ierr, message) + + if (ierr /= 0) then + message = 'Lower parameter bounds: '//trim(message) + return + end if + + call vector_to_search_space(bu_phys_msp, transform_codes, bu_search, ierr, message) + + if (ierr /= 0) then + message = 'Upper parameter bounds: '//trim(message) + return + end if + + end subroutine setup_parameter_transforms + end module sce_driver_MODULE diff --git a/build/FUSE_SRC/util/parameter_transform.f90 b/build/FUSE_SRC/util/parameter_transform.f90 index c3150fd..50c6880 100644 --- a/build/FUSE_SRC/util/parameter_transform.f90 +++ b/build/FUSE_SRC/util/parameter_transform.f90 @@ -92,8 +92,7 @@ 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) + 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(:) @@ -119,30 +118,16 @@ subroutine vector_to_search_space(physical_values, transform_codes, search_value 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 (TRANS_NONE, TRANS_LOG10, TRANS_LN) + search_values(i) = to_search_space(physical_values(i), transform_codes(i)) case default ierr = 12 - write(message,'(a,i0,a,i0)') & - 'Unknown transformation code ', transform_codes(i), & - ' for parameter ', i + 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 @@ -151,8 +136,7 @@ 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) + 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(:) @@ -179,14 +163,11 @@ subroutine vector_to_physical_space(search_values, transform_codes, physical_val select case (transform_codes(i)) case (TRANS_NONE, TRANS_LOG10, TRANS_LN) - physical_values(i) = & - to_physical_space(search_values(i), transform_codes(i)) + 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 + write(message,'(a,i0,a,i0)') 'Unknown transformation code ', transform_codes(i), ' for parameter ', i return end select @@ -199,8 +180,7 @@ 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) + subroutine validate_transform(parname, lower, default_value, upper, transform_code, ierr, message) character(len=*), intent(in) :: parname real(MSP), intent(in) :: lower @@ -215,15 +195,13 @@ subroutine validate_transform(parname, lower, default_value, upper, & if (lower > upper) then ierr = 1 - write(message,'(a,a)') & - 'Invalid parameter bounds for ', trim(parname) + 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) + write(message,'(a,a)') 'Default value outside bounds for ', trim(parname) return end if @@ -234,13 +212,10 @@ subroutine validate_transform(parname, lower, default_value, upper, & case (TRANS_LOG10, TRANS_LN) - if (lower <= 0.0_MSP .or. default_value <= 0.0_MSP .or. & - upper <= 0.0_MSP) then + 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), '.' + write(message,'(a,a,a)') 'Logarithmic transformation requires strictly positive values for ', trim(parname) return end if @@ -248,9 +223,7 @@ subroutine validate_transform(parname, lower, default_value, upper, & case default ierr = 4 - write(message,'(a,i0,a,a)') & - 'Unknown transformation code ', transform_code, & - ' for parameter ', trim(parname) + write(message,'(a,i0,a,a)') 'Unknown transformation code ', transform_code, ' for parameter ', trim(parname) return end select