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

Sun Mar 20 13:13:44 2005 UTC (14 years, 11 months ago) by dpvc
File size: 13040 byte(s)
```Implements an experimental answer checker for implicitly defined
curves and surfaces.
```

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