Java Code Examples for org.apache.commons.math3.analysis.differentiation.DerivativeStructure#compose()

The following examples show how to use org.apache.commons.math3.analysis.differentiation.DerivativeStructure#compose() . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
Source File: HarmonicOscillator.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t) {
    final double x = t.getValue();
    double[] f = new double[t.getOrder() + 1];

    final double alpha = omega * x + phase;
    f[0] = amplitude * FastMath.cos(alpha);
    if (f.length > 1) {
        f[1] = -amplitude * omega * FastMath.sin(alpha);
        final double mo2 = - omega * omega;
        for (int i = 2; i < f.length; ++i) {
            f[i] = mo2 * f[i - 2];
        }
    }

    return t.compose(f);

}
 
Example 2
Source File: HarmonicOscillator.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {
    final double x = t.getValue();
    double[] f = new double[t.getOrder() + 1];

    final double alpha = omega * x + phase;
    f[0] = amplitude * FastMath.cos(alpha);
    if (f.length > 1) {
        f[1] = -amplitude * omega * FastMath.sin(alpha);
        final double mo2 = - omega * omega;
        for (int i = 2; i < f.length; ++i) {
            f[i] = mo2 * f[i - 2];
        }
    }

    return t.compose(f);

}
 
Example 3
Source File: Nopol2017_0064_s.java    From coming with MIT License 5 votes vote down vote up
/** Convert a {@link DifferentiableUnivariateFunction} into a {@link UnivariateDifferentiable}.
 * <p>
 * Note that the converted function is able to handle {@link DerivativeStructure} with
 * <em>only</em> one parameter and up to order one. If the function is called with
 * more parameters or higher order, a {@link DimensionMismatchException} will be thrown.
 * </p>
 * @param f function to convert
 * @return converted function
 * @deprecated this conversion method is temporary in version 3.1, as the {@link
 * DifferentiableUnivariateFunction} interface itself is deprecated
 */
@Deprecated
public static UnivariateDifferentiable toUnivariateDifferential(final DifferentiableUnivariateFunction f) {
    return new UnivariateDifferentiable() {

        /** {@inheritDoc} */
        public double value(final double x) {
            return f.value(x);
        }

        /** {@inheritDoc}
         * @exception DimensionMismatchException if number of parameters or derivation
         * order are higher than 1
         */
        public DerivativeStructure value(final DerivativeStructure t)
            throws DimensionMismatchException {
            if (t.getFreeParameters() != 1) {
                throw new DimensionMismatchException(t.getFreeParameters(), 1);
            }
            if (t.getOrder() > 1) {
                throw new DimensionMismatchException(t.getOrder(), 1);
            }
            return t.compose(new double[] {
                f.value(t.getValue()),
                f.derivative().value(t.getValue())
            });
        }

    };
}
 
Example 4
Source File: Nopol2017_0064_t.java    From coming with MIT License 5 votes vote down vote up
/** Convert a {@link DifferentiableUnivariateFunction} into a {@link UnivariateDifferentiable}.
 * <p>
 * Note that the converted function is able to handle {@link DerivativeStructure} with
 * <em>only</em> one parameter and up to order one. If the function is called with
 * more parameters or higher order, a {@link DimensionMismatchException} will be thrown.
 * </p>
 * @param f function to convert
 * @return converted function
 * @deprecated this conversion method is temporary in version 3.1, as the {@link
 * DifferentiableUnivariateFunction} interface itself is deprecated
 */
@Deprecated
public static UnivariateDifferentiable toUnivariateDifferential(final DifferentiableUnivariateFunction f) {
    return new UnivariateDifferentiable() {

        /** {@inheritDoc} */
        public double value(final double x) {
            return f.value(x);
        }

        /** {@inheritDoc}
         * @exception DimensionMismatchException if number of parameters or derivation
         * order are higher than 1
         */
        public DerivativeStructure value(final DerivativeStructure t)
            throws DimensionMismatchException {
            if (t.getFreeParameters() != 1) {
                throw new DimensionMismatchException(t.getFreeParameters(), 1);
            }
            if (t.getOrder() > 1) {
                throw new DimensionMismatchException(t.getOrder(), 1);
            }
            return t.compose(new double[] {
                f.value(t.getValue()),
                f.derivative().value(t.getValue())
            });
        }

    };
}
 
