Line data Source code
1 : !> \file mo_cost.f90
2 : !> \brief \copybrief mo_cost
3 : !> \details \copydetails mo_cost
4 :
5 : !> \brief Added for testing purposes of test_mo_anneal
6 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
7 : !! FORCES is released under the LGPLv3+ license \license_note
8 : Module mo_cost
9 :
10 : use mo_kind, only: sp, dp
11 :
12 : IMPLICIT NONE
13 :
14 : PRIVATE
15 :
16 : PUBLIC :: cost_sp, cost_dp
17 : PUBLIC :: cost_valid_sp, cost_valid_dp
18 : PUBLIC :: range_sp, range_dp
19 : public :: cost_objective
20 :
21 : CONTAINS
22 :
23 : !> \brief function: `f(x) = ax^3 + bx^2 + cx + d`
24 0 : FUNCTION cost_sp (paraset)
25 :
26 : implicit none
27 :
28 : REAL(SP), DIMENSION(:), INTENT(IN) :: paraset
29 : REAL(SP) :: cost_sp
30 0 : REAL(SP), DIMENSION(6,2) :: meas
31 0 : REAL(SP), DIMENSION(6) :: calc
32 :
33 : ! function: f(x) = ax^3 + bx^2 + cx + d
34 : ! measurements: (0.5,5.725), (1.0, 21.7), (1.5, 49.175), (2.0, 88.9), (2.5, 141.625), (3.0, 208.1)
35 : ! --> a=1.0, b=20.0, c=0.2, d=0.5
36 :
37 0 : meas(:,1) = (/0.5_sp, 1.0_sp, 1.5_sp, 2.0_sp, 2.5_sp, 3.0_sp/)
38 0 : meas(:,2) = (/5.7250_sp, 21.7000_sp, 49.1750_sp, 88.9000_sp, 141.6250_sp, 208.1000_sp/)
39 :
40 0 : calc(:) = paraset(1)*meas(:,1)**3+paraset(2)*meas(:,1)**2+paraset(3)*meas(:,1)+paraset(4)
41 :
42 : ! MAE Mean Absolute Error
43 0 : cost_sp = sum(abs( meas(:,2)-calc(:) ))/size(meas,1)
44 :
45 : RETURN
46 0 : END FUNCTION cost_sp
47 :
48 : !> \brief function: `f(x) = ax^3 + bx^2 + cx + d`
49 0 : FUNCTION cost_dp (paraset)
50 :
51 : implicit none
52 :
53 : REAL(DP), DIMENSION(:), INTENT(IN) :: paraset
54 : REAL(DP) :: cost_dp
55 0 : REAL(DP), DIMENSION(6,2) :: meas
56 0 : REAL(DP), DIMENSION(6) :: calc
57 :
58 : ! function: f(x) = ax^3 + bx^2 + cx + d
59 : ! measurements: (0.5,5.725), (1.0, 21.7), (1.5, 49.175), (2.0, 88.9), (2.5, 141.625), (3.0, 208.1)
60 : ! --> a=1.0, b=20.0, c=0.2, d=0.5
61 :
62 0 : meas(:,1) = (/0.5_dp, 1.0_dp, 1.5_dp, 2.0_dp, 2.5_dp, 3.0_dp/)
63 0 : meas(:,2) = (/5.7250_dp, 21.7000_dp, 49.1750_dp, 88.9000_dp, 141.6250_dp, 208.1000_dp/)
64 :
65 0 : calc(:) = paraset(1)*meas(:,1)**3+paraset(2)*meas(:,1)**2+paraset(3)*meas(:,1)+paraset(4)
66 :
67 : ! MAE Mean Absolute Error
68 0 : cost_dp = sum(abs( meas(:,2)-calc(:) ))/size(meas,1)
69 :
70 : RETURN
71 0 : END FUNCTION cost_dp
72 :
73 : !> \brief function: `f(x) = ax^3 + bx^2 + cx + d`
74 0 : FUNCTION cost_valid_sp (paraset,status_in)
75 :
76 : implicit none
77 :
78 : REAL(SP), DIMENSION(:), INTENT(IN) :: paraset
79 : LOGICAL, OPTIONAL, INTENT(OUT) :: status_in
80 : REAL(SP) :: cost_valid_sp
81 0 : REAL(SP), DIMENSION(6,2) :: meas
82 0 : REAL(SP), DIMENSION(6) :: calc
83 :
84 : ! function: f(x) = ax^3 + bx^2 + cx + d
85 : ! measurements: (0.5,5.725), (1.0, 21.7), (1.5, 49.175), (2.0, 88.9), (2.5, 141.625), (3.0, 208.1)
86 : ! --> a=1.0, b=20.0, c=0.2, d=0.5
87 :
88 0 : meas(:,1) = (/0.5_sp, 1.0_sp, 1.5_sp, 2.0_sp, 2.5_sp, 3.0_sp/)
89 0 : meas(:,2) = (/5.7250_sp, 21.7000_sp, 49.1750_sp, 88.9000_sp, 141.6250_sp, 208.1000_sp/)
90 :
91 0 : calc(:) = paraset(1)*meas(:,1)**3+paraset(2)*meas(:,1)**2+paraset(3)*meas(:,1)+paraset(4)
92 :
93 0 : if (present(status_in)) then
94 0 : status_in = .true.
95 : ! Define a status .false. if calculation of "calc" was not successful
96 : end if
97 :
98 : ! MAE Mean Absolute Error
99 0 : cost_valid_sp = sum(abs( meas(:,2)-calc(:) ))/size(meas,1)
100 :
101 : RETURN
102 0 : END FUNCTION cost_valid_sp
103 :
104 : !> \brief function: `f(x) = ax^3 + bx^2 + cx + d`
105 0 : FUNCTION cost_valid_dp (paraset,status_in)
106 :
107 : implicit none
108 :
109 : REAL(DP), DIMENSION(:), INTENT(IN) :: paraset
110 : LOGICAL, OPTIONAL, INTENT(OUT) :: status_in
111 : REAL(DP) :: cost_valid_dp
112 0 : REAL(DP), DIMENSION(6,2) :: meas
113 0 : REAL(DP), DIMENSION(6) :: calc
114 :
115 : ! function: f(x) = ax^3 + bx^2 + cx + d
116 : ! measurements: (0.5,5.725), (1.0, 21.7), (1.5, 49.175), (2.0, 88.9), (2.5, 141.625), (3.0, 208.1)
117 : ! --> a=1.0, b=20.0, c=0.2, d=0.5
118 :
119 0 : meas(:,1) = (/0.5_dp, 1.0_dp, 1.5_dp, 2.0_dp, 2.5_dp, 3.0_dp/)
120 0 : meas(:,2) = (/5.7250_dp, 21.7000_dp, 49.1750_dp, 88.9000_dp, 141.6250_dp, 208.1000_dp/)
121 :
122 0 : calc(:) = paraset(1)*meas(:,1)**3+paraset(2)*meas(:,1)**2+paraset(3)*meas(:,1)+paraset(4)
123 :
124 0 : if (present(status_in)) then
125 0 : status_in = .true.
126 : ! Define a status .false. if calculation of "calc" was not successful
127 : end if
128 :
129 : ! MAE Mean Absolute Error
130 0 : cost_valid_dp = sum(abs( meas(:,2)-calc(:) ))/size(meas,1)
131 :
132 : RETURN
133 0 : END FUNCTION cost_valid_dp
134 :
135 : !> \brief dummy range
136 1287500 : SUBROUTINE range_dp(paraset, iPar, rangePar)
137 0 : use mo_kind
138 : REAL(DP), DIMENSION(:), INTENT(IN) :: paraset
139 : INTEGER(I4), INTENT(IN) :: iPar
140 : REAL(DP), DIMENSION(2), INTENT(OUT) :: rangePar
141 :
142 : ! Range does not depend on parameter set
143 : ! select case(iPar)
144 : ! case(1_i4)
145 : ! rangePar(1) = 0.0_dp
146 : ! rangePar(2) = 10.0_dp
147 : ! case(2_i4)
148 : ! rangePar(1) = 0.0_dp
149 : ! rangePar(2) = 40.0_dp
150 : ! case(3_i4)
151 : ! rangePar(1) = 0.0_dp
152 : ! rangePar(2) = 10.0_dp
153 : ! case(4_i4)
154 : ! rangePar(1) = 0.0_dp
155 : ! rangePar(2) = 5.0_dp
156 : ! end select
157 :
158 : ! Range of parameter 2 depends on value of parameter 1:
159 : ! parameter 2 at most 40* parameter 1 :
160 : ! 0 <= p2 <= 40p1
161 : ! 0 <= p1 <= 0.025p2
162 1610002 : select case(iPar)
163 : case(1_i4)
164 322502 : rangePar(1) = 0.025_dp*paraset(2)
165 322502 : rangePar(2) = 10.0_dp
166 : case(2_i4)
167 321474 : rangePar(1) = 0.0_dp
168 321474 : rangePar(2) = 40.0_dp*paraset(1)
169 : case(3_i4)
170 321610 : rangePar(1) = 0.0_dp
171 321610 : rangePar(2) = 10.0_dp
172 : case(4_i4)
173 321914 : rangePar(1) = 0.0_dp
174 1287500 : rangePar(2) = 5.0_dp
175 : end select
176 :
177 1287500 : END SUBROUTINE range_dp
178 :
179 : !> \brief dummy range
180 0 : SUBROUTINE range_sp(paraset, iPar, rangePar)
181 1287500 : use mo_kind
182 : REAL(SP), DIMENSION(:), INTENT(IN) :: paraset
183 : INTEGER(I4), INTENT(IN) :: iPar
184 : REAL(SP), DIMENSION(2), INTENT(OUT) :: rangePar
185 :
186 : ! Range does not depend on parameter set
187 : ! select case(iPar)
188 : ! case(1_i4)
189 : ! rangePar(1) = 0.0_sp
190 : ! rangePar(2) = 10.0_sp
191 : ! case(2_i4)
192 : ! rangePar(1) = 0.0_sp
193 : ! rangePar(2) = 40.0_sp
194 : ! case(3_i4)
195 : ! rangePar(1) = 0.0_sp
196 : ! rangePar(2) = 10.0_sp
197 : ! case(4_i4)
198 : ! rangePar(1) = 0.0_sp
199 : ! rangePar(2) = 5.0_sp
200 : ! end select
201 :
202 : ! Range of parameter 2 depends on value of parameter 1:
203 : ! parameter 2 at most 4* parameter 1 :
204 : ! 0 <= p2 <= 4p1
205 : ! 0.25p2 <= p1 <= 10.0
206 0 : select case(iPar)
207 : case(1_i4)
208 0 : rangePar(1) = 0.025_sp*paraset(2)
209 0 : rangePar(2) = 10.0_sp
210 : case(2_i4)
211 0 : rangePar(1) = 0.0_sp
212 0 : rangePar(2) = 40.0_sp*paraset(1)
213 : case(3_i4)
214 0 : rangePar(1) = 0.0_sp
215 0 : rangePar(2) = 10.0_sp
216 : case(4_i4)
217 0 : rangePar(1) = 0.0_sp
218 0 : rangePar(2) = 5.0_sp
219 : end select
220 :
221 0 : END SUBROUTINE range_sp
222 :
223 : !> \brief dummy cost objective function
224 1287511 : FUNCTION cost_objective(parameterset, eval, arg1, arg2, arg3)
225 :
226 0 : use mo_kind, only: dp
227 : use mo_optimization_utils, only: eval_interface
228 : use mo_optimization_types, only : optidata_sim
229 :
230 : implicit none
231 :
232 : real(dp), intent(in), dimension(:) :: parameterset
233 : procedure(eval_interface), INTENT(IN), pointer :: eval
234 : real(dp), optional, intent(in) :: arg1
235 : real(dp), optional, intent(out) :: arg2
236 : real(dp), optional, intent(out) :: arg3
237 : real(dp) :: cost_objective
238 :
239 1287511 : type(optidata_sim), dimension(:), allocatable :: et_opti
240 19312665 : REAL(DP), DIMENSION(6,2) :: meas
241 9012577 : REAL(DP), DIMENSION(6) :: calc
242 :
243 1287511 : call eval(parameterset, etOptiSim=et_opti)
244 :
245 : ! function: f(x) = ax^3 + bx^2 + cx + d
246 : ! measurements: (0.5,5.725), (1.0, 21.7), (1.5, 49.175), (2.0, 88.9), (2.5, 141.625), (3.0, 208.1)
247 : ! --> a=1.0, b=20.0, c=0.2, d=0.5
248 :
249 9012577 : meas(:,1) = (/0.5_dp, 1.0_dp, 1.5_dp, 2.0_dp, 2.5_dp, 3.0_dp/)
250 9012577 : meas(:,2) = (/5.7250_dp, 21.7000_dp, 49.1750_dp, 88.9000_dp, 141.6250_dp, 208.1000_dp/)
251 :
252 9012577 : calc(:) = parameterset(1)*meas(:,1)**3+parameterset(2)*meas(:,1)**2+parameterset(3)*meas(:,1)+parameterset(4)
253 :
254 : ! MAE Mean Absolute Error
255 9012577 : cost_objective = sum(abs( meas(:,2)-calc(:) ))/size(meas,1)
256 :
257 : RETURN
258 1287511 : END FUNCTION cost_objective
259 :
260 :
261 : END MODULE mo_cost
|
