Java Code Examples for org.ejml.ops.CommonOps#scale()

private MatrixFormulation() {
    int numRows = response.length;
    int numCols = predictors.length + ((hasIntercept) ? 1 : 0);
    this.X = createMatrixA(numRows, numCols);
    this.Xt = new DenseMatrix64F(numCols, numRows);
    CommonOps.transpose(X, Xt);
    this.XtXInv = new DenseMatrix64F(numCols, numCols);
    this.b = new DenseMatrix64F(numCols, 1);
    this.y = new DenseMatrix64F(numRows, 1);
    solveSystem(numRows, numCols);
    this.fitted = computeFittedValues();
    this.residuals = computeResiduals();
    this.sigma2 = estimateSigma2(numCols);
    this.covarianceMatrix = new DenseMatrix64F(numCols, numCols);
    CommonOps.scale(sigma2, XtXInv, covarianceMatrix);
}
 

private void computeIOUActualizedDisplacement(final double[] optimalRates,
                                              final int offset,
                                              final int pio,
                                              double branchLength,
                                              DenseMatrix64F inverseSelectionStrength) {
    DenseMatrix64F displacementOU = wrap(displacements, pio, dimProcess, 1);
    DenseMatrix64F optVal = wrap(optimalRates, offset, dimProcess, 1);
    DenseMatrix64F displacement = new DenseMatrix64F(dimProcess, 1);

    CommonOps.mult(inverseSelectionStrength, displacementOU, displacement);
    CommonOps.scale(-1.0, displacement);

    CommonOps.addEquals(displacement, branchLength, optVal);

    unwrap(displacement, displacements, pio + dimProcess);
}
 

private static void rotateNd(double[] x, int dim) {

        // Get first `dim` locations
        DenseMatrix64F matrix = new DenseMatrix64F(dim, dim);
        for (int row = 0; row < dim; ++row) {
            for (int col = 0; col < dim; ++col) {
                matrix.set(row, col, x[col * dim + row]);
            }
        }

        // Do a QR decomposition
        QRDecomposition<DenseMatrix64F> qr = DecompositionFactory.qr(dim, dim);
        qr.decompose(matrix);
        DenseMatrix64F qm = qr.getQ(null, true);
        DenseMatrix64F rm = qr.getR(null, true);

        // Reflection invariance
        if (rm.get(0,0) < 0) {
            CommonOps.scale(-1, rm);
            CommonOps.scale(-1, qm);
        }

        // Compute Q^{-1}
        DenseMatrix64F qInv = new DenseMatrix64F(dim, dim);
        CommonOps.transpose(qm, qInv);

        // Apply to each location
        for (int location = 0; location < x.length / dim; ++location) {
            WrappedVector locationVector = new WrappedVector.Raw(x, location * dim, dim);
            MissingOps.matrixVectorMultiple(qInv, locationVector, locationVector, dim);
        }
    }
 

DenseMatrix64F getPrecisionBranch(double branchPrecision){
    if (!hasDrift) {
        DenseMatrix64F P1 = new DenseMatrix64F(dimTrait, dimTrait);
        CommonOps.scale(branchPrecision, Pd, P1);
        return P1;
    } else {
        return DenseMatrix64F.wrap(dimTrait, dimTrait, precisionBuffer);
    }
}
 

DenseMatrix64F getVarianceBranch(double branchPrecision){
    if (!hasDrift) {
        final DenseMatrix64F V1 = new DenseMatrix64F(dimTrait, dimTrait);
        CommonOps.scale(1.0 / branchPrecision, Vd, V1);
        return V1;
    } else {
        DenseMatrix64F P = getPrecisionBranch(branchPrecision);
        DenseMatrix64F V = new DenseMatrix64F(dimTrait, dimTrait);
        CommonOps.invert(P, V);
        return V;
    }
}
 

private DenseMatrix64F scaleGradient(double scalar, DenseMatrix64F gradient) {
    DenseMatrix64F result = gradient.copy();
    if (scalar == 0.0) {
        CommonOps.fill(result, 0.0);
    } else {
        CommonOps.scale(scalar, result);
    }
    return result;
}
 