Example 5
Source File: Sigmoid.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    double[] f = new double[t.getOrder() + 1];
    final double exp = FastMath.exp(-t.getValue());
    if (Double.isInfinite(exp)) {

        // special handling near lower boundary, to avoid NaN
        f[0] = lo;
        Arrays.fill(f, 1, f.length, 0.0);

    } else {

        // the nth order derivative of sigmoid has the form:
        // dn(sigmoid(x)/dxn = P_n(exp(-x)) / (1+exp(-x))^(n+1)
        // where P_n(t) is a degree n polynomial with normalized higher term
        // P_0(t) = 1, P_1(t) = t, P_2(t) = t^2 - t, P_3(t) = t^3 - 4 t^2 + t...
        // the general recurrence relation for P_n is:
        // P_n(x) = n t P_(n-1)(t) - t (1 + t) P_(n-1)'(t)
        final double[] p = new double[f.length];

        final double inv   = 1 / (1 + exp);
        double coeff = hi - lo;
        for (int n = 0; n < f.length; ++n) {

            // update and evaluate polynomial P_n(t)
            double v = 0;
            p[n] = 1;
            for (int k = n; k >= 0; --k) {
                v = v * exp + p[k];
                if (k > 1) {
                    p[k - 1] = (n - k + 2) * p[k - 2] - (k - 1) * p[k - 1];
                } else {
                    p[0] = 0;
                }
            }

            coeff *= inv;
            f[n]   = coeff * v;

        }

        // fix function value
        f[0] += lo;

    }

    return t.compose(f);

}
 
Example 6
Source File: Logit.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 * @exception OutOfRangeException if parameter is outside of function domain
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws OutOfRangeException {
    final double x = t.getValue();
    if (x < lo || x > hi) {
        throw new OutOfRangeException(x, lo, hi);
    }
    double[] f = new double[t.getOrder() + 1];

    // function value
    f[0] = FastMath.log((x - lo) / (hi - x));

    if (Double.isInfinite(f[0])) {

        if (f.length > 1) {
            f[1] = Double.POSITIVE_INFINITY;
        }
        // fill the array with infinities
        // (for x close to lo the signs will flip between -inf and +inf,
        //  for x close to hi the signs will always be +inf)
        // this is probably overkill, since the call to compose at the end
        // of the method will transform most infinities into NaN ...
        for (int i = 2; i < f.length; ++i) {
            f[i] = f[i - 2];
        }

    } else {

        // function derivatives
        final double invL = 1.0 / (x - lo);
        double xL = invL;
        final double invH = 1.0 / (hi - x);
        double xH = invH;
        for (int i = 1; i < f.length; ++i) {
            f[i] = xL + xH;
            xL  *= -i * invL;
            xH  *=  i * invH;
        }
    }

    return t.compose(f);
}
 
Example 7
Source File: Gaussian.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    final double u = is * (t.getValue() - mean);
    double[] f = new double[t.getOrder() + 1];

    // the nth order derivative of the Gaussian has the form:
    // dn(g(x)/dxn = (norm / s^n) P_n(u) exp(-u^2/2) with u=(x-m)/s
    // where P_n(u) is a degree n polynomial with same parity as n
    // P_0(u) = 1, P_1(u) = -u, P_2(u) = u^2 - 1, P_3(u) = -u^3 + 3 u...
    // the general recurrence relation for P_n is:
    // P_n(u) = P_(n-1)'(u) - u P_(n-1)(u)
    // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
    final double[] p = new double[f.length];
    p[0] = 1;
    final double u2 = u * u;
    double coeff = norm * FastMath.exp(-0.5 * u2);
    if (coeff <= Precision.SAFE_MIN) {
        Arrays.fill(f, 0.0);
    } else {
        f[0] = coeff;
        for (int n = 1; n < f.length; ++n) {

            // update and evaluate polynomial P_n(x)
            double v = 0;
            p[n] = -p[n - 1];
            for (int k = n; k >= 0; k -= 2) {
                v = v * u2 + p[k];
                if (k > 2) {
                    p[k - 2] = (k - 1) * p[k - 1] - p[k - 3];
                } else if (k == 2) {
                    p[0] = p[1];
                }
            }
            if ((n & 0x1) == 1) {
                v *= u;
            }

            coeff *= is;
            f[n] = coeff * v;

        }
    }

    return t.compose(f);

}
 
