Examples with DoubleMatrix1D cern.colt.matrix.DoubleMatrix1D used on opensource projects

Example 66 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project tetrad by cmu-phil.

the class MGM method smoothValue.

/**
 * non-penalized -log(pseudolikelihood) this is the smooth function g(x) in prox gradient
 *
 * @param parIn
 * @return
 */
public double smoothValue(DoubleMatrix1D parIn) {
    // work with copy
    MGMParams par = new MGMParams(parIn, p, lsum);
    for (int i = 0; i < par.betad.size(); i++) {
        if (par.betad.get(i) < 0)
            return Double.POSITIVE_INFINITY;
    }
    // double nll = 0;
    // int n = xDat.rows();
    // beta=beta+beta';
    // phi=phi+phi';
    upperTri(par.beta, 1);
    par.beta.assign(alg.transpose(par.beta), Functions.plus);
    for (int i = 0; i < q; i++) {
        par.phi.viewPart(lcumsum[i], lcumsum[i], l[i], l[i]).assign(0);
    }
    // ensure mats are upper triangular
    upperTri(par.phi, 0);
    par.phi.assign(alg.transpose(par.phi), Functions.plus);
    // Xbeta=X*beta*diag(1./betad);
    DoubleMatrix2D divBetaD = factory2D.diagonal(factory1D.make(p, 1.0).assign(par.betad, Functions.div));
    DoubleMatrix2D xBeta = alg.mult(xDat, alg.mult(par.beta, divBetaD));
    // Dtheta=D*theta*diag(1./betad);
    DoubleMatrix2D dTheta = alg.mult(alg.mult(dDat, par.theta), divBetaD);
    // Squared loss
    // sqloss=-n/2*sum(log(betad))+...
    // .5*norm((X-e*alpha1'-Xbeta-Dtheta)*diag(sqrt(betad)),'fro')^2;
    DoubleMatrix2D tempLoss = factory2D.make(n, xDat.columns());
    // wxprod=X*(theta')+D*phi+e*alpha2';
    DoubleMatrix2D wxProd = alg.mult(xDat, alg.transpose(par.theta));
    wxProd.assign(alg.mult(dDat, par.phi), Functions.plus);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < xDat.columns(); j++) {
            tempLoss.set(i, j, xDat.get(i, j) - par.alpha1.get(j) - xBeta.get(i, j) - dTheta.get(i, j));
        }
        for (int j = 0; j < dDat.columns(); j++) {
            wxProd.set(i, j, wxProd.get(i, j) + par.alpha2.get(j));
        }
    }
    double sqloss = -n / 2.0 * par.betad.copy().assign(Functions.log).zSum() + .5 * Math.pow(alg.normF(alg.mult(tempLoss, factory2D.diagonal(par.betad.copy().assign(Functions.sqrt)))), 2);
    // categorical loss
    /*catloss=0;
        wxprod=X*(theta')+D*phi+e*alpha2'; %this is n by Ltot
        for r=1:q
            wxtemp=wxprod(:,Lsum(r)+1:Lsum(r)+L(r));
            denom= logsumexp(wxtemp,2); %this is n by 1
            catloss=catloss-sum(wxtemp(sub2ind([n L(r)],(1:n)',Y(:,r))));
            catloss=catloss+sum(denom);
        end
        */
    double catloss = 0;
    for (int i = 0; i < yDat.columns(); i++) {
        DoubleMatrix2D wxTemp = wxProd.viewPart(0, lcumsum[i], n, l[i]);
        for (int k = 0; k < n; k++) {
            DoubleMatrix1D curRow = wxTemp.viewRow(k);
            catloss -= curRow.get((int) yDat.get(k, i) - 1);
            catloss += logsumexp(curRow);
        }
    }
    return (sqloss + catloss) / ((double) n);
}

Example 67 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project tetrad by cmu-phil.

the class MGM method smoothGradient.

/**
 * Gradient of the pseudolikelihood
 *
 * @param parIn
 * @return
 */
