ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation.
pm_sampleMean::setMean Interface Reference

Return the (weighted) mean of a pair of time series or of an input sample of nsam observations with ndim = 1 or 2 attributes, optionally weighted by the input weight, optionally also sum(weight) and optionally, sum(weight**2).
More...

Detailed Description

Return the (weighted) mean of a pair of time series or of an input sample of nsam observations with ndim = 1 or 2 attributes, optionally weighted by the input weight, optionally also sum(weight) and optionally, sum(weight**2).

This generic interface is developed to specifically improve the performance of other procedures where sum(weight) or sum(weight**2) are also needed in addition to the mean of the weighted sample.
For example, this situation occurs frequently in the computation of the variance or covariance matrix of a weighted sample.

Parameters
[out]mean: The output object of,
  1. type complex of kind any supported by the processor (e.g., CK, CK32, CK64, or CK128),
  2. type real of kind any supported by the processor (e.g., RK, RK32, RK64, or RK128),
containing the data mean.
  1. When the input sample has the shape (nsam), the output mean must be a scalar.
  2. When the input x and y are present, the output mean must be of shape mean(1:2).
  3. When the input sample has the shape (nsam, ndim) or (ndim, nsam), the output mean must be a vector of length ndim.
[in]x: The contiguous array of shape (nsam) of the same type and kind as the output mean, containing the first of the two data series whose mean is to be computed.
(optional. It must be present *if and only if** the input argument sample is missing and y is present.)
[in]y: The contiguous of the same type and kind and shape as the input x, containing the second of the two data series whose mean is to be computed.
(optional. It must be present *if and only if** the input argument sample is missing and x is present.)
[in]sample: The contiguous array of shape (nsam), (ndim, nsam) or (nsam, ndim) of the same type and kind as the output mean, containing the sample whose mean is to be computed.
(optional. It must be present *if and only if** the input arguments x and y are missing.)
[in]dim: The input scalar integer of default kind IK representing the dimension (1 or 2) of the input sample along which the mean must be computed.
  1. If dim = 1, the input sample of rank 2 is assumed to have the shape (nsam, ndim).
  2. If dim = 2, the input sample of rank 2 is assumed to have the shape (ndim, nsam).
The input dim must always be 1 or missing for an input sample of rank 1.
(optional. It can be present only if sample is also present. If missing, the mean of the whole input sample is computed.)
[in]weight: The contiguous vector of length nsam of
  1. type integer of default kind IK, or
  2. type real of the same kind as that of the output mean,
containing the corresponding weights of the data points in the input sample or the input pair x and y.
[out]weisum: The output scalar of the same type and kind as the input weight.
On output, it will contain sum(weight).
This quantity is frequently needed in computing the weighted sample variance.
(optional. It must be present if and only if the input argument weight is also present.)


Possible calling interfaces