Example 8
Source File: Sinc.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    final double scaledX  = (normalized ? FastMath.PI : 1) * t.getValue();
    final double scaledX2 = scaledX * scaledX;

    double[] f = new double[t.getOrder() + 1];

    if (FastMath.abs(scaledX) <= SHORTCUT) {

        for (int i = 0; i < f.length; ++i) {
            final int k = i / 2;
            if ((i & 0x1) == 0) {
                // even derivation order
                f[i] = (((k & 0x1) == 0) ? 1 : -1) *
                       (1.0 / (i + 1) - scaledX2 * (1.0 / (2 * i + 6) - scaledX2 / (24 * i + 120)));
            } else {
                // odd derivation order
                f[i] = (((k & 0x1) == 0) ? -scaledX : scaledX) *
                       (1.0 / (i + 2) - scaledX2 * (1.0 / (6 * i + 24) - scaledX2 / (120 * i + 720)));
            }
        }

    } else {

        final double inv = 1 / scaledX;
        final double cos = FastMath.cos(scaledX);
        final double sin = FastMath.sin(scaledX);

        f[0] = inv * sin;

        // the nth order derivative of sinc has the form:
        // dn(sinc(x)/dxn = [S_n(x) sin(x) + C_n(x) cos(x)] / x^(n+1)
        // where S_n(x) is an even polynomial with degree n-1 or n (depending on parity)
        // and C_n(x) is an odd polynomial with degree n-1 or n (depending on parity)
        // S_0(x) = 1, S_1(x) = -1, S_2(x) = -x^2 + 2, S_3(x) = 3x^2 - 6...
        // C_0(x) = 0, C_1(x) = x, C_2(x) = -2x, C_3(x) = -x^3 + 6x...
        // the general recurrence relations for S_n and C_n are:
        // S_n(x) = x S_(n-1)'(x) - n S_(n-1)(x) - x C_(n-1)(x)
        // C_n(x) = x C_(n-1)'(x) - n C_(n-1)(x) + x S_(n-1)(x)
        // as per polynomials parity, we can store both S_n and C_n in the same array
        final double[] sc = new double[f.length];
        sc[0] = 1;

        double coeff = inv;
        for (int n = 1; n < f.length; ++n) {

            double s = 0;
            double c = 0;

            // update and evaluate polynomials S_n(x) and C_n(x)
            final int kStart;
            if ((n & 0x1) == 0) {
                // even derivation order, S_n is degree n and C_n is degree n-1
                sc[n] = 0;
                kStart = n;
            } else {
                // odd derivation order, S_n is degree n-1 and C_n is degree n
                sc[n] = sc[n - 1];
                c = sc[n];
                kStart = n - 1;
            }

            // in this loop, k is always even
            for (int k = kStart; k > 1; k -= 2) {

                // sine part
                sc[k]     = (k - n) * sc[k] - sc[k - 1];
                s         = s * scaledX2 + sc[k];

                // cosine part
                sc[k - 1] = (k - 1 - n) * sc[k - 1] + sc[k -2];
                c         = c * scaledX2 + sc[k - 1];

            }
            sc[0] *= -n;
            s      = s * scaledX2 + sc[0];

            coeff *= inv;
            f[n]   = coeff * (s * sin + c * scaledX * cos);

        }

    }

    if (normalized) {
        double scale = FastMath.PI;
        for (int i = 1; i < f.length; ++i) {
            f[i]  *= scale;
            scale *= FastMath.PI;
        }
    }

    return t.compose(f);

}
 
