← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/odr/__odrpack.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/odr/__odrpack.c
 
1
/* Anti-Copyright
 
2
 *
 
3
 * I hereby release this code into the PUBLIC DOMAIN AS IS.  There is no
 
4
 * support, warranty, or guarantee.  I will gladly accept comments, bug 
 
5
 * reports, and patches, however.
 
6
 *
 
7
 * Robert Kern
 
8
 * kern@caltech.edu
 
9
 *
 
10
 */
 
11
 
 
12
#include "odrpack.h"
 
13
 
 
14
 
 
15
void F_FUNC(dodrc,DODRC)(void (*fcn)(int *n, int *m, int *np, int *nq, int *ldn, int *ldm, 
 
16
            int *ldnp, double *beta, double *xplusd, int *ifixb, int *ifixx, 
 
17
            int *ldifx, int *ideval, double *f, double *fjacb, double *fjacd, 
 
18
            int *istop), 
 
19
           int *n, int *m, int *np, int *nq, double *beta, double *y, int *ldy,
 
20
           double *x, int *ldx, double *we, int *ldwe, int *ld2we, double *wd,
 
21
           int *ldwd, int *ld2wd, int *ifixb, int *ifixx, int *ldifx, int *job,
 
22
           int *ndigit, double *taufac, double *sstol, double *partol, 
 
23
           int *maxit, int *iprint, int *lunerr, int *lunrpt, double *stpb,
 
24
           double *stpd, int *ldstpd, double *sclb, double *scld, int *ldscld,
 
25
           double *work, int *lwork, int *iwork, int *liwork, int *info);
 
26
void F_FUNC(dwinf,DWINF)(int *n, int *m, int *np, int *nq, int *ldwe, int *ld2we, int *isodr,
 
27
        int *delta, int *eps, int *xplus, int *fn, int *sd, int *vcv, int *rvar,
 
28
        int *wss, int *wssde, int *wssep, int *rcond, int *eta, int *olmav, 
 
29
        int *tau, int *alpha, int *actrs, int *pnorm, int *rnors, int *prers,
 
30
        int *partl, int *sstol, int *taufc, int *apsma, int *betao, int *betac,
 
31
        int *betas, int *betan, int *s, int *ss, int *ssf, int *qraux, int *u,
 
32
        int *fs, int *fjacb, int *we1, int *diff, int *delts, int *deltn, 
 
33
        int *t, int *tt, int *omega, int *fjacd, int *wrk1, int *wrk2, 
 
34
        int *wrk3, int *wrk4, int *wrk5, int *wrk6, int *wrk7, int *lwkmn);
 
35
void F_FUNC(dluno,DLUNO)(int *lun, char *fn, int fnlen);
 
36
void F_FUNC(dlunc,DLUNC)(int *lun);
 
37
 
 
38
 
 
39
 
 
40
/* callback to pass to DODRC; calls the Python function in the global structure |odr_global| */
 
41
void fcn_callback(int *n, int *m, int *np, int *nq, int *ldn, int *ldm,
 
42
                  int *ldnp, double *beta, double *xplusd, int *ifixb,
 
43
                  int *ifixx, int *ldfix, int *ideval, double *f,
 
44
                  double *fjacb, double *fjacd, int *istop)
 
