# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import numpy as np
from scipy.optimize import linprog
import matplotlib.pyplot as plt

v_t = np.loadtxt('../data/quadrature_debug/v_t.txt')
b = np.loadtxt('../data/quadrature_debug/b.txt')
w0 = np.loadtxt('../data/quadrature_debug/w0.txt')
#w = np.loadtxt('../w.txt')
#points = np.loadtxt('../points.txt')
# r = v_t.dot(w) - b
# np.linalg.norm(r)

import gurobipy as gp
from gurobipy import GRB

# +
try:
    num_vars = v_t.shape[1]
    num_constraints = v_t.shape[0]
    V = sum(w0)
    
    # Create a new model
    model = gp.Model("milp")

    w_vars = [ model.addVar(vtype=GRB.CONTINUOUS, name="w_{}".format(i)) for i in range(0, num_vars)]
    b_vars = [ model.addVar(vtype=GRB.BINARY, obj=1, name="b_{}".format(i)) for i in range(0, num_vars)]
    
    model.addMConstrs(v_t, w_vars, '=', b)
    
    for i in range(0, num_vars):
        model.addConstr(w_vars[i] <= V * b_vars[i])
    

#     # TODO: Add P constraints
#     def create_V_t_row_constraint(V_t, b, i):
#         constraint = solver.Constraint(b[i], b[i], 'V_t_{}'.format(i))
#         for j in range(0, num_vars):
#             constraint.SetCoefficient(w_vars[j], V_t[i, j])

#     for j in range(0, num_vars):
#         constraint = solver.Constraint(-solver.infinity(), 0, "w_b_{}".format(j))
#         constraint.SetCoefficient(w_vars[j], 1)
#         constraint.SetCoefficient(b_vars[j], -V)

#     P_constraints = [ create_V_t_row_constraint(v_t, b, i) for i in range(0, num_constraints) ]
#     objective = solver.Objective()
#     for b_i in b_vars:
#         objective.SetCoefficient(b_i, 1)

#     solver.EnableOutput()
#     print(solver.Solve())
#     objective.Value()

#     w_glop = np.array([ w_i.solution_value() for w_i in w_vars ])
#     b_glop = np.array([ b_i.solution_value() for b_i in b_vars ])

#     r = v_t.dot(w_glop) - b
#     print(np.linalg.norm(r))


#     # Create variables
#     x = m.addVar(vtype=GRB.BINARY, name="x")
#     y = m.addVar(vtype=GRB.BINARY, name="y")
#     z = m.addVar(vtype=GRB.BINARY, name="z")

#     # Set objective
#     m.setObjective(x + y + 2 * z, GRB.MAXIMIZE)

#     # Add constraint: x + 2 y + 3 z <= 4
#     m.addConstr(x + 2 * y + 3 * z <= 4, "c0")

#     # Add constraint: x + y >= 1
#     m.addConstr(x + y >= 1, "c1")

    # Optimize model
    model.optimize()
    
    w_gurobi = np.array([w.x for w in w_vars])
    b_gurobi = np.array([b.x for b in b_vars])

#     for v in m.getVars():
#         print('%s %g' % (v.varName, v.x))

#     print('Obj: %g' % m.objVal)

except gp.GurobiError as e:
    print('Error code ' + str(e.errno) + ': ' + str(e))

except AttributeError as e:
    print('Encountered an attribute error: ' + str(e))
# -

num_nonzero_gurobi = len([w_i for w_i in w_gurobi if w_i > 1e-14])
print("Num nonzero gurobi: {}".format(num_nonzero_gurobi))

plt.figure(figsize=(20,20))
plt.plot(w_gurobi, color="green",  linewidth=2.5, linestyle="-", label="w_gurobi")
# plt.plot(b_gurobi, color="blue",  linewidth=2.5, linestyle="-", label="b_gurobi")
plt.legend(loc='upper left', frameon=False)
plt.show()

# Try the same thing with Google's GLOP solver
from ortools.linear_solver import pywraplp

# +
# num_vars = v_t.shape[1]
# num_constraints = v_t.shape[0]
# V = sum(w0)
# solver = pywraplp.Solver('quadrature_lp', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

# w_vars = [ solver.NumVar(0.0, float("inf"), "w_{}".format(i)) for i in range(0, num_vars)]
# b_vars = [ solver.IntVar(0.0, 1.0, "b_{}".format(i)) for i in range(0, num_vars)]

# # TODO: Add P constraints
# def create_V_t_row_constraint(V_t, b, i):
#     constraint = solver.Constraint(b[i], b[i], 'V_t_{}'.format(i))
#     for j in range(0, num_vars):
#         constraint.SetCoefficient(w_vars[j], V_t[i, j])
        
# for j in range(0, num_vars):
#     constraint = solver.Constraint(-solver.infinity(), 0, "w_b_{}".format(j))
#     constraint.SetCoefficient(w_vars[j], 1)
#     constraint.SetCoefficient(b_vars[j], -V)
        
# P_constraints = [ create_V_t_row_constraint(v_t, b, i) for i in range(0, num_constraints) ]
# objective = solver.Objective()
# for b_i in b_vars:
#     objective.SetCoefficient(b_i, 1)

# solver.EnableOutput()
# print(solver.Solve())
# objective.Value()

# w_glop = np.array([ w_i.solution_value() for w_i in w_vars ])
# b_glop = np.array([ b_i.solution_value() for b_i in b_vars ])

# r = v_t.dot(w_glop) - b
# print(np.linalg.norm(r))
# -

plt.figure(figsize=(20,20))
plt.plot(w_glop, color="green",  linewidth=2.5, linestyle="-", label="w_glop")
plt.plot(b_glop, color="blue",  linewidth=2.5, linestyle="-", label="b_glop")
plt.legend(loc='upper left', frameon=False)
plt.show()

num_nonzero_glop = len([w_i for w_i in w_glop if w_i > 1e-14])
print("Num nonzero glop: {}".format(num_nonzero_glop))

# +
# from bokeh.plotting import figure, show, output_file, output_notebook
# from math import sqrt, pi

# output_notebook()

# def plot_points(points, weights):
#     x = points[:, 0]
#     y = points[:, 1]
#     print(sum(weights))
#     radii = [ sqrt(w_i / pi) for w_i in weights ]
#     colors = [
#         "#%02x%02x%02x" % (int(r), int(g), 150) for r, g in zip(50+2*x, 30+2*y)
#     ]

#     TOOLS="hover,crosshair,pan,wheel_zoom,zoom_in,zoom_out,box_zoom,undo,redo,reset,tap,save,box_select,poly_select,lasso_select,"

#     p = figure(tools=TOOLS)

#     p.scatter(x, y, radius=radii,
#               fill_color=colors, fill_alpha=0.6,
#               line_color=None)

#     show(p)
# -



StarcoderdataJupyterScriptsDedupFiltered