Examples with GRBVar gurobi.GRBVar used on opensource projects

Example 1 with GRBVar

use of gurobi.GRBVar in project chordatlas by twak.

the class GurobiSkelSolver method dumpModelStats.

private void dumpModelStats(List<HalfEdge> edges) throws GRBException {
    print("profile count " + (globalProfs == null ? 0 : globalProfs.size()));
    print("acute corner constraints " + countBadCorners);
    print("parallel edge constraints " + countBadEdges);
    print("nearby profile terms " + countNearbyProfiles);
    print("total half edges " + edges.size());
    print("total faces " + faceInfo.size());
    if (profFit.get(edges.get(0)) != null)
        print("total profiles " + profFit.get(edges.get(0)).length);
    print("total edge length " + totalEdgeLength);
    print("\nobjective fn breakdown follows...\n");
    Cache<String, Double> penalties = new Cach<String, Double>(x -> 0.);
    Cache<String, Integer> count = new Cach<String, Integer>(x -> 0);
    for (GRBVar v : model.getVars()) {
        penalties.cache.put(v.get(StringAttr.VarName), penalties.get(v.get(StringAttr.VarName)) + v.get(GRB.DoubleAttr.Obj) * v.get(GRB.DoubleAttr.X));
        count.cache.put(v.get(StringAttr.VarName), count.get(v.get(StringAttr.VarName)) + 1);
    }
    List<String> vals = new ArrayList<>(penalties.cache.keySet());
    Collections.sort(vals);
    double obj = penalties.cache.values().stream().mapToDouble(x -> x).sum();
    for (String k : vals) print(k + "(" + count.get(k) + ") :: " + penalties.get(k) + " ");
    print("\ntotal obj " + obj + "\n\n");
}

Example 2 with GRBVar

use of gurobi.GRBVar in project chordatlas by twak.

the class GurobiSkelSolver method buildMini.

private void buildMini() throws GRBException {
    Cache<HalfEdge, GRBLinExpr> isMiniExpr = new Cach<HalfMesh2.HalfEdge, GRBLinExpr>(he -> new GRBLinExpr());
    if (minis == null)
        return;
    for (MegaFeatures mf : minis.keySet()) {
        List<MiniPtCluster> clusterVars = new ArrayList<>();
        mega2clusters.put(mf, clusterVars);
        double mLen = mf.megafacade.length();
        DumbCluster1D<MFPoint> clusters = clusterMinis(mf, minis);
        for (DumbCluster1D.Cluster<MFPoint> d : clusters) {
            // for each cluster, we pick a single MFPoint as the boundary
            MiniPtCluster miniPtVar = new MiniPtCluster();
            Cache<HalfEdge, GRBLinExpr> isEdgeBind = new Cach<HalfMesh2.HalfEdge, GRBLinExpr>(he -> new GRBLinExpr());
            GRBLinExpr selectHE = new GRBLinExpr();
            boolean one = false;
            for (MFPoint pt : d.things) {
                MiniSelectEdge miniSelectEdge = new MiniSelectEdge();
                miniPtVar.put(pt, miniSelectEdge);
                List<HalfEdge> nearCorners = findNear(mf.megafacade, pt, mesh);
                try {
                    for (HalfEdge he : nearCorners) {
                        GRBVar heVar = model.addVar(0.0, 1.0, 0, GRB.BINARY, MINI_TO_EDGE);
                        miniSelectEdge.edge.put(he, heVar);
                        selectHE.addTerm(1, heVar);
                        double cost = he.end.distance(pt);
                        if (he.over != null) {
                            if (pt.right != null)
                                cost += Math.abs(pt.right.height - ((SuperFace) he.face).height);
                            if (pt.left != null)
                                cost += Math.abs(pt.left.height - ((SuperFace) he.over.face).height);
                            isEdgeBind.get(he).addTerm(1, heVar);
                        } else
                            // bonus for being on a corner;
                            cost -= totalEdgeLength * 0.1;
                        target.addTerm(cost, heVar);
                        isMiniExpr.get(he).addTerm(1, heVar);
                        one = true;
                    }
                } catch (Throwable th) {
                    th.printStackTrace();
                }
            }
            if (one) {
                clusterVars.add(miniPtVar);
                model.addConstr(selectHE, GRB.EQUAL, 1, "pick one near " + d.things.iterator().next());
            } else
                print("warning skipping minifacade loction " + d.things.size());
            for (HalfEdge he : isEdgeBind.cache.keySet()) model.addConstr(isEdgeBind.get(he), GRB.EQUAL, edgeInfo.get(he).isEdge, "minifacade boundary must terminate on edge " + he);
            miniPtVar.mean = mf.megafacade.fromPPram(d.mean / mLen);
        }
    }
    double penalty = totalEdgeLength * 0.1;
    for (HalfEdge he : edges) {
        if (he.over != null && he.next.over == null) {
            // edge comes to boundary without minifacade --> penalty
            OptionalDouble miniDist = minis.keySet().stream().map(mf -> mf.megafacade).mapToDouble(line -> line.distance(he.end, true)).min();
            if (!miniDist.isPresent() || miniDist.getAsDouble() > 4)
                continue;
            EdgeVars ei = edgeInfo.get(he);
            ei.edgeNoMini = model.addVar(0.0, 1.0, 0, GRB.BINARY, EDGE_NO_MINI);
            if (isMiniExpr.cache.containsKey(he)) {
                ei.isMini = model.addVar(0.0, 1.0, 0, GRB.BINARY, IS_MINI);
                GRBLinExpr is = isMiniExpr.get(he);
                /* ei.isMini might take 0 or positive integer... assume it's below 10 (0.1 = 1/10) */
                model.addConstr(ei.isMini, GRB.LESS_EQUAL, is, "is minifacade on edge " + he);
                model.addConstr(scale(0.1, is), GRB.LESS_EQUAL, ei.isMini, "is minifacade on edge " + he);
                GRBLinExpr expr = new GRBLinExpr();
                expr.addTerm(1, ei.isEdge);
                expr.addTerm(1, ei.isMini);
                model.addConstr(ei.edgeNoMini, GRB.LESS_EQUAL, expr, null);
                expr = new GRBLinExpr();
                expr.addTerm(1, ei.edgeNoMini);
                expr.addTerm(1, ei.isMini);
                model.addConstr(expr, GRB.LESS_EQUAL, 1, null);
                expr = new GRBLinExpr();
                expr.addTerm(1, ei.isEdge);
                expr.addTerm(-1, ei.isMini);
                model.addConstr(ei.edgeNoMini, GRB.GREATER_EQUAL, expr, null);
            } else {
                // no mini, but easier debug
                model.addConstr(ei.edgeNoMini, GRB.EQUAL, ei.isEdge, null);
            }
            target.addTerm(penalty, ei.edgeNoMini);
        }
    }
}