public DoubleMatrix1D smoothGradient(DoubleMatrix1D parIn) {
    int n = xDat.rows();
    MGMParams grad = new MGMParams();
    // 
    MGMParams par = new MGMParams(parIn, p, lsum);
    upperTri(par.beta, 1);
    par.beta.assign(alg.transpose(par.beta), Functions.plus);
    for (int i = 0; i < q; i++) {
        par.phi.viewPart(lcumsum[i], lcumsum[i], l[i], l[i]).assign(0);
    }
    upperTri(par.phi, 0);
    par.phi.assign(alg.transpose(par.phi), Functions.plus);
    // Xbeta=X*beta*diag(1./betad);
    // Dtheta=D*theta*diag(1./betad);
    DoubleMatrix2D divBetaD = factory2D.diagonal(factory1D.make(p, 1.0).assign(par.betad, Functions.div));
    DoubleMatrix2D xBeta = alg.mult(alg.mult(xDat, par.beta), divBetaD);
    DoubleMatrix2D dTheta = alg.mult(alg.mult(dDat, par.theta), divBetaD);
    // res=Xbeta-X+e*alpha1'+Dtheta;
    DoubleMatrix2D negLoss = factory2D.make(n, xDat.columns());
    // wxprod=X*(theta')+D*phi+e*alpha2';
    DoubleMatrix2D wxProd = alg.mult(xDat, alg.transpose(par.theta));
    wxProd.assign(alg.mult(dDat, par.phi), Functions.plus);
    for (int i = 0; i < n; i++) {
        if (Thread.currentThread().isInterrupted()) {
            break;
        }
        for (int j = 0; j < p; j++) {
            if (Thread.currentThread().isInterrupted()) {
                break;
            }
            negLoss.set(i, j, xBeta.get(i, j) - xDat.get(i, j) + par.alpha1.get(j) + dTheta.get(i, j));
        }
        for (int j = 0; j < dDat.columns(); j++) {
            if (Thread.currentThread().isInterrupted()) {
                break;
            }
            wxProd.set(i, j, wxProd.get(i, j) + par.alpha2.get(j));
        }
    }
    // gradbeta=X'*(res);
    grad.beta = alg.mult(alg.transpose(xDat), negLoss);
    // gradbeta=gradbeta-diag(diag(gradbeta)); % zero out diag
    // gradbeta=tril(gradbeta)'+triu(gradbeta);
    DoubleMatrix2D lowerBeta = alg.transpose(lowerTri(grad.beta.copy(), -1));
    upperTri(grad.beta, 1).assign(lowerBeta, Functions.plus);
    // gradalpha1=diag(betad)*sum(res,1)';
    grad.alpha1 = alg.mult(factory2D.diagonal(par.betad), margSum(negLoss, 1));
    // gradtheta=D'*(res);
    grad.theta = alg.mult(alg.transpose(dDat), negLoss);
    for (int i = 0; i < yDat.columns(); i++) {
        DoubleMatrix2D wxTemp = wxProd.viewPart(0, lcumsum[i], n, l[i]);
        // does this need to be done in log space??
        wxTemp.assign(Functions.exp);
        DoubleMatrix1D invDenom = factory1D.make(n, 1.0).assign(margSum(wxTemp, 2), Functions.div);
        wxTemp.assign(alg.mult(factory2D.diagonal(invDenom), wxTemp));
        for (int k = 0; k < n; k++) {
            DoubleMatrix1D curRow = wxTemp.viewRow(k);
            // wxtemp(sub2ind(size(wxtemp),(1:n)',Y(:,r)))=wxtemp(sub2ind(size(wxtemp),(1:n)',Y(:,r)))-1;
            curRow.set((int) yDat.get(k, i) - 1, curRow.get((int) yDat.get(k, i) - 1) - 1);
        }
    }
    // gradalpha2=sum(wxprod,1)';
    grad.alpha2 = margSum(wxProd, 1);
    // gradw=X'*wxprod;
    DoubleMatrix2D gradW = alg.mult(alg.transpose(xDat), wxProd);
    // gradtheta=gradtheta+gradw';
    grad.theta.assign(alg.transpose(gradW), Functions.plus);
    // gradphi=D'*wxprod;
    grad.phi = alg.mult(alg.transpose(dDat), wxProd);
    // end
    for (int i = 0; i < q; i++) {
        if (Thread.currentThread().isInterrupted()) {
            break;
        }
        grad.phi.viewPart(lcumsum[i], lcumsum[i], l[i], l[i]).assign(0);
    }
    // gradphi=tril(gradphi)'+triu(gradphi);
    DoubleMatrix2D lowerPhi = alg.transpose(lowerTri(grad.phi.copy(), 0));
    upperTri(grad.phi, 0).assign(lowerPhi, Functions.plus);
    /*
        for s=1:p
            gradbetad(s)=-n/(2*betad(s))+1/2*norm(res(:,s))^2-res(:,s)'*(Xbeta(:,s)+Dtheta(:,s));
        end
         */
    grad.betad = factory1D.make(xDat.columns());
    for (int i = 0; i < p; i++) {
        if (Thread.currentThread().isInterrupted()) {
            break;
        }
        grad.betad.set(i, -n / (2.0 * par.betad.get(i)) + alg.norm2(negLoss.viewColumn(i)) / 2.0 - alg.mult(negLoss.viewColumn(i), xBeta.viewColumn(i).copy().assign(dTheta.viewColumn(i), Functions.plus)));
    }
    grad.alpha1.assign(Functions.div((double) n));
    grad.alpha2.assign(Functions.div((double) n));
    grad.betad.assign(Functions.div((double) n));
    grad.beta.assign(Functions.div((double) n));
    grad.theta.assign(Functions.div((double) n));
    grad.phi.assign(Functions.div((double) n));
    return grad.toMatrix1D();
}

