[system] / trunk / pg / lib / Value / Formula.pm Repository: Repository Listing bbplugincoursesdistsnplrochestersystemwww

# View of /trunk/pg/lib/Value/Formula.pm

Sun Oct 16 03:37:17 2005 UTC (14 years, 1 month ago) by dpvc
File size: 16143 byte(s)
```In the past, when Value objects were inserted into strings, they would
automatically include parentheses so that if you had \$f equal to 1+x
and \$g equal to 1-x, then Formula("\$f/\$g") would mean (1+x)/(1-x)
rather than 1+(x/1)-x, which is what would happen as a straing string
substitution.

The problem is that this would also happen for real numbers, vectors,
and everything else, even when it wasn't necessary.  So if \$x=Real(3),
then "Let x = \$x" would be "Let x = (3)".

I have changed the behavior of the string concatenation for Value
objects so that parentheses are only added in a few cases: for
Formulas, Complex numbers, and Unions.  This makes the other Value
objects work more like regular variables in strings, but might cause
some problems with strings that are used as formulas.  For example, if
\$a = Real(-3), then "x + 2 \$a" will become "x + 2 -3", or "x-1" rather
than the expected "x - 6".  (The old approach would have made it "x +
2 (-3)" which would have worked properly).  For the most part, it is
easier to use something like "x + 2*\$a" or even "x" + 2*\$a in this
case, so the extra trouble of having to avoid parentheses when you
really meant to substitute the value into a string didn't seem worth
it.
```

