Merge pull request #2940 from GEOS-ESM/feature/mathomp4/2926-update-q…

…uantize

Fixes #2926. Update quantization properties

Apps/Regrid_Util.F90

@@ -96,7 +96,7 @@ subroutine process_command_line(this,rc)
     this%lat_range=uninit
     this%shave=64
     this%deflate=0
-    this%quantize_algorithm=1
+    this%quantize_algorithm=0
     this%quantize_level=0
     this%use_weights = .false.
     nargs = command_argument_count()

CHANGELOG.md

@@ -37,6 +37,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
     - Update `esma_add_fortran_submodules` function
     - Move MPI detection out of FindBaselibs
     - Add SMOD to submodule generator
+- Add support for preliminary CF Conventions quantization properties
+  - Add new quantization keyword `granular_bitround` to History. This will be the preferred keyword for quantization in the future
+    replacing `GranularBR`
 
 ### Fixed
 
@@ -47,6 +50,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Deprecated
 
+- Deprecate `GranularBR` as a quantization method keyword in History. We will prefer `granular_bitround` in the future to match
+  draft CF conventions. This will be removed in MAPL 3.
+
 ## [2.47.1] - 2024-07-17
 
 ### Fixed

gridcomps/History/MAPL_HistoryGridComp.F90

@@ -425,6 +425,7 @@ subroutine Initialize ( gc, import, dumexport, clock, rc )
     logical :: has_conservative_keyword, has_regrid_keyword
     integer :: create_mode
     character(len=:), allocatable :: uppercase_algorithm
+    character(len=2) :: tmpchar
 
 ! Begin
 !------