public void makeGradientMatrices0(DenseMatrix64F matrix1, DenseMatrix64F logDetComponent,
                                  BranchSufficientStatistics statistics, double differentialScaling) {

    final NormalSufficientStatistics above = statistics.getAbove();
    final NormalSufficientStatistics branch = statistics.getBranch();

    DenseMatrix64F Qi = above.getRawPrecision();
    CommonOps.scale(differentialScaling, branch.getRawVariance(), matrix1); //matrix1 = Si
    CommonOps.mult(Qi, matrix1, logDetComponent); //matrix0 = logDetComponent
    CommonOps.mult(logDetComponent, Qi, matrix1); //matrix1 = QuadraticComponent

}
 

@Override
public double[] chainRule(BranchSufficientStatistics statistics, NodeRef node,
                          DenseMatrix64F gradQInv, DenseMatrix64F gradN) {

    final double rate = branchRateModel.getBranchRate(tree, node);
    final double differential = branchRateModel.getBranchRateDifferential(tree, node);
    final double scaling = differential / rate;

    // Q_i w.r.t. rate
    DenseMatrix64F gradMatQInv = matrixJacobianQInv;
    CommonOps.scale(scaling, statistics.getBranch().getRawVariance(), gradMatQInv);

    double[] gradient = new double[1];
    for (int i = 0; i < gradMatQInv.getNumElements(); i++) {
        gradient[0] += gradMatQInv.get(i) * gradQInv.get(i);
    }

    // n_i w.r.t. rate
    // TODO: Fix delegate to (possibly) un-link drift from arbitrary rate
    DenseMatrix64F gradMatN = matrixJacobianN;
    CommonOps.scale(scaling, statistics.getBranch().getRawDisplacement(), gradMatN);
    for (int i = 0; i < gradMatN.numRows; i++) {
        gradient[0] += gradMatN.get(i) * gradN.get(i);
    }

    return gradient;

}
 

private void actualizeDisplacementGradient(ContinuousDiffusionIntegrator cdi,
                                               int nodeIndex, DenseMatrix64F gradient) {
        // q_i
        double[] qi = new double[dim * dim];
        cdi.getBranch1mActualization(getMatrixBufferOffsetIndex(nodeIndex), qi);
        DenseMatrix64F Actu = wrap(qi, 0, dim, dim);
        CommonOps.scale(-1.0, Actu);
//        for (int i = 0; i < dim; i++) {
//            Actu.unsafe_set(i, i, Actu.unsafe_get(i, i) - 1.0);
//        }
        DenseMatrix64F tmp = new DenseMatrix64F(dim, 1);
        CommonOps.mult(Actu, gradient, tmp);
        CommonOps.scale(-1.0, tmp, gradient);
    }
 

private static void idMinusA(DenseMatrix64F A) {
    CommonOps.scale(-1.0, A);
    for (int i = 0; i < A.numCols; i++) {
        A.set(i, i, 1.0 + A.get(i, i));
    }
}
 

void actualizePrecision(DenseMatrix64F P, DenseMatrix64F QP, int jbo, int jmo, int jdo) {
    CommonOps.scale(1.0, P, QP);
}
 

static void getGradientQInvForBranch(DenseMatrix64F Qi, DenseMatrix64F Wi,
                                     DenseMatrix64F Vi, DenseMatrix64F delta,
                                     DenseMatrix64F grad) {

    CommonOps.scale(0.5, Wi, grad);

    CommonOps.multAddTransB(-0.5, delta, delta, grad);

    CommonOps.addEquals(grad, -0.5, Vi);

    if (DEBUG) {
        System.err.println("\tgradientQi = " + NormalSufficientStatistics.toVectorizedString(grad));
    }

    MultivariateChainRule ruleI = new MultivariateChainRule.InverseGeneral(Qi);
    ruleI.chainGradient(grad);

    if (DEBUG) {
        System.err.println("\tgradientQiInv = " + NormalSufficientStatistics.toVectorizedString(grad));
    }

}
 

