1 /* $Id: thrdrsmp.c,v 1.4 2000/12/21 14:14:41 beloslyu Exp $
2 *===========================================================================
3 *
4 * PUBLIC DOMAIN NOTICE
5 * National Center for Biotechnology Information
6 *
7 * This software/database is a "United States Government Work" under the
8 * terms of the United States Copyright Act. It was written as part of
9 * the author's official duties as a United States Government employee and
10 * thus cannot be copyrighted. This software/database is freely available
11 * to the public for use. The National Library of Medicine and the U.S.
12 * Government have not placed any restriction on its use or reproduction.
13 *
14 * Although all reasonable efforts have been taken to ensure the accuracy
15 * and reliability of the software and data, the NLM and the U.S.
16 * Government do not and cannot warrant the performance or results that
17 * may be obtained by using this software or data. The NLM and the U.S.
18 * Government disclaim all warranties, express or implied, including
19 * warranties of performance, merchantability or fitness for any particular
20 * purpose.
21 *
22 * Please cite the author in any work or product based on this material.
23 *
24 * ===========================================================================
25 *
26 * File Name: thrdrsmp.c
27 *
28 * Author: Stephen Bryant
29 *
30 * Initial Version Creation Date: 08/16/2000
31 *
32 * $Revision: 1.4 $
33 *
34 * File Description: threader
35 *
36 * Modifications:
37 * --------------------------------------------------------------------------
38 * $Log: thrdrsmp.c,v $
39 * Revision 1.4 2000/12/21 14:14:41 beloslyu
40 * c++ comments are not allowed in c code
41 *
42 * Revision 1.3 2000/12/20 18:56:41 hurwitz
43 * new random num gen, more debug printing
44 *
45 * Revision 1.2 2000/08/16 21:18:57 hurwitz
46 * fix dangerous warnings found by MS Visual C++, replace rand48 functions with toolkit functions
47 *
48 * Revision 1.1 2000/08/16 20:45:21 hurwitz
49 * initial check in of threading routines
50 *
51 * ==========================================================================
52 */
53
54
55
56 /* Choose randomly from a multinomial distribution */
57
58 #include <thrdatd.h>
59 #include <thrddecl.h>
60 #include <ncbimath.h>
61 #include <stdio.h>
62 #include <math.h>
63
rsmp(Rnd_Smp * pvl)64 int rsmp(Rnd_Smp* pvl) {
65 /*----------------------------------------------------*/
66 /* pvl: Number and probabilities of parameter values */
67 /*----------------------------------------------------*/
68
69 int i; /* The i-th offset parameter value will be the choice */
70 float c; /* Cumulative probabilites across parameter values */
71 float r; /* A uniform random number on the interval 0 - 1 */
72 int idum = 1;
73
74 /* for(i=0;i<pvl->n;i++) printf("%.4f ",pvl->p[i]); printf("pvl->p\n"); */
75
76 /* r=drand48(); c=0.; */
77 r = Rand01(&idum);
78 c = 0.0;
79
80 /* printf("r: %.4f\n",r); */
81
82 for(i=0; i<pvl->n; i++) {
83 c=c+pvl->p[i];
84 /* printf("c: %.4f\n",c); */
85 if(c>=r) return(i); }
86
87 return(i-1);
88 }
89
90
91 #define M1 259200
92 #define IA1 7141
93 #define IC1 54773
94 #define RM1 (1.0/M1)
95 #define M2 134456
96 #define IA2 8121
97 #define IC2 28411
98 #define RM2 (1.0/M2)
99 #define M3 243000
100 #define IA3 4561
101 #define IC3 51349
102
Rand01(int * idum)103 float Rand01(int* idum) {
104 /*----------------------------------------------------------------------------
105 // Returns a uniform deviate between 0.0 and 1.0.
106 // Set idum to any negative value to initialize or reinitialize the sequence.
107 //
108 // This is copied from Numerical Recipes, chapter 7.1, 1988.
109 // (ISBN 0-521-35465-X)
110 //--------------------------------------------------------------------------*/
111 static long ix1,ix2,ix3;
112 static float r[98];
113 float temp;
114 static int iff=0;
115 int j;
116
117 if (*idum<0 || iff==0) {
118 iff=1;
119 ix1=(IC1-*idum) % M1;
120 ix1=(IA1*ix1+IC1) % M1;
121 ix2=ix1 % M2;
122 ix1=(IA1*ix1+IC1) % M1;
123 ix3=ix1 % M3;
124 for (j=1; j<=97; j++) {
125 ix1=(IA1*ix1+IC1) % M1;
126 ix2=(IA2*ix2+IC2) % M2;
127 r[j]=(ix1+ix2*RM2)*RM1;
128 }
129 *idum=1;
130 }
131 ix1=(IA1*ix1+IC1) % M1;
132 ix2=(IA2*ix2+IC2) % M2;
133 ix3=(IA3*ix3+IC3) % M3;
134 j=1 + ((97*ix3)/M3);
135 temp=r[j];
136 r[j]=(ix1+ix2*RM2)*RM1;
137 return(temp);
138 }
139