List of usage examples for org.apache.commons.math3.linear FieldLUDecomposition FieldLUDecomposition
public FieldLUDecomposition(FieldMatrix<T> matrix)
From source file:model.LP.java
/** * Find a leaving variable index that is the most * bounding on the given entering variable index. * * @param entering/* w w w .j av a 2s . com*/ * an entering variable index. * @param dual * If true, find a leaving variable index for the dual dictionary. * Otherwise, find one for the primal dictionary. * @return * A leaving variable index. */ private int leaving(int entering, boolean dual) { FieldVector<BigFraction> check; FieldVector<BigFraction> sd; FieldMatrix<BigFraction> bin = new FieldLUDecomposition<BigFraction>(B_).getSolver().getInverse() .multiply(N_); if (dual) { check = c_; FieldVector<BigFraction> unit = new ArrayFieldVector<BigFraction>(bin.getRowDimension(), BigFraction.ZERO); unit.setEntry(entering, BigFraction.ONE); sd = bin.transpose().scalarMultiply(BigFraction.MINUS_ONE).operate(unit); } else { check = b_; FieldVector<BigFraction> unit = new ArrayFieldVector<BigFraction>(bin.getColumnDimension(), BigFraction.ZERO); unit.setEntry(entering, BigFraction.ONE); sd = bin.operate(unit); } boolean unbounded = true; int index = -1; /* Check for unboundedness and find first non-zero element in check */ for (int i = 0; i < sd.getDimension(); i++) { if (!check.getEntry(i).equals(BigFraction.ZERO) && index == -1) { index = i; } if (sd.getEntry(i).compareTo(BigFraction.ZERO) > 0) { unbounded = false; } } String e = "Program is unbounded."; if (unbounded) throw new RuntimeException(e); /* Set temporarily max value as ratio of the first divisible pair. */ BigFraction max = sd.getEntry(index).divide(check.getEntry(index)); for (int i = 0; i < sd.getDimension(); i++) { BigFraction num = sd.getEntry(i); BigFraction denom = check.getEntry(i); if (!denom.equals(BigFraction.ZERO)) { BigFraction val = num.divide(denom); if (val.compareTo(max) > 0) { max = val; index = i; } } else { if (num.compareTo(BigFraction.ZERO) > 0) return i; } } return index; }
From source file:controller.VisLP.java
private static Point2D[] getFeasibleIntersections(FieldMatrix<BigFraction> cons) { FieldMatrix<BigFraction> N = cons.getSubMatrix(0, cons.getRowDimension() - 1, 0, cons.getColumnDimension() - 2); FieldVector<BigFraction> b = cons.getColumnVector(cons.getColumnDimension() - 1); HashSet<Point2D> points = new HashSet<Point2D>(); unb = new ArrayList<Point2D>(); /* Find all intersections */ for (int i = 0; i < N.getRowDimension(); i++) { for (int j = 0; j < N.getRowDimension(); j++) { if (i == j) continue; FieldMatrix<BigFraction> line1 = N.getRowMatrix(i); FieldMatrix<BigFraction> line2 = N.getRowMatrix(j); BigFraction[] bval = new BigFraction[] { b.getEntry(i), b.getEntry(j) }; FieldVector<BigFraction> bsys = new ArrayFieldVector<BigFraction>(bval); FieldMatrix<BigFraction> sys = LP.addBlock(line1, line2, LP.UNDER); try { FieldVector<BigFraction> point = new FieldLUDecomposition<BigFraction>(sys).getSolver() .getInverse().operate(bsys); double x = point.getEntry(0).doubleValue(); double y = point.getEntry(1).doubleValue(); Point2D p2d = new Point2D.Double(x, y); /* Only add feasible points */ if (feasible(p2d, N, b)) { if (i >= N.getRowDimension() - 1) unb.add(p2d);/* w w w.jav a2s . c om*/ points.add(p2d); } } catch (IllegalArgumentException e) { /* * Two lines that don't intersect forms an invertible * matrix. Skip these points. */ } } } return points.toArray(new Point2D[0]); }
From source file:model.LP.java
/** * Do one iteration of the simplex method. * * @param entering/* w w w .j ava2s .c o m*/ * Index of variable to enter the basis. * @param leaving * Index of variable to leave the basis. * @return * A linear program after one iteration. */ public LP pivot(int entering, int leaving) { FieldMatrix<BigFraction> bin = new FieldLUDecomposition<BigFraction>(B_).getSolver().getInverse() .multiply(N_); // Step 1: Check for optimpivality // Step 2: Select entering variable. // Naive method. Does not check for optimality. Assumes feasibility. // Entering variable is given. // Step 3: Compute primal step direction. FieldVector<BigFraction> ej = new ArrayFieldVector<BigFraction>(bin.getColumnDimension(), BigFraction.ZERO); ej.setEntry(entering, BigFraction.ONE); FieldVector<BigFraction> psd = bin.operate(ej); // Step 4: Compute primal step length. // Step 5: Select leaving variable. // Leaving variable is given. BigFraction t = b_.getEntry(leaving).divide(psd.getEntry(leaving)); // Step 6: Compute dual step direction. FieldVector<BigFraction> ei = new ArrayFieldVector<BigFraction>(bin.getRowDimension(), BigFraction.ZERO); ei.setEntry(leaving, BigFraction.ONE); FieldVector<BigFraction> dsd = bin.transpose().scalarMultiply(BigFraction.MINUS_ONE).operate(ei); // Step 7: Compute dual step length. BigFraction s = c_.getEntry(entering).divide(dsd.getEntry(entering)); // Step 8: Update current primal and dual solutions. FieldVector<BigFraction> nb_ = b_.subtract(psd.mapMultiply(t)); nb_.setEntry(leaving, t); FieldVector<BigFraction> nc_ = c_.subtract(dsd.mapMultiply(s)); nc_.setEntry(entering, s); // Step 9: Update basis. FieldVector<BigFraction> temp = B_.getColumnVector(leaving); FieldMatrix<BigFraction> nB_ = B_.copy(); nB_.setColumn(leaving, N_.getColumn(entering)); FieldMatrix<BigFraction> nN_ = N_.copy(); nN_.setColumnVector(entering, temp); int[] nBi = Bi.clone(); int[] nNi = Ni.clone(); nBi[leaving] = Ni[entering]; nNi[entering] = Bi[leaving]; return new LP(B, N, b, c, nB_, nN_, nb_, nc_, x, nBi, nNi); }
From source file:model.LP.java
public LP reinstate() { FieldVector<BigFraction> nc_ = new ArrayFieldVector<BigFraction>(c_.getDimension(), BigFraction.ZERO); FieldMatrix<BigFraction> bin = new FieldLUDecomposition<BigFraction>(B_).getSolver().getInverse() .multiply(N_);/*w w w . ja v a 2 s . c om*/ for (int i = 0; i < Bi.length; i++) { int k = Bi[i]; if (k < Ni.length) { for (int j = 0; j < Ni.length; j++) { BigFraction bf = c.getEntry(k).multiply(bin.getEntry(i, j)); nc_.setEntry(j, nc_.getEntry(j).add(bf)); } } } for (int i = 0; i < Ni.length; i++) { int k = Ni[i]; if (k < Ni.length) { nc_.setEntry(i, nc_.getEntry(i).subtract(c.getEntry(i))); } } return new LP(B, N, b, c, B_, N_, b_, nc_, x, Bi, Ni); }
From source file:model.LP.java
/** * @return/*from w w w. j a v a 2 s . com*/ * A {@code Matrix} of double precision numbers representing * the dictionary of the current Linear Program. */ public FieldMatrix<BigFraction> dictionary() { BigFraction[][] data = new BigFraction[Bi.length + 1][Ni.length + 1]; for (int i = 0; i < Ni.length; i++) { data[0][i + 1] = c_.getEntry(i).negate(); } for (int i = 0; i < Bi.length; i++) { data[i + 1][0] = b_.getEntry(i); } data[0][0] = objVal(); FieldMatrix<BigFraction> values = new FieldLUDecomposition<BigFraction>(B_).getSolver().getInverse() .multiply(N_); for (int i = 0; i < Bi.length; i++) { for (int j = 0; j < Ni.length; j++) { data[i + 1][j + 1] = values.getEntry(i, j).negate(); } } return new Array2DRowFieldMatrix<BigFraction>(data); }
From source file:org.rhwlab.BHCnotused.GaussianGIWPrior.java
public void init() { int n = data.size(); int d = m.getDimension(); Dfp rP = r.add(n);//from w ww. java2 s .com // System.out.printf("rP=%s\n",rP.toString()); double nuP = nu + n; // System.out.printf("nuP=%e\n", nuP); FieldMatrix C = new Array2DRowFieldMatrix(field, d, d); for (int row = 0; row < C.getRowDimension(); ++row) { for (int col = 0; col < C.getColumnDimension(); ++col) { C.setEntry(row, col, field.getZero()); } } FieldVector X = new ArrayFieldVector(field, d); // a vector of zeros for (FieldVector v : data) { X = X.add(v); FieldMatrix v2 = v.outerProduct(v); C = C.add(v2); } FieldVector mP = (m.mapMultiply(r).add(X)).mapDivide(r.add(n)); FieldMatrix Sp = C.add(S); FieldMatrix rmmP = mP.outerProduct(mP).scalarMultiply(rP); Sp = Sp.add(rmm).subtract(rmmP); FieldLUDecomposition ed = new FieldLUDecomposition(Sp); Dfp det = (Dfp) ed.getDeterminant(); Dfp detSp = det.pow(field.newDfp(nuP / 2.0)); Dfp gamma = field.getOne(); Dfp gammaP = field.getOne(); for (int i = 1; i <= d; ++i) { gamma = gamma.multiply(Gamma.gamma((nu + 1 - i) / 2.0)); gammaP = gammaP.multiply(Gamma.gamma((nuP + 1 - i) / 2.0)); } Dfp t1 = field.getPi().pow(-n * d / 2.0); Dfp t2 = r.divide(rP).pow(d / 2.0); Dfp t3 = detS.divide(detSp); Dfp t4 = gammaP.divide(gamma); Dfp t34 = t3.multiply(t4); /* System.out.printf("detSp=%s\n", detSp.toString()); System.out.printf("det=%s\n", det.toString()); System.out.printf("gamma=%s\n", gamma.toString()); System.out.printf("gammaP=%s\n", gammaP.toString()); System.out.printf("t1=%s\n", t1.toString()); System.out.printf("t2=%s\n", t2.toString()); System.out.printf("t3=%s\n", t3.toString()); System.out.printf("t4=%s\n", t4.toString()); */ likelihood = t2.multiply(t34).multiply(t1); double realLike = likelihood.getReal(); // System.out.printf("Likelihood=%e\n", realLike); int uhfd = 0; }
From source file:org.rhwlab.BHCnotused.GaussianGIWPrior.java
static void setParameters(double n, double beta, Dfp[] mu, double s) { GaussianGIWPrior.nu = n;//from w w w .j a va 2s .c o m GaussianGIWPrior.S = new Array2DRowFieldMatrix(field, mu.length, mu.length); for (int i = 0; i < mu.length; ++i) { S.setEntry(i, i, field.newDfp(s)); } GaussianGIWPrior.r = field.newDfp(beta); GaussianGIWPrior.m = new ArrayFieldVector(field, mu); // EigenDecomposition ed = new EigenDecomposition(S); FieldLUDecomposition ed = new FieldLUDecomposition(S); detS = ((Dfp) ed.getDeterminant()); System.out.printf("detS=%s\n", detS.toString()); detS = detS.pow(nu / 2.0); System.out.printf("detSnu=%s\n", detS.toString()); rmm = m.outerProduct(m).scalarMultiply(field.newDfp(r)); }