hocprax.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
/*
Hoc interface to praxis.
 
See praxis.cpp in the scopmath library about
tolerance (t0), maxstepsize (h0), and printmode (prin).
These are set with
    attr_praxis(t0, h0, prin)
 
Minimise an interpreted hoc function or compiled hoc function with
    double x[n]
    fit_praxis(n, "funname", &x[0], "afterquad statement")
where n is the number of variables on which funname depends and
x is a vector containing on entry a guess of the point of minimum
and on exit contains the estimated point of minimum.
(The third arg may be a Vector)
Funname will be called with 2 args to be retrieved by $1 and $&2[i].
afterquad statement (if the arg exists , will be executed
at the end of each main loop iteration (complete conjugate gradient search
followed by calculation of quadratic form)
eg.
 
double x[2]
attr_praxis(1e-5, 1, -1)
func f() {
    return $&0[0]^2 + 2*$&0[1]^4
}
fit_praxis(2, "f", &x)
 
After fit_praxis exits, (or on execution of afterquad stmt)
one can retrieve the i'th principal value
with pval_praxis(i). If you also want the i'th principal axis then use
double a[n]
pval = pval_praxis(i, &a)
also can use
pval = pval_praxis(i, Vector)
*/
 
#include <stdlib.h>
#include "hocdec.h"
#include "parse.hpp"
 
#if defined(__cplusplus)
extern "C" {
#endif
 
extern double praxis(double* t0,
                     double* machep,
                     double* h0,
                     long int nval,
                     long int* prin,
                     double* x,
                     double (*f)(double*, long int),
                     double* fmin,
                     char* after_quad);
extern double praxis_pval(int), *praxis_paxis(int);
extern int praxis_stop(int);
 
#if defined(__cplusplus)
}
#endif
 
extern int stoprun;
extern double chkarg(int, double, double);
 
extern void vector_resize(IvocVect*, int);
extern double* vector_vec(IvocVect*);
extern IvocVect* vector_arg(int);  // TODO: IvocVect?
extern IvocVect* vector_new2(IvocVect* vec);
extern void vector_delete(IvocVect* vec);
extern int vector_capacity(IvocVect*);
extern Object** vector_pobj(IvocVect* v);
extern int nrn_praxis_ran_index;
extern Object** hoc_objgetarg(int);
 
static double efun(double*, long int);
static Symbol* hoc_efun_sym;
 
static double tolerance, machep, maxstepsize;
static long int printmode;
/* for prax_pval to know the proper size, following has to be the value
at return of previous prax call
*/
static long int nvar;
 
double (*nrnpy_praxis_efun)(Object* pycallable, Object* hvec);
static Object* efun_py;
static Object* efun_py_arg;
static void* vec_py_save;
 
void stop_praxis(void) {
    int i = 1;
    if (ifarg(1)) {
        i = (int) chkarg(1, 0., 1e5);
    }
    i = praxis_stop(i);
    hoc_retpushx((double) i);
}
 
/* want to get best result if user does a stoprun */
static double minerr, *minarg; /* too bad this is not recursive */
 
