feat(app): optionally plot unnormalized histograms

cloud-microphysics/app/gnuplot.inp → cloud-microphysics/app/plot-normalized-histograms.gnu

@@ -1,11 +1,11 @@
 set logscale y 10
 set key font ",20"
-set title "Normalized Input Tensor Histograms" font "Arial,20"
+set title "Input-Tensor Histograms Mapped to the Interval [0,1]" font "Arial,20"
 
 plot for [col=2:8] 'training_inputs_stats.plt' using 1:col with linespoints title columnheader
 pause mouse any "Plot 1 of 2. Press any key to continue.\n"
 
-set title "Normalized Output Tensor Histograms" font "Arial,20"
+set title "Output-Tensor Histograms Mapped to the Interval [0,1]" font "Arial,20"
 
 plot for [col=2:6] 'training_outputs_stats.plt' using 1:col with linespoints title columnheader
 pause mouse any "Plot 2 of 2. Press any key to exit.\n"

cloud-microphysics/app/plot-raw-histograms.gnu

@@ -0,0 +1,22 @@
+set logscale y 10
+set term png size 2600, 1200
+
+set output "input-tensor-histograms.png"
+set multiplot layout 3,3
+do for [name in "potential-temperature pressure qc qr qs qv temperature"] {
+  filename = name . ".plt"
+  set title sprintf("Histogram for %s",name)
+  plot filename using 1:2 with linespoints title columnheader
+}
+print "\n Output-tensor histograms have been written to the graphics file 'output-tensor-histograms.png'."
+
+unset multiplot
+
+set output "output-tensor-histograms.png"
+set multiplot layout 3,2
+do for [name in "dptdt dqcdt dqrdt dqsdt dqvdt"] {
+  filename = name . ".plt"
+  set title sprintf("Histogram for %s",name)
+  plot filename using 1:2 with linespoints title columnheader
+}
+print " Input-tensor histograms have been written to the graphics file 'input-tensor-histograms.png'.\n"

cloud-microphysics/app/tensor-statistics.f90

@@ -23,35 +23,41 @@ program tensor_statistics
     'Usage: ' // new_line('a') // new_line('a') // &
     './build/run-fpm.sh run tensor-statistics -- \' // new_line('a') // &
     '  --base <string> --bins <integer> \' // new_line('a') // &
-    '  [--start <integer>] [--end <integer>] [--stride <integer>]' // &
+    '  [--raw] [--start <integer>] [--end <integer>] [--stride <integer>]' // &
     new_line('a') // new_line('a') // &
-    'where angular brackets denote user-provided values and square brackets denote optional arguments.' // new_line('a') // &
-    'The presence of a file named "stop" halts execution gracefully.'
+    'where angular brackets denote user-provided values and square brackets denote optional arguments.' // new_line('a')
 
   integer(int64) t_start, t_finish, clock_rate
   integer num_bins, start_step, stride
   integer, allocatable :: end_step
   character(len=:), allocatable :: base_name
+  logical raw
 
   call system_clock(t_start, clock_rate)
-  call get_command_line_arguments(base_name, num_bins, start_step, end_step, stride)
-  call compute_histograms(base_name)
+  call get_command_line_arguments(base_name, num_bins, start_step, end_step, stride, raw)
+  call compute_histograms(base_name, raw)
   call system_clock(t_finish)
 
   print *,"System clock time: ", real(t_finish - t_start, real64)/real(clock_rate, real64)
-  print *,"---"
   print *
-  print *,"---> Please see the *.plt files for the tensor ranges and histograms. <---"
-  print *,"---> Execute `gnuplot app/gnuplot.inp` to graph the histograms.       <---"
+  print *,"______________________________________________________"
+  print *,"The *.plt files contain tensor ranges and histograms."
+  print *,"If you have gnuplot installed, please execute the"
+  print *,"following command to produce histogram visualizations:" 
   print *