Example 9
Source File: Logit.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 * @exception OutOfRangeException if parameter is outside of function domain
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws OutOfRangeException {
    final double x = t.getValue();
    if (x < lo || x > hi) {
        throw new OutOfRangeException(x, lo, hi);
    }
    double[] f = new double[t.getOrder() + 1];

    // function value
    f[0] = FastMath.log((x - lo) / (hi - x));

    if (Double.isInfinite(f[0])) {

        if (f.length > 1) {
            f[1] = Double.POSITIVE_INFINITY;
        }
        // fill the array with infinities
        // (for x close to lo the signs will flip between -inf and +inf,
        //  for x close to hi the signs will always be +inf)
        // this is probably overkill, since the call to compose at the end
        // of the method will transform most infinities into NaN ...
        for (int i = 2; i < f.length; ++i) {
            f[i] = f[i - 2];
        }

    } else {

        // function derivatives
        final double invL = 1.0 / (x - lo);
        double xL = invL;
        final double invH = 1.0 / (hi - x);
        double xH = invH;
        for (int i = 1; i < f.length; ++i) {
            f[i] = xL + xH;
            xL  *= -i * invL;
            xH  *=  i * invH;
        }
    }

    return t.compose(f);
}
 
Example 10
Source File: Gaussian.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    final double u = is * (t.getValue() - mean);
    double[] f = new double[t.getOrder() + 1];

    // the nth order derivative of the Gaussian has the form:
    // dn(g(x)/dxn = (norm / s^n) P_n(u) exp(-u^2/2) with u=(x-m)/s
    // where P_n(u) is a degree n polynomial with same parity as n
    // P_0(u) = 1, P_1(u) = -u, P_2(u) = u^2 - 1, P_3(u) = -u^3 + 3 u...
    // the general recurrence relation for P_n is:
    // P_n(u) = P_(n-1)'(u) - u P_(n-1)(u)
    // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
    final double[] p = new double[f.length];
    p[0] = 1;
    final double u2 = u * u;
    double coeff = norm * FastMath.exp(-0.5 * u2);
    if (coeff <= Precision.SAFE_MIN) {
        Arrays.fill(f, 0.0);
    } else {
        f[0] = coeff;
        for (int n = 1; n < f.length; ++n) {

            // update and evaluate polynomial P_n(x)
            double v = 0;
            p[n] = -p[n - 1];
            for (int k = n; k >= 0; k -= 2) {
                v = v * u2 + p[k];
                if (k > 2) {
                    p[k - 2] = (k - 1) * p[k - 1] - p[k - 3];
                } else if (k == 2) {
                    p[0] = p[1];
                }
            }
            if ((n & 0x1) == 1) {
                v *= u;
            }

            coeff *= is;
            f[n] = coeff * v;

        }
    }

    return t.compose(f);

}
 
Example 11
Source File: Gaussian.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    final double u = is * (t.getValue() - mean);
    double[] f = new double[t.getOrder() + 1];

    // the nth order derivative of the Gaussian has the form:
    // dn(g(x)/dxn = (norm / s^n) P_n(u) exp(-u^2/2) with u=(x-m)/s
    // where P_n(u) is a degree n polynomial with same parity as n
    // P_0(u) = 1, P_1(u) = -u, P_2(u) = u^2 - 1, P_3(u) = -u^3 + 3 u...
    // the general recurrence relation for P_n is:
    // P_n(u) = P_(n-1)'(u) - u P_(n-1)(u)
    // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
    final double[] p = new double[f.length];
    p[0] = 1;
    final double u2 = u * u;
    double coeff = norm * FastMath.exp(-0.5 * u2);
    if (coeff <= Precision.SAFE_MIN) {
        Arrays.fill(f, 0.0);
    } else {
        f[0] = coeff;
        for (int n = 1; n < f.length; ++n) {

            // update and evaluate polynomial P_n(x)
            double v = 0;
            p[n] = -p[n - 1];
            for (int k = n; k >= 0; k -= 2) {
                v = v * u2 + p[k];
                if (k > 2) {
                    p[k - 2] = (k - 1) * p[k - 1] - p[k - 3];
                } else if (k == 2) {
                    p[0] = p[1];
                }
            }
            if ((n & 0x1) == 1) {
                v *= u;
            }

            coeff *= is;
            f[n] = coeff * v;

        }
    }

    return t.compose(f);

}
 