Example 68 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project tetrad by cmu-phil.

the class MGM method proximalOperator.

/**
 * A proximal operator for the MGM
 *
 * @param t parameter for operator, must be positive
 * @param X input vector to operator
 * @return output vector, same dimension as X
 */
public DoubleMatrix1D proximalOperator(double t, DoubleMatrix1D X) {
    // System.out.println("PROX with t = " + t);
    if (t <= 0)
        throw new IllegalArgumentException("t must be positive: " + t);
    DoubleMatrix1D tlam = lambda.copy().assign(Functions.mult(t));
    // Constructor copies and checks dimension
    // par is a copy so we can update it
    MGMParams par = new MGMParams(X.copy(), p, lsum);
    // penbeta = t(1).*(wv(1:p)'*wv(1:p));
    // betascale=zeros(size(beta));
    // betascale=max(0,1-penbeta./abs(beta));
    DoubleMatrix2D weightMat = alg.multOuter(weights, weights, null);
    DoubleMatrix2D betaWeight = weightMat.viewPart(0, 0, p, p);
    DoubleMatrix2D betascale = betaWeight.copy().assign(Functions.mult(-tlam.get(0)));
    betascale.assign(par.beta.copy().assign(Functions.abs), Functions.div);
    betascale.assign(Functions.plus(1));
    betascale.assign(Functions.max(0));
    // par.beta.assign(betascale, Functions.mult);
    for (int i = 0; i < p; i++) {
        if (Thread.currentThread().isInterrupted()) {
            break;
        }
        for (int j = 0; j < p; j++) {
            if (Thread.currentThread().isInterrupted()) {
                break;
            }
            double curVal = par.beta.get(i, j);
            if (curVal != 0) {
                par.beta.set(i, j, curVal * betascale.get(i, j));
            }
        }
    }
    /*
        thetanorms=0;
        for s=1:p
            for j=1:q
                tempvec=theta(Lsums(j)+1:Lsums(j+1),s);
                tempvec=max(0,1-t(2)*(wv(s)*wv(p+j))/norm(tempvec))*tempvec;
                thetanorms=thetanorms+(wv(s)*wv(p+j))*norm(tempvec);
                theta(Lsums(j)+1:Lsums(j+1),s)=tempvec(1:L(j));
            end
        end
        */
    for (int i = 0; i < p; i++) {
        if (Thread.currentThread().isInterrupted()) {
            break;
        }
        for (int j = 0; j < lcumsum.length - 1; j++) {
            if (Thread.currentThread().isInterrupted()) {
                break;
            }
            DoubleMatrix1D tempVec = par.theta.viewColumn(i).viewPart(lcumsum[j], l[j]);
            // double thetaScale = Math.max(0, 1 - tlam.get(1)*weightMat.get(i, p+j)/Math.sqrt(alg.norm2(tempVec)));
            double foo = norm2(tempVec);
            double thetaScale = Math.max(0, 1 - tlam.get(1) * weightMat.get(i, p + j) / norm2(tempVec));
            tempVec.assign(Functions.mult(thetaScale));
        }
    }
    /*
        for r=1:q
            for j=1:q
                if r<j
                    tempmat=phi(Lsums(r)+1:Lsums(r+1),Lsums(j)+1:Lsums(j+1));
                    tempmat=max(0,1-t(3)*(wv(p+r)*wv(p+j))/norm(tempmat))*tempmat; % Lj by 2*Lr
                    phinorms=phinorms+(wv(p+r)*wv(p+j))*norm(tempmat,'fro');
                    phi( Lsums(r)+1:Lsums(r+1),Lsums(j)+1:Lsums(j+1) )=tempmat;
                end
            end
        end
         */
    for (int i = 0; i < lcumsum.length - 1; i++) {
        if (Thread.currentThread().isInterrupted()) {
            break;
        }
        for (int j = i + 1; j < lcumsum.length - 1; j++) {
            if (Thread.currentThread().isInterrupted()) {
                break;
            }
            DoubleMatrix2D tempMat = par.phi.viewPart(lcumsum[i], lcumsum[j], l[i], l[j]);
            // Not sure why this isnt Frobenius norm...
            // double phiScale = Math.max(0, 1-tlam.get(2)*weightMat.get(p+i,p+j)/alg.norm2(tempMat));
            double phiScale = Math.max(0, 1 - tlam.get(2) * weightMat.get(p + i, p + j) / norm2(tempMat));
            // double phiScale = Math.max(0, 1-tlam.get(2)*weightMat.get(p+i,p+j)/alg.normF(tempMat));
            tempMat.assign(Functions.mult(phiScale));
        }
    }
    return par.toMatrix1D();
}

