chem/KPP/kpp/kpp-2.1/int/tau_leap.f90

   1 MODULE KPP_ROOT_Random
   2 ! A module for random number generation from the following distributions:
   3 !
   4 !     Distribution                    Function/subroutine name
   5 !
   6 !     Normal (Gaussian)               random_normal
   7 !     Gamma                           random_gamma
   8 !     Chi-squared                     random_chisq
   9 !     Exponential                     random_exponential
  10 !     Weibull                         random_Weibull
  11 !     Beta                            random_beta
  12 !     t                               random_t
  13 !     Multivariate normal             random_mvnorm
  14 !     Generalized inverse Gaussian    random_inv_gauss
  15 !     Poisson                         random_Poisson
  16 !     Binomial                        random_binomial1   *
  17 !                                     random_binomial2   *
  18 !     Negative binomial               random_neg_binomial
  19 !     von Mises                       random_von_Mises
  20 !     Cauchy                          random_Cauchy
  21 !
  22 !  Generate a random ordering of the integers 1 .. N
  23 !                                     random_order
  24 !     Initialize (seed) the uniform random number generator for ANY compiler
  25 !                                     seed_random_number
  26
  27 !     Lognormal - see note below.
  28
  29 !  ** Two functions are provided for the binomial distribution.
  30 !  If the parameter values remain constant, it is recommended that the
  31 !  first function is used (random_binomial1).   If one or both of the
  32 !  parameters change, use the second function (random_binomial2).
  33
  34 ! The compilers own random number generator, SUBROUTINE RANDOM_NUMBER(r),
  35 ! is used to provide a source of uniformly distributed random numbers.
  36
  37 ! N.B. At this stage, only one random number is generated at each call to
  38 !      one of the functions above.
  39
  40 ! The module uses the following functions which are included here:
  41 ! bin_prob to calculate a single binomial probability
  42 ! lngamma  to calculate the logarithm to base e of the gamma function
  43
  44 ! Some of the code is adapted from Dagpunar's book:
  45 !     Dagpunar, J. 'Principles of random variate generation'
  46 !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  47 !
  48 ! In most of Dagpunar's routines, there is a test to see whether the value
  49 ! of one or two floating-point parameters has changed since the last call.
  50 ! These tests have been replaced by using a logical variable FIRST.
  51 ! This should be set to .TRUE. on the first call using new values of the
  52 ! parameters, and .FALSE. if the parameter values are the same as for the
  53 ! previous call.
  54
  55 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  56 ! Lognormal distribution
  57 ! If X has a lognormal distribution, then log(X) is normally distributed.
  58 ! Here the logarithm is the natural logarithm, that is to base e, sometimes
  59 ! denoted as ln.  To generate random variates from this distribution, generate
  60 ! a random deviate from the normal distribution with mean and variance equal
  61 ! to the mean and variance of the logarithms of X, then take its exponential.
  62
  63 ! Relationship between the mean & variance of log(X) and the mean & variance
  64 ! of X, when X has a lognormal distribution.
  65 ! Let m = mean of log(X), and s^2 = variance of log(X)
  66 ! Then
  67 ! mean of X     = exp(m + 0.5s^2)
  68 ! variance of X = (mean(X))^2.[exp(s^2) - 1]
  69
  70 ! In the reverse direction (rarely used)
  71 ! variance of log(X) = log[1 + var(X)/(mean(X))^2]
  72 ! mean of log(X)     = log(mean(X) - 0.5var(log(X))
  73
  74 ! N.B. The above formulae relate to population parameters; they will only be
  75 !      approximate if applied to sample values.
  76 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  77
  78 ! Version 1.13, 2 October 2000
  79 ! Changes from version 1.01
  80 ! 1. The random_order, random_Poisson & random_binomial routines have been
  81 !    replaced with more efficient routines.
  82 ! 2. A routine, seed_random_number, has been added to seed the uniform random
  83 !    number generator.   This requires input of the required number of seeds
  84 !    for the particular compiler from a specified I/O unit such as a keyboard.
  85 ! 3. Made compatible with Lahey's ELF90.
  86 ! 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1.
  87 ! 5. INTENT for array f corrected in random_mvnorm.
  88
  89 !     Author: Alan Miller
  90 !             CSIRO Division of Mathematical & Information Sciences
  91 !             Private Bag 10, Clayton South MDC
  92 !             Clayton 3169, Victoria, Australia
  93 !     Phone: (+61) 3 9545-8016      Fax: (+61) 3 9545-8080
  94 !     e-mail: amiller @ bigpond.net.au
  95
  96 IMPLICIT NONE
  97 REAL, PRIVATE      :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0,   &
  98                       vsmall = TINY(1.0), vlarge = HUGE(1.0)
  99 PRIVATE            :: integral
 100 INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
 101
 102
 103 CONTAINS
 104
 105
 106 FUNCTION random_normal() RESULT(fn_val)
 107
 108 ! Adapted from the following Fortran 77 code
 109 !      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
 110 !      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
 111 !      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
 112
 113 !  The function random_normal() returns a normally distributed pseudo-random
 114 !  number with zero mean and unit variance.
 115
 116 !  The algorithm uses the ratio of uniforms method of A.J. Kinderman
 117 !  and J.F. Monahan augmented with quadratic bounding curves.
 118
 119 REAL :: fn_val
 120
 121 !     Local variables
 122 REAL     :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472,           &
 123             r1 = 0.27597, r2 = 0.27846, u, v, x, y, q
 124
 125 !     Generate P = (u,v) uniform in rectangle enclosing acceptance region
 126
 127 DO
 128   CALL RANDOM_NUMBER(u)
 129   CALL RANDOM_NUMBER(v)
 130   v = 1.7156 * (v - half)
 131
 132 !     Evaluate the quadratic form
 133   x = u - s
 134   y = ABS(v) - t
 135   q = x**2 + y*(a*y - b*x)
 136
 137 !     Accept P if inside inner ellipse
 138   IF (q < r1) EXIT
 139 !     Reject P if outside outer ellipse
 140   IF (q > r2) CYCLE
 141 !     Reject P if outside acceptance region
 142   IF (v**2 < -4.0*LOG(u)*u**2) EXIT
 143 END DO
 144
 145 !     Return ratio of P's coordinates as the normal deviate
 146 fn_val = v/u
 147 RETURN
 148
 149 END FUNCTION random_normal
 150
 151
 152
 153 FUNCTION random_gamma(s, first) RESULT(fn_val)
 154
 155 ! Adapted from Fortran 77 code from the book:
 156 !     Dagpunar, J. 'Principles of random variate generation'
 157 !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
 158
 159 !     FUNCTION GENERATES A RANDOM GAMMA VARIATE.
 160 !     CALLS EITHER random_gamma1 (S > 1.0)
 161 !     OR random_exponential (S = 1.0)
 162 !     OR random_gamma2 (S < 1.0).
 163
 164 !     S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL).
 165
 166 REAL, INTENT(IN)    :: s
 167 LOGICAL, INTENT(IN) :: first
 168 REAL                :: fn_val
 169
 170 IF (s <= zero) THEN
 171   WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE'
 172   STOP
 173 END IF
 174
 175 IF (s > one) THEN
 176   fn_val = random_gamma1(s, first)
 177 ELSE IF (s < one) THEN
 178   fn_val = random_gamma2(s, first)
 179 ELSE
 180   fn_val = random_exponential()
 181 END IF
 182
 183 RETURN
 184 END FUNCTION random_gamma
 185
 186
 187
 188 FUNCTION random_gamma1(s, first) RESULT(fn_val)
 189
 190 ! Uses the algorithm in
 191 ! Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating
 192 ! gamma variables', Trans. om Math. Software (TOMS), vol.26(3), pp.363-372.
 193
 194 ! Generates a random gamma deviate for shape parameter s >= 1.
 195
 196 REAL, INTENT(IN)    :: s
 197 LOGICAL, INTENT(IN) :: first
 198 REAL                :: fn_val
 199
 200 ! Local variables
 201 REAL, SAVE  :: c, d
 202 REAL        :: u, v, x
 203
 204 IF (first) THEN
 205   d = s - one/3.
 206   c = one/SQRT(9.0*d)
 207 END IF
 208
 209 ! Start of main loop
 210 DO
 211
 212 ! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0.
 213
 214   DO
 215     x = random_normal()
 216     v = (one + c*x)**3
 217     IF (v > zero) EXIT
 218   END DO
 219
 220 ! Generate uniform variable U
 221
 222   CALL RANDOM_NUMBER(u)
 223   IF (u < one - 0.0331*x**4) THEN
 224     fn_val = d*v
 225     EXIT
 226   ELSE IF (LOG(u) < half*x**2 + d*(one - v + LOG(v))) THEN
 227     fn_val = d*v
 228     EXIT
 229   END IF
 230 END DO
 231
 232 RETURN
 233 END FUNCTION random_gamma1
 234
 235
 236
 237 FUNCTION random_gamma2(s, first) RESULT(fn_val)
 238
 239 ! Adapted from Fortran 77 code from the book:
 240 !     Dagpunar, J. 'Principles of random variate generation'
 241 !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
 242
 243 ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
 244 ! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO
 245 ! GAMMA2**(S-1) * EXP(-GAMMA2),
 246 ! USING A SWITCHING METHOD.
 247
 248 !    S = SHAPE PARAMETER OF DISTRIBUTION
 249 !          (REAL < 1.0)
 250
 251 REAL, INTENT(IN)    :: s
 252 LOGICAL, INTENT(IN) :: first
 253 REAL                :: fn_val
 254
 255 !     Local variables
 256 REAL       :: r, x, w
 257 REAL, SAVE :: a, p, c, uf, vr, d
 258
 259 IF (s <= zero .OR. s >= one) THEN
 260   WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE'
 261   STOP
 262 END IF
 263
 264 IF (first) THEN                        ! Initialization, if necessary
 265   a = one - s
 266   p = a/(a + s*EXP(-a))
 267   IF (s < vsmall) THEN
 268     WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL'
 269     STOP
 270   END IF
 271   c = one/s
 272   uf = p*(vsmall/a)**s
 273   vr = one - vsmall
 274   d = a*LOG(a)
 275 END IF
 276
 277 DO
 278   CALL RANDOM_NUMBER(r)
 279   IF (r >= vr) THEN
 280     CYCLE
 281   ELSE IF (r > p) THEN
 282     x = a - LOG((one - r)/(one - p))
 283     w = a*LOG(x)-d
 284   ELSE IF (r > uf) THEN
 285     x = a*(r/p)**c
 286     w = x
 287   ELSE
 288     fn_val = zero
 289     RETURN
 290   END IF
 291
 292   CALL RANDOM_NUMBER(r)
 293   IF (one-r <= w .AND. r > zero) THEN
 294     IF (r*(w + one) >= one) CYCLE
 295     IF (-LOG(r) <= w) CYCLE
 296   END IF
 297   EXIT
 298 END DO
 299
 300 fn_val = x
 301 RETURN
 302
 303 END FUNCTION random_gamma2
 304
 305
 306
 307 FUNCTION random_chisq(ndf, first) RESULT(fn_val)
 308
 309 !     Generates a random variate from the chi-squared distribution with
 310 !     ndf degrees of freedom
 311
 312 INTEGER, INTENT(IN) :: ndf
 313 LOGICAL, INTENT(IN) :: first
 314 REAL                :: fn_val
 315
 316 fn_val = two * random_gamma(half*ndf, first)
 317 RETURN
 318
 319 END FUNCTION random_chisq
 320
 321
 322
 323 FUNCTION random_exponential() RESULT(fn_val)
 324
 325 ! Adapted from Fortran 77 code from the book:
 326 !     Dagpunar, J. 'Principles of random variate generation'
 327 !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
 328
 329 ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
 330 ! A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL
 331 ! TO EXP(-random_exponential), USING INVERSION.
 332
 333 REAL  :: fn_val
 334
 335 !     Local variable
 336 REAL  :: r
 337
 338 DO
 339   CALL RANDOM_NUMBER(r)
 340   IF (r > zero) EXIT
 341 END DO
 342
 343 fn_val = -LOG(r)
 344 RETURN
 345
 346 END FUNCTION random_exponential
 347
 348
 349
 350 FUNCTION random_Weibull(a) RESULT(fn_val)
 351
 352 !     Generates a random variate from the Weibull distribution with
 353 !     probability density:
 354 !                      a
 355 !               a-1  -x
 356 !     f(x) = a.x    e
 357
 358 REAL, INTENT(IN) :: a
 359 REAL             :: fn_val
 360
 361 !     For speed, there is no checking that a is not zero or very small.
 362
 363 fn_val = random_exponential() ** (one/a)
 364 RETURN
 365
 366 END FUNCTION random_Weibull
 367
 368
 369
 370 FUNCTION random_beta(aa, bb, first) RESULT(fn_val)
 371
 372 ! Adapted from Fortran 77 code from the book:
 373 !     Dagpunar, J. 'Principles of random variate generation'
 374 !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
 375
 376 ! FUNCTION GENERATES A RANDOM VARIATE IN [0,1]
 377 ! FROM A BETA DISTRIBUTION WITH DENSITY
 378 ! PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1).
 379 ! USING CHENG'S LOG LOGISTIC METHOD.
 380
 381 !     AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
 382 !     BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
 383
 384 REAL, INTENT(IN)    :: aa, bb
 385 LOGICAL, INTENT(IN) :: first
 386 REAL                :: fn_val
 387
 388 !     Local variables
 389 REAL, PARAMETER  :: aln4 = 1.3862944
 390 REAL             :: a, b, g, r, s, x, y, z
 391 REAL, SAVE       :: d, f, h, t, c
 392 LOGICAL, SAVE    :: swap
 393
 394 IF (aa <= zero .OR. bb <= zero) THEN
 395   WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE(S)'
 396   STOP
 397 END IF
 398
 399 IF (first) THEN                        ! Initialization, if necessary
 400   a = aa
 401   b = bb
 402   swap = b > a
 403   IF (swap) THEN
 404     g = b
 405     b = a
 406     a = g
 407   END IF
 408   d = a/b
 409   f = a+b
 410   IF (b > one) THEN
 411     h = SQRT((two*a*b - f)/(f - two))
 412     t = one
 413   ELSE
 414     h = b
 415     t = one/(one + (a/(vlarge*b))**b)
 416   END IF
 417   c = a+h
 418 END IF
 419
 420 DO
 421   CALL RANDOM_NUMBER(r)
 422   CALL RANDOM_NUMBER(x)
 423   s = r*r*x
 424   IF (r < vsmall .OR. s <= zero) CYCLE
 425   IF (r < t) THEN
 426     x = LOG(r/(one - r))/h
 427     y = d*EXP(x)
 428     z = c*x + f*LOG((one + d)/(one + y)) - aln4
 429     IF (s - one > z) THEN
 430       IF (s - s*z > one) CYCLE
 431       IF (LOG(s) > z) CYCLE
 432     END IF
 433     fn_val = y/(one + y)
 434   ELSE
 435     IF (4.0*s > (one + one/d)**f) CYCLE
 436     fn_val = one
 437   END IF
 438   EXIT
 439 END DO
 440
 441 IF (swap) fn_val = one - fn_val
 442 RETURN
 443 END FUNCTION random_beta
 444
 445
 446
 447 FUNCTION random_t(m) RESULT(fn_val)
 448
 449 ! Adapted from Fortran 77 code from the book:
 450 !     Dagpunar, J. 'Principles of random variate generation'
 451 !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
 452
 453 ! FUNCTION GENERATES A RANDOM VARIATE FROM A
 454 ! T DISTRIBUTION USING KINDERMAN AND MONAHAN'S RATIO METHOD.
 455
 456 !     M = DEGREES OF FREEDOM OF DISTRIBUTION
 457 !           (1 <= 1NTEGER)
 458
 459 INTEGER, INTENT(IN) :: m
 460 REAL                :: fn_val
 461
 462 !     Local variables
 463 REAL, SAVE      :: s, c, a, f, g
 464 REAL            :: r, x, v
 465
 466 REAL, PARAMETER :: three = 3.0, four = 4.0, quart = 0.25,   &
 467                    five = 5.0, sixteen = 16.0
 468 INTEGER         :: mm = 0
 469
 470 IF (m < 1) THEN
 471   WRITE(*, *) 'IMPERMISSIBLE DEGREES OF FREEDOM'
 472   STOP
 473 END IF
 474
 475 IF (m /= mm) THEN                    ! Initialization, if necessary
 476   s = m
 477   c = -quart*(s + one)
 478   a = four/(one + one/s)**c
 479   f = sixteen/a
 480   IF (m > 1) THEN
 481     g = s - one
 482     g = ((s + one)/g)**c*SQRT((s+s)/g)
 483   ELSE
 484     g = one
 485   END IF
 486   mm = m
 487 END IF
 488
 489 DO
 490   CALL RANDOM_NUMBER(r)
 491   IF (r <= zero) CYCLE
 492   CALL RANDOM_NUMBER(v)
 493   x = (two*v - one)*g/r
 494   v = x*x
 495   IF (v > five - a*r) THEN
 496     IF (m >= 1 .AND. r*(v + three) > f) CYCLE
 497     IF (r > (one + v/s)**c) CYCLE
 498   END IF
 499   EXIT
 500 END DO
 501
 502 fn_val = x
 503 RETURN
 504 END FUNCTION random_t
 505
 506
 507
 508 SUBROUTINE random_mvnorm(n, h, d, f, first, x, ier)
 509
 510 ! Adapted from Fortran 77 code from the book:
 511 !     Dagpunar, J. 'Principles of random variate generation'
 512 !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
 513
 514 ! N.B. An extra argument, ier, has been added to Dagpunar's routine
 515
 516 !     SUBROUTINE GENERATES AN N VARIATE RANDOM NORMAL
 517 !     VECTOR USING A CHOLESKY DECOMPOSITION.
 518
 519 ! ARGUMENTS:
 520 !        N = NUMBER OF VARIATES IN VECTOR
 521 !           (INPUT,INTEGER >= 1)
 522 !     H(J) = J'TH ELEMENT OF VECTOR OF MEANS
 523 !           (INPUT,REAL)
 524 !     X(J) = J'TH ELEMENT OF DELIVERED VECTOR
 525 !           (OUTPUT,REAL)
 526 !
 527 !    D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I)
 528 !            (INPUT,REAL)
 529 !    F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR
 530 !           DECOMPOSITION OF VARIANCE MATRIX (J <= I)
 531 !            (OUTPUT,REAL)
 532
 533 !    FIRST = .TRUE. IF THIS IS THE FIRST CALL OF THE ROUTINE
 534 !    OR IF THE DISTRIBUTION HAS CHANGED SINCE THE LAST CALL OF THE ROUTINE.
 535 !    OTHERWISE SET TO .FALSE.
 536 !            (INPUT,LOGICAL)
 537
 538 !    ier = 1 if the input covariance matrix is not +ve definite
 539 !        = 0 otherwise
 540
 541 INTEGER, INTENT(IN)   :: n
 542 REAL, INTENT(IN)      :: h(:), d(:)   ! d(n*(n+1)/2)
 543 REAL, INTENT(IN OUT)  :: f(:)         ! f(n*(n+1)/2)
 544 REAL, INTENT(OUT)     :: x(:)
 545 LOGICAL, INTENT(IN)   :: first
 546 INTEGER, INTENT(OUT)  :: ier
 547
 548 !     Local variables
 549 INTEGER       :: j, i, m
 550 REAL          :: y, v
 551 INTEGER, SAVE :: n2
 552
 553 IF (n < 1) THEN
 554   WRITE(*, *) 'SIZE OF VECTOR IS NON POSITIVE'
 555   STOP
 556 END IF
 557
 558 ier = 0
 559 IF (first) THEN                        ! Initialization, if necessary
 560   n2 = 2*n
 561   IF (d(1) < zero) THEN
 562     ier = 1
 563     RETURN
 564   END IF
 565
 566   f(1) = SQRT(d(1))
 567   y = one/f(1)
 568   DO j = 2,n
 569     f(j) = d(1+j*(j-1)/2) * y
 570   END DO
 571
 572   DO i = 2,n
 573     v = d(i*(i-1)/2+i)
 574     DO m = 1,i-1
 575       v = v - f((m-1)*(n2-m)/2+i)**2
 576     END DO
 577
 578     IF (v < zero) THEN
 579       ier = 1
 580       RETURN
 581     END IF
 582
 583     v = SQRT(v)
 584     y = one/v
 585     f((i-1)*(n2-i)/2+i) = v
 586     DO j = i+1,n
 587       v = d(j*(j-1)/2+i)
 588       DO m = 1,i-1
 589         v = v - f((m-1)*(n2-m)/2+i)*f((m-1)*(n2-m)/2 + j)
 590       END DO ! m = 1,i-1
 591       f((i-1)*(n2-i)/2 + j) = v*y
 592     END DO ! j = i+1,n
 593   END DO ! i = 2,n
 594 END IF
 595
 596 x(1:n) = h(1:n)
 597 DO j = 1,n
 598   y = random_normal()
 599   DO i = j,n
 600     x(i) = x(i) + f((j-1)*(n2-j)/2 + i) * y
 601   END DO ! i = j,n
 602 END DO ! j = 1,n
 603
 604 RETURN
 605 END SUBROUTINE random_mvnorm
 606
 607
 608
 609 FUNCTION random_inv_gauss(h, b, first) RESULT(fn_val)
 610
 611 ! Adapted from Fortran 77 code from the book:
 612 !     Dagpunar, J. 'Principles of random variate generation'
 613 !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
 614
 615 ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY] FROM
 616 ! A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION
 617 ! WITH DENSITY PROPORTIONAL TO  GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG))
 618 ! USING A RATIO METHOD.
 619
 620 !     H = PARAMETER OF DISTRIBUTION (0 <= REAL)
 621 !     B = PARAMETER OF DISTRIBUTION (0 < REAL)
 622
 623 REAL, INTENT(IN)    :: h, b
 624 LOGICAL, INTENT(IN) :: first
 625 REAL                :: fn_val
 626
 627 !     Local variables
 628 REAL            :: ym, xm, r, w, r1, r2, x
 629 REAL, SAVE      :: a, c, d, e
 630 REAL, PARAMETER :: quart = 0.25
 631
 632 IF (h < zero .OR. b <= zero) THEN
 633   WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'
 634   STOP
 635 END IF
 636
 637 IF (first) THEN                        ! Initialization, if necessary
 638   IF (h > quart*b*SQRT(vlarge)) THEN
 639     WRITE(*, *) 'THE RATIO H:B IS TOO SMALL'
 640     STOP
 641   END IF
 642   e = b*b
 643   d = h + one
 644   ym = (-d + SQRT(d*d + e))/b
 645   IF (ym < vsmall) THEN
 646     WRITE(*, *) 'THE VALUE OF B IS TOO SMALL'
 647     STOP
 648   END IF
 649
 650   d = h - one
 651   xm = (d + SQRT(d*d + e))/b
 652   d = half*d
 653   e = -quart*b
 654   r = xm + one/xm
 655   w = xm*ym
 656   a = w**(-half*h) * SQRT(xm/ym) * EXP(-e*(r - ym - one/ym))
 657   IF (a < vsmall) THEN
 658     WRITE(*, *) 'THE VALUE OF H IS TOO LARGE'
 659     STOP
 660   END IF
 661   c = -d*LOG(xm) - e*r
 662 END IF
 663
 664 DO
 665   CALL RANDOM_NUMBER(r1)
 666   IF (r1 <= zero) CYCLE
 667   CALL RANDOM_NUMBER(r2)
 668   x = a*r2/r1
 669   IF (x <= zero) CYCLE
 670   IF (LOG(r1) < d*LOG(x) + e*(x + one/x) + c) EXIT
 671 END DO
 672
 673 fn_val = x
 674
 675 RETURN
 676 END FUNCTION random_inv_gauss
 677
 678
 679
 680 FUNCTION random_Poisson(mu, first) RESULT(ival)
 681 !**********************************************************************
 682 !     Translated to Fortran 90 by Alan Miller from:
 683 !                           RANLIB
 684 !
 685 !     Library of Fortran Routines for Random Number Generation
 686 !
 687 !                    Compiled and Written by:
 688 !
 689 !                         Barry W. Brown
 690 !                          James Lovato
 691 !
 692 !             Department of Biomathematics, Box 237
 693 !             The University of Texas, M.D. Anderson Cancer Center
 694 !             1515 Holcombe Boulevard
 695 !             Houston, TX      77030
 696 !
 697 ! This work was supported by grant CA-16672 from the National Cancer Institute.
 698
 699 !                    GENerate POIsson random deviate
 700
 701 !                            Function
 702
 703 ! Generates a single random deviate from a Poisson distribution with mean mu.
 704
 705 !                            Arguments
 706
 707 !     mu --> The mean of the Poisson distribution from which
 708 !            a random deviate is to be generated.
 709 !                              REAL mu
 710
 711 !                              Method
 712
 713 !     For details see:
 714
 715 !               Ahrens, J.H. and Dieter, U.
 716 !               Computer Generation of Poisson Deviates
 717 !               From Modified Normal Distributions.
 718 !               ACM Trans. Math. Software, 8, 2
 719 !               (June 1982),163-179
 720
 721 !     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
 722 !     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
 723
 724 !     SEPARATION OF CASES A AND B
 725
 726 !     .. Scalar Arguments ..
 727 REAL, INTENT(IN)    :: mu
 728 LOGICAL, INTENT(IN) :: first
 729 INTEGER             :: ival
 730 !     ..
 731 !     .. Local Scalars ..
 732 REAL          :: b1, b2, c, c0, c1, c2, c3, del, difmuk, e, fk, fx, fy, g,  &
 733                  omega, px, py, t, u, v, x, xx
 734 REAL, SAVE    :: s, d, p, q, p0
 735 INTEGER       :: j, k, kflag
 736 LOGICAL, SAVE :: full_init
 737 INTEGER, SAVE :: l, m
 738 !     ..
 739 !     .. Local Arrays ..
 740 REAL, SAVE    :: pp(35)
 741 !     ..
 742 !     .. Data statements ..
 743 REAL, PARAMETER :: a0 = -.5, a1 = .3333333, a2 = -.2500068, a3 = .2000118,  &
 744                    a4 = -.1661269, a5 = .1421878, a6 = -.1384794,   &
 745                    a7 = .1250060
 746
 747 REAL, PARAMETER :: fact(10) = (/ 1., 1., 2., 6., 24., 120., 720., 5040.,  &
 748                                  40320., 362880. /)
 749
 750 !     ..
 751 !     .. Executable Statements ..
 752 IF (mu > 10.0) THEN
 753 !     C A S E  A. (RECALCULATION OF S, D, L IF MU HAS CHANGED)
 754
 755   IF (first) THEN
 756     s = SQRT(mu)
 757     d = 6.0*mu*mu
 758
 759 !             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
 760 !             PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
 761 !             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
 762
 763     l = mu - 1.1484
 764     full_init = .false.
 765   END IF
 766
 767
 768 !     STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE
 769
 770   g = mu + s*random_normal()
 771   IF (g > 0.0) THEN
 772     ival = g
 773
 774 !     STEP I. IMMEDIATE ACCEPTANCE IF ival IS LARGE ENOUGH
 775
 776     IF (ival>=l) RETURN
 777
 778 !     STEP S. SQUEEZE ACCEPTANCE - SAMPLE U
 779
 780     fk = ival
 781     difmuk = mu - fk
 782     CALL RANDOM_NUMBER(u)
 783     IF (d*u >= difmuk*difmuk*difmuk) RETURN
 784   END IF
 785
 786 !     STEP P. PREPARATIONS FOR STEPS Q AND H.
 787 !             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
 788 !             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
 789 !             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
 790 !             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
 791 !             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
 792
 793   IF (.NOT. full_init) THEN
 794     omega = .3989423/s
 795     b1 = .4166667E-1/mu
 796     b2 = .3*b1*b1
 797     c3 = .1428571*b1*b2
 798     c2 = b2 - 15.*c3
 799     c1 = b1 - 6.*b2 + 45.*c3
 800     c0 = 1. - b1 + 3.*b2 - 15.*c3
 801     c = .1069/mu
 802     full_init = .true.
 803   END IF
 804
 805   IF (g < 0.0) GO TO 50
 806
 807 !             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
 808
 809   kflag = 0
 810   GO TO 70
 811
 812 !     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
 813
 814   40 IF (fy-u*fy <= py*EXP(px-fx)) RETURN
 815
 816 !     STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL
 817 !             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
 818 !             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
 819
 820   50 e = random_exponential()
 821   CALL RANDOM_NUMBER(u)
 822   u = u + u - one
 823   t = 1.8 + SIGN(e, u)
 824   IF (t <= (-.6744)) GO TO 50
 825   ival = mu + s*t
 826   fk = ival
 827   difmuk = mu - fk
 828
 829 !             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
 830
 831   kflag = 1
 832   GO TO 70
 833
 834 !     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
 835
 836   60 IF (c*ABS(u) > py*EXP(px+e) - fy*EXP(fx+e)) GO TO 50
 837   RETURN
 838
 839 !     STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY.
 840 !             CASE ival < 10 USES FACTORIALS FROM TABLE FACT
 841
 842   70 IF (ival>=10) GO TO 80
 843   px = -mu
 844   py = mu**ival/fact(ival+1)
 845   GO TO 110
 846
 847 !             CASE ival >= 10 USES POLYNOMIAL APPROXIMATION
 848 !             A0-A7 FOR ACCURACY WHEN ADVISABLE
 849 !             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
 850
 851   80 del = .8333333E-1/fk
 852   del = del - 4.8*del*del*del
 853   v = difmuk/fk
 854   IF (ABS(v)>0.25) THEN
 855     px = fk*LOG(one + v) - difmuk - del
 856   ELSE
 857     px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - del
 858   END IF
 859   py = .3989423/SQRT(fk)
 860   110 x = (half - difmuk)/s
 861   xx = x*x
 862   fx = -half*xx
 863   fy = omega* (((c3*xx + c2)*xx + c1)*xx + c0)
 864   IF (kflag <= 0) GO TO 40
 865   GO TO 60
 866
 867 !---------------------------------------------------------------------------
 868 !     C A S E  B.    mu < 10
 869 !     START NEW TABLE AND CALCULATE P0 IF NECESSARY
 870
 871 ELSE
 872   IF (first) THEN
 873     m = MAX(1, INT(mu))
 874     l = 0
 875     p = EXP(-mu)
 876     q = p
 877     p0 = p
 878   END IF
 879
 880 !     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
 881
 882   DO
 883     CALL RANDOM_NUMBER(u)
 884     ival = 0
 885     IF (u <= p0) RETURN
 886
 887 !     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
 888 !             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
 889 !             (0.458=PP(9) FOR MU=10)
 890
 891     IF (l == 0) GO TO 150
 892     j = 1
 893     IF (u > 0.458) j = MIN(l, m)
 894     DO k = j, l
 895       IF (u <= pp(k)) GO TO 180
 896     END DO
 897     IF (l == 35) CYCLE
 898
 899 !     STEP C. CREATION OF NEW POISSON PROBABILITIES P
 900 !             AND THEIR CUMULATIVES Q=PP(K)
 901
 902     150 l = l + 1
 903     DO k = l, 35
 904       p = p*mu / k
 905       q = q + p
 906       pp(k) = q
 907       IF (u <= q) GO TO 170
 908     END DO
 909     l = 35
 910   END DO
 911
 912   170 l = k
 913   180 ival = k
 914   RETURN
 915 END IF
 916
 917 RETURN
 918 END FUNCTION random_Poisson
 919
 920
 921
 922 FUNCTION random_binomial1(n, p, first) RESULT(ival)
 923
 924 ! FUNCTION GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method.
 925 ! This algorithm is suitable when many random variates are required
 926 ! with the SAME parameter values for n & p.
 927
 928 !    P = BERNOULLI SUCCESS PROBABILITY
 929 !           (0 <= REAL <= 1)
 930 !    N = NUMBER OF BERNOULLI TRIALS
 931 !           (1 <= INTEGER)
 932 !    FIRST = .TRUE. for the first call using the current parameter values
 933 !          = .FALSE. if the values of (n,p) are unchanged from last call
 934
 935 ! Reference: Kemp, C.D. (1986). `A modal method for generating binomial
 936 !            variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.
 937
 938 INTEGER, INTENT(IN) :: n
 939 REAL, INTENT(IN)    :: p
 940 LOGICAL, INTENT(IN) :: first
 941 INTEGER             :: ival
 942
 943 !     Local variables
 944
 945 INTEGER         :: ru, rd
 946 INTEGER, SAVE   :: r0
 947 REAL            :: u, pd, pu
 948 REAL, SAVE      :: odds_ratio, p_r
 949 REAL, PARAMETER :: zero = 0.0, one = 1.0
 950
 951 IF (first) THEN
 952   r0 = (n+1)*p
 953   p_r = bin_prob(n, p, r0)
 954   odds_ratio = p / (one - p)
 955 END IF
 956
 957 CALL RANDOM_NUMBER(u)
 958 u = u - p_r
 959 IF (u < zero) THEN
 960   ival = r0
 961   RETURN
 962 END IF
 963
 964 pu = p_r
 965 ru = r0
 966 pd = p_r
 967 rd = r0
 968 DO
 969   rd = rd - 1
 970   IF (rd >= 0) THEN
 971     pd = pd * (rd+1) / (odds_ratio * (n-rd))
 972     u = u - pd
 973     IF (u < zero) THEN
 974       ival = rd
 975       RETURN
 976     END IF
 977   END IF
 978
 979   ru = ru + 1
 980   IF (ru <= n) THEN
 981     pu = pu * (n-ru+1) * odds_ratio / ru
 982     u = u - pu
 983     IF (u < zero) THEN
 984       ival = ru
 985       RETURN
 986     END IF
 987   END IF
 988 END DO
 989
 990 !     This point should not be reached, but just in case:
 991
 992 ival = r0
 993 RETURN
 994
 995 END FUNCTION random_binomial1
 996
 997
 998
 999 FUNCTION bin_prob(n, p, r) RESULT(fn_val)
1000 !     Calculate a binomial probability
1001
1002 INTEGER, INTENT(IN) :: n, r
1003 REAL, INTENT(IN)    :: p
1004 REAL                :: fn_val
1005
1006 !     Local variable
1007 REAL                :: one = 1.0
1008
1009 fn_val = EXP( lngamma(DBLE(n+1)) - lngamma(DBLE(r+1)) - lngamma(DBLE(n-r+1)) &
1010               + r*LOG(p) + (n-r)*LOG(one - p) )
1011 RETURN
1012
1013 END FUNCTION bin_prob
1014
1015
1016
1017 FUNCTION lngamma(x) RESULT(fn_val)
1018
1019 ! Logarithm to base e of the gamma function.
1020 !
1021 ! Accurate to about 1.e-14.
1022 ! Programmer: Alan Miller
1023
1024 ! Latest revision of Fortran 77 version - 28 February 1988
1025
1026 REAL (dp), INTENT(IN) :: x
1027 REAL (dp)             :: fn_val
1028
1029 !       Local variables
1030
1031 REAL (dp) :: a1 = -4.166666666554424D-02, a2 = 2.430554511376954D-03,  &
1032              a3 = -7.685928044064347D-04, a4 = 5.660478426014386D-04,  &
1033              temp, arg, product, lnrt2pi = 9.189385332046727D-1,       &
1034              pi = 3.141592653589793D0
1035 LOGICAL   :: reflect
1036
1037 !       lngamma is not defined if x = 0 or a negative integer.
1038
1039 IF (x > 0.d0) GO TO 10
1040 IF (x /= INT(x)) GO TO 10
1041 fn_val = 0.d0
1042 RETURN
1043
1044 !       If x < 0, use the reflection formula:
1045 !               gamma(x) * gamma(1-x) = pi * cosec(pi.x)
1046
1047 10 reflect = (x < 0.d0)
1048 IF (reflect) THEN
1049   arg = 1.d0 - x
1050 ELSE
1051   arg = x
1052 END IF
1053
1054 !       Increase the argument, if necessary, to make it > 10.
1055
1056 product = 1.d0
1057 20 IF (arg <= 10.d0) THEN
1058   product = product * arg
1059   arg = arg + 1.d0
1060   GO TO 20
1061 END IF
1062
1063 !  Use a polynomial approximation to Stirling's formula.
1064 !  N.B. The real Stirling's formula is used here, not the simpler, but less
1065 !       accurate formula given by De Moivre in a letter to Stirling, which
1066 !       is the one usually quoted.
1067
1068 arg = arg - 0.5D0
1069 temp = 1.d0/arg**2
1070 fn_val = lnrt2pi + arg * (LOG(arg) - 1.d0 + &
1071                   (((a4*temp + a3)*temp + a2)*temp + a1)*temp) - LOG(product)
1072 IF (reflect) THEN
1073   temp = SIN(pi * x)
1074   fn_val = LOG(pi/temp) - fn_val
1075 END IF
1076 RETURN
1077 END FUNCTION lngamma
1078
1079
1080
1081 FUNCTION random_binomial2(n, pp, first) RESULT(ival)
1082 !**********************************************************************
1083 !     Translated to Fortran 90 by Alan Miller from:
1084 !                              RANLIB
1085 !
1086 !     Library of Fortran Routines for Random Number Generation
1087 !
1088 !                      Compiled and Written by:
1089 !
1090 !                           Barry W. Brown
1091 !                            James Lovato
1092 !
1093 !               Department of Biomathematics, Box 237
1094 !               The University of Texas, M.D. Anderson Cancer Center
1095 !               1515 Holcombe Boulevard
1096 !               Houston, TX      77030
1097 !
1098 ! This work was supported by grant CA-16672 from the National Cancer Institute.
1099
1100 !                    GENerate BINomial random deviate
1101
1102 !                              Function
1103
1104 !     Generates a single random deviate from a binomial
1105 !     distribution whose number of trials is N and whose
1106 !     probability of an event in each trial is P.
1107
1108 !                              Arguments
1109
1110 !     N  --> The number of trials in the binomial distribution
1111 !            from which a random deviate is to be generated.
1112 !                              INTEGER N
1113
1114 !     P  --> The probability of an event in each trial of the
1115 !            binomial distribution from which a random deviate
1116 !            is to be generated.
1117 !                              REAL P
1118
1119 !     FIRST --> Set FIRST = .TRUE. for the first call to perform initialization
1120 !               the set FIRST = .FALSE. for further calls using the same pair
1121 !               of parameter values (N, P).
1122 !                              LOGICAL FIRST
1123
1124 !     random_binomial2 <-- A random deviate yielding the number of events
1125 !                from N independent trials, each of which has
1126 !                a probability of event P.
1127 !                              INTEGER random_binomial
1128
1129 !                              Method
1130
1131 !     This is algorithm BTPE from:
1132
1133 !         Kachitvichyanukul, V. and Schmeiser, B. W.
1134 !         Binomial Random Variate Generation.
1135 !         Communications of the ACM, 31, 2 (February, 1988) 216.
1136
1137 !**********************************************************************
1138
1139 !*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
1140
1141 !     ..
1142 !     .. Scalar Arguments ..
1143 REAL, INTENT(IN)    :: pp
1144 INTEGER, INTENT(IN) :: n
1145 LOGICAL, INTENT(IN) :: first
1146 INTEGER             :: ival
1147 !     ..
1148 !     .. Local Scalars ..
1149 REAL            :: alv, amaxp, f, f1, f2, u, v, w, w2, x, x1, x2, ynorm, z, z2
1150 REAL, PARAMETER :: zero = 0.0, half = 0.5, one = 1.0
1151 INTEGER         :: i, ix, ix1, k, mp
1152 INTEGER, SAVE   :: m
1153 REAL, SAVE      :: p, q, xnp, ffm, fm, xnpq, p1, xm, xl, xr, c, al, xll,  &
1154                    xlr, p2, p3, p4, qn, r, g
1155
1156 !     ..
1157 !     .. Executable Statements ..
1158
1159 !*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
1160
1161 IF (first) THEN
1162   p = MIN(pp, one-pp)
1163   q = one - p
1164   xnp = n * p
1165 END IF
1166
1167 IF (xnp > 30.) THEN
1168   IF (first) THEN
1169     ffm = xnp + p
1170     m = ffm
1171     fm = m
1172     xnpq = xnp * q
1173     p1 = INT(2.195*SQRT(xnpq) - 4.6*q) + half
1174     xm = fm + half
1175     xl = xm - p1
1176     xr = xm + p1
1177     c = 0.134 + 20.5 / (15.3 + fm)
1178     al = (ffm-xl) / (ffm - xl*p)
1179     xll = al * (one + half*al)
1180     al = (xr - ffm) / (xr*q)
1181     xlr = al * (one + half*al)
1182     p2 = p1 * (one + c + c)
1183     p3 = p2 + c / xll
1184     p4 = p3 + c / xlr
1185   END IF
1186
1187 !*****GENERATE VARIATE, Binomial mean at least 30.
1188
1189   20 CALL RANDOM_NUMBER(u)
1190   u = u * p4
1191   CALL RANDOM_NUMBER(v)
1192
1193 !     TRIANGULAR REGION
1194
1195   IF (u <= p1) THEN
1196     ix = xm - p1 * v + u
1197     GO TO 110
1198   END IF
1199
1200 !     PARALLELOGRAM REGION
1201
1202   IF (u <= p2) THEN
1203     x = xl + (u-p1) / c
1204     v = v * c + one - ABS(xm-x) / p1
1205     IF (v > one .OR. v <= zero) GO TO 20
1206     ix = x
1207   ELSE
1208
1209 !     LEFT TAIL
1210
1211     IF (u <= p3) THEN
1212       ix = xl + LOG(v) / xll
1213       IF (ix < 0) GO TO 20
1214       v = v * (u-p2) * xll
1215     ELSE
1216
1217 !     RIGHT TAIL
1218
1219       ix = xr - LOG(v) / xlr
1220       IF (ix > n) GO TO 20
1221       v = v * (u-p3) * xlr
1222     END IF
1223   END IF
1224
1225 !*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
1226
1227   k = ABS(ix-m)
1228   IF (k <= 20 .OR. k >= xnpq/2-1) THEN
1229
1230 !     EXPLICIT EVALUATION
1231
1232     f = one
1233     r = p / q
1234     g = (n+1) * r
1235     IF (m < ix) THEN
1236       mp = m + 1
1237       DO i = mp, ix
1238         f = f * (g/i-r)
1239       END DO
1240
1241     ELSE IF (m > ix) THEN
1242       ix1 = ix + 1
1243       DO i = ix1, m
1244         f = f / (g/i-r)
1245       END DO
1246     END IF
1247
1248     IF (v > f) THEN
1249       GO TO 20
1250     ELSE
1251       GO TO 110
1252     END IF
1253   END IF
1254
1255 !     SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X))
1256
1257   amaxp = (k/xnpq) * ((k*(k/3. + .625) + .1666666666666)/xnpq + half)
1258   ynorm = -k * k / (2.*xnpq)
1259   alv = LOG(v)
1260   IF (alv<ynorm - amaxp) GO TO 110
1261   IF (alv>ynorm + amaxp) GO TO 20
1262
1263 !     STIRLING'S (actually de Moivre's) FORMULA TO MACHINE ACCURACY FOR
1264 !     THE FINAL ACCEPTANCE/REJECTION TEST
1265
1266   x1 = ix + 1
1267   f1 = fm + one
1268   z = n + 1 - fm
1269   w = n - ix + one
1270   z2 = z * z
1271   x2 = x1 * x1
1272   f2 = f1 * f1
1273   w2 = w * w
1274   IF (alv - (xm*LOG(f1/x1) + (n-m+half)*LOG(z/w) + (ix-m)*LOG(w*p/(x1*q)) +    &
1275       (13860.-(462.-(132.-(99.-140./f2)/f2)/f2)/f2)/f1/166320. +               &
1276       (13860.-(462.-(132.-(99.-140./z2)/z2)/z2)/z2)/z/166320. +                &
1277       (13860.-(462.-(132.-(99.-140./x2)/x2)/x2)/x2)/x1/166320. +               &
1278       (13860.-(462.-(132.-(99.-140./w2)/w2)/w2)/w2)/w/166320.) > zero) THEN
1279     GO TO 20
1280   ELSE
1281     GO TO 110
1282   END IF
1283
1284 ELSE
1285 !     INVERSE CDF LOGIC FOR MEAN LESS THAN 30
1286   IF (first) THEN
1287     qn = q ** n
1288     r = p / q
1289     g = r * (n+1)
1290   END IF
1291
1292   90 ix = 0
1293   f = qn
1294   CALL RANDOM_NUMBER(u)
1295   100 IF (u >= f) THEN
1296     IF (ix > 110) GO TO 90
1297     u = u - f
1298     ix = ix + 1
1299     f = f * (g/ix - r)
1300     GO TO 100
1301   END IF
1302 END IF
1303
1304 110 IF (pp > half) ix = n - ix
1305 ival = ix
1306 RETURN
1307
1308 END FUNCTION random_binomial2
1309
1310
1311
1312
1313 FUNCTION random_neg_binomial(sk, p) RESULT(ival)
1314
1315 ! Adapted from Fortran 77 code from the book:
1316 !     Dagpunar, J. 'Principles of random variate generation'
1317 !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
1318
1319 ! FUNCTION GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED
1320 ! INVERSION AND/OR THE REPRODUCTIVE PROPERTY.
1321
1322 !    SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words!)
1323 !       = the `power' parameter of the negative binomial
1324 !           (0 < REAL)
1325 !    P = BERNOULLI SUCCESS PROBABILITY
1326 !           (0 < REAL < 1)
1327
1328 ! THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H,
1329 ! OTHERWISE A COMBINATION OF UNSTORED INVERSION AND
1330 ! THE REPRODUCTIVE PROPERTY IS USED.
1331
1332 REAL, INTENT(IN)   :: sk, p
1333 INTEGER            :: ival
1334
1335 !     Local variables
1336 ! THE PARAMETER ULN = -LOG(MACHINE'S SMALLEST REAL NUMBER).
1337
1338 REAL, PARAMETER    :: h = 0.7
1339 REAL               :: q, x, st, uln, v, r, s, y, g
1340 INTEGER            :: k, i, n
1341
1342 IF (sk <= zero .OR. p <= zero .OR. p >= one) THEN
1343   WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'
1344   STOP
1345 END IF
1346
1347 q = one - p
1348 x = zero
1349 st = sk
1350 IF (p > h) THEN
1351   v = one/LOG(p)
1352   k = st
1353   DO i = 1,k
1354     DO
1355       CALL RANDOM_NUMBER(r)
1356       IF (r > zero) EXIT
1357     END DO
1358     n = v*LOG(r)
1359     x = x + n
1360   END DO
1361   st = st - k
1362 END IF
1363
1364 s = zero
1365 uln = -LOG(vsmall)
1366 IF (st > -uln/LOG(q)) THEN
1367   WRITE(*, *) ' P IS TOO LARGE FOR THIS VALUE OF SK'
1368   STOP
1369 END IF
1370
1371 y = q**st
1372 g = st
1373 CALL RANDOM_NUMBER(r)
1374 DO
1375   IF (y > r) EXIT
1376   r = r - y
1377   s = s + one
1378   y = y*p*g/s
1379   g = g + one
1380 END DO
1381
1382 ival = x + s + half
1383 RETURN
1384 END FUNCTION random_neg_binomial
1385
1386
1387
1388 FUNCTION random_von_Mises(k, first) RESULT(fn_val)
1389
1390 !     Algorithm VMD from:
1391 !     Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a
1392 !     comparison of random numbers', J. of Appl. Statist., 17, 165-168.
1393
1394 !     Fortran 90 code by Alan Miller
1395 !     CSIRO Division of Mathematical & Information Sciences
1396
1397 !     Arguments:
1398 !     k (real)        parameter of the von Mises distribution.
1399 !     first (logical) set to .TRUE. the first time that the function
1400 !                     is called, or the first time with a new value
1401 !                     for k.   When first = .TRUE., the function sets
1402 !                     up starting values and may be very much slower.
1403
1404 REAL, INTENT(IN)     :: k
1405 LOGICAL, INTENT(IN)  :: first
1406 REAL                 :: fn_val
1407
1408 !     Local variables
1409
1410 INTEGER          :: j, n
1411 INTEGER, SAVE    :: nk
1412 REAL, PARAMETER  :: pi = 3.14159265
1413 REAL, SAVE       :: p(20), theta(0:20)
1414 REAL             :: sump, r, th, lambda, rlast
1415 REAL (dp)        :: dk
1416
1417 IF (first) THEN                        ! Initialization, if necessary
1418   IF (k < zero) THEN
1419     WRITE(*, *) '** Error: argument k for random_von_Mises = ', k
1420     RETURN
1421   END IF
1422
1423   nk = k + k + one
1424   IF (nk > 20) THEN
1425     WRITE(*, *) '** Error: argument k for random_von_Mises = ', k
1426     RETURN
1427   END IF
1428
1429   dk = k
1430   theta(0) = zero
1431   IF (k > half) THEN
1432
1433 !     Set up array p of probabilities.
1434
1435     sump = zero
1436     DO j = 1, nk
1437       IF (j < nk) THEN
1438         theta(j) = ACOS(one - j/k)
1439       ELSE
1440         theta(nk) = pi
1441       END IF
1442
1443 !     Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j)
1444
1445       CALL integral(theta(j-1), theta(j), p(j), dk)
1446       sump = sump + p(j)
1447     END DO
1448     p(1:nk) = p(1:nk) / sump
1449   ELSE
1450     p(1) = one
1451     theta(1) = pi
1452   END IF                         ! if k > 0.5
1453 END IF                           ! if first
1454
1455 CALL RANDOM_NUMBER(r)
1456 DO j = 1, nk
1457   r = r - p(j)
1458   IF (r < zero) EXIT
1459 END DO
1460 r = -r/p(j)
1461
1462 DO
1463   th = theta(j-1) + r*(theta(j) - theta(j-1))
1464   lambda = k - j + one - k*COS(th)
1465   n = 1
1466   rlast = lambda
1467
1468   DO
1469     CALL RANDOM_NUMBER(r)
1470     IF (r > rlast) EXIT
1471     n = n + 1
1472     rlast = r
1473   END DO
1474
1475   IF (n .NE. 2*(n/2)) EXIT         ! is n even?
1476   CALL RANDOM_NUMBER(r)
1477 END DO
1478
1479 fn_val = SIGN(th, (r - rlast)/(one - rlast) - half)
1480 RETURN
1481 END FUNCTION random_von_Mises
1482
1483
1484
1485 SUBROUTINE integral(a, b, result, dk)
1486
1487 !     Gaussian integration of exp(k.cosx) from a to b.
1488
1489 REAL (dp), INTENT(IN) :: dk
1490 REAL, INTENT(IN)      :: a, b
1491 REAL, INTENT(OUT)     :: result
1492
1493 !     Local variables
1494
1495 REAL (dp)  :: xmid, range, x1, x2,                                    &
1496   x(3) = (/0.238619186083197_dp, 0.661209386466265_dp, 0.932469514203152_dp/), &
1497   w(3) = (/0.467913934572691_dp, 0.360761573048139_dp, 0.171324492379170_dp/)
1498 INTEGER    :: i
1499
1500 xmid = (a + b)/2._dp
1501 range = (b - a)/2._dp
1502
1503 result = 0._dp
1504 DO i = 1, 3
1505   x1 = xmid + x(i)*range
1506   x2 = xmid - x(i)*range
1507   result = result + w(i)*(EXP(dk*COS(x1)) + EXP(dk*COS(x2)))
1508 END DO
1509
1510 result = result * range
1511 RETURN
1512 END SUBROUTINE integral
1513
1514
1515
1516 FUNCTION random_Cauchy() RESULT(fn_val)
1517
1518 !     Generate a random deviate from the standard Cauchy distribution
1519
1520 REAL     :: fn_val
1521
1522 !     Local variables
1523 REAL     :: v(2)
1524
1525 DO
1526   CALL RANDOM_NUMBER(v)
1527   v = two*(v - half)
1528   IF (ABS(v(2)) < vsmall) CYCLE               ! Test for zero
1529   IF (v(1)**2 + v(2)**2 < one) EXIT
1530 END DO
1531 fn_val = v(1) / v(2)
1532
1533 RETURN
1534 END FUNCTION random_Cauchy
1535
1536
1537
1538 SUBROUTINE random_order(order, n)
1539
1540 !     Generate a random ordering of the integers 1 ... n.
1541
1542 INTEGER, INTENT(IN)  :: n
1543 INTEGER, INTENT(OUT) :: order(n)
1544
1545 !     Local variables
1546
1547 INTEGER :: i, j, k
1548 REAL    :: wk
1549
1550 DO i = 1, n
1551   order(i) = i
1552 END DO
1553
1554 !     Starting at the end, swap the current last indicator with one
1555 !     randomly chosen from those preceeding it.
1556
1557 DO i = n, 2, -1
1558   CALL RANDOM_NUMBER(wk)
1559   j = 1 + i * wk
1560   IF (j < i) THEN
1561     k = order(i)
1562     order(i) = order(j)
1563     order(j) = k
1564   END IF
1565 END DO
1566
1567 RETURN
1568 END SUBROUTINE random_order
1569
1570
1571
1572 SUBROUTINE seed_random_number(iounit)
1573
1574 INTEGER, INTENT(IN)  :: iounit
1575
1576 ! Local variables
1577
1578 INTEGER              :: k
1579 INTEGER, ALLOCATABLE :: seed(:)
1580
1581 CALL RANDOM_SEED(SIZE=k)
1582 ALLOCATE( seed(k) )
1583
1584 WRITE(*, '(a, i2, a)')' Enter ', k, ' integers for random no. seeds: '
1585 READ(*, *) seed
1586 WRITE(iounit, '(a, (7i10))') ' Random no. seeds: ', seed
1587 CALL RANDOM_SEED(PUT=seed)
1588
1589 DEALLOCATE( seed )
1590
1591 RETURN
1592 END SUBROUTINE seed_random_number
1593
1594
1595 END MODULE KPP_ROOT_Random
1596
1597 MODULE KPP_ROOT_Integrator
1598   USE KPP_ROOT_Random
1599   USE KPP_ROOT_Parameters, ONLY : NVAR, NFIX, NREACT
1600   USE KPP_ROOT_Global, ONLY : TIME, RCONST, Volume
1601   USE KPP_ROOT_Stoichiom
1602   USE KPP_ROOT_Stochastic
1603   USE KPP_ROOT_Rates
1604   USE KPP_ROOT_Random, ddp => dp
1605   IMPLICIT NONE
1606
1607 CONTAINS
1608
1609 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1610 SUBROUTINE TauLeap(Nsteps, Tau, T, SCT, NmlcV, NmlcF)
1611 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1612 !
1613 !   Tau-leap stochastic integration
1614 !   INPUT:
1615 !       Nsteps  = no. of tau-leap steps to be simulated
1616 !       Tau     = time step length
1617 !       T       = time
1618 !       SCT     = stochastic rate constants
1619 !       NmlcV, NmlcF = no. of molecules for variable and fixed species
1620 !   OUTPUT:
1621 !       T       = updated time (after Nsteps)
1622 !       NmlcV   = updated no. of molecules for variable species
1623 !
1624 !
1625 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1626   IMPLICIT NONE
1627 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1628
1629       KPP_REAL:: T, Tau
1630       INTEGER :: Nsteps
1631       INTEGER :: NmlcV(NVAR), NmlcF(NFIX)
1632       INTEGER :: i, j, irct, id, istep, Nfirings(NREACT)
1633       REAL   :: mu
1634       KPP_REAL :: A(NREACT), SCT(NREACT), x
1635       LOGICAL, SAVE :: First = .TRUE.
1636
1637       DO istep = 1, Nsteps
1638
1639           ! Propensity vector
1640           CALL  Propensity ( NmlcV, NmlcF, SCT, A )
1641
1642           ! Index of next reaction
1643           DO irct = 1, NREACT
1644             mu = A(irct)*Tau
1645             Nfirings(irct) = random_Poisson(mu, First)
1646             First = .TRUE.
1647           END DO
1648
1649           ! Update time with the leap interval
1650           T = T + Tau;
1651
1652           ! Directly update state vector
1653           DO irct = 1, NREACT
1654              DO i = CCOL_STOICM(irct), CCOL_STOICM(irct+1)-1
1655                 id = IROW_STOICM(i)
1656                 NmlcV(id) = MAX(0, NmlcV(id) + Nfirings(irct)*INT(STOICM(i)))
1657              END DO
1658           END DO
1659
1660           ! Update state vector
1661           ! DO irct = 1, NREACT
1662           !   DO j = 1, Nfirings(irct)
1663           !    CALL MoleculeChange( irct, NmlcV )
1664           !   END DO
1665           ! END DO
1666
1667      END DO
1668
1669 CONTAINS
1670
1671      SUBROUTINE PropensityTemplate( T, NmlcV, NmlcF, Prop )
1672       KPP_REAL, INTENT(IN)  :: T
1673       INTEGER, INTENT(IN)   :: NmlcV(NVAR), NmlcF(NFIX)
1674       KPP_REAL, INTENT(OUT) :: Prop(NREACT)
1675       KPP_REAL :: Tsave
1676 ! Update the stochastic reaction rates, which may be time dependent
1677       Tsave = TIME
1678       TIME = T
1679       CALL Update_RCONST()
1680       CALL StochasticRates( RCONST, Volume, SCT )
1681       CALL Propensity ( NmlcV, NmlcF, SCT, Prop )
1682       TIME = Tsave
1683      END SUBROUTINE PropensityTemplate
1684
1685 END SUBROUTINE TauLeap
1686 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1687
1688 END MODULE KPP_ROOT_Integrator