Example 12
Source File: Gaussian.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t) {

    final double u = is * (t.getValue() - mean);
    double[] f = new double[t.getOrder() + 1];

    // the nth order derivative of the Gaussian has the form:
    // dn(g(x)/dxn = (norm / s^n) P_n(u) exp(-u^2/2) with u=(x-m)/s
    // where P_n(u) is a degree n polynomial with same parity as n
    // P_0(u) = 1, P_1(u) = -u, P_2(u) = u^2 - 1, P_3(u) = -u^3 + 3 u...
    // the general recurrence relation for P_n is:
    // P_n(u) = P_(n-1)'(u) - u P_(n-1)(u)
    // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
    final double[] p = new double[f.length];
    p[0] = 1;
    final double u2 = u * u;
    double coeff = norm * FastMath.exp(-0.5 * u2);
    if (coeff <= Precision.SAFE_MIN) {
        Arrays.fill(f, 0.0);
    } else {
        f[0] = coeff;
        for (int n = 1; n < f.length; ++n) {

            // update and evaluate polynomial P_n(x)
            double v = 0;
            p[n] = -p[n - 1];
            for (int k = n; k >= 0; k -= 2) {
                v = v * u2 + p[k];
                if (k > 2) {
                    p[k - 2] = (k - 1) * p[k - 1] - p[k - 3];
                } else if (k == 2) {
                    p[0] = p[1];
                }
            }
            if ((n & 0x1) == 1) {
                v *= u;
            }

            coeff *= is;
            f[n] = coeff * v;

        }
    }

    return t.compose(f);

}
 
Example 13
Source File: Sigmoid.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    double[] f = new double[t.getOrder() + 1];
    final double exp = FastMath.exp(-t.getValue());
    if (Double.isInfinite(exp)) {

        // special handling near lower boundary, to avoid NaN
        f[0] = lo;
        Arrays.fill(f, 1, f.length, 0.0);

    } else {

        // the nth order derivative of sigmoid has the form:
        // dn(sigmoid(x)/dxn = P_n(exp(-x)) / (1+exp(-x))^(n+1)
        // where P_n(t) is a degree n polynomial with normalized higher term
        // P_0(t) = 1, P_1(t) = t, P_2(t) = t^2 - t, P_3(t) = t^3 - 4 t^2 + t...
        // the general recurrence relation for P_n is:
        // P_n(x) = n t P_(n-1)(t) - t (1 + t) P_(n-1)'(t)
        final double[] p = new double[f.length];

        final double inv   = 1 / (1 + exp);
        double coeff = hi - lo;
        for (int n = 0; n < f.length; ++n) {

            // update and evaluate polynomial P_n(t)
            double v = 0;
            p[n] = 1;
            for (int k = n; k >= 0; --k) {
                v = v * exp + p[k];
                if (k > 1) {
                    p[k - 1] = (n - k + 2) * p[k - 2] - (k - 1) * p[k - 1];
                } else {
                    p[0] = 0;
                }
            }

            coeff *= inv;
            f[n]   = coeff * v;

        }

        // fix function value
        f[0] += lo;

    }

    return t.compose(f);

}
 
Example 14
Source File: Logit.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 * @exception OutOfRangeException if parameter is outside of function domain
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws OutOfRangeException {
    final double x = t.getValue();
    if (x < lo || x > hi) {
        throw new OutOfRangeException(x, lo, hi);
    }
    double[] f = new double[t.getOrder() + 1];

    // function value
    f[0] = FastMath.log((x - lo) / (hi - x));

    if (Double.isInfinite(f[0])) {

        if (f.length > 1) {
            f[1] = Double.POSITIVE_INFINITY;
        }
        // fill the array with infinities
        // (for x close to lo the signs will flip between -inf and +inf,
        //  for x close to hi the signs will always be +inf)
        // this is probably overkill, since the call to compose at the end
        // of the method will transform most infinities into NaN ...
        for (int i = 2; i < f.length; ++i) {
            f[i] = f[i - 2];
        }

    } else {

        // function derivatives
        final double invL = 1.0 / (x - lo);
        double xL = invL;
        final double invH = 1.0 / (hi - x);
        double xH = invH;
        for (int i = 1; i < f.length; ++i) {
            f[i] = xL + xH;
            xL  *= -i * invL;
            xH  *=  i * invH;
        }
    }
    
    return t.compose(f);

}
 