Example 3 with GRBVar

use of gurobi.GRBVar in project chordatlas by twak.

the class SliceSolver method sliceOptimize.

public static int[] sliceOptimize(int numSlices, double[][] data, double[][] alignment) {
    try {
        GRBEnv env = new GRBEnv("mip1.log");
        GRBModel model = new GRBModel(env);
        GRBVar[][] xij = new GRBVar[numSlices][numSlices];
        for (int i = 0; i < numSlices; i++) for (int j = 0; j < numSlices; j++) xij[i][j] = model.addVar(0.0, 1.0, 1.0, GRB.BINARY, "x" + i + "_" + j);
        // overlap term n^2 terms
        for (int i = 0; i < numSlices; i++) {
            GRBLinExpr expr = new GRBLinExpr();
            for (int j = 0; j < numSlices; j++) expr.addTerm(1, xij[i][j]);
            model.addConstr(expr, GRB.EQUAL, 1, "overlap_" + i);
        }
        // for ( int i = 0; i < numSlices-1; i++ ) {
        // for ( int j = 0; j < numSlices; j++ ) {
        // GRBLinExpr expr = new GRBLinExpr();
        // expr.addTerm(  1, xij[ i   ] [ j ] );
        // expr.addTerm( -1, xij[ i+1 ] [ j ] );
        // model.addConstr( expr, GRB.EQUAL, 0, "super_adjacency" + j );
        // }
        // }
        // data fitting term n^2 terms
        GRBQuadExpr target = new GRBQuadExpr();
        for (int i = 0; i < numSlices; i++) {
            for (int j = 0; j < numSlices; j++) {
                if (data[i][j] != Double.MAX_VALUE)
                    target.addTerm(data[i][j] * 30, xij[i][j]);
                else
                    target.addTerm(1e3, xij[i][j]);
            }
        }
        // neighbour alignment n^3 quadratic terms
        for (int i = 0; i < numSlices - 1; i++) for (int ja = 0; ja < numSlices; ja++) for (int jb = 0; jb < numSlices; jb++) target.addTerm(alignment[ja][jb] == Double.MAX_VALUE ? 1e2 : alignment[ja][jb], xij[i][ja], xij[i + 1][jb]);
        model.setObjective(target, GRB.MINIMIZE);
        model.getEnv().set(GRB.DoubleParam.TimeLimit, 10.0);
        model.optimize();
        // xij[0][0].set( DoubleAttr.X, 1 );
        // target.getValue();
        System.out.println("Obj: " + model.get(GRB.DoubleAttr.ObjVal));
        int[] out = new int[numSlices];
        for (int i = 0; i < numSlices; i++) {
            System.out.print("i: " + i + " ");
            for (int j = 0; j < numSlices; j++) {
                System.out.print(xij[i][j].get(GRB.DoubleAttr.X) + " ");
                if (xij[i][j].get(GRB.DoubleAttr.X) == 1)
                    out[i] = j;
            }
            System.out.println();
        }
        // Dispose of model and environment
        model.dispose();
        env.dispose();
        return out;
    } catch (GRBException e) {
        System.out.println("Error code: " + e.getErrorCode() + ". " + e.getMessage());
        e.printStackTrace();
    }
    return null;
}