Example 69 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project tetrad by cmu-phil.

the class ProximalGradient method learnBackTrack.

// run FISTA with step size backtracking attempt to speed up
public DoubleMatrix1D learnBackTrack(ConvexProximal cp, DoubleMatrix1D Xin, double epsilon, int iterLimit) {
    DoubleMatrix1D X = cp.proximalOperator(1.0, Xin.copy());
    DoubleMatrix1D Y = X.copy();
    DoubleMatrix1D Z = X.copy();
    DoubleMatrix1D GrY = cp.smoothGradient(Y);
    DoubleMatrix1D GrX = cp.smoothGradient(X);
    int iterCount = 0;
    int noEdgeChangeCount = 0;
    double theta = Double.POSITIVE_INFINITY;
    double thetaOld = theta;
    double L = 1.0;
    double Lold = L;
    boolean backtrackSwitch = true;
    double dx;
    double Fx = Double.POSITIVE_INFINITY;
    double Gx = Double.POSITIVE_INFINITY;
    double Fy;
    double obj;
    while (true) {
        Lold = L;
        L = L * alpha;
        thetaOld = theta;
        DoubleMatrix1D Xold = X.copy();
        obj = Fx + Gx;
        while (true) {
            theta = 2.0 / (1.0 + Math.sqrt(1.0 + (4.0 * L) / (Lold * Math.pow(thetaOld, 2))));
            if (theta < 1) {
                Y.assign(Xold.copy().assign(Functions.mult(1 - theta)));
                Y.assign(Z.copy().assign(Functions.mult(theta)), Functions.plus);
            }
            Fy = cp.smooth(Y, GrY);
            DoubleMatrix1D temp = Y.copy().assign(GrY.copy().assign(Functions.mult(1.0 / L)), Functions.minus);
            Gx = cp.nonSmooth(1.0 / L, temp, X);
            if (backtrackSwitch) {
                Fx = cp.smoothValue(X);
            } else {
                // tempPar = new MGMParams();
                Fx = cp.smooth(X, GrX);
            // GrX.assign(factory1D.make(tempPar.toVector()[0]));
            }
            DoubleMatrix1D XmY = X.copy().assign(Y, Functions.minus);
            double normXY = alg.norm2(XmY);
            if (normXY == 0)
                break;
            double Qx;
            double LocalL;
            if (backtrackSwitch) {
                // System.out.println("Back Norm");
                Qx = Fy + alg.mult(XmY, GrY) + (L / 2.0) * normXY;
                LocalL = L + 2 * Math.max(Fx - Qx, 0) / normXY;
                backtrackSwitch = Math.abs(Fy - Fx) >= backtrackTol * Math.max(Math.abs(Fx), Math.abs(Fy));
            } else {
                // System.out.println("Close Rule");
                // it shouldn't be possible for GrX to be null here...
                // if(GrX==null)
                // GrX = factory1D.make(gradient(Xpar).toVector()[0]);
                // Fx = alg.mult(YmX, Gx.assign(G, Functions.minus));
                // Qx = (L / 2.0) * alg.norm2(YmX);
                LocalL = 2 * alg.mult(XmY, GrX.assign(GrY, Functions.minus)) / normXY;
            }
            // System.out.println("LocalL: " + LocalL + " L: " + L);
            if (LocalL <= L) {
                break;
            } else if (LocalL != Double.POSITIVE_INFINITY) {
                L = LocalL;
            } else {
                LocalL = L;
            }
            L = Math.max(LocalL, L / beta);
        }
        int diffEdges = 0;
        for (int i = 0; i < X.size(); i++) {
            double a = X.get(i);
            double b = Xold.get(i);
            if (a != 0 & b == 0) {
                diffEdges++;
            } else if (a == 0 & b != 0) {
                diffEdges++;
            }
        }
        dx = norm2(X.copy().assign(Xold, Functions.minus)) / Math.max(1, norm2(X));
        // sometimes there are more edge changes after initial 0, so may want to do two zeros in a row...
        if (diffEdges == 0 && edgeConverge) {
            noEdgeChangeCount++;
            if (noEdgeChangeCount >= noEdgeChangeTol) {
                System.out.println("Edges converged at iter: " + iterCount + " with |dx|/|x|: " + dx);
                System.out.println("Iter: " + iterCount + " |dx|/|x|: " + dx + " normX: " + norm2(X) + " nll: " + Fx + " reg: " + Gx + " DiffEdges: " + diffEdges + " L: " + L);
                break;
            }
        // negative noEdgeChangeTol stops when diffEdges <= |noEdgeChangeTol|
        } else if (noEdgeChangeTol < 0 && diffEdges <= Math.abs(noEdgeChangeTol)) {
            System.out.println("Edges converged at iter: " + iterCount + " with |dx|/|x|: " + dx);
            System.out.println("Iter: " + iterCount + " |dx|/|x|: " + dx + " normX: " + norm2(X) + " nll: " + Fx + " reg: " + Gx + " DiffEdges: " + diffEdges + " L: " + L);
            break;
        } else {
            noEdgeChangeCount = 0;
        }
        // edge converge should happen before params converge, unless epsilon is big
        if (dx < epsilon && !edgeConverge) {
            System.out.println("Converged at iter: " + iterCount + " with |dx|/|x|: " + dx + " < epsilon: " + epsilon);
            System.out.println("Iter: " + iterCount + " |dx|/|x|: " + dx + " normX: " + norm2(X) + " nll: " + Fx + " reg: " + Gx + " DiffEdges: " + diffEdges + " L: " + L);
            break;
        }
        // restart acceleration if objective got worse
        if (Fx + Gx > obj) {
            theta = Double.POSITIVE_INFINITY;
            Y.assign(X.copy());
            // Ypar = new MGMParams(Xpar);
            Z.assign(X.copy());
        // Fy = Fx;
        // GrY.assign(GrX.copy());
        } else if (theta == 1) {
            Z.assign(X.copy());
        } else {
            Z.assign(X.copy().assign(Functions.mult(1 / theta)));
            Z.assign(Xold.copy().assign(Functions.mult(1 - (1.0 / theta))), Functions.plus);
        }
        if (iterCount % printIter == 0) {
            System.out.println("Iter: " + iterCount + " |dx|/|x|: " + dx + " normX: " + norm2(X) + " nll: " + Fx + " reg: " + Gx + " DiffEdges: " + diffEdges + " L: " + L);
        // System.out.println("Iter: " + iterCount + " |dx|/|x|: " + dx + " nll: " + negLogLikelihood(params) + " reg: " + regTerm(params));
        }
        // System.out.println("t: " + t);
        // System.out.println("Parameters: " + params);
        iterCount++;
        if (iterCount >= iterLimit) {
            System.out.println("Iter limit reached");
            break;
        }
    }
    return X;
}