Example 15
Source File: Gaussian.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    final double u = is * (t.getValue() - mean);
    double[] f = new double[t.getOrder() + 1];

    // the nth order derivative of the Gaussian has the form:
    // dn(g(x)/dxn = (norm / s^n) P_n(u) exp(-u^2/2) with u=(x-m)/s
    // where P_n(u) is a degree n polynomial with same parity as n
    // P_0(u) = 1, P_1(u) = -u, P_2(u) = u^2 - 1, P_3(u) = -u^3 + 3 u...
    // the general recurrence relation for P_n is:
    // P_n(u) = P_(n-1)'(u) - u P_(n-1)(u)
    // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
    final double[] p = new double[f.length];
    p[0] = 1;
    final double u2 = u * u;
    double coeff = norm * FastMath.exp(-0.5 * u2);
    if (coeff <= Precision.SAFE_MIN) {
        Arrays.fill(f, 0.0);
    } else {
        f[0] = coeff;
        for (int n = 1; n < f.length; ++n) {

            // update and evaluate polynomial P_n(x)
            double v = 0;
            p[n] = -p[n - 1];
            for (int k = n; k >= 0; k -= 2) {
                v = v * u2 + p[k];
                if (k > 2) {
                    p[k - 2] = (k - 1) * p[k - 1] - p[k - 3];
                } else if (k == 2) {
                    p[0] = p[1];
                }
            }
            if ((n & 0x1) == 1) {
                v *= u;
            }

            coeff *= is;
            f[n] = coeff * v;

        }
    }

    return t.compose(f);

}
 
Example 16
Source File: Gaussian.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    final double u = is * (t.getValue() - mean);
    double[] f = new double[t.getOrder() + 1];

    // the nth order derivative of the Gaussian has the form:
    // dn(g(x)/dxn = (norm / s^n) P_n(u) exp(-u^2/2) with u=(x-m)/s
    // where P_n(u) is a degree n polynomial with same parity as n
    // P_0(u) = 1, P_1(u) = -u, P_2(u) = u^2 - 1, P_3(u) = -u^3 + 3 u...
    // the general recurrence relation for P_n is:
    // P_n(u) = P_(n-1)'(u) - u P_(n-1)(u)
    // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
    final double[] p = new double[f.length];
    p[0] = 1;
    final double u2 = u * u;
    double coeff = norm * FastMath.exp(-0.5 * u2);
    if (coeff <= Precision.SAFE_MIN) {
        Arrays.fill(f, 0.0);
    } else {
        f[0] = coeff;
        for (int n = 1; n < f.length; ++n) {

            // update and evaluate polynomial P_n(x)
            double v = 0;
            p[n] = -p[n - 1];
            for (int k = n; k >= 0; k -= 2) {
                v = v * u2 + p[k];
                if (k > 2) {
                    p[k - 2] = (k - 1) * p[k - 1] - p[k - 3];
                } else if (k == 2) {
                    p[0] = p[1];
                }
            }
            if ((n & 0x1) == 1) {
                v *= u;
            }

            coeff *= is;
            f[n] = coeff * v;

        }
    }

    return t.compose(f);

}
 
