Java Code Examples for org.apache.commons.math3.linear.RealVector#dotProduct()

/**
 * @return a comparator for sorting the optima.
 */
private Comparator<PointVectorValuePair> getPairComparator() {
    return new Comparator<PointVectorValuePair>() {
        private final RealVector target = new ArrayRealVector(optimizer.getTarget(), false);
        private final RealMatrix weight = optimizer.getWeight();

        public int compare(final PointVectorValuePair o1,
                           final PointVectorValuePair o2) {
            if (o1 == null) {
                return (o2 == null) ? 0 : 1;
            } else if (o2 == null) {
                return -1;
            }
            return Double.compare(weightedResidual(o1),
                                  weightedResidual(o2));
        }

        private double weightedResidual(final PointVectorValuePair pv) {
            final RealVector v = new ArrayRealVector(pv.getValueRef(), false);
            final RealVector r = target.subtract(v);
            return r.dotProduct(weight.operate(r));
        }
    };
}
 

public double[] predict(LearningInstance instance) {
    RealMatrix A = variableSpace.getMatrixVarByName(LinearUCBKey.A.get());
    RealVector B = variableSpace.getScalarVarByName(LinearUCBKey.B.get());
    RealMatrix invA = new LUDecomposition(A).getSolver().getInverse();
    RealVector theta = invA.operate(B);
    RealVector x = extractDenseVector(theta.getDimension(), instance);
    double bound = Math.sqrt(x.dotProduct(invA.operate(x)));
    double mean = x.dotProduct(theta);
    double pred = mean + alpha * bound;
    if (Double.isNaN(pred)) {
        logger.error("Prediction is NaN, model parameter A probably goes wrong.");
        pred = 0.0;
    }
    double[] preds = new double[1];
    preds[0] = pred;
    return preds;
}
 

/**
 * @return a comparator for sorting the optima.
 */
private Comparator<PointVectorValuePair> getPairComparator() {
    return new Comparator<PointVectorValuePair>() {
        private final RealVector target = new ArrayRealVector(optimizer.getTarget(), false);
        private final RealMatrix weight = optimizer.getWeight();

        public int compare(final PointVectorValuePair o1,
                           final PointVectorValuePair o2) {
            if (o1 == null) {
                return (o2 == null) ? 0 : 1;
            } else if (o2 == null) {
                return -1;
            }
            return Double.compare(weightedResidual(o1),
                                  weightedResidual(o2));
        }

        private double weightedResidual(final PointVectorValuePair pv) {
            final RealVector v = new ArrayRealVector(pv.getValueRef(), false);
            final RealVector r = target.subtract(v);
            return r.dotProduct(weight.operate(r));
        }
    };
}
 

/**
 * @return a comparator for sorting the optima.
 */
private Comparator<PointVectorValuePair> getPairComparator() {
    return new Comparator<PointVectorValuePair>() {
        private final RealVector target = new ArrayRealVector(optimizer.getTarget(), false);
        private final RealMatrix weight = optimizer.getWeight();

        public int compare(final PointVectorValuePair o1,
                           final PointVectorValuePair o2) {
            if (o1 == null) {
                return (o2 == null) ? 0 : 1;
            } else if (o2 == null) {
                return -1;
            }
            return Double.compare(weightedResidual(o1),
                                  weightedResidual(o2));
        }

        private double weightedResidual(final PointVectorValuePair pv) {
            final RealVector v = new ArrayRealVector(pv.getValueRef(), false);
            final RealVector r = target.subtract(v);
            return r.dotProduct(weight.operate(r));
        }
    };
}
 

@Override
protected double computeTSS(RealVector y) {
    if (!hasIntercept()) {
        return y.dotProduct(y);
    } else {
        final C regressand = getRegressand();
        final double sumOfWeights = weights.stats().sum().doubleValue();
        final Array<Double> yValues = Array.of(frame().col(regressand).toDoubleStream().toArray());
        final double weightedAvg = yValues.mapToDoubles(v -> v.getDouble() * weights.getDouble(v.index())).stats().sum().doubleValue() / sumOfWeights;
        final Array<Double> diffSquared = yValues.mapToDoubles(v -> weights.getDouble(v.index()) * Math.pow(v.getDouble() - weightedAvg, 2d));
        return diffSquared.stats().sum().doubleValue();
    }
}
 