private void chainRuleActualizationWrtAttenuationDiagonal(double ti, DenseMatrix64F grad) {
//        MissingOps.diagMult(actualization, grad);
        CommonOps.scale(-ti, grad);
    }
 

private DenseMatrix64F getGradientVarianceWrtActualizationDiagonal(ContinuousDiffusionIntegrator cdi, BranchSufficientStatistics statistics,
                                                                       int nodeIndex, DenseMatrix64F gradient) {
        // q_i
//        double[] qi = new double[dim];
//        cdi.getBranchActualization(getMatrixBufferOffsetIndex(nodeIndex), qi);
//        DenseMatrix64F qi = wrapDiagonal(actualization, 0, dim);
        // q_i^-1
//        DenseMatrix64F qiInv = wrapDiagonalInverse(actualization, 0, dim);
        // Q_i^-
        DenseMatrix64F Wi = statistics.getAbove().getRawVarianceCopy();
        // Gamma
//        DenseMatrix64F Gamma = wrap(
//                ((SafeMultivariateDiagonalActualizedWithDriftIntegrator) cdi).getStationaryVariance(getEigenBufferOffsetIndex(0)),
//                0, dim, dim);
        // Branch variance
        double[] branchVariance = new double[dim * dim];
        cdi.getBranchVariance(getMatrixBufferOffsetIndex(nodeIndex), getEigenBufferOffsetIndex(0) , branchVariance);
        DenseMatrix64F Sigma_i = wrap(branchVariance, 0, dim, dim);

        // factor
        DenseMatrix64F res = new DenseMatrix64F(dim, dim);
        CommonOps.addEquals(Wi, -1, Sigma_i);
//        MissingOps.diagDiv(qi, Wi);
//        CommonOps.multTransA(qiInv, tmp, res);
        CommonOps.multTransB(Wi, gradient, res);

        // temp + temp^T
//        MissingOps.addTransEquals(res); No need for diagonal case
        CommonOps.scale(2.0, res);

        // Get diagonal
        DenseMatrix64F resDiag = new DenseMatrix64F(dim, 1);
        CommonOps.extractDiag(res, resDiag);

        // Wrt attenuation
        double ti = cdi.getBranchLength(getMatrixBufferOffsetIndex(nodeIndex));
        chainRuleActualizationWrtAttenuationDiagonal(ti, resDiag);

        return resDiag;

    }
 

