2
! Module : DiscreteSampler
3
! Author : Valentin Hirschi
6
! A relatively simple and flexible module to do
7
! sampling of discrete dimensions for Monte-Carlo
10
! List of public subroutines and usage :
12
! DS_initialize(tolerate_zero_grid, verbose)
13
! :: Initializes the general options of the DiscreteSampler.
14
! :: tolerate_zero_grid is a logical that specifies whether
15
! :: DiscreteSampler should crash when asked to sample a dimension
17
! :: Verbose is another logical setting the verbosity of the
20
! DS_register_dimension(name, n_bins,(all_grid|void))
21
! :: Register a new dimension with its name and number of bins
22
! :: If all_grid is specified and set to False, then only a
23
! :: running grid is specified, it is useful for registering
24
! :: grid intended for convolution only.
26
! DS_remove_dimension(name)
27
! :: Removes and clear the dimension of input name
29
! DS_print_global_info(name|index|void)
30
! :: Print global information on a registered information, using
31
! :: either its name or index. or all if none is selected
34
! :: Reinitializes the module and all grid data
36
! DS_binID(integer_ID)
37
! :: Creates an object of binID type from an integer ID. Notice
38
! :: that you can also use the assignment operator directly.
40
! DS_add_bin(dimension_name, (binID|integerID|void))
41
! :: Add one bin to dimension of name dimension_name.
42
! :: The user can add an ID to the added bin or just use the
43
! :: default sequential labeling
45
! DS_remove_bin(dimension_name, (binID|integerID))
46
! :: Remove a bin in a dimension by either specfiying its index
47
! :: in this dimension or its binID (remember you can create a
48
! :: binID on the flight from an integer using the function
51
! DS_add_entry(dimension_name, (binID|integerID), weight, (reset|void))
52
! :: Add a new weight to a certan bin (characterized by either
53
! :: its binID or the integer of the binID). If the dimension or
54
! :: the binID do not exist at the time of adding them, they will
55
! :: be registered on the flight. By setting the option 'reset'
56
! :: to '.True.' (default is '.False.'), you can chose to hardset
57
! :: the weight of this bin (with n_entries=min_bin_probing_points
58
! :: then) instead of adding it.
60
! DS_update_grid((dim_name|void))
61
! :: Update the reference grid of the dimension dim_name or
62
! :: update all of them at the same time without argument.
63
! :: It uses the running grid for this and it reinitilizes it.
65
! DS_write_grid((file_name|stream_id), (dim_name|void), (grid_type|void))
66
! :: Append to file 'file_name' or a given stream the data for
67
! :: the current reference grid for dimension dim_name or
68
! :: or all of them if called without dim_name.
69
! :: You can specify grid_type to be written out to be either 'ref'
70
! :: or 'run', and it is 'ref' by default.
72
! DS_load_grid((file_name|stream_id), (dim_name|void))
73
! :: Reset the run_grid for dimension dim_name and loads in
74
! :: the data obtained from file 'file_name' or stream_id for this
75
! :: dimension or all of them if called without dim_name.
76
! :: The data is loaded in the running grid only (which is
77
! :: re-initialized before this), so you have to call
78
! :: 'DS_update_grid' if you want it pased to the ref_grid.
80
! DS_get_point(dim_name, random_variable,
81
! (binIDPicked|integerIDPicked), jacobian_weight, (mode|void),
82
! (convoluted_grid_names|void))
83
! :: From a given random variable in [0.0,1.0] and a dimension
84
! :: name, this subroutine returns the picked bin or index and
85
! :: the Jacobian weight associated.
86
! :: The mode is by default 'norm' but can also be 'variance'
87
! :: which means that the sampling is done either on the
88
! :: variance in each bin or the absolute value of the weights.
90
! :: Jacobian = 1.0d0 / normalized_abs_wgt_in_selected_bin
91
! :: mode == 'variance'
92
! :: Jacobian = 1.0d0 / normalized_variance_in_selected_bin
93
! :: mode == 'uniform'
94
! :: Jacobian = 1.0d0 / n_bins_in_dimension
95
! :: Setting the option 'convoluted_grid_names' to an array of
96
! :: dimension names already registered in the module will have
97
! :: the subroutine DS_get_point return a point sampled as the
98
! :: grid 'dim_name' *convoluted with all grids specified in
99
! :: convoluted_grid_names*.
101
! :: call DS_get_point('MyDim',0.02,out_binPicked,out_jac,mode='norm',
102
! :: & convoluted_grid_names = (/'ConvolutionDim1','ConvolutionDim2'/))
104
! DS_set_min_points(min_point, (dim_name|void))
105
! :: Sets the minimum number of points that must be used to probe
106
! :: each bin of a particular dimension (or all if not specified)
107
! :: before DS_get_point uses a uniform sampling on the bins
108
! :: (possibly with convolution) when the reference grid is empty.
109
! :: By default it is 10.
111
! DS_get_dim_status(dim_name)
112
! :: Return an integer that specifies the status of a given dimension
113
! :: dim_status = -1 : Dimension is not registered.
114
! :: dim_status = 0 : Dimension exists but the reference grid is
115
! :: empty and the running grid does not yet
116
! :: have all its bins filled with min_point.
117
! :: dim_status = 1 : Dimension exists with a non-empty
118
! :: reference grid and a running grid with all
119
! :: bins filled with more than min_points entries.
121
! DS_set_grid_mode(grid_name,grid_mode)
122
! :: Sets the kind of reference grid which is currently active.
123
! :: grid_mode = 'default' : This means that the reference grid holds
124
! :: the same kind of weights than the running grid. When the reference
125
! :: grid will be updated, the running grid will be *combined* with
126
! :: the reference grid, and not overwritten by it.
127
! :: grid_mode = 'init' : This means that the reference grid is used for
128
! :: initialisation, and its weights do not compare with those put
129
! :: in the running grid. When updated, the reference grid will
130
! :: therefore be *overwritten* by the running grid.
133
! DS_get_damping_for_grid(grid_name, small_contrib, damping_power)
134
! :: Returns the current value stored in the run grid for
135
! :: dimension grid_name of what are the parameter damping the
136
! :: bin with small contributions and whose Jacobian can
137
! :: potentially be very large. See the definition of these
138
! :: parameters for a description of the procedure.
140
! DS_set_damping_for_grid(grid_name, small_contrib, damping_power)
141
! :: Sets the value for both the ref and running grid of the
142
! :: dimension grid_name of what are the parameter damping the
143
! :: bin with small contributions and whose Jacobian can
144
! :: potentially be very large. See the definition of these
145
! :: parameters for a description of the procedure.
147
module DiscreteSampler
151
! Global options for the module
154
logical DS_tolerate_zero_norm
155
save DS_tolerate_zero_norm
156
! An allocatable to mark initialization
157
logical, dimension(:), allocatable :: DS_isInitialized
159
! This parameter sets how large must be the sampling bar when
160
! displaying information about a dimension
161
integer samplingBarWidth
162
parameter (samplingBarWidth=80)
164
! Attributes identifying a bin
165
! For now just an integer
169
! And an easy way to create a binIDs
170
interface assignment (=)
171
module procedure binID_from_binID
172
module procedure binID_from_integer
173
end interface assignment (=)
174
! Define and easy way of comparing binIDs
175
interface operator (==)
176
module procedure equal_binID
177
end interface operator (==)
179
! Information relevant to a bin
182
! Sum of the squared weights, to compute the variance
184
! Sum of the absolute value of the weights put in this bin,
185
! necessary when in presence of negative weights.
188
! Practical to be able to identify a bin by its id
192
! Define and easy way of adding Bins
193
interface operator (+)
194
module procedure DS_combine_two_bins
195
end interface operator (+)
197
type sampledDimension
198
! The grid_mode can take the following values
199
! grid_mode = 1 : This means that the reference grid holds the
200
! same kind of weights than the running grid. When the reference
201
! grid will be updated, the running grid will be *combined* with
202
! the reference grid, and not overwrite it.
203
! grid_mode = 2 : This means that the reference grid is used for
204
! initialisation, and its weights do not compare with those put
205
! in the running grid. When updated, the reference grid will
206
! therefore be *overwritten* by the running grid.
209
! Treat specially bin with a contribution (i.e. weight) worth less than
210
! 'small_contrib_threshold' of the averaged contributionover all bins.
211
! For those, we sample according to the square root (or the specified power
212
! 'damping power' of the difference between the reference value corresponding
213
! to the chosen mode and the small_contrib_threshold.
214
! In this way, we are less sensitive to possible large fluctuations
215
! of very suppressed contributions for which the Jacobian would be
216
! really big. However, the square-root is such that a really
217
! suppressed contribution at the level of numerical precision
218
! would still never be probed.
219
! Notice that this procedure does *not* change the weight in the
220
! bin, but only how it is used for bin picking.
221
real*8 :: small_contrib_threshold
222
real*8 :: damping_power
223
! Minimum number of points to probe each bin with when the reference
224
! grid is empty. Once each bin has been probed that many times, the
225
! subroutine DS_get_point will use a uniform distribution
226
integer :: min_bin_probing_points
227
! Keep track of the norm (i.e. sum of all weights) and the total
228
! number of points for ease and optimisation purpose
230
! The sum of the absolute value of the weight in each bin
232
! The sum of the variance of the weight in each bin
233
real*8 :: variance_norm
234
! The sum of the squared weights in each bin
236
integer :: n_tot_entries
237
! A handy way of referring to the dimension by its name rather than
239
character, dimension(:), allocatable :: dimension_name
241
type(bin) , dimension(:), allocatable :: bins
242
endtype sampledDimension
244
! This stores the overall discrete reference grid
245
type(sampledDimension), dimension(:), allocatable :: ref_grid
247
! This is the running grid, whose weights are being updated for each point
248
! but not yet used for the sampling. The user must call the 'update'
249
! function for the running grid to be merged to the reference one.
250
type(sampledDimension), dimension(:), allocatable :: run_grid
252
interface DS_add_entry
253
module procedure DS_add_entry_with_BinID
254
module procedure DS_add_entry_with_BinIntID
255
end interface DS_add_entry ! DS_add_entry
257
interface DS_print_global_info
258
module procedure DS_print_dim_global_info_from_name
259
module procedure DS_print_dim_global_info_from_void
260
end interface ! DS_print_dim_global_info
263
module procedure DS_add_bin_with_binID
264
module procedure DS_add_bin_with_IntegerID
265
module procedure DS_add_bin_with_void
266
end interface ! DS_add_bin
268
interface DS_remove_bin
269
module procedure DS_remove_bin_withIntegerID
270
module procedure DS_remove_bin_withBinID
271
end interface ! DS_remove_bin
274
module procedure DS_get_bin_from_binID
275
module procedure DS_get_bin_from_binID_and_dimName
276
end interface ! DS_get_bin
278
interface DS_update_grid
279
module procedure DS_update_grid_with_dim_name
280
module procedure DS_update_all_grids
281
module procedure DS_update_grid_with_dim_index
282
end interface ! DS_update_grid
284
interface DS_write_grid
285
module procedure DS_write_grid_with_filename
286
module procedure DS_write_grid_with_streamID
287
end interface ! DS_write_grid
289
interface DS_load_grid
290
module procedure DS_load_grid_with_filename
291
module procedure DS_load_grid_with_streamID
292
end interface ! DS_load_grid
294
interface DS_dim_index
295
module procedure DS_dim_index_default
296
module procedure DS_dim_index_with_force
297
module procedure DS_dim_index_default_with_chararray
298
module procedure DS_dim_index_with_force_with_chararray
299
end interface ! DS_dim_index
301
interface DS_bin_index
302
module procedure DS_bin_index_default
303
module procedure DS_bin_index_with_force
304
end interface ! DS_bin_index
306
interface DS_get_point
307
module procedure DS_get_point_with_integerBinID
308
module procedure DS_get_point_with_BinID
309
end interface ! DS_get_point
313
! ---------------------------------------------------------------
314
! This subroutine is simply the logger of this module
315
! ---------------------------------------------------------------
316
subroutine DS_initialize(tolerate_zero_norm, verbose)
319
! Subroutine arguments
321
logical, optional, intent(in) :: tolerate_zero_norm
322
logical, optional, intent(in) :: verbose
324
if (allocated(DS_isInitialized)) then
325
write(*,*) "DiscreteSampler:: Error: The DiscreteSampler"//
326
& " module can only be initialized once."
329
allocate(DS_isInitialized(1))
331
if (present(verbose)) then
337
if (present(tolerate_zero_norm)) then
338
DS_tolerate_zero_norm = tolerate_zero_norm
340
DS_tolerate_zero_norm = .True.
343
! Re-instore the if statement below if too annoying
344
! if(DS_verbose) then
346
write(*,*) '********************************************'
347
write(*,*) '* You are using the DiscreteSampler module *'
348
write(*,*) '* part of the MG5_aMC framework *'
349
write(*,*) '* Author: Valentin Hirschi *'
350
write(*,*) '********************************************'
354
end subroutine DS_initialize
356
! ---------------------------------------------------------------
357
! This subroutine is simply the logger of this module
358
! ---------------------------------------------------------------
360
subroutine DS_Logger(msg)
363
! Subroutine arguments
365
character(len=*), intent(in) :: msg
367
if (DS_verbose) write(*,*) msg
369
end subroutine DS_Logger
371
! ---------------------------------------------------------------
372
! This subroutine clears the module and reinitialize all data
373
! ---------------------------------------------------------------
374
subroutine DS_clear()
375
call DS_deallocate_grid(ref_grid)
376
call DS_deallocate_grid(run_grid)
377
end subroutine DS_clear
379
subroutine DS_deallocate_grid(grid)
381
type(sampledDimension), dimension(:), allocatable,
382
& intent(inout) :: grid
383
if (allocated(grid)) then
385
if (allocated(grid(i)%bins)) then
386
deallocate(grid(i)%bins)
388
if (allocated(grid(i)%dimension_name)) then
389
deallocate(grid(i)%dimension_name)
394
end subroutine DS_deallocate_grid
396
! ---------------------------------------------------------------
397
! This subroutine takes care of registering a new dimension in
398
! the DSampler module by characterizin it by its name and number
400
! ---------------------------------------------------------------
401
subroutine DS_register_dimension(dim_name,n_bins,all_grids)
404
! Subroutine arguments
406
integer , intent(in) :: n_bins
407
character(len=*), intent(in) :: dim_name
408
logical , optional :: all_grids
412
logical :: do_all_grids
416
! Make sure the module is initialized
417
if (.not.allocated(DS_isInitialized)) then
420
if (present(all_grids)) then
421
do_all_grids = all_grids
423
do_all_grids = .True.
425
if (do_all_grids) then
426
call DS_add_dimension_to_grid(ref_grid, dim_name, n_bins)
428
call DS_add_dimension_to_grid(run_grid, dim_name, n_bins)
430
call DS_Logger("DiscreteSampler:: Successfully registered "//
431
$ "dimension '"//dim_name//"' ("//TRIM(toStr(n_bins))//' bins)')
433
end subroutine DS_register_dimension
435
! ---------------------------------------------------------------
436
! This subroutine registers a dimension to a given grid
437
! ---------------------------------------------------------------
438
subroutine DS_add_dimension_to_grid(grid, dim_name, n_bins)
441
! Subroutine arguments
443
type(sampledDimension), dimension(:), allocatable,
444
& intent(inout) :: grid
445
integer , intent(in) :: n_bins
446
character(len=*), intent(in) :: dim_name
451
type(sampledDimension), dimension(:), allocatable :: tmp
456
! Make sure the module is initialized
457
if (.not.allocated(DS_isInitialized)) then
460
if(allocated(grid)) then
461
dim_index = DS_dim_index(grid,dim_name,.True.)
462
if (dim_index.ne.-1) then
463
write(*,*) 'DiscreteSampler:: Error, the dimension'//
464
$ " with name '"//dim_name//"' is already registered."
469
! Either allocate the discrete grids or append a dimension
470
if (allocated(grid)) then
471
allocate(tmp(size(grid)))
473
call DS_copy_dimension(grid(i), tmp(i))
475
call DS_deallocate_grid(grid)
476
allocate(grid(size(tmp)+1))
478
call DS_copy_dimension(tmp(i), grid(i))
480
call DS_deallocate_grid(tmp)
484
! Now we can fill in the appended element with the
485
! characteristics of the dimension which must be added
486
allocate(grid(size(grid))%bins(n_bins))
487
allocate(grid(size(grid))%dimension_name(len(dim_name)))
488
! Initialize the values of the grid with default
489
call DS_initialize_dimension(grid(size(grid)))
490
! For the string assignation, I have to it character by
492
do i=1, len(dim_name)
493
grid(size(grid))%dimension_name(i) = dim_name(i:i)
496
end subroutine DS_add_dimension_to_grid
498
! ----------------------------------------------------------------------
499
! Copy a dimension from source to target, making sure to allocate
500
! ----------------------------------------------------------------------
501
subroutine DS_copy_dimension(source, trget)
502
type(sampledDimension), intent(out) :: trget
503
type(sampledDimension), intent(in) :: source
506
if (allocated(trget%bins)) then
507
deallocate(trget%bins)
509
allocate(trget%bins(size(source%bins)))
510
do i=1,size(source%bins)
511
call DS_copy_bin(source%bins(i),trget%bins(i))
513
if (allocated(trget%dimension_name)) then
514
deallocate(trget%dimension_name)
516
allocate(trget%dimension_name(size(source%dimension_name)))
517
do i=1, size(source%dimension_name)
518
trget%dimension_name(i) = source%dimension_name(i)
520
trget%norm = source%norm
521
trget%abs_norm = source%abs_norm
522
trget%variance_norm = source%variance_norm
523
trget%norm_sqr = source%norm_sqr
524
trget%n_tot_entries = source%n_tot_entries
525
trget%min_bin_probing_points = source%min_bin_probing_points
526
trget%grid_mode = source%grid_mode
527
trget%damping_power = source%damping_power
528
trget%small_contrib_threshold = source%small_contrib_threshold
529
end subroutine DS_copy_dimension
531
! ----------------------------------------------------------------------
532
! This subroutine removes a dimension at index dim_index from a given grid
533
! ----------------------------------------------------------------------
534
subroutine DS_remove_dimension(dim_name)
537
! Subroutine arguments
539
character(len=*), intent(in) :: dim_name
543
integer :: ref_dim_index, run_dim_index
547
ref_dim_index = DS_dim_index(ref_grid, dim_name)
548
run_dim_index = DS_dim_index(run_grid, dim_name)
549
call DS_remove_dimension_from_grid(ref_grid, ref_dim_index)
550
call DS_remove_dimension_from_grid(run_grid, run_dim_index)
551
end subroutine DS_remove_dimension
554
! ----------------------------------------------------------------------
555
! This subroutine removes a dimension at index dim_index from a given grid
556
! ----------------------------------------------------------------------
557
subroutine DS_remove_dimension_from_grid(grid, dim_index)
560
! Subroutine arguments
562
type(sampledDimension), dimension(:), allocatable,
563
& intent(inout) :: grid
564
integer, intent(in) :: dim_index
568
type(sampledDimension), dimension(:), allocatable :: tmp
574
allocate(tmp(size(grid)-1))
576
call DS_copy_dimension(grid(i), tmp(i))
578
do i=dim_index+1,size(grid)
579
call DS_copy_dimension(grid(i), tmp(i-1))
581
call DS_deallocate_grid(grid)
582
allocate(grid(size(tmp)))
584
call DS_copy_dimension(tmp(i), grid(i))
586
call DS_deallocate_grid(tmp)
587
end subroutine DS_remove_dimension_from_grid
589
! ---------------------------------------------------------------
590
! This subroutine takes care of reinitializing a given dimension
591
! with default values
592
! ---------------------------------------------------------------
593
subroutine DS_reinitialize_dimension(d_dim)
596
! Subroutine arguments
598
type(sampledDimension), intent(inout) :: d_dim
607
do i=1, size(d_dim%bins)
608
call DS_reinitialize_bin(d_dim%bins(i))
610
d_dim%norm_sqr = 0.0d0
611
d_dim%abs_norm = 0.0d0
612
d_dim%variance_norm = 0.0d0
614
d_dim%n_tot_entries = 0
616
end subroutine DS_reinitialize_dimension
618
! ---------------------------------------------------------------
619
! This subroutine takes care of initializing a given dimension
620
! with default values
621
! ---------------------------------------------------------------
622
subroutine DS_initialize_dimension(d_dim)
625
! Subroutine arguments
627
type(sampledDimension), intent(inout) :: d_dim
636
call DS_reinitialize_dimension(d_dim)
637
do i=1, size(d_dim%bins)
638
call DS_initialize_bin(d_dim%bins(i))
640
do i= 1, len(d_dim%dimension_name)
641
d_dim%dimension_name(i:i) = ' '
643
! By default require each bin to be probed by 10 points
644
! before a uniform distribution is used when the reference grid
646
d_dim%min_bin_probing_points = 10
648
! Turn off the damping of small contributions by default
649
d_dim%small_contrib_threshold = 0.0d0
650
d_dim%damping_power = 0.5d0
651
! By default give sequential ids to the bins
652
do i=1, size(d_dim%bins)
653
d_dim%bins(i)%bid = i
655
end subroutine DS_initialize_dimension
657
! ---------------------------------------------------------------
658
! This subroutine takes care of reinitializing a given bin
659
! ---------------------------------------------------------------
660
subroutine DS_initialize_bin(d_bin)
663
! Subroutine arguments
665
type(bin), intent(inout) :: d_bin
669
call DS_reinitialize_bin(d_bin)
671
end subroutine DS_initialize_bin
673
! ---------------------------------------------------------------
674
! This subroutine takes care of initializing a given bin
675
! ---------------------------------------------------------------
676
subroutine DS_reinitialize_bin(d_bin)
679
! Subroutine arguments
681
type(bin), intent(inout) :: d_bin
685
d_bin%weight_sqr = 0.0d0
686
d_bin%abs_weight = 0.0d0
689
end subroutine DS_reinitialize_bin
691
! ---------------------------------------------------------------
692
! Set the minimum number of point for which the bins must be
693
! probed before a uniform distribution is used when the reference
695
! ---------------------------------------------------------------
696
subroutine DS_set_min_points(min_points, dim_name)
699
! Subroutine arguments
702
integer, intent(in) :: min_points
703
character(len=*), intent(in), optional :: dim_name
711
if(present(dim_name)) then
712
ref_grid(DS_dim_index(ref_grid,dim_name))%
713
& min_bin_probing_points = min_points
714
run_grid(DS_dim_index(ref_grid,dim_name))%
715
& min_bin_probing_points = min_points
717
do i=1,size(ref_grid)
718
ref_grid(i)%min_bin_probing_points = min_points
719
run_grid(i)%min_bin_probing_points = min_points
722
end subroutine DS_set_min_points
724
! ---------------------------------------------------------------
725
! Returns an integer that specifies the status of a given dimension
726
! dim_status = -1 : Dimension is not registered.
727
! dim_status = 0 : Dimension exists but the reference grid is
728
! empty and the running grid does not yet
729
! have all its bins filled with min_point.
730
! dim_status = 1 : Dimension exists with a non-empty
731
! reference grid and a running grid with all
732
! bins filled with more than min_points entries.
733
! ---------------------------------------------------------------
734
function DS_get_dim_status(grid_name)
739
character(len=*), intent(in) :: grid_name
740
integer :: DS_get_dim_status
744
integer :: ref_grid_index
745
integer :: run_grid_index
746
integer :: int_grid_mode
752
ref_grid_index = DS_dim_index(ref_grid, grid_name, .True.)
753
run_grid_index = DS_dim_index(run_grid, grid_name, .True.)
754
if (ref_grid_index.eq.-1.or.run_grid_index.eq.-1) then
755
DS_get_dim_status = -1
758
if (ref_grid(ref_grid_index)%n_tot_entries.ne.0) then
759
DS_get_dim_status = 1
763
! If the running grid has zero entries, then consider the grid
765
if(size(run_grid(run_grid_index)%bins).eq.0) then
766
DS_get_dim_status = 0
770
do i=1,size(ref_grid(ref_grid_index)%bins)
771
mRunBin = DS_get_bin(run_grid(run_grid_index)%bins,
772
& ref_grid(ref_grid_index)%bins(i)%bid)
773
if (mRunBin%n_entries.lt.ref_grid(ref_grid_index)
774
& %min_bin_probing_points) then
775
DS_get_dim_status = 0
780
DS_get_dim_status = 1
782
end function DS_get_dim_status
784
! ---------------------------------------------------------------
785
! Access function to modify the damping parameters of small
787
! ---------------------------------------------------------------
788
subroutine DS_set_damping_for_grid(grid_name, in_small_contrib,
792
! Subroutine arguments
794
character(len=*), intent(in) :: grid_name
795
real*8, intent(in) :: in_small_contrib
796
real*8, intent(in) :: in_damping_power
800
integer :: ref_grid_index
801
integer :: run_grid_index
805
ref_grid_index = DS_dim_index(ref_grid, grid_name, .True.)
806
if (ref_grid_index.eq.-1) then
807
write(*,*) "DiscreteSampler:: Error in 'DS_set_damping_"//
808
& "for_grid', dimension '"//grid_name//"' could not be"//
809
& " found in the reference grid."
812
run_grid_index = DS_dim_index(run_grid, grid_name, .True.)
813
if (run_grid_index.eq.-1) then
814
write(*,*) "DiscreteSampler:: Error in 'DS_set_damping_"//
815
& "for_grid', dimension '"//grid_name//"' could not be"//
816
& " found in the running grid."
820
! Limit arbitrarily at 50% because anything above that really
821
! breaks the assumption of a small grid deformation not
822
! significantly affecting the averaged contribution taked as
824
if (in_small_contrib.lt.0.0d0.or.
825
& in_small_contrib.gt.0.5d0) then
826
write(*,*) "The small relative contribution threshold "//
827
& toStr_real_with_ndig(in_small_contrib,3)
828
& //") given in argument of the function 'DS_set_damping_"//
829
& "for_grid' must be >=0.0 and <= 0.5."
833
if (in_damping_power.lt.0.0d0.or.
834
& in_damping_power.gt.1.0d0) then
835
write(*,*) "The damping power ("//
836
& toStr_real_with_ndig(in_damping_power,3)
837
& //") given in argument of the function 'DS_set_damping_"//
838
& "for_grid' must be >= 0.0 and <= 1.0."
842
ref_grid(ref_grid_index)%small_contrib_threshold =
844
ref_grid(ref_grid_index)%damping_power = in_damping_power
845
run_grid(run_grid_index)%small_contrib_threshold =
847
run_grid(run_grid_index)%damping_power = in_damping_power
848
end subroutine DS_set_damping_for_grid
850
! ---------------------------------------------------------------
851
! Access function to access the damping parameters for small
852
! contributions stored in the reference grid
853
! ---------------------------------------------------------------
854
subroutine DS_get_damping_for_grid(grid_name, out_small_contrib,
858
! Subroutine arguments
860
character(len=*), intent(in) :: grid_name
861
real*8, intent(out) :: out_small_contrib
862
real*8, intent(out) :: out_damping_power
866
integer :: run_grid_index
870
run_grid_index = DS_dim_index(run_grid, grid_name, .True.)
871
if (run_grid_index.eq.-1) then
872
write(*,*) "DiscreteSampler:: Error in 'DS_get_damping_"//
873
& "for_grid', dimension '"//grid_name//"' could not be"//
874
& " found in the running grid."
878
out_small_contrib = run_grid(run_grid_index)%
879
& small_contrib_threshold
880
out_damping_power = run_grid(run_grid_index)%damping_power
882
end subroutine DS_get_damping_for_grid
884
! ---------------------------------------------------------------
885
! Access function to modify the mode of the reference grid:
886
! grid_mode = 'default' : This means that the reference grid holds
887
! the same kind of weights than the running grid. When the reference
888
! grid will be updated, the running grid will be *combined* with
889
! the reference grid, and not overwritten by it.
890
! grid_mode = 'init' : This means that the reference grid is used for
891
! initialisation, and its weights do not compare with those put
892
! in the running grid. When updated, the reference grid will
893
! therefore be *overwritten* by the running grid.
894
! ---------------------------------------------------------------
895
subroutine DS_set_grid_mode(grid_name, grid_mode)
898
! Subroutine arguments
900
character(len=*), intent(in) :: grid_mode
901
character(len=*), intent(in) :: grid_name
905
integer :: ref_grid_index
906
integer :: int_grid_mode
910
ref_grid_index = DS_dim_index(ref_grid, grid_name, .True.)
911
if (ref_grid_index.eq.-1) then
912
write(*,*) 'DiscreteSampler:: Error in DS_set_grid_mode, '//
913
& "dimension '"//grid_name//"' could not be found in the "//
917
if (grid_mode.eq.'init') then
919
elseif (grid_mode.eq.'default') then
922
write(*,*) 'DiscreteSampler:: Error in DS_set_grid_mode, '//
923
& " grid_mode '"//grid_mode//"' not recognized. It must "//
924
& " be one of the following: 'default', 'init'."
928
! Notice that we don't change the mode of the running_grid
929
! because in this way, after any DS_update() is done, the
930
! ref_grid will automatically turn its mode to 'default' because
931
! it inherits the attribute of the running grid.
932
! However, if the running grid was loaded from a saved grid file
933
! then it might be that the run_grid also has the grid_mode set
934
! to 'initialization' which will then correctly be copied to the
935
! ref_grid after the DS_update() performed at the end of
936
! DS_load() which correctly reproduce the state of the
937
! DiscreteSampler module at the time it wrote the grids.
938
ref_grid(ref_grid_index)%grid_mode = int_grid_mode
939
end subroutine DS_set_grid_mode
941
! ---------------------------------------------------------------
942
! Dictionary access-like subroutine to obtain a grid from its name
943
! ---------------------------------------------------------------
945
function DS_get_dimension(grid, dim_name)
950
type(sampledDimension), dimension(:), intent(in), allocatable
952
character(len=*), intent(in) :: dim_name
953
type(sampledDimension) :: DS_get_dimension
957
DS_get_dimension = grid(DS_dim_index(grid,dim_name))
958
end function DS_get_dimension
960
! ---------------------------------------------------------------
961
! Returns the index of a bin with mBinID in the list bins
962
! ---------------------------------------------------------------
963
function DS_bin_index_default(bins, mBinID)
968
type(Bin), dimension(:), intent(in)
970
type(BinID) :: mBinID
971
integer :: DS_bin_index_default
975
DS_bin_index_default = DS_bin_index_with_force(bins,mBinID,
977
end function DS_bin_index_default
979
function DS_bin_index_with_force(bins, mBinID,force)
984
type(Bin), dimension(:), intent(in)
986
type(BinID) :: mBinID
987
integer :: DS_bin_index_with_force
996
! For efficiency first look at index mBinID%id
997
if (size(bins).ge.mBinID%id) then
998
if (bins(mBinID%id)%bid==mBinID) then
999
DS_bin_index_with_force = mBinID%id
1004
DS_bin_index_with_force = -1
1005
do i = 1, size(bins)
1006
if (bins(i)%bid==mBinID) then
1007
DS_bin_index_with_force = i
1011
if (DS_bin_index_with_force.eq.-1.and.(.not.Force)) then
1012
write(*,*) 'DiscreteSampler:: Error in function bin_index'//
1013
& "(), bin with BinID '"//trim(DS_toStr(mBinID))
1017
end function DS_bin_index_with_force
1019
! ---------------------------------------------------------------
1020
! Functions of the interface get_bin facilitating the access to a
1022
! ---------------------------------------------------------------
1024
function DS_get_bin_from_binID(bins, mBinID)
1027
! Function arguments
1029
type(Bin), dimension(:), intent(in)
1031
type(BinID) :: mBinID
1032
type(Bin) :: DS_get_bin_from_binID
1040
DS_get_bin_from_binID = bins(DS_bin_index(bins,mBinID))
1041
end function DS_get_bin_from_binID
1043
function DS_get_bin_from_binID_and_dimName(grid, dim_name,
1047
! Function arguments
1049
type(sampledDimension), dimension(:), intent(in), allocatable
1051
character(len=*), intent(in) :: dim_name
1052
type(BinID) :: mBinID
1053
type(Bin) :: DS_get_bin_from_binID_and_dimName
1058
type(SampledDimension) :: m_dim
1062
m_dim = DS_get_dimension(grid,dim_name)
1063
DS_get_bin_from_binID_and_dimName = DS_get_bin_from_binID(
1064
& m_dim%bins,mBinID)
1065
end function DS_get_bin_from_binID_and_dimName
1068
! ---------------------------------------------------------------
1069
! Add a new weight to a certan bin (characterized by either its
1071
! ---------------------------------------------------------------
1072
subroutine DS_add_entry_with_BinID(dim_name, mBinID, weight,
1076
! Subroutine arguments
1078
character(len=*), intent(in) :: dim_name
1079
type(BinID) :: mBinID
1081
logical, optional :: reset
1085
integer dim_index, bin_index
1087
integer :: n_entries
1088
logical :: opt_reset
1092
if (present(reset)) then
1098
dim_index = DS_dim_index(run_grid, dim_name, .TRUE.)
1099
if (dim_index.eq.-1) then
1100
call DS_Logger('Dimension '//dim_name//
1101
& ' does not exist in the run grid. Creating it now.')
1102
call DS_register_dimension(dim_name,0)
1103
dim_index = DS_dim_index(run_grid, dim_name)
1106
bin_index = DS_bin_index(
1107
& run_grid(dim_index)%bins,mBinID,.TRUE.)
1108
if (bin_index.eq.-1) then
1109
call DS_Logger('Bin with binID '//trim(DS_toStr(mBinID))//
1110
& ' does not exist in the run grid. Creating it now.')
1111
call DS_reinitialize_bin(newBin)
1113
call DS_add_bin_to_bins(run_grid(dim_index)%bins,newBin)
1114
bin_index = DS_bin_index(run_grid(dim_index)%bins,mBinID)
1117
! First remove bin from global cumulative information in the grid
1118
run_grid(dim_index)%norm = run_grid(dim_index)%norm -
1119
& run_grid(dim_index)%bins(bin_index)%weight
1120
run_grid(dim_index)%norm_sqr = run_grid(dim_index)%norm_sqr -
1121
& run_grid(dim_index)%bins(bin_index)%weight_sqr
1122
run_grid(dim_index)%abs_norm = run_grid(dim_index)%abs_norm -
1123
& run_grid(dim_index)%bins(bin_index)%abs_weight
1124
run_grid(dim_index)%variance_norm =
1125
& run_grid(dim_index)%variance_norm -
1126
& DS_bin_variance(run_grid(dim_index)%bins(bin_index))
1127
run_grid(dim_index)%n_tot_entries =
1128
& run_grid(dim_index)%n_tot_entries -
1129
& run_grid(dim_index)%bins(bin_index)%n_entries
1130
! Update the information directly stored in the bin
1131
if(.not.opt_reset) then
1132
n_entries = run_grid(dim_index)%bins(bin_index)%n_entries
1133
run_grid(dim_index)%bins(bin_index)%weight =
1134
& (run_grid(dim_index)%bins(bin_index)%weight*n_entries
1135
& + weight)/(n_entries+1)
1136
run_grid(dim_index)%bins(bin_index)%weight_sqr =
1137
& (run_grid(dim_index)%bins(bin_index)%weight_sqr*n_entries
1138
& + weight**2)/(n_entries+1)
1139
run_grid(dim_index)%bins(bin_index)%abs_weight =
1140
& (run_grid(dim_index)%bins(bin_index)%abs_weight*n_entries
1141
& + abs(weight))/(n_entries+1)
1142
run_grid(dim_index)%bins(bin_index)%n_entries = n_entries+1
1144
run_grid(dim_index)%bins(bin_index)%weight = weight
1145
run_grid(dim_index)%bins(bin_index)%weight_sqr = weight**2
1146
run_grid(dim_index)%bins(bin_index)%abs_weight = abs(weight)
1147
run_grid(dim_index)%bins(bin_index)%n_entries =
1148
& run_grid(dim_index)%min_bin_probing_points
1150
! Now add the bin information back to the info in the grid
1151
run_grid(dim_index)%norm = run_grid(dim_index)%norm +
1152
& run_grid(dim_index)%bins(bin_index)%weight
1153
run_grid(dim_index)%norm_sqr = run_grid(dim_index)%norm_sqr +
1154
& run_grid(dim_index)%bins(bin_index)%weight_sqr
1155
run_grid(dim_index)%abs_norm = run_grid(dim_index)%abs_norm +
1156
& run_grid(dim_index)%bins(bin_index)%abs_weight
1157
run_grid(dim_index)%variance_norm =
1158
& run_grid(dim_index)%variance_norm +
1159
& DS_bin_variance(run_grid(dim_index)%bins(bin_index))
1160
run_grid(dim_index)%n_tot_entries =
1161
& run_grid(dim_index)%n_tot_entries +
1162
& run_grid(dim_index)%bins(bin_index)%n_entries
1164
end subroutine DS_add_entry_with_BinID
1166
subroutine DS_add_entry_with_BinIntID(dim_name, BinIntID,
1170
! Subroutine arguments
1172
character(len=*), intent(in) :: dim_name
1175
logical, optional :: reset
1179
if (present(reset)) then
1180
call DS_add_entry_with_BinID(dim_name, DS_BinID(BinIntID),
1183
call DS_add_entry_with_BinID(dim_name, DS_BinID(BinIntID),
1186
end subroutine DS_add_entry_with_BinIntID
1188
! ---------------------------------------------------------------
1189
! Prints out all informations for dimension of index d_index, or
1191
! ---------------------------------------------------------------
1192
subroutine DS_print_dim_global_info_from_void()
1194
if(allocated(ref_grid).and.allocated(run_grid)) then
1195
do i = 1, size(ref_grid)
1196
call DS_print_dim_global_info_from_name(
1197
& trim(toStr(ref_grid(i)%dimension_name)))
1200
write(*,*) 'DiscreteSampler:: No dimension setup yet.'
1202
end subroutine DS_print_dim_global_info_from_void
1204
subroutine DS_print_dim_global_info_from_name(d_name)
1207
! Function arguments
1209
character(len=*), intent(in) :: d_name
1213
integer n_bins, ref_dim_index, run_dim_index
1217
! The running grid and ref grid must have the same number of
1218
! bins at this stage
1220
if(.not.(allocated(ref_grid).and.allocated(run_grid))) then
1221
write(*,*) 'DiscreteSampler:: No dimension setup yet.'
1225
ref_dim_index = DS_dim_index(ref_grid,d_name,.TRUE.)
1226
run_dim_index = DS_dim_index(run_grid,d_name,.TRUE.)
1228
if (ref_dim_index.ne.-1) then
1229
n_bins = size(ref_grid(DS_dim_index(ref_grid,d_name))%bins)
1230
elseif (run_dim_index.ne.-1) then
1231
n_bins = size(run_grid(DS_dim_index(run_grid,d_name))%bins)
1233
write(*,*) 'DiscreteSampler:: No grid registered for name'//
1234
& " '"//d_name//"'."
1238
write(*,*) "DiscreteSampler:: ========================"//
1239
& "=========================="
1240
write(*,*) "DiscreteSampler:: Information for dimension '"//
1241
& d_name//"' ("//trim(toStr(n_bins))//" bins):"
1242
write(*,*) "DiscreteSampler:: -> Grids status ID : "//
1243
& trim(toStr(DS_get_dim_status(d_name)))
1244
if (ref_dim_index.ne.-1) then
1245
write(*,*) "DiscreteSampler:: || Reference grid "
1246
select case(ref_grid(ref_dim_index)%grid_mode)
1248
write(*,*) "DiscreteSampler:: -> Grid mode : default"
1250
write(*,*) "DiscreteSampler:: -> Grid mode : "//
1253
call DS_print_dim_info(ref_grid(ref_dim_index))
1255
write(*,*) "DiscreteSampler:: || No reference grid for "//
1258
if (run_dim_index.ne.-1) then
1259
write(*,*) "DiscreteSampler:: || Running grid "
1260
write(*,*) "DiscreteSampler:: -> Initialization "//
1261
& "minimum points : "//trim(toStr(run_grid(
1262
& run_dim_index)%min_bin_probing_points))
1263
call DS_print_dim_info(run_grid(run_dim_index))
1265
write(*,*) "DiscreteSampler:: || No running grid for "//
1268
write(*,*) "DiscreteSampler:: ========================"//
1269
& "=========================="
1270
end subroutine DS_print_dim_global_info_from_name
1272
! ---------------------------------------------------------------
1273
! Print all informations related to a specific sampled dimension
1275
! ---------------------------------------------------------------
1276
subroutine DS_print_dim_info(d_dim)
1279
! Function arguments
1281
type(sampledDimension), intent(in) :: d_dim
1285
integer i,j, curr_pos1, curr_pos2, curr_pos3
1286
integer n_bins, bin_width
1287
! Adding the minimum size for the separators '|' and binID assumed
1288
! of being of length 2 at most, so 10*2+11 and + 20 security :)
1290
character(samplingBarWidth+10*2+11+20) :: samplingBar1
1291
character(samplingBarWidth+10*2+11+20) :: samplingBar2
1292
character(samplingBarWidth+10*2+11+20) :: samplingBar3
1293
real*8 :: tot_entries, tot_variance, tot_abs_weight
1298
! Setup the sampling bars
1301
tot_variance = 0.0d0
1302
tot_abs_weight = 0.0d0
1303
do i=1,min(size(d_dim%bins),10)
1304
tot_entries = tot_entries + d_dim%bins(i)%n_entries
1305
tot_variance = tot_variance + DS_bin_variance(d_dim%bins(i))
1306
tot_abs_weight = tot_abs_weight + d_dim%bins(i)%abs_weight
1308
if (d_dim%n_tot_entries.eq.0) then
1309
samplingBar1 = "| Empty grid |"
1310
samplingBar2 = "| Empty grid |"
1311
samplingBar3 = "| Empty grid |"
1313
do i=1,len(samplingBar1)
1314
samplingBar1(i:i)=' '
1316
do i=1,len(samplingBar2)
1317
samplingBar2(i:i)=' '
1319
do i=1,len(samplingBar3)
1320
samplingBar3(i:i)=' '
1322
samplingBar1(1:1) = '|'
1323
samplingBar2(1:1) = '|'
1324
samplingBar3(1:1) = '|'
1328
do i=1,min(10,size(d_dim%bins))
1329
samplingBar1(curr_pos1:curr_pos1+1) =
1330
& trim(DS_toStr(d_dim%bins(i)%bid))
1331
samplingBar2(curr_pos2:curr_pos2+1) =
1332
& trim(DS_toStr(d_dim%bins(i)%bid))
1333
samplingBar3(curr_pos3:curr_pos3+1) =
1334
& trim(DS_toStr(d_dim%bins(i)%bid))
1335
curr_pos1 = curr_pos1+2
1336
curr_pos2 = curr_pos2+2
1337
curr_pos3 = curr_pos3+2
1339
if (tot_abs_weight.ne.0.0d0) then
1340
bin_width = int((d_dim%bins(i)%abs_weight/
1341
& tot_abs_weight)*samplingBarWidth)
1343
samplingBar1(curr_pos1+j:curr_pos1+j) = ' '
1345
curr_pos1 = curr_pos1+bin_width+1
1346
samplingBar1(curr_pos1:curr_pos1) = '|'
1347
curr_pos1 = curr_pos1+1
1350
if (tot_entries.ne.0) then
1351
bin_width = int((float(d_dim%bins(i)%n_entries)/
1352
& tot_entries)*samplingBarWidth)
1354
samplingBar2(curr_pos2+j:curr_pos2+j) = ' '
1356
curr_pos2 = curr_pos2+bin_width+1
1357
samplingBar2(curr_pos2:curr_pos2) = '|'
1358
curr_pos2 = curr_pos2+1
1361
if (tot_variance.ne.0.0d0) then
1362
bin_width = int((DS_bin_variance(d_dim%bins(i))/
1363
& tot_variance)*samplingBarWidth)
1365
samplingBar3(curr_pos3+j:curr_pos3+j) = ' '
1367
curr_pos3 = curr_pos3+bin_width+1
1368
samplingBar3(curr_pos3:curr_pos3) = '|'
1369
curr_pos3 = curr_pos3+1
1372
if (tot_abs_weight.eq.0.0d0) then
1373
samplingBar1 = "| All considered bins have zero weight |"
1375
if (tot_entries.eq.0) then
1376
samplingBar2 = "| All considered bins have no entries |"
1378
if (tot_variance.eq.0.0d0) then
1379
samplingBar3 = "| All variances are zeros in considered"//
1380
& " bins. Maybe not enough entries (need at least one bin"//
1381
& " with >=2 entries). |"
1387
n_bins = size(d_dim%bins)
1389
write(*,*) "DiscreteSampler:: -> Total number of "//
1390
& "entries : "//trim(toStr(d_dim%n_tot_entries))
1391
if (n_bins.gt.10) then
1392
write(*,*) "DiscreteSampler:: -> Sampled as"//
1393
& " (first 10 bins):"
1395
write(*,*) "DiscreteSampler:: -> Sampled as:"
1397
write(*,*) "DiscreteSampler:: "//trim(samplingBar2)
1398
write(*,*) "DiscreteSampler:: -> (norm_sqr , "//
1399
& "abs_norm , norm , variance ) :"
1400
write(*,*) "DiscreteSampler:: ("//
1401
& trim(toStr(d_dim%norm_sqr,'Ew.3'))//", "//
1402
& trim(toStr(d_dim%abs_norm,'Ew.3'))//", "//
1403
& trim(toStr(d_dim%norm,'Ew.3'))//", "//
1404
& trim(toStr(d_dim%variance_norm,'Ew.3'))//")"
1405
if (n_bins.gt.10) then
1406
write(*,*) "DiscreteSampler:: -> Abs weights sampled as"//
1407
& " (first 10 bins):"
1409
write(*,*) "DiscreteSampler:: -> Abs weights sampled as:"
1411
write(*,*) "DiscreteSampler:: "//trim(samplingBar1)
1412
if (n_bins.gt.10) then
1413
write(*,*) "DiscreteSampler:: -> Variance sampled as"//
1414
& " (first 10 bins):"
1416
write(*,*) "DiscreteSampler:: -> Variance sampled as:"
1418
write(*,*) "DiscreteSampler:: "//trim(samplingBar3)
1420
end subroutine DS_print_dim_info
1422
! ---------------------------------------------------------------
1423
! Functions to add a bin with different binID specifier
1424
! ---------------------------------------------------------------
1425
subroutine DS_add_bin_with_IntegerID(dim_name,intID)
1428
! Subroutine arguments
1430
integer, intent(in) :: intID
1431
character(len=*) :: dim_name
1435
call DS_add_bin_with_binID(dim_name,DS_binID(intID))
1436
end subroutine DS_add_bin_with_IntegerID
1438
subroutine DS_add_bin_with_void(dim_name)
1441
! Subroutine arguments
1443
character(len=*) :: dim_name
1447
integer :: ref_size, run_size
1451
ref_size=size(ref_grid(DS_dim_index(ref_grid,dim_name))%bins)
1452
run_size=size(run_grid(DS_dim_index(run_grid,dim_name))%bins)
1453
call DS_add_bin_with_binID(dim_name,DS_binID(
1454
& max(ref_size, run_size)+1))
1455
end subroutine DS_add_bin_with_void
1457
subroutine DS_add_bin_with_binID(dim_name,mBinID)
1460
! Subroutine arguments
1462
type(binID), intent(in) :: mBinID
1463
character(len=*) :: dim_name
1467
type(Bin) :: new_bin
1471
call DS_reinitialize_bin(new_bin)
1472
new_bin%bid = mBinID
1473
call DS_add_bin_to_bins(ref_grid(DS_dim_index(ref_grid,
1474
& dim_name))%bins,new_bin)
1475
call DS_add_bin_to_bins(run_grid(DS_dim_index(run_grid,
1476
& dim_name))%bins,new_bin)
1477
end subroutine DS_add_bin_with_binID
1479
subroutine DS_add_bin_to_bins(bins,new_bin)
1482
! Subroutine arguments
1484
type(Bin), dimension(:), allocatable, intent(inout)
1486
type(Bin) :: new_bin
1490
type(Bin), dimension(:), allocatable :: tmp
1491
integer :: i, bin_index
1495
bin_index = DS_bin_index(bins,new_bin%bid,.True.)
1496
if (bin_index.ne.-1) then
1497
write(*,*)"DiscreteSampler:: Error, the bin with binID '"//
1498
& trim(DS_toStr(new_bin%bid))//"' cannot be added "//
1499
& "be added because it already exists."
1504
allocate(tmp(size(bins)+1))
1506
call DS_copy_bin(bins(i),tmp(i))
1508
tmp(size(bins)+1) = new_bin
1510
allocate(bins(size(tmp)))
1512
call DS_copy_bin(tmp(i),bins(i))
1515
end subroutine DS_add_bin_to_bins
1517
subroutine DS_copy_bin(source, trget)
1519
type(Bin), intent(out) :: trget
1520
type(Bin), intent(in) :: source
1521
trget%weight = source%weight
1522
trget%weight_sqr = source%weight_sqr
1523
trget%abs_weight = source%abs_weight
1524
trget%n_entries = source%n_entries
1525
trget%bid = DS_binID(source%bid%id)
1526
end subroutine DS_copy_bin
1528
! ---------------------------------------------------------------
1529
! Functions to remove a bin from a dimension
1530
! ---------------------------------------------------------------
1531
subroutine DS_remove_bin_withIndex(dim_name, binIndex)
1534
! Subroutine arguments
1536
character(len=*), intent(in) :: dim_name
1537
integer, intent(in) :: binIndex
1542
call DS_remove_bin_from_grid(run_grid(
1543
& DS_dim_index(run_grid, dim_name)),binIndex)
1544
end subroutine DS_remove_bin_withIndex
1546
subroutine DS_remove_bin_withBinID(dim_name, mbinID)
1549
! Subroutine arguments
1551
character(len=*), intent(in) :: dim_name
1552
type(binID), intent(in) :: mbinID
1556
integer :: ref_dim_index,run_dim_index
1557
integer :: ref_bin_index,run_bin_index
1561
ref_dim_index = DS_dim_index(ref_grid, dim_name)
1562
ref_bin_index = DS_bin_index(ref_grid(ref_dim_index)%bins,
1564
call DS_remove_bin_from_grid(ref_grid(ref_dim_index),
1566
run_dim_index = DS_dim_index(run_grid, dim_name)
1567
run_bin_index = DS_bin_index(run_grid(run_dim_index)%bins,
1569
call DS_remove_bin_from_grid(run_grid(run_dim_index),
1571
end subroutine DS_remove_bin_withBinID
1573
subroutine DS_remove_bin_withIntegerID(dim_name, mBinIntID)
1576
! Subroutine arguments
1578
character(len=*), intent(in) :: dim_name
1579
integer, intent(in) :: mBinIntID
1583
call DS_remove_bin_withBinID(dim_name,DS_binID(mBinIntID))
1584
end subroutine DS_remove_bin_withIntegerID
1586
subroutine DS_remove_bin_from_grid(grid, bin_index)
1589
! Subroutine arguments
1591
type(SampledDimension), intent(inout) :: grid
1592
integer, intent(in) :: bin_index
1596
type(Bin), dimension(:), allocatable :: tmp
1602
! Update the norm, norm_sqr and the number of entries in
1603
! the corresponding dimension
1604
grid%norm = grid%norm - grid%bins(bin_index)%weight
1605
grid%norm_sqr = grid%norm_sqr -
1606
& grid%bins(bin_index)%weight_sqr
1607
grid%abs_norm = grid%abs_norm -
1608
& grid%bins(bin_index)%abs_weight
1609
grid%variance_norm = grid%variance_norm
1610
& - DS_bin_variance(grid%bins(bin_index))
1611
grid%n_tot_entries = grid%n_tot_entries
1612
& - grid%bins(bin_index)%n_entries
1613
allocate(tmp(size(grid%bins)-1))
1615
tmp(i) = grid%bins(i)
1617
do i=bin_index+1,size(grid%bins)
1618
tmp(i-1) = grid%bins(i)
1620
deallocate(grid%bins)
1621
allocate(grid%bins(size(tmp)))
1626
end subroutine DS_remove_bin_from_grid
1629
! ---------------------------------------------------------------
1630
! Function to update the reference grid with the running one
1631
! ---------------------------------------------------------------
1632
subroutine DS_update_all_grids(filterZeros)
1635
! Subroutine arguments
1637
logical, optional :: filterZeros
1642
logical :: do_filterZeros
1646
if (.not.allocated(run_grid)) then
1649
if(present(filterZeros)) then
1650
do_filterZeros = filterZeros
1652
do_filterZeros = .False.
1654
do i=1, size(run_grid)
1655
call DS_update_grid_with_dim_index(i,do_filterZeros)
1657
end subroutine DS_update_all_grids
1659
subroutine DS_update_grid_with_dim_name(dim_name, filterZeros)
1662
! Subroutine arguments
1664
character(len=*) :: dim_name
1665
logical, optional :: filterZeros
1670
logical :: do_filterZeros
1674
if(present(filterZeros)) then
1675
do_filterZeros = filterZeros
1677
do_filterZeros = .False.
1679
call DS_update_grid_with_dim_index(
1680
& DS_dim_index(run_grid,dim_name),do_filterZeros)
1682
end subroutine DS_update_grid_with_dim_name
1684
subroutine DS_update_grid_with_dim_index(d_index,filterOutZeros)
1687
! Subroutine arguments
1690
logical :: filterOutZeros
1694
integer :: i, ref_d_index
1695
integer :: ref_bin_index
1697
character, dimension(:), allocatable :: dim_name
1698
type(BinID) :: mBinID
1699
type(Bin) :: new_bin, ref_bin, run_bin
1700
logical :: empty_ref_grid
1704
allocate(dim_name(size(run_grid(d_index)%dimension_name)))
1705
dim_name = run_grid(d_index)%dimension_name
1706
call DS_Logger("Updating dimension '"//
1707
& trim(toStr(dim_name))//"'.")
1709
! Start by making sure that the dimension exists in the
1710
! reference grid. If not, then create it.
1711
if (DS_dim_index(ref_grid,
1712
& run_grid(d_index)%dimension_name,.True.).eq.-1) then
1713
call DS_Logger('Reference grid does not have dimension '//
1714
& trim(toStr(dim_name))//'. Adding it now')
1715
call DS_add_dimension_to_grid(ref_grid,
1716
& trim(toStr(dim_name)) , 0)
1718
ref_d_index = DS_dim_index(ref_grid, dim_name)
1720
empty_ref_grid = (ref_grid(ref_d_index)%n_tot_entries.eq.0)
1722
do i=1,size(run_grid(d_index)%bins)
1723
mBinID = run_grid(d_index)%bins(i)%bid
1724
ref_bin_index = DS_bin_index(
1725
& ref_grid(ref_d_index)%bins,mBinID,.True.)
1726
if (ref_bin_index.eq.-1) then
1727
call DS_Logger('Bin with binID '//trim(DS_toStr(mBinID))//
1728
& ' is missing in the reference grid. Adding it now.')
1729
call DS_reinitialize_bin(new_bin)
1730
new_bin%bid = mBinID
1731
call DS_add_bin_to_bins(ref_grid(ref_d_index)%bins,
1733
ref_bin_index = DS_bin_index(
1734
& ref_grid(ref_d_index)%bins,mBinID)
1737
run_bin = run_grid(d_index)%bins(i)
1738
if ((run_bin%n_entries.lt.ref_grid(ref_d_index)%
1739
& min_bin_probing_points).and.empty_ref_grid) then
1740
write(*,*) "DiscreteSampler:: WARNING, the bin '"//
1741
& trim(DS_toStr(run_bin%bid))//"' of dimension '"//
1742
& trim(toStr(dim_name))//"' will be used for reference"//
1743
& " even though it has been probed only "//
1744
& trim(toStr(run_bin%n_entries))//" times (minimum "//
1745
& "requested is "//trim(toStr(ref_grid(ref_d_index)%
1746
& min_bin_probing_points))//" times)."
1749
ref_bin = ref_grid(ref_d_index)%bins(ref_bin_index)
1750
if (ref_grid(ref_d_index)%grid_mode.eq.2) then
1751
! This means that the reference grid is in 'initialization'
1752
! mode and should be overwritten by the running grid (instead
1753
! of being combined with it) when updated except for the
1754
! bins with not enough entries in the run_grid.
1755
if (run_bin%n_entries.ge.ref_grid(ref_d_index)%
1756
& min_bin_probing_points) then
1757
call DS_reinitialize_bin(ref_bin)
1759
! Then we combine the run_bin and the ref_bin by weighting
1760
! the ref_bin with the ratio of the corresponding norms
1761
ref_bin%weight = ref_bin%weight * (run_grid(
1762
& d_index)%abs_norm / ref_grid(ref_d_index)%abs_norm)
1763
ref_bin%abs_weight = ref_bin%abs_weight * (run_grid(
1764
& d_index)%abs_norm / ref_grid(ref_d_index)%abs_norm)
1765
ref_bin%weight_sqr = ref_bin%weight_sqr * (run_grid(
1766
& d_index)%norm_sqr / ref_grid(ref_d_index)%norm_sqr)
1770
new_bin = ref_bin + run_bin
1772
! Now update the ref grid bin
1773
ref_grid(ref_d_index)%bins(ref_bin_index) = new_bin
1776
call DS_synchronize_grid_with_bins(ref_grid(ref_d_index))
1778
! Now we set the global attribute of the reference_grid to be
1779
! the ones of the running grid.
1780
ref_grid(ref_d_index)%min_bin_probing_points =
1781
& run_grid(d_index)%min_bin_probing_points
1782
ref_grid(ref_d_index)%grid_mode = run_grid(d_index)%grid_mode
1783
ref_grid(ref_d_index)%small_contrib_threshold =
1784
& run_grid(d_index)%small_contrib_threshold
1785
ref_grid(ref_d_index)%damping_power =
1786
& run_grid(d_index)%damping_power
1788
! Now filter all bins in ref_grid that have 0.0 weight and
1789
! remove them! They will not be probed anyway.
1790
if (filterOutZeros) then
1792
do j=1,size(ref_grid(ref_d_index)%bins)
1794
if ((ref_grid(ref_d_index)%bins(i)%weight.eq.0.0d0).and.
1795
& (ref_grid(ref_d_index)%bins(i)%abs_weight.eq.0.0d0).and.
1796
& (ref_grid(ref_d_index)%bins(i)%weight_sqr.eq.0.0d0)) then
1797
call DS_Logger('Bin with binID '//
1798
& trim(DS_toStr(ref_grid(ref_d_index)%bins(i)%bid))//
1799
& ' is zero and will be filtered out. Removing it now.')
1800
call DS_remove_bin_from_grid(ref_grid(ref_d_index),i)
1806
! Clear the running grid now
1807
call DS_reinitialize_dimension(run_grid(d_index))
1809
deallocate(dim_name)
1811
end subroutine DS_update_grid_with_dim_index
1814
function DS_combine_two_bins(BinA, BinB) result(CombinedBin)
1817
! Function arguments
1820
Type(Bin), intent(in) :: BinA, BinB
1821
Type(Bin) :: CombinedBin
1825
call DS_reinitialize_bin(CombinedBin)
1826
if(.not.(BinA%bid==BinB%bid)) then
1827
write(*,*) 'DiscreteSampler:: Error in function '//
1828
& 'DS_combine_two_bins, cannot combine two bins '//
1829
& ' with different bin IDs : '//trim(DS_toStr(BinA%bid))//
1830
& ', '//trim(DS_toStr(BinB%bid))
1833
CombinedBin%bid = BinA%bid
1834
CombinedBin%n_entries = BinA%n_entries + BinB%n_entries
1835
if (CombinedBin%n_entries.eq.0) then
1836
CombinedBin%weight = 0.0d0
1837
CombinedBin%abs_weight = 0.0d0
1838
CombinedBin%weight_sqr = 0.0d0
1840
CombinedBin%weight = (BinA%weight*BinA%n_entries +
1841
& BinB%weight*BinB%n_entries)/CombinedBin%n_entries
1842
CombinedBin%abs_weight = (BinA%abs_weight*BinA%n_entries +
1843
& BinB%abs_weight*BinB%n_entries)/CombinedBin%n_entries
1844
CombinedBin%weight_sqr = (BinA%weight_sqr*BinA%n_entries +
1845
& BinB%weight_sqr*BinB%n_entries)/CombinedBin%n_entries
1847
end function DS_combine_two_bins
1849
! ================================================
1850
! Main function to pick a point
1851
! ================================================
1853
subroutine DS_get_point_with_integerBinID(dim_name,
1854
& random_variable, integerIDPicked, jacobian_weight,mode,
1855
& convoluted_grid_names)
1857
! Subroutine arguments
1859
character(len=*), intent(in) :: dim_name
1860
real*8, intent(in) :: random_variable
1861
integer, intent(out) :: integerIDPicked
1862
real*8, intent(out) :: jacobian_weight
1863
character(len=*), intent(in), optional :: mode
1864
character(len=*), dimension(:), intent(in), optional ::
1865
& convoluted_grid_names
1869
type(BinID) :: mBinID
1873
if (present(mode)) then
1874
if (present(convoluted_grid_names)) then
1875
call DS_get_point_with_BinID(dim_name,random_variable,
1876
& mBinID,jacobian_weight,mode=mode,
1877
& convoluted_grid_names=convoluted_grid_names)
1879
call DS_get_point_with_BinID(dim_name,random_variable,
1880
& mBinID,jacobian_weight,mode=mode)
1883
if (present(convoluted_grid_names)) then
1884
call DS_get_point_with_BinID(dim_name,random_variable,
1885
& mBinID,jacobian_weight,
1886
& convoluted_grid_names=convoluted_grid_names)
1888
call DS_get_point_with_BinID(dim_name,random_variable,
1889
& mBinID,jacobian_weight)
1892
integerIDPicked = mBinID%id
1893
end subroutine DS_get_point_with_integerBinID
1895
subroutine DS_get_point_with_BinID(dim_name,
1896
& random_variable, mBinID, jacobian_weight, mode,
1897
& convoluted_grid_names)
1899
! Subroutine arguments
1901
character(len=*), intent(in) :: dim_name
1902
real*8, intent(in) :: random_variable
1903
type(BinID), intent(out) :: mBinID
1904
real*8, intent(out) :: jacobian_weight
1905
character(len=*), intent(in), optional :: mode
1906
character(len=*), dimension(:), intent(in), optional ::
1907
& convoluted_grid_names
1911
! chose_mode = 1 : Sampling accoridng to variance
1912
! chose_mode = 2 : Sampling according to norm
1913
! chose_mode = 3 : Uniform sampling
1914
integer :: chosen_mode
1915
type(SampledDimension) :: mGrid, runGrid
1916
type(Bin) :: mBin, mRunBin
1917
integer :: ref_grid_index, run_grid_index
1919
real*8 :: running_bound
1920
real*8 :: normalized_bin_bound
1921
logical, dimension(:), allocatable :: bin_indices_to_fill
1922
logical :: initialization_done
1923
real*8 :: sampling_norm
1924
! Local variables related to convolution
1925
real*8, dimension(:), allocatable :: convolution_factors
1926
integer :: conv_bin_index
1927
type(SampledDimension) :: conv_dim
1928
logical :: one_norm_is_zero
1929
real*8 :: small_contrib_thres
1933
if (present(mode)) then
1934
if (mode.eq.'variance') then
1936
elseif (mode.eq.'norm') then
1938
elseif (mode.eq.'uniform') then
1941
write(*,*) "DiscreteSampler:: Error in subroutine"//
1942
& " DS_get_point, mode '"//mode//"' is not recognized."
1949
if (.not.allocated(ref_grid)) then
1950
write(*,*) "DiscreteSampler:: Error, dimensions"//
1951
& " must first be registered with 'DS_register_dimension'"//
1952
& " before the module can be used to pick a point."
1956
ref_grid_index = DS_dim_index(ref_grid, dim_name,.True.)
1957
if (ref_grid_index.eq.-1) then
1958
write(*,*) "DiscreteSampler:: Error in subroutine"//
1959
& " DS_get_point, dimension '"//dim_name//"' not found."
1962
mGrid = ref_grid(ref_grid_index)
1963
run_grid_index = DS_dim_index(run_grid, dim_name,.True.)
1964
if (run_grid_index.eq.-1) then
1965
write(*,*) "DiscreteSampler:: Error in subroutine"//
1966
& " DS_get_point, dimension '"//dim_name//"' not found"//
1967
& " in the running grid."
1970
runGrid = run_grid(run_grid_index)
1972
! If the reference grid is empty, force the use of uniform
1974
if (mGrid%n_tot_entries.eq.0) then
1978
! Pick the right norm for the chosen mode
1979
if (chosen_mode.eq.1) then
1980
sampling_norm = mGrid%variance_norm
1981
elseif (chosen_mode.eq.2) then
1982
sampling_norm = mGrid%abs_norm
1983
elseif (chosen_mode.eq.3) then
1984
sampling_norm = float(size(mGrid%bins))
1987
! If the grid is empty we must first make sure that each bin was
1988
! probed with min_bin_probing_points before using a uniform grid
1989
allocate(bin_indices_to_fill(size(mGrid%bins)))
1990
initialization_done = .True.
1991
if(mGrid%n_tot_entries.eq.0) then
1993
do i=1,size(mGrid%bins)
1994
mRunBin = DS_get_bin(runGrid%bins,mGrid%bins(i)%bid)
1995
if (mRunBin%n_entries.lt.mGrid%min_bin_probing_points) then
1996
bin_indices_to_fill(i) = .True.
1997
initialization_done = .False.
1999
bin_indices_to_fill(i) = .False.
2002
if(.not.initialization_done) then
2003
! In this case, we will only fill in bins which do not have
2004
! have enough entries (and select them uniformly) and veto the
2005
! others. The jacobian returned is still the one corresponding
2006
! to a uniform distributions over the whole set of bins.
2007
! Possible convolutions are ignored
2008
sampling_norm = 0.0d0
2009
do i=1,size(bin_indices_to_fill)
2010
if (bin_indices_to_fill(i)) then
2011
sampling_norm = sampling_norm + 1.0d0
2017
if (initialization_done) then
2018
do i=1,size(mGrid%bins)
2019
bin_indices_to_fill(i) = .True.
2023
! Pick the right reference bin value for the chosen mode. Note
2024
! that this reference value is stored in the %weight attribute
2025
! of the reference grid local copy mGrid
2026
do i=1,size(mGrid%bins)
2027
if (.not.bin_indices_to_fill(i)) then
2028
mGrid%bins(i)%weight = 0.0d0
2029
elseif (chosen_mode.eq.1) then
2030
mGrid%bins(i)%weight = DS_bin_variance(mGrid%bins(i))
2031
elseif (chosen_mode.eq.2) then
2032
mGrid%bins(i)%weight = mGrid%bins(i)%abs_weight
2033
elseif (chosen_mode.eq.3) then
2034
mGrid%bins(i)%weight = 1.0d0
2039
! Treat specially contributions worth less than 5% of the
2040
! contribution averaged over all bins. For those, we sample
2041
! according to the square root (or the specified power 'pow'
2042
! of the reference value corresponding to the chosen mode.
2043
! In this way, we are less sensitive to possible large fluctuations
2044
! of very suppressed contributions for which the Jacobian would be
2045
! really big. However, the square-root is such that a really
2046
! suppressed contribution at the level of numerical precision
2047
! would still never be probed.
2049
average_contrib = sampling_norm / size(mGrid%bins)
2050
! Ignore this if the average contribution is zero
2051
if (average_contrib.gt.0.0d0) then
2052
do i=1,size(mGrid%bins)
2053
mBin = mGrid%bins(i)
2054
if ( (mBin%weight/average_contrib) .lt.
2055
& runGrid%small_contrib_threshold) then
2056
sampling_norm = sampling_norm - mGrid%bins(i)%weight
2057
mGrid%bins(i)%weight =
2058
& ((mBin%weight/(runGrid%small_contrib_threshold
2059
& *average_contrib))**runGrid%damping_power)*
2060
& runGrid%small_contrib_threshold*average_contrib
2061
sampling_norm = sampling_norm + mGrid%bins(i)%weight
2066
! Now appropriately set the convolution factors
2068
allocate(convolution_factors(size(mGrid%bins)))
2069
if (present(convoluted_grid_names).and.initialization_done) then
2071
do j=1,size(convoluted_grid_names)
2072
if (DS_dim_index(run_grid,convoluted_grid_names(j),
2073
& .True.).eq.-1) then
2074
write(*,*) "DiscreteSampler:: Error, dimension '"//
2075
& convoluted_grid_names(j)//"' for convolut"//
2076
& "ion could not be found in the running grid."
2080
sampling_norm = 0.0d0
2081
do i=1,size(mGrid%bins)
2082
convolution_factors(i) = 1.0d0
2083
do j=1,size(convoluted_grid_names)
2084
conv_dim = DS_get_dimension(
2085
& run_grid,convoluted_grid_names(j))
2086
conv_bin_index = DS_bin_index(conv_dim%bins,
2087
& mGrid%bins(i)%bid,.True.)
2088
if (conv_bin_index.eq.-1) then
2089
write(*,*) "DiscreteSampler:: Error, bin '"//
2090
& trim(DS_toStr(mGrid%bins(i)%bid))//"' could not be fo"//
2091
& "und in convoluted dimension '"//
2092
& convoluted_grid_names(j)//"'."
2095
! Notice that for the convolution we always use the
2096
! absolute value of the weight because we assume the user
2097
! has edited this grid by himself for with a single entry.
2098
convolution_factors(i) = convolution_factors(i)*
2099
& conv_dim%bins(conv_bin_index)%abs_weight
2101
sampling_norm = sampling_norm +
2102
& convolution_factors(i)*mGrid%bins(i)%weight
2105
do i=1,size(mGrid%bins)
2106
convolution_factors(i) = 1.0d0
2110
! Now crash nicely on zero norm grid
2111
if (sampling_norm.eq.0d0.and..not.DS_tolerate_zero_norm) then
2112
one_norm_is_zero = .FALSE.
2113
write(*,*) 'DiscreteSampler:: Error, all bins'//
2114
& " of sampled dimension '"//dim_name//"' or of the"//
2115
& " following convoluted dimensions have zero weight:"
2116
if (chosen_mode.eq.2) then
2117
write(*,*) "DiscreteSampler:: Sampled dimension "//
2118
& " : '"//trim(toStr(mGrid%dimension_name))//"' with norm "//
2119
& trim(toStr(mGrid%abs_norm,'ENw.3'))//"."
2120
one_norm_is_zero = (one_norm_is_zero.or.
2121
& mGrid%abs_norm.eq.0.0d0)
2122
elseif (chosen_mode.eq.1) then
2123
write(*,*) "DiscreteSampler:: Sampled dimension "//
2124
& " : '"//trim(toStr(mGrid%dimension_name))//"' with norm "//
2125
& trim(toStr(mGrid%variance_norm,'ENw.3'))//"."
2126
one_norm_is_zero = (one_norm_is_zero.or.
2127
& mGrid%variance_norm.eq.0.0d0)
2128
elseif (chosen_mode.eq.3) then
2129
write(*,*) "DiscreteSampler:: Norm of sampled dimension '"//
2130
& trim(toStr(mGrid%dimension_name))//"' irrelevant since"//
2131
& " uniform sampling was selected."
2133
if(present(convoluted_grid_names).and.initialization_done)then
2134
do i=1,size(convoluted_grid_names)
2135
conv_dim = DS_get_dimension(run_grid,
2136
& convoluted_grid_names(i))
2137
write(*,*) "DiscreteSampler:: Convoluted dimension "//
2138
& trim(toStr(i))//": '"//convoluted_grid_names(i)//
2139
& "' with norm "//trim(toStr(conv_dim%abs_norm,'ENw.3'))//"."
2140
one_norm_is_zero = (one_norm_is_zero.or.
2141
& conv_dim%abs_norm.eq.0.0d0)
2144
if(present(convoluted_grid_names).and.initialization_done
2145
& .and.(.not.one_norm_is_zero))then
2146
write(*,*) "DiscreteSampler:: None of the norm above"
2147
& //" is zero, this means that the convolution (product)"
2148
& //" of the grids yields zero for each bin, even though"
2149
& //" they are not zero separately."
2150
write(*,*) "DiscreteSampler:: Use DS_print_global_info()"//
2151
& " to investigate further."
2153
write(*,*) "DiscreteSampler:: One norm is zero, no sampling"//
2154
& " can be done in these conditions. Set 'tolerate_zero_norm"//
2155
& "' to .True. when initializating the module to proceed wi"//
2156
& "th a uniform distribution for the grids of zero norm."
2160
! Or make it pure random if DS_tolerate_zero_norm is True.
2161
if (sampling_norm.eq.0d0) then
2162
do i=1,size(mGrid%bins)
2163
bin_indices_to_fill(i) = .True.
2164
if(chosen_mode.eq.2.and.mGrid%abs_norm.eq.0.0d0.or.
2165
& chosen_mode.eq.1.and.mGrid%variance_norm.eq.0.0d0) then
2166
mGrid%bins(i)%weight = 1.0d0
2168
if (present(convoluted_grid_names).and.
2169
& initialization_done.and.conv_dim%abs_norm.eq.0.0d0) then
2170
conv_dim = DS_get_dimension(run_grid,
2171
& convoluted_grid_names(i))
2172
if (conv_dim%abs_norm.eq.0.0d0) then
2173
convolution_factors(i) = 1.0d0
2176
sampling_norm = sampling_norm +
2177
& mGrid%bins(i)%weight*convolution_factors(i)
2179
! If sampling_norm is again zero it means that the two grids
2180
! are "orthogonal" so that we have no choice but to randomize
2182
if (sampling_norm.eq.0.0d0) then
2183
do i=1,size(mGrid%bins)
2184
mGrid%bins(i)%weight = 1.0d0
2185
convolution_factors(i) = 1.0d0
2186
sampling_norm = sampling_norm + 1.0d0
2192
! Now come the usual sampling method
2194
running_bound = 0.0d0
2195
do i=1,size(mGrid%bins)
2196
if (.not.bin_indices_to_fill(i)) then
2199
mBin = mGrid%bins(i)
2200
normalized_bin_bound = mBin%weight *
2201
& ( convolution_factors(i) / sampling_norm )
2202
running_bound = running_bound + normalized_bin_bound
2203
if (random_variable.le.running_bound) then
2204
mBinID = mGrid%bins(i)%bid
2205
jacobian_weight = 1.0d0 / normalized_bin_bound
2206
deallocate(convolution_factors)
2207
deallocate(bin_indices_to_fill)
2211
! If no point was picked at this stage, there was a problem
2212
write(*,*) 'DiscreteSampler:: Error, no point could be '//
2213
& 'picked with random variable '//trim(toStr(random_variable))//
2214
& ' using upper bound found of '//trim(toStr(running_bound))//'.'
2216
end subroutine DS_get_point_with_BinID
2218
function DS_bin_variance(mBin)
2220
! Function arguments
2222
type(Bin), intent(in) :: mBin
2223
real*8 :: DS_bin_variance
2227
DS_bin_variance = ((mBin%weight_sqr - mBin%weight**2) *
2228
& (mBin%n_entries))/(mBin%n_entries+1)
2229
end function DS_bin_variance
2230
! ================================================
2231
! Grid I/O functions
2232
! ================================================
2234
! ---------------------------------------------------------------
2235
! This function writes the ref_grid to a file specified by its
2237
! ---------------------------------------------------------------
2238
subroutine DS_write_grid_with_filename(filename, dim_name,
2242
! Subroutine arguments
2244
character(len=*), intent(in) :: filename
2245
character(len=*), intent(in), optional :: dim_name
2246
character(len=*), intent(in), optional :: grid_type
2254
inquire(file=filename, exist=fileExist)
2256
call DS_Logger('DiscreteSampler:: The file '
2257
& //filename//' already exists, so beware that '//
2258
& ' the grid information will be appended to it.')
2260
open(123, file=filename, err=11, access='append',
2264
write(*,*) 'DiscreteSampler :: Error, file '//filename//
2265
& ' could not be opened for writing.'
2268
if (present(dim_name)) then
2269
if (present(grid_type)) then
2270
call DS_write_grid_with_streamID(123, dim_name, grid_type)
2272
call DS_write_grid_with_streamID(123, dim_name)
2275
if (present(grid_type)) then
2276
call DS_write_grid_with_streamID(123, grid_type=grid_type)
2278
call DS_write_grid_with_streamID(123)
2282
end subroutine DS_write_grid_with_filename
2284
! ---------------------------------------------------------------
2285
! This function writes the ref_grid or all grids to a file
2286
! specified by its stream ID.
2287
! ---------------------------------------------------------------
2288
subroutine DS_write_grid_with_streamID(streamID, dim_name,
2292
! Subroutine arguments
2294
integer, intent(in) :: streamID
2295
character(len=*), intent(in), optional :: dim_name
2296
character(len=*), intent(in), optional :: grid_type
2300
type(SampledDimension) :: grid
2302
integer :: chosen_grid
2306
if (present(grid_type)) then
2307
if (grid_type.eq.'ref') then
2309
elseif (grid_type.eq.'run') then
2311
elseif (grid_type.eq.'all') then
2314
write(*,*) 'DiscreteSampler:: Error in'//
2315
& " subroutine 'DS_write_grid_with_streamID',"//
2316
& " argument grid_type='"//grid_type//"' not"//
2323
if ((chosen_grid.eq.1.or.chosen_grid.eq.3)
2324
& .and..not.allocated(ref_grid)) then
2327
if ((chosen_grid.eq.2..or.chosen_grid.eq.3)
2328
& .and..not.allocated(run_grid)) then
2331
if (present(dim_name)) then
2332
if (chosen_grid.eq.1.or.chosen_grid.eq.3) then
2333
grid = ref_grid(DS_dim_index(ref_grid, dim_name))
2334
call DS_write_grid_from_grid(grid, streamID,'ref')
2336
if (chosen_grid.eq.2.or.chosen_grid.eq.3) then
2337
grid = run_grid(DS_dim_index(run_grid, dim_name))
2338
call DS_write_grid_from_grid(grid, streamID,'run')
2341
if (chosen_grid.eq.1.or.chosen_grid.eq.3) then
2342
do i=1,size(ref_grid)
2344
call DS_write_grid_from_grid(grid, streamID,'ref')
2347
if (chosen_grid.eq.2.or.chosen_grid.eq.3) then
2348
do i=1,size(run_grid)
2350
call DS_write_grid_from_grid(grid, streamID,'run')
2354
end subroutine DS_write_grid_with_streamID
2356
! ---------------------------------------------------------------
2357
! This function writes a given grid to a file.
2358
! ---------------------------------------------------------------
2359
subroutine DS_write_grid_from_grid(grid, streamID, grid_type)
2362
! Subroutine arguments
2364
integer, intent(in) :: streamID
2365
type(SampledDimension), intent(in) :: grid
2366
character(len=*), intent(in) :: grid_type
2375
write(streamID,*) ' <DiscreteSampler_grid>'
2376
write(streamID,*) ' '//trim(toStr(grid%dimension_name))
2377
if (grid_type.eq.'ref') then
2378
write(streamID,*) ' '//trim(toStr(1))
2379
& //" # 1 for a reference and 2 for a running grid."
2380
elseif (grid_type.eq.'run') then
2381
write(streamID,*) ' '//trim(toStr(2))
2382
& //" # 1 for a reference and 2 for a running grid."
2384
write(*,*) "DiscreteSampler:: Error, grid_type'"//
2385
& grid_type//"' not recognized."
2388
write(streamID,*) ' '//trim(toStr(grid%min_bin_probing_points
2389
& ))//" # Attribute 'min_bin_probing_points' of the grid."
2390
write(streamID,*) ' '//trim(toStr(grid%grid_mode
2391
& ))//" # Attribute 'grid_mode' of the grid. 1=='default',"
2392
& //"2=='initialization'"
2393
write(streamID,*) ' '//trim(toStr(grid%small_contrib_threshold
2394
& ))//" # Attribute 'small_contrib_threshold' of the grid."
2395
write(streamID,*) ' '//trim(toStr(grid%damping_power
2396
& ))//" # Attribute 'damping_power' of the grid."
2397
write(streamID,*) '# binID n_entries weight weight_sqr'//
2399
do i=1,size(grid%bins)
2401
& ' '//trim(DS_toStr(grid%bins(i)%bid))//
2402
& ' '//trim(toStr(grid%bins(i)%n_entries))//
2403
& ' '//trim(toStr(grid%bins(i)%weight,'ESw.15E3'))//
2404
& ' '//trim(toStr(grid%bins(i)%weight_sqr,'ESw.15E3'))//
2405
& ' '//trim(toStr(grid%bins(i)%abs_weight,'ESw.15E3'))
2407
write(streamID,*) ' </DiscreteSampler_grid>'
2409
end subroutine DS_write_grid_from_grid
2411
! ---------------------------------------------------------------
2412
! This function loads the grid specified in a file specified by its
2413
! stream ID into the run_grid.
2414
! ---------------------------------------------------------------
2415
subroutine DS_load_grid_with_filename(filename, dim_name)
2418
! Subroutine arguments
2420
character(len=*), intent(in) :: filename
2421
character(len=*), intent(in), optional :: dim_name
2429
! Make sure the module is initialized
2430
if (.not.allocated(DS_isInitialized)) then
2431
call DS_initialize()
2433
inquire(file=filename, exist=fileExist)
2434
if (.not.fileExist) then
2435
write(*,*) 'DiscreteSampler:: Error, the file '//filename//
2436
& ' could not be found.'
2439
open(124, file=filename, err=13, action='read')
2442
write(*,*) 'DiscreteSampler :: Error, file '//filename//
2443
& ' exists but could not be read.'
2445
if (present(dim_name)) then
2446
call DS_load_grid_with_streamID(124, dim_name)
2448
call DS_load_grid_with_streamID(124)
2451
end subroutine DS_load_grid_with_filename
2453
! ---------------------------------------------------------------
2454
! This function loads the grid specified in a file specified by its
2455
! stream ID into the run_grid.
2456
! ---------------------------------------------------------------
2457
subroutine DS_load_grid_with_streamID(streamID, dim_name)
2460
! Subroutine arguments
2462
integer, intent(in) :: streamID
2463
character(len=*), intent(in), optional :: dim_name
2468
character(512) :: buff
2469
character(2) :: TwoBuff
2470
character(3) :: ThreeBuff
2471
logical :: startedGrid
2472
real*8 :: weight, abs_weight, weight_sqr
2473
integer :: n_entries, bid
2474
type(Bin) :: new_bin
2475
integer :: char_size
2476
integer :: read_position
2477
integer :: run_dim_index
2478
integer :: grid_mode
2479
real*8 :: small_contrib_threshold
2480
real*8 :: damping_power
2484
! Make sure the module is initialized
2485
if (.not.allocated(DS_isInitialized)) then
2486
call DS_initialize()
2488
! Now start reading the file
2489
startedGrid = .False.
2493
read(streamID, "(A)", size=char_size, eor=998,
2494
& end=999, advance='no') TwoBuff
2497
if (char_size.le.1) then
2500
if (TwoBuff(1:1).eq.'#'.or.TwoBuff(2:2).eq.'#') then
2501
! Advance the stream
2502
read(streamID,*,end=990) buff
2505
if (startedGrid) then
2506
read(streamID, "(A)", size=char_size,
2507
& end=999, advance='no') TwoBuff
2508
if (TwoBuff(1:2).eq.'</') then
2509
! Advance the stream
2510
read(streamID,*,end=990) buff
2511
startedGrid = .False.
2515
read(streamID,*,end=990) bid, n_entries, weight,
2516
& weight_sqr, abs_weight
2518
new_bin%n_entries = n_entries
2519
new_bin%weight = weight
2520
new_bin%weight_sqr = weight_sqr
2521
new_bin%abs_weight = abs_weight
2522
call DS_add_bin_to_bins(run_grid(size(run_grid))%bins,
2525
! Advance the stream
2526
if (read_position.eq.0) then
2527
read(streamID,*,end=990) buff
2528
if (buff(1:22).eq.'<DiscreteSampler_grid>') then
2529
read_position = read_position + 1
2532
select case(read_position)
2534
read(streamID,*,end=990) buff
2535
run_dim_index = DS_dim_index(run_grid,
2536
& trim(buff),.True.)
2537
if (run_dim_index.ne.-1) then
2538
call DS_remove_dimension_from_grid(run_grid,
2541
call DS_register_dimension(trim(buff),0,.False.)
2543
read(streamID,*,end=990) grid_mode
2544
if (grid_mode.ne.1) then
2545
write(*,*) 'DiscreteSampler:: Warning, the '//
2546
& "grid read is not of type 'reference'."//
2547
& " It will be skipped."
2548
call DS_remove_dimension_from_grid(run_grid,
2551
startedGrid = .False.
2555
read(streamID,*,end=990)
2556
& run_grid(size(run_grid))%min_bin_probing_points
2558
read(streamID,*,end=990)
2559
& run_grid(size(run_grid))%grid_mode
2561
read(streamID,*,end=990) small_contrib_threshold
2562
if (small_contrib_threshold.lt.0.0d0.or.
2563
& small_contrib_threshold.gt.0.5d0) then
2564
write(*,*) 'DiscreteSampler:: The '//
2565
& 'small_contrib_threshold must be >= 0.0 and '//
2566
& '< 0.5 to be meaningful.'
2569
run_grid(size(run_grid))%small_contrib_threshold
2570
& = small_contrib_threshold
2572
read(streamID,*,end=990) damping_power
2573
if (damping_power.lt.0.0d0.or.
2574
& damping_power.gt.1.0d0) then
2575
write(*,*) 'DiscreteSampler:: The damping power'//
2576
& ' must be >= 0.0 and <= 1.0.'
2579
run_grid(size(run_grid))%damping_power
2581
! Make sure that the last info read before reading the
2582
! bin content (here the info with read_position=6)
2583
! sets startedGrid to .True. to start the bin readout
2584
startedGrid = .True.
2586
write(*,*) 'DiscreteSampler:: Number of entries'//
2587
& ' before reaching bin lists exceeded.'
2590
read_position = read_position + 1
2596
write(*,*) 'DiscreteSampler:: Error, when loading grids'//
2601
! Now update the running grid into the reference one
2602
call DS_update_grid()
2603
end subroutine DS_load_grid_with_streamID
2606
! ---------------------------------------------------------------
2607
! Synchronizes the cumulative information in a given grid from
2609
! ---------------------------------------------------------------
2610
subroutine DS_synchronize_grid_with_bins(grid)
2613
! Subroutine argument
2615
type(sampledDimension), intent(inout) :: grid
2619
real*8 :: norm, abs_norm, norm_sqr, variance_norm
2620
integer :: i, n_tot_entries
2627
variance_norm = 0.0d0
2629
do i=1,size(grid%bins)
2630
n_tot_entries = n_tot_entries + grid%bins(i)%n_entries
2631
norm_sqr = norm_sqr + grid%bins(i)%weight_sqr
2632
abs_norm = abs_norm + grid%bins(i)%abs_weight
2633
norm = norm + grid%bins(i)%weight
2634
variance_norm = variance_norm +
2635
& DS_bin_variance(grid%bins(i))
2637
grid%n_tot_entries = n_tot_entries
2638
grid%norm_sqr = norm_sqr
2639
grid%abs_norm = abs_norm
2641
grid%variance_norm = variance_norm
2642
end subroutine DS_synchronize_grid_with_bins
2644
! ================================================
2645
! Functions and subroutine handling derived types
2646
! ================================================
2648
! ---------------------------------------------------------------
2649
! Specify how bin idea should be compared
2650
! ---------------------------------------------------------------
2651
function equal_binID(binID1,binID2)
2654
! Function arguments
2656
type(binID), intent(in) :: binID1, binID2
2657
logical :: equal_binID
2661
if(binID1%id.ne.binID2%id) then
2662
equal_binID = .False.
2665
equal_binID = .True.
2667
end function equal_binID
2669
! ---------------------------------------------------------------
2670
! BinIDs constructors
2671
! ---------------------------------------------------------------
2672
pure elemental subroutine binID_from_binID(binID1,binID2)
2675
! Function arguments
2677
type(binID), intent(out) :: binID1
2678
type(binID), intent(in) :: binID2
2682
binID1%id = binID2%id
2683
end subroutine binID_from_binID
2685
pure elemental subroutine binID_from_integer(binID1,binIDInt)
2688
! Function arguments
2690
type(binID), intent(out) :: binID1
2691
integer, intent(in) :: binIDInt
2695
binID1%id = binIDInt
2696
end subroutine binID_from_integer
2698
! Provide a constructor-like way of creating a binID
2699
function DS_binID(binIDInt)
2702
! Function arguments
2704
type(binID) :: DS_binID
2705
integer, intent(in) :: binIDInt
2710
end function DS_binID
2711
! ---------------------------------------------------------------
2712
! String representation of a binID
2713
! ---------------------------------------------------------------
2714
function DS_toStr(mBinID)
2717
! Function arguments
2719
type(binID), intent(in) :: mBinID
2720
character(100) :: DS_toStr
2724
DS_toStr = trim(toStr(mBinID%id))
2725
end function DS_toStr
2728
! ================================================
2729
! Access routines emulating a dictionary
2730
! ================================================
2732
! ---------------------------------------------------------------
2733
! Returns the index of the discrete dimension with name dim_name
2734
! ---------------------------------------------------------------
2735
function DS_dim_index_default(grid, dim_name)
2738
! Function arguments
2740
type(sampledDimension), dimension(:), intent(in), allocatable
2742
character(len=*), intent(in) :: dim_name
2743
integer :: DS_dim_index_default
2747
DS_dim_index_default =
2748
& DS_dim_index_with_force(grid, dim_name, .False.)
2749
end function DS_dim_index_default
2751
function DS_dim_index_with_force(grid, dim_name, force)
2754
! Function arguments
2756
type(sampledDimension), dimension(:), intent(in), allocatable
2758
character(len=*), intent(in) :: dim_name
2759
integer :: DS_dim_index_with_force
2769
DS_dim_index_with_force = -1
2770
if (.not.allocated(grid)) then
2773
do i = 1, size(grid)
2774
if (len(dim_name).ne.size(grid(i)%dimension_name)) cycle
2775
do j =1, len(dim_name)
2776
if(grid(i)%dimension_name(j).ne.dim_name(j:j)) then
2780
DS_dim_index_with_force = i
2784
if (DS_dim_index_with_force.eq.-1.and.(.not.force)) then
2785
write(*,*) 'DiscreteSampler:: Error in function dim_index'//
2786
& "(), dimension name '"//dim_name//"' not found."
2789
end function DS_dim_index_with_force
2791
function DS_dim_index_default_with_chararray(grid, dim_name)
2794
! Function arguments
2796
type(sampledDimension), dimension(:), intent(in), allocatable
2798
character, dimension(:), intent(in) :: dim_name
2799
integer :: DS_dim_index_default_with_chararray
2803
DS_dim_index_default_with_chararray =
2804
& DS_dim_index_with_force_with_chararray(
2805
& grid, dim_name, .False.)
2806
end function DS_dim_index_default_with_chararray
2808
function DS_dim_index_with_force_with_chararray(
2809
& grid, dim_name, force)
2812
! Function arguments
2814
type(sampledDimension), dimension(:), intent(in), allocatable
2816
character, dimension(:), intent(in) :: dim_name
2817
integer :: DS_dim_index_with_force_with_chararray
2827
DS_dim_index_with_force_with_chararray = -1
2828
if (.not.allocated(grid)) then
2831
do i = 1, size(grid)
2832
if (size(dim_name).ne.size(grid(i)%dimension_name)) cycle
2833
do j =1, size(dim_name)
2834
if(grid(i)%dimension_name(j).ne.dim_name(j)) then
2838
DS_dim_index_with_force_with_chararray = i
2842
if (DS_dim_index_with_force_with_chararray.eq.-1.and.
2843
& (.not.force)) then
2844
write(*,*) 'DiscreteSampler:: Error in function dim_index'//
2845
& "(), dimension name '"//dim_name//"' not found."
2848
end function DS_dim_index_with_force_with_chararray
2851
end module DiscreteSampler
added
removed
vendor/DiscreteSampler/DiscreteSampler.f
