Line data Source code
1 : /**
2 : * WARNING: Do Not Remove This Section
3 : *
4 : * $LastChangedRevision: 415 $
5 : * $LastChangedDate: 2010-12-18 17:21:05 +0100 (Sa, 18. Dez 2010) $
6 : * $LastChangedBy: davidjansen $
7 : *
8 : * Modified by Benoît Barbot on 25/01/2014 to integrate it into Cosmos
9 : *
10 : * MRMC is a model checker for discrete-time and continuous-time Markov
11 : * reward models. It supports reward extensions of PCTL and CSL (PRCTL
12 : * and CSRL), and allows for the automated verification of properties
13 : * concerning long-run and instantaneous rewards as well as cumulative
14 : * rewards.
15 : *
16 : * Copyright (C) The University of Twente, 2004-2008.
17 : * Copyright (C) RWTH Aachen, 2008-2009.
18 : * Authors: Ivan Zapreev
19 : *
20 : * This program is free software; you can redistribute it and/or
21 : * modify it under the terms of the GNU General Public License
22 : * as published by the Free Software Foundation; either version 2
23 : * of the License, or (at your option) any later version.
24 : *
25 : * This program is distributed in the hope that it will be useful,
26 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
27 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28 : * GNU General Public License for more details.
29 : *
30 : * You should have received a copy of the GNU General Public License
31 : * along with this program; if not, write to the Free Software
32 : * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
33 : * MA 02110-1301, USA.
34 : *
35 : * Main contact:
36 : * Lehrstuhl für Informatik 2, RWTH Aachen University
37 : * Ahornstrasse 55, 52074 Aachen, Germany
38 : * E-mail: info@mrmc-tool.org
39 : *
40 : * Old contact:
41 : * Formal Methods and Tools Group, University of Twente,
42 : * P.O. Box 217, 7500 AE Enschede, The Netherlands,
43 : * Phone: +31 53 4893767, Fax: +31 53 4893247,
44 : * E-mail: mrmc@cs.utwente.nl
45 : *
46 : * Source description: Calculate left/right truncation pts. and
47 : * weights based on Fox&Glynn algorithm.
48 : * Uses: Definition of foxglynn - foxglynn.h
49 : * Reference: "Computing Poisson Probabilities", Bennett L. Fox and
50 : * Peter W. Glynn,
51 : * Communications of the ACM , vol. 31, pp. 440-445
52 : * (1988 ).
53 : */
54 :
55 : #include "foxglynn.hpp"
56 :
57 : #include <stdio.h>
58 : #include <stdlib.h>
59 : #include <math.h>
60 :
61 : /*The pi constant*/
62 : static const double pi = 3.14159265358979323846264;
63 : static const double lambda_25 = 25.0;
64 : static const double lambda_400 = 400.0;
65 :
66 : /*****************************************************************************
67 : Name : finder
68 : Role : The FINDER function from the Fox-Glynn algorithm
69 : @param : int m: the middle point from weighter.
70 : @param : double lambda: (rate of uniformization)*(mission time)
71 : @param : double tau: underflow
72 : @param : double omega: overflow
73 : @param : double epsilon: error bound
74 : @param : double * pw_m: pointer to the initial value in the array of sets, will be set here
75 : @param : FoxGlynn *: return by reference
76 : @return : TRUE if everything is fine, otherwise FALSE.
77 : This is the F parameter of Fox-Glynn finder function.
78 : remark : This implementation is based on the original Fox-Glynn paper with
79 : small modifications.
80 : The most remarkable alterations are the removal of terminating indices for the search for the
81 : truncation points. This will result in higher accuracy when needed.
82 : It can be shown that the introduced error converges to a value below the requested error level
83 : (for epsilon > 0). Hence the algorithm will always terminate.
84 : ******************************************************************************/
85 0 : bool FoxGlynn::finder(const int m, const double lambda, const double tau, const double omega,
86 : const double epsilon, double * pw_m)
87 : {
88 0 : const double sqrt_2_pi = sqrt( 2.0 * pi );
89 0 : const double sqrt_2 = sqrt(2.0);
90 0 : const double sqrt_lambda = sqrt(lambda);
91 0 : double lambda_max, k, k_rtp = HUGE_VAL, k_prime, c_m_inf, result, al, dkl, bl;
92 :
93 : /*Simple bad cases, when we quit*/
94 0 : if( lambda == 0.0 )
95 : {
96 0 : fprintf(stderr,"ERROR: Fox-Glynn: lambda = 0, terminating the algorithm\n");
97 0 : return false;
98 : }
99 : /* The requested error level must not be smaller than the minimum machine precision
100 : (needed to guarantee the convergence of the error conditions) */
101 0 : if( epsilon < tau)
102 : {
103 0 : fprintf(stderr,"ERROR: Fox-Glynn: epsilon < tau, invalid error level, terminating the algorithm\n");
104 0 : fprintf(stderr,"epsilon %e, tau %e\n",epsilon,tau);
105 0 : return false;
106 : }
107 : /* zero is used as left truncation point for lambda <= 25 */
108 0 : left = 0;
109 0 : lambda_max = lambda;
110 :
111 : /* for lambda below 25 the exponential can be smaller than tau */
112 : /* if that is the case we expect underflows and warn the user */
113 0 : if( 0.0 < lambda && lambda <= lambda_25 )
114 : {
115 0 : if( exp( -lambda ) <= tau )
116 : {
117 0 : fprintf(stderr,"ERROR: Fox-Glynn: 0 < lambda < 25, underflow. The results are UNRELIABLE.\n");
118 : }
119 : }
120 :
121 0 : bl = (1.0 + 1.0/lambda) * exp(1.0 / (8.0 * lambda));
122 :
123 : /****Compute pFG->right truncation point****/
124 : /*According to Fox-Glynn, if lambda < 400 we should take lambda = 400,
125 : otherwise use the original value. This is for computing the right truncation point*/
126 0 : if(lambda < lambda_400)
127 0 : lambda_max = lambda_400;
128 0 : k = 4;
129 0 : al = (1.0+1.0/lambda_max) * exp(1.0/16.0) * sqrt_2;
130 0 : dkl = exp(-2.0/9.0 * (k * sqrt(2.0 * lambda_max) + 1.5 ));
131 0 : dkl = 1.0 / (1.0 - dkl);
132 : /* find right truncation point */
133 :
134 : /* This loop is a modification to the original Fox-Glynn paper.
135 : The search for the right truncation point is only terminated by
136 : the error condition and not by the stop index from the FG paper.
137 : This can yield more accurate results if neccesary.*/
138 0 : while((epsilon/2.0) < ((al * dkl * exp(-(k*k)/2.0))/(k*sqrt_2_pi)))
139 : {
140 0 : k++;
141 0 : dkl = exp(-2.0/9.0 * (k * sqrt_2 * sqrt(lambda_max) + 1.5 ));
142 0 : dkl = 1.0 / (1.0 - dkl);
143 : }
144 0 : k_rtp = k;
145 0 : right = (int)ceil(m + k_rtp * sqrt_2 * sqrt(lambda_max) + 1.5 );
146 :
147 :
148 : /****Compute pFG->left truncation point****/
149 : /* compute the left truncation point for lambda > 25 */
150 : /* for lambda <= 25 we use zero as left truncation point */
151 0 : if(lambda > lambda_25)
152 : {
153 : /*Start looking for the left truncation point*/
154 : /* start search at k=4 (taken from original Fox-Glynn paper */
155 0 : k = 4;
156 : /* increase the left truncation point as long as we fulfill the error condition */
157 :
158 : /* This loop is a modification to the original Fox-Glynn paper.
159 : The search for the left truncation point is only terminated by
160 : the error condition and not by the stop index from the FG paper.
161 : This can yield more accurate results if neccesary.*/
162 0 : while((epsilon/2.0) < ((bl * exp(-(k*k)/2.0))/(k * sqrt_2_pi)))
163 0 : k++;
164 : /*Finally the left truncation point is found*/
165 0 : left = (int)floor(m - k*sqrt(lambda)- 1.5 );
166 : /* for small lambda the above calculation can yield negative truncation points, crop them here */
167 0 : if(left < 0)left = 0;
168 : /* perform underflow check */
169 0 : k_prime = k + 3.0 / (2.0 * sqrt_lambda);
170 : /*We take the c_m_inf = 0.02935 / sqrt( m ), as for lambda >= 25
171 : c_m = 1 / ( sqrt( 2.0 * pi * m ) ) * exp( m - lambda - 1 / ( 12.0 * m ) ) => c_m_inf*/
172 0 : c_m_inf = 0.02935 / sqrt((double) m);
173 0 : result = 0.0;
174 0 : if( 0.0 < k_prime && k_prime <= sqrt_lambda / 2.0 )
175 : {
176 0 : result = c_m_inf * exp( - pow(k_prime,2.0) / 2.0 - pow(k_prime, 3.0) / (3.0 * sqrt_lambda) );
177 : }
178 : else
179 : {
180 0 : if( k_prime <= sqrt( m + 1.0 ) / m )
181 : {
182 0 : double result_1 = c_m_inf * pow(
183 0 : 1.0 - k_prime / sqrt((double) (m + 1)),
184 0 : k_prime * sqrt((double) (m + 1)));
185 0 : double result_2 = exp( - lambda );
186 : /*Take the maximum*/
187 0 : result = ( result_1 > result_2 ? result_1 : result_2);
188 : }
189 : else
190 : {
191 : /*NOTE: It will be an underflow error*/;
192 0 : fprintf(stderr,"ERROR: Fox-Glynn: lambda >= 25, underflow. The results are UNRELIABLE.\n");
193 : }
194 : }
195 0 : if ( result * omega / ( 1.0e+10 * ( right - left ) ) <= tau )
196 : {
197 0 : fprintf(stderr,"ERROR: Fox-Glynn: lambda >= 25, underflow. The results are UNRELIABLE.\n");
198 : }
199 : }
200 :
201 :
202 :
203 : /*We still have to perform an underflow check for the right truncation point when lambda >= 400*/
204 0 : if( lambda >= lambda_400 )
205 : {
206 0 : k_prime = k_rtp * sqrt_2 + 3.0 / (2.0 * sqrt_lambda);
207 : /*We take the c_m_inf = 0.02935 / sqrt( m ), as for lambda >= 25
208 : c_m = 1 / ( sqrt( 2.0 * pi * m ) ) * exp( m - lambda - 1 / ( 12.0 * m ) ) => c_m_inf*/
209 0 : c_m_inf = 0.02935 / sqrt((double) m);
210 0 : result = c_m_inf * exp( - pow( k_prime + 1.0 , 2.0 ) / 2.0 );
211 0 : if( result * omega / ( 1.0e+10 * ( right - left ) ) <= tau)
212 : {
213 0 : fprintf(stderr,"ERROR: Fox-Glynn: lambda >= 400, underflow. The results are UNRELIABLE.\n");
214 : }
215 : }
216 : /*Time to set the initial value for weights*/
217 0 : *pw_m = omega / ( 1.0e+10 * ( right - left ) );
218 :
219 0 : return true;
220 : }
221 :
222 : /*****************************************************************************
223 : Name : weighter
224 : Role : The WEIGHTER function from the Fox-Glynn algorithm
225 : @param : double lambda: (rate of uniformization)*(mission time)
226 : @param : double tau: underflow
227 : @param : double omega: overflow
228 : @param : double epsilon: error bound
229 : @param : FoxGlynn *: return by reference
230 : @return : TRUE if everything is fine, otherwise FALSE.
231 : This is the F parameter of Fox-Glynn finder function.
232 : remark :
233 : ******************************************************************************/
234 0 : bool FoxGlynn::weighter(const double lambda, const double tau, const double omega, const double epsilon)
235 : {
236 : /*The magic m point*/
237 0 : const int m = (int)floor(lambda);
238 0 : double w_m = 0;
239 : int j, s, t;
240 :
241 0 : if( ! finder( m, lambda, tau, omega, epsilon, &w_m ) )
242 0 : return false;
243 :
244 : /*Allocate space for weights*/
245 0 : weights = (double *) calloc((size_t) (right - left + 1),
246 : sizeof(double));
247 : /*Set an initial weight*/
248 0 : weights[ m - left ] = w_m;
249 :
250 : /*Fill the left side of the array*/
251 0 : for( j = m; j > left; j-- )
252 0 : weights[ ( j - left ) - 1 ] = ( j / lambda ) * weights[ j - left ];
253 :
254 : /*Fill the right side of the array, have two cases lambda < 400 & lambda >= 400*/
255 0 : if( lambda < lambda_400 )
256 : {
257 : /*Perform the underflow check, according to Fox-Glynn*/
258 0 : if( right > 600 )
259 : {
260 0 : fprintf(stderr,"ERROR: Fox-Glynn: pFG->right > 600, underflow is possible\n");
261 0 : return false;
262 : }
263 : /*Compute weights*/
264 0 : for( j = m; j < right; j++ )
265 : {
266 0 : double q = lambda / ( j + 1 );
267 0 : if( weights[ j - left ] > tau / q )
268 : {
269 0 : weights[ ( j - left ) + 1 ] = q * weights[ j - left ];
270 : }else{
271 0 : right = j;
272 0 : break; /*It's time to compute W*/
273 : }
274 : }
275 : }else{
276 : /*Compute weights*/
277 0 : for( j = m; j < right; j++ )
278 0 : weights[ ( j - left ) + 1 ] = ( lambda / ( j + 1 ) ) * weights[ j - left ];
279 : }
280 :
281 : /*It is time to compute the normalization weight W*/
282 0 : total_weight = 0.0;
283 0 : s = left;
284 0 : t = right;
285 :
286 0 : while( s < t )
287 : {
288 0 : if( weights[ s - left ] <= weights[ t - left ] )
289 : {
290 0 : total_weight += weights[ s - left ];
291 0 : s++;
292 : }else{
293 0 : total_weight += weights[ t - left ];
294 0 : t--;
295 : }
296 : }
297 0 : total_weight += weights[ s - left ];
298 :
299 : /* printf("Fox-Glynn: ltp = %d, rtp = %d, w = %10.15le \n", pFG->left, pFG->right, pFG->total_weight); */
300 :
301 0 : return true;
302 : }
303 :
304 : /*****************************************************************************
305 : Name : fox_glynn
306 : Role : get poisson probabilities.
307 : @param : double lambda: (rate of uniformization)*(mission time)
308 : @param : double tau: underflow
309 : @param : double omega: overflow
310 : @param : double epsilon: error bound
311 : @param : FoxGlynn **: return a new FoxGlynn structure by reference
312 : @return : TRUE if it worked fine, otherwise false
313 : remark :
314 : ******************************************************************************/
315 0 : FoxGlynn::FoxGlynn(const double lambda, const double tau, const double omega, const double epsilon):
316 0 : isValid(false),left(0),right(0),total_weight(0.0),weights(NULL)
317 : {
318 : /* printf("Fox-Glynn: lambda = %3.3le, epsilon = %1.8le\n",lambda, epsilon); */
319 0 : isValid = weighter(lambda, tau, omega, epsilon);
320 0 : }
321 :
322 :
323 : /**
324 : * Fries the memory allocated for the FoxGlynn structure
325 : * @param fg the structure to free
326 : */
327 0 : FoxGlynn::~FoxGlynn()
328 : {
329 0 : if( weights )free(weights);
330 0 : }
|
