Skip to content

Support parameter transformations in SCE calibration #43

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
CyrilThebault wants to merge 2 commits into
base: staging
Choose a base branch
from
Open

Support parameter transformations in SCE calibration #43

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7a04203
Add parameter transformation utilities
CyrilThebault Jul 29, 2026
30166e8
Integrate parameter transformations into SCE calibration
CyrilThebault Jul 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 31 additions & 7 deletions build/FUSE_SRC/driver/functn.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,13 +4,15 @@ 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:
! --------
! Wrapper for SCE (used to compute the objective function)
! ---------------------------------------------------------------------------------------
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
Expand All @@ -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

Expand All @@ -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)

Expand Down
14 changes: 13 additions & 1 deletion build/FUSE_SRC/driver/sce_callback_context.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
118 changes: 103 additions & 15 deletions build/FUSE_SRC/driver/sce_driver.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Comment on lines +91 to +152

@martynpclark martynpclark Aug 1, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be better if these lines are moved to a separate subroutine. You can include the new subroutine within this module.


! 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)
Expand All @@ -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
Loading