```    1 ###########################################################################
2 #
3 #  Implements the Formula class.
4 #
5 package Value::Formula;
6 my \$pkg = 'Value::Formula';
7
8 my \$UNDEF = bless {}, "UNDEF";
9
10 use strict;
11 use vars qw(@ISA);
12 @ISA = qw(Parser Value);
13
16        '-'    => sub {shift->sub(@_)},
17        '*'    => sub {shift->mult(@_)},
18        '/'    => sub {shift->div(@_)},
19        '**'   => sub {shift->power(@_)},
20        '.'    => sub {shift->_dot(@_)},
21        'x'    => sub {shift->cross(@_)},
22        '<=>'  => sub {shift->compare(@_)},
23        'cmp'  => sub {shift->compare_string(@_)},
24        '~'    => sub {shift->call('conj',@_)},
25        'neg'  => sub {shift->neg},
26        'sin'  => sub {shift->call('sin',@_)},
27        'cos'  => sub {shift->call('cos',@_)},
28        'exp'  => sub {shift->call('exp',@_)},
29        'abs'  => sub {shift->call('abs',@_)},
30        'log'  => sub {shift->call('log',@_)},
31        'sqrt' => sub {shift->call('sqrt',@_)},
32       'atan2' => sub {shift->atan2(@_)},
33    'nomethod' => sub {shift->nomethod(@_)},
34          '""' => sub {shift->stringify(@_)};
35
36 #
37 #  Call Parser to make the new item, copying important
38 #    fields from the tree.  The Context can override the
39 #    Context()->{value}{Formula} setting to substitue a
40 #    different class to call for creating the formula.
41 #
42 sub new {
43   shift; my \$self = \$\$Value::context->{value}{Formula}->create(@_);
44   foreach my \$id ('open','close') {\$self->{\$id} = \$self->{tree}{\$id}}
45   return \$self;
46 }
47
48 #
49 #  Call Parser to creat the formula
50 #
51 sub create {shift; \$pkg->SUPER::new(@_)}
52
53 #
54 #  Create the new parser with no string
55 #    (we'll fill in its tree by hand)
56 #
57 sub blank {\$pkg->SUPER::new('')}
58
59 #
60 #  with() changes tree element as well
61 #    as the formula itself.
62 #
63 sub with {
64   my \$self = shift; my %hash = @_;
65   foreach my \$id (keys(%hash)) {
66     \$self->{tree}{\$id} = \$hash{\$id};
67     \$self->{\$id} = \$hash{\$id};
68   }
69   return \$self;
70 }
71
72 #
73 #  Get the type from the tree
74 #
75 sub typeRef {(shift)->{tree}->typeRef}
76 sub length {(shift)->{tree}->typeRef->{length}}
77
78 sub isZero {(shift)->{tree}{isZero}}
79 sub isOne {(shift)->{tree}{isOne}}
80
81 sub isSetOfReals {(shift)->{tree}->isSetOfReals}
82 sub canBeInUnion {(shift)->{tree}->canBeInUnion}
83
84 ############################################
85 #
86 #  Create a BOP from two operands
87 #
88 #  Get the context and variables from the left and right operands
89 #    if they are formulas
90 #  Make them into Value objects if they aren't already.
91 #  Convert '+' to union for intervals or unions.
92 #  Make a new BOP with the two operands.
93 #  Record the variables.
94 #  Evaluate the formula if it is constant.
95 #
96 sub bop {
97   my (\$bop,\$l,\$r,\$flag) = @_;
98   my \$call = \$\$Value::context->{method}{\$bop};
99   if (\$l->promotePrecedence(\$r)) {return \$r->\$call(\$l,!\$flag)}
100   if (\$flag) {my \$tmp = \$l; \$l = \$r; \$r = \$tmp}
101   my \$formula = \$pkg->blank; my \$parser = \$formula->{context}{parser};
102   if (ref(\$r) eq \$pkg) {
103     \$formula->{context} = \$r->{context};
104     \$r = \$r->{tree}->copy(\$formula);
105   }
106   if (ref(\$l) eq \$pkg) {
107     \$formula->{context} = \$l->{context};
108     \$l = \$l->{tree}->copy(\$formula);
109   }
110   \$l = \$pkg->new(\$l) if (!ref(\$l) && Value::getType(\$formula,\$l) eq "unknown");
111   \$r = \$pkg->new(\$r) if (!ref(\$r) && Value::getType(\$formula,\$r) eq "unknown");
112   \$l = \$parser->{Value}->new(\$formula,\$l) unless ref(\$l) =~ m/^Parser::/;
113   \$r = \$parser->{Value}->new(\$formula,\$r) unless ref(\$r) =~ m/^Parser::/;
114   \$bop = 'U' if \$bop eq '+' &&
115     (\$l->type =~ m/Interval|Set|Union/ || \$r->type =~ m/Interval|Set|Union/);
116   \$formula->{tree} = \$parser->{BOP}->new(\$formula,\$bop,\$l,\$r);
117   \$formula->{variables} = \$formula->{tree}->getVariables;
118 #  return \$formula->eval if \$formula->{isConstant};
119   return \$formula;
120 }
121
123 sub sub   {bop('-',@_)}
124 sub mult  {bop('*',@_)}
125 sub div   {bop('/',@_)}
126 sub power {bop('**',@_)}
127 sub cross {bop('><',@_)}
128
129 #
130 #  Make dot work for vector operands
131 #
132 sub _dot   {
133   my (\$l,\$r,\$flag) = @_;
134   if (\$l->promotePrecedence(\$r)) {return \$r->_dot(\$l,!\$flag)}
135   return bop('.',@_) if \$l->type eq 'Vector' &&
136      Value::isValue(\$r) && \$r->type eq 'Vector';
137   \$l->SUPER::_dot(\$r,\$flag);
138 }
139
140 sub pdot {'('.(shift->stringify).')'}
141
142 #
143 #  Call the Parser::Function call function
144 #
145 sub call {
146   my \$self = shift; my \$name = shift;
147   Parser::Function->call(\$name,\$self);
148 }
149
150 ############################################
151 #
152 #  Form the negation of a formula
153 #
154 sub neg {
155   my \$self = shift;
156   my \$formula = \$self->blank;
157   \$formula->{context} = \$self->{context};
158   \$formula->{variables} = \$self->{variables};
159   \$formula->{tree} = \$formula->{context}{parser}{UOP}->new(\$formula,'u-',\$self->{tree}->copy(\$formula));
160   return \$formula;
161 }
162
163 #
164 #  Form the function atan2 function call on two operands
165 #
166 sub atan2 {
167   my (\$l,\$r,\$flag) = @_;
168   if (\$l->promotePrecedence(\$r)) {return \$r->atan2(\$l,!\$flag)}
169   if (\$flag) {my \$tmp = \$l; \$l = \$r; \$r = \$tmp}
170   Parser::Function->call('atan2',\$l,\$r);
171 }
172
173 ############################################
174 #
175 #  Compare two functions for equality
176 #
177 sub compare {
178   my (\$l,\$r,\$flag) = @_; my \$self = \$l;
179   if (\$l->promotePrecedence(\$r)) {return \$r->compare(\$l,!\$flag)}
180   \$r = Value::Formula->new(\$r) unless Value::isFormula(\$r);
181   Value::Error("Functions from different contexts can't be compared")
182     unless \$l->{context} == \$r->{context};
183
184   #
185   #  Get the test points and evaluate the functions at those points
186   #
187   ##  FIXME: Check given points for consistency
188   my \$points  = \$l->{test_points} || \$l->createRandomPoints(undef,\$l->{test_at});
189   my \$lvalues = \$l->{test_values} || \$l->createPointValues(\$points,1,1);
190   my \$rvalues = \$r->createPointValues(\$points,0,1,\$l->{checkUndefinedPoints});
191   #
192   # Note: \$l is bigger if \$r can't be evaluated at one of the points
193   return 1 unless \$rvalues;
194
195   my (\$i, \$cmp);
196
197   #
199   #    Get the tolerances, and check each adapted point relative
200   #    to the ORIGINAL correct answer.  (This will have to be
201   #    fixed if we ever do adaptive parameters for non-real formulas)
202   #
205     my \$tolerance  = \$self->getFlag('tolerance',1E-4);
206     my \$isRelative = \$self->getFlag('tolType','relative') eq 'relative';
207     my \$zeroLevel  = \$self->getFlag('zeroLevel',1E-14);
208     foreach \$i (0..scalar(@{\$lvalues})-1) {
209       my \$tol = \$tolerance;
210       my (\$lv,\$rv,\$av) = (\$lvalues->[\$i]->value,\$rvalues->[\$i]->value,\$avalues->[\$i]->value);
211       \$tol *= abs(\$lv) if \$isRelative && abs(\$lv) > \$zeroLevel;
212       return \$rv <=> \$av unless abs(\$rv - \$av) < \$tol;
213     }
214     return 0;
215   }
216
217   #
218   #  Look through the two lists of values to see if they are equal.
219   #  If not, return the comparison of the first unequal value
220   #    (not good for < and >, but OK for ==).
221   #
222   my \$domainError = 0;
223   foreach \$i (0..scalar(@{\$lvalues})-1) {
224     if (ref(\$lvalues->[\$i]) eq 'UNDEF' ^ ref(\$rvalues->[\$i]) eq 'UNDEF') {\$domainError = 1; next}
225     \$cmp = \$lvalues->[\$i] <=> \$rvalues->[\$i];
226     return \$cmp if \$cmp;
227   }
228   \$l->{domainMismatch} = \$domainError;  # return this value
229 }
230
231 #
232 #  Create the value list from a given set of test points
233 #
234 sub createPointValues {
235   my \$self = shift;
236   my \$points = shift || \$self->{test_points} || \$self->createRandomPoints;
237   my \$showError = shift; my \$cacheResults = shift;
238   my @vars   = \$self->{context}->variables->variables;
239   my @params = \$self->{context}->variables->parameters;
240   my @zeros  = (0) x scalar(@params);
241   my \$f = \$self->{f}; \$f = \$self->{f} = \$self->perlFunction(undef,[@vars,@params]) unless \$f;
242   my \$checkUndef = scalar(@params) == 0 && (shift || \$self->getFlag('checkUndefinedPoints',0));
243
244   my \$values = []; my \$v;
245   foreach my \$p (@{\$points}) {
246     \$v = eval {&\$f(@{\$p},@zeros)};
247     if (!defined(\$v) && !\$checkUndef) {
248       return unless \$showError;
249       Value::Error("Can't evaluate formula on test point (%s)",join(',',@{\$p}));
250     }
251     push @{\$values}, (defined(\$v)? Value::makeValue(\$v): \$UNDEF);
252   }
253   if (\$cacheResults) {
254     \$self->{test_points} = \$points;
255     \$self->{test_values} = \$values;
256   }
257   return \$values;
258 }
259
260 #
261 #  Create the adapted value list for the test points
262 #
264   my \$self = shift;
265   my \$points = shift || \$self->{test_points} || \$self->createRandomPoints;
266   my \$showError = shift;
267   my @vars   = \$self->{context}->variables->variables;
268   my @params = \$self->{context}->variables->parameters;
269   my \$f = \$self->{f}; \$f = \$self->{f} = \$self->perlFunction(undef,[@vars,@params]) unless \$f;
270
271   my \$values = []; my \$v;
273   foreach my \$p (@{\$points}) {
275     if (!defined(\$v)) {
276       return unless \$showError;
277       Value::Error("Can't evaluate formula on test point (%s) with parameters (%s)",
279     }
280     push @{\$values}, Value::makeValue(\$v);
281   }
283 }
284
285 #
286 #  Create a list of random points, making sure that the function
287 #  is defined at the given points.  Error if we can't find enough.
288 #
289 sub createRandomPoints {
290   my \$self = shift;
291   my (\$num_points,\$include) = @_; my \$cacheResults = !defined(\$num_points);
292   \$num_points = int(\$self->getFlag('num_points',5)) unless defined(\$num_points);
293   \$num_points = 1 if \$num_points < 1;
294
295   my @vars   = \$self->{context}->variables->variables;
296   my @params = \$self->{context}->variables->parameters;
297   my @limits = \$self->getVariableLimits(@vars);
298   my @make   = \$self->getVariableTypes(@vars);
299   my @zeros  = (0) x scalar(@params);
300   my \$f = \$self->{f}; \$f = \$self->{f} = \$self->perlFunction(undef,[@vars,@params]) unless \$f;
301   my \$seedRandom = \$self->{context}->flag('random_seed')? 'PGseedRandom' : 'seedRandom';
302   my \$getRandom  = \$self->{context}->flag('random_seed')? 'PGgetRandom'  : 'getRandom';
303   my \$checkUndef = scalar(@params) == 0 && \$self->getFlag('checkUndefinedPoints',0);
304   my \$max_undef  = \$self->getFlag('max_undefined',\$num_points);
305
306   \$self->\$seedRandom;
307   my \$points = []; my \$values = []; my \$num_undef = 0;
308   if (\$include) {
309     push(@{\$points},@{\$include});
310     push(@{\$values},@{\$self->createPointValues(\$include,1,\$cacheResults,\$self->{checkundefinedPoints})});
311   }
312   my (@P,@p,\$v,\$i); my \$k = 0;
313   while (scalar(@{\$points}) < \$num_points+\$num_undef && \$k < 10) {
314     @P = (); \$i = 0;
315     foreach my \$limit (@limits) {
316       @p = (); foreach my \$I (@{\$limit}) {push @p, \$self->\$getRandom(@{\$I})}
317       push @P, \$make[\$i++]->make(@p);
318     }
319     \$v = eval {&\$f(@P,@zeros)};
320     if (!defined(\$v)) {
321       if (\$checkUndef && \$num_undef < \$max_undef) {
322   push @{\$points}, [@P];
323   push @{\$values}, \$UNDEF;
324   \$num_undef++;
325       }
326       \$k++;
327     } else {
328       push @{\$points}, [@P];
329       push @{\$values}, Value::makeValue(\$v);
330       \$k = 0; # reset count when we find a point
331     }
332   }
333
334   Value::Error("Can't generate enough valid points for comparison") if \$k;
335   return (\$points,\$values) unless \$cacheResults;
336   \$self->{test_values} = \$values;
337   \$self->{test_points} = \$points;
338 }
339
340 #
341 #  Get the array of variable limits
342 #
343 sub getVariableLimits {
344   my \$self = shift;
345   my \$userlimits = \$self->{limits};
346   if (defined(\$userlimits)) {
347     \$userlimits = [[[-\$userlimits,\$userlimits]]] unless ref(\$userlimits) eq 'ARRAY';
348     \$userlimits = [\$userlimits] unless ref(\$userlimits->[0]) eq 'ARRAY';
349     \$userlimits = [\$userlimits] if scalar(@_) == 1 && ref(\$userlimits->[0][0]) ne 'ARRAY';
350     foreach my \$I (@{\$userlimits}) {\$I = [\$I] unless ref(\$I->[0]) eq 'ARRAY'};
351   }
352   \$userlimits = [] unless \$userlimits; my @limits;
353   my \$default;  \$default = \$userlimits->[0][0] if defined(\$userlimits->[0]);
354   \$default = \$default || \$self->{context}{flags}{limits} || [-2,2];
355   my \$granularity = \$self->getFlag('granularity',1000);
356   my \$resolution = \$self->getFlag('resolution');
357   my \$i = 0;
358   foreach my \$x (@_) {
359     my \$def = \$self->{context}{variables}{\$x};
360     my \$limit = \$userlimits->[\$i++] || \$def->{limits} || [];
361     \$limit = [\$limit] if defined(\$limit->[0]) && ref(\$limit->[0]) ne 'ARRAY';
362     push(@{\$limit},\$limit->[0] || \$default) while (scalar(@{\$limit}) < \$def->{type}{length});
363     pop(@{\$limit}) while (scalar(@{\$limit}) > \$def->{type}{length});
365   }
366   return @limits;
367 }
368
369 #
370 #  Add the granularity to the limit intervals
371 #
373   my \$self = shift; my \$limit = shift; my \$def = shift;
374   my \$granularity = shift; my \$resolution = shift;
375   \$granularity = \$def->{granularity} || \$granularity;
376   \$resolution = \$def->{resolution} || \$resolution;
377   foreach my \$I (@{\$limit}) {
378     my (\$a,\$b,\$n) = @{\$I}; \$b = -\$a unless defined \$b;
379     \$I = [\$a,\$b,(\$n || \$resolution || abs(\$b-\$a)/\$granularity)];
380   }
381   return \$limit;
382 }
383
384 #
385 #  Get the routines to make the coordinates of the points
386 #
387 sub getVariableTypes {
388   my \$self = shift;
389   my @make;
390   foreach my \$x (@_) {
391     my \$type = \$self->{context}{variables}{\$x}{type};
392     if (\$type->{name} eq 'Number') {
393       push @make,(\$type->{length} == 1)? 'Value::Formula::number': 'Value::Complex';
394     } else {
395       push @make, "Value::\$type->{name}";
396     }
397   }
398   return @make;
399 }
400
401 #
402 #  Fake object for making reals (rather than use overhead of Value::Real)
403 #
404 sub Value::Formula::number::make {shift; shift}
405
406 #
407 #  Find adaptive parameters, if any
408 #
410   my \$l = shift; my \$r = shift;
411   my @params = @_; my \$d = scalar(@params);
412   return 0 if \$d == 0; return 0 unless \$l->usesOneOf(@params);
413   \$l->Error("Adaptive parameters can only be used for real-valued functions")
414     unless \$l->{tree}->isRealNumber;
415   #
416   #  Get coefficient matrix of adaptive parameters
417   #  and value vector for linear system
418   #
419   my (\$p,\$v) = \$l->createRandomPoints(\$d);
420   my @P = (0) x \$d; my (\$f,\$F) = (\$l->{f},\$r->{f});
421   my @A = (); my @b = ();
422   foreach my \$i (0..\$d-1) {
423     my @a = (); my @p = @{\$p->[\$i]};
424     foreach my \$j (0..\$d-1) {
425       \$P[\$j] = 1; push(@a,&\$f(@p,@P)-\$v->[\$i]);
426       \$P[\$j] = 0;
427     }
428     push @A, [@a]; push @b, [&\$F(@p,@P)-\$v->[\$i]];
429   }
430   #
431   #  Use MatrixReal1.pm to solve system of linear equations
432   #
433   my \$M = MatrixReal1->new(\$d,\$d); \$M->[0] = \@A;
434   my \$B = MatrixReal1->new(\$d,1);  \$B->[0] = \@b;
435   (\$M,\$B) = \$M->normalize(\$B);
436   \$M = \$M->decompose_LR;
437   if ((\$d,\$B,\$M) = \$M->solve_LR(\$B)) {
438     if (\$d == 0) {
439       #
440       #  Get parameter values and recompute the points using them
441       #
442       my @a; my \$i = 0; my \$max = Value::Real->new(\$l->getFlag('max_adapt',1E8));
443       foreach my \$row (@{\$B->[0]}) {
444   if (abs(\$row->[0]) > \$max) {
445     \$l->Error("Constant of integration is too large: %s\n(maximum allowed is %s)",
446         \$row->[0]->string,\$max->string) if (\$params[\$i] eq 'C0');
447     \$l->Error("Adaptive constant is too large: %s = %s\n(maximum allowed is %s)",
448         \$params[\$i],\$row->[0]->string,\$max->string);
449   }
450   push @a, \$row->[0]; \$i++;
451       }
452       \$l->{parameters} = [@a];
454       return 1;
455     }
456   }
457   \$l->Error("Can't solve for adaptive parameters");
458 }
459
460 sub usesOneOf {
461   my \$self = shift;
462   foreach my \$x (@_) {return 1 if \$self->{variables}{\$x}}
463   return 0;
464 }
465
466 ##
467 ##  debugging routine
468 ##
469 #sub main::Format {
470 #  my \$v = scalar(@_) > 1? [@_]: shift;
471 #  \$v = [%{\$v}] if ref(\$v) eq 'HASH';
472 #  return \$v unless ref(\$v) eq 'ARRAY';
473 #  my @V; foreach my \$x (@{\$v}) {push @V, main::Format(\$x)}
474 #  return '['.join(",",@V).']';
475 #}
476
477 #
478 #  Random number generator  (replaced by Value::WeBWorK.pm)
479 #
480 sub seedRandom {srand}
481 sub getRandom {
482   my \$self = shift;
483   my (\$m,\$M,\$n) = @_; \$n = 1 unless \$n;
484   return \$m + \$n*int(rand()*(int((\$M-\$m)/\$n)+1));
485 }
486
487 ############################################
488 #
489 #  Check if the value of a formula is constant
490 #    (could use shift->{tree}{isConstant}, but I don't trust it)
491 #
492 sub isConstant {
493   my @vars = (%{shift->{variables}});
494   return scalar(@vars) == 0;
495 }
496
497 ############################################
498 #
499 #  Provide output formats
500 #
501 sub stringify {
502   my \$self = shift;
503   return \$self->TeX if \$\$Value::context->flag('StringifyAsTeX');
504   \$self->string;
505 }
506
507 ###########################################################################
508
509 1;
```