void fit_praxis(void) {
    char* after_quad;
    int i;
    double err, fmin;
    double* px;
    /* allow nested calls to fit_praxis. I.e. store all the needed
       statics specific to this invocation with proper reference
       counting and then unref/destoy on exit from this invocation.
       Before the prax call save the statics from earlier invocation
       without increasing the
       ref count and on exit restore without decreasing the ref count.
    */
 
    /* save before setting statics, restore after prax */
    double minerrsav, *minargsav, maxstepsizesav, tolerancesav;
    long int printmodesav, nvarsav;
    Symbol* efun_sym_sav;
    Object *efun_py_save, *efun_py_arg_save;
    void* vec_py_save_save;
 
    /* store statics specified by this invocation */
    /* will be copied just before calling prax */
    double minerr_, *minarg_;
    long int nvar_;
    Symbol* efun_sym_;
    Object *efun_py_, *efun_py_arg_;
    void* vec_py_save_;
 
    minerr_ = 0.0;
    nvar_ = 0;
    minarg_ = NULL;
    efun_sym_ = NULL;
    efun_py_ = NULL;
    efun_py_arg_ = NULL;
    vec_py_save_ = NULL;
 
    fmin = 0.;
 
    if (hoc_is_object_arg(1)) {
        assert(nrnpy_praxis_efun);
        efun_py_ = *hoc_objgetarg(1);
        hoc_obj_ref(efun_py_);
        efun_py_arg_ = *vector_pobj(vector_arg(2));
        hoc_obj_ref(efun_py_arg_);
        vec_py_save_ = vector_new2(static_cast<IvocVect*>(efun_py_arg_->u.this_pointer));
        nvar_ = vector_capacity(static_cast<IvocVect*>(vec_py_save_));
        px = vector_vec(static_cast<IvocVect*>(vec_py_save_));
    } else {
        nvar_ = (int) chkarg(1, 0., 1e6);
        efun_sym_ = hoc_lookup(gargstr(2));
        if (!efun_sym_ || (efun_sym_->type != FUNCTION && efun_sym_->type != FUN_BLTIN)) {
            hoc_execerror(gargstr(2), "not a function name");
        }
 
        if (!hoc_is_pdouble_arg(3)) {
            IvocVect* vec = vector_arg(3);
            if (vector_capacity(vec) != nvar_) {
                hoc_execerror("first arg not equal to size of Vector", 0);
            }
            px = vector_vec(vec);
        } else {
            px = hoc_pgetarg(3);
        }
    }
    minarg_ = (double*) ecalloc(nvar_, sizeof(double));
 
    if (maxstepsize == 0.) {
        hoc_execerror("call attr_praxis first to set attributes", 0);
    }
    machep = 1e-15;
 
    if (ifarg(4)) {
        after_quad = gargstr(4);
    } else {
        after_quad = (char*) 0;
    }
 
    /* save the values set by earlier invocation */
    minerrsav = minerr;
    minargsav = minarg;
    tolerancesav = tolerance;
    maxstepsizesav = maxstepsize;
    printmodesav = printmode;
    nvarsav = nvar;
    efun_sym_sav = hoc_efun_sym;
    efun_py_save = efun_py;
    efun_py_arg_save = efun_py_arg;
    vec_py_save_save = vec_py_save;
 
 
    /* copy this invocation values to the statics */
    minerr = minerr_;
    minarg = minarg_;
    nvar = nvar_;
    hoc_efun_sym = efun_sym_;
    efun_py = efun_py_;
    efun_py_arg = efun_py_arg_;
    vec_py_save = vec_py_save_;
 
 
    minerr = 1e9;
    err = praxis(&tolerance, &machep, &maxstepsize, nvar, &printmode, px, efun, &fmin, after_quad);
    err = minerr;
    if (minerr < 1e9) {
        for (i = 0; i < nvar; ++i) {
            px[i] = minarg[i];
        }
    }
 
    /* restore the values set by earlier invocation */
    minerr = minerrsav;
    minarg = minargsav;
    tolerance = tolerancesav;
    maxstepsize = maxstepsizesav;
    printmode = printmodesav;
    nvar = nvar_; /* in case one calls prax_pval */
    hoc_efun_sym = efun_sym_sav;
    efun_py = efun_py_save;
    efun_py_arg = efun_py_arg_save;
    vec_py_save = vec_py_save_save;
 
    if (efun_py_) {
        double* px = vector_vec(static_cast<IvocVect*>(efun_py_arg_->u.this_pointer));
        for (i = 0; i < nvar_; ++i) {
            px[i] = minarg_[i];
        }
        hoc_obj_unref(efun_py_);
        hoc_obj_unref(efun_py_arg_);
        vector_delete(static_cast<IvocVect*>(vec_py_save_));
    }
    if (minarg_) {
        free(minarg_);
    }
    hoc_retpushx(err);
}
 
extern "C" void hoc_after_prax_quad(char* s) {
    efun(minarg, nvar);
    hoc_obj_run(s, hoc_thisobject);
}
 
void attr_praxis(void) {
    if (ifarg(2)) {
        tolerance = *getarg(1);
        maxstepsize = *getarg(2);
        printmode = (int) chkarg(3, 0., 3.);
        hoc_retpushx(0.);
    } else {
        int old = nrn_praxis_ran_index;
        if (ifarg(1)) {
            nrn_praxis_ran_index = (int) chkarg(1, 0., 1e9);
        }
        hoc_retpushx((double) old);
    }
}
 
void pval_praxis(void) {
    int i;
    i = (int) chkarg(1, 0., (double) (nvar - 1));
    if (ifarg(2)) {
        int j;
        double *p1, *p2;
        p1 = praxis_paxis(i);
        if (!hoc_is_pdouble_arg(2)) {
            IvocVect* vec = vector_arg(2);
            vector_resize(vec, nvar);
            p2 = vector_vec(vec);
        } else {
            p2 = hoc_pgetarg(2);
        }
        for (j = 0; j < nvar; ++j) {
            p2[j] = p1[j];
        }
    }
    hoc_retpushx(praxis_pval(i));
}
 
static double efun(double* v, long int n) {
    int i;
    double err;
    if (efun_py) {
        double* px = vector_vec(static_cast<IvocVect*>(efun_py_arg->u.this_pointer));
        for (i = 0; i < n; ++i) {
            px[i] = v[i];
        }
        err = nrnpy_praxis_efun(efun_py, efun_py_arg);
        for (i = 0; i < n; ++i) {
            v[i] = px[i];
        }
    } else {
        int i;
        hoc_pushx((double) n);
        hoc_pushpx(v);
        err = hoc_call_func(hoc_efun_sym, 2);
    }
    if (!stoprun && err < minerr) {
        minerr = err;
        for (i = 0; i < n; ++i) {
            minarg[i] = v[i];
        }
    }
    return err;
}