Example usage for org.apache.commons.math3.linear BlockRealMatrix operate

Introduction

In this page you can find the example usage for org.apache.commons.math3.linear BlockRealMatrix operate.

Prototype

@Override
public double[] operate(final double[] v) throws DimensionMismatchException 

Source Link

Usage

From source file:gamlss.smoothing.PB.java

/**
 * The main fitting method, initiate appropriate smoothing 
 * function according to incoming parameters.
 * @param additiveFit -object of AdditiveFoit class
 * @return reiduals/*from  w  w w. j a  v a 2  s.c o  m*/
 */
//gamlss.pb <- function(x, y, w, xeval = NULL, ...)      
public ArrayRealVector solve(final AdditiveFit additiveFit) {

    Double lambda = Controls.LAMBDAS_USER_DEFINED;
    colNum = additiveFit.getColNum();
    y = additiveFit.getZ();
    w = additiveFit.getW();
    whichDistParameter = additiveFit.getWhichDistParameter();
    basisM = (BlockRealMatrix) getBasisM().get(colNum);
    penaltyM = (BlockRealMatrix) getPenaltyM().get(colNum);
    df = (Integer) getDfValues().get(colNum);

    //n <- nrow(X) # the no of observations
    n = basisM.getRowDimension();

    //lambdaS <-  get(startLambdaName, envir=gamlss.env)
    //geting the starting value
    lambdaS = getLambdas().get(whichDistParameter);

    //if (lambdaS>=1e+07) lambda <- 1e+07
    if (lambdaS >= 1e+07) {
        lambda = 1e+07;
    }
    //if (lambdaS<=1e-07) lambda <- 1e-07
    if (lambdaS <= 1e-07) {
        lambda = 1e-07;
    }

    //if (is.null(df)&&!is.null(lambda)||!is.null(df)&&!is.null(lambda))
    if (lambda != null) {

        //fit <- regpen(y, X, w, lambda, D)
        beta = regpen(lambda);

        //fv <- X %*% fit$beta
        fv = (ArrayRealVector) basisM.operate(beta);

        //else if (is.null(df)&&is.null(lambda))
    } else if (df == null) {

        //lambda <- lambdaS
        lambda = lambdaS;

        //switch(control$c," ML"={
        switch (Controls.SMOOTH_METHOD) {
        case DistributionSettings.GAIC:
            fv = functionGAIC(lambda);
            break;
        case DistributionSettings.ML:
            fv = functionML(lambda);
            break;
        case DistributionSettings.ML1:
            fv = functionML1(lambda);
            break;
        case DistributionSettings.EM:
            System.err.println("EM has not been implemented yet");
            break;
        case DistributionSettings.GCV:
            fv = functionGCV(lambda);
            break;
        default:
            System.err.println(
                    " Cannot recognise the " + "smoothing method or it has" + "not been implemented yet");
        }
    } else {

        //QR <- qr(sqrt(w)*X)
        //Rinv <- solve(qr.R(QR))
        rM = (BlockRealMatrix) new QRDecomposition(
                MatrixFunctions.multVectorMatrix(MatrixFunctions.sqrtVec(w), basisM)).getR();
        rM = rM.getSubMatrix(0, rM.getColumnDimension() - 1, 0, rM.getColumnDimension() - 1);
        rMinv = (BlockRealMatrix) new QRDecomposition(rM).getSolver().getInverse();

        //S   <- t(D)%*%D
        sM = penaltyM.transpose().multiply(penaltyM);

        //UDU <- eigen(t(Rinv)%*%S%*%Rinv)
        uduV = new EigenDecomposition(rMinv.transpose().multiply(sM).multiply(rMinv)).getRealEigenvalues();

        //lambda <- if (sign(edf1_df(0))==sign(edf1_df(100000))) 100000
        //in case they have the some sign
        edfTemp1 = edf1_df(0, df);
        edfTemp2 = edf1_df(100000.0, df);

        if (FastMath.signum(edfTemp1) == FastMath.signum(edfTemp2)) {
            lambda = 100000.0;
        } else {

            //else uniroot(edf1_df, c(0,100000))$root
            uniRootObj.setDf(df);
            final double relativeAccuracy = 1.0e-12;
            final double absoluteAccuracy = 1.0e-8;
            UnivariateSolver uniRootSolver = new BrentSolver(relativeAccuracy, absoluteAccuracy);
            lambda = uniRootSolver.solve(1000, uniRootObj, 0.0, 100000.0);
        }

        //fit <- regpen(y, X, w, lambda, D)
        beta = regpen(lambda);
        fv = (ArrayRealVector) basisM.operate(beta);
    }
    if (Controls.IF_NEED_THIS) {

        //waug <- as.vector(c(w, rep(1,nrow(D))))
        waug = w.append(MatrixFunctions.repV(1.0, penaltyM.getRowDimension()));

        //xaug <- as.matrix(rbind(X,sqrt(lambda)*D))
        xaug = MatrixFunctions.appendMatricesRows(basisM,
                (BlockRealMatrix) penaltyM.scalarMultiply(FastMath.sqrt(lambda)));

        //lev <- hat(sqrt(waug)*xaug,intercept=FALSE)[1:n]
        //get the hat matrix
        lev = (ArrayRealVector) MatrixFunctions.getMainDiagonal(new BlockRealMatrix(MatrixFunctions
                .calculateHat(MatrixFunctions.multVectorMatrix(MatrixFunctions.sqrtVec(waug), xaug)).getData()))
                .getSubVector(0, n);

        //lev <- (lev-.hat.WX(w,x))
        //subtract  the linear since is already fitted 
        lev = lev.subtract(
                hatWX((ArrayRealVector) getSmootherMatrices().get(whichDistParameter).getColumnVector(colNum)));

        // var <- lev/w
        //the variance of the smootherz
        var = lev.ebeDivide(w);
    }

    //if (is.null(xeval)) # if no prediction
    if (Controls.XEVAL_USER_DEFINED == null) {
        // Residuals
        return y.subtract(fv);

        //else # for prediction
    } else {

        //ll <- dim(as.matrix(attr(x,"X")))[1]                 
        //nx <- as.matrix(attr(x,"X"))[seq(length(y)+1,ll),]
        tempArr = ArithmeticSeries.getSeries(y.getDimension() + 1, basisM.getRowDimension(), 1);
        BlockRealMatrix nx = new BlockRealMatrix(tempArr.length, basisM.getColumnDimension());
        for (int i = 0; i < tempArr.length; i++) {
            nx.setRowVector(i, basisM.getRowVector((int) tempArr[i]));
        }

        //pred <- drop(nx %*% fit$beta)
        ArrayRealVector pred = (ArrayRealVector) nx.operate(beta);

        // Residuals
        return y.subtract(fv);
    }
}

---