1 /******************************************************************************
2
3 # # ## ##### # # # # ##### ####
4 ## ## # # # # # # # # # # #
5 # ## # # # # ###### # # ##### #
6 # # ###### # # # # # # # ### #
7 # # # # # # # # # # # ### # #
8 # # # # # # # ###### # ##### ### ####
9
10 ******************************************************************************/
11 /* This file is part of MAPMAKER 3.0b, Copyright 1987-1992, Whitehead Institute
12 for Biomedical Research. All rights reserved. See READ.ME for license. */
13
14 #define INC_MATH
15 #define INC_MEM
16 #define INC_MSG
17 #define INC_STR
18 #define INC_IO
19
20 #define INC_MISC
21 #define INC_HELP_DEFS
22 #include "system.h"
23
sq(r)24 real sq(r) real r; { return(r*r); }
rmaxf(r,s)25 real rmaxf(r,s) real r,s; { return(r>=s ? r : s); }
rminf(r,s)26 real rminf(r,s) real r,s; { return(r<=s ? r : s); }
27
lpow2(i)28 long lpow2(i) /* Return 2 to the i-th power for 0<=i<=LONGBITS (31) */
29 int i;
30 {
31 long result; int j;
32
33 if (i<0 || i>LONGBITS) send(CRASH);
34 for(result=(long)1, j=0; j<i; j++) result <<= 1;
35 return(result);
36 }
37
ipow2(i)38 int ipow2(i) /* Return 2 to the i-th power 0<=i<=INTBITS (15) */
39 int i;
40 {
41 int result; int j;
42
43 if (i<0 || i>INTBITS) send(CRASH);
44 for(result=1, j=0; j<i; j++) result <<= 1;
45 return(result);
46 }
47
48
lpow(x,i)49 long lpow(x,i) /* Return x to the i-th power */
50 int x,i;
51 {
52 long result; int j;
53
54 if (i<0 || i>LONGBITS) send(CRASH);
55 for(result=(long)1, j=0; j<i; j++) result*=(long)x;
56 return(result);
57 }
58
ipow(x,i)59 int ipow(x,i) /* Return 2 to the i-th power 0<=i<=INTBITS (15) */
60 int i;
61 {
62 int result; int j;
63
64 if (i<0 || i>INTBITS) send(CRASH);
65 for(result=1, j=0; j<i; j++) result*=x;
66 return(result);
67 }
68
ichoose(n,k)69 int ichoose(n,k) /* Best algorithm for small k */
70 int n, k;
71 {
72 long a, b, Ln, Lk, L1, Lx;
73
74 if (k>n || n<1 || k<1) send(CRASH);
75 Ln= (long)n; Lk= (long)k; Lx= (long)(n-k+1); L1= (long)1;
76
77 for (a=L1; Ln>=Lx; Ln--) a*=Ln; /* is product of n...n-k+1, or n!/(n-k)! */
78 for (b=L1; Lk>L1; Lk--) b*=Lk;
79 return((int)(a/b));
80 }
81
imaxf(r,s)82 int imaxf(r,s) int r,s; { return(r>=s ? r : s); }
iminf(r,s)83 int iminf(r,s) int r,s; { return(r<=s ? r : s); }
lmaxf(r,s)84 long lmaxf(r,s) long r,s; { return(r>=s ? r : s); }
lminf(r,s)85 long lminf(r,s) long r,s; { return(r<=s ? r : s); }
86
87
88 /*** Compare Functions for Array Sort Routines ***/
89
90 int icomp(x,y)
91 QSORT_COMPARE_PTR_TO(int) *x, *y;
92 { return( (*(int*)y) - (*(int*)x) ); }
93
94 int lcomp(x,y)
95 QSORT_COMPARE_PTR_TO(long) *x, *y;
96 { long a, b; a= *(long*)x; b= *(long*)y;
97 if (a<b) return(-1); else if (a==b) return(0); else return(1); }
98
99 int rcomp(x,y)
100 QSORT_COMPARE_PTR_TO(real) *x, *y;
101 { real a, b; a= *(real*)x; b= *(real*)y;
102 if (a<b) return(-1); else if (a==b) return(0); else return(1); }
103
104 int scomp(x,y)
105 QSORT_COMPARE_PTR_TO(char) **x, **y; /* does this work? */
106 { return((int)strcmp(*(char**)x,*(char**)y)); }
107
108 int inv_icomp(x,y)
109 QSORT_COMPARE_PTR_TO(int) *x, *y;
110 { return( (*(int*)x) - (*(int*)y) ); }
111
112 int inv_rcomp(x,y)
113 QSORT_COMPARE_PTR_TO(real) *x, *y;
114 { real a, b; a= *(real*)x; b= *(real*)y;
115 if (a>b) return(-1); else if (a==b) return(0); else return(1); }
116
117 /* this needs work */
rhistogram(data,length,min_num_buckets,scale_quantization,scale_limit_quantization)118 int rhistogram(data,length,min_num_buckets,
119 scale_quantization,scale_limit_quantization) /* return T or F */
120 real *data;
121 int length, min_num_buckets;
122 real scale_quantization; /* for present, ignored */
123 real scale_limit_quantization; /* for present, ignored */
124 {
125 int i, j, index, stars, max_stars;
126 int *bucket, num_buckets, largest_bucket;
127 real max_data, scale, min_data;
128 real stars_per;
129
130 num_buckets= min_num_buckets;
131 for (i=0, max_data= -VERY_BIG, min_data= VERY_BIG; i<length; i++) {
132 if (data[i]>max_data) max_data=data[i];
133 if (data[i]<min_data) min_data=data[i];
134 }
135 scale= (max_data - min_data)/REAL(num_buckets); /* range per bucket */
136
137 run { array(bucket, num_buckets, int); }
138 except_when (NOMEMORY) return(FALSE);
139 for (j=0; j<num_buckets; j++) bucket[j]=0;
140
141 for (i=0; i<length; i++) {
142 if (data[i]==max_data) index=num_buckets-1; else
143 index= (int)((data[i]-min_data)/scale);
144 if (!irange(&index,0,num_buckets-1)) send(CRASH);
145 ++(bucket[index]);
146 }
147
148 max_stars= LINE-10;
149 for (j=0, largest_bucket=0; j<num_buckets; j++)
150 if (bucket[j]>largest_bucket) largest_bucket=bucket[j];
151 if (largest_bucket<=max_stars) stars_per=REAL(max_stars/largest_bucket);
152 else stars_per=REAL(max_stars)/REAL(largest_bucket);
153
154 for (j=0; j<num_buckets; j++) {
155 sf(ps,"%6.3lf | ",((real)j)*scale+(scale/2.0)+min_data); print(ps);
156 stars= INT(REAL(bucket[j]) * stars_per);
157 irange(&stars,0,max_stars);
158 if (bucket[j]>0 && stars<1) print(":");
159 else for (i=0; i<INT(stars); i++) print("*");
160 nl();
161 }
162 unarray(bucket, int);
163 return(TRUE);
164 }
165
166
rmean(data,length)167 real rmean(data,length)
168 real *data;
169 int length;
170 {
171 int i;
172 real total;
173
174 if (length==0) return(0.0);
175 for (i=0, total=0.0; i<length; i++) total+= data[i];
176 return(total/REAL(length));
177 }
178
179
rmaxin(data,length)180 real rmaxin(data,length)
181 real *data;
182 int length;
183 {
184 int i;
185 real largest;
186
187 if (length==0) return(0.0);
188 for (i=0, largest= -VERY_BIG; i<length; i++)
189 if (data[i]>largest) largest=data[i];
190 return(largest);
191 }
192
193
rmedian(data,length)194 real rmedian(data,length)
195 real *data;
196 int length;
197 {
198 real *copy, median;
199
200 if (length==0) return(0.0);
201 array(copy, length, real);
202 rcopy(copy, data, length);
203 rsort(copy,length);
204 median= rmiddle(copy,length);
205 unarray(copy, real);
206 return(median);
207 }
208
209
rmiddle(data,length)210 real rmiddle(data,length)
211 real *data;
212 int length;
213 {
214 if (length==0) return(0.0);
215 if (length<3) return(data[0]);
216 return(data[(length-1)/2]);
217 }
218
219
rcopy(to,from,length)220 void rcopy(to,from,length)
221 real *to, *from;
222 int length;
223 { int i; for (i=0; i<length; i++) to[i]=from[i]; }
224
225
dummy_math_calls()226 void dummy_math_calls()
227 {
228 real x,y,z; x=y=z=0.5;
229
230 x=log10(y);
231 x=pow(y,z);
232 x=log(y);
233 x=exp(y);
234
235 x=floor(y);
236 x=ceil(y);
237 x=fabs(y);
238
239 x=sin(y);
240 x=cos(y);
241 x=tan(y);
242 x=sinh(y);
243 x=cosh(y);
244 x=tanh(y);
245 x=asin(y);
246 x=acos(y);
247 x=atan(y);
248 }
249
250
251 real two_sigma_sq;
252
normal_func(deviation)253 real normal_func(deviation) /* return P(deviation) */
254 real deviation;
255 { return(exp(-sq(deviation)/two_sigma_sq)); }
256
make_normal_dist(mu,sigma,inc,limit)257 DISTRIBUTION *make_normal_dist(mu,sigma,inc,limit)
258 real mu, sigma, inc, limit; /* inc and limit are fractions of sigma */
259 {
260 two_sigma_sq= 2.0 * sq(sigma);
261 return(make_distribution(inc*sigma,inc,limit*sigma,mu,normal_func));
262 }
263
264
265
make_distribution(sigma,inc,limit,mean,prob_func)266 DISTRIBUTION *make_distribution(sigma,inc,limit,mean,prob_func)
267 real sigma, inc, limit, mean; /* inc and limit here are offsets from mean */
268 real (*prob_func)(); /* a ptr to a function */
269 {
270 real d, total_prob, cum_prob, end;
271 real *prob, *deviation;
272 int length, i, j;
273 DISTRIBUTION *dist;
274
275 prob=NULL; deviation=NULL; dist=NULL;
276 run {
277 end= limit + VERY_SMALL;
278 length= ((int)(2.0*limit/inc + 0.99999)) + 1;
279
280 array(prob, length, real);
281 array(deviation, length, real);
282 single(dist, DISTRIBUTION);
283
284 total_prob= 0.0;
285 for (d= -limit, i=0; d<end; d+=inc, i++) {
286 if (i>=length) send(CRASH);
287 deviation[i]= d+mean;
288 total_prob+= prob[i]= (*prob_func)(d);
289 }
290
291 cum_prob= 0.0;
292 for (j=0; j<i; j++) {
293 cum_prob+= prob[j];
294 prob[j]= cum_prob/total_prob;
295 }
296
297 dist->entries= i;
298 dist->start= deviation[0];
299 dist->increment= inc;
300 dist->mean= mean;
301 dist->deviation= deviation;
302 dist->prob= prob;
303
304 } when_aborting {
305 unsingle(dist, DISTRIBUTION);
306 unarray(prob, real);
307 unarray(deviation, real);
308 relay;
309 return((DISTRIBUTION*)NULL); /* never reached */
310 }
311 return(dist);
312 }
313
314
pick_from_distribution(dist,prob)315 real pick_from_distribution(dist,prob)
316 DISTRIBUTION *dist;
317 real *prob; /* may be NULL */
318 {
319 real p, *probs;
320 int i, last;
321
322 p= randnum();
323 probs= dist->prob; last= dist->entries-1;
324 for (i=0; probs[i]<p; i++) if (i==last) break;
325 if (prob!=NULL) *prob=p;
326 return(dist->deviation[i]);
327 }
328
329
330 void eliminate(); /* defined below */
331
mat_invert(m,size,m_inverse)332 void mat_invert(m,size,m_inverse) /* Invert square matrix by Gauss' method */
333 real **m;
334 int size;
335 real **m_inverse;
336 /* m_inverse should be 2*size columns (2nd index) by size rows (first index)
337 its left square will be left with the result (ignore the right side) */
338 {
339 int row, col, c, twice_size;
340 real diff, value;
341
342 twice_size= 2 * size;
343 for (row=0; row<size; row++)
344 for (col=0; col<size; col++) {
345 m_inverse[row][col]= m[row][col];
346 m_inverse[row][col+size]= (row==col ? 1.0 : 0.0);
347 }
348
349 for (col=0; col<size; col++) {
350 if (m_inverse[col][col]==0.0) {
351 for (row=0; row<size && m_inverse[row][col]==0.0; row++);
352 if (row==size) send(CRASH);
353 for (c=0; c<twice_size; c++)
354 m_inverse[col][c]+= m_inverse[row][c];
355 }
356 for (row=0; row<size; row++)
357 if (row!=col) eliminate(row,col,col,m_inverse,twice_size);
358 }
359
360 for (row=0; row<size; row++) {
361 value= m_inverse[row][row];
362 for (col=0; col<twice_size; col++)
363 m_inverse[row][col]/= value;
364 }
365
366 for (row=0; row<size; row++)
367 for (col=0; col<size; col++)
368 m_inverse[row][col]= m_inverse[row][col+size];
369
370 /* if ((diff=test_mult(m,m_inverse,size))>mat_tolerance)
371 { MATINV_diff= diff; send(MATINV); } OLD */
372 }
373
374
eliminate(row,col,row_to_sub,m,n_cols)375 void eliminate(row,col,row_to_sub,m,n_cols)
376 int row, col, row_to_sub;
377 real **m;
378 int n_cols;
379 {
380 int j;
381 real value;
382
383 /* Get m[row][col]=0 by subtracting the row with a 1 in col
384 from m[row], where the row is scaled appropriately */
385
386 if (m[row][col]==0.0) return;
387 if (m[row_to_sub][col]==0.0) send(CRASH); /* used to send(SINGMAT) */
388 /* but SINGMAT was undefined */
389 value= m[row][col]/m[row_to_sub][col];
390 for (j=0; j<n_cols; j++) m[row_to_sub][j]*= value;
391 for (j=0; j<n_cols; j++) m[row][j]-= (m[row_to_sub][j]);
392 }
393
394
mat_mult(m,m2,size,result)395 void mat_mult(m,m2,size,result) /* NEEDS TESTING */
396 real **m, **m2;
397 int size;
398 real **result;
399 {
400 int i, j, k;
401 real entry;
402
403 for (i=0; i<size; i++)
404 for (j=0; j<size; j++) {
405 for (k=0, entry= 0.0; k<size; k++) entry+= m[i][k] * m2[k][j];
406 result[i][j]= entry;
407 }
408 }
409
410
array_times_matrix(a,b,rows,columns,c)411 void array_times_matrix(a,b,rows,columns,c)
412 real *a, **b;
413 int rows,columns;
414 real *c;
415 {
416 int i, j;
417
418 for (i=0; i <columns; i++)
419 for (j=0, c[i] = 0.0; j < rows; j++) c[i] += a[j] * b[j][i];
420 }
421
422
423
424