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

Tue Aug 28 22:40:15 2007 UTC (12 years, 5 months ago) by dpvc
File size: 13530 byte(s)
```Add context names for the context(s) created here.
```

```    1 loadMacros("MathObjects.pl");
2
3 sub _parserImplicitEquation_init {ImplicitEquation::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 #  Create the ImplicitEquation package
132 #
133 package ImplicitEquation;
134 our @ISA = qw(Value::Formula);
135
136 sub Init {
137   my \$context = \$main::context{ImplicitEquation} = Parser::Context->getCopy("Numeric");
138   \$context->{name} = "ImplicitEquation";
139   \$context->variables->are(x=>'Real',y=>'Real');
140   \$context{precedence}{ImplicitEquation} = \$context->{precedence}{special};
141   Parser::BOP::equality->Allow(\$context);
142   \$context->flags->set(
143     ImplicitPoints => 10,
144     ImplicitTolerance => 1E-6,
145     ImplicitAbsoluteMinTolerance => 1E-3,
146     ImplicitAbsoluteMaxTolerance => 1,
147     ImplicitPointTolerance => 1E-9,
148     BisectionTolerance => .01,
149     BisectionCutoff => 40,
150   );
151
152   main::Context("ImplicitEquation");  ### FIXEME:  probably should require author to set this explicitly
153
154   main::PG_restricted_eval('sub ImplicitEquation {ImplicitEquation->new(@_)}');
155 }
156
157 sub new {
158   my \$self = shift; my \$class = ref(\$self) || \$self;
159   my \$context = (Value::isContext(\$_) ? shift : \$self->context);
160   my \$f = shift; return \$f if ref(\$f) eq \$class;
161   \$f = main::Formula(\$f);
162   Value::Error("Your formula doesn't look like an implicit equation")
163     unless \$f->type eq 'Equality';
164   my \$F = (\$context->Package("Formula")->new(\$context,\$f->{tree}{lop}) -
165      \$context->Package("Formula")->new(\$context,\$f->{tree}{rop}))->reduce;
166   \$F = \$context->Package("Formula")->new(\$F) unless Value::isFormula(\$F);
167   Value::Error("Your equation must be real-valued") unless \$F->isRealNumber;
168   Value::Error("Your equation should not be constant") if \$F->isConstant;
170     if (\$F->usesOneOf(\$context->variables->parameters));
171   \$F = bless \$F, \$class;
172   my %options = (@_);  # user can supply limits, tolerance, etc.
173   foreach my \$id (keys %options) {\$F->{\$id} = \$options{\$id}}
174   \$F->{F} = \$f; \$F->{isValue} = \$F->{isFormula} = 1;
175   \$F->createPoints unless \$F->{solutions};
176   return \$F;
177 }
178
179 #
180 #  Override the comparison method.
181 #
182 #  Turn the right-hand equation into an ImplicitEquation.  This creates
183 #  the test points (i.e., finds the solution points).  Then check
184 #  the professor's function on the student's test points and the
185 #  student's function on the professor's test points.
186 #
187
188 sub compare {
189   my (\$l,\$r) = @_; my \$self = \$l; my \$tolerance;
190   my @params; @params = (limits=>\$l->{limits}) if \$l->{limits};
191   \$r = ImplicitEquation->new(\$r,@params);
192   Value::Error("Functions from different contexts can't be compared")
193     unless \$l->{context} == \$r->{context};
194
195   #
196   #  They are not equal if couldn't get solutions for one of them
197   #
198   return 1 unless \$l->{solutions} && \$r->{solutions};
199
200   #
201   #  Test the right-hand function on the solutions of the left-hand one
202   #  and vice-versa
203   #
204   my \$rzeros = \$r->createPointValues(\$l->{solutions});
205   my \$lzeros = \$l->createPointValues(\$r->{solutions});
206   return 1 unless \$lzeros && \$rzeros;
207
208   #
209   #  Check that the values are, in fact, zeros
210   #
211   \$tolerance = \$r->getFlag('tolerance',1E-3);
212   foreach my \$v (@{\$rzeros}) {return 1 unless abs(\$v) < \$tolerance}
213   \$tolerance = \$l->getFlag('tolerance',1E-3);
214   foreach my \$v (@{\$lzeros}) {return 1 unless abs(\$v) < \$tolerance}
215
216   return 0; # equal
217 }
218
219 #
220 #  Use the original equation for these (but not for perl(), since we
221 #  need that to use perlFunction).
222 #
223 sub string {shift->{F}->string}
224 sub TeX    {shift->{F}->TeX}
225
226 sub cmp_class {'an Implicit Equation'}
227 sub showClass {shift->cmp_class}
228
229 #
230 #  Locate points that satisfy the equation
231 #
232 sub createPoints {
233   my \$self = shift;
234   my \$num_points = int(\$self->getFlag('ImplicitPoints',10));
235   \$num_points = 1 if \$num_points < 1;
236
237   #
238   #  Get some positive and negative test points (try up to 5 times)
239   #
240   my (\$OK,\$p,\$n,\$z,\$f,@zero); my \$k = 5;
241   while (!\$OK && --\$k) {(\$OK,\$p,\$n,\$z,\$f) = \$self->getPositiveNegativeZero(\$num_points)}
242   Value::Error("Can't find any solutions to your equation") unless \$OK;
243   my (\$P,@intervals) = \$self->getIntervals(\$p,\$n);
244
245   #
246   #  Get relative tolerance values and make them absolute
247   #
248   my \$minTolerance = \$self->getFlag('ImplicitAbsoluteMinTolerance');
249   my \$maxTolerance = \$self->getFlag('ImplicitAbsoluteMaxTolerance');
250   my \$tolerance = \$f * \$self->getFlag('ImplicitTolerance');
251   \$tolerance = \$minTolerance if \$tolerance < \$minTolerance;
252   \$tolerance = \$maxTolerance if \$tolerance > \$maxTolerance;
253   \$self->{tolerance} = \$tolerance unless \$self->{tolerance};
254   \$self->{bisect_tolerance} = \$self->{tolerance} * \$self->getFlag('BisectionTolerance')
255     unless \$self->{bisect_tolerance};
256   \$self->{point_tolerance} = \$P * \$self->getFlag('ImplicitPointTolerance')
257     unless \$self->{point_tolerance};
258
259   #
260   #  Locate solutions to be used for comparison test
261   #
262   @zero = @{\$z}; @zero = \$zero[0..\$num_points-1] if (scalar(@zero) > \$num_points);
263   for (\$i = 0; scalar(@zero) < \$num_points && \$i < scalar(@intervals); \$i++) {
264     my \$Q = \$self->Bisect(\$intervals[\$i],\$intervals[\$i]);
265     push(@zero,\$Q) if \$Q;
266   }
267   Value::Error("Can't find enough solutions for an effective test")
268     unless scalar(@zero) == \$num_points;
269
270   #
271   #  Save the solutions to the equation
272   #
273   \$self->{solutions} = [@zero];
274 }
275
276 #
277 #  Get random points and sort them by sign of the function.
278 #  Also, get the average function value, and indicate if
279 #  we actually did find both positive and negative values.
280 #
281 sub getPositiveNegativeZero {
282   my \$self = shift; my \$n = shift;
283   my (\$p,\$v) = \$self->SUPER::createRandomPoints(3*\$n);
284   my (@pos,@neg,@zero);
285   my \$f = 0; my \$k = 0;
286   foreach my \$i (0..scalar(@{\$v})-1) {
287     if (\$v->[\$i] == 0) {push(@zero,\$p->[\$i])} else {
288       \$f += abs(\$v->[\$i]); \$k++;
289       if (\$v->[\$i] > 0) {push(@pos,\$p->[\$i])} else {push(@neg,\$p->[\$i])}
290     }
291   }
292   \$f /= \$k if \$k;
293   return (scalar(@pos) && scalar(@neg),[@pos],[@neg],[@zero],\$f);
294 }
295
296 #
297 #  Make a list of positive and negative points sorted by
298 #  the distance between them.  Also return the average distance
299 #  between points.
300 #
301 sub getIntervals {
302   my \$self = shift; my \$pos = shift; my \$neg = shift;
303   my @intervals = (); my \$D = 0;
304   my \$point = \$self->Package("Point");
305   my \$context = \$self->context;
306   foreach my \$p (@{\$pos}) {
307     foreach my \$n (@{\$neg}) {
308       my \$d = abs(\$point->make(\$context,@{\$p}) - \$point->make(\$context,@{\$n}));
309       push(@intervals,[\$p,\$n,\$d]); \$D += \$d;
310     }
311   }
312   @intervals = main::PGsort(sub {\$_-> < \$_->},@intervals);
313   return(\$D/scalar(@intervals),@intervals);
314 }
315
316 #
317 #  Use bisection algorithm to find a point where the function is zero
318 #  (i.e., where the original equation is an equality)
319 #  If we can't find a point (e.g., we are at a discontinuity),
320 #  return an undefined value.
321 #
322 sub Bisect {
323   my \$self = shift;
324   my \$tolerance  = \$self->getFlag('bisect_tolerance',1E-5);
325   my \$ptolerance = \$self->getFlag('point_tolerance',1E-9);
326   my \$m = \$self->getFlag('BisectionCutoff',30); my (\$P,\$f);
327   my \$point = \$self->Package("Point"); my \$context = \$self->context;
328   my \$P0 = \$point->make(\$context,@{\$_}); my \$P1 = \$point->make(\$context,@{\$_});
329   my (\$f0,\$f1) = @{\$self->createPointValues([\$P0->data,\$P1->data],1)};
330   for (my \$i = 0; \$i < \$m; \$i++) {
331     \$P = (\$P0+\$P1)/2; \$f = \$self->createPointValues([\$P->data]);
332     return unless ref(\$f);
333     \$f = \$f->;
334     return [\$P->value] if abs(\$f) < \$tolerance && abs(\$P1-\$P0) < \$ptolerance;
335     if (\$f > 0) {\$P0 = \$P; \$f0 = \$f} else {\$P1 = \$P; \$f1 = \$f}
336   }
337   return;
338 }
339
340 1;
```