-  print *,"______ tensor_statistics done_______"
+  associate(file_name => trim(merge("plot-raw-histograms       ", "plot-normalized-histograms", raw)) // ".gnu")
+   print*,"  gnuplot app/" // file_name
+  end associate
+  print *
+  print *,"_______________ tensor_statistics done________________"
 
 contains
 
-  subroutine get_command_line_arguments(base_name, num_bins, start_step, end_step, stride)
+  subroutine get_command_line_arguments(base_name, num_bins, start_step, end_step, stride, raw)
     character(len=:), allocatable, intent(out) :: base_name
     integer, intent(out) :: num_bins, start_step, stride
     integer, intent(out), allocatable :: end_step
+    logical, intent(out) :: raw
 
     ! local variables
     type(command_line_t) command_line
@@ -62,6 +68,7 @@ subroutine get_command_line_arguments(base_name, num_bins, start_step, end_step,
     start_string = command_line%flag_value("--start")
     end_string = command_line%flag_value("--end")
     stride_string = command_line%flag_value("--stride")
+    raw = command_line%argument_present(["--raw"])
 
     associate(required_arguments => len(base_name)/=0 .and. len(bins_string)/=0)
        if (.not. required_arguments) error stop usage 
@@ -88,8 +95,9 @@ subroutine get_command_line_arguments(base_name, num_bins, start_step, end_step,
 
   end subroutine get_command_line_arguments
 
-  subroutine compute_histograms(base_name)
+  subroutine compute_histograms(base_name, raw)
     character(len=*), intent(in) :: base_name
+    logical, intent(in) :: raw
 
     real, allocatable, dimension(:,:,:,:) :: &
       pressure_in , potential_temperature_in , temperature_in , &
@@ -99,6 +107,7 @@ subroutine compute_histograms(base_name)
       dpt_dt, dqv_dt, dqc_dt, dqr_dt, dqs_dt
     double precision, allocatable, dimension(:) :: time_in, time_out
     double precision, parameter :: tolerance = 1.E-07
+    type(histogram_t), allocatable :: histograms(:)
     type(ubounds_t), allocatable :: ubounds(:)
     integer, allocatable :: lbounds(:)
     integer t, t_end
@@ -125,24 +134,31 @@ subroutine compute_histograms(base_name)
            ubounds_t(ubound(pressure_in)), ubounds_t(ubound(temperature_in)) &
           ]
 
-        print *,"Calculating input tensor histograms and write to file"
-
-        associate(histograms => [ &
-          histogram_t(pressure_in, "pressure", num_bins) &
-         ,histogram_t(potential_temperature_in, '"potential temperature"', num_bins) &
-         ,histogram_t(temperature_in, "temperature", num_bins) &
-         ,histogram_t(qv_in, "qv", num_bins) &
-         ,histogram_t(qc_in, "qc", num_bins) &
-         ,histogram_t(qr_in, "qr", num_bins) &
-         ,histogram_t(qs_in, "qs", num_bins) &
-        ])
-          block
-            type(file_t) histograms_file
-
-            histograms_file = to_file(histograms)
-            call histograms_file%write_lines(string_t(base_name // "_inputs_stats.plt"))
-          end block
-        end associate
+        print *,"Calculating input tensor histograms and writing to file"
+
+        histograms = [ &
+          histogram_t(pressure_in, "pressure", num_bins, raw) &
+         ,histogram_t(potential_temperature_in, 'potential-temperature', num_bins, raw) &
+         ,histogram_t(temperature_in, "temperature", num_bins, raw) &
+         ,histogram_t(qv_in, "qv", num_bins, raw) &
+         ,histogram_t(qc_in, "qc", num_bins, raw) &
+         ,histogram_t(qr_in, "qr", num_bins, raw) &
+         ,histogram_t(qs_in, "qs", num_bins, raw) &
+        ]
+        block
+          type(file_t) histograms_file
+          integer h
+
+          if (raw) then
+            do h = 1, size(histograms)
+              histograms_file = to_file(histograms(h))
+              call histograms_file%write_lines(string_t(histograms(h)%variable_name() // ".plt"))
+            end do
+          else
+              histograms_file = to_file(histograms)
+              call histograms_file%write_lines(string_t(base_name // "_inputs_stats.plt"))
+          end if
+        end block
       end associate
     end associate
 
@@ -190,22 +206,29 @@ subroutine compute_histograms(base_name)
       call assert(.not. any(ieee_is_nan(dqr_dt)), ".not. any(ieee_is_nan(dqr_dt)")
       call assert(.not. any(ieee_is_nan(dqs_dt)), ".not. any(ieee_is_nan(dqs_dt)")
 
-      print *,"Calculating output tensor histograms"
-
-      associate(histograms => [ &
-         histogram_t(dpt_dt, "d(pt)/dt", num_bins) &
-        ,histogram_t(dqv_dt, "d(qv)/dt", num_bins) &
-        ,histogram_t(dqc_dt, "d(qc)/dt", num_bins) &
-        ,histogram_t(dqr_dt, "d(qr)/dt", num_bins) &
-        ,histogram_t(dqs_dt, "d(qs)/dt", num_bins) &
-      ])
-        block
-          type(file_t) histograms_file
-
-          histograms_file = to_file(histograms)
-          call histograms_file%write_lines(string_t(base_name // "_outputs_stats.plt"))
-        end block
-      end associate
+      print *,"Calculating output tensor histograms and writing to file"
+
+      histograms = [ &
+         histogram_t(dpt_dt, "dptdt", num_bins, raw) &
+        ,histogram_t(dqv_dt, "dqvdt", num_bins, raw) &
+        ,histogram_t(dqc_dt, "dqcdt", num_bins, raw) &
+        ,histogram_t(dqr_dt, "dqrdt", num_bins, raw) &
+        ,histogram_t(dqs_dt, "dqsdt", num_bins, raw) &
+      ]
+      block
+        type(file_t) histograms_file
+        integer h
+
+        if (raw) then
+          do h = 1, size(histograms)
+            histograms_file = to_file(histograms(h))
+            call histograms_file%write_lines(string_t(histograms(h)%variable_name() // ".plt"))
+          end do
+        else
+            histograms_file = to_file(histograms)
+            call histograms_file%write_lines(string_t(base_name // "_outputs_stats.plt"))
+        end if
+      end block
     end associate
   end subroutine
 

cloud-microphysics/src/histogram_m.f90

@@ -26,25 +26,35 @@ module histogram_m
 
   interface histogram_t
 
-    pure module function histogram_on_unit_interval(v, variable_name, num_bins) result(histogram)
+    pure module function construct(v, variable_name, num_bins, raw) result(histogram)
       implicit none
       real, intent(in) :: v(:,:,:,:)
       character(len=*), intent(in) :: variable_name
       integer, intent(in) :: num_bins
+      logical, intent(in) :: raw
       type(histogram_t) histogram
     end function
 
   end interface
 
-  interface
+  interface to_file
+
+    pure module function to_separate_file(histogram) result(file)
+      implicit none
+      type(histogram_t), intent(in) :: histogram
+      type(file_t) file
+    end function
 
-    !pure module function to_file(histograms) result(file)
-    module function to_file(histograms) result(file)
+    pure module function to_common_file(histograms) result(file)
       implicit none
       type(histogram_t), intent(in) :: histograms(:)
       type(file_t) file
     end function
 
+  end interface
+
+  interface
+
     pure module function num_bins(self) result(bins)
       implicit none
       class(histogram_t), intent(in) :: self

cloud-microphysics/src/histogram_s.f90

@@ -36,7 +36,11 @@
     frequency = self%frequency_(bin)
   end procedure
 
-  module procedure to_file
+  module procedure to_separate_file
+     file = to_common_file([histogram])
+  end procedure
+
+  module procedure to_common_file
     type(string_t), allocatable :: comments(:), columns(:)
     type(string_t) column_headings
 
@@ -54,7 +58,7 @@
 
       block
         integer h
-        column_headings = "bin" // .cat. [(string_t(histograms(h)%variable_name()) // "   ", h=1,num_histograms)]
+        column_headings = "    bin    " // .cat. [("    " // string_t(histograms(h)%variable_name()) // "    ", h=1,num_histograms)]
       end block
 
       associate(num_bins =>  histograms(1)%num_bins())
@@ -89,9 +93,8 @@ pure function normalize(x, x_min, x_max) result(x_normalized)
     x_normalized = (x - x_min)/(x_max - x_min)
   end function
 
-  module procedure histogram_on_unit_interval
+  module procedure construct
 
-    real, parameter :: v_mapped_min = 0., v_mapped_max = 1.
     integer, allocatable :: in_bin(:)
     integer i
 
@@ -104,22 +107,34 @@ pure function normalize(x, x_min, x_max) result(x_normalized)
     allocate(                in_bin(num_bins))
 
     associate(v_min => (histogram%unmapped_min_), v_max => (histogram%unmapped_max_), cardinality => size(v))
-      associate(v_mapped => normalize(v, v_min, v_max), dv => (v_mapped_max - v_mapped_min)/real(num_bins))
-        associate(v_bin_min => [(v_mapped_min + (i-1)*dv, i=1,num_bins)])
-          associate(smidgen => .0001*abs(dv)) ! use to make the high end of the bin range inclusive of the max value
-            associate(v_bin_max => [v_bin_min(2:), v_mapped_max + smidgen])
-              do concurrent(i = 1:num_bins)
-                in_bin(i) = count(v_mapped >= v_bin_min(i) .and. v_mapped < v_bin_max(i)) ! replace with Fortran 2023 reduction
-                histogram%frequency_(i) = real(in_bin(i)) / real(cardinality)
-                histogram%bin_midpoint_(i) = v_bin_min(i) + 0.5*dv
-              end do
+      block
+        real, allocatable :: v_mapped(:,:,:,:)
+
+        if (raw)  then
+          v_mapped = v
+        else
+          v_mapped = normalize(v, v_min, v_max)
+        end if
+
+        associate(v_mapped_min => merge(v_min, 0., raw), v_mapped_max => merge(v_max, 1., raw))
+          associate(dv => (v_mapped_max - v_mapped_min)/real(num_bins))
+            associate(v_bin_min => [(v_mapped_min + (i-1)*dv, i=1,num_bins)])
+              associate(smidgen => .0001*abs(dv)) ! use to make the high end of the bin range inclusive of the max value
+                associate(v_bin_max => [v_bin_min(2:), v_mapped_max + smidgen])
+                  do concurrent(i = 1:num_bins)
+                    in_bin(i) = count(v_mapped >= v_bin_min(i) .and. v_mapped < v_bin_max(i)) ! replace with Fortran 2023 reduction
+                    histogram%frequency_(i) = real(in_bin(i)) / real(cardinality)
+                    histogram%bin_midpoint_(i) = v_bin_min(i) + 0.5*dv
+                  end do
+                end associate
+              end associate
             end associate
           end associate
         end associate
-      end associate
-      associate(binned => sum(in_bin))
-        call assert(cardinality == binned, "histogram_m(normalize): lossless binning", intrinsic_array_t([cardinality, binned]))
-      end associate
+        associate(binned => sum(in_bin))
+          call assert(cardinality == binned, "histogram_m(normalize): lossless binning", intrinsic_array_t([cardinality, binned]))
+        end associate
+      end block
     end associate
 
   end procedure