View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.optim.nonlinear.scalar.gradient;
24  
25  import org.hipparchus.analysis.MultivariateFunction;
26  import org.hipparchus.analysis.MultivariateVectorFunction;
27  import org.hipparchus.exception.MathRuntimeException;
28  import org.hipparchus.geometry.euclidean.twod.Vector2D;
29  import org.hipparchus.linear.BlockRealMatrix;
30  import org.hipparchus.linear.RealMatrix;
31  import org.hipparchus.optim.InitialGuess;
32  import org.hipparchus.optim.MaxEval;
33  import org.hipparchus.optim.PointValuePair;
34  import org.hipparchus.optim.SimpleBounds;
35  import org.hipparchus.optim.SimpleValueChecker;
36  import org.hipparchus.optim.nonlinear.scalar.GoalType;
37  import org.hipparchus.optim.nonlinear.scalar.ObjectiveFunction;
38  import org.hipparchus.optim.nonlinear.scalar.ObjectiveFunctionGradient;
39  import org.junit.Assert;
40  import org.junit.Test;
41  
42  /**
43   * <p>Some of the unit tests are re-implementations of the MINPACK <a
44   * href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
45   * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files.
46   * The redistribution policy for MINPACK is available <a
47   * href="http://www.netlib.org/minpack/disclaimer">here</a>, for
48   * convenience, it is reproduced below.</p>
49   *
50   * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
51   * <tr><td>
52   *    Minpack Copyright Notice (1999) University of Chicago.
53   *    All rights reserved
54   * </td></tr>
55   * <tr><td>
56   * Redistribution and use in source and binary forms, with or without
57   * modification, are permitted provided that the following conditions
58   * are met:
59   * <ol>
60   *  <li>Redistributions of source code must retain the above copyright
61   *      notice, this list of conditions and the following disclaimer.</li>
62   * <li>Redistributions in binary form must reproduce the above
63   *     copyright notice, this list of conditions and the following
64   *     disclaimer in the documentation and/or other materials provided
65   *     with the distribution.</li>
66   * <li>The end-user documentation included with the redistribution, if any,
67   *     must include the following acknowledgment:
68   *     <code>This product includes software developed by the University of
69   *           Chicago, as Operator of Argonne National Laboratory.</code>
70   *     Alternately, this acknowledgment may appear in the software itself,
71   *     if and wherever such third-party acknowledgments normally appear.</li>
72   * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
73   *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
74   *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
75   *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
76   *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
77   *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
78   *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
79   *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
80   *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
81   *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
82   *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
83   *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
84   *     BE CORRECTED.</strong></li>
85   * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
86   *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
87   *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
88   *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
89   *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
90   *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
91   *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
92   *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
93   *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
94   *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
95   * <ol></td></tr>
96   * </table>
97   *
98   * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
99   * @author Burton S. Garbow (original fortran minpack tests)
100  * @author Kenneth E. Hillstrom (original fortran minpack tests)
101  * @author Jorge J. More (original fortran minpack tests)
102  * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
103  */
104 public class NonLinearConjugateGradientOptimizerTest {
105     @Test(expected=MathRuntimeException.class)
106     public void testBoundsUnsupported() {
107         LinearProblem problem
108             = new LinearProblem(new double[][] { { 2 } }, new double[] { 3 });
109         NonLinearConjugateGradientOptimizer optimizer
110             = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
111                                                       new SimpleValueChecker(1e-6, 1e-6),
112                                                       1e-3, 1e-3, 1);
113         optimizer.optimize(new MaxEval(100),
114                            problem.getObjectiveFunction(),
115                            problem.getObjectiveFunctionGradient(),
116                            GoalType.MINIMIZE,
117                            new InitialGuess(new double[] { 0 }),
118                            new SimpleBounds(new double[] { -1 },
119                                             new double[] { 1 }));
120     }
121 
122     @Test
123     public void testTrivial() {
124         LinearProblem problem
125             = new LinearProblem(new double[][] { { 2 } }, new double[] { 3 });
126         NonLinearConjugateGradientOptimizer optimizer
127             = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
128                                                       new SimpleValueChecker(1e-6, 1e-6),
129                                                       1e-3, 1e-3, 1);
130         PointValuePair optimum
131             = optimizer.optimize(new MaxEval(100),
132                                  problem.getObjectiveFunction(),
133                                  problem.getObjectiveFunctionGradient(),
134                                  GoalType.MINIMIZE,
135                                  new InitialGuess(new double[] { 0 }));
136         Assert.assertEquals(1.5, optimum.getPoint()[0], 1.0e-10);
137         Assert.assertEquals(0.0, optimum.getValue(), 1.0e-10);
138 
139         // Check that the number of iterations is updated (MATH-949).
140         Assert.assertTrue(optimizer.getIterations() > 0);
141     }
142 
143     @Test
144     public void testColumnsPermutation() {
145         LinearProblem problem
146             = new LinearProblem(new double[][] { { 1.0, -1.0 }, { 0.0, 2.0 }, { 1.0, -2.0 } },
147                                 new double[] { 4.0, 6.0, 1.0 });
148 
149         NonLinearConjugateGradientOptimizer optimizer
150             = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
151                                                       new SimpleValueChecker(1e-6, 1e-6),
152                                                       1e-3, 1e-3, 1);
153         PointValuePair optimum
154             = optimizer.optimize(new MaxEval(100),
155                                  problem.getObjectiveFunction(),
156                                  problem.getObjectiveFunctionGradient(),
157                                  GoalType.MINIMIZE,
158                                  new InitialGuess(new double[] { 0, 0 }));
159         Assert.assertEquals(7.0, optimum.getPoint()[0], 1.0e-10);
160         Assert.assertEquals(3.0, optimum.getPoint()[1], 1.0e-10);
161         Assert.assertEquals(0.0, optimum.getValue(), 1.0e-10);
162 
163     }
164 
165     @Test
166     public void testNoDependency() {
167         LinearProblem problem = new LinearProblem(new double[][] {
168                 { 2, 0, 0, 0, 0, 0 },
169                 { 0, 2, 0, 0, 0, 0 },
170                 { 0, 0, 2, 0, 0, 0 },
171                 { 0, 0, 0, 2, 0, 0 },
172                 { 0, 0, 0, 0, 2, 0 },
173                 { 0, 0, 0, 0, 0, 2 }
174         }, new double[] { 0.0, 1.1, 2.2, 3.3, 4.4, 5.5 });
175         NonLinearConjugateGradientOptimizer optimizer
176             = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
177                                                       new SimpleValueChecker(1e-6, 1e-6),
178                                                       1e-3, 1e-3, 1);
179         PointValuePair optimum
180             = optimizer.optimize(new MaxEval(100),
181                                  problem.getObjectiveFunction(),
182                                  problem.getObjectiveFunctionGradient(),
183                                  GoalType.MINIMIZE,
184                                  new InitialGuess(new double[] { 0, 0, 0, 0, 0, 0 }));
185         for (int i = 0; i < problem.target.length; ++i) {
186             Assert.assertEquals(0.55 * i, optimum.getPoint()[i], 1.0e-10);
187         }
188     }
189 
190     @Test
191     public void testOneSet() {
192         LinearProblem problem = new LinearProblem(new double[][] {
193                 {  1,  0, 0 },
194                 { -1,  1, 0 },
195                 {  0, -1, 1 }
196         }, new double[] { 1, 1, 1});
197         NonLinearConjugateGradientOptimizer optimizer
198             = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
199                                                       new SimpleValueChecker(1e-6, 1e-6),
200                                                       1e-3, 1e-3, 1);
201         PointValuePair optimum
202             = optimizer.optimize(new MaxEval(100),
203                                  problem.getObjectiveFunction(),
204                                  problem.getObjectiveFunctionGradient(),
205                                  GoalType.MINIMIZE,
206                                  new InitialGuess(new double[] { 0, 0, 0 }));
207         Assert.assertEquals(1.0, optimum.getPoint()[0], 1.0e-10);
208         Assert.assertEquals(2.0, optimum.getPoint()[1], 1.0e-10);
209         Assert.assertEquals(3.0, optimum.getPoint()[2], 1.0e-10);
210 
211     }
212 
213     @Test
214     public void testTwoSets() {
215         final double epsilon = 1.0e-7;
216         LinearProblem problem = new LinearProblem(new double[][] {
217                 {  2,  1,   0,  4,       0, 0 },
218                 { -4, -2,   3, -7,       0, 0 },
219                 {  4,  1,  -2,  8,       0, 0 },
220                 {  0, -3, -12, -1,       0, 0 },
221                 {  0,  0,   0,  0, epsilon, 1 },
222                 {  0,  0,   0,  0,       1, 1 }
223         }, new double[] { 2, -9, 2, 2, 1 + epsilon * epsilon, 2});
224 
225         final Preconditioner preconditioner
226             = new Preconditioner() {
227                     public double[] precondition(double[] point, double[] r) {
228                         double[] d = r.clone();
229                         d[0] /=  72.0;
230                         d[1] /=  30.0;
231                         d[2] /= 314.0;
232                         d[3] /= 260.0;
233                         d[4] /= 2 * (1 + epsilon * epsilon);
234                         d[5] /= 4.0;
235                         return d;
236                     }
237                 };
238 
239         NonLinearConjugateGradientOptimizer optimizer
240            = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
241                                                      new SimpleValueChecker(1e-13, 1e-13),
242                                                      1e-7, 1e-7, 1,
243                                                      preconditioner);
244 
245         PointValuePair optimum
246             = optimizer.optimize(new MaxEval(100),
247                                  problem.getObjectiveFunction(),
248                                  problem.getObjectiveFunctionGradient(),
249                                  GoalType.MINIMIZE,
250                                  new InitialGuess(new double[] { 0, 0, 0, 0, 0, 0 }));
251 
252         final double[] result = optimum.getPoint();
253         final double[] expected = {3, 4, -1, -2, 1 + epsilon, 1 - epsilon};
254 
255         Assert.assertEquals(expected[0], result[0], 1.0e-7);
256         Assert.assertEquals(expected[1], result[1], 1.0e-7);
257         Assert.assertEquals(expected[2], result[2], 1.0e-9);
258         Assert.assertEquals(expected[3], result[3], 1.0e-8);
259         Assert.assertEquals(expected[4] + epsilon, result[4], 1.0e-6);
260         Assert.assertEquals(expected[5] - epsilon, result[5], 1.0e-6);
261 
262     }
263 
264     @Test
265     public void testNonInversible() {
266         LinearProblem problem = new LinearProblem(new double[][] {
267                 {  1, 2, -3 },
268                 {  2, 1,  3 },
269                 { -3, 0, -9 }
270         }, new double[] { 1, 1, 1 });
271         NonLinearConjugateGradientOptimizer optimizer
272             = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
273                                                       new SimpleValueChecker(1e-6, 1e-6),
274                                                       1e-3, 1e-3, 1);
275         PointValuePair optimum
276             = optimizer.optimize(new MaxEval(100),
277                                  problem.getObjectiveFunction(),
278                                  problem.getObjectiveFunctionGradient(),
279                                  GoalType.MINIMIZE,
280                                  new InitialGuess(new double[] { 0, 0, 0 }));
281         Assert.assertTrue(optimum.getValue() > 0.5);
282     }
283 
284     @Test
285     public void testIllConditioned() {
286         LinearProblem problem1 = new LinearProblem(new double[][] {
287                 { 10.0, 7.0,  8.0,  7.0 },
288                 {  7.0, 5.0,  6.0,  5.0 },
289                 {  8.0, 6.0, 10.0,  9.0 },
290                 {  7.0, 5.0,  9.0, 10.0 }
291         }, new double[] { 32, 23, 33, 31 });
292         NonLinearConjugateGradientOptimizer optimizer
293             = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
294                                                       new SimpleValueChecker(1e-13, 1e-13),
295                                                       1e-15, 1e-15, 1);
296         PointValuePair optimum1
297             = optimizer.optimize(new MaxEval(200),
298                                  problem1.getObjectiveFunction(),
299                                  problem1.getObjectiveFunctionGradient(),
300                                  GoalType.MINIMIZE,
301                                  new InitialGuess(new double[] { 0, 1, 2, 3 }));
302         Assert.assertEquals(1.0, optimum1.getPoint()[0], 1.0e-4);
303         Assert.assertEquals(1.0, optimum1.getPoint()[1], 1.0e-3);
304         Assert.assertEquals(1.0, optimum1.getPoint()[2], 1.0e-4);
305         Assert.assertEquals(1.0, optimum1.getPoint()[3], 1.0e-4);
306 
307         LinearProblem problem2 = new LinearProblem(new double[][] {
308                 { 10.00, 7.00, 8.10, 7.20 },
309                 {  7.08, 5.04, 6.00, 5.00 },
310                 {  8.00, 5.98, 9.89, 9.00 },
311                 {  6.99, 4.99, 9.00, 9.98 }
312         }, new double[] { 32, 23, 33, 31 });
313         PointValuePair optimum2
314             = optimizer.optimize(new MaxEval(200),
315                                  problem2.getObjectiveFunction(),
316                                  problem2.getObjectiveFunctionGradient(),
317                                  GoalType.MINIMIZE,
318                                  new InitialGuess(new double[] { 0, 1, 2, 3 }));
319 
320         final double[] result2 = optimum2.getPoint();
321         final double[] expected2 = {-81, 137, -34, 22};
322 
323         Assert.assertEquals(expected2[0], result2[0], 2);
324         Assert.assertEquals(expected2[1], result2[1], 4);
325         Assert.assertEquals(expected2[2], result2[2], 1);
326         Assert.assertEquals(expected2[3], result2[3], 1);
327     }
328 
329     @Test
330     public void testMoreEstimatedParametersSimple() {
331         LinearProblem problem = new LinearProblem(new double[][] {
332                 { 3.0, 2.0,  0.0, 0.0 },
333                 { 0.0, 1.0, -1.0, 1.0 },
334                 { 2.0, 0.0,  1.0, 0.0 }
335         }, new double[] { 7.0, 3.0, 5.0 });
336 
337         NonLinearConjugateGradientOptimizer optimizer
338             = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
339                                                       new SimpleValueChecker(1e-6, 1e-6),
340                                                       1e-3, 1e-3, 1);
341         PointValuePair optimum
342             = optimizer.optimize(new MaxEval(100),
343                                  problem.getObjectiveFunction(),
344                                  problem.getObjectiveFunctionGradient(),
345                                  GoalType.MINIMIZE,
346                                  new InitialGuess(new double[] { 7, 6, 5, 4 }));
347         Assert.assertEquals(0, optimum.getValue(), 1.0e-10);
348 
349     }
350 
351     @Test
352     public void testMoreEstimatedParametersUnsorted() {
353         LinearProblem problem = new LinearProblem(new double[][] {
354                  { 1.0, 1.0,  0.0,  0.0, 0.0,  0.0 },
355                  { 0.0, 0.0,  1.0,  1.0, 1.0,  0.0 },
356                  { 0.0, 0.0,  0.0,  0.0, 1.0, -1.0 },
357                  { 0.0, 0.0, -1.0,  1.0, 0.0,  1.0 },
358                  { 0.0, 0.0,  0.0, -1.0, 1.0,  0.0 }
359         }, new double[] { 3.0, 12.0, -1.0, 7.0, 1.0 });
360         NonLinearConjugateGradientOptimizer optimizer
361            = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
362                                                      new SimpleValueChecker(1e-6, 1e-6),
363                                                      1e-3, 1e-3, 1);
364         PointValuePair optimum
365             = optimizer.optimize(new MaxEval(100),
366                                  problem.getObjectiveFunction(),
367                                  problem.getObjectiveFunctionGradient(),
368                                  GoalType.MINIMIZE,
369                                  new InitialGuess(new double[] { 2, 2, 2, 2, 2, 2 }));
370         Assert.assertEquals(0, optimum.getValue(), 1.0e-10);
371     }
372 
373     @Test
374     public void testRedundantEquations() {
375         LinearProblem problem = new LinearProblem(new double[][] {
376                 { 1.0,  1.0 },
377                 { 1.0, -1.0 },
378                 { 1.0,  3.0 }
379         }, new double[] { 3.0, 1.0, 5.0 });
380 
381         NonLinearConjugateGradientOptimizer optimizer
382             = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
383                                                       new SimpleValueChecker(1e-6, 1e-6),
384                                                       1e-3, 1e-3, 1);
385         PointValuePair optimum
386             = optimizer.optimize(new MaxEval(100),
387                                  problem.getObjectiveFunction(),
388                                  problem.getObjectiveFunctionGradient(),
389                                  GoalType.MINIMIZE,
390                                  new InitialGuess(new double[] { 1, 1 }));
391         Assert.assertEquals(2.0, optimum.getPoint()[0], 1.0e-8);
392         Assert.assertEquals(1.0, optimum.getPoint()[1], 1.0e-8);
393 
394     }
395 
396     @Test
397     public void testInconsistentEquations() {
398         LinearProblem problem = new LinearProblem(new double[][] {
399                 { 1.0,  1.0 },
400                 { 1.0, -1.0 },
401                 { 1.0,  3.0 }
402         }, new double[] { 3.0, 1.0, 4.0 });
403 
404         NonLinearConjugateGradientOptimizer optimizer
405             = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
406                                                       new SimpleValueChecker(1e-6, 1e-6),
407                                                       1e-3, 1e-3, 1);
408         PointValuePair optimum
409             = optimizer.optimize(new MaxEval(100),
410                                  problem.getObjectiveFunction(),
411                                  problem.getObjectiveFunctionGradient(),
412                                  GoalType.MINIMIZE,
413                                  new InitialGuess(new double[] { 1, 1 }));
414         Assert.assertTrue(optimum.getValue() > 0.1);
415 
416     }
417 
418     @Test
419     public void testCircleFitting() {
420         CircleScalar problem = new CircleScalar();
421         problem.addPoint( 30.0,  68.0);
422         problem.addPoint( 50.0,  -6.0);
423         problem.addPoint(110.0, -20.0);
424         problem.addPoint( 35.0,  15.0);
425         problem.addPoint( 45.0,  97.0);
426         NonLinearConjugateGradientOptimizer optimizer
427            = new NonLinearConjugateGradientOptimizer(NonLinearConjugateGradientOptimizer.Formula.POLAK_RIBIERE,
428                                                      new SimpleValueChecker(1e-30, 1e-30),
429                                                      1e-15, 1e-13, 1);
430         PointValuePair optimum
431             = optimizer.optimize(new MaxEval(100),
432                                  problem.getObjectiveFunction(),
433                                  problem.getObjectiveFunctionGradient(),
434                                  GoalType.MINIMIZE,
435                                  new InitialGuess(new double[] { 98.680, 47.345 }));
436         Vector2D center = new Vector2D(optimum.getPointRef()[0], optimum.getPointRef()[1]);
437         Assert.assertEquals(69.960161753, problem.getRadius(center), 1.0e-8);
438         Assert.assertEquals(96.075902096, center.getX(), 1.0e-7);
439         Assert.assertEquals(48.135167894, center.getY(), 1.0e-6);
440     }
441 
442     private static class LinearProblem {
443         final RealMatrix factors;
444         final double[] target;
445 
446         public LinearProblem(double[][] factors,
447                              double[] target) {
448             this.factors = new BlockRealMatrix(factors);
449             this.target  = target;
450         }
451 
452         public ObjectiveFunction getObjectiveFunction() {
453             return new ObjectiveFunction(new MultivariateFunction() {
454                     public double value(double[] point) {
455                         double[] y = factors.operate(point);
456                         double sum = 0;
457                         for (int i = 0; i < y.length; ++i) {
458                             double ri = y[i] - target[i];
459                             sum += ri * ri;
460                         }
461                         return sum;
462                     }
463                 });
464         }
465 
466         public ObjectiveFunctionGradient getObjectiveFunctionGradient() {
467             return new ObjectiveFunctionGradient(new MultivariateVectorFunction() {
468                     public double[] value(double[] point) {
469                         double[] r = factors.operate(point);
470                         for (int i = 0; i < r.length; ++i) {
471                             r[i] -= target[i];
472                         }
473                         double[] p = factors.transpose().operate(r);
474                         for (int i = 0; i < p.length; ++i) {
475                             p[i] *= 2;
476                         }
477                         return p;
478                     }
479                 });
480         }
481     }
482 }