1 /**
2 * @file evaluateMath.c
3 * @brief Evaluates and outputs infix expressions
4 * @author Rainer Machne <raim@tbi.univie.ac.at>
5 * @author Ben Bornstein
6 * @author Michael Hucka
7 *
8 * Copyright 2003 Rainer Machne
9 *
10 * This library is free software; you can redistribute it and/or modify it
11 * under the terms of the GNU Lesser General Public License as published by
12 * the Free Software Foundation. A copy of the license agreement is
13 * provided in the file named "LICENSE.txt" included with this software
14 * distribution. It is also available online at
15 * http://sbml.org/software/libsbml/license.html
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with this library; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
20 */
21
22 #include <math.h>
23 #include <stdlib.h>
24 #include <stdio.h>
25 #include <stddef.h>
26 #include <string.h>
27
28 #include <sbml/util/util.h>
29 #include "util.h"
30
31 #include <sbml/math/FormulaFormatter.h>
32 #include <sbml/math/FormulaParser.h>
33
34
35 #define SQR(x) ((x)*(x))
36 #define SQRT(x) pow((x),(.5))
37
38 double
39 evalAST(ASTNode_t *n);
40
41
42 /**
43 * This program asks the user to enter an infix formula, translates it to
44 * an Abstract Syntax tree using the function:
45 *
46 * ASTNode_t *SBML_parseFormula(char *)
47 *
48 * evaluates the formula and returns the result. See comments for
49 * double evalAST(ASTNode_t *n) for further information.
50 */
51 int
main()52 main()
53 {
54 char *line;
55 #ifdef __BORLANDC__
56 unsigned long start, stop;
57 #else
58 unsigned long long start, stop;
59 #endif
60
61 ASTNode_t *n;
62 double result;
63
64 printf( "\n" );
65 printf( "This program evaluates math formulas in infix notation.\n" );
66 printf( "Typing 'enter' triggers evaluation.\n" );
67 printf( "Typing 'enter' on an empty line stops the program.\n" );
68 printf( "\n" );
69
70 while (1)
71 {
72 printf( "Enter an infix formula (empty line to quit):\n\n" );
73 printf( "> " );
74
75 if ( strlen(line = trim_whitespace(get_line( stdin ))) == 0 ) break;
76
77 n = SBML_parseFormula(line);
78
79 start = getCurrentMillis();
80 result = evalAST(n);
81 stop = getCurrentMillis();
82
83 printf("\n%s\n= %.10g\n\n", SBML_formulaToString(n), result);
84 printf("evaluation time: %llu ms\n\n", stop - start);
85
86 free(line);
87 ASTNode_free(n);
88 }
89
90 return 0;
91 }
92
93
94 /**
95 * The function evalAST(ASTNode_t) evaluates the formula of an
96 * Abstract Syntax Tree by simple recursion and returns the result
97 * as a double value.
98 *
99 * If variables (ASTNodeType_t AST_NAME) occur in the formula the user is
100 * asked to provide a numerical value. When evaluating ASTs within an SBML
101 * document or simulating an SBML model this node type includes parameters
102 * and variables of the model. Parameters should be retrieved from the
103 * SBML file, time and variables from current values of the simulation.
104 *
105 * Not implemented:
106 *
107 * - PIECEWISE, LAMBDA, and the SBML model specific functions DELAY and
108 * TIME and user-defined functions.
109 *
110 * - Complex numbers and/or checking for domains of trigonometric and root
111 * functions.
112 *
113 * - Checking for precision and rounding errors.
114 *
115 * The Nodetypes AST_TIME, AST_DELAY and AST_PIECEWISE default to 0. The
116 * SBML DELAY function and unknown functions (SBML user-defined functions)
117 * use the value of the left child (first argument to function) or 0 if the
118 * node has no children.
119 */
120 double
evalAST(ASTNode_t * n)121 evalAST(ASTNode_t *n)
122 {
123 int i;
124 double result;
125
126 int childnum = ASTNode_getNumChildren(n);
127 ASTNode_t **child = (ASTNode_t **) malloc(childnum * sizeof(ASTNode_t*));
128
129
130 for (i = 0; i < childnum; i++)
131 {
132 child[i] = ASTNode_getChild(n, i);
133 }
134
135 switch (ASTNode_getType(n))
136 {
137 case AST_INTEGER:
138 result = (double) ASTNode_getInteger(n);
139 break;
140
141 case AST_REAL:
142 result = ASTNode_getReal(n);
143 break;
144
145 case AST_REAL_E:
146 result = ASTNode_getReal(n);
147 break;
148
149 case AST_RATIONAL:
150 result = ASTNode_getReal(n);
151 break;
152
153 case AST_NAME:
154 {
155 char *l;
156 double var;
157 printf("\n-------------MESSAGE FROM EVALUATION FUNCTION-------------\n");
158 printf("Please enter a number for the variable!\n");
159 printf("If you do not enter a valid number (empty or characters), the \n");
160 printf("evaluation will proceed with a current internal value and the \n");
161 printf("result will make no sense.\n");
162 printf("%s=",ASTNode_getName(n));
163 l = get_line(stdin);
164 sscanf(l, "%lf", &var);
165 free(l);
166 printf("%s = %f\n", ASTNode_getName(n), var);
167 printf("-----------------------END MESSAGE--------------------------\n\n");
168 result = var;
169 }
170 break;
171
172 case AST_FUNCTION_DELAY:
173 printf("\n-------------MESSAGE FROM EVALUATION FUNCTION-------------\n");
174 printf("Delays can only be evaluated during a time series simulation.\n");
175 printf("The value of the first child (ie. the first argument to the function)\n");
176 printf("is used for this evaluation. If the function node has no children the\n");
177 printf("value defaults to 0.\n");
178 printf("-----------------------END MESSAGE--------------------------\n\n");
179 if(i>0)
180 result = evalAST(child[0]);
181 else
182 result = 0.0;
183 break;
184
185 case AST_NAME_TIME:
186 printf("\n-------------MESSAGE FROM EVALUATION FUNCTION-------------\n");
187 printf("The time can only be evaluated during a time series simulation.\n");
188 printf("The value of defaults to 0\n");
189 printf("-----------------------END MESSAGE--------------------------\n\n");
190 result = 0.0;
191 break;
192
193 case AST_CONSTANT_E:
194 /* exp(1) is used to adjust exponentiale to machine precision */
195 result = exp(1);
196 break;
197
198 case AST_CONSTANT_FALSE:
199 result = 0.0;
200 break;
201
202 case AST_CONSTANT_PI:
203 /* pi = 4 * atan 1 is used to adjust Pi to machine precision */
204 result = 4.*atan(1.);
205 break;
206
207 case AST_CONSTANT_TRUE:
208 result = 1.0;
209 break;
210
211 case AST_PLUS:
212 result = evalAST(child[0]) + evalAST(child[1]);
213 break;
214
215 case AST_MINUS:
216 if(childnum==1)
217 result = - (evalAST(child[0]));
218 else
219 result = evalAST(child[0]) - evalAST(child[1]);
220 break;
221
222 case AST_TIMES:
223 result = evalAST(child[0]) * evalAST(child[1]) ;
224 break;
225
226 case AST_DIVIDE:
227 result = evalAST(child[0]) / evalAST(child[1]);
228 break;
229
230 case AST_POWER:
231 result = pow(evalAST(child[0]),evalAST(child[1]));
232 break;
233
234 case AST_LAMBDA:
235 printf("\n-------------MESSAGE FROM EVALUATION FUNCTION-------------\n");
236 printf("This function is not implemented yet.\n");
237 printf("The value defaults to 0.\n");
238 printf("-----------------------END MESSAGE--------------------------\n\n");
239 result = 0.0;
240 break;
241
242 case AST_FUNCTION:
243 printf("\n-------------MESSAGE FROM EVALUATION FUNCTION-------------\n");
244 printf("This function is not known.\n");
245 printf("Within an SBML document new functions can be defined by the user or \n");
246 printf("application. The value of the first child (ie. the first argument to \n");
247 printf("the function) is used for this evaluation. If the function node has\n");
248 printf("no children the value defaults to 0.\n");
249 printf("-----------------------END MESSAGE--------------------------\n\n");
250 if(childnum>0)
251 result = evalAST(child[0]);
252 else
253 result = 0.0;
254 break;
255
256 case AST_FUNCTION_ABS:
257 result = (double) fabs(evalAST(child[0]));
258 break;
259
260 case AST_FUNCTION_ARCCOS:
261 result = acos(evalAST(child[0])) ;
262 break;
263
264 case AST_FUNCTION_ARCCOSH:
265 #ifndef WIN32
266 result = acosh(evalAST(child[0]));
267 #else
268 result = log(evalAST(child[0]) + SQR(evalAST(child[0]) * evalAST(child[0]) - 1.));
269 #endif
270 break;
271 case AST_FUNCTION_ARCCOT:
272 /* arccot x = arctan (1 / x) */
273 result = atan(1./ evalAST(child[0]));
274 break;
275
276 case AST_FUNCTION_ARCCOTH:
277 /* arccoth x = 1/2 * ln((x+1)/(x-1)) */
278 result = ((1./2.)*log((evalAST(child[0])+1.)/(evalAST(child[0])-1.)) );
279 break;
280
281 case AST_FUNCTION_ARCCSC:
282 /* arccsc(x) = Arctan(1 / sqrt((x - 1)(x + 1))) */
283 result = atan( 1. / SQRT( (evalAST(child[0])-1.)*(evalAST(child[0])+1.) ) );
284 break;
285
286 case AST_FUNCTION_ARCCSCH:
287 /* arccsch(x) = ln((1 + sqrt(1 + x^2)) / x) */
288 result = log((1.+SQRT((1+SQR(evalAST(child[0]))))) /evalAST(child[0]));
289 break;
290
291 case AST_FUNCTION_ARCSEC:
292 /* arcsec(x) = arctan(sqrt((x - 1)(x + 1))) */
293 result = atan( SQRT( (evalAST(child[0])-1.)*(evalAST(child[0])+1.) ) );
294 break;
295
296 case AST_FUNCTION_ARCSECH:
297 /* arcsech(x) = ln((1 + sqrt(1 - x^2)) / x) */
298 result = log((1.+pow((1-SQR(evalAST(child[0]))),0.5))/evalAST(child[0]));
299 break;
300
301 case AST_FUNCTION_ARCSIN:
302 result = asin(evalAST(child[0]));
303 break;
304 case AST_FUNCTION_ARCSINH:
305 #ifndef WIN32
306 result = asinh(evalAST(child[0]));
307 #else
308 result = log(evalAST(child[0]) + SQR(evalAST(child[0]) * evalAST(child[0]) + 1.));
309 #endif
310 break;
311
312 case AST_FUNCTION_ARCTAN:
313 result = atan(evalAST(child[0]));
314 break;
315 case AST_FUNCTION_ARCTANH:
316 #ifndef WIN32
317 result = atanh(evalAST(child[0]));
318 #else
319 result = log((1. / evalAST(child[0]) + 1.) / (1. / evalAST(child[0]) - 1.)) / 2.;
320 #endif
321 break;
322
323 case AST_FUNCTION_CEILING:
324 result = ceil(evalAST(child[0]));
325 break;
326
327 case AST_FUNCTION_COS:
328 result = cos(evalAST(child[0]));
329 break;
330 case AST_FUNCTION_COSH:
331 result = cosh(evalAST(child[0]));
332 break;
333
334 case AST_FUNCTION_COT:
335 /* cot x = 1 / tan x */
336 result = (1./tan(evalAST(child[0])));
337 break;
338 case AST_FUNCTION_COTH:
339 /* coth x = cosh x / sinh x */
340 result = cosh(evalAST(child[0])) / sinh(evalAST(child[0]));
341 break;
342 case AST_FUNCTION_CSC:
343 /* csc x = 1 / sin x */
344 result = (1./sin(evalAST(child[0])));
345 break;
346
347 case AST_FUNCTION_CSCH:
348 /* csch x = 1 / cosh x */
349 result = (1./cosh(evalAST(child[0])));
350 break;
351
352 case AST_FUNCTION_EXP:
353 result = exp(evalAST(child[0]));
354 break;
355
356 case AST_FUNCTION_FACTORIAL:
357 {
358 printf("\n-------------MESSAGE FROM EVALUATION FUNCTION-------------\n");
359 printf("The factorial is only implemented for integer values. If a floating\n");
360 printf("point number is passed, the floor value is used for calculation!\n");
361 printf("-----------------------END MESSAGE--------------------------\n\n");
362 i = (int)floor(evalAST(child[0]));
363 for(result=1;i>1;--i)
364 result *= i;
365 }
366 break;
367
368 case AST_FUNCTION_FLOOR:
369 result = floor(evalAST(child[0]));
370 break;
371
372 case AST_FUNCTION_LN:
373 result = log(evalAST(child[0]));
374 break;
375
376 case AST_FUNCTION_LOG:
377 result = log10(evalAST(child[0]));
378 break;
379
380 case AST_FUNCTION_PIECEWISE:
381 printf("\n-------------MESSAGE FROM EVALUATION FUNCTION-------------\n");
382 printf("This function is not implemented yet.\n");
383 printf("The value defaults to 0.\n");
384 printf("-----------------------END MESSAGE--------------------------\n\n");
385 result = 0.0;
386 break;
387
388 case AST_FUNCTION_POWER:
389 result = pow(evalAST(child[0]),evalAST(child[1]));
390 break;
391
392 case AST_FUNCTION_ROOT:
393 result = pow(evalAST(child[1]),(1./evalAST(child[0])));
394 break;
395
396 case AST_FUNCTION_SEC:
397 /* sec x = 1 / cos x */
398 result = 1./cos(evalAST(child[0]));
399 break;
400
401 case AST_FUNCTION_SECH:
402 /* sech x = 1 / sinh x */
403 result = 1./sinh(evalAST(child[0]));
404 break;
405
406 case AST_FUNCTION_SIN:
407 result = sin(evalAST(child[0]));
408 break;
409
410 case AST_FUNCTION_SINH:
411 result = sinh(evalAST(child[0]));
412 break;
413
414 case AST_FUNCTION_TAN:
415 result = tan(evalAST(child[0]));
416 break;
417
418 case AST_FUNCTION_TANH:
419 result = tanh(evalAST(child[0]));
420 break;
421
422 case AST_LOGICAL_AND:
423 result = (double) ((evalAST(child[0])) && (evalAST(child[1])));
424 break;
425
426 case AST_LOGICAL_NOT:
427 result = (double) (!(evalAST(child[0])));
428 break;
429
430 case AST_LOGICAL_OR:
431 result = (double) ((evalAST(child[0])) || (evalAST(child[1])));
432 break;
433
434 case AST_LOGICAL_XOR:
435 result = (double) ((!(evalAST(child[0])) && (evalAST(child[1])))
436 || ((evalAST(child[0])) && !(evalAST(child[1]))));
437 break;
438
439 case AST_RELATIONAL_EQ :
440 result = (double) ((evalAST(child[0])) == (evalAST(child[1])));
441 break;
442
443 case AST_RELATIONAL_GEQ:
444 result = (double) ((evalAST(child[0])) >= (evalAST(child[1])));
445 break;
446
447 case AST_RELATIONAL_GT:
448 result = (double) ((evalAST(child[0])) > (evalAST(child[1])));
449 break;
450
451 case AST_RELATIONAL_LEQ:
452 result = (double) ((evalAST(child[0])) <= (evalAST(child[1])));
453 break;
454
455 case AST_RELATIONAL_LT :
456 result = (double) ((evalAST(child[0])) < (evalAST(child[1])));
457 break;
458
459 default:
460 result = 0;
461 break;
462 }
463
464 free(child);
465
466 return result;
467 }
468