Example 17
Source File: Sinc.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    final double scaledX  = (normalized ? FastMath.PI : 1) * t.getValue();
    final double scaledX2 = scaledX * scaledX;

    double[] f = new double[t.getOrder() + 1];

    if (FastMath.abs(scaledX) <= SHORTCUT) {

        for (int i = 0; i < f.length; ++i) {
            final int k = i / 2;
            if ((i & 0x1) == 0) {
                // even derivation order
                f[i] = (((k & 0x1) == 0) ? 1 : -1) *
                       (1.0 / (i + 1) - scaledX2 * (1.0 / (2 * i + 6) - scaledX2 / (24 * i + 120)));
            } else {
                // odd derivation order
                f[i] = (((k & 0x1) == 0) ? -scaledX : scaledX) *
                       (1.0 / (i + 2) - scaledX2 * (1.0 / (6 * i + 24) - scaledX2 / (120 * i + 720)));
            }
        }

    } else {

        final double inv = 1 / scaledX;
        final double cos = FastMath.cos(scaledX);
        final double sin = FastMath.sin(scaledX);

        f[0] = inv * sin;

        // the nth order derivative of sinc has the form:
        // dn(sinc(x)/dxn = [S_n(x) sin(x) + C_n(x) cos(x)] / x^(n+1)
        // where S_n(x) is an even polynomial with degree n-1 or n (depending on parity)
        // and C_n(x) is an odd polynomial with degree n-1 or n (depending on parity)
        // S_0(x) = 1, S_1(x) = -1, S_2(x) = -x^2 + 2, S_3(x) = 3x^2 - 6...
        // C_0(x) = 0, C_1(x) = x, C_2(x) = -2x, C_3(x) = -x^3 + 6x...
        // the general recurrence relations for S_n and C_n are:
        // S_n(x) = x S_(n-1)'(x) - n S_(n-1)(x) - x C_(n-1)(x)
        // C_n(x) = x C_(n-1)'(x) - n C_(n-1)(x) + x S_(n-1)(x)
        // as per polynomials parity, we can store both S_n and C_n in the same array
        final double[] sc = new double[f.length];
        sc[0] = 1;

        double coeff = inv;
        for (int n = 1; n < f.length; ++n) {

            double s = 0;
            double c = 0;

            // update and evaluate polynomials S_n(x) and C_n(x)
            final int kStart;
            if ((n & 0x1) == 0) {
                // even derivation order, S_n is degree n and C_n is degree n-1
                sc[n] = 0;
                kStart = n;
            } else {
                // odd derivation order, S_n is degree n-1 and C_n is degree n
                sc[n] = sc[n - 1];
                c = sc[n];
                kStart = n - 1;
            }

            // in this loop, k is always even
            for (int k = kStart; k > 1; k -= 2) {

                // sine part
                sc[k]     = (k - n) * sc[k] - sc[k - 1];
                s         = s * scaledX2 + sc[k];

                // cosine part
                sc[k - 1] = (k - 1 - n) * sc[k - 1] + sc[k -2];
                c         = c * scaledX2 + sc[k - 1];

            }
            sc[0] *= -n;
            s      = s * scaledX2 + sc[0];

            coeff *= inv;
            f[n]   = coeff * (s * sin + c * scaledX * cos);

        }

    }

    if (normalized) {
        double scale = FastMath.PI;
        for (int i = 1; i < f.length; ++i) {
            f[i]  *= scale;
            scale *= FastMath.PI;
        }
    }

    return t.compose(f);

}
 
Example 18
Source File: Logit.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 * @exception OutOfRangeException if parameter is outside of function domain
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws OutOfRangeException {
    final double x = t.getValue();
    if (x < lo || x > hi) {
        throw new OutOfRangeException(x, lo, hi);
    }
    double[] f = new double[t.getOrder() + 1];

    // function value
    f[0] = FastMath.log((x - lo) / (hi - x));

    if (Double.isInfinite(f[0])) {

        if (f.length > 1) {
            f[1] = Double.POSITIVE_INFINITY;
        }
        // fill the array with infinities
        // (for x close to lo the signs will flip between -inf and +inf,
        //  for x close to hi the signs will always be +inf)
        // this is probably overkill, since the call to compose at the end
        // of the method will transform most infinities into NaN ...
        for (int i = 2; i < f.length; ++i) {
            f[i] = f[i - 2];
        }

    } else {

        // function derivatives
        final double invL = 1.0 / (x - lo);
        double xL = invL;
        final double invH = 1.0 / (hi - x);
        double xH = invH;
        for (int i = 1; i < f.length; ++i) {
            f[i] = xL + xH;
            xL  *= -i * invL;
            xH  *=  i * invH;
        }
    }

    return t.compose(f);
}
 
Example 19
Source File: Sigmoid.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    double[] f = new double[t.getOrder() + 1];
    final double exp = FastMath.exp(-t.getValue());
    if (Double.isInfinite(exp)) {

        // special handling near lower boundary, to avoid NaN
        f[0] = lo;
        Arrays.fill(f, 1, f.length, 0.0);

    } else {

        // the nth order derivative of sigmoid has the form:
        // dn(sigmoid(x)/dxn = P_n(exp(-x)) / (1+exp(-x))^(n+1)
        // where P_n(t) is a degree n polynomial with normalized higher term
        // P_0(t) = 1, P_1(t) = t, P_2(t) = t^2 - t, P_3(t) = t^3 - 4 t^2 + t...
        // the general recurrence relation for P_n is:
        // P_n(x) = n t P_(n-1)(t) - t (1 + t) P_(n-1)'(t)
        final double[] p = new double[f.length];

        final double inv   = 1 / (1 + exp);
        double coeff = hi - lo;
        for (int n = 0; n < f.length; ++n) {

            // update and evaluate polynomial P_n(t)
            double v = 0;
            p[n] = 1;
            for (int k = n; k >= 0; --k) {
                v = v * exp + p[k];
                if (k > 1) {
                    p[k - 1] = (n - k + 2) * p[k - 2] - (k - 1) * p[k - 1];
                } else {
                    p[0] = 0;
                }
            }

            coeff *= inv;
            f[n]   = coeff * v;

        }

        // fix function value
        f[0] += lo;

    }

    return t.compose(f);

}
 
