public MGMParams(DoubleMatrix1D vec, int p, int ltot){ int[] lens = {p*p, p, p*ltot, ltot*ltot, p, ltot}; int[] lenSums = new int[lens.length]; lenSums[0] = lens[0]; for(int i = 1; i < lenSums.length; i++){ lenSums[i] = lens[i] + lenSums[i-1]; } if(vec.size() != lenSums[5]) throw new IllegalArgumentException("Param vector dimension doesn't match: Found " + vec.size() + " need " + lenSums[5]); beta = DoubleFactory2D.dense.make(vec.viewPart(0, lens[0]).toArray(), p); betad = vec.viewPart(lenSums[0], lens[1]).copy(); theta = DoubleFactory2D.dense.make(vec.viewPart(lenSums[1], lens[2]).toArray(), ltot); phi = DoubleFactory2D.dense.make(vec.viewPart(lenSums[2], lens[3]).toArray(), ltot); alpha1 = vec.viewPart(lenSums[3], lens[4]).copy(); alpha2 = vec.viewPart(lenSums[4], lens[5]).copy(); }
/** */ public static void doubleTest31(int size) { System.out.println("\ninit"); DoubleMatrix1D a = Factory1D.dense.descending(size); DoubleMatrix1D b = new WrapperDoubleMatrix1D(a); DoubleMatrix1D c = b.viewPart(2,3); DoubleMatrix1D d = c.viewFlip(); //DoubleMatrix1D c = b.viewFlip(); //DoubleMatrix1D d = c.viewFlip(); d.set(0,99); b = b.viewSorted(); System.out.println("a = "+a); System.out.println("b = "+b); System.out.println("c = "+c); System.out.println("d = "+d); System.out.println("done"); } /**
/** */ public static void doubleTest31(int size) { System.out.println("\ninit"); DoubleMatrix1D a = Factory1D.dense.descending(size); DoubleMatrix1D b = new WrapperDoubleMatrix1D(a); DoubleMatrix1D c = b.viewPart(2,3); DoubleMatrix1D d = c.viewFlip(); //DoubleMatrix1D c = b.viewFlip(); //DoubleMatrix1D d = c.viewFlip(); d.set(0,99); b = b.viewSorted(); System.out.println("a = "+a); System.out.println("b = "+b); System.out.println("c = "+c); System.out.println("d = "+d); System.out.println("done"); } /**
/** C = A||A||..||A; Constructs a new matrix which is concatenated <tt>repeat</tt> times. Example: <pre> 0 1 repeat(3) --> 0 1 0 1 0 1 </pre> */ public DoubleMatrix1D repeat(DoubleMatrix1D A, int repeat) { int size = A.size(); DoubleMatrix1D matrix = make(repeat * size); for (int i=repeat; --i >= 0; ) { matrix.viewPart(size*i,size).assign(A); } return matrix; } /**
/** C = A||A||..||A; Constructs a new matrix which is concatenated <tt>repeat</tt> times. Example: <pre> 0 1 repeat(3) --> 0 1 0 1 0 1 </pre> */ public DoubleMatrix1D repeat(DoubleMatrix1D A, int repeat) { int size = A.size(); DoubleMatrix1D matrix = make(repeat * size); for (int i=repeat; --i >= 0; ) { matrix.viewPart(size*i,size).assign(A); } return matrix; } /**
@Override public double value(double[] Y) { DoubleMatrix1D y = DoubleFactory1D.dense.make(Y); DoubleMatrix1D X = y.viewPart(0, originalDim); return originalFi.value(X.toArray()) - y.get(dim-1); }
@Override public double value(double[] Y) { DoubleMatrix1D y = DoubleFactory1D.dense.make(Y); DoubleMatrix1D X = y.viewPart(0, dim); return originalFi.value(X.toArray()) - y.get(dimPh1-1); }
/** C = A||B; Constructs a new matrix which is the concatenation of two other matrices. Example: <tt>0 1</tt> append <tt>3 4</tt> --> <tt>0 1 3 4</tt>. */ public DoubleMatrix1D append(DoubleMatrix1D A, DoubleMatrix1D B) { // concatenate DoubleMatrix1D matrix = make(A.size()+B.size()); matrix.viewPart(0,A.size()).assign(A); matrix.viewPart(A.size(),B.size()).assign(B); return matrix; } /**
/** C = A||B; Constructs a new matrix which is the concatenation of two other matrices. Example: <tt>0 1</tt> append <tt>3 4</tt> --> <tt>0 1 3 4</tt>. */ public DoubleMatrix1D append(DoubleMatrix1D A, DoubleMatrix1D B) { // concatenate DoubleMatrix1D matrix = make(A.size()+B.size()); matrix.viewPart(0,A.size()).assign(A); matrix.viewPart(A.size(),B.size()).assign(B); return matrix; } /**
/** Generates and returns the (economy-sized) orthogonal factor <tt>Q</tt>. @return <tt>Q</tt> */ public DoubleMatrix2D getQ () { cern.jet.math.Functions F = cern.jet.math.Functions.functions; DoubleMatrix2D Q = QR.like(); //double[][] Q = X.getArray(); for (int k = n-1; k >= 0; k--) { DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k); Q.setQuick(k,k, 1); for (int j = k; j < n; j++) { if (QR.getQuick(k,k) != 0) { DoubleMatrix1D Qcolj = Q.viewColumn(j).viewPart(k,m-k); double s = QRcolk.zDotProduct(Qcolj); s = -s / QR.getQuick(k,k); Qcolj.assign(QRcolk, F.plusMult(s)); } } } return Q; } /**
/** Generates and returns the (economy-sized) orthogonal factor <tt>Q</tt>. @return <tt>Q</tt> */ public DoubleMatrix2D getQ () { cern.jet.math.Functions F = cern.jet.math.Functions.functions; DoubleMatrix2D Q = QR.like(); //double[][] Q = X.getArray(); for (int k = n-1; k >= 0; k--) { DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k); Q.setQuick(k,k, 1); for (int j = k; j < n; j++) { if (QR.getQuick(k,k) != 0) { DoubleMatrix1D Qcolj = Q.viewColumn(j).viewPart(k,m-k); double s = QRcolk.zDotProduct(Qcolj); s = -s / QR.getQuick(k,k); Qcolj.assign(QRcolk, F.plusMult(s)); } } } return Q; } /**
@Override protected boolean checkCustomExitConditions(DoubleMatrix1D Y){ DoubleMatrix1D X = Y.viewPart(0, getDim()-1); //equalities DoubleMatrix1D originalRPriX = F1.make(0); if(getA()!=null){ //originalRPriX = originalProblem.getA().zMult(X, originalProblem.getB().copy(), 1., -1., false); originalRPriX = ColtUtils.zMult(originalProblem.getA(), X, originalProblem.getB(), -1.); } boolean b2 = Math.sqrt(ALG.norm2(originalRPriX)) < originalProblem.getToleranceFeas(); //inequalities boolean b1 = !Double.isNaN(originalProblem.getBarrierFunction().value(X.toArray())) || Y.get(Y.size()-1)<0; log.debug("checkCustomExitConditions: " + (b1 && b2)); return b1 && b2; } }
@Override public double[] gradient(double[] Y) { DoubleMatrix1D y = DoubleFactory1D.dense.make(Y); DoubleMatrix1D X = y.viewPart(0, originalDim); DoubleMatrix1D origGrad = F1.make(originalFi.gradient(X.toArray())); DoubleMatrix1D ret = F1.make(1, -1); ret = F1.append(origGrad, ret); return ret.toArray(); }
@Override public double[] gradient(double[] Y) { DoubleMatrix1D y = DoubleFactory1D.dense.make(Y); DoubleMatrix1D X = y.viewPart(0, dim); DoubleMatrix1D origGrad = F1.make(originalFi.gradient(X.toArray())); DoubleMatrix1D ret = F1.make(1, -1); ret = F1.append(origGrad, ret); return ret.toArray(); }
/** Constructs a matrix which is the concatenation of all given parts. Cells are copied. */ public DoubleMatrix1D make(DoubleMatrix1D[] parts) { if (parts.length==0) return make(0); int size = 0; for (int i=0; i < parts.length; i++) size += parts[i].size(); DoubleMatrix1D vector = make(size); size = 0; for (int i=0; i < parts.length; i++) { vector.viewPart(size,parts[i].size()).assign(parts[i]); size += parts[i].size(); } return vector; } /**
@Override protected boolean checkCustomExitConditions(DoubleMatrix1D Y){ DoubleMatrix1D X = Y.viewPart(0, getDim()-1); DoubleMatrix1D ineqX = originalProblem.getFi(X); int ineqMaxIndex = Utils.getMaxIndex(ineqX); boolean isInternal = (ineqX.get(ineqMaxIndex) + getTolerance() <0) || Y.get(Y.size()-1)<0; log.info("isInternal : " + isInternal); if(!isInternal){ return false; } DoubleMatrix1D originalRPriX = F1.make(0); if(getA()!=null){ //originalRPriX = originalProblem.getA().zMult(X, originalProblem.getB().copy(), 1., -1., false); originalRPriX = ColtUtils.zMult(originalProblem.getA(), X, originalProblem.getB(), -1); } boolean isPrimalFeas = Math.sqrt(ALG.norm2(originalRPriX)) < originalProblem.getToleranceFeas(); log.info("isPrimalFeas: " + isPrimalFeas); log.info("checkCustomExitConditions: " + (isInternal && isPrimalFeas)); return isInternal && isPrimalFeas; } }
/** Constructs a matrix which is the concatenation of all given parts. Cells are copied. */ public DoubleMatrix1D make(DoubleMatrix1D[] parts) { if (parts.length==0) return make(0); int size = 0; for (int i=0; i < parts.length; i++) size += parts[i].size(); DoubleMatrix1D vector = make(size); size = 0; for (int i=0; i < parts.length; i++) { vector.viewPart(size,parts[i].size()).assign(parts[i]); size += parts[i].size(); } return vector; } /**
@Override protected boolean checkCustomExitConditions(DoubleMatrix1D Xs){ DoubleMatrix1D X = Xs.viewPart(0, getDim()-1); DoubleMatrix1D ineqX = originalProblem.getFi(X); int ineqMaxIndex = Utils.getMaxIndex(ineqX); //log.debug("ineqMaxIndex: " + ineqMaxIndex); //log.debug("ineqMaxValue: " + ineqX.get(ineqMaxIndex)); boolean isInternal = (ineqX.get(ineqMaxIndex) + getTolerance() <0) || Xs.get(Xs.size()-1)<0; log.info("isInternal : " + isInternal); if(!isInternal){ return false; } DoubleMatrix1D originalRPriX = F1.make(0); if(getA()!=null){ //originalRPriX = originalProblem.getA().zMult(X, originalProblem.getB().copy(), 1., -1., false); originalRPriX = ColtUtils.zMult(originalProblem.getA(), X, originalProblem.getB(), -1); } boolean isPrimalFeas = Math.sqrt(ALG.norm2(originalRPriX)) < originalProblem.getToleranceFeas(); log.info("isPrimalFeas: " + isPrimalFeas); log.info("checkCustomExitConditions: " + (isInternal && isPrimalFeas)); return isInternal && isPrimalFeas; }
@Override public double[][] hessian(double[] Y) { DoubleMatrix1D y = DoubleFactory1D.dense.make(Y); DoubleMatrix1D X = y.viewPart(0, dim); DoubleMatrix2D origHess; double[][] origFiHessX = originalFi.hessian(X.toArray()); if(origFiHessX == FunctionsUtils.ZEROES_2D_ARRAY_PLACEHOLDER){ return FunctionsUtils.ZEROES_2D_ARRAY_PLACEHOLDER; }else{ origHess = F2.make(origFiHessX); DoubleMatrix2D[][] parts = new DoubleMatrix2D[][]{{origHess, null},{null,F2.make(1, 1)}}; return F2.compose(parts).toArray(); } }
@Override public double[][] hessian(double[] Y) { DoubleMatrix1D y = DoubleFactory1D.dense.make(Y); DoubleMatrix1D X = y.viewPart(0, originalDim); double[][] originalFiHess = originalFi.hessian(X.toArray()); if(originalFiHess==FunctionsUtils.ZEROES_2D_ARRAY_PLACEHOLDER){ return FunctionsUtils.ZEROES_2D_ARRAY_PLACEHOLDER; }else{ DoubleMatrix2D origHess = (originalFiHess!=FunctionsUtils.ZEROES_2D_ARRAY_PLACEHOLDER)? F2.make(originalFi.hessian(X.toArray())) : F2.make(X.size(), X.size()); DoubleMatrix2D[][] parts = new DoubleMatrix2D[][]{{origHess, null},{null,F2.make(1, 1)}}; return F2.compose(parts).toArray(); } }