@@ -838,40 +839,51 @@ subroutine Initialize ( gc, import, dumexport, clock, rc )
        call ESMF_ConfigGetAttribute ( cfg, list(n)%quantize_algorithm_string, default='NONE', &
                                       label=trim(string) // 'quantize_algorithm:' ,_RC )
 
+       call ESMF_ConfigGetAttribute ( cfg, list(n)%quantize_level, default=0, &
+                                      label=trim(string) // 'quantize_level:' ,_RC )
+
        ! Uppercase the algorithm string just to allow for any case
+       ! CF Conventions will prefer 'bitgroom', 'bitround', and 'granular_bitround'
+       ! but we will allow 'GranularBR' in MAPL2, deprecate it, and remove it in MAPL3
        uppercase_algorithm = ESMF_UtilStringUpperCase(list(n)%quantize_algorithm_string,_RC)
        select case (trim(uppercase_algorithm))
        case ('NONE')
-          list(n)%quantize_algorithm = MAPL_Quantize_Disabled
+          list(n)%quantize_algorithm = MAPL_NOQUANTIZE
+          ! If quantize_algorithm is 0, then quantize_level must be 0
+          _ASSERT( list(n)%quantize_level == 0, 'quantize_algorithm is none, so quantize_level must be "none"')
        case ('BITGROOM')
-          list(n)%quantize_algorithm = MAPL_Quantize_BitGroom
-       case ('GRANULARBR')
-          list(n)%quantize_algorithm = MAPL_Quantize_GranularBR
+          list(n)%quantize_algorithm = MAPL_QUANTIZE_BITGROOM
+       case ('GRANULARBR', 'GRANULAR_BITROUND')
+          list(n)%quantize_algorithm = MAPL_QUANTIZE_GRANULAR_BITROUND
        case ('BITROUND')
-          list(n)%quantize_algorithm = MAPL_Quantize_BitRound
+          list(n)%quantize_algorithm = MAPL_QUANTIZE_BITROUND
        case default
-          _FAIL('Invalid quantize_algorithm. Allowed values are NONE, BitGroom, GranularBR, BitRound')
+          _FAIL('Invalid quantize_algorithm. Allowed values are none, bitgroom, granular_bitround, granularbr (deprecated), and bitround')
        end select
 
-       call ESMF_ConfigGetAttribute ( cfg, list(n)%quantize_level, default=0, &
-                                      label=trim(string) // 'quantize_level:' ,_RC )
-
        ! If nbits_to_keep < MAPL_NBITS_UPPER_LIMIT (24) and quantize_algorithm greater than 0, then a user might be doing different
        ! shaving algorithms. We do not allow this
-       _ASSERT( .not. ( (list(n)%nbits_to_keep < MAPL_NBITS_UPPER_LIMIT) .and. (list(n)%quantize_algorithm > 0) ), 'nbits < 24 and quantize_algorithm > 0 is not allowed. Choose one bit grooming method.')
-
-       ! quantize_algorithm must be between 0 and 3 where 0 means not enabled
-       _ASSERT( (list(n)%quantize_algorithm >= 0) .and. (list(n)%quantize_algorithm <= 3), 'quantize_algorithm must be between 0 and 3, where 0 means not enabled')
+       _ASSERT( .not. ( (list(n)%nbits_to_keep < MAPL_NBITS_UPPER_LIMIT) .and. (list(n)%quantize_algorithm > MAPL_NOQUANTIZE) ), 'nbits < 24 and quantize_algorithm not "none" is not allowed. Choose a supported quantization method.')
 
        ! Now we test in the case that a valid quantize algorithm is chosen
-       if (list(n)%quantize_algorithm == 0) then
-         ! If quantize_algorithm is 0, then quantize_level must be 0
-          _ASSERT( list(n)%quantize_level == 0, 'quantize_algorithm is 0, so quantize_level must be 0')
-       else
+       if (list(n)%quantize_algorithm /= MAPL_NOQUANTIZE) then
          ! If quantize_algorithm is greater than 0, then quantize_level must be greater than or equal to 0
          _ASSERT( list(n)%quantize_level >= 0, 'netCDF quantize has been enabled, so quantize_level must be greater than or equal to 0')
        end if
 
+       ! If a user has chosen MAPL_QUANTIZE_BITROUND, then we allow a maximum of 23 bits to be kept
+       if (list(n)%quantize_algorithm == MAPL_QUANTIZE_BITROUND) then
+          write(tmpchar, '(I2)') MAPL_QUANTIZE_MAX_NSB
+         _ASSERT( list(n)%quantize_level <= MAPL_QUANTIZE_MAX_NSB, 'netCDF bitround has been enabled, so number of significant bits (quantize_level) must be less than or equal to ' // trim(tmpchar))
+       end if
+
+       ! For MAPL_QUANTIZE_GRANULAR_BITROUND and MAPL_QUANTIZE_BITGROOM, these use number of
+       ! significant digits, so for single precision, we allow a maximum of 7 digits to be kept
+       if (list(n)%quantize_algorithm == MAPL_QUANTIZE_GRANULAR_BITROUND .or. list(n)%quantize_algorithm == MAPL_QUANTIZE_BITGROOM) then
+          write(tmpchar, '(I2)') MAPL_QUANTIZE_MAX_NSD
+         _ASSERT( list(n)%quantize_level <= MAPL_QUANTIZE_MAX_NSD, 'netCDF granular bitround or bitgroom has been enabled, so number of significant digits (quantize_level) must be less than or equal to ' // trim(tmpchar))
+       end if
+
        tm_default = -1
        call ESMF_ConfigGetAttribute ( cfg, list(n)%tm, default=tm_default, &
                                       label=trim(string) // 'tm:', _RC )

griddedio/GriddedIO.F90

@@ -26,6 +26,7 @@ module MAPL_GriddedIOMod
   use, intrinsic :: ISO_C_BINDING
   use, intrinsic :: iso_fortran_env, only: REAL64
   use ieee_arithmetic, only: isnan => ieee_is_nan
+  use netcdf, only: nf90_inq_libvers
   implicit none
 
   private
@@ -51,7 +52,7 @@ module MAPL_GriddedIOMod
      type(VerticalData) :: vdata
      type(GriddedIOitemVector) :: items
      integer :: deflateLevel = 0
-     integer :: quantizeAlgorithm = 1
+     integer :: quantizeAlgorithm = MAPL_NOQUANTIZE
      integer :: quantizeLevel = 0
      integer, allocatable :: chunking(:)
      logical :: itemOrderAlphabetical = .true.
@@ -60,6 +61,7 @@ module MAPL_GriddedIOMod
      contains
         procedure :: CreateFileMetaData
         procedure :: CreateVariable
+        procedure :: CreateQuantizationInfo
         procedure :: modifyTime
         procedure :: modifyTimeIncrement
         procedure :: bundlePost
@@ -190,19 +192,23 @@ subroutine CreateFileMetaData(this,items,bundle,timeInfo,vdata,ogrid,global_attr
         call this%timeInfo%add_time_to_metadata(this%metadata,rc=status)
         _VERIFY(status)
 
-        iter = this%items%begin()
         if (.not.allocated(this%chunking)) then
            call this%set_default_chunking(rc=status)
            _VERIFY(status)
         else
            call this%check_chunking(this%vdata%lm,_RC)
         end if
 
-
         order = this%metadata%get_order(rc=status)
         _VERIFY(status)
         metadataVarsSize = order%size()
 
+        ! If quantize algorithm is set, create a quantization_info variable
+        if (this%quantizeAlgorithm /= MAPL_NOQUANTIZE) then
+           call this%CreateQuantizationInfo(_RC)
+        end if
+
+        iter = this%items%begin()
         do while (iter /= this%items%end())
            item => iter%get()
            if (item%itemType == ItemTypeScalar) then
@@ -423,6 +429,31 @@ subroutine CreateVariable(this,itemName,rc)
 #else
         call v%add_attribute('regrid_method', regrid_method_int_to_string(this%regrid_method))
 #endif
+        ! The CF Convention will soon support quantization. This requires three new attributes
+        ! if enabled:
+        ! 1. quantization --> Will point to a quantization_info container with the quantization algorithm
+        !                     (NOTE: this will need to be programmatic when per-variable quantization is enabled)
+        ! 2a. quantization_nsb --> Number of significant bits (only for bitround)
+        ! 2b. quantization_nsd --> Number of significant digits (only for bitgroom and granular_bitround)
+        ! 3. quantization_maximum_relative_error --> Maximum relative error (defined as 2^(-nsb) for bitround, and UNDEFINED? for bitgroom and granular_bitround)
+
+        ! Bitround
+        if (this%quantizeAlgorithm == MAPL_QUANTIZE_BITROUND) then
+           call v%add_attribute('quantization', 'quantization_info')
+           call v%add_attribute('quantization_nsb', this%quantizeLevel)
+           call v%add_attribute('quantization_maximum_relative_error', 0.5 * 2.0**(-this%quantizeLevel))
+        end if
+        ! granular_bitround and bitgroom
+        if (this%quantizeAlgorithm == MAPL_QUANTIZE_BITGROOM .or. this%quantizeAlgorithm == MAPL_QUANTIZE_GRANULAR_BITROUND) then
+           call v%add_attribute('quantization', 'quantization_info')
+           call v%add_attribute('quantization_nsd', this%quantizeLevel)
+           ! Per czender, these have maximum_absolute_error. We use the calculate_mae function below
+           ! which replicates a table in doi 10.5194/gmd-12-4099-2019
+           ! NOTE: This might not be the right formula. As the CF Convention draft is updated,
+           ! we will update this code.
+           call v%add_attribute('quantization_maximum_absolute_error', calculate_mae(this%quantizeLevel))
+        end if
+
         call factory%append_variable_metadata(v)
         call this%metadata%add_variable(trim(varName),v,rc=status)
         _VERIFY(status)
@@ -442,6 +473,81 @@ subroutine CreateVariable(this,itemName,rc)
 
      end subroutine CreateVariable
 
+     function calculate_mae(nsd) result(mae)
+
+        ! This function is based on Table 3 of doi 10.5194/gmd-12-4099-2019
+        ! The algorithm is weird, but it does duplicate the table
+
+        integer, intent(in) :: nsd
+        real(kind=REAL32) :: mae
+        real(kind=REAL32) :: mae_base
+        integer :: correction
+
+        mae_base = 4.0 * (1.0/16.0)**floor(real(nsd)/2.0) * (1.0/8.0)**ceiling(real(nsd)/2.0)
+
+        correction = 1
+        if ( (nsd > 2 .and. mod(nsd, 2) == 0) .or. nsd == 7 ) then
+           correction = 2
+        end if
+
+        mae = mae_base * correction
+     end function calculate_mae
+
+     subroutine CreateQuantizationInfo(this,rc)
+        class (MAPL_GriddedIO), intent(inout) :: this
+        integer, optional, intent(out) :: rc
+
+        integer :: status
+
+        class (AbstractGridFactory), pointer :: factory
+        character(len=:), allocatable :: varName, netcdf_version
+        type(Variable) :: v
+
+        factory => get_factory(this%output_grid,_RC)
+
+        v = Variable(type=PFIO_CHAR)
+
+        ! In the future when we can do per variable quantization, we will need
+        ! to do things like quantization_info1, quantization_info2, etc.
+        ! For now, we will just use quantization_info as it is per collection
+        varName = "quantization_info"
+
+        ! We need to convert the quantization algorithm to a string
+        select case (this%quantizeAlgorithm)
+        case (MAPL_QUANTIZE_BITGROOM)
+           call v%add_attribute('algorithm', 'bitgroom')
+        case (MAPL_QUANTIZE_BITROUND)
+           call v%add_attribute('algorithm', 'bitround')
+        case (MAPL_QUANTIZE_GRANULAR_BITROUND)
+           call v%add_attribute('algorithm', 'granular_bitround')
+        case default
+           _FAIL('Unknown quantization algorithm')
+        end select
+
+        ! Next add the implementation details
+        ! 3. implementation: This property contains free-form text
+        ! that concisely conveys the algorithm provenance, including the
+        ! name of the library or client that performed the quantization,
+        ! the software version, and the name of the author(s) if deemed
+        ! relevant.
+        !
+        ! In the current case, all algorithms are from libnetcdf
+        ! we make a string using nf90_inq_libvers()
+
+        netcdf_version = 'libnetcdf ' // nf90_inq_libvers()
+        call v%add_attribute('implementation', netcdf_version)
+
+        ! NOTE: In the future if we add the MAPL bit-shaving
+        ! to use the quantization parts of the code, it will
+        ! need a different implementation string
+
+        call factory%append_variable_metadata(v)
+        call this%metadata%add_variable(trim(varName),v,_RC)
+
+        _RETURN(_SUCCESS)
+
+     end subroutine CreateQuantizationInfo
+
      subroutine modifyTime(this, oClients, rc)
         class(MAPL_GriddedIO), intent(inout) :: this
         type (ClientManager), optional, intent(inout) :: oClients

pfio/Variable.F90

@@ -28,7 +28,7 @@ module pFIO_VariableMod
       type (StringAttributeMap) :: attributes
       type (UnlimitedEntity) :: const_value
       integer :: deflation = 0 ! default no compression
-      integer :: quantize_algorithm = 1 ! default bitgroom
+      integer :: quantize_algorithm = 0 ! default no quantization
       integer :: quantize_level = 0 ! default no quantize_level
       integer, allocatable :: chunksizes(:)
    contains
@@ -85,7 +85,7 @@ function new_Variable(unusable, type, dimensions, chunksizes,const_value, deflat
 
       var%type = -1
       var%deflation = 0
-      var%quantize_algorithm = 1
+      var%quantize_algorithm = 0
       var%quantize_level = 0
       var%chunksizes = empty
       var%dimensions = StringVector()

shared/Constants/InternalConstants.F90

@@ -166,13 +166,17 @@ module MAPL_InternalConstantsMod
 
    integer, parameter :: MAPL_NBITS_NOT_SET = 1000
    integer, parameter :: MAPL_NBITS_UPPER_LIMIT = 24
-   ! Constants for netCDF quantize
+   ! Constants for netCDF quantize (these echo the values in the netcdf-fortran library)
    enum, bind(c)
-      enumerator MAPL_Quantize_Disabled
-      enumerator MAPL_Quantize_BitGroom
-      enumerator MAPL_Quantize_GranularBR
-      enumerator MAPL_Quantize_BitRound
+      enumerator MAPL_NOQUANTIZE
+      enumerator MAPL_QUANTIZE_BITGROOM
+      enumerator MAPL_QUANTIZE_GRANULAR_BITROUND
+      enumerator MAPL_QUANTIZE_BITROUND
    endenum
+   ! Maximum number of significant digits for quantization (bitgroom, granular_bitround)
+   integer, parameter :: MAPL_QUANTIZE_MAX_NSD = 7
+   ! Maximum number of significant bits for quantization (bitround)
+   integer, parameter :: MAPL_QUANTIZE_MAX_NSB = 23
    ! Constant masking
    enum, bind(c)
       enumerator MAPL_MASK_OUT