@Override
            public double getGradientForBranch(BranchSufficientStatistics statistics, double differentialScaling) {

                final NormalSufficientStatistics below = statistics.getBelow();
                final NormalSufficientStatistics branch = statistics.getBranch();
                final NormalSufficientStatistics above = statistics.getAbove();

                if (DEBUG) {
                    System.err.println("B = " + statistics.toVectorizedString());
                }

//                if (DEBUG) {
//                    System.err.println("\tQi = " + NormalSufficientStatistics.toVectorizedString(Qi));
//                    System.err.println("\tSi = " + NormalSufficientStatistics.toVectorizedString(Si));
//                }

                DenseMatrix64F Qi = above.getRawPrecision();

                DenseMatrix64F logDetComponent = matrix0;
                DenseMatrix64F quadraticComponent = matrix1;
                makeGradientMatrices0(quadraticComponent, logDetComponent, statistics, differentialScaling);

                double grad1 = 0.0;
                for (int row = 0; row < dim; ++row) {
                    grad1 -= 0.5 * logDetComponent.unsafe_get(row, row);
                }

                if (DEBUG) {
                    System.err.println("grad1 = " + grad1);
                }


//                if (DEBUG) {
//                    System.err.println("\tQi = " + NormalSufficientStatistics.toVectorizedString(Qi));
//                    System.err.println("\tQQ = " + NormalSufficientStatistics.toVectorizedString(quadraticComponent));
//                }

                NormalSufficientStatistics jointStatistics = computeJointStatistics(below, above, dim);

                DenseMatrix64F additionalVariance = matrix0; //new DenseMatrix64F(dim, dim);
                makeGradientMatrices1(additionalVariance, quadraticComponent,
                        jointStatistics);


                DenseMatrix64F delta = vector0;
                makeDeltaVector(delta, jointStatistics, above);

                double grad2 = 0.0;
                for (int row = 0; row < dim; ++row) {
                    for (int col = 0; col < dim; ++col) {

                        grad2 += 0.5 * delta.unsafe_get(row, 0)
                                * quadraticComponent.unsafe_get(row, col)
                                * delta.unsafe_get(col, 0);

                    }

                    grad2 += 0.5 * additionalVariance.unsafe_get(row, row);
                }

                if (DEBUG) {
                    System.err.println("\tjoint mean  = " + NormalSufficientStatistics.toVectorizedString(jointStatistics.getRawMean()));
                    System.err.println("\tabove mean = " + NormalSufficientStatistics.toVectorizedString(above.getRawMean()));
                    System.err.println("\tbelow mean  = " + NormalSufficientStatistics.toVectorizedString(below.getRawMean()));
                    System.err.println("\tjoint precision  = " + NormalSufficientStatistics.toVectorizedString(jointStatistics.getRawPrecision()));
                    System.err.println("\tabove precision = " + NormalSufficientStatistics.toVectorizedString(above.getRawPrecision()));
                    System.err.println("\tbelow precision  = " + NormalSufficientStatistics.toVectorizedString(below.getRawPrecision()));
                    System.err.println("\tbelow variance   = " + NormalSufficientStatistics.toVectorizedString(below.getRawVariance()));
                    System.err.println("\tquadratic      = " + NormalSufficientStatistics.toVectorizedString(quadraticComponent));
                    System.err.println("\tadditional     = " + NormalSufficientStatistics.toVectorizedString(additionalVariance));
                    System.err.println("delta: " + NormalSufficientStatistics.toVectorizedString(delta));
                    System.err.println("grad2 = " + grad2);
                }

                // W.r.t. drift
                // TODO: Fix delegate to (possibly) un-link drift from arbitrary rate
                DenseMatrix64F Di = new DenseMatrix64F(dim, 1);
                CommonOps.scale(differentialScaling, branch.getRawMean(), Di);

                double grad3 = 0.0;
                for (int row = 0; row < dim; ++row) {
                    for (int col = 0; col < dim; ++col) {

                        grad3 += delta.unsafe_get(row, 0)
                                * Qi.unsafe_get(row, col)
                                * Di.unsafe_get(col, 0);

                    }
                }

                if (DEBUG) {
                    System.err.println("\tDi     = " + NormalSufficientStatistics.toVectorizedString(branch.getRawMean()));
                    System.err.println("grad3 = " + grad3);
                }

                return grad1 + grad2 + grad3;

            }
 

private void computeVariance(DenseMatrix64F matrix, double[] data, int numRows, int numCols) {

        double[] buffer = new double[numRows];
        DenseMatrix64F sumVec = new DenseMatrix64F(numCols, 1);
        DenseMatrix64F matrixBuffer = new DenseMatrix64F(numCols, numCols);

        Arrays.fill(matrix.getData(), 0);

        for (int i = 0; i < numRows; i++) {
            int offset = numCols * i;

            DenseMatrix64F wrapper = MissingOps.wrap(data, offset, numCols, 1, buffer);

            CommonOps.multTransB(wrapper, wrapper, matrixBuffer);
            CommonOps.addEquals(matrix, matrixBuffer);
            CommonOps.addEquals(sumVec, wrapper);

        }

        CommonOps.multTransB(sumVec, sumVec, matrixBuffer);
        CommonOps.addEquals(matrix, -1.0 / numRows, matrixBuffer);
        CommonOps.scale(1.0 / numRows, matrix);


    }
 

public static void forceSymmetric(DenseMatrix64F P) {
    DenseMatrix64F Ptrans = new DenseMatrix64F(P);
    CommonOps.transpose(P, Ptrans);
    CommonOps.addEquals(P, Ptrans);
    CommonOps.scale(0.5, P);
}