1 /*
2 * This file is a collection (more can be added) of wrappers around some
3 * CDFLIB Fortran algorithms.
4 */
5
6 /*
7 * Notice q and p are used in reverse from their meanings in distributions.py
8 */
9
10 #include "cdf_wrappers.h"
11
12 #if defined(NO_APPEND_FORTRAN)
13 #if defined(UPPERCASE_FORTRAN)
14 #define F_FUNC(f,F) F
15 #else
16 #define F_FUNC(f,F) f
17 #endif
18 #else
19 #if defined(UPPERCASE_FORTRAN)
20 #define F_FUNC(f,F) F##_
21 #else
22 #define F_FUNC(f,F) f##_
23 #endif
24 #endif
25
26 /*
27 * Call macros for different numbers of distribution parameters.
28 */
29
30 #define CDFLIB_CALL2(func, name, a, b, result, return_bound) \
31 if (npy_isnan(p) || npy_isnan(q) || npy_isnan(a) || npy_isnan(b) || \
32 npy_isnan(bound)) { \
33 return NPY_NAN; \
34 } \
35 func(&which, &p, &q, &a, &b, &status, &bound); \
36 return get_result(name, status, bound, result, return_bound)
37
38 #define CDFLIB_CALL3(func, name, a, b, c, result, return_bound) \
39 if (npy_isnan(p) || npy_isnan(q) || npy_isnan(a) || npy_isnan(b) || \
40 npy_isnan(c) || npy_isnan(bound)) { \
41 return NPY_NAN; \
42 } \
43 func(&which, &p, &q, &a, &b, &c, &status, &bound); \
44 return get_result(name, status, bound, result, return_bound)
45
46 #define CDFLIB_CALL4(func, name, a, b, c, d, result, return_bound) \
47 if (npy_isnan(p) || npy_isnan(q) || npy_isnan(a) || npy_isnan(b) || \
48 npy_isnan(c) || npy_isnan(d) || npy_isnan(bound)) { \
49 return NPY_NAN; \
50 } \
51 func(&which, &p, &q, &a, &b, &c, &d, &status, &bound); \
52 return get_result(name, status, bound, result, return_bound)
53
54
55 /* Return nan on status==1,2 */
56 #define NO_RETURN_BOUND 0
57 /* Return bound on status==1,2 */
58 #define RETURN_BOUND 1
59
60
61 /*
62 * Return status checking function
63 */
64
get_result(char * name,int status,double bound,double result,int return_bound)65 static double get_result(char *name, int status, double bound, double result, int return_bound) {
66 if (status < 0) {
67 sf_error(name, SF_ERROR_ARG, "(Fortran) input parameter %d is out of range", (-status));
68 }
69 else {
70 switch (status) {
71 case 0:
72 /* no error */
73 return result;
74 case 1:
75 sf_error(name, SF_ERROR_OTHER, "Answer appears to be lower than lowest search bound (%g)", bound);
76 if (return_bound) {
77 return bound;
78 }
79 break;
80 case 2:
81 sf_error(name, SF_ERROR_OTHER, "Answer appears to be higher than highest search bound (%g)", bound);
82 if (return_bound) {
83 return bound;
84 }
85 break;
86 case 3:
87 case 4:
88 sf_error(name, SF_ERROR_OTHER, "Two parameters that should sum to 1.0 do not");
89 break;
90 case 10:
91 sf_error(name, SF_ERROR_OTHER, "Computational error");
92 break;
93 default:
94 sf_error(name, SF_ERROR_OTHER, "Unknown error");
95 }
96 }
97 return NPY_NAN;
98 }
99
100
101 extern void F_FUNC(cdfbet,CDFBET)(int*,double*,double*,double*,double*,double*,double*,int*,double*);
102
cdfbet3_wrap(double p,double b,double x)103 double cdfbet3_wrap(double p, double b, double x) {
104 int which=3;
105 double q=1.0-p, y=1.0-x, a=0, bound=0;
106 int status=10;
107
108 CDFLIB_CALL4(F_FUNC(cdfbet,CDFBET), "btdtria",
109 x, y, a, b,
110 a, RETURN_BOUND);
111 }
112
cdfbet4_wrap(double a,double p,double x)113 double cdfbet4_wrap(double a, double p, double x) {
114 int which=4;
115 double q=1.0-p, y=1.0-x, b=0, bound=0;
116 int status=10;
117
118 CDFLIB_CALL4(F_FUNC(cdfbet,CDFBET), "btdtrib",
119 x, y, a, b,
120 b, RETURN_BOUND);
121 }
122
123
124 extern void F_FUNC(cdfbin,CDFBIN)(int*,double*,double*,double*,double*,double*,double*,int*,double*);
125
cdfbin2_wrap(double p,double xn,double pr)126 double cdfbin2_wrap(double p, double xn, double pr) {
127 int which=2;
128 double q=1.0-p, s=0, ompr=1.0-pr, bound=0;
129 int status=10;
130
131 CDFLIB_CALL4(F_FUNC(cdfbin,CDFBIN), "bdtrik",
132 s, xn, pr, ompr,
133 s, RETURN_BOUND);
134 }
135
cdfbin3_wrap(double s,double p,double pr)136 double cdfbin3_wrap(double s, double p, double pr) {
137 int which=3;
138 double q=1.0-p, xn=0, ompr=1.0-pr, bound=0;
139 int status=10;
140
141 CDFLIB_CALL4(F_FUNC(cdfbin,CDFBIN), "bdtrin",
142 s, xn, pr, ompr,
143 xn, RETURN_BOUND);
144 }
145
146 extern void F_FUNC(cdfchi,CDFCHI)(int*,double*,double*,double*,double*,int*,double*);
cdfchi3_wrap(double p,double x)147 double cdfchi3_wrap(double p, double x){
148 int which=3;
149 double q=1.0-p, df=0, bound=0;
150 int status=10;
151
152 CDFLIB_CALL2(F_FUNC(cdfchi,CDFCHI), "chdtriv",
153 x, df,
154 df, RETURN_BOUND);
155 }
156
157 extern void F_FUNC(cdfchn,CDFCHN)(int*,double*,double*,double*,double*,double*,int*,double*);
cdfchn1_wrap(double x,double df,double nc)158 double cdfchn1_wrap(double x, double df, double nc) {
159 int which=1;
160 double q=0, p=0, bound=0;
161 int status=10;
162
163 CDFLIB_CALL3(F_FUNC(cdfchn,CDFCHN), "chndtr",
164 x, df, nc,
165 p, RETURN_BOUND);
166 }
167
cdfchn2_wrap(double p,double df,double nc)168 double cdfchn2_wrap(double p, double df, double nc) {
169 int which=2;
170 double q=1.0-p, x=0, bound=0;
171 int status=10;
172
173 CDFLIB_CALL3(F_FUNC(cdfchn,CDFCHN), "chndtrix",
174 x, df, nc,
175 x, NO_RETURN_BOUND);
176 }
177
cdfchn3_wrap(double x,double p,double nc)178 double cdfchn3_wrap(double x, double p, double nc) {
179 int which=3;
180 double q=1.0-p, df=0, bound=0;
181 int status=10;
182
183 CDFLIB_CALL3(F_FUNC(cdfchn,CDFCHN), "chndtridf",
184 x, df, nc,
185 df, RETURN_BOUND);
186 }
187
cdfchn4_wrap(double x,double df,double p)188 double cdfchn4_wrap(double x, double df, double p) {
189 int which=4;
190 double q=1.0-p, nc=0, bound=0;
191 int status=10;
192
193 CDFLIB_CALL3(F_FUNC(cdfchn,CDFCHN), "chndtrinc",
194 x, df, nc,
195 nc, RETURN_BOUND);
196 }
197
198 extern void F_FUNC(cdff,CDFF)(int*,double*,double*,double*,double*,double*,int*,double*);
199 /*
200 double cdff1_wrap(double dfn, double dfd, double f) {
201 int which=1;
202 double q, p, bound;
203 int status=10;
204
205 CDFLIB_CALL3(F_FUNC(cdff,CDFF), "fdtr",
206 f, dfn, dfd,
207 p, NO_RETURN_BOUND);
208 }
209
210 double cdff2_wrap(double dfn, double dfd, double p) {
211 int which=2;
212 double q=1.0-p, f, bound;
213 int status=10;
214
215 CDFLIB_CALL3(F_FUNC(cdff,CDFF), "fdtri",
216 f, dfn, dfd,
217 f, NO_RETURN_BOUND);
218 }
219 */
220
221 /* This seem to give some trouble. No idea why... */
cdff3_wrap(double p,double dfd,double f)222 double cdff3_wrap(double p, double dfd, double f) {
223 int which=3;
224 double q=1.0-p, dfn=0, bound=0;
225 int status=10;
226
227 CDFLIB_CALL3(F_FUNC(cdff,CDFF), "fdtridfn",
228 f, dfn, dfd,
229 dfn, RETURN_BOUND);
230 }
231
cdff4_wrap(double dfn,double p,double f)232 double cdff4_wrap(double dfn, double p, double f) {
233 int which=4;
234 double q=1.0-p, dfd=0, bound=0;
235 int status=10;
236
237 CDFLIB_CALL3(F_FUNC(cdff,CDFF), "fdtridfd",
238 f, dfn, dfd,
239 dfd, RETURN_BOUND);
240 }
241
242
243 extern void F_FUNC(cdffnc,CDFFNC)(int*,double*,double*,double*,double*,double*,double*,int*,double*);
cdffnc1_wrap(double dfn,double dfd,double nc,double f)244 double cdffnc1_wrap(double dfn, double dfd, double nc, double f) {
245 int which=1;
246 double q=0, p=0, bound=0;
247 int status=10;
248
249 CDFLIB_CALL4(F_FUNC(cdffnc,CDFFNC), "ncfdtr",
250 f, dfn, dfd, nc,
251 p, NO_RETURN_BOUND);
252 }
253
cdffnc2_wrap(double dfn,double dfd,double nc,double p)254 double cdffnc2_wrap(double dfn, double dfd, double nc, double p) {
255 int which=2;
256 double q=1.0-p, f=0, bound=0;
257 int status=10;
258
259 CDFLIB_CALL4(F_FUNC(cdffnc,CDFFNC), "ncfdtri",
260 f, dfn, dfd, nc,
261 f, RETURN_BOUND);
262 }
263
264
cdffnc3_wrap(double p,double dfd,double nc,double f)265 double cdffnc3_wrap(double p, double dfd, double nc, double f) {
266 int which=3;
267 double q=1.0-p, dfn=0, bound=0;
268 int status=10;
269
270 CDFLIB_CALL4(F_FUNC(cdffnc,CDFFNC), "ncfdtridfn",
271 f, dfn, dfd, nc,
272 dfn, RETURN_BOUND);
273 }
274
cdffnc4_wrap(double dfn,double p,double nc,double f)275 double cdffnc4_wrap(double dfn, double p, double nc, double f) {
276 int which=4;
277 double q=1.0-p, dfd=0, bound=0;
278 int status=10;
279
280 CDFLIB_CALL4(F_FUNC(cdffnc,CDFFNC), "ncfdtridfd",
281 f, dfn, dfd, nc,
282 dfd, RETURN_BOUND);
283 }
284
cdffnc5_wrap(double dfn,double dfd,double p,double f)285 double cdffnc5_wrap(double dfn, double dfd, double p, double f) {
286 int which=5;
287 double q=1.0-p, nc=0, bound=0;
288 int status=10;
289
290 CDFLIB_CALL4(F_FUNC(cdffnc,CDFFNC), "ncfdtrinc",
291 f, dfn, dfd, nc,
292 nc, RETURN_BOUND);
293 }
294
295 /* scl == a in gdtr
296 shp == b in gdtr
297 */
298 extern void F_FUNC(cdfgam,CDFGAM)(int*,double*,double*,double*,double*,double*,int*,double*);
cdfgam1_wrap(double scl,double shp,double x)299 double cdfgam1_wrap(double scl, double shp, double x) {
300 int which=1;
301 double q=0, p=0, bound=0;
302 int status=10;
303
304 CDFLIB_CALL3(F_FUNC(cdfgam,CDFGAM), "gdtr",
305 x, shp, scl,
306 p, NO_RETURN_BOUND);
307 }
308
cdfgam2_wrap(double scl,double shp,double p)309 double cdfgam2_wrap(double scl, double shp, double p) {
310 int which=2;
311 double q=1.0-p, x=0, bound=0;
312 int status=10;
313
314 CDFLIB_CALL3(F_FUNC(cdfgam,CDFGAM), "gdtrix",
315 x, shp, scl,
316 x, RETURN_BOUND);
317 }
318
cdfgam3_wrap(double scl,double p,double x)319 double cdfgam3_wrap(double scl, double p, double x) {
320 int which=3;
321 double q=1.0-p, shp=0, bound=0;
322 int status=10;
323
324 CDFLIB_CALL3(F_FUNC(cdfgam,CDFGAM), "gdtrib",
325 x, shp, scl,
326 shp, RETURN_BOUND);
327 }
328
cdfgam4_wrap(double p,double shp,double x)329 double cdfgam4_wrap(double p, double shp, double x) {
330 int which=4;
331 double q=1.0-p, scl=0, bound=0;
332 int status=10;
333
334 CDFLIB_CALL3(F_FUNC(cdfgam,CDFGAM), "gdtria",
335 x, shp, scl,
336 scl, RETURN_BOUND);
337 }
338
339 extern void F_FUNC(cdfnbn,CDFNBN)(int*,double*,double*,double*,double*,double*,double*,int*,double*);
cdfnbn2_wrap(double p,double xn,double pr)340 double cdfnbn2_wrap(double p, double xn, double pr) {
341 int which=2;
342 double q=1.0-p, s=0, ompr=1.0-pr, bound=0;
343 int status=10;
344
345 CDFLIB_CALL4(F_FUNC(cdfnbn,CDFNBN), "nbdtrik",
346 s, xn, pr, ompr,
347 s, RETURN_BOUND);
348 }
349
cdfnbn3_wrap(double s,double p,double pr)350 double cdfnbn3_wrap(double s, double p, double pr) {
351 int which=3;
352 double q=1.0-p, xn=0, ompr=1.0-pr, bound=0;
353 int status=10;
354
355 CDFLIB_CALL4(F_FUNC(cdfnbn,CDFNBN), "nbdtrin",
356 s, xn, pr, ompr,
357 xn, RETURN_BOUND);
358 }
359
360 extern void F_FUNC(cdfnor,CDFNOR)(int*,double*,double*,double*,double*,double*,int*,double*);
cdfnor3_wrap(double p,double std,double x)361 double cdfnor3_wrap(double p, double std, double x) {
362 int which=3;
363 double q=1.0-p, mn=0, bound=0;
364 int status=10;
365
366 CDFLIB_CALL3(F_FUNC(cdfnor,CDFNOR), "nrdtrimn",
367 x, mn, std,
368 mn, RETURN_BOUND);
369 }
370
cdfnor4_wrap(double mn,double p,double x)371 double cdfnor4_wrap(double mn, double p, double x) {
372 int which=4;
373 double q=1.0-p, std=0, bound=0;
374 int status=10;
375
376 CDFLIB_CALL3(F_FUNC(cdfnor,CDFNOR), "nrdtrisd",
377 x, mn, std,
378 std, RETURN_BOUND);
379 }
380
381 extern void F_FUNC(cdfpoi,CDFPOI)(int*,double*,double*,double*,double*,int*,double*);
cdfpoi2_wrap(double p,double xlam)382 double cdfpoi2_wrap(double p, double xlam){
383 int which=2;
384 double q=1.0-p, s=0, bound=0;
385 int status=10;
386
387 CDFLIB_CALL2(F_FUNC(cdfpoi,CDFPOI), "pdtrik",
388 s, xlam,
389 s, RETURN_BOUND);
390 }
391
392 extern void F_FUNC(cdft,CDFT)(int*,double*,double*,double*,double*,int*,double*);
cdft1_wrap(double df,double t)393 double cdft1_wrap(double df, double t){
394 int which=1;
395 double q=0, p=0, bound=0;
396 int status=10;
397
398 CDFLIB_CALL2(F_FUNC(cdft,CDFT), "stdtr",
399 t, df,
400 p, NO_RETURN_BOUND);
401 }
402
cdft2_wrap(double df,double p)403 double cdft2_wrap(double df, double p){
404 int which=2;
405 double q=1.0-p, t=0, bound=0;
406 int status=10;
407
408 CDFLIB_CALL2(F_FUNC(cdft,CDFT), "stdtrit",
409 t, df,
410 t, RETURN_BOUND);
411 }
412
cdft3_wrap(double p,double t)413 double cdft3_wrap(double p, double t){
414 int which=3;
415 double q=1.0-p, df=0, bound=0;
416 int status=10;
417
418 CDFLIB_CALL2(F_FUNC(cdft,CDFT), "stdtridf",
419 t, df,
420 df, RETURN_BOUND);
421 }
422
423 extern void F_FUNC(cdftnc,CDFTNC)(int*,double*,double*,double*,double*,double*,int*,double*);
cdftnc1_wrap(double df,double nc,double t)424 double cdftnc1_wrap(double df, double nc, double t) {
425 int which=1;
426 double q=0, p=0, bound=0;
427 int status=10;
428
429 CDFLIB_CALL3(F_FUNC(cdftnc,CDFTNC), "nctdtr",
430 t, df, nc,
431 p, RETURN_BOUND);
432 }
433
cdftnc2_wrap(double df,double nc,double p)434 double cdftnc2_wrap(double df, double nc, double p) {
435 int which=2;
436 double q=1.0-p, t=0, bound=0;
437 int status=10;
438
439 CDFLIB_CALL3(F_FUNC(cdftnc,CDFTNC), "nctdtrit",
440 t, df, nc,
441 t, RETURN_BOUND);
442 }
443
cdftnc3_wrap(double p,double nc,double t)444 double cdftnc3_wrap(double p, double nc, double t) {
445 int which=3;
446 double q=1.0-p, df=0, bound=0;
447 int status=10;
448
449 CDFLIB_CALL3(F_FUNC(cdftnc,CDFTNC), "nctdtridf",
450 t, df, nc,
451 df, RETURN_BOUND);
452 }
453
cdftnc4_wrap(double df,double p,double t)454 double cdftnc4_wrap(double df, double p, double t) {
455 int which=4;
456 double q=1.0-p, nc=0, bound=0;
457 int status=10;
458
459 CDFLIB_CALL3(F_FUNC(cdftnc,CDFTNC), "nctdtrinc",
460 t, df, nc,
461 nc, RETURN_BOUND);
462 }
463