Source code

001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math3.analysis.differentiation;
018
019import java.io.Serializable;
020import java.util.Collections;
021import java.util.HashMap;
022import java.util.Map;
023
024import org.apache.commons.math3.Field;
025import org.apache.commons.math3.FieldElement;
026import org.apache.commons.math3.RealFieldElement;
027import org.apache.commons.math3.exception.DimensionMismatchException;
028import org.apache.commons.math3.util.FastMath;
029import org.apache.commons.math3.util.MathArrays;
030import org.apache.commons.math3.util.MathUtils;
031import org.apache.commons.math3.util.Precision;
032
033/**
034 * First derivative computation with large number of variables.
035 * <p>
036 * This class plays a similar role to {@link DerivativeStructure}, with
037 * a focus on efficiency when dealing with large number of independent variables
038 * and most computation depend only on a few of them, and when only first derivative
039 * is desired. When these conditions are met, this class should be much faster than
040 * {@link DerivativeStructure} and use less memory.
041 * </p>
042 *
043 * @since 3.3
044 */
045public class SparseGradient implements RealFieldElement<SparseGradient>, Serializable {
046
047    /** Serializable UID. */
048    private static final long serialVersionUID = 20131025L;
049
050    /** Value of the calculation. */
051    private double value;
052
053    /** Stored derivative, each key representing a different independent variable. */
054    private final Map<Integer, Double> derivatives;
055
056    /** Internal constructor.
057     * @param value value of the function
058     * @param derivatives derivatives map, a deep copy will be performed,
059     * so the map given here will remain safe from changes in the new instance,
060     * may be null to create an empty derivatives map, i.e. a constant value
061     */
062    private SparseGradient(final double value, final Map<Integer, Double> derivatives) {
063        this.value = value;
064        this.derivatives = new HashMap<Integer, Double>();
065        if (derivatives != null) {
066            this.derivatives.putAll(derivatives);
067        }
068    }
069
070    /** Internal constructor.
071     * @param value value of the function
072     * @param scale scaling factor to apply to all derivatives
073     * @param derivatives derivatives map, a deep copy will be performed,
074     * so the map given here will remain safe from changes in the new instance,
075     * may be null to create an empty derivatives map, i.e. a constant value
076     */
077    private SparseGradient(final double value, final double scale,
078                             final Map<Integer, Double> derivatives) {
079        this.value = value;
080        this.derivatives = new HashMap<Integer, Double>();
081        if (derivatives != null) {
082            for (final Map.Entry<Integer, Double> entry : derivatives.entrySet()) {
083                this.derivatives.put(entry.getKey(), scale * entry.getValue());
084            }
085        }
086    }
087
088    /** Factory method creating a constant.
089     * @param value value of the constant
090     * @return a new instance
091     */
092    public static SparseGradient createConstant(final double value) {
093        return new SparseGradient(value, Collections.<Integer, Double> emptyMap());
094    }
095
096    /** Factory method creating an independent variable.
097     * @param idx index of the variable
098     * @param value value of the variable
099     * @return a new instance
100     */
101    public static SparseGradient createVariable(final int idx, final double value) {
102        return new SparseGradient(value, Collections.singletonMap(idx, 1.0));
103    }
104
105    /**
106     * Find the number of variables.
107     * @return number of variables
108     */
109    public int numVars() {
110        return derivatives.size();
111    }
112
113    /**
114     * Get the derivative with respect to a particular index variable.
115     *
116     * @param index index to differentiate with.
117     * @return derivative with respect to a particular index variable
118     */
119    public double getDerivative(final int index) {
120        final Double out = derivatives.get(index);
121        return (out == null) ? 0.0 : out;
122    }
123
124    /**
125     * Get the value of the function.
126     * @return value of the function.
127     */
128    public double getValue() {
129        return value;
130    }
131
132    /** {@inheritDoc} */
133    public double getReal() {
134        return value;
135    }
136
137    /** {@inheritDoc} */
138    public SparseGradient add(final SparseGradient a) {
139        final SparseGradient out = new SparseGradient(value + a.value, derivatives);
140        for (Map.Entry<Integer, Double> entry : a.derivatives.entrySet()) {
141            final int id = entry.getKey();
142            final Double old = out.derivatives.get(id);
143            if (old == null) {
144                out.derivatives.put(id, entry.getValue());
145            } else {
146                out.derivatives.put(id, old + entry.getValue());
147            }
148        }
149
150        return out;
151    }
152
153    /**
154     * Add in place.
155     * <p>
156     * This method is designed to be faster when used multiple times in a loop.
157     * </p>
158     * <p>
159     * The instance is changed here, in order to not change the
160     * instance the {@link #add(SparseGradient)} method should
161     * be used.
162     * </p>
163     * @param a instance to add
164     */
165    public void addInPlace(final SparseGradient a) {
166        value += a.value;
167        for (final Map.Entry<Integer, Double> entry : a.derivatives.entrySet()) {
168            final int id = entry.getKey();
169            final Double old = derivatives.get(id);
170            if (old == null) {
171                derivatives.put(id, entry.getValue());
172            } else {
173                derivatives.put(id, old + entry.getValue());
174            }
175        }
176    }
177
178    /** {@inheritDoc} */
179    public SparseGradient add(final double c) {
180        final SparseGradient out = new SparseGradient(value + c, derivatives);
181        return out;
182    }
183
184    /** {@inheritDoc} */
185    public SparseGradient subtract(final SparseGradient a) {
186        final SparseGradient out = new SparseGradient(value - a.value, derivatives);
187        for (Map.Entry<Integer, Double> entry : a.derivatives.entrySet()) {
188            final int id = entry.getKey();
189            final Double old = out.derivatives.get(id);
190            if (old == null) {
191                out.derivatives.put(id, -entry.getValue());
192            } else {
193                out.derivatives.put(id, old - entry.getValue());
194            }
195        }
196        return out;
197    }
198
199    /** {@inheritDoc} */
200    public SparseGradient subtract(double c) {
201        return new SparseGradient(value - c, derivatives);
202    }
203
204    /** {@inheritDoc} */
205    public SparseGradient multiply(final SparseGradient a) {
206        final SparseGradient out =
207            new SparseGradient(value * a.value, Collections.<Integer, Double> emptyMap());
208
209        // Derivatives.
210        for (Map.Entry<Integer, Double> entry : derivatives.entrySet()) {
211            out.derivatives.put(entry.getKey(), a.value * entry.getValue());
212        }
213        for (Map.Entry<Integer, Double> entry : a.derivatives.entrySet()) {
214            final int id = entry.getKey();
215            final Double old = out.derivatives.get(id);
216            if (old == null) {
217                out.derivatives.put(id, value * entry.getValue());
218            } else {
219                out.derivatives.put(id, old + value * entry.getValue());
220            }
221        }
222        return out;
223    }
224
225    /**
226     * Multiply in place.
227     * <p>
228     * This method is designed to be faster when used multiple times in a loop.
229     * </p>
230     * <p>
231     * The instance is changed here, in order to not change the
232     * instance the {@link #add(SparseGradient)} method should
233     * be used.
234     * </p>
235     * @param a instance to multiply
236     */
237    public void multiplyInPlace(final SparseGradient a) {
238        // Derivatives.
239        for (Map.Entry<Integer, Double> entry : derivatives.entrySet()) {
240            derivatives.put(entry.getKey(), a.value * entry.getValue());
241        }
242        for (Map.Entry<Integer, Double> entry : a.derivatives.entrySet()) {
243            final int id = entry.getKey();
244            final Double old = derivatives.get(id);
245            if (old == null) {
246                derivatives.put(id, value * entry.getValue());
247            } else {
248                derivatives.put(id, old + value * entry.getValue());
249            }
250        }
251        value *= a.value;
252    }
253
254    /** {@inheritDoc} */
255    public SparseGradient multiply(final double c) {
256        return new SparseGradient(value * c, c, derivatives);
257    }
258
259    /** {@inheritDoc} */
260    public SparseGradient multiply(final int n) {
261        return new SparseGradient(value * n, n, derivatives);
262    }
263
264    /** {@inheritDoc} */
265    public SparseGradient divide(final SparseGradient a) {
266        final SparseGradient out = new SparseGradient(value / a.value, Collections.<Integer, Double> emptyMap());
267
268        // Derivatives.
269        for (Map.Entry<Integer, Double> entry : derivatives.entrySet()) {
270            out.derivatives.put(entry.getKey(), entry.getValue() / a.value);
271        }
272        for (Map.Entry<Integer, Double> entry : a.derivatives.entrySet()) {
273            final int id = entry.getKey();
274            final Double old = out.derivatives.get(id);
275            if (old == null) {
276                out.derivatives.put(id, -out.value / a.value * entry.getValue());
277            } else {
278                out.derivatives.put(id, old - out.value / a.value * entry.getValue());
279            }
280        }
281        return out;
282    }
283
284    /** {@inheritDoc} */
285    public SparseGradient divide(final double c) {
286        return new SparseGradient(value / c, 1.0 / c, derivatives);
287    }
288
289    /** {@inheritDoc} */
290    public SparseGradient negate() {
291        return new SparseGradient(-value, -1.0, derivatives);
292    }
293
294    /** {@inheritDoc} */
295    public Field<SparseGradient> getField() {
296        return new Field<SparseGradient>() {
297
298            /** {@inheritDoc} */
299            public SparseGradient getZero() {
300                return createConstant(0);
301            }
302
303            /** {@inheritDoc} */
304            public SparseGradient getOne() {
305                return createConstant(1);
306            }
307
308            /** {@inheritDoc} */
309            public Class<? extends FieldElement<SparseGradient>> getRuntimeClass() {
310                return SparseGradient.class;
311            }
312
313        };
314    }
315
316    /** {@inheritDoc} */
317    public SparseGradient remainder(final double a) {
318        return new SparseGradient(FastMath.IEEEremainder(value, a), derivatives);
319    }
320
321    /** {@inheritDoc} */
322    public SparseGradient remainder(final SparseGradient a) {
323
324        // compute k such that lhs % rhs = lhs - k rhs
325        final double rem = FastMath.IEEEremainder(value, a.value);
326        final double k   = FastMath.rint((value - rem) / a.value);
327
328        return subtract(a.multiply(k));
329
330    }
331
332    /** {@inheritDoc} */
333    public SparseGradient abs() {
334        if (Double.doubleToLongBits(value) < 0) {
335            // we use the bits representation to also handle -0.0
336            return negate();
337        } else {
338            return this;
339        }
340    }
341
342    /** {@inheritDoc} */
343    public SparseGradient ceil() {
344        return createConstant(FastMath.ceil(value));
345    }
346
347    /** {@inheritDoc} */
348    public SparseGradient floor() {
349        return createConstant(FastMath.floor(value));
350    }
351
352    /** {@inheritDoc} */
353    public SparseGradient rint() {
354        return createConstant(FastMath.rint(value));
355    }
356
357    /** {@inheritDoc} */
358    public long round() {
359        return FastMath.round(value);
360    }
361
362    /** {@inheritDoc} */
363    public SparseGradient signum() {
364        return createConstant(FastMath.signum(value));
365    }
366
367    /** {@inheritDoc} */
368    public SparseGradient copySign(final SparseGradient sign) {
369        final long m = Double.doubleToLongBits(value);
370        final long s = Double.doubleToLongBits(sign.value);
371        if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
372            return this;
373        }
374        return negate(); // flip sign
375    }
376
377    /** {@inheritDoc} */
378    public SparseGradient copySign(final double sign) {
379        final long m = Double.doubleToLongBits(value);
380        final long s = Double.doubleToLongBits(sign);
381        if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
382            return this;
383        }
384        return negate(); // flip sign
385    }
386
387    /** {@inheritDoc} */
388    public SparseGradient scalb(final int n) {
389        final SparseGradient out = new SparseGradient(FastMath.scalb(value, n), Collections.<Integer, Double> emptyMap());
390        for (Map.Entry<Integer, Double> entry : derivatives.entrySet()) {
391            out.derivatives.put(entry.getKey(), FastMath.scalb(entry.getValue(), n));
392        }
393        return out;
394    }
395
396    /** {@inheritDoc} */
397    public SparseGradient hypot(final SparseGradient y) {
398        if (Double.isInfinite(value) || Double.isInfinite(y.value)) {
399            return createConstant(Double.POSITIVE_INFINITY);
400        } else if (Double.isNaN(value) || Double.isNaN(y.value)) {
401            return createConstant(Double.NaN);
402        } else {
403
404            final int expX = FastMath.getExponent(value);
405            final int expY = FastMath.getExponent(y.value);
406            if (expX > expY + 27) {
407                // y is negligible with respect to x
408                return abs();
409            } else if (expY > expX + 27) {
410                // x is negligible with respect to y
411                return y.abs();
412            } else {
413
414                // find an intermediate scale to avoid both overflow and underflow
415                final int middleExp = (expX + expY) / 2;
416
417                // scale parameters without losing precision
418                final SparseGradient scaledX = scalb(-middleExp);
419                final SparseGradient scaledY = y.scalb(-middleExp);
420
421                // compute scaled hypotenuse
422                final SparseGradient scaledH =
423                        scaledX.multiply(scaledX).add(scaledY.multiply(scaledY)).sqrt();
424
425                // remove scaling
426                return scaledH.scalb(middleExp);
427
428            }
429
430        }
431    }
432
433    /**
434     * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
435     * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
436     * avoiding intermediate overflow or underflow.
437     *
438     * <ul>
439     * <li> If either argument is infinite, then the result is positive infinity.</li>
440     * <li> else, if either argument is NaN then the result is NaN.</li>
441     * </ul>
442     *
443     * @param x a value
444     * @param y a value
445     * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
446     */
447    public static SparseGradient hypot(final SparseGradient x, final SparseGradient y) {
448        return x.hypot(y);
449    }
450
451    /** {@inheritDoc} */
452    public SparseGradient reciprocal() {
453        return new SparseGradient(1.0 / value, -1.0 / (value * value), derivatives);
454    }
455
456    /** {@inheritDoc} */
457    public SparseGradient sqrt() {
458        final double sqrt = FastMath.sqrt(value);
459        return new SparseGradient(sqrt, 0.5 / sqrt, derivatives);
460    }
461
462    /** {@inheritDoc} */
463    public SparseGradient cbrt() {
464        final double cbrt = FastMath.cbrt(value);
465        return new SparseGradient(cbrt, 1.0 / (3 * cbrt * cbrt), derivatives);
466    }
467
468    /** {@inheritDoc} */
469    public SparseGradient rootN(final int n) {
470        if (n == 2) {
471            return sqrt();
472        } else if (n == 3) {
473            return cbrt();
474        } else {
475            final double root = FastMath.pow(value, 1.0 / n);
476            return new SparseGradient(root, 1.0 / (n * FastMath.pow(root, n - 1)), derivatives);
477        }
478    }
479
480    /** {@inheritDoc} */
481    public SparseGradient pow(final double p) {
482        return new SparseGradient(FastMath.pow(value,  p), p * FastMath.pow(value,  p - 1), derivatives);
483    }
484
485    /** {@inheritDoc} */
486    public SparseGradient pow(final int n) {
487        if (n == 0) {
488            return getField().getOne();
489        } else {
490            final double valueNm1 = FastMath.pow(value,  n - 1);
491            return new SparseGradient(value * valueNm1, n * valueNm1, derivatives);
492        }
493    }
494
495    /** {@inheritDoc} */
496    public SparseGradient pow(final SparseGradient e) {
497        return log().multiply(e).exp();
498    }
499
500    /** Compute a<sup>x</sup> where a is a double and x a {@link SparseGradient}
501     * @param a number to exponentiate
502     * @param x power to apply
503     * @return a<sup>x</sup>
504     */
505    public static SparseGradient pow(final double a, final SparseGradient x) {
506        if (a == 0) {
507            if (x.value == 0) {
508                return x.compose(1.0, Double.NEGATIVE_INFINITY);
509            } else if (x.value < 0) {
510                return x.compose(Double.NaN, Double.NaN);
511            } else {
512                return x.getField().getZero();
513            }
514        } else {
515            final double ax = FastMath.pow(a, x.value);
516            return new SparseGradient(ax, ax * FastMath.log(a), x.derivatives);
517        }
518    }
519
520    /** {@inheritDoc} */
521    public SparseGradient exp() {
522        final double e = FastMath.exp(value);
523        return new SparseGradient(e, e, derivatives);
524    }
525
526    /** {@inheritDoc} */
527    public SparseGradient expm1() {
528        return new SparseGradient(FastMath.expm1(value), FastMath.exp(value), derivatives);
529    }
530
531    /** {@inheritDoc} */
532    public SparseGradient log() {
533        return new SparseGradient(FastMath.log(value), 1.0 / value, derivatives);
534    }
535
536    /** Base 10 logarithm.
537     * @return base 10 logarithm of the instance
538     */
539    public SparseGradient log10() {
540        return new SparseGradient(FastMath.log10(value), 1.0 / (FastMath.log(10.0) * value), derivatives);
541    }
542
543    /** {@inheritDoc} */
544    public SparseGradient log1p() {
545        return new SparseGradient(FastMath.log1p(value), 1.0 / (1.0 + value), derivatives);
546    }
547
548    /** {@inheritDoc} */
549    public SparseGradient cos() {
550        return new SparseGradient(FastMath.cos(value), -FastMath.sin(value), derivatives);
551    }
552
553    /** {@inheritDoc} */
554    public SparseGradient sin() {
555        return new SparseGradient(FastMath.sin(value), FastMath.cos(value), derivatives);
556    }
557
558    /** {@inheritDoc} */
559    public SparseGradient tan() {
560        final double t = FastMath.tan(value);
561        return new SparseGradient(t, 1 + t * t, derivatives);
562    }
563
564    /** {@inheritDoc} */
565    public SparseGradient acos() {
566        return new SparseGradient(FastMath.acos(value), -1.0 / FastMath.sqrt(1 - value * value), derivatives);
567    }
568
569    /** {@inheritDoc} */
570    public SparseGradient asin() {
571        return new SparseGradient(FastMath.asin(value), 1.0 / FastMath.sqrt(1 - value * value), derivatives);
572    }
573
574    /** {@inheritDoc} */
575    public SparseGradient atan() {
576        return new SparseGradient(FastMath.atan(value), 1.0 / (1 + value * value), derivatives);
577    }
578
579    /** {@inheritDoc} */
580    public SparseGradient atan2(final SparseGradient x) {
581
582        // compute r = sqrt(x^2+y^2)
583        final SparseGradient r = multiply(this).add(x.multiply(x)).sqrt();
584
585        final SparseGradient a;
586        if (x.value >= 0) {
587
588            // compute atan2(y, x) = 2 atan(y / (r + x))
589            a = divide(r.add(x)).atan().multiply(2);
590
591        } else {
592
593            // compute atan2(y, x) = +/- pi - 2 atan(y / (r - x))
594            final SparseGradient tmp = divide(r.subtract(x)).atan().multiply(-2);
595            a = tmp.add(tmp.value <= 0 ? -FastMath.PI : FastMath.PI);
596
597        }
598
599        // fix value to take special cases (+0/+0, +0/-0, -0/+0, -0/-0, +/-infinity) correctly
600        a.value = FastMath.atan2(value, x.value);
601
602        return a;
603
604    }
605
606    /** Two arguments arc tangent operation.
607     * @param y first argument of the arc tangent
608     * @param x second argument of the arc tangent
609     * @return atan2(y, x)
610     */
611    public static SparseGradient atan2(final SparseGradient y, final SparseGradient x) {
612        return y.atan2(x);
613    }
614
615    /** {@inheritDoc} */
616    public SparseGradient cosh() {
617        return new SparseGradient(FastMath.cosh(value), FastMath.sinh(value), derivatives);
618    }
619
620    /** {@inheritDoc} */
621    public SparseGradient sinh() {
622        return new SparseGradient(FastMath.sinh(value), FastMath.cosh(value), derivatives);
623    }
624
625    /** {@inheritDoc} */
626    public SparseGradient tanh() {
627        final double t = FastMath.tanh(value);
628        return new SparseGradient(t, 1 - t * t, derivatives);
629    }
630
631    /** {@inheritDoc} */
632    public SparseGradient acosh() {
633        return new SparseGradient(FastMath.acosh(value), 1.0 / FastMath.sqrt(value * value - 1.0), derivatives);
634    }
635
636    /** {@inheritDoc} */
637    public SparseGradient asinh() {
638        return new SparseGradient(FastMath.asinh(value), 1.0 / FastMath.sqrt(value * value + 1.0), derivatives);
639    }
640
641    /** {@inheritDoc} */
642    public SparseGradient atanh() {
643        return new SparseGradient(FastMath.atanh(value), 1.0 / (1.0 - value * value), derivatives);
644    }
645
646    /** Convert radians to degrees, with error of less than 0.5 ULP
647     *  @return instance converted into degrees
648     */
649    public SparseGradient toDegrees() {
650        return new SparseGradient(FastMath.toDegrees(value), FastMath.toDegrees(1.0), derivatives);
651    }
652
653    /** Convert degrees to radians, with error of less than 0.5 ULP
654     *  @return instance converted into radians
655     */
656    public SparseGradient toRadians() {
657        return new SparseGradient(FastMath.toRadians(value), FastMath.toRadians(1.0), derivatives);
658    }
659
660    /** Evaluate Taylor expansion of a sparse gradient.
661     * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
662     * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
663     */
664    public double taylor(final double ... delta) {
665        double y = value;
666        for (int i = 0; i < delta.length; ++i) {
667            y += delta[i] * getDerivative(i);
668        }
669        return y;
670    }
671
672    /** Compute composition of the instance by a univariate function.
673     * @param f0 value of the function at (i.e. f({@link #getValue()}))
674     * @param f1 first derivative of the function at
675     * the current point (i.e. f'({@link #getValue()}))
676     * @return f(this)
677    */
678    public SparseGradient compose(final double f0, final double f1) {
679        return new SparseGradient(f0, f1, derivatives);
680    }
681
682    /** {@inheritDoc} */
683    public SparseGradient linearCombination(final SparseGradient[] a,
684                                              final SparseGradient[] b)
685        throws DimensionMismatchException {
686
687        // compute a simple value, with all partial derivatives
688        SparseGradient out = a[0].getField().getZero();
689        for (int i = 0; i < a.length; ++i) {
690            out = out.add(a[i].multiply(b[i]));
691        }
692
693        // recompute an accurate value, taking care of cancellations
694        final double[] aDouble = new double[a.length];
695        for (int i = 0; i < a.length; ++i) {
696            aDouble[i] = a[i].getValue();
697        }
698        final double[] bDouble = new double[b.length];
699        for (int i = 0; i < b.length; ++i) {
700            bDouble[i] = b[i].getValue();
701        }
702        out.value = MathArrays.linearCombination(aDouble, bDouble);
703
704        return out;
705
706    }
707
708    /** {@inheritDoc} */
709    public SparseGradient linearCombination(final double[] a, final SparseGradient[] b) {
710
711        // compute a simple value, with all partial derivatives
712        SparseGradient out = b[0].getField().getZero();
713        for (int i = 0; i < a.length; ++i) {
714            out = out.add(b[i].multiply(a[i]));
715        }
716
717        // recompute an accurate value, taking care of cancellations
718        final double[] bDouble = new double[b.length];
719        for (int i = 0; i < b.length; ++i) {
720            bDouble[i] = b[i].getValue();
721        }
722        out.value = MathArrays.linearCombination(a, bDouble);
723
724        return out;
725
726    }
727
728    /** {@inheritDoc} */
729    public SparseGradient linearCombination(final SparseGradient a1, final SparseGradient b1,
730                                              final SparseGradient a2, final SparseGradient b2) {
731
732        // compute a simple value, with all partial derivatives
733        SparseGradient out = a1.multiply(b1).add(a2.multiply(b2));
734
735        // recompute an accurate value, taking care of cancellations
736        out.value = MathArrays.linearCombination(a1.value, b1.value, a2.value, b2.value);
737
738        return out;
739
740    }
741
742    /** {@inheritDoc} */
743    public SparseGradient linearCombination(final double a1, final SparseGradient b1,
744                                              final double a2, final SparseGradient b2) {
745
746        // compute a simple value, with all partial derivatives
747        SparseGradient out = b1.multiply(a1).add(b2.multiply(a2));
748
749        // recompute an accurate value, taking care of cancellations
750        out.value = MathArrays.linearCombination(a1, b1.value, a2, b2.value);
751
752        return out;
753
754    }
755
756    /** {@inheritDoc} */
757    public SparseGradient linearCombination(final SparseGradient a1, final SparseGradient b1,
758                                              final SparseGradient a2, final SparseGradient b2,
759                                              final SparseGradient a3, final SparseGradient b3) {
760
761        // compute a simple value, with all partial derivatives
762        SparseGradient out = a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3));
763
764        // recompute an accurate value, taking care of cancellations
765        out.value = MathArrays.linearCombination(a1.value, b1.value,
766                                                 a2.value, b2.value,
767                                                 a3.value, b3.value);
768
769        return out;
770
771    }
772
773    /** {@inheritDoc} */
774    public SparseGradient linearCombination(final double a1, final SparseGradient b1,
775                                              final double a2, final SparseGradient b2,
776                                              final double a3, final SparseGradient b3) {
777
778        // compute a simple value, with all partial derivatives
779        SparseGradient out = b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3));
780
781        // recompute an accurate value, taking care of cancellations
782        out.value = MathArrays.linearCombination(a1, b1.value,
783                                                 a2, b2.value,
784                                                 a3, b3.value);
785
786        return out;
787
788    }
789
790    /** {@inheritDoc} */
791    public SparseGradient linearCombination(final SparseGradient a1, final SparseGradient b1,
792                                              final SparseGradient a2, final SparseGradient b2,
793                                              final SparseGradient a3, final SparseGradient b3,
794                                              final SparseGradient a4, final SparseGradient b4) {
795
796        // compute a simple value, with all partial derivatives
797        SparseGradient out = a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3)).add(a4.multiply(b4));
798
799        // recompute an accurate value, taking care of cancellations
800        out.value = MathArrays.linearCombination(a1.value, b1.value,
801                                                 a2.value, b2.value,
802                                                 a3.value, b3.value,
803                                                 a4.value, b4.value);
804
805        return out;
806
807    }
808
809    /** {@inheritDoc} */
810    public SparseGradient linearCombination(final double a1, final SparseGradient b1,
811                                              final double a2, final SparseGradient b2,
812                                              final double a3, final SparseGradient b3,
813                                              final double a4, final SparseGradient b4) {
814
815        // compute a simple value, with all partial derivatives
816        SparseGradient out = b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3)).add(b4.multiply(a4));
817
818        // recompute an accurate value, taking care of cancellations
819        out.value = MathArrays.linearCombination(a1, b1.value,
820                                                 a2, b2.value,
821                                                 a3, b3.value,
822                                                 a4, b4.value);
823
824        return out;
825
826    }
827
828    /**
829     * Test for the equality of two sparse gradients.
830     * <p>
831     * Sparse gradients are considered equal if they have the same value
832     * and the same derivatives.
833     * </p>
834     * @param other Object to test for equality to this
835     * @return true if two sparse gradients are equal
836     */
837    @Override
838    public boolean equals(Object other) {
839
840        if (this == other) {
841            return true;
842        }
843
844        if (other instanceof SparseGradient) {
845            final SparseGradient rhs = (SparseGradient)other;
846            if (!Precision.equals(value, rhs.value, 1)) {
847                return false;
848            }
849            if (derivatives.size() != rhs.derivatives.size()) {
850                return false;
851            }
852            for (final Map.Entry<Integer, Double> entry : derivatives.entrySet()) {
853                if (!rhs.derivatives.containsKey(entry.getKey())) {
854                    return false;
855                }
856                if (!Precision.equals(entry.getValue(), rhs.derivatives.get(entry.getKey()), 1)) {
857                    return false;
858                }
859            }
860            return true;
861        }
862
863        return false;
864
865    }
866
867    /**
868     * Get a hashCode for the derivative structure.
869     * @return a hash code value for this object
870     * @since 3.2
871     */
872    @Override
873    public int hashCode() {
874        return 743 + 809 * MathUtils.hash(value) + 167 * derivatives.hashCode();
875    }
876
877}