[system] / trunk / pg / macros / parserImplicitEquation.pl Repository: Repository Listing bbplugincoursesdistsnplrochestersystemwww

View of /trunk/pg/macros/parserImplicitEquation.pl

Thu Jun 28 22:31:59 2007 UTC (12 years, 7 months ago) by dpvc
File size: 13371 byte(s)
```Updated to use the new context methods.
```

```    1 loadMacros("Parser.pl");
2
3 sub _parserImplicitEquation_init {}; # don't reload this file
4
6
7 ######################################################################
8 #
9 #  This is a MathObject class that implements an answer checker for
10 #  implicitly defined equations.  The checker looks for the zeros of
11 #  the equation and tests that the student and professor equations
12 #  both have the same solutions.
13 #
14 #  This type of check is very subtle, and there are important issues
15 #  that you may have to take into account.  The solutions to the
16 #  equations are found numerically, and so they will not be exact;
17 #  that means that there are tolerances that may need to be adjusted
18 #  for your particular equation.  Also, it is always possible for the
19 #  student to represent the function in a form that will exceed those
20 #  tolerances, and so be marked as incorrect.  The answer checker
21 #  attempts to set the parameters based on the values of the functions
22 #  involved, but this may not work perfectly.
23 #
24 #  The method used to locate the solutions of A=B is by finding zeros
25 #  of A-B, and it requires this function to take on both positive and
26 #  negative values (that is, it can only find transverse intersections
27 #  of the surface A-B=0 and the plane at height 0).  For example, even
28 #  though the solutions of (x^2+y^1-1)^2=0 form a circle, the
29 #  algorithm will fail to find any solutions for this equation.
30 #
31 #  In order to locate the zeros, you may need to change the limits so
32 #  that they include regions where the function is both positive and
33 #  negative (see below).  The algorithm will avoid discontinuities, so
34 #  you can specify things like x-y=1/(x+y) rather than x^2-y^2=1.
35 #
36 #  Because the solutions are found using a random search, it is
37 #  possible the randomly chosen starting points don't locate enough
38 #  zeros for the function, in which case the check will fail.  This
39 #  can happen for both the professor's function and the student's,
40 #  since zeros are found for both.  This means that a correct answer
41 #  can sometimes be marked incorrect depending on the random points
42 #  chosen initially.  These points also affect the values selected for
43 #  the tolerances used to determine when a function's value is zero,
44 #  and so can affect whether the student's function is marked as
45 #  correct or not.
46 #
47 #  If an equation has several components or branches, it is possible
48 #  that the random location of solutions will not find zeros on some
49 #  of the branches, and so might incorrectly mark as correct an
50 #  equation that only is zero on one of the components.  For example,
51 #  x^2-y^2=0 has solutions along the lines y=x and y=-x, so it is
52 #  possible that x-y=0 or x+y=0 will be marked as correct if the
53 #  random points are unluckily chosen.  One way to reduce this problem
54 #  is to increase the number of solutions that are required (by
55 #  setting the ImplicitPoints flag in the Context).  Another is to
56 #  specify the solutions yourself, so that you are sure there are
57 #  points on each component.
58 #
59 #  These problems should be rare, and the values for the various
60 #  parameters have been set in an attempt to minimize the possibility
61 #  of these errors, but they can occur, and you should be aware of
62 #  them, and their possible solutions.
63 #
64 #
65 #  Usage examples:
66 #
67 #     Context("ImplicitEquation");
68 #     \$f = ImplicitEquation("x^2 = cos(y)");
69 #     \$f = ImplicitEquation("x^2 - 2y^2 = 5",limits=>[[-3,3],[-2,2]]);
70 #     \$f = ImplicitEquation("x=1/y",tolerance=>.0001);
71 #
72 #  Then use
73 #
74 #     ANS(\$f->cmp);
75 #
76 #  to get the answer checker for \$f.
77 #
78 #  There are a number of Context flags that control the answer checker.
79 #  These include:
80 #
81 #    ImplicitPoints => 7                   (the number of solutions to test)
82 #    ImplicitTolerance => 1E-6             (relative tolerance value for when
83 #                                           the tested function is zero)
84 #    ImplicitAbsoluteMinTolerance => 1E-3  (the minimum tolerance allowed)
85 #    ImplicitAbsoluteMaxTolerance => 1E-3  (the maximum tolerance allowed)
86 #    ImplicitPointTolerance => 1E-9        (relative tolerance for how close
87 #                                           the solution point must be to an
88 #                                           actual solution)
89 #    BisectionTolerance => .01             (extra factor used for the tolerance
90 #                                           when finding the solutions)
91 #    BisectionCutoff => 40                 (maximum number of bisections to
92 #                                           perform when looking for a solution)
93 #
94 #  You may set any of these using Context()->flags->set(...).
95 #
96 #  In addition to the Context flags, you can set some values within
97 #  the ImplicitEquation itself:
98 #
99 #    tolerance                (the absolute tolerance for zeros of the function)
100 #    bisect_tolerance         (the tolerance used when searching for zeros)
101 #    point_tolerance          (the absolute tolerance for how close to an
102 #                              actual solution the located solution must be)
103 #    limits                   (the domain to use for the function; see the
104 #                              documentation for the Formula object)
105 #    solutions                (a reference to an array of references to arrays
106 #                              that contain the coordinates of the points
107 #                              that are the solutions of the equation)
108 #
109 #  These can be set in the in the ImplicitEquation() call that creates
110 #  the object, as in the examples below:
111 #
112 #  For example:
113 #
114 #    \$f = ImplicitEquation("x^2-y^2=0",
115 #            solutions => [[0,0],[1,1],[-1,1],[-1,-1],[1,-1]],
116 #            tolerance => .001
117 #         );
118 #
119 #
120 #    \$f = ImplicitEquation("xy=5",limits=>[-3,3]);
121 #
122 #  The limits value can be set globally within the Context, if you wish,
123 #  and the others can be controlled by the Context flags discussed
124 #  above.
125 #
126 ######################################################################
127
128 =cut
129
130 #
131 #  Set up the context for ImplicitEquations and activate it
132 #
133 \$context{ImplicitEquation} = Parser::Context->getCopy(undef,"Numeric");
134 \$context{ImplicitEquation}->variables->are(x=>'Real',y=>'Real');
135 \$context{ImplicitEquation}{precedence}{ImplicitEquation} = Context()->{precedence}{special};
136 Parser::BOP::equality->Allow(\$context{ImplicitEquation});
137 \$context{ImplicitEquation}->flags->set(
138   ImplicitPoints => 10,
139   ImplicitTolerance => 1E-6,
140   ImplicitAbsoluteMinTolerance => 1E-3,
141   ImplicitAbsoluteMaxTolerance => 1,
142   ImplicitPointTolerance => 1E-9,
143   BisectionTolerance => .01,
144   BisectionCutoff => 40,
145 );
146
147 Context("ImplicitEquation");
148
149 #
150 #  Syntactic sugar for creating ImplicitEquations
151 #
152 sub ImplicitEquation {ImplicitEquation->new(@_)}
153
154 #
155 #  Create the ImplicitEquation package
156 #
157 package ImplicitEquation;
158 our @ISA = qw(Value::Formula);
159
160 sub new {
161   my \$self = shift; my \$class = ref(\$self) || \$self;
162   my \$context = (Value::isContext(\$_[0]) ? shift : \$self->context);
163   my \$f = shift; return \$f if ref(\$f) eq \$class;
164   \$f = main::Formula(\$f);
165   Value::Error("Your formula doesn't look like an implicit equation")
166     unless \$f->type eq 'Equality';
167   my \$F = (\$context->Package("Formula")->new(\$context,\$f->{tree}{lop}) -
168      \$context->Package("Formula")->new(\$context,\$f->{tree}{rop}))->reduce;
169   \$F = \$context->Package("Formula")->new(\$F) unless Value::isFormula(\$F);
170   Value::Error("Your equation must be real-valued") unless \$F->isRealNumber;
171   Value::Error("Your equation should not be constant") if \$F->isConstant;
173     if (\$F->usesOneOf(\$context->variables->parameters));
174   \$F = bless \$F, \$class;
175   my %options = (@_);  # user can supply limits, tolerance, etc.
176   foreach my \$id (keys %options) {\$F->{\$id} = \$options{\$id}}
177   \$F->{F} = \$f; \$F->{isValue} = \$F->{isFormula} = 1;
178   \$F->createPoints unless \$F->{solutions};
179   return \$F;
180 }
181
182 #
183 #  Override the comparison method.
184 #
185 #  Turn the right-hand equation into an ImplicitEquation.  This creates
186 #  the test points (i.e., finds the solution points).  Then check
187 #  the professor's function on the student's test points and the
188 #  student's function on the professor's test points.
189 #
190
191 sub compare {
192   my (\$l,\$r) = @_; my \$self = \$l; my \$tolerance;
193   my @params; @params = (limits=>\$l->{limits}) if \$l->{limits};
194   \$r = ImplicitEquation->new(\$r,@params);
195   Value::Error("Functions from different contexts can't be compared")
196     unless \$l->{context} == \$r->{context};
197
198   #
199   #  They are not equal if couldn't get solutions for one of them
200   #
201   return 1 unless \$l->{solutions} && \$r->{solutions};
202
203   #
204   #  Test the right-hand function on the solutions of the left-hand one
205   #  and vice-versa
206   #
207   my \$rzeros = \$r->createPointValues(\$l->{solutions});
208   my \$lzeros = \$l->createPointValues(\$r->{solutions});
209   return 1 unless \$lzeros && \$rzeros;
210
211   #
212   #  Check that the values are, in fact, zeros
213   #
214   \$tolerance = \$r->getFlag('tolerance',1E-3);
215   foreach my \$v (@{\$rzeros}) {return 1 unless abs(\$v) < \$tolerance}
216   \$tolerance = \$l->getFlag('tolerance',1E-3);
217   foreach my \$v (@{\$lzeros}) {return 1 unless abs(\$v) < \$tolerance}
218
219   return 0; # equal
220 }
221
222 #
223 #  Use the original equation for these (but not for perl(), since we
224 #  need that to use perlFunction).
225 #
226 sub string {shift->{F}->string}
227 sub TeX    {shift->{F}->TeX}
228
229 sub cmp_class {'an Implicit Equation'}
230 sub showClass {shift->cmp_class}
231
232 #
233 #  Locate points that satisfy the equation
234 #
235 sub createPoints {
236   my \$self = shift;
237   my \$num_points = int(\$self->getFlag('ImplicitPoints',10));
238   \$num_points = 1 if \$num_points < 1;
239
240   #
241   #  Get some positive and negative test points (try up to 5 times)
242   #
243   my (\$OK,\$p,\$n,\$z,\$f,@zero); my \$k = 5;
244   while (!\$OK && --\$k) {(\$OK,\$p,\$n,\$z,\$f) = \$self->getPositiveNegativeZero(\$num_points)}
245   Value::Error("Can't find any solutions to your equation") unless \$OK;
246   my (\$P,@intervals) = \$self->getIntervals(\$p,\$n);
247
248   #
249   #  Get relative tolerance values and make them absolute
250   #
251   my \$minTolerance = \$self->getFlag('ImplicitAbsoluteMinTolerance');
252   my \$maxTolerance = \$self->getFlag('ImplicitAbsoluteMaxTolerance');
253   my \$tolerance = \$f * \$self->getFlag('ImplicitTolerance');
254   \$tolerance = \$minTolerance if \$tolerance < \$minTolerance;
255   \$tolerance = \$maxTolerance if \$tolerance > \$maxTolerance;
256   \$self->{tolerance} = \$tolerance unless \$self->{tolerance};
257   \$self->{bisect_tolerance} = \$self->{tolerance} * \$self->getFlag('BisectionTolerance')
258     unless \$self->{bisect_tolerance};
259   \$self->{point_tolerance} = \$P * \$self->getFlag('ImplicitPointTolerance')
260     unless \$self->{point_tolerance};
261
262   #
263   #  Locate solutions to be used for comparison test
264   #
265   @zero = @{\$z}; @zero = \$zero[0..\$num_points-1] if (scalar(@zero) > \$num_points);
266   for (\$i = 0; scalar(@zero) < \$num_points && \$i < scalar(@intervals); \$i++) {
267     my \$Q = \$self->Bisect(\$intervals[\$i][0],\$intervals[\$i][1]);
268     push(@zero,\$Q) if \$Q;
269   }
270   Value::Error("Can't find enough solutions for an effective test")
271     unless scalar(@zero) == \$num_points;
272
273   #
274   #  Save the solutions to the equation
275   #
276   \$self->{solutions} = [@zero];
277 }
278
279 #
280 #  Get random points and sort them by sign of the function.
281 #  Also, get the average function value, and indicate if
282 #  we actually did find both positive and negative values.
283 #
284 sub getPositiveNegativeZero {
285   my \$self = shift; my \$n = shift;
286   my (\$p,\$v) = \$self->SUPER::createRandomPoints(3*\$n);
287   my (@pos,@neg,@zero);
288   my \$f = 0; my \$k = 0;
289   foreach my \$i (0..scalar(@{\$v})-1) {
290     if (\$v->[\$i] == 0) {push(@zero,\$p->[\$i])} else {
291       \$f += abs(\$v->[\$i]); \$k++;
292       if (\$v->[\$i] > 0) {push(@pos,\$p->[\$i])} else {push(@neg,\$p->[\$i])}
293     }
294   }
295   \$f /= \$k if \$k;
296   return (scalar(@pos) && scalar(@neg),[@pos],[@neg],[@zero],\$f);
297 }
298
299 #
300 #  Make a list of positive and negative points sorted by
301 #  the distance between them.  Also return the average distance
302 #  between points.
303 #
304 sub getIntervals {
305   my \$self = shift; my \$pos = shift; my \$neg = shift;
306   my @intervals = (); my \$D = 0;
307   my \$point = \$self->Package("Point");
308   my \$context = \$self->context;
309   foreach my \$p (@{\$pos}) {
310     foreach my \$n (@{\$neg}) {
311       my \$d = abs(\$point->make(\$context,@{\$p}) - \$point->make(\$context,@{\$n}));
312       push(@intervals,[\$p,\$n,\$d]); \$D += \$d;
313     }
314   }
315   @intervals = main::PGsort(sub {\$_[0]->[2] < \$_[1]->[2]},@intervals);
316   return(\$D/scalar(@intervals),@intervals);
317 }
318
319 #
320 #  Use bisection algorithm to find a point where the function is zero
321 #  (i.e., where the original equation is an equality)
322 #  If we can't find a point (e.g., we are at a discontinuity),
323 #  return an undefined value.
324 #
325 sub Bisect {
326   my \$self = shift;
327   my \$tolerance  = \$self->getFlag('bisect_tolerance',1E-5);
328   my \$ptolerance = \$self->getFlag('point_tolerance',1E-9);
329   my \$m = \$self->getFlag('BisectionCutoff',30); my (\$P,\$f);
330   my \$point = \$self->Package("Point"); my \$context = \$self->context;
331   my \$P0 = \$point->make(\$context,@{\$_[0]}); my \$P1 = \$point->make(\$context,@{\$_[1]});
332   my (\$f0,\$f1) = @{\$self->createPointValues([\$P0->data,\$P1->data],1)};
333   for (my \$i = 0; \$i < \$m; \$i++) {
334     \$P = (\$P0+\$P1)/2; \$f = \$self->createPointValues([\$P->data]);
335     return unless ref(\$f);
336     \$f = \$f->[0];
337     return [\$P->value] if abs(\$f) < \$tolerance && abs(\$P1-\$P0) < \$ptolerance;
338     if (\$f > 0) {\$P0 = \$P; \$f0 = \$f} else {\$P1 = \$P; \$f1 = \$f}
339   }
340   return;
341 }
342
343 1;
```