45
{
 
46
  PyObject *arg01, *arglist;
 
47
  PyObject *result;
 
48
  PyArrayObject *result_array = NULL;
 
49
  PyArrayObject *pyXplusD;
 
50
 
 
51
  arg01 = PyTuple_New(2);
 
52
 
 
53
  if (*m != 1)
 
54
    {
 
55
      int dim2[2];
 
56
      dim2[0] = *m;
 
57
      dim2[1] = *n;
 
58
      pyXplusD = (PyArrayObject *) PyArray_FromDims(2, dim2, PyArray_DOUBLE);
 
59
      memcpy(pyXplusD->data, (void *)xplusd, (*m) * (*n) * sizeof(double));
 
60
    }
 
61
  else
 
62
    {
 
63
      int dim1[1];
 
64
      dim1[0] = *n;
 
65
      pyXplusD = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
66
      memcpy(pyXplusD->data, (void *)xplusd, (*n) * sizeof(double));
 
67
    }
 
68
 
 
69
  PyTuple_SetItem(arg01, 0, odr_global.pyBeta);
 
70
  Py_INCREF(odr_global.pyBeta);
 
71
  PyTuple_SetItem(arg01, 1, (PyObject *) pyXplusD);
 
72
  Py_INCREF((PyObject *) pyXplusD);
 
73
 
 
74
  if (odr_global.extra_args != NULL)
 
75
    {
 
76
      arglist = PySequence_Concat(arg01, odr_global.extra_args);
 
77
    }
 
78
  else
 
79
    {
 
80
      arglist = PySequence_Tuple(arg01);        /* make a copy */
 
81
    }
 
82
 
 
83
  Py_DECREF(arg01);
 
84
  *istop = 0;
 
85
 
 
86
  memcpy(((PyArrayObject *) (odr_global.pyBeta))->data, (void *)beta,
 
87
         (*np) * sizeof(double));
 
88
 
 
89
  if ((*ideval % 10) >= 1)
 
90
    {
 
91
      /* compute f with odr_global.fcn */
 
92
 
 
93
      if (odr_global.fcn == NULL)
 
94
        {
 
95
          /* we don't have a function to call */
 
96
          PYERR2(odr_error, "Function has not been initialized");
 
97
        }
 
98
 
 
99
      if ((result = PyEval_CallObject(odr_global.fcn, arglist)) == NULL)
 
100
        {
 
101
          PyObject *tmpobj, *str1;
 
102
 
 
103
          if (PyErr_ExceptionMatches(odr_stop))
 
104
            {
 
105
              /* stop, don't fail */
 
106
              *istop = 1;
 
107
 
 
108
              Py_DECREF(arglist);
 
109
              return;
 
110
            }
 
111
 
 
112
          PyErr_Print();
 
113
          tmpobj = PyObject_GetAttrString(odr_global.fcn, "func_name");
 
114
          if (tmpobj == NULL)
 
115
            goto fail;
 
116
 
 
117
          str1 =
 
118
            PyString_FromString
 
119
            ("Error occured while calling the Python function named ");
 
120
          if (str1 == NULL)
 
121
            {
 
122
              Py_DECREF(tmpobj);
 
123
              goto fail;
 
124
            }
 
125
          PyString_ConcatAndDel(&str1, tmpobj);
 
126
          PyErr_SetString(odr_error, PyString_AsString(str1));
 
127
          Py_DECREF(str1);
 
128
          goto fail;
 
129
        }
 
130
 
 
131
      if ((result_array =
 
132
           (PyArrayObject *) PyArray_ContiguousFromObject(result,
 
133
                                                          PyArray_DOUBLE, 0,
 
134
                                                          2)) == NULL)
 
135
        {
 
136
          PYERR2(odr_error,
 
137
                 "Result from function call is not a proper array of floats.");
 
138
        }
 
139
 
 
140
      memcpy((void *)f, result_array->data, (*n) * (*nq) * sizeof(double));
 
141
      Py_DECREF(result_array);
 
142
    }
 
143
 
 
144
  if (((*ideval) / 10) % 10 >= 1)
 
145
    {
 
146
      /* compute fjacb with odr_global.fjacb */
 
147
 
 
148
      if (odr_global.fjacb == NULL)
 
149
        {
 
150
          /* we don't have a function to call */
 
151
          PYERR2(odr_error, "Function has not been initialized");
 
152
        }
 
153
 
 
154
      if ((result = PyEval_CallObject(odr_global.fjacb, arglist)) == NULL)
 
155
        {
 
156
          PyObject *tmpobj, *str1;
 
157
 
 
158
          if (PyErr_ExceptionMatches(odr_stop))
 
159
            {
 
160
              /* stop, don't fail */
 
161
              *istop = 1;
 
162
 
 
163
              Py_DECREF(arglist);
 
164
              return;
 
165
            }
 
166
 
 
167
          PyErr_Print();
 
168
          tmpobj = PyObject_GetAttrString(odr_global.fjacb, "func_name");
 
169
          if (tmpobj == NULL)
 
170
            goto fail;
 
171
 
 
172
          str1 =
 
173
            PyString_FromString
 
174
            ("Error occured while calling the Python function named ");
 
175
          if (str1 == NULL)
 
176
            {
 
177
              Py_DECREF(tmpobj);
 
178
              goto fail;
 
179
            }
 
180
          PyString_ConcatAndDel(&str1, tmpobj);
 
181
          PyErr_SetString(odr_error, PyString_AsString(str1));
 
182
          Py_DECREF(str1);
 
183
          goto fail;
 
184
        }
 
185
 
 
186
      if ((result_array =
 
187
           (PyArrayObject *) PyArray_ContiguousFromObject(result,
 
188
                                                          PyArray_DOUBLE, 0,
 
189
                                                          2)) == NULL)
 
190
        {
 
191
          PYERR2(odr_error,
 
192
                 "Result from function call is not a proper array of floats.");
 
193
        }
 
194
 
 
195
      if (*nq != 1 && *np != 1)
 
196
        {
 
197
          /* result_array should be rank-3 */
 
198
 
 
199
          if (result_array->nd != 3)
 
200
            {
 
201
              Py_DECREF(result_array);
 
202
              PYERR2(odr_error, "Beta Jacobian is not rank-3");
 
203
            }
 
204
        }
 
205
      else if (*nq == 1)
 
206
        {
 
207
          /* result_array should be rank-2 */
 
208
 
 
209
          if (result_array->nd != 2)
 
210
            {
 
211
              Py_DECREF(result_array);
 
212
              PYERR2(odr_error, "Beta Jacobian is not rank-2");
 
213
            }
 
214
        }
 
215
 
 
216
      memcpy((void *)fjacb, result_array->data,
 
217
             (*n) * (*nq) * (*np) * sizeof(double));
 
218
      Py_DECREF(result_array);
 
219
 
 
220
    }
 
221
 
 
222
  if (((*ideval) / 100) % 10 >= 1)
 
223
    {
 
224
      /* compute fjacd with odr_global.fjacd */
 
225
 
 
226
      if (odr_global.fjacd == NULL)
 
227
        {
 
228
          /* we don't have a function to call */
 
229
          PYERR2(odr_error, "fjcad has not been initialized");
 
230
        }
 
231
 
 
232
      if ((result = PyEval_CallObject(odr_global.fjacd, arglist)) == NULL)
 
233
        {
 
234
          PyObject *tmpobj, *str1;
 
235
 
 
236
          if (PyErr_ExceptionMatches(odr_stop))
 
237
            {
 
238
              /* stop, don't fail */
 
239
              *istop = 1;
 
240
 
 
241
              Py_DECREF(arglist);
 
242
              return;
 
243
            }
 
244
 
 
245
          PyErr_Print();
 
246
          tmpobj = PyObject_GetAttrString(odr_global.fjacd, "func_name");
 
247
          if (tmpobj == NULL)
 
248
            goto fail;
 
249
 
 
250
          str1 =
 
251
            PyString_FromString
 
252
            ("Error occured while calling the Python function named ");
 
253
          if (str1 == NULL)
 
254
            {
 
255
              Py_DECREF(tmpobj);
 
256
              goto fail;
 
257
            }
 
258
          PyString_ConcatAndDel(&str1, tmpobj);
 
259
          PyErr_SetString(odr_error, PyString_AsString(str1));
 
260
          Py_DECREF(str1);
 
261
          goto fail;
 
262
        }
 
263
 
 
264
      if ((result_array =
 
265
           (PyArrayObject *) PyArray_ContiguousFromObject(result,
 
266
                                                          PyArray_DOUBLE, 0,
 
267
                                                          2)) == NULL)
 
268
        {
 
269
          PYERR2(odr_error,
 
270
                 "Result from function call is not a proper array of floats.");
 
271
        }
 
272
 
 
273
      if (*nq != 1 && *m != 1)
 
274
        {
 
275
          /* result_array should be rank-3 */
 
276
 
 
277
          if (result_array->nd != 3)
 
278
            {
 
279
              Py_DECREF(result_array);
 
280
              PYERR2(odr_error, "xplusd Jacobian is not rank-3");
 
281
            }
 
282
        }
 
283
      else if (*nq == 1 && *m != 1)
 
284
        {
 
285
          /* result_array should be rank-2 */
 
286
 
 
287
          if (result_array->nd != 2)
 
288
            {
 
289
              Py_DECREF(result_array);
 
290
              PYERR2(odr_error, "xplusd Jacobian is not rank-2");
 
291
            }
 
292
        }
 
293
      else if (*nq == 1 && *m == 1)
 
294
        {
 
295
          /* result_array should be rank-1 */
 
296
 
 
297
          if (result_array->nd != 1)
 
298
            {
 
299
              Py_DECREF(result_array);
 
300
              PYERR2(odr_error, "xplusd Jacobian is not rank-1");
 
301
            }
 
302
        }
 
303
 
 
304
      memcpy((void *)fjacd, result_array->data,
 
305
             (*n) * (*nq) * (*m) * sizeof(double));
 
306
      Py_DECREF(result_array);
 
307
    }
 
308
 
 
309
  Py_DECREF(result);
 
310
  Py_DECREF(arglist);
 
311
  Py_DECREF(pyXplusD);
 
312
 
 
313
  return;
 
314
 
 
315
fail:
 
316
  Py_XDECREF(result);
 
317
  Py_XDECREF(arglist);
 
318
  Py_XDECREF(pyXplusD);
 
319
  *istop = -1;
 
320
  return;
 
321
}
 
322
 
 
323
 
 
324
/* generates Python output from the raw output from DODRC */
 
325
PyObject *gen_output(int n, int m, int np, int nq, int ldwe, int ld2we,
 
326
                     PyArrayObject * beta, PyArrayObject * work,
 
327
                     PyArrayObject * iwork, int isodr, int info,
 
328
                     int full_output)
 