Example 4 with GRBVar

use of gurobi.GRBVar in project chordatlas by twak.

the class GurobiSkelSolver method solve.

public void solve() {
    try {
        for (HalfFace f : mesh.faces) for (HalfEdge e : f.edges()) {
            edges.add(e);
            totalEdgeLength += e.length();
        }
        print("total edge length " + totalEdgeLength);
        progress.setProgress((int) (Math.random() * 90));
        buildProblem();
        model.getEnv().set(GRB.DoubleParam.TimeLimit, GRB.INFINITY);
        long time = System.currentTimeMillis();
        new Thread(() -> {
            while (!gurobiDone) {
                try {
                    Thread.sleep(300);
                    if ((System.currentTimeMillis() - time) / 1000 > maxTimeSecs) {
                        print("time's up");
                        model.terminate();
                        progress.close();
                        return;
                    }
                } catch (InterruptedException e1) {
                }
                progress.setProgress((int) (Math.random() * 90));
                if (progress.isCanceled()) {
                    print("user cancelled");
                    model.terminate();
                }
            }
            progress.setProgress(100);
        }).start();
        try {
            String file = Tweed.SCRATCH + "problem_" + System.currentTimeMillis() + ".mps";
            new File(file).getParentFile().mkdirs();
            model.write(file);
            model.optimize();
        } finally {
            gurobiDone = true;
            progress.setProgress(100);
        }
        print("time " + (System.currentTimeMillis() - time) / 1000. + " seconds");
        if (model.get(GRB.IntAttr.Status) == GRB.Status.INFEASIBLE) {
            print("Can't solve; won't solve");
            return;
        }
        for (HalfFace f : mesh.faces) ((SuperFace) f).classification = index(faceInfo.get(f).color);
        if (minis != null)
            for (MegaFeatures mf : minis.keySet()) {
                if (mega2clusters.containsKey(mf))
                    for (MiniPtCluster pt : mega2clusters.get(mf)) {
                        for (Map.Entry<MFPoint, MiniSelectEdge> selectPt : pt.entrySet()) {
                            for (Map.Entry<HalfEdge, GRBVar> selectEdge : selectPt.getValue().edge.entrySet()) {
                                if (selectEdge.getValue().get(GRB.DoubleAttr.X) > 0.5) {
                                    HalfEdge e = selectEdge.getKey();
                                    selectPt.getKey().selectedEdge = e;
                                    selectPt.getKey().sameCluster = pt.keySet();
                                    if (e.over != null && e.next.over == null) {
                                        for (MFPoint p : pt.keySet()) ((SuperEdge) e.next).addMini(p.right);
                                        for (MFPoint p : pt.keySet()) ((SuperEdge) e.over.findBefore()).addMini(p.left);
                                    }
                                    if (e.over == null && e.next.over == null) {
                                        Line parallel = selectPt.getKey().mega.megafacade;
                                        if (parallel.absAngle(e.line()) < parallel.absAngle(e.next.line()))
                                            for (MFPoint p : pt.keySet()) ((SuperEdge) e).addMini(p.left);
                                        else
                                            for (MFPoint p : pt.keySet()) ((SuperEdge) e.next).addMini(p.right);
                                    }
                                }
                            }
                        }
                    }
            }
        if (false)
            for (HalfEdge he : edges) {
                if (he.over != null && he.next.over == null) {
                    // face comes to boundary
                    boolean viz = false;
                    EdgeVars ei = edgeInfo.get(he);
                    if (ei.edgeNoMini != null) {
                        viz = true;
                    }
                    if (viz)
                        PaintThing.debug.put("dd", new Point2d(he.end));
                }
            }
        if (globalProfs != null)
            for (HalfEdge e : edges) {
                EdgeVars ei = edgeInfo.get(e);
                GRBVar plot = ei.debug;
                if (plot != null) {
                    // && plot.get(GRB.DoubleAttr.X) > 0.5 ) {
                    Point2d o = new Point2d(e.end);
                    PaintThing.debug.put(1, o);
                }
                int p = index(ei.profile);
                ((SuperEdge) e).profI = p;
                ((SuperEdge) e).prof = p < 0 ? null : globalProfs.get(p);
            }
        dumpModelStats(edges);
        model.getEnv().dispose();
        model.dispose();
    } catch (GRBException e) {
        print("Error code: " + e.getErrorCode() + ". " + e.getMessage());
        e.printStackTrace();
    }
}