private TrcWaypoint interpolatePoints(TrcWaypoint prev, TrcWaypoint point, double robotX, double robotY)
{
    // Find intersection of path segment with circle with radius followingDistance and center at robot
    RealVector start = new ArrayRealVector(new double[] { prev.x, prev.y });
    RealVector end = new ArrayRealVector(new double[] { point.x, point.y });
    RealVector robot = new ArrayRealVector(new double[] { robotX, robotY });

    RealVector startToEnd = end.subtract(start);
    RealVector robotToStart = start.subtract(robot);
    // Solve quadratic formula
    double a = startToEnd.dotProduct(startToEnd);
    double b = 2 * robotToStart.dotProduct(startToEnd);
    double c = robotToStart.dotProduct(robotToStart) - followingDistance * followingDistance;

    double discriminant = b * b - 4 * a * c;
    if (discriminant < 0)
    {
        // No valid intersection.
        return null;
    }
    else
    {
        // line is a parametric equation, where t=0 is start, t=1 is end.
        discriminant = Math.sqrt(discriminant);
        double t1 = (-b - discriminant) / (2 * a);
        double t2 = (-b + discriminant) / (2 * a);
        double t = Math.max(t1, t2); // We want the furthest intersection
        // If the intersection is not on the line segment, it's invalid.
        if (!TrcUtil.inRange(t, 0.0, 1.0))
        {
            return null;
        }
        return interpolate(prev, point, t);
    }
}
 

/**
 * pdf(x, x_hat) = exp(-0.5 * (x-x_hat) * inv(Σ) * (x-x_hat)T) / ( 2π^0.5d * det(Σ)^0.5)
 * 
 * @return value of probabilistic density function
 * @link https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Density_function
 */
public static double pdf(@Nonnull final RealVector x, @Nonnull final RealVector x_hat,
        @Nonnull final RealMatrix sigma) {
    final int dim = x.getDimension();
    Preconditions.checkArgument(x_hat.getDimension() == dim,
        "|x| != |x_hat|, |x|=" + dim + ", |x_hat|=" + x_hat.getDimension());
    Preconditions.checkArgument(sigma.getRowDimension() == dim,
        "|x| != |sigma|, |x|=" + dim + ", |sigma|=" + sigma.getRowDimension());
    Preconditions.checkArgument(sigma.isSquare(), "Sigma is not square matrix");

    LUDecomposition LU = new LUDecomposition(sigma);
    final double detSigma = LU.getDeterminant();
    double denominator = Math.pow(2.d * Math.PI, 0.5d * dim) * Math.pow(detSigma, 0.5d);
    if (denominator == 0.d) { // avoid divide by zero
        return 0.d;
    }

    final RealMatrix invSigma;
    DecompositionSolver solver = LU.getSolver();
    if (solver.isNonSingular() == false) {
        SingularValueDecomposition svd = new SingularValueDecomposition(sigma);
        invSigma = svd.getSolver().getInverse(); // least square solution
    } else {
        invSigma = solver.getInverse();
    }
    //EigenDecomposition eigen = new EigenDecomposition(sigma);
    //double detSigma = eigen.getDeterminant();
    //RealMatrix invSigma = eigen.getSolver().getInverse();

    RealVector diff = x.subtract(x_hat);
    RealVector premultiplied = invSigma.preMultiply(diff);
    double sum = premultiplied.dotProduct(diff);
    double numerator = Math.exp(-0.5d * sum);

    return numerator / denominator;
}
 

/**
 * Lanczos tridiagonalization for a symmetric matrix C to make s * s tridiagonal matrix T.
 *
 * @see http://www.cas.mcmaster.ca/~qiao/publications/spie05.pdf
 * @param C target symmetric matrix
 * @param a initial vector
 * @param T result is stored here
 */
public static void lanczosTridiagonalization(@Nonnull final RealMatrix C,
        @Nonnull final double[] a, @Nonnull final RealMatrix T) {
    Preconditions.checkArgument(Arrays.deepEquals(C.getData(), C.transpose().getData()),
        "Target matrix C must be a symmetric matrix");
    Preconditions.checkArgument(C.getColumnDimension() == a.length,
        "Column size of A and length of a should be same");
    Preconditions.checkArgument(T.getRowDimension() == T.getColumnDimension(),
        "T must be a square matrix");

    int s = T.getRowDimension();

    // initialize T with zeros
    T.setSubMatrix(new double[s][s], 0, 0);

    RealVector a0 = new ArrayRealVector(a.length);
    RealVector r = new ArrayRealVector(a);

    double beta0 = 1.d;

    for (int i = 0; i < s; i++) {
        RealVector a1 = r.mapDivide(beta0);
        RealVector Ca1 = C.operate(a1);

        double alpha1 = a1.dotProduct(Ca1);

        r = Ca1.add(a1.mapMultiply(-1.d * alpha1)).add(a0.mapMultiply(-1.d * beta0));

        double beta1 = r.getNorm();

        T.setEntry(i, i, alpha1);
        if (i - 1 >= 0) {
            T.setEntry(i, i - 1, beta0);
        }
        if (i + 1 < s) {
            T.setEntry(i, i + 1, beta1);
        }

        a0 = a1.copy();
        beta0 = beta1;
    }
}
 

/**
 * Calculates the estimated variance of the error term using the formula
 * <pre>
 *  Var(u) = Tr(u' Omega^-1 u)/(n-k)
 * </pre>
 * where n and k are the row and column dimensions of the design
 * matrix X.
 *
 * @return error variance
 * @since 2.2
 */