329
{
 
330
  PyArrayObject *sd_beta, *cov_beta;
 
331
 
 
332
  int delta, eps, xplus, fn, sd, vcv, rvar, wss, wssde, wssep, rcond;
 
333
  int eta, olmav, tau, alpha, actrs, pnorm, rnors, prers, partl, sstol;
 
334
  int taufc, apsma, betao, betac, betas, betan, s, ss, ssf, qraux, u;
 
335
  int fs, fjacb, we1, diff, delts, deltn, t, tt, omega, fjacd;
 
336
  int wrk1, wrk2, wrk3, wrk4, wrk5, wrk6, wrk7, lwkmn;
 
337
 
 
338
  PyObject *retobj;
 
339
 
 
340
  int dim1[1], dim2[2];
 
341
 
 
342
  if (info == 50005) {
 
343
      /* fatal error in fcn call; return NULL to propogate the exception */
 
344
 
 
345
      return NULL;
 
346
  }
 
347
 
 
348
  lwkmn = work->dimensions[0];
 
349
 
 
350
  F_FUNC(dwinf,DWINF)(&n, &m, &np, &nq, &ldwe, &ld2we, &isodr,
 
351
        &delta, &eps, &xplus, &fn, &sd, &vcv, &rvar, &wss, &wssde,
 
352
        &wssep, &rcond, &eta, &olmav, &tau, &alpha, &actrs, &pnorm,
 
353
        &rnors, &prers, &partl, &sstol, &taufc, &apsma, &betao, &betac,
 
354
        &betas, &betan, &s, &ss, &ssf, &qraux, &u, &fs, &fjacb, &we1,
 
355
        &diff, &delts, &deltn, &t, &tt, &omega, &fjacd, &wrk1, &wrk2,
 
356
        &wrk3, &wrk4, &wrk5, &wrk6, &wrk7, &lwkmn);
 
357
 
 
358
  /* convert FORTRAN indices to C indices */
 
359
  delta--;
 
360
  eps--;
 
361
  xplus--;
 
362
  fn--;
 
363
  sd--;
 
364
  vcv--;
 
365
  rvar--;
 
366
  wss--;
 
367
  wssde--;
 
368
  wssep--;
 
369
  rcond--;
 
370
  eta--;
 
371
  olmav--;
 
372
  tau--;
 
373
  alpha--;
 
374
  actrs--;
 
375
  pnorm--;
 
376
  rnors--;
 
377
  prers--;
 
378
  partl--;
 
379
  sstol--;
 
380
  taufc--;
 
381
  apsma--;
 
382
  betao--;
 
383
  betac--;
 
384
  betas--;
 
385
  betan--;
 
386
  s--;
 
387
  ss--;
 
388
  ssf--;
 
389
  qraux--;
 
390
  u--;
 
391
  fs--;
 
392
  fjacb--;
 
393
  we1--;
 
394
  diff--;
 
395
  delts--;
 
396
  deltn--;
 
397
  t--;
 
398
  tt--;
 
399
  omega--;
 
400
  fjacd--;
 
401
  wrk1--;
 
402
  wrk2--;
 
403
  wrk3--;
 
404
  wrk4--;
 
405
  wrk5--;
 
406
  wrk6--;
 
407
  wrk7--;
 
408
 
 
409
  dim1[0] = beta->dimensions[0];
 
410
  sd_beta = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
411
  dim2[0] = beta->dimensions[0];
 
412
  dim2[1] = beta->dimensions[0];
 
413
  cov_beta = (PyArrayObject *) PyArray_FromDims(2, dim2, PyArray_DOUBLE);
 
414
 
 
415
  memcpy(sd_beta->data, (void *)((double *)(work->data) + sd),
 
416
         np * sizeof(double));
 
417
  memcpy(cov_beta->data, (void *)((double *)(work->data) + vcv),
 
418
         np * np * sizeof(double));
 
419
 
 
420
  if (!full_output)
 
421
    {
 
422
      retobj = Py_BuildValue("OOO", PyArray_Return(beta),
 
423
                             PyArray_Return(sd_beta),
 
424
                             PyArray_Return(cov_beta));
 
425
      Py_DECREF((PyObject *) sd_beta);
 
426
      Py_DECREF((PyObject *) cov_beta);
 
427
 
 
428
      return retobj;
 
429
    }
 
430
  else
 
431
    {
 
432
      PyArrayObject *deltaA, *epsA, *xplusA, *fnA;
 
433
      double res_var, sum_square, sum_square_delta, sum_square_eps;
 
434
      double inv_condnum, rel_error;
 
435
      PyObject *work_ind;
 
436
 
 
437
      work_ind =
 
438
        Py_BuildValue
 
439
        ("{s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l,s:l}",
 
440
         "delta", delta, "eps", eps, "xplus", xplus, "fn", fn, "sd", sd, "sd",
 
441
         vcv, "rvar", rvar, "wss", wss, "wssde", wssde, "wssep", wssep,
 
442
         "rcond", rcond, "eta", eta, "olmav", olmav, "tau", tau, "alpha",
 
443
         alpha, "actrs", actrs, "pnorm", pnorm, "rnors", rnors, "prers",
 
444
         prers, "partl", partl, "sstol", sstol, "taufc", taufc, "apsma",
 
445
         apsma, "betao", betao, "betac", betac, "betas", betas, "betan",
 
446
         betan, "s", s, "ss", ss, "ssf", ssf, "qraux", qraux, "u", u, "fs",
 
447
         fs, "fjacb", fjacb, "we1", we1, "diff", diff, "delts", delts,
 
448
         "deltn", deltn, "t", t, "tt", tt, "omega", omega, "fjacd", fjacd,
 
449
         "wrk1", wrk1, "wrk2", wrk2, "wrk3", wrk3, "wrk4", wrk4, "wrk5", wrk5,
 
450
         "wrk6", wrk6, "wrk7", wrk7);
 
451
 
 
452
      if (m == 1)
 
453
        {
 
454
          dim1[0] = n;
 
455
          deltaA =
 
456
            (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
457
          xplusA =
 
458
            (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
459
        }
 
460
      else
 
461
        {
 
462
          dim2[0] = m;
 
463
          dim2[1] = n;
 
464
          deltaA =
 
465
            (PyArrayObject *) PyArray_FromDims(2, dim2, PyArray_DOUBLE);
 
466
          xplusA =
 
467
            (PyArrayObject *) PyArray_FromDims(2, dim2, PyArray_DOUBLE);
 
468
        }
 
469
 
 
470
      if (nq == 1)
 
471
        {
 
472
          dim1[0] = n;
 
473
          epsA = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
474
          fnA = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
475
        }
 
476
      else
 
477
        {
 
478
          dim2[0] = nq;
 
479
          dim2[1] = n;
 
480
          epsA = (PyArrayObject *) PyArray_FromDims(2, dim2, PyArray_DOUBLE);
 
481
          fnA = (PyArrayObject *) PyArray_FromDims(2, dim2, PyArray_DOUBLE);
 
482
        }
 
483
 
 
484
      memcpy(deltaA->data, (void *)((double *)(work->data) + delta),
 
485
             m * n * sizeof(double));
 
486
      memcpy(epsA->data, (void *)((double *)(work->data) + eps),
 
487
             nq * n * sizeof(double));
 
488
      memcpy(xplusA->data, (void *)((double *)(work->data) + xplus),
 
489
             m * n * sizeof(double));
 
490
      memcpy(fnA->data, (void *)((double *)(work->data) + fn),
 
491
             nq * n * sizeof(double));
 
492
 
 
493
      res_var = *((double *)(work->data) + rvar);
 
494
      sum_square = *((double *)(work->data) + wss);
 
495
      sum_square_delta = *((double *)(work->data) + wssde);
 
496
      sum_square_eps = *((double *)(work->data) + wssep);
 
497
      inv_condnum = *((double *)(work->data) + rcond);
 
498
      rel_error = *((double *)(work->data) + eta);
 
499
 
 
500
      retobj =
 
501
        Py_BuildValue
 
502
        ("OOO{s:O,s:O,s:O,s:O,s:d,s:d,s:d,s:d,s:d,s:d,s:O,s:O,s:O,s:l}",
 
503
         PyArray_Return(beta), PyArray_Return(sd_beta),
 
504
         PyArray_Return(cov_beta), "delta", PyArray_Return(deltaA), "eps",
 
505
         PyArray_Return(epsA), "xplus", PyArray_Return(xplusA), "y",
 
506
         PyArray_Return(fnA), "res_var", res_var, "sum_square", sum_square,
 
507
         "sum_square_delta", sum_square_delta, "sum_square_eps",
 
508
         sum_square_eps, "inv_condnum", inv_condnum, "rel_error", rel_error,
 
509
         "work", PyArray_Return(work), "work_ind", work_ind, "iwork",
 
510
         PyArray_Return(iwork), "info", info);
 
511
      Py_DECREF((PyObject *) sd_beta);
 
512
      Py_DECREF((PyObject *) cov_beta);
 
513
      Py_DECREF((PyObject *) deltaA);
 
514
      Py_DECREF((PyObject *) epsA);
 
515
      Py_DECREF((PyObject *) xplusA);
 
516
      Py_DECREF((PyObject *) fnA);
 
517
      Py_DECREF((PyObject *) work_ind);
 
518
 
 
519
      return retobj;
 
520
    }
 
521
}
 
522
 
 
523
PyObject *odr(PyObject * self, PyObject * args, PyObject * kwds)
 
524
{
 
525
  PyObject *fcn, *initbeta, *py, *px, *pwe = NULL, *pwd = NULL, *fjacb = NULL;
 
526
  PyObject *fjacd = NULL, *pifixb = NULL, *pifixx = NULL;
 
527
  PyObject *pstpb = NULL, *pstpd = NULL, *psclb = NULL, *pscld = NULL;
 
528
  PyObject *pwork = NULL, *piwork = NULL, *extra_args = NULL;
 
529
  int job = 0, ndigit = 0, maxit = -1, iprint = 0;
 
530
  int full_output = 0;
 
531
  double taufac = 0.0, sstol = -1.0, partol = -1.0;
 
532
  char *errfile = NULL, *rptfile = NULL;
 
533
  int lerrfile = 0, lrptfile = 0;
 
534
  PyArrayObject *beta = NULL, *y = NULL, *x = NULL, *we = NULL, *wd = NULL;
 
535
  PyArrayObject *ifixb = NULL, *ifixx = NULL;
 
536
  PyArrayObject *stpb = NULL, *stpd = NULL, *sclb = NULL, *scld = NULL;
 
537
  PyArrayObject *work = NULL, *iwork = NULL;
 
538
  int n, m, np, nq, ldy, ldx, ldwe, ld2we, ldwd, ld2wd, ldifx;
 
539
  int lunerr = -1, lunrpt = -1, ldstpd, ldscld, lwork, liwork, info = 0;
 
540
  static char *kw_list[] = { "fcn", "initbeta", "y", "x", "we", "wd", "fjacb",
 
541
    "fjacd", "extra_args", "ifixb", "ifixx", "job", "iprint", "errfile",
 
542
    "rptfile", "ndigit", "taufac", "sstol", "partol",
 
543
    "maxit", "stpb", "stpd", "sclb", "scld", "work",
 
544
    "iwork", "full_output", NULL
 
545
  };
 
546
  int isodr = 1;
 
547
  PyObject *result;
 
548
  int dim1[1], dim2[2], dim3[3];
 
549
  int implicit;                 /* flag for implicit model */
 
550
 
 
551
 
 
552
  if (kwds == NULL)
 
553
    {
 
554
      if (!PyArg_ParseTuple(args, "OOOO|OOOOOOOllz#z#ldddlOOOOOOi:odr",
 
555
                            &fcn, &initbeta, &py, &px, &pwe, &pwd,
 
556
                            &fjacb, &fjacd, &extra_args, &pifixb, &pifixx,
 
557
                            &job, &iprint, &errfile, &lerrfile, &rptfile,
 
558
                            &lrptfile, &ndigit, &taufac, &sstol, &partol,
 
559
                            &maxit, &pstpb, &pstpd, &psclb, &pscld, &pwork,
 
560
                            &piwork, &full_output))
 
561
        {
 
562
          return NULL;
 
563
        }
 
564
    }
 
565
  else
 
566
    {
 
567
      if (!PyArg_ParseTupleAndKeywords(args, kwds,
 
568
                                       "OOOO|OOOOOOOllz#z#ldddlOOOOOOi:odr",
 
569
                                       kw_list, &fcn, &initbeta, &py, &px,
 
570
                                       &pwe, &pwd, &fjacb, &fjacd,
 
571
                                       &extra_args, &pifixb, &pifixx, &job,
 
572
                                       &iprint, &errfile, &lerrfile, &rptfile,
 
573
                                       &lrptfile, &ndigit, &taufac, &sstol,
 
574
                                       &partol, &maxit, &pstpb, &pstpd,
 
575
                                       &psclb, &pscld, &pwork, &piwork,
 
576
                                       &full_output))
 
577
        {
 
578
          return NULL;
 
579
        }
 
580
    }
 
581
 
 
582
  /* Check the validity of all arguments */
 
583
 
 
584
  if (!PyCallable_Check(fcn))
 
585
    {
 
586
      PYERR(PyExc_TypeError, "fcn must be callable");
 
587
    }
 
588
  if (!PySequence_Check(initbeta))
 
589
    {
 
590
      PYERR(PyExc_TypeError, "initbeta must be a sequence");
 
591
    }
 
592
  if (!PySequence_Check(py) && !PyInt_Check(py))
 
593
    {
 
594
      PYERR(PyExc_TypeError,
 
595
            "y must be a sequence or integer (if model is implicit)");
 
596
    }
 
597
  if (!PySequence_Check(px))
 
598
    {
 
599
      PYERR(PyExc_TypeError, "x must be a sequence");
 
600
    }
 
601
  if (pwe != NULL && !PySequence_Check(pwe) && !PyNumber_Check(pwe))
 
602
    {
 
603
      PYERR(PyExc_TypeError, "we must be a sequence or a number");
 
604
    }
 
605
  if (pwd != NULL && !PySequence_Check(pwd) && !PyNumber_Check(pwd))
 
606
    {
 
607
      PYERR(PyExc_TypeError, "wd must be a sequence or a number");
 
608
    }
 
609
  if (fjacb != NULL && !PyCallable_Check(fjacb))
 
610
    {
 
611
      PYERR(PyExc_TypeError, "fjacb must be callable");
 
612
    }
 
613
  if (fjacd != NULL && !PyCallable_Check(fjacd))
 
614
    {
 
615
      PYERR(PyExc_TypeError, "fjacd must be callable");
 
616
    }
 
617
  if (extra_args != NULL && !PySequence_Check(extra_args))
 
618
    {
 
619
      PYERR(PyExc_TypeError, "extra_args must be a sequence");
 
620
    }
 
621
  if (pifixx != NULL && !PySequence_Check(pifixx))
 
622
    {
 
623
      PYERR(PyExc_TypeError, "ifixx must be a sequence");
 
624
    }
 
625
  if (pifixb != NULL && !PySequence_Check(pifixb))
 
626
    {
 
627
      PYERR(PyExc_TypeError, "ifixb must be a sequence");
 
628
    }
 
629
  if (pstpb != NULL && !PySequence_Check(pstpb))
 
630
    {
 
631
      PYERR(PyExc_TypeError, "stpb must be a sequence");
 
632
    }
 
633
  if (pstpd != NULL && !PySequence_Check(pstpd))
 
634
    {
 
635
      PYERR(PyExc_TypeError, "stpd must be a sequence");
 
636
    }
 
637
  if (psclb != NULL && !PySequence_Check(psclb))
 
638
    {
 
639
      PYERR(PyExc_TypeError, "sclb must be a sequence");
 
640
    }
 
641
  if (pscld != NULL && !PySequence_Check(pscld))
 
642
    {
 
643
      PYERR(PyExc_TypeError, "scld must be a sequence");
 
644
    }
 
645
  if (pwork != NULL && !PyArray_Check(pwork))
 
646
    {
 
647
      PYERR(PyExc_TypeError, "work must be an array");
 
648
    }
 
649
  if (piwork != NULL && !PyArray_Check(piwork))
 
650
    {
 
651
      PYERR(PyExc_TypeError, "iwork must be an array");
 
652
    }
 
653
 
 
654
  /* start processing the arguments and check for errors on the way */
 
655
 
 
656
  /* check for implicit model */
 
657
 
 
658
  implicit = (job % 10 == 1);
 
659
 
 
660
  if (!implicit)
 
661
    {
 
662
      if ((y =
 
663
           (PyArrayObject *) PyArray_CopyFromObject(py, PyArray_DOUBLE, 1,
 
664
                                                    2)) == NULL)
 
665
        {
 
666
          PYERR(PyExc_ValueError,
 
667
                "y could not be made into a suitable array");
 
668
        }
 
669
      n = y->dimensions[y->nd - 1];     /* pick the last dimension */
 
670
      if ((x =
 
671
           (PyArrayObject *) PyArray_CopyFromObject(px, PyArray_DOUBLE, 1,
 
672
                                                    2)) == NULL)
 
673
        {
 
674
          PYERR(PyExc_ValueError,
 
675
                "x could not be made into a suitable array");
 
676
        }
 
677
      if (n != x->dimensions[x->nd - 1])
 
678
        {
 
679
          PYERR(PyExc_ValueError,
 
680
                "x and y don't have matching numbers of observations");
 
681
        }
 
682
      if (y->nd == 1)
 
683
        {
 
684
          nq = 1;
 
685
        }
 
686
      else
 
687
        {
 
688
          nq = y->dimensions[0];
 
689
        }
 
690
 
 
691
      ldx = ldy = n;
 
692
    }
 
693
  else
 
694
    {                           /* we *do* have an implicit model */
 
695
      ldy = 1;
 
696
      nq = (int)PyInt_AsLong(py);
 
697
      dim1[0] = 1;
 
698
 
 
699
      /* initialize y to a dummy array; never referenced */
 
700
      y = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
701
 
 
702
      if ((x =
 
703
           (PyArrayObject *) PyArray_CopyFromObject(px, PyArray_DOUBLE, 1,
 
704
                                                    2)) == NULL)
 
705
        {
 
706
          PYERR(PyExc_ValueError,
 
707
                "x could not be made into a suitable array");
 
708
        }
 
709
 
 
710
      n = x->dimensions[x->nd - 1];
 
711
      ldx = n;
 
712
    }
 
713
 
 
714
  if (x->nd == 1)
 
715
    {
 
716
      m = 1;
 
717
    }
 
718
  else
 
719
    {
 
720
      m = x->dimensions[0];
 
721
    }                           /* x, y */
 
722
 
 
723
  if ((beta =
 
724
       (PyArrayObject *) PyArray_CopyFromObject(initbeta, PyArray_DOUBLE, 1,
 
725
                                                1)) == NULL)
 
726
    {
 
727
      PYERR(PyExc_ValueError,
 
728
            "initbeta could not be made into a suitable array");
 
729
    }
 
730
  np = beta->dimensions[0];
 
731
 
 
732
  if (pwe == NULL)
 
733
    {
 
734
      ldwe = ld2we = 1;
 
735
      dim1[0] = n;
 
736
      we = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
737
      ((double *)(we->data))[0] = -1.0;
 
738
    }
 
739
  else if (PyNumber_Check(pwe) && !PyArray_Check(pwe))
 
740
    {
 
741
      /* we is a single weight, set the first value of we to -pwe */
 
742
      PyObject *tmp;
 
743
      double val;
 
744
 
 
745
      tmp = PyNumber_Float(pwe);
 
746
      if (tmp == NULL)
 
747
        PYERR(PyExc_ValueError, "could not convert we to a suitable array");
 
748
      val = PyFloat_AsDouble(tmp);
 
749
      Py_DECREF(tmp);
 
750
 
 
751
      dim3[0] = nq;
 
752
      dim3[1] = 1;
 
753
      dim3[2] = 1;
 
754
      we = (PyArrayObject *) PyArray_FromDims(3, dim3, PyArray_DOUBLE);
 
755
      if (implicit)
 
756
        {
 
757
          ((double *)(we->data))[0] = val;
 
758
        }
 
759
      else
 
760
        {
 
761
          ((double *)(we->data))[0] = -val;
 
762
        }
 
763
      ldwe = ld2we = 1;
 
764
    }
 
765
  else if (PySequence_Check(pwe))
 
766
    {
 
767
      /* we needs to be turned into an array */
 
768
 
 
769
      if ((we =
 
770
           (PyArrayObject *) PyArray_CopyFromObject(pwe, PyArray_DOUBLE, 1,
 
771
                                                    3)) == NULL)
 
772
        {
 
773
          PYERR(PyExc_ValueError, "could not convert we to a suitable array");
 
774
        }
 
775
 
 
776
      if (we->nd == 1 && nq == 1)
 
777
        {
 
778
 
 
779
          ldwe = n;
 
780
          ld2we = 1;
 
781
        }
 
782
      else if (we->nd == 1 && we->dimensions[0] == nq)
 
783
        {
 
784
          /* we is a rank-1 array with diagonal weightings to be broadcast 
 
785
           * to all observations */
 
786
          ldwe = 1;
 
787
          ld2we = 1;
 
788
        }
 
789
      else if (we->nd == 3 && we->dimensions[0] == nq
 
790
               && we->dimensions[1] == nq && we->dimensions[2] == 1)
 
791
        {
 
792
          /* we is a rank-3 array with the covariant weightings 
 
793
             to be broadcast to all observations */
 
794
          ldwe = 1;
 
795
          ld2we = nq;
 
796
        }
 
797
      else if (we->nd == 2 && we->dimensions[0] == nq
 
798
               && we->dimensions[1] == nq)
 
799
        {
 
800
          /* we is a rank-2 array with the full covariant weightings 
 
801
             to be broadcast to all observations */
 
802
          ldwe = 1;
 
803
          ld2we = nq;
 
804
        }
 
805
 
 
806
      else if (we->nd == 2 && we->dimensions[0] == nq
 
807
               && we->dimensions[1] == n)
 
808
        {
 
809
          /* we is a rank-2 array with the diagonal elements of the 
 
810
             covariant weightings for each observation */
 
811
          ldwe = n;
 
812
          ld2we = 1;
 
813
        }
 
814
      else if (we->nd == 3 && we->dimensions[0] == nq
 
815
               && we->dimensions[1] == nq && we->dimensions[2] == n)
 
816
        {
 
817
          /* we is the full specification of the covariant weights
 
818
             for each observation */
 
819
          ldwe = n;
 
820
          ld2we = nq;
 
821
        }
 
822
      else
 
823
        {
 
824
          PYERR(PyExc_ValueError, "could not convert we to a suitable array");
 
825
        }
 
826
    }                           /* we */
 
827
 
 
828
  if (pwd == NULL)
 
829
    {
 
830
      ldwd = ld2wd = 1;
 
831
 
 
832
      dim1[0] = m;
 
833
      wd = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
834
      ((double *)(wd->data))[0] = -1.0;
 
835
    }
 
836
  else if (PyNumber_Check(pwd) && !PyArray_Check(pwd))
 
837
    {
 
838
      /* wd is a single weight, set the first value of wd to -pwd */
 
839
      PyObject *tmp;
 
840
      double val;
 
841
 
 
842
      tmp = PyNumber_Float(pwd);
 
843
      if (tmp == NULL)
 
844
        PYERR(PyExc_ValueError, "could not convert wd to a suitable array");
 
845
      val = PyFloat_AsDouble(tmp);
 
846
      Py_DECREF(tmp);
 
847
 
 
848
      dim3[0] = 1;
 
849
      dim3[1] = 1;
 
850
      dim3[2] = m;
 
851
      wd = (PyArrayObject *) PyArray_FromDims(3, dim3, PyArray_DOUBLE);
 
852
      ((double *)(wd->data))[0] = -val;
 
853
      ldwd = ld2wd = 1;
 
854
    }
 
855
  else if (PySequence_Check(pwd))
 
856
    {
 
857
      /* wd needs to be turned into an array */
 
858
 
 
859
      if ((wd =
 
860
           (PyArrayObject *) PyArray_CopyFromObject(pwd, PyArray_DOUBLE, 1,
 
861
                                                    3)) == NULL)
 
862
        {
 
863
          PYERR(PyExc_ValueError, "could not convert wd to a suitable array");
 
864
        }
 
865
 
 
866
      if (wd->nd == 1 && m == 1)
 
867
        {
 
868
          ldwd = n;
 
869
          ld2wd = 1;
 
870
        }
 
871
      else if (wd->nd == 1 && wd->dimensions[0] == m)
 
872
        {
 
873
          /* wd is a rank-1 array with diagonal weightings to be broadcast 
 
874
           * to all observations */
 
875
          ldwd = 1;
 
876
          ld2wd = 1;
 
877
        }
 
878
 
 
879
      else if (wd->nd == 3 && wd->dimensions[0] == m
 
880
               && wd->dimensions[1] == m && wd->dimensions[2] == 1)
 
881
        {
 
882
          /* wd is a rank-3 array with the covariant wdightings 
 
883
             to be broadcast to all observations */
 
884
          ldwd = 1;
 
885
          ld2wd = m;
 
886
        }
 
887
      else if (wd->nd == 2 && wd->dimensions[0] == m
 
888
               && wd->dimensions[1] == m)
 
889
        {
 
890
          /* wd is a rank-2 array with the full covariant weightings 
 
891
             to be broadcast to all observations */
 
892
          ldwd = 1;
 
893
          ld2wd = m;
 
894
        }
 
895
 
 
896
      else if (wd->nd == 2 && wd->dimensions[0] == m
 
897
               && wd->dimensions[1] == n)
 
898
        {
 
899
          /* wd is a rank-2 array with the diagonal elements of the 
 
900
             covariant weightings for each observation */
 
901
          ldwd = n;
 
902
          ld2wd = 1;
 
903
        }
 
904
      else if (wd->nd == 3 && wd->dimensions[0] == m
 
905
               && wd->dimensions[1] == m && wd->dimensions[2] == n)
 
906
        {
 
907
          /* wd is the full specification of the covariant weights
 
908
             for each observation */
 
909
          ldwd = n;
 
910
          ld2wd = m;
 
911
        }
 
912
      else
 
913
        {
 
914
          PYERR(PyExc_ValueError, "could not convert wd to a suitable array");
 
915
        }
 
916
 
 
917
    }                           /* wd */
 
918
 
 
919
 
 
920
  if (pifixb == NULL)
 
921
    {
 
922
      dim1[0] = np;
 
923
      ifixb = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_LONG);
 
924
      *(int *)(ifixb->data) = -1;      /* set first element negative */
 
925
    }
 
926
  else
 
927
    {
 
928
      /* pifixb is a sequence as checked before */
 
929
 
 
930
      if ((ifixb =
 
931
           (PyArrayObject *) PyArray_CopyFromObject(pifixb, PyArray_LONG, 1,
 
932
                                                    1)) == NULL)
 
933
        {
 
934
          PYERR(PyExc_ValueError,
 
935
                "could not convert ifixb to a suitable array");
 
936
        }
 
937
 
 
938
      if (ifixb->dimensions[0] != np)
 
939
        {
 
940
          PYERR(PyExc_ValueError,
 
941
                "could not convert ifixb to a suitable array");
 
942
        }
 
943
    }                           /* ifixb */
 
944
 
 
945
  if (pifixx == NULL)
 
946
    {
 
947
      dim2[0] = m;
 
948
      dim2[1] = 1;
 
949
      ifixx = (PyArrayObject *) PyArray_FromDims(2, dim2, PyArray_LONG);
 
950
      *(int *)(ifixx->data) = -1;      /* set first element negative */
 
951
      ldifx = 1;
 
952
    }
 
953
  else
 
954
    {
 
955
      /* pifixx is a sequence as checked before */
 
956
 
 
957
      if ((ifixx =
 
958
           (PyArrayObject *) PyArray_CopyFromObject(pifixx, PyArray_LONG, 1,
 
959
                                                    2)) == NULL)
 
960
        {
 
961
          PYERR(PyExc_ValueError,
 
962
                "could not convert ifixx to a suitable array");
 
963
        }
 
964
 
 
965
      if (ifixx->nd == 1 && ifixx->dimensions[0] == m)
 
966
        {
 
967
          ldifx = 1;
 
968
        }
 
969
      else if (ifixx->nd == 1 && ifixx->dimensions[0] == n && m == 1)
 
970
        {
 
971
          ldifx = n;
 
972
        }
 
973
      else if (ifixx->nd == 2 && ifixx->dimensions[0] == m
 
974
               && ifixx->dimensions[1] == n)
 
975
        {
 
976
          ldifx = n;
 
977
        }
 
978
      else
 
979
        {
 
980
          PYERR(PyExc_ValueError,
 
981
                "could not convert ifixx to a suitable array");
 
982
        }
 
983
    }                           /* ifixx */
 
984
 
 
985
  if (errfile != NULL)
 
986
    {
 
987
      /* call FORTRAN's OPEN to open the file with a logical unit of 18 */
 
988
      lunerr = 18;
 
989
      F_FUNC(dluno,DLUNO)(&lunerr, errfile, lerrfile);
 
990
    }
 
991
 
 
992
  if (rptfile != NULL)
 
993
    {
 
994
      /* call FORTRAN's OPEN to open the file with a logical unit of 19 */
 
995
      lunrpt = 19;
 
996
      F_FUNC(dluno,DLUNO)(&lunrpt, rptfile, lrptfile);
 
997
    }
 
998
 
 
999
  if (pstpb == NULL)
 
1000
    {
 
1001
      dim1[0] = np;
 
1002
      stpb = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
1003
      *(double *)(stpb->data) = 0.0;
 
1004
    }
 
1005
  else                          /* pstpb is a sequence */
 
1006
    {
 
1007
      if ((stpb =
 
1008
           (PyArrayObject *) PyArray_CopyFromObject(pstpb, PyArray_DOUBLE, 1,
 
1009
                                                    1)) == NULL
 
1010
          || stpb->dimensions[0] != np)
 
1011
        {
 
1012
          PYERR(PyExc_ValueError,
 
1013
                "could not convert stpb to a suitable array");
 
1014
        }
 
1015
    }                           /* stpb */
 
1016
 
 
1017
  if (pstpd == NULL)
 
1018
    {
 
1019
      dim2[0] = 1;
 
1020
      dim2[1] = m;
 
1021
      stpd = (PyArrayObject *) PyArray_FromDims(2, dim2, PyArray_DOUBLE);
 
1022
      *(double *)(stpd->data) = 0.0;
 
1023
      ldstpd = 1;
 
1024
    }
 
1025
  else
 
1026
    {
 
1027
      if ((stpd =
 
1028
           (PyArrayObject *) PyArray_CopyFromObject(pstpd, PyArray_DOUBLE, 1,
 
1029
                                                    2)) == NULL)
 
1030
        {
 
1031
          PYERR(PyExc_ValueError,
 
1032
                "could not convert stpb to a suitable array");
 
1033
        }
 
1034
 
 
1035
      if (stpd->nd == 1 && stpd->dimensions[0] == m)
 
1036
        {
 
1037
          ldstpd = 1;
 
1038
        }
 
1039
      else if (stpd->nd == 1 && stpd->dimensions[0] == n && m == 1)
 
1040
        {
 
1041
          ldstpd = n;
 
1042
        }
 
1043
      else if (stpd->nd == 2 && stpd->dimensions[0] == n &&
 
1044
               stpd->dimensions[1] == m)
 
1045
        {
 
1046
          ldstpd = n;
 
1047
        }
 
1048
    }                           /* stpd */
 
1049
 
 
1050
  if (psclb == NULL)
 
1051
    {
 
1052
      dim1[0] = np;
 
1053
      sclb = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
1054
      *(double *)(sclb->data) = 0.0;
 
1055
    }
 
1056
  else                          /* psclb is a sequence */
 
1057
    {
 
1058
      if ((sclb =
 
1059
           (PyArrayObject *) PyArray_CopyFromObject(psclb, PyArray_DOUBLE, 1,
 
1060
                                                    1)) == NULL
 
1061
          || sclb->dimensions[0] != np)
 
1062
        {
 
1063
          PYERR(PyExc_ValueError,
 
1064
                "could not convert sclb to a suitable array");
 
1065
        }
 
1066
    }                           /* sclb */
 
1067
 
 
1068
  if (pscld == NULL)
 
1069
    {
 
1070
      dim2[0] = 1;
 
1071
      dim2[1] = n;
 
1072
      scld = (PyArrayObject *) PyArray_FromDims(2, dim2, PyArray_DOUBLE);
 
1073
      *(double *)(scld->data) = 0.0;
 
1074
      ldscld = 1;
 
1075
    }
 
1076
  else
 
1077
    {
 
1078
      if ((scld =
 
1079
           (PyArrayObject *) PyArray_CopyFromObject(pscld, PyArray_DOUBLE, 1,
 
1080
                                                    2)) == NULL)
 
1081
        {
 
1082
          PYERR(PyExc_ValueError,
 
1083
                "could not convert stpb to a suitable array");
 
1084
        }
 
1085
 
 
1086
      if (scld->nd == 1 && scld->dimensions[0] == m)
 
1087
        {
 
1088
          ldscld = 1;
 
1089
        }
 
1090
      else if (scld->nd == 1 && scld->dimensions[0] == n && m == 1)
 
1091
        {
 
1092
          ldscld = n;
 
1093
        }
 
1094
      else if (scld->nd == 2 && scld->dimensions[0] == n &&
 
1095
               scld->dimensions[1] == m)
 
1096
        {
 
1097
          ldscld = n;
 
1098
        }
 
1099
    }                           /* scld */
 
1100
 
 
1101
  if (job % 10 < 2)
 
1102
    {
 
1103
      /* ODR, not OLS */
 
1104
 
 
1105
      lwork =
 
1106
        18 + 11 * np + np * np + m + m * m + 4 * n * nq + 6 * n * m +
 
1107
        2 * n * nq * np + 2 * n * nq * m + nq * nq + 5 * nq + nq * (np + m) +
 
1108
        ldwe * ld2we * nq;
 
1109
 
 
1110
      isodr = 1;
 
1111
    }
 
1112
  else
 
1113
    {
 
1114
      /* OLS, not ODR */
 
1115
 
 
1116
      lwork =
 
1117
        18 + 11 * np + np * np + m + m * m + 4 * n * nq + 2 * n * m +
 
1118
        2 * n * nq * np + 5 * nq + nq * (np + m) + ldwe * ld2we * nq;
 
1119
 
 
1120
      isodr = 0;
 
1121
    }
 
1122
 
 
1123
  liwork = 20 + np + nq * (np + m);
 
1124
 
 
1125
  if ((job / 10000) % 10 >= 1)
 
1126
    {
 
1127
      /* fit is a restart, make sure work and iwork are input */
 
1128
 
 
1129
      if (pwork == NULL || piwork == NULL)
 
1130
        {
 
1131
          PYERR(PyExc_ValueError,
 
1132
                "need to input work and iwork arrays to restart");
 
1133
        }
 
1134
    }
 
1135
 
 
1136
  if ((job / 1000) % 10 >= 1)
 
1137
    {
 
1138
      /* delta should be supplied, make sure the user does */
 
1139
 
 
1140
      if (pwork == NULL)
 
1141
        {
 
1142
          PYERR(PyExc_ValueError,
 
1143
                "need to input work array for delta initialization");
 
1144
        }
 
1145
    }
 
1146
 
 
1147
  if (pwork != NULL)
 
1148
    {
 
1149
      if ((work =
 
1150
           (PyArrayObject *) PyArray_CopyFromObject(pwork, PyArray_DOUBLE, 1,
 
1151
                                                    1)) == NULL)
 
1152
        {
 
1153
          PYERR(PyExc_ValueError,
 
1154
                "could not convert work to a suitable array");
 
1155
        }
 
1156
      if (work->dimensions[0] < lwork)
 
1157
        {
 
1158
            printf("%d %d\n", work->dimensions[0], lwork);
 
1159
          PYERR(PyExc_ValueError, "work is too small");
 
1160
        }
 
1161
    }
 
1162
  else
 
1163
    {
 
1164
      /* initialize our own work array */
 
1165
      dim1[0] = lwork;
 
1166
      work = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_DOUBLE);
 
1167
    }                           /* work */
 
1168
 
 
1169
  if (piwork != NULL)
 
1170
    {
 
1171
      if ((iwork =
 
1172
           (PyArrayObject *) PyArray_CopyFromObject(piwork, PyArray_LONG, 1,
 
1173
                                                    1)) == NULL)
 
1174
        {
 
1175
          PYERR(PyExc_ValueError,
 
1176
                "could not convert iwork to a suitable array");
 
1177
        }
 
1178
 
 
1179
      if (iwork->dimensions[0] < liwork)
 
1180
        {
 
1181
          PYERR(PyExc_ValueError, "iwork is too small");
 
1182
        }
 
1183
    }
 
1184
  else
 
1185
    {
 
1186
      /* initialize our own iwork array */
 
1187
      dim1[0] = liwork;
 
1188
      iwork = (PyArrayObject *) PyArray_FromDims(1, dim1, PyArray_LONG);
 
1189
    }                           /* iwork */
 
1190
 
 
1191
  /* check if what JOB requests can be done with what the user has 
 
1192
     input into the function */
 
1193
 
 
1194
  if ((job / 10) % 10 >= 2)
 
1195
    {
 
1196
      /* derivatives are supposed to be supplied */
 
1197
 
 
1198
      if (fjacb == NULL || fjacd == NULL)
 
1199
        {
 
1200
          PYERR(PyExc_ValueError,
 
1201
                "need fjacb and fjacd to calculate derivatives");
 
1202
        }
 
1203
    }
 
1204
 
 
1205
  /* setup the global data for the callback */
 
1206
  odr_global.fcn = fcn;
 
1207
  Py_INCREF(fcn);
 
1208
  odr_global.fjacb = fjacb;
 
1209
  Py_XINCREF(fjacb);
 
1210
  odr_global.fjacd = fjacd;
 
1211
  Py_XINCREF(fjacd);
 
1212
  odr_global.pyBeta = (PyObject *) beta;
 
1213
  Py_INCREF(beta);
 
1214
  odr_global.extra_args = extra_args;
 
1215
  Py_XINCREF(extra_args);
 
1216
 
 
1217
  /* now call DODRC */
 
1218
  F_FUNC(dodrc,DODRC)(fcn_callback, &n, &m, &np, &nq, (double *)(beta->data),
 
1219
        (double *)(y->data), &ldy, (double *)(x->data), &ldx,
 
1220
        (double *)(we->data), &ldwe, &ld2we,
 
1221
        (double *)(wd->data), &ldwd, &ld2wd,
 
1222
        (int *)(ifixb->data), (int *)(ifixx->data), &ldifx,
 
1223
        &job, &ndigit, &taufac, &sstol, &partol, &maxit,
 
1224
        &iprint, &lunerr, &lunrpt,
 
1225
        (double *)(stpb->data), (double *)(stpd->data), &ldstpd,
 
1226
        (double *)(sclb->data), (double *)(scld->data), &ldscld,
 
1227
        (double *)(work->data), &lwork, (int *)(iwork->data), &liwork,
 
1228
        &info);
 
1229
 
 
1230
  result = gen_output(n, m, np, nq, ldwe, ld2we,
 
1231
                      beta, work, iwork, isodr, info, full_output);
 
1232
 
 
1233
  if (result == NULL)
 
1234
    PYERR(PyExc_RuntimeError, "could not generate output");
 
1235
 
 
1236
  if (lunerr != -1)
 
1237
    {
 
1238
      F_FUNC(dlunc,DLUNC)(&lunerr);
 
1239
    }
 
1240
  if (lunrpt != -1)
 
1241
    {
 
1242
      F_FUNC(dlunc,DLUNC)(&lunrpt);
 
1243
    }
 
1244
 
 
1245
  Py_DECREF(odr_global.fcn);
 
1246
  Py_XDECREF(odr_global.fjacb);
 
1247
  Py_XDECREF(odr_global.fjacd);
 
1248
  Py_DECREF(odr_global.pyBeta);
 
1249
  Py_XDECREF(odr_global.extra_args);
 
1250
 
 
1251
  odr_global.fcn = odr_global.fjacb = odr_global.fjacd = odr_global.pyBeta =
 
1252
    odr_global.extra_args = NULL;
 
1253
 
 
1254
  Py_DECREF(beta);
 
1255
  Py_DECREF(y);
 
1256
  Py_DECREF(x);
 
1257
  Py_DECREF(we);
 
1258
  Py_DECREF(wd);
 
1259
  Py_DECREF(ifixb);
 
1260
  Py_DECREF(ifixx);
 
1261
  Py_DECREF(stpb);
 
1262
  Py_DECREF(stpd);
 
1263
  Py_DECREF(sclb);
 
1264
  Py_DECREF(scld);
 
1265
  Py_DECREF(work);
 
1266
  Py_DECREF(iwork);
 
1267
 
 
1268
  return result;
 
1269
 
 
1270
fail:
 
1271
 
 
1272
 
 
1273
  if (lunerr != -1)
 
1274
    {
 
1275
      F_FUNC(dlunc,DLUNC)(&lunerr);
 
1276
    }
 
1277
  if (lunrpt != -1)
 
1278
    {
 
1279
      F_FUNC(dlunc,DLUNC)(&lunrpt);
 
1280
    }
 
1281
 
 
1282
  Py_XDECREF(beta);
 
1283
  Py_XDECREF(y);
 
1284
  Py_XDECREF(x);
 
1285
  Py_XDECREF(we);
 
1286
  Py_XDECREF(wd);
 
1287
  Py_XDECREF(ifixb);
 
1288
  Py_XDECREF(ifixx);
 
1289
  Py_XDECREF(stpb);
 
1290
  Py_XDECREF(stpd);
 
1291
  Py_XDECREF(sclb);
 
1292
  Py_XDECREF(scld);
 
1293
  Py_XDECREF(work);
 
1294
  Py_XDECREF(iwork);
 
1295
 
 
1296
  return NULL;
 
1297
}
 
