Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 : !!****m* ABINIT/m_Stat
6 : !! NAME
7 : !! m_Stat
8 : !!
9 : !! FUNCTION
10 : !! FIXME: add description.
11 : !!
12 : !! COPYRIGHT
13 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
14 : !! This file is distributed under the terms of the
15 : !! GNU General Public License, see ~abinit/COPYING
16 : !! or http://www.gnu.org/copyleft/gpl.txt .
17 : !!
18 : !! NOTES
19 : !!
20 : !! SOURCE
21 :
22 : #include "defs.h"
23 : MODULE m_Stat
24 : USE m_global
25 : IMPLICIT NONE
26 :
27 : PRIVATE
28 :
29 : PUBLIC :: Stat_average
30 : PUBLIC :: Stat_variance
31 : PUBLIC :: Stat_coVariance
32 : PUBLIC :: Stat_deviation
33 : PUBLIC :: Stat_linearReg
34 : PUBLIC :: Stat_powerReg
35 :
36 : CONTAINS
37 : !!***
38 :
39 272 : DOUBLE PRECISION FUNCTION Stat_average(tab)
40 :
41 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: tab
42 : !INTEGER :: sizet
43 :
44 : !sizet = SIZE(tab)
45 272 : Stat_average = 0.d0
46 : !IF ( sizet .GT. 0 ) &
47 344064 : Stat_average = SUM(tab)/DBLE(SIZE(tab))
48 272 : END FUNCTION Stat_average
49 : !!***
50 :
51 136 : DOUBLE PRECISION FUNCTION Stat_variance(tab)
52 :
53 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: tab
54 : INTEGER :: sizet
55 : INTEGER :: i
56 : DOUBLE PRECISION :: average
57 : DOUBLE PRECISION :: tab2
58 :
59 136 : sizet = SIZE(tab)
60 136 : Stat_variance = 0.d0
61 :
62 : !IF ( sizet .GT. 0 ) THEN
63 136 : average = Stat_average(tab)
64 343656 : DO i = 1, sizet
65 343520 : tab2 = tab(i)-average
66 343520 : tab2 = tab2 * tab2
67 343656 : Stat_variance = Stat_variance + tab2
68 : END DO
69 136 : Stat_variance = Stat_variance / DBLE(sizet)
70 : !ELSE
71 : ! Stat_variance = 0.d-16
72 : !END IF
73 136 : END FUNCTION Stat_variance
74 : !!***
75 :
76 34 : DOUBLE PRECISION FUNCTION Stat_coVariance(tab1, tab2)
77 :
78 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: tab1
79 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: tab2
80 : INTEGER :: size1
81 : INTEGER :: size2
82 : INTEGER :: i
83 : DOUBLE PRECISION :: average1
84 : DOUBLE PRECISION :: average2
85 : DOUBLE PRECISION :: tmp
86 :
87 34 : size1 = SIZE(tab1)
88 34 : size2 = SIZE(tab2)
89 :
90 34 : IF ( size1 .NE. size2 ) &
91 0 : CALL ERROR("Stat_coVariance : Array sizes mismatch ")
92 :
93 34 : average1 = Stat_average(tab1)
94 34 : average2 = Stat_average(tab2)
95 34 : Stat_coVariance = 0.d0
96 :
97 : !IF ( size1 .GT. 0 ) THEN
98 102 : DO i = 1, size1
99 68 : tmp = (tab1(i)-average1)*(tab2(i)-average2)
100 102 : Stat_coVariance = Stat_coVariance + tmp
101 : END DO
102 34 : Stat_coVariance = Stat_coVariance / DBLE(size1)
103 : !ELSE
104 : ! Stat_coVariance = 0.d-16
105 : !END IF
106 34 : END FUNCTION Stat_coVariance
107 : !!***
108 :
109 102 : DOUBLE PRECISION FUNCTION Stat_deviation(tab1)
110 :
111 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: tab1
112 :
113 102 : Stat_deviation = SQRT(Stat_variance(tab1))
114 102 : END FUNCTION Stat_deviation
115 : !!***
116 :
117 : !!****f* ABINIT/m_Stat/Stat_linearReg
118 : !! NAME
119 : !! Stat_linearReg
120 : !!
121 : !! FUNCTION
122 : !! FIXME: add description.
123 : !!
124 : !! COPYRIGHT
125 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
126 : !! This file is distributed under the terms of the
127 : !! GNU General Public License, see ~abinit/COPYING
128 : !! or http://www.gnu.org/copyleft/gpl.txt .
129 : !!
130 : !! INPUTS
131 : !! argin(sizein)=description
132 : !!
133 : !! OUTPUT
134 : !! argout(sizeout)=description
135 : !!
136 : !! SIDE EFFECTS
137 : !!
138 : !! NOTES
139 : !!
140 : !! SOURCE
141 :
142 34 : SUBROUTINE Stat_linearReg(tabX, tabY, a, b, R)
143 : !Arguments ------------------------------------
144 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN ) :: tabX
145 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN ) :: tabY
146 : DOUBLE PRECISION , INTENT(OUT) :: a
147 : DOUBLE PRECISION , INTENT(OUT) :: b
148 : DOUBLE PRECISION , INTENT(OUT) :: R
149 : DOUBLE PRECISION :: coVar
150 : DOUBLE PRECISION :: Var
151 :
152 34 : coVar = Stat_coVariance(tabX, tabY)
153 34 : Var = Stat_variance(tabX)
154 34 : a = coVar / var
155 34 : b = Stat_average(tabY) - a* Stat_average(tabX)
156 34 : R = ABS(a * SQRT(var) / Stat_deviation(tabY))
157 34 : END SUBROUTINE Stat_linearReg
158 : !!***
159 :
160 : !!****f* ABINIT/m_Stat/Stat_powerReg
161 : !! NAME
162 : !! Stat_powerReg
163 : !!
164 : !! FUNCTION
165 : !! FIXME: add description.
166 : !!
167 : !! COPYRIGHT
168 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
169 : !! This file is distributed under the terms of the
170 : !! GNU General Public License, see ~abinit/COPYING
171 : !! or http://www.gnu.org/copyleft/gpl.txt .
172 : !!
173 : !! INPUTS
174 : !! argin(sizein)=description
175 : !!
176 : !! OUTPUT
177 : !! argout(sizeout)=description
178 : !!
179 : !! SIDE EFFECTS
180 : !!
181 : !! NOTES
182 : !!
183 : !! SOURCE
184 :
185 34 : SUBROUTINE Stat_powerReg(tabX, tabY, a, b, R)
186 : !Arguments ------------------------------------
187 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN ) :: tabX
188 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN ) :: tabY
189 : DOUBLE PRECISION , INTENT(OUT) :: a
190 : DOUBLE PRECISION , INTENT(OUT) :: b
191 : DOUBLE PRECISION , INTENT(OUT) :: R
192 : INTEGER :: size1
193 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: tab1
194 : DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: tab2
195 : INTEGER :: i
196 :
197 34 : size1 = SIZE(tabX)
198 34 : IF ( size1 .NE. SIZE(tabY) ) &
199 0 : CALL ERROR("Stat_powerReg : Array sizes mismatch ")
200 :
201 : FREEIF(tab1)
202 : FREEIF(tab2)
203 102 : MALLOC(tab1,(1:size1))
204 68 : MALLOC(tab2,(1:size1))
205 :
206 :
207 : ! IF ( ISNAN(a) .OR. ISNAN(b) ) THEN
208 102 : DO i = 1, size1
209 68 : tab1(i) = LOG(tabX(i))
210 102 : IF ( tabY(i) .NE. 0.d0 ) THEN
211 68 : tab2(i) = LOG(tabY(i))
212 : ELSE
213 0 : tab2(i) = LOG(1.d-16)
214 : END IF
215 : !WRITE(92,'(4E22.14)') tabX(i), tab1(i), tabY(i), tab2(i)
216 : END DO
217 : !WRITE(92,*)
218 : !END IF
219 34 : CALL Stat_linearReg(tab1, tab2, b, a, R)
220 34 : a = EXP(a)
221 34 : FREE(tab1)
222 34 : FREE(tab2)
223 34 : END SUBROUTINE Stat_powerReg
224 : !!***
225 :
226 : DOUBLE PRECISION FUNCTION Stat_simpson(func, a, b, N)
227 :
228 : DOUBLE PRECISION, DIMENSION(:), POINTER :: func !vz_i
229 : DOUBLE PRECISION , INTENT(IN) :: a
230 : DOUBLE PRECISION , INTENT(IN) :: b
231 : INTEGER , INTENT(IN) :: N
232 : INTEGER :: i
233 : INTEGER :: x0
234 : INTEGER :: x1
235 : INTEGER :: x2
236 : DOUBLE PRECISION :: dtau
237 : DOUBLE PRECISION :: J
238 :
239 : dtau = (b-a)/DBLE(N-1)
240 : J=0.d0
241 : x2 = 1
242 : x1 = 0
243 : DO i = 1, N-2, 2
244 : x0 = x2
245 : x1 = x1 + 2
246 : x2 = x0 + 2
247 : J = J + func(x0) + 4.d0*func(x1) + func(x2)
248 : END DO
249 : J = J * dtau / 3.d0
250 : Stat_simpson = J
251 : END FUNCTION Stat_simpson
252 : !!***
253 :
254 : END MODULE m_Stat
|
