PWL functionality for B sources

This commit is contained in:
dwarning 2009-10-11 08:50:54 +00:00
parent 05116f0e0e
commit 8e8ecc784c
4 changed files with 175 additions and 3 deletions

View File

@ -51,6 +51,7 @@ typedef struct INPparseNode {
char *funcname; /* If INP_FUNCTION, name of function, */
int funcnum; /* ... one of PTF_*, */
double (*function)(); /* ... and pointer to the function. */
void *data; /* private data for certain functions, currently PTF_PWL */
} INPparseNode;
/* These are the possible types of nodes we can have in the parse tree. The
@ -94,6 +95,8 @@ typedef struct INPparseNode {
#define PTF_URAMP 20
/* MW. PTF_CIF - next function */
#define PTF_USTEP2 21
#define PTF_PWL 22
#define PTF_PWL_DERIVATIVE 23
/* The following things are used by the parser -- these are the token types the
@ -160,6 +163,8 @@ extern double PTtanh();
extern double PTustep();
/* MW. PTcif declaration */
extern double PTustep2();
extern double PTpwl();
extern double PTpwl_derivative();
extern double PTuramp();
extern double PTuminus();

View File

@ -74,7 +74,10 @@ PTeval(INPparseNode * tree, double gmin, double *res, double *vals)
err = PTeval(tree->left, gmin, &r1, vals);
if (err != OK)
return (err);
*res = (*tree->function) (r1);
if(tree->data == NULL)
*res = (*tree->function) (r1);
else
*res = (*tree->function) (r1, tree->data);
if (*res == HUGE) {
fprintf(stderr, "Error: %g out of range for %s\n",
r1, tree->funcname);

View File

@ -82,7 +82,9 @@ static struct func {
{ "uramp", PTF_URAMP, PTuramp } ,
{ "-", PTF_UMINUS, PTuminus },
/* MW. cif function added */
{ "u2", PTF_USTEP2, PTustep2}
{ "u2", PTF_USTEP2, PTustep2},
{ "pwl", PTF_PWL, PTpwl},
{ "pwl_derivative", PTF_PWL_DERIVATIVE, PTpwl_derivative}
} ;
#define NUM_FUNCS (sizeof (funcs) / sizeof (struct func))
@ -385,6 +387,15 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum)
arg1 = mkcon((double) - 1.0);
break;
case PTF_PWL: /* PWL(var, x1, y1, x2, y2, ... a const list) */
arg1 = mkf(PTF_PWL_DERIVATIVE, p->left);
arg1->data = p->data;
break;
case PTF_PWL_DERIVATIVE: /* d/dvar PWL(var, ...) */
arg1 = mkcon((double) 0.0);
break;
default:
fprintf(stderr, "Internal Error: bad function # %d\n",
p->funcnum);
@ -528,6 +539,8 @@ static INPparseNode *mkf(int type, INPparseNode * arg)
p->function = funcs[i].funcptr;
p->funcname = funcs[i].name;
p->data = NULL;
return (p);
}
@ -551,6 +564,7 @@ static int PTcheck(INPparseNode * p)
case PT_TIMES:
case PT_DIVIDE:
case PT_POWER:
case PT_COMMA:
return (PTcheck(p->left) && PTcheck(p->right));
default:
@ -574,7 +588,7 @@ static char prectable[11][11] = {
/* * */ {G, G, G, G, G, L, L, L, G, L, G},
/* / */ {G, G, G, G, G, L, L, L, G, L, G},
/* ^ */ {G, G, G, G, G, L, L, L, G, L, G},
/* u-*/ {G, G, G, G, G, G, G, L, G, L, R},
/* u-*/ {G, G, G, G, G, G, G, L, G, L, G},
/* ( */ {R, L, L, L, L, L, L, L, E, L, L},
/* ) */ {G, G, G, G, G, G, G, R, G, R, G},
/* v */ {G, G, G, G, G, G, G, G, G, R, G},
@ -743,6 +757,89 @@ static INPparseNode *mkbnode(int opnum, INPparseNode * arg1,
return (p);
}
/*
* prepare_PTF_PWL()
* for the PWL(expr, points...) function
* collect the given points, which are expected to be given
* literal constant
* strip them from the INPparseNode
* and pass them as an opaque struct alongside the
* INPparseNode for the PWL(expr) function call
*
* Note:
* the functionINPgetTree() is missing a recursive decending simplifier
* as a consequence we can meet a PTF_UMINUS->PTF_CONSTANT
* instead of a plain PTF_CONSTANT here
* of course this can get arbitrarily more complex
* for example PTF_TIMES -> PTF_CONSTANT, PTF_CONSTANT
* etc.
* currently we support only PFT_CONST and PTF_UMINUS->PTF_CONST
*/
#define Breakpoint do { __asm__ __volatile__ ("int $03"); } while(0)
static INPparseNode *prepare_PTF_PWL(INPparseNode *p)
{
INPparseNode *w;
struct pwldata { int n; double vals[0]; } *data;
int i;
if (p->funcnum != PTF_PWL) {
fprintf(stderr, "MY-INFO: %s, very unexpected\n", __FUNCTION__);
exit(1);
}
fprintf(stderr, "MY-INFO: %s building a PTF_PWL\n", __FUNCTION__);
i = 0;
for(w = p->left; w->type == PT_COMMA; w = w->left)
i++;
if (i<2 || (i%1)) {
fprintf(stderr, "Error: PWL(expr, points...) needs an even and >=2 number of constant args\n");
return (NULL);
}
data = (struct pwldata *)
MALLOC(sizeof(struct pwldata) + i*sizeof(double));
data->n = i;
for (w = p->left ; --i >= 0 ; w = w->left)
if (w->right->type == PT_CONSTANT) {
data->vals[i] = w->right->constant;
} else if (w->right->type == PT_FUNCTION &&
w->right->funcnum == PTF_UMINUS &&
w->right->left->type == PT_CONSTANT) {
data->vals[i] = - w->right->left->constant;
} else {
fprintf(stderr, "MY-ERROR: %s, not a constant\n", __FUNCTION__);
fprintf(stderr, " type = %d\n", w->right->type);
//Breakpoint;
fprintf(stderr, "Error: PWL(expr, points...) only *literal* points are supported\n");
return (NULL);
}
for (i = 0 ; i < data->n ; i += 2)
fprintf(stderr, " (%lf %lf)\n", data->vals[i], data->vals[i+1]);
for (i = 2 ; i < data->n ; i += 2)
if(data->vals[i-2] >= data->vals[i]) {
fprintf(stderr, "Error: PWL(expr, points...) the abscissa of points must be ascending\n");
return (NULL);
}
/* strip all but the first arg,
* and attach the rest as opaque data to the INPparseNode
*/
p->left = w;
p->data = (void *) data;
return (p);
}
static INPparseNode *mkfnode(char *fname, INPparseNode * arg)
{
int i;
@ -845,6 +942,10 @@ static INPparseNode *mkfnode(char *fname, INPparseNode * arg)
p->funcname = funcs[i].name;
p->funcnum = funcs[i].number;
p->function = funcs[i].funcptr;
p->data = NULL;
if(p->funcnum == PTF_PWL)
p = prepare_PTF_PWL(p);
}
return (p);

View File

@ -244,3 +244,66 @@ PTuminus(double arg)
return (- arg);
}
double
PTpwl(double arg, void *data)
{
struct pwldata { int n; double vals[0]; } *thing = data;
// fprintf(stderr, "%s called, data=%p, n=%d\n", __FUNCTION__, data, thing->n);
// for(i=0; i<thing->n; i++)
// fprintf(stderr, " %lf", thing->vals[i]);
// fprintf(stderr, "\n");
double y;
int k0 = 0;
int k1 = thing->n/2 - 1;
while(k1-k0 > 1) {
int k = (k0+k1)/2;
if(thing->vals[2*k] > arg)
k1 = k;
else
k0 = k;
}
y = thing->vals[2*k0+1] +
(thing->vals[2*k1+1] - thing->vals[2*k0+1]) *
(arg - thing->vals[2*k0]) / (thing->vals[2*k1] - thing->vals[2*k0]);
// fprintf(stderr, "foo: interval %lf %lf %lf -> %lf\n", thing->vals[2*k0], arg, thing->vals[2*k1], y);
return y;
}
double
PTpwl_derivative(double arg, void *data)
{
struct pwldata { int n; double vals[0]; } *thing = data;
// fprintf(stderr, "%s called, data=%p, n=%d\n", __FUNCTION__, data, thing->n);
// for(i=0; i<thing->n; i++)
// fprintf(stderr, " %lf", thing->vals[i]);
// fprintf(stderr, "\n");
double y;
int k0 = 0;
int k1 = thing->n/2 - 1;
while(k1-k0 > 1) {
int k = (k0+k1)/2;
if(thing->vals[2*k] > arg)
k1 = k;
else
k0 = k;
}
y =
(thing->vals[2*k1+1] - thing->vals[2*k0+1]) /
(thing->vals[2*k1] - thing->vals[2*k0]);
// fprintf(stderr, "bar: interval %lf %lf %lf -> %lf\n", thing->vals[2*k0], arg, thing->vals[2*k1], y);
return y;
}