1298
 
 
1299
static void check_args(int n, int m, int np, int nq,
 
1300
                       PyArrayObject * beta,
 
1301
                       PyArrayObject * y, int ldy,
 
1302
                       PyArrayObject * x, int ldx,
 
1303
                       PyArrayObject * we, int ldwe, int ld2we,
 
1304
                       PyArrayObject * wd, int ldwd, int ld2wd,
 
1305
                       PyArrayObject * ifixb, PyArrayObject * ifixx,
 
1306
                       int ldifx, int job, int ndigit, double taufac,
 
1307
                       double sstol, double partol, int maxit,
 
1308
                       PyArrayObject * stpb, PyArrayObject * stpd,
 
1309
                       int ldstpd, PyArrayObject * sclb,
 
1310
                       PyArrayObject * scld, int ldscld,
 
1311
                       PyArrayObject * work, int lwork,
 
1312
                       PyArrayObject * iwork, int liwork, int info)
 
1313
{
 
1314
  PyObject *printdict;
 
1315
 
 
1316
  printdict =
 
1317
    Py_BuildValue
 
1318
    ("{s:l,s:l,s:l,s:l,s:O,s:O,s:l,s:O,s:l,s:O,s:l,s:l,s:O,s:l,s:l,s:O,s:O,s:l,s:l,s:l,s:d,s:d,s:d,s:l,s:O,s:O,s:l,s:O,s:O,s:l,s:O,s:l,s:O,s:l,s:l}",
 
1319
     "n", n, "m", m, "np", np, "nq", nq, "beta", (PyObject *) beta, "y",
 
1320
     (PyObject *) y, "ldy", ldy, "x", (PyObject *) x, "ldx", ldx, "we",
 
1321
     (PyObject *) we, "ldwe", ldwe, "ld2we", ld2we, "wd", (PyObject *) wd,
 
1322
     "ldwd", ldwd, "ld2wd", ld2wd, "ifixb", (PyObject *) ifixb, "ifixx",
 
1323
     (PyObject *) ifixx, "ldifx", ldifx, "job", job, "ndigit", ndigit,
 
1324
     "taufac", taufac, "sstol", sstol, "partol", partol, "maxit", maxit,
 
1325
     "stpb", (PyObject *) stpb, "stpd", (PyObject *) stpd, "ldstpd", ldstpd,
 
1326
     "sclb", (PyObject *) sclb, "scld", (PyObject *) scld, "ldscld", ldscld,
 
1327
     "work", (PyObject *) work, "lwork", lwork, "iwork", (PyObject *) iwork,
 
1328
     "liwork", liwork, "info", info);
 
1329
  if (printdict == NULL)
 
1330
    {
 
1331
      PyErr_Print();
 
1332
      return;
 
1333
    }
 
1334
 
 
1335
  PyObject_Print(printdict, stdout, Py_PRINT_RAW);
 
1336
  printf("\n");
 
1337
  Py_XDECREF(printdict);
 
1338
}
 