! XY sample.
call setMean(mean(1:2), x(1 : nsam), y(1 : nsam))
call setMean(mean(1:2), x(1 : nsam), y(1 : nsam), weight(1 : nsam), weisum)
! 1D sample.
call setMean(mean, sample(1 : nsam))
call setMean(mean, sample(1 : nsam), weight(1 : nsam), weisum)
call setMean(mean, sample(1 : nsam), dim)
call setMean(mean, sample(1 : nsam), dim, weight(1 : nsam), weisum)
! 2D sample.
call setMean(mean(1 : ndim), sample(:,:))
call setMean(mean(1 : ndim), sample(:,:), weight(1 : nsam), weisum)
call setMean(mean(1 : ndim), sample(:,:), dim)
call setMean(mean(1 : ndim), sample(:,:), dim, weight(1 : nsam), weisum)
Return the (weighted) mean of a pair of time series or of an input sample of nsam observations with n...
This module contains classes and procedures for computing the first moment (i.e., the statistical mea...
Warning
The condition all(0. <= weight) must hold for the corresponding input arguments.
The condition 1 <= dim .and. dim <= rank(sample) must hold for the corresponding input arguments.
The condition size(sample, dim) == size(weight, 1) must hold for the corresponding input arguments.
The condition (dim == 1 .and. size(mean, 1) == size(sample, 2)) .or. (dim == 2 .and. size(mean, 1) == size(sample, 1)) must hold for the corresponding input arguments.
The condition size(x) == size(weight) must hold for the corresponding input arguments.
The condition size(x) == size(y) must hold for the corresponding input arguments.
These conditions are verified only if the library is built with the preprocessor macro CHECK_ENABLED=1.
The pure procedure(s) documented herein become impure when the ParaMonte library is compiled with preprocessor macro CHECK_ENABLED=1.
By default, these procedures are pure in release build and impure in debug and testing builds.
Note
If the input sample is to be an array of type integer, simply convert the sample to an array of type real of the desired kind for the output real mean of the sample.
There is no point in accepting an input sample of type integer since it will be inevitably converted to an array of type real within the procedure to avoid potential integer overflow.
Furthermore, an sample of type integer creates ambiguity about the kind of the real-valued returned mean by the procedure.
See the notes in the description of the pm_sampleMean.
Note that the mean of any one or two-dimensional sample can be simply computed via the Fortran intrinsic routine sum():
integer :: i
integer , parameter :: NDIM = 3_IK
integer , parameter :: NSAM = 1000_IK
real , parameter :: sample(NDIM,NSAM) = reshape([( real(i,RK), i = 1, NSAM )], shape = shape(sample))
real , allocatable :: mean(:)
mean = sum(sample, dim = 1) / size(transpose(sample), dim = 1) ! assuming the first dimension represents the observations.
mean = sum(sample, dim = 2) / size(sample, dim = 2) ! assuming the second dimension represents the observations.
The mean of a whole multidimensional array can be obtained by either,
  1. reshaping the array to a vector form and passing it to this procedure, or
  2. mapping the array to a 1-dimensional pointer array of the same size as the ndim dimensional array.
See the examples below.
Developer Remark:
The logic behind allowing a pair of time series data x and y is to allow fast computation of the mean of the pair and the weight sums in one loop.
This pattern occurs frequently in the computation of correlation coefficients.
See also
getVar
getMean
setVarMean


Example usage

1program example
2
3 use pm_kind, only: SK, IK
4 use pm_kind, only: TKG => RKS ! All other real types are also supported.
6 use pm_distUnif, only: getUnifRand
7 use pm_sampleMean, only: setMean
8 use pm_io, only: display_type
9
10 implicit none
11
12 integer(IK) :: idim, ndim, nsam
13 type(display_type) :: disp
14 disp = display_type(file = "main.out.F90")
15
16 call disp%skip()
17 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
18 call disp%show("!Compute the mean of a 1-D array.")
19 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
20 call disp%skip()
21
22 block
23 real(TKG) :: mean
24 real(TKG), allocatable :: sample(:)
25 call disp%skip()
26 call disp%show("nsam = getUnifRand(1, 5)")
27 nsam = getUnifRand(1, 5)
28 call disp%show("sample = getUnifRand(0., 1., nsam)")
29 sample = getUnifRand(0., 1., nsam)
30 call disp%show("sample")
31 call disp%show( sample )
32 call disp%show("call setMean(mean, sample)")
33 call setMean(mean, sample)
34 call disp%show("mean")
35 call disp%show( mean )
36 call disp%show("call setMean(mean, sample, dim = 1_IK)")
37 call setMean(mean, sample, dim = 1_IK)
38 call disp%show("mean")
39 call disp%show( mean )
40 call disp%skip()
41 end block
42
43 call disp%skip()
44 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
45 call disp%show("!Compute the mean of a pair of data series.")
46 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
47 call disp%skip()
48
49 block
50 real(TKG) :: mean(2)
51 real(TKG), allocatable :: sample(:,:)
52 real(TKG), allocatable :: weight(:)
53 real(TKG) :: weisum
54 call disp%skip()
55 call disp%show("nsam = getUnifRand(1, 5)")
56 nsam = getUnifRand(1, 5)
57 call disp%show("sample = getUnifRand(0., 1., nsam, 2_IK)")
58 sample = getUnifRand(0., 1., nsam, 2_IK)
59 call disp%show("sample")
60 call disp%show( sample )
61 call disp%show("call setMean(mean, sample(:,1), sample(:,2))")
62 call setMean(mean, sample(:,1), sample(:,2))
63 call disp%show("mean")
64 call disp%show( mean )
65 call disp%show("call setMean(mean, sample, dim = 1_IK) ! for comparison.")
66 call setMean(mean, sample, dim = 1_IK) ! for comparison.
67 call disp%show("mean")
68 call disp%show( mean )
69 call disp%skip()
70 call disp%show("weight = getUnifRand(1, 100, nsam)")
71 weight = getUnifRand(1, 100, nsam)
72 call disp%show("weight")
73 call disp%show( weight )
74 call disp%show("call setMean(mean, sample(:,1), sample(:,2), weight, weisum)")
75 call setMean(mean, sample(:,1), sample(:,2), weight, weisum)
76 call disp%show("weisum")
77 call disp%show( weisum )
78 call disp%show("mean")
79 call disp%show( mean )
80 call disp%show("call setMean(mean, sample, 1_IK, weight, weisum) ! for comparison.")
81 call setMean(mean, sample, 1_IK, weight, weisum) ! for comparison.
82 call disp%show("weisum")
83 call disp%show( weisum )
84 call disp%show("mean")
85 call disp%show( mean )
86 call disp%skip()
87 call disp%show("call setMean(mean, sample(:,1), sample(:,2), weight, weisum)")
88 call setMean(mean, sample(:,1), sample(:,2), weight, weisum)
89 call disp%show("weisum")
90 call disp%show( weisum )
91 call disp%show("mean")
92 call disp%show( mean )
93 call disp%show("call setMean(mean, sample, 1_IK, weight, weisum) ! for comparison.")
94 call setMean(mean, sample, 1_IK, weight, weisum) ! for comparison.
95 call disp%show("weisum")
96 call disp%show( weisum )
97 call disp%show("mean")
98 call disp%show( mean )
99 call disp%skip()
100 end block
101
102 call disp%skip()
103 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
104 call disp%show("!Compute the mean of a 2-D array along a specific dimension.")
105 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
106 call disp%skip()
107
108 block
109 real(TKG), allocatable :: mean(:)
110 real(TKG), allocatable :: sample(:,:)
111 call disp%skip()
112 call disp%show("ndim = getUnifRand(1, 3); nsam = getUnifRand(1, 5)")
113 ndim = getUnifRand(1, 3); nsam = getUnifRand(1, 5)
114 call disp%show("sample = getUnifRand(0., 1., ndim, nsam)")
115 sample = getUnifRand(0., 1., ndim, nsam)
116 call disp%show("sample")
117 call disp%show( sample )
118 call disp%show("call setResized(mean, ndim)")
119 call setResized(mean, ndim)
120 call disp%show("call setMean(mean(1), sample)")
121 call setMean(mean(1), sample)
122 call disp%show("mean(1)")
123 call disp%show( mean(1) )
124 call disp%show("call setMean(mean, sample, dim = 2_IK)")
125 call setMean(mean, sample, dim = 2_IK)
126 call disp%show("mean")
127 call disp%show( mean )
128 call disp%show("call setMean(mean, transpose(sample), dim = 1_IK)")
129 call setMean(mean, transpose(sample), dim = 1_IK)
130 call disp%show("mean")
131 call disp%show( mean )
132 call disp%skip()
133 end block
134
135 call disp%skip()
136 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
137 call disp%show("!Compute the mean of a 1-D weighted array.")
138 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
139 call disp%skip()
140
141 block
142 real(TKG) :: mean
143 real(TKG), allocatable :: sample(:)
144 real(TKG), allocatable :: weight(:)
145 real(TKG) :: weisum
146 call disp%skip()
147 call disp%show("nsam = getUnifRand(1, 5)")
148 nsam = getUnifRand(1, 5)
149 call disp%show("sample = getUnifRand(0., 1., nsam)")
150 sample = getUnifRand(0., 1., nsam)
151 call disp%show("sample")
152 call disp%show( sample )
153 call disp%show("weight = getUnifRand(1, 100, nsam)")
154 weight = getUnifRand(1, 100, nsam)
155 call disp%show("weight")
156 call disp%show( weight )
157 call disp%show("call setMean(mean, sample)")
158 call setMean(mean, sample)
159 call disp%show("mean")
160 call disp%show( mean )
161 call disp%show("call setMean(mean, sample, dim = 1_IK)")
162 call setMean(mean, sample, dim = 1_IK)
163 call disp%show("mean")
164 call disp%show( mean )
165 call disp%show("call setMean(mean, sample, weight, weisum)")
166 call setMean(mean, sample, weight, weisum)
167 call disp%show("[weisum, sum(weight)]")
168 call disp%show( [weisum, sum(weight)] )
169 call disp%show("mean")
170 call disp%show( mean )
171 call disp%show("call setMean(mean, sample, 1_IK, weight, weisum)")
172 call setMean(mean, sample, 1_IK, weight, weisum)
173 call disp%show("[weisum, sum(weight)]")
174 call disp%show( [weisum, sum(weight)] )
175 call disp%show("mean")
176 call disp%show( mean )
177 call disp%skip()
178 end block
179
180 call disp%skip()
181 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
182 call disp%show("!Compute the mean of a 2-D weighted array along a specific dimension.")
183 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
184 call disp%skip()
185
186 block
187 real(TKG), allocatable :: mean(:)
188 real(TKG), allocatable :: sample(:,:)
189 real(TKG), allocatable :: weight(:)
190 real(TKG) :: weisum
191 call disp%skip()
192 call disp%show("ndim = getUnifRand(1, 3); nsam = getUnifRand(1, 5)")
193 ndim = getUnifRand(1, 3); nsam = getUnifRand(1, 5)
194 call disp%show("sample = getUnifRand(0., 1., ndim, nsam)")
195 sample = getUnifRand(0., 1., ndim, nsam)
196 call disp%show("sample")
197 call disp%show( sample )
198 call disp%show("weight = getUnifRand(1, 100, nsam)")
199 weight = getUnifRand(1, 100, nsam)
200 call disp%show("weight")
201 call disp%show( weight )
202 call disp%show("call setResized(mean, ndim)")
203 call setResized(mean, ndim)
204 call disp%show("call setMean(mean(1), sample, [(weight, idim = 1, ndim)], weisum)")
205 call setMean(mean(1), sample, [(weight, idim = 1, ndim)], weisum)
206 call disp%show("mean(1)")
207 call disp%show( mean(1) )
208 call disp%show("call setMean(mean, sample, 2_IK, weight, weisum)")
209 call setMean(mean, sample, 2_IK, weight, weisum)
210 call disp%show("[weisum, sum(weight)]")
211 call disp%show( [weisum, sum(weight)] )
212 call disp%show("mean")
213 call disp%show( mean )
214 call disp%show("call setMean(mean, transpose(sample), 1_IK, weight, weisum)")
215 call setMean(mean, transpose(sample), 1_IK, weight, weisum)
216 call disp%show("[weisum, sum(weight)]")
217 call disp%show( [weisum, sum(weight)] )
218 call disp%show("mean")
219 call disp%show( mean )
220 call disp%skip()
221 end block
222
223 call disp%skip()
224 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
225 call disp%show("!Compute the mean of a multidimensional array by associating it with a 1D pointer.")
226 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
227 call disp%skip()
228
229 block
230 real(TKG) :: mean
231 integer(IK) :: nslice
232 real(TKG), allocatable, target :: sample(:,:,:)
233 real(TKG), pointer :: samptr(:)
234 call disp%skip()
235 call disp%show("ndim = getUnifRand(1, 2); nsam = getUnifRand(1, 3); nslice = getUnifRand(1, 4)")
236 ndim = getUnifRand(1, 2); nsam = getUnifRand(1, 3); nslice = getUnifRand(1, 4)
237 call disp%show("sample = getUnifRand(0., 1., ndim, nsam, nslice)")
238 sample = getUnifRand(0., 1., ndim, nsam, nslice)
239 call disp%show("sample")
240 call disp%show( sample )
241 call disp%show("shape(sample)")
242 call disp%show( shape(sample) )
243 call disp%show("samptr(1:product(shape(sample))) => sample")
244 samptr(1:product(shape(sample))) => sample
245 call disp%show("call setMean(mean, samptr)")
246 call setMean(mean, samptr)
247 call disp%show("nullify(samptr)")
248 nullify(samptr)
249 call disp%show("mean")
250 call disp%show( mean )
251 call disp%show("call setMean(mean, reshape(sample, [product(shape(sample))]))")
252 call setMean(mean, reshape(sample, [product(shape(sample))]))
253 call disp%show("mean")
254 call disp%show( mean )
255 call disp%skip()
256 end block
257
258end program example
Allocate or resize (shrink or expand) an input allocatable scalar string or array of rank 1....
Generate and return a scalar or a contiguous array of rank 1 of length s1 of randomly uniformly distr...
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11726
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11508
This module contains procedures and generic interfaces for resizing allocatable arrays of various typ...
This module contains classes and procedures for computing various statistical quantities related to t...
This module contains classes and procedures for input/output (IO) or generic display operations on st...
Definition: pm_io.F90:252
type(display_type) disp
This is a scalar module variable an object of type display_type for general display.
Definition: pm_io.F90:11393
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
Definition: pm_kind.F90:540
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
Definition: pm_kind.F90:539
integer, parameter RKS
The single-precision real kind in Fortran mode. On most platforms, this is an 32-bit real kind.
Definition: pm_kind.F90:567
Generate and return an object of type display_type.
Definition: pm_io.F90:10282

Example Unix compile command via Intel ifort compiler
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example output
1
2!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3!Compute the mean of a 1-D array.
4!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5
6
7nsam = getUnifRand(1, 5)
8sample = getUnifRand(0., 1., nsam)
9sample
10+0.842041314, +0.839977026
11call setMean(mean, sample)
12mean
13+0.841009140
14call setMean(mean, sample, dim = 1_IK)
15mean
16+0.841009140
17
18
19!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
20!Compute the mean of a pair of data series.
21!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22
23
24nsam = getUnifRand(1, 5)
25sample = getUnifRand(0., 1., nsam, 2_IK)
26sample
27+0.765661240, +0.417786241E-1
28+0.293720424, +0.342870712
29+0.694622219, +0.966399550
30+0.335063875, +0.886310935
31+0.640322268, +0.760899365
32call setMean(mean, sample(:,1), sample(:,2))
33mean
34+0.545878053, +0.599651873
35call setMean(mean, sample, dim = 1_IK) ! for comparison.
36mean
37+0.545878053, +0.599651873
38
39weight = getUnifRand(1, 100, nsam)
40weight
41+18.0000000, +62.0000000, +48.0000000, +86.0000000, +85.0000000
42call setMean(mean, sample(:,1), sample(:,2), weight, weisum)
43weisum
44+299.000000
45mean
46+0.496914089, +0.699987829
47call setMean(mean, sample, 1_IK, weight, weisum) ! for comparison.
48weisum
49+299.000000
50mean
51+0.496914089, +0.699987829
52
53call setMean(mean, sample(:,1), sample(:,2), weight, weisum)
54weisum
55+299.000000
56mean
57+0.496914089, +0.699987829
58call setMean(mean, sample, 1_IK, weight, weisum) ! for comparison.
59weisum
60+299.000000
61mean
62+0.496914089, +0.699987829
63
64
65!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66!Compute the mean of a 2-D array along a specific dimension.
67!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
68
69
70ndim = getUnifRand(1, 3); nsam = getUnifRand(1, 5)
71sample = getUnifRand(0., 1., ndim, nsam)
72sample
73+0.214186609, +0.162420273E-1, +0.383274555, +0.664523840
74call setResized(mean, ndim)
75call setMean(mean(1), sample)
76mean(1)
77+0.319556773
78call setMean(mean, sample, dim = 2_IK)
79mean
80+0.319556773
81call setMean(mean, transpose(sample), dim = 1_IK)
82mean
83+0.319556773
84
85
86!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
87!Compute the mean of a 1-D weighted array.
88!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
89
90
91nsam = getUnifRand(1, 5)
92sample = getUnifRand(0., 1., nsam)
93sample
94+0.894065917, +0.508470118, +0.556979775, +0.781553388, +0.955277503
95weight = getUnifRand(1, 100, nsam)
96weight
97+44.0000000, +13.0000000, +30.0000000, +1.00000000, +91.0000000
98call setMean(mean, sample)
99mean
100+0.739269376
101call setMean(mean, sample, dim = 1_IK)
102mean
103+0.739269376
104call setMean(mean, sample, weight, weisum)
105[weisum, sum(weight)]
106+179.000000, +179.000000
107mean
108+0.840057015
109call setMean(mean, sample, 1_IK, weight, weisum)
110[weisum, sum(weight)]
111+179.000000, +179.000000
112mean
113+0.840057015
114
115
116!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117!Compute the mean of a 2-D weighted array along a specific dimension.
118!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119
120
121ndim = getUnifRand(1, 3); nsam = getUnifRand(1, 5)
122sample = getUnifRand(0., 1., ndim, nsam)
123sample
124+0.540724993
125weight = getUnifRand(1, 100, nsam)
126weight
127+95.0000000
128call setResized(mean, ndim)
129call setMean(mean(1), sample, [(weight, idim = 1, ndim)], weisum)
130mean(1)
131+0.540724993
132call setMean(mean, sample, 2_IK, weight, weisum)
133[weisum, sum(weight)]
134+95.0000000, +95.0000000
135mean
136+0.540724993
137call setMean(mean, transpose(sample), 1_IK, weight, weisum)
138[weisum, sum(weight)]
139+95.0000000, +95.0000000
140mean
141+0.540724993
142
143
144!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145!Compute the mean of a multidimensional array by associating it with a 1D pointer.
146!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147
148
149ndim = getUnifRand(1, 2); nsam = getUnifRand(1, 3); nslice = getUnifRand(1, 4)
150sample = getUnifRand(0., 1., ndim, nsam, nslice)
151sample
152slice(:,:,1) =
153+0.342382669
154+0.412743688E-1
155slice(:,:,2) =
156+0.914745092
157+0.980952859
158shape(sample)
159+2, +1, +2
160samptr(1:product(shape(sample))) => sample
161call setMean(mean, samptr)
162nullify(samptr)
163mean
164+0.569838762
165call setMean(mean, reshape(sample, [product(shape(sample))]))
166mean
167+0.569838762
168
169
Benchmarks:


Benchmark :: The runtime performance of setMean along different sample dimensions.

1! Test the performance of Cholesky factorization computation using an assumed-shape interface vs. explicit-shape interface.
2program benchmark
3
4 use pm_kind, only: IK, LK, RKG => RKD, SK
5 use pm_bench, only: bench_type
6
7 implicit none
8
9 integer(IK) :: itry, ntry
10 integer(IK) :: i
11 integer(IK) :: iarr
12 integer(IK) :: fileUnit
13 integer(IK) , parameter :: NARR = 18_IK
14 real(RKG) , allocatable :: samdim1(:,:)
15 real(RKG) , allocatable :: samdim2(:,:)
16 type(bench_type), allocatable :: bench(:)
17 integer(IK) , parameter :: nsammax = 2**NARR
18 integer(IK) , parameter :: ndim = 5_IK
19 real(RKG) :: mean(ndim)
20 integer(IK) :: isam, nsam
21 real(RKG) :: dumm
22
23 bench = [ bench_type(name = SK_"intrinsicDIM1", exec = intrinsicDIM1, overhead = setOverhead) &
24 , bench_type(name = SK_"intrinsicDIM2", exec = intrinsicDIM2, overhead = setOverhead) &
25 , bench_type(name = SK_"setMeanDIM1", exec = setMeanDIM1, overhead = setOverhead) &
26 , bench_type(name = SK_"setMeanDIM2", exec = setMeanDIM2, overhead = setOverhead) &
27 ]
28
29 write(*,"(*(g0,:,' '))")
30 write(*,"(*(g0,:,' '))") "sample covariance benchmarking..."
31 write(*,"(*(g0,:,' '))")
32
33 open(newunit = fileUnit, file = "main.out", status = "replace")
34
35 write(fileUnit, "(*(g0,:,','))") "nsam", (bench(i)%name, i = 1, size(bench))
36
37 dumm = 0._RKG
38 loopOverMatrixSize: do iarr = 1, NARR - 1
39
40 nsam = 2**iarr
41 ntry = nsammax / nsam
42 allocate(samdim1(nsam, ndim), samdim2(ndim, nsam))
43 write(*,"(*(g0,:,' '))") "Benchmarking setMeanDIM1() vs. setMeanDIM2()", nsam, ntry
44
45 do i = 1, size(bench)
46 bench(i)%timing = bench(i)%getTiming()
47 end do
48
49 write(fileUnit,"(*(g0,:,','))") nsam, (bench(i)%timing%mean / ntry, i = 1, size(bench))
50 deallocate(samdim1, samdim2)
51
52 end do loopOverMatrixSize
53 write(*,"(*(g0,:,' '))") dumm
54
55 close(fileUnit)
56
57contains
58
59 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
60 ! procedure wrappers.
61 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62
63 subroutine setOverhead()
64 do itry = 1, ntry
65 call random_number(mean)
66 dumm = dumm - mean(1)
67 end do
68 end subroutine
69
70 subroutine setSamDIM1()
71 do isam = 1, nsam
72 call random_number(samdim1(isam, 1 : ndim))
73 end do
74 end subroutine
75
76 subroutine setSamDIM2()
77 do isam = 1, nsam
78 call random_number(samdim2(1 : ndim, isam))
79 end do
80 end subroutine
81
82 subroutine setMeanDIM1()
83 use pm_sampleMean, only: setMean
84 do itry = 1, ntry
85 call setSamDIM1()
86 call setMean(mean, samdim1, dim = 1_IK)
87 dumm = dumm + mean(1)
88 end do
89 end subroutine
90
91 subroutine setMeanDIM2()
92 use pm_sampleMean, only: setMean
93 do itry = 1, ntry
94 call setSamDIM2()
95 call setMean(mean, samdim2, dim = 2_IK)
96 dumm = dumm + mean(1)
97 end do
98 end subroutine
99
100 subroutine intrinsicDIM1()
101 do itry = 1, ntry
102 call setSamDIM1()
103 mean = sum(samdim1, dim = 1) / size(samdim1, dim = 1)
104 dumm = dumm + mean(1)
105 end do
106 end subroutine
107
108 subroutine intrinsicDIM2()
109 do itry = 1, ntry
110 call setSamDIM2()
111 mean = sum(samdim2, dim = 2) / size(samdim1, dim = 2)
112 dumm = dumm + mean(1)
113 end do
114 end subroutine
115
116end program benchmark
Generate and return an object of type timing_type containing the benchmark timing information and sta...
Definition: pm_bench.F90:574
This module contains abstract interfaces and types that facilitate benchmarking of different procedur...
Definition: pm_bench.F90:41
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541
integer, parameter RKD
The double precision real kind in Fortran mode. On most platforms, this is an 64-bit real kind.
Definition: pm_kind.F90:568
This is the class for creating benchmark and performance-profiling objects.
Definition: pm_bench.F90:386
subroutine bench(sort, arraySize)

Example Unix compile command via Intel ifort compiler
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Postprocessing of the benchmark output
1#!/usr/bin/env python
2
3import matplotlib.pyplot as plt
4import pandas as pd
5import numpy as np
6
7import os
8dirname = os.path.basename(os.getcwd())
9
10fontsize = 14
11
12df = pd.read_csv("main.out", delimiter = ",")
13colnames = list(df.columns.values)
14
15
18
19ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
20ax = plt.subplot()
21
22for colname in colnames[1:]:
23 plt.plot( df[colnames[0]].values
24 , df[colname].values
25 , linewidth = 2
26 )
27
28plt.xticks(fontsize = fontsize)
29plt.yticks(fontsize = fontsize)
30ax.set_xlabel(colnames[0], fontsize = fontsize)
31ax.set_ylabel("Runtime [ seconds ]", fontsize = fontsize)
32ax.set_title(" vs. ".join(colnames[1:])+"\nLower is better.", fontsize = fontsize)
33ax.set_xscale("log")
34ax.set_yscale("log")
35plt.minorticks_on()
36plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
37ax.tick_params(axis = "y", which = "minor")
38ax.tick_params(axis = "x", which = "minor")
39ax.legend ( colnames[1:]
40 #, loc='center left'
41 #, bbox_to_anchor=(1, 0.5)
42 , fontsize = fontsize
43 )
44
45plt.tight_layout()
46plt.savefig("benchmark." + dirname + ".runtime.png")
47
48
51
52ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
53ax = plt.subplot()
54
55plt.plot( df[colnames[0]].values
56 , np.ones(len(df[colnames[0]].values))
57 , linestyle = "--"
58 #, color = "black"
59 , linewidth = 2
60 )
61for colname in colnames[2:]:
62 plt.plot( df[colnames[0]].values
63 , df[colname].values / df[colnames[1]].values
64 , linewidth = 2
65 )
66
67plt.xticks(fontsize = fontsize)
68plt.yticks(fontsize = fontsize)
69ax.set_xlabel(colnames[0], fontsize = fontsize)
70ax.set_ylabel("Runtime compared to {}".format(colnames[1]), fontsize = fontsize)
71ax.set_title("Runtime Ratio Comparison. Lower means faster.\nLower than 1 means faster than {}().".format(colnames[1]), fontsize = fontsize)
72ax.set_xscale("log")
73ax.set_yscale("log")
74plt.minorticks_on()
75plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
76ax.tick_params(axis = "y", which = "minor")
77ax.tick_params(axis = "x", which = "minor")
78ax.legend ( colnames[1:]
79 #, bbox_to_anchor = (1, 0.5)
80 #, loc = "center left"
81 , fontsize = fontsize
82 )
83
84plt.tight_layout()
85plt.savefig("benchmark." + dirname + ".runtime.ratio.png")

Visualization of the benchmark output

Benchmark moral
  1. The procedures under the generic interface setMean can compute the covariance under two different sample axes.
  2. Recall that C is a row-major language while Fortran is a column-major language.
  3. As such, one would expect the computations for a sample whose observations are stored along the second axis would be faster in the Fortran programming language.
  4. However, such an expectation does not appear to hold at all times and appears to depend significantly on the computing architecture and the number of data attributes involved.
  5. The higher the number of data attributes, the more likely the computations along the second axis of sample will be faster.
  6. Note that for small number of data attributes, the computations along the second data axis involve a small loop that has significant computational cost due to the implicit branching involved in the loop.
Test:
test_pm_sampleMean


Final Remarks


If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

  1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
  2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Author:
Amir Shahmoradi, April 21, 2017, 1:54 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin

Definition at line 2008 of file pm_sampleMean.F90.


The documentation for this interface was generated from the following file: