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

Example 6 with DoubleMatrix1D

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

the class MGM method logsumexp.

// avoid underflow in log(sum(exp(x))) calculation
private double logsumexp(DoubleMatrix1D x) {
    DoubleMatrix1D myX = x.copy();
    double maxX = StatUtils.max(myX.toArray());
    return Math.log(myX.assign(Functions.minus(maxX)).assign(Functions.exp).zSum()) + maxX;
}

Example 7 with DoubleMatrix1D

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

the class MGM method margSum.

/*
     * PRIVATE UTILS
     */
// Utils
// sum rows together if marg == 1 and cols together if marg == 2
// Using row-major speeds up marg=1 5x
private static DoubleMatrix1D margSum(DoubleMatrix2D mat, int marg) {
    int n = 0;
    DoubleMatrix1D vec = null;
    DoubleFactory1D fac = DoubleFactory1D.dense;
    if (marg == 1) {
        n = mat.columns();
        vec = fac.make(n);
        for (int j = 0; j < mat.rows(); j++) {
            if (Thread.currentThread().isInterrupted()) {
                break;
            }
            for (int i = 0; i < n; i++) {
                vec.setQuick(i, vec.getQuick(i) + mat.getQuick(j, i));
            }
        }
    } else if (marg == 2) {
        n = mat.rows();
        vec = fac.make(n);
        for (int i = 0; i < n; i++) {
            if (Thread.currentThread().isInterrupted()) {
                break;
            }
            vec.setQuick(i, mat.viewRow(i).zSum());
        }
    }
    return vec;
}

Example 8 with DoubleMatrix1D

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

the class MGM method initParameters.

// init all parameters to zeros except for betad which is set to 1s
private void initParameters() {
    lcumsum = new int[l.length + 1];
    lcumsum[0] = 0;
    for (int i = 0; i < l.length; i++) {
        lcumsum[i + 1] = lcumsum[i] + l[i];
    }
    lsum = lcumsum[l.length];
    // LH init to zeros, maybe should be random init?
    // continuous-continuous
    DoubleMatrix2D beta = factory2D.make(xDat.columns(), xDat.columns());
    // cont squared node pot
    DoubleMatrix1D betad = factory1D.make(xDat.columns(), 1.0);
    // continuous-discrete
    DoubleMatrix2D theta = factory2D.make(lsum, xDat.columns());
    // continuous-discrete
    ;
    // discrete-discrete
    DoubleMatrix2D phi = factory2D.make(lsum, lsum);
    // cont linear node pot
    DoubleMatrix1D alpha1 = factory1D.make(xDat.columns());
    // disc node potbeta =
    DoubleMatrix1D alpha2 = factory1D.make(lsum);
    params = new MGMParams(beta, betad, theta, phi, alpha1, alpha2);
// separate lambda for each type of edge, [cc, cd, dd]
// lambda = factory1D.make(3);
}

Example 9 with DoubleMatrix1D

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

the class MGM method smooth.

/**
 * non-penalized -log(pseudolikelihood) this is the smooth function g(x) in prox gradient
 * this overloaded version calculates both nll and the smooth gradient at the same time
 * any value in gradOut will be replaced by the new calculations
 *
 * @param parIn
 * @param gradOutVec
 * @return
 */
public double smooth(DoubleMatrix1D parIn, DoubleMatrix1D gradOutVec) {
    // work with copy
    MGMParams par = new MGMParams(parIn, p, lsum);
    MGMParams gradOut = new MGMParams();
    for (int i = 0; i < par.betad.size(); i++) {
        if (par.betad.get(i) < 0)
            return Double.POSITIVE_INFINITY;
    }
    // beta=beta-diag(diag(beta));
    // for r=1:q
    // phi(Lsum(r)+1:Lsum(r+1),Lsum(r)+1:Lsum(r+1))=0;
    // end
    // beta=triu(beta); phi=triu(phi);
    // 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 matrix is 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
    // tempLoss =  (X-e*alpha1'-Xbeta-Dtheta) = -res (in gradient code)
    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++) {
        if (Thread.currentThread().isInterrupted()) {
            break;
        }
        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));
        }
    }
    // sqloss=-n/2*sum(log(betad))+...
    // .5*norm((X-e*alpha1'-Xbeta-Dtheta)*diag(sqrt(betad)),'fro')^2;
    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);
    // ok now tempLoss = res
    tempLoss.assign(Functions.mult(-1));
    // gradbeta=X'*(res);
    gradOut.beta = alg.mult(alg.transpose(xDat), tempLoss);
    // gradbeta=gradbeta-diag(diag(gradbeta)); % zero out diag
    // gradbeta=tril(gradbeta)'+triu(gradbeta);
    DoubleMatrix2D lowerBeta = alg.transpose(lowerTri(gradOut.beta.copy(), -1));
    upperTri(gradOut.beta, 1).assign(lowerBeta, Functions.plus);
    // gradalpha1=diag(betad)*sum(res,1)';
    gradOut.alpha1 = alg.mult(factory2D.diagonal(par.betad), margSum(tempLoss, 1));
    // gradtheta=D'*(res);
    gradOut.theta = alg.mult(alg.transpose(dDat), tempLoss);
    // 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++) {
        if (Thread.currentThread().isInterrupted()) {
            break;
        }
        DoubleMatrix2D wxTemp = wxProd.viewPart(0, lcumsum[i], n, l[i]);
        // need to copy init values for calculating nll
        DoubleMatrix2D wxTemp0 = wxTemp.copy();
        // 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++) {
            if (Thread.currentThread().isInterrupted()) {
                break;
            }
            DoubleMatrix1D curRow = wxTemp.viewRow(k);
            DoubleMatrix1D curRow0 = wxTemp0.viewRow(k);
            catloss -= curRow0.get((int) yDat.get(k, i) - 1);
            catloss += logsumexp(curRow0);
            // 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)';
    gradOut.alpha2 = margSum(wxProd, 1);
    // gradw=X'*wxprod;
    DoubleMatrix2D gradW = alg.mult(alg.transpose(xDat), wxProd);
    // gradtheta=gradtheta+gradw';
    gradOut.theta.assign(alg.transpose(gradW), Functions.plus);
    // gradphi=D'*wxprod;
    gradOut.phi = alg.mult(alg.transpose(dDat), wxProd);
    // end
    for (int i = 0; i < q; i++) {
        gradOut.phi.viewPart(lcumsum[i], lcumsum[i], l[i], l[i]).assign(0);
    }
    // gradphi=tril(gradphi)'+triu(gradphi);
    DoubleMatrix2D lowerPhi = alg.transpose(lowerTri(gradOut.phi.copy(), 0));
    upperTri(gradOut.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
         */
    gradOut.betad = factory1D.make(xDat.columns());
    for (int i = 0; i < p; i++) {
        gradOut.betad.set(i, -n / (2.0 * par.betad.get(i)) + alg.norm2(tempLoss.viewColumn(i)) / 2.0 - alg.mult(tempLoss.viewColumn(i), xBeta.viewColumn(i).copy().assign(dTheta.viewColumn(i), Functions.plus)));
    }
    gradOut.alpha1.assign(Functions.div((double) n));
    gradOut.alpha2.assign(Functions.div((double) n));
    gradOut.betad.assign(Functions.div((double) n));
    gradOut.beta.assign(Functions.div((double) n));
    gradOut.theta.assign(Functions.div((double) n));
    gradOut.phi.assign(Functions.div((double) n));
    gradOutVec.assign(gradOut.toMatrix1D());
    return (sqloss + catloss) / ((double) n);
}

Example 10 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project tdq-studio-se by Talend.

the class WrapperDoubleMatrix1D method viewSelection.

/**
 *Constructs and returns a new <i>selection view</i> that is a matrix holding the indicated cells.
 *There holds <tt>view.size() == indexes.length</tt> and <tt>view.get(i) == this.get(indexes[i])</tt>.
 *Indexes can occur multiple times and can be in arbitrary order.
 *<p>
 *<b>Example:</b>
 *<br>
 *<pre>
 *this     = (0,0,8,0,7)
 *indexes  = (0,2,4,2)
 *-->
 *view     = (0,8,7,8)
 *</pre>
 *Note that modifying <tt>indexes</tt> after this call has returned has no effect on the view.
 *The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
 *
 *@param  indexes   The indexes of the cells that shall be visible in the new view. To indicate that <i>all</i> cells shall be visible, simply set this parameter to <tt>null</tt>.
 *@return the new view.
 *@throws IndexOutOfBoundsException if <tt>!(0 <= indexes[i] < size())</tt> for any <tt>i=0..indexes.length()-1</tt>.
 */
public DoubleMatrix1D viewSelection(int[] indexes) {
    // check for "all"
    if (indexes == null) {
        indexes = new int[size];
        for (int i = size; --i >= 0; ) indexes[i] = i;
    }
    checkIndexes(indexes);
    final int[] idx = indexes;
    DoubleMatrix1D view = new WrapperDoubleMatrix1D(this) {

        public double getQuick(int i) {
            return content.get(idx[i]);
        }

        public void setQuick(int i, double value) {
            content.set(idx[i], value);
        }
    };
    view.size = indexes.length;
    return view;
}