@Override
protected double calculateErrorVariance() {
    RealVector residuals = calculateResiduals();
    double t = residuals.dotProduct(getOmegaInverse().operate(residuals));
    return t / (getX().getRowDimension() - getX().getColumnDimension());

}
 

/**
 * Calculates the estimated variance of the error term using the formula
 * <pre>
 *  Var(u) = Tr(u' Omega^-1 u)/(n-k)
 * </pre>
 * where n and k are the row and column dimensions of the design
 * matrix X.
 *
 * @return error variance
 * @since 2.2
 */
@Override
protected double calculateErrorVariance() {
    RealVector residuals = calculateResiduals();
    double t = residuals.dotProduct(getOmegaInverse().operate(residuals));
    return t / (getX().getRowDimension() - getX().getColumnDimension());

}
 

/**
 * Calculates the estimated variance of the error term using the formula
 * <pre>
 *  Var(u) = Tr(u' Omega^-1 u)/(n-k)
 * </pre>
 * where n and k are the row and column dimensions of the design
 * matrix X.
 *
 * @return error variance
 * @since 2.2
 */
@Override
protected double calculateErrorVariance() {
    RealVector residuals = calculateResiduals();
    double t = residuals.dotProduct(getOmegaInverse().operate(residuals));
    return t / (getX().getRowDimension() - getX().getColumnDimension());

}
 

/**
 * <p>Calculates the variance of the error term.</p>
 * Uses the formula <pre>
 * var(u) = u &middot; u / (n - k)
 * </pre>
 * where n and k are the row and column dimensions of the design
 * matrix X.
 *
 * @return error variance estimate
 * @since 2.2
 */
protected double calculateErrorVariance() {
    RealVector residuals = calculateResiduals();
    return residuals.dotProduct(residuals) /
           (xMatrix.getRowDimension() - xMatrix.getColumnDimension());
}
 

/**
 * Returns the sum of squared residuals.
 *
 * @return residual sum of squares
 * @since 2.2
 */
public double calculateResidualSumOfSquares() {
    final RealVector residuals = calculateResiduals();
    // No advertised DME, args are valid
    return residuals.dotProduct(residuals);
}
 

/**
 * <p>Calculates the variance of the error term.</p>
 * Uses the formula <pre>
 * var(u) = u &middot; u / (n - k)
 * </pre>
 * where n and k are the row and column dimensions of the design
 * matrix X.
 *
 * @return error variance estimate
 * @since 2.2
 */
protected double calculateErrorVariance() {
    RealVector residuals = calculateResiduals();
    return residuals.dotProduct(residuals) /
           (xMatrix.getRowDimension() - xMatrix.getColumnDimension());
}
 

/**
 * Compute the cosine of the angle between 2 vectors.
 *
 * @param v Vector.
 * @param w Vector.
 * @return the cosine of the angle between {@code v} and {@code w}.
 */
private double cosAngle(final RealVector v, final RealVector w) {
    return v.dotProduct(w) / (v.getNorm() * w.getNorm());
}
 

/**
 * Returns the sum of squared residuals.
 *
 * @return residual sum of squares
 * @since 2.2
 */
public double calculateResidualSumOfSquares() {
    final RealVector residuals = calculateResiduals();
    return residuals.dotProduct(residuals);
}
 

/**
 * Compute the cosine of the angle between 2 vectors.
 *
 * @param v Vector.
 * @param w Vector.
 * @return the cosine of the angle between {@code v} and {@code w}.
 */
private double cosAngle(final RealVector v, final RealVector w) {
    return v.dotProduct(w) / (v.getNorm() * w.getNorm());
}
 

/**
 * <p>Calculates the variance of the error term.</p>
 * Uses the formula <pre>
 * var(u) = u &middot; u / (n - k)
 * </pre>
 * where n and k are the row and column dimensions of the design
 * matrix X.
 *
 * @return error variance estimate
 * @since 2.2
 */
protected double calculateErrorVariance() {
    RealVector residuals = calculateResiduals();
    return residuals.dotProduct(residuals) /
           (xMatrix.getRowDimension() - xMatrix.getColumnDimension());
}
 

/**
 * <p>Calculates the variance of the error term.</p>
 * Uses the formula <pre>
 * var(u) = u &middot; u / (n - k)
 * </pre>
 * where n and k are the row and column dimensions of the design
 * matrix X.
 *
 * @return error variance estimate
 * @since 2.2
 */
protected double calculateErrorVariance() {
    RealVector residuals = calculateResiduals();
    return residuals.dotProduct(residuals) /
           (xMatrix.getRowDimension() - xMatrix.getColumnDimension());
}
 

/**
 * Compute the cosine of the angle between 2 vectors.
 *
 * @param v Vector.
 * @param w Vector.
 * @return the cosine of the angle between {@code v} and {@code w}.
 */
private double cosAngle(final RealVector v, final RealVector w) {
    return v.dotProduct(w) / (v.getNorm() * w.getNorm());
}