Example 5 with GRBVar

use of gurobi.GRBVar in project chordatlas by twak.

the class GurobiSkelSolver method buildProfiles.

private void buildProfiles() throws GRBException {
    for (HalfEdge e : edges) {
        if (// when a profile ends, we assume it can't start again...
        ((SuperEdge) e).profLine == null)
            continue;
        EdgeVars ev = edgeInfo.get(e);
        ev.profile = new GRBVar[globalProfs.size()];
        for (int p = 0; p < globalProfs.size(); p++) {
            ev.profile[p] = model.addVar(0.0, 1.0, 0, GRB.BINARY, PROFILE_SELECT);
            set(ev.profile[p], p == 0 ? 1 : 0);
        }
    }
    Cache2<HalfEdge, HalfEdge, GRBVar> isProfileDifferent = new Cache2<HalfMesh2.HalfEdge, HalfMesh2.HalfEdge, GRBVar>() {

        @Override
        public GRBVar create(HalfEdge e1, HalfEdge e2) {
            try {
                GRBVar s = model.addVar(0.0, 1.0, 0.0, GRB.BINARY, PROFILE_DIFFERENT);
                buildIsDifferentColor(s, edgeInfo.get(e1).profile, edgeInfo.get(e2).profile, "profile");
                s.set(DoubleAttr.Start, 0);
                return s;
            } catch (GRBException e) {
                e.printStackTrace();
            }
            return null;
        }
    };
    for (HalfEdge e1 : edges) {
        // if (e1.over != null)
        // continue;
        EdgeVars ev = edgeInfo.get(e1);
        if (ev.profile == null)
            continue;
        double[] fit = profFit.get(e1);
        GRBLinExpr pickOneProfile = new GRBLinExpr();
        for (int p = 0; p < globalProfs.size(); p++) {
            GRBVar a = ev.profile[p];
            pickOneProfile.addTerm(1, a);
            target.addTerm(.001 * fit[p] * e1.length(), a);
        }
        model.addConstr(pickOneProfile, GRB.EQUAL, 1, /*ev.isEdge*/
        "pick only one profile for " + e1);
        List<HalfEdge> atEnd = e1.collectAroundEnd();
        for (HalfEdge e2 : atEnd) {
            if (e1 != e2 && edgeInfo.get(e2).profile != null && // only constrain if ~parallel
            e1.line().absAngle(e2.line()) < 0.1) {
                GRBVar s = isProfileDifferent.get(e1, e2);
                // ev.debug = s;
                // Point2d dbg = new Point2d(e1.end);
                // dbg.add(new Point2d(Math.random() * 0.1, Math.random() * 0.1));
                GRBLinExpr perpIsEdge = new GRBLinExpr();
                for (HalfEdge e3 : atEnd) {
                    if (e3 == e1 || e3 == e2 || e3 == e1.over || e3 == e2.over || (e3.over != null && (e3.over == e1.over || e3.over == e2.over)))
                        continue;
                    if (!e1.line().isOnLeft(e1.end.distanceSquared(e3.start) > e1.end.distanceSquared(e3.end) ? e3.start : e3.end))
                        continue;
                    perpIsEdge.addTerm(1, edgeInfo.get(e3).isEdge);
                // PaintThing.debug.put(e2, dbg);
                // Point2d d2 = new Point2d(e3.line().dir());
                // d2.scale (0.1/new Vector2d(d2).length());
                // d2.add(dbg);
                // PaintThing.debug.put(e2, new Line (dbg, d2));
                // model.addConstr( s, GRB.LESS_EQUAL, edgeInfo.get(e3).isEdge, "only change profile if no adjacent edge "+e1 );
                }
                model.addConstr(s, GRB.LESS_EQUAL, perpIsEdge, "dont' change profile over " + e1);
            // PaintThing.debug.put(e2, dbg);
            // Point2d d2 = new Point2d(e2.line().dir());
            // d2.scale (0.2/new Vector2d(d2).length());
            // d2.add(dbg);
            // PaintThing.debug.put(e2, new Line (dbg, d2));
            }
        }
    }
    countNearbyProfiles = 0;
    if (false)
        for (int i = 0; i < edges.size(); i++) {
            HalfEdge e1 = edges.get(i);
            if (edgeInfo.get(e1).profile == null)
                continue;
            print("building edge locality term " + i + " / " + edges.size());
            for (HalfEdge e2 : edges) if (lt(e1, e2) && edgeInfo.get(e2).profile != null) {
                if (e1.line().distance(e2.line()) < 2) {
                    GRBVar s = isProfileDifferent.get(e1, e2);
                    target.addTerm(0.1 * (e1.length() + e2.length()), s);
                    countNearbyProfiles++;
                }
            }
        }
}