1339
 
 
1340
static char odr__doc__[] =
 
1341
  "odr(fcn, beta0, y, x,\nwe=None, wd=None, fjacb=None, fjacd=None,\nextra_args=None, ifixx=None, ifixb=None, job=0, iprint=0,\nerrfile=None, rptfile=None, ndigit=0,\ntaufac=0.0, sstol=-1.0, partol=-1.0,\nmaxit=-1, stpb=None, stpd=None,\nsclb=None, scld=None, work=None, iwork=None,\nfull_output=0)";
 
1342
 
 
1343
static PyMethodDef methods[] = {
 
1344
  {"odr", (PyCFunction) odr, METH_VARARGS | METH_KEYWORDS, odr__doc__},
 
1345
  {NULL, NULL},
 
1346
};
 
1347
 
 
1348
PyMODINIT_FUNC init__odrpack(void)
 
1349
{
 
1350
  PyObject *m, *d;
 
1351
 
 
1352
  import_array();
 
1353
 
 
1354
  m = Py_InitModule("__odrpack", methods);
 
1355
  d = PyModule_GetDict(m);
 
1356
  odr_error = PyErr_NewException("odr.odrpack.odr_error", NULL, NULL);
 
1357
  odr_stop = PyErr_NewException("odr.odrpack.odr_stop", NULL, NULL);
 
1358
  PyDict_SetItemString(d, "odr_error", odr_error);
 
1359
  PyDict_SetItemString(d, "odr_stop", odr_stop);
 
1360
}