Example 20
Source File: Sinc.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc}
 * @since 3.1
 */
public DerivativeStructure value(final DerivativeStructure t)
    throws DimensionMismatchException {

    final double scaledX  = (normalized ? FastMath.PI : 1) * t.getValue();
    final double scaledX2 = scaledX * scaledX;

    double[] f = new double[t.getOrder() + 1];

    if (FastMath.abs(scaledX) <= SHORTCUT) {

        for (int i = 0; i < f.length; ++i) {
            final int k = i / 2;
            if ((i & 0x1) == 0) {
                // even derivation order
                f[i] = (((k & 0x1) == 0) ? 1 : -1) *
                       (1.0 / (i + 1) - scaledX2 * (1.0 / (2 * i + 6) - scaledX2 / (24 * i + 120)));
            } else {
                // odd derivation order
                f[i] = (((k & 0x1) == 0) ? -scaledX : scaledX) *
                       (1.0 / (i + 2) - scaledX2 * (1.0 / (6 * i + 24) - scaledX2 / (120 * i + 720)));
            }
        }

    } else {

        final double inv = 1 / scaledX;
        final double cos = FastMath.cos(scaledX);
        final double sin = FastMath.sin(scaledX);

        f[0] = inv * sin;

        // the nth order derivative of sinc has the form:
        // dn(sinc(x)/dxn = [S_n(x) sin(x) + C_n(x) cos(x)] / x^(n+1)
        // where S_n(x) is an even polynomial with degree n-1 or n (depending on parity)
        // and C_n(x) is an odd polynomial with degree n-1 or n (depending on parity)
        // S_0(x) = 1, S_1(x) = -1, S_2(x) = -x^2 + 2, S_3(x) = 3x^2 - 6...
        // C_0(x) = 0, C_1(x) = x, C_2(x) = -2x, C_3(x) = -x^3 + 6x...
        // the general recurrence relations for S_n and C_n are:
        // S_n(x) = x S_(n-1)'(x) - n S_(n-1)(x) - x C_(n-1)(x)
        // C_n(x) = x C_(n-1)'(x) - n C_(n-1)(x) + x S_(n-1)(x)
        // as per polynomials parity, we can store both S_n and C_n in the same array
        final double[] sc = new double[f.length];
        sc[0] = 1;

        double coeff = inv;
        for (int n = 1; n < f.length; ++n) {

            double s = 0;
            double c = 0;

            // update and evaluate polynomials S_n(x) and C_n(x)
            final int kStart;
            if ((n & 0x1) == 0) {
                // even derivation order, S_n is degree n and C_n is degree n-1
                sc[n] = 0;
                kStart = n;
            } else {
                // odd derivation order, S_n is degree n-1 and C_n is degree n
                sc[n] = sc[n - 1];
                c = sc[n];
                kStart = n - 1;
            }

            // in this loop, k is always even
            for (int k = kStart; k > 1; k -= 2) {

                // sine part
                sc[k]     = (k - n) * sc[k] - sc[k - 1];
                s         = s * scaledX2 + sc[k];

                // cosine part
                sc[k - 1] = (k - 1 - n) * sc[k - 1] + sc[k -2];
                c         = c * scaledX2 + sc[k - 1];

            }
            sc[0] *= -n;
            s      = s * scaledX2 + sc[0];

            coeff *= inv;
            f[n]   = coeff * (s * sin + c * scaledX * cos);

        }

    }

    if (normalized) {
        double scale = FastMath.PI;
        for (int i = 1; i < f.length; ++i) {
            f[i]  *= scale;
            scale *= FastMath.PI;
        }
    }

    return t.compose(f);

}