Example 70 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project Gemma by PavlidisLab.

the class ComBat method plot.

/**
 * Make diagnostic plots.
 * FIXME: As in the original ComBat, this only graphs the first batch's statistics. In principle we can (and perhaps
 * should) examine these plots for all the batches.
 *
 * @param filePrefix file prefix
 */
public void plot(String filePrefix) {
    if (this.gammaHat == null)
        throw new IllegalArgumentException("You must call 'run' first");
    /*
         * View the distribution of gammaHat, which we assume will have a normal distribution
         */
    DoubleMatrix1D ghr = gammaHat.viewRow(0);
    int NUM_HIST_BINS = 100;
    Histogram gammaHatHist = new Histogram("GammaHat", NUM_HIST_BINS, ghr);
    XYSeries ghplot = gammaHatHist.plot();
    Normal rn = new Normal(this.gammaBar.get(0), Math.sqrt(this.t2.get(0)), new MersenneTwister());
    Histogram ghtheoryT = new Histogram("Gamma", NUM_HIST_BINS, gammaHatHist.min(), gammaHatHist.max());
    for (int i = 0; i < 10000; i++) {
        double n = rn.nextDouble();
        ghtheoryT.fill(n);
    }
    XYSeries ghtheory = ghtheoryT.plot();
    File tmpfile;
    try {
        tmpfile = File.createTempFile(filePrefix + ".gammahat.histogram.", ".png");
        ComBat.log.info(tmpfile);
    } catch (IOException e) {
        throw new RuntimeException(e);
    }
    try (OutputStream os = new FileOutputStream(tmpfile)) {
        this.writePlot(os, ghplot, ghtheory);
        /*
             * View the distribution of deltaHat, which we assume has an inverse gamma distribution
             */
        DoubleMatrix1D dhr = deltaHat.viewRow(0);
        Histogram deltaHatHist = new Histogram("DeltaHat", NUM_HIST_BINS, dhr);
        XYSeries dhplot = deltaHatHist.plot();
        Gamma g = new Gamma(aPrior.get(0), bPrior.get(0), new MersenneTwister());
        Histogram deltaHatT = new Histogram("Delta", NUM_HIST_BINS, deltaHatHist.min(), deltaHatHist.max());
        for (int i = 0; i < 10000; i++) {
            double invg = 1.0 / g.nextDouble();
            deltaHatT.fill(invg);
        }
        XYSeries dhtheory = deltaHatT.plot();
        tmpfile = File.createTempFile(filePrefix + ".deltahat.histogram.", ".png");
        ComBat.log.info(tmpfile);
        try (OutputStream os2 = new FileOutputStream(tmpfile)) {
            this.writePlot(os2, dhplot, dhtheory);
        }
    } catch (IOException e) {
        throw new RuntimeException(e);
    }
}