linopy example

test.py

import time

import numpy as np
from numpy import arange
import pyoptinterface as poi
from linopy import Model
from pyoptinterface import gurobi


def to_poi(m):
    m.constraints.sanitize_missings()
    M = m.matrices

    model = gurobi.Model()

    M_vlabels = M.vlabels
    names = "x" + M_vlabels.astype(str).astype(object)
    Nvariables = M_vlabels.shape[0]
    x = np.empty(Nvariables, dtype=object)
    M_vtypes = M.vtypes
    M_lb = M.lb
    M_ub = M.ub
    for i in range(Nvariables):
        vtype = M_vtypes[i]
        domain = poi.VariableDomain.Continuous
        if vtype == "B":
            domain = poi.VariableDomain.Binary
        elif vtype == "I":
            domain = poi.VariableDomain.Integer
        x[i] = model.add_variable(lb=M_lb[i], ub=M_ub[i], name=names[i], domain=domain)

    obj = poi.ExprBuilder()
    if m.is_quadratic:
        Q = M.Q
        for col in range(Q.shape[1]):
            # For each column, slice the relevant section of the indices and data arrays
            start_index = Q.indptr[col]
            end_index = Q.indptr[col + 1]
            rows = Q.indices[start_index:end_index]
            values = Q.data[start_index:end_index]

            # Iterate through rows and values
            for row, value in zip(rows, values):
                obj.add_quadratic_term(x[col], x[row], value)
    M_c = M.c
    for i in range(Nvariables):
        obj.add_affine_term(x[i], M_c[i])

    if m.objective.sense == "max":
        obj = -1.0 * obj

    model.set_objective(obj)

    M_A = M.A
    M_clabels = M.clabels
    if len(m.constraints):
        Nconstraints = M_A.shape[0]
        names = "c" + M_clabels.astype(str).astype(object)

        # convert A to csr_matrix
        A = M_A.tocsr()
        # iterate through rows
        M_sense = M.sense
        M_b = M.b
        for row in range(Nconstraints):
            # For each row, slice the relevant section of the indices and data arrays
            start_index = A.indptr[row]
            end_index = A.indptr[row + 1]
            cols = A.indices[start_index:end_index]
            values = A.data[start_index:end_index]

            # Iterate through columns and values
            expr = poi.ExprBuilder()
            for col, value in zip(cols, values):
                expr.add_affine_term(x[col], value)
            char_sense = M_sense[row]
            rhs = M_b[row]
            sense = {"<": poi.Leq, ">": poi.Geq, "=": poi.Eq, }[char_sense]
            c = model.add_linear_constraint(expr, sense, rhs, name=names[row])

    return model


def create_model(N):
    m = Model()
    x = m.add_variables(coords=[arange(N), arange(N)])
    y = m.add_variables(coords=[arange(N), arange(N)])
    m.add_constraints(x - y >= arange(N))
    m.add_constraints(x + y >= 0)
    m.add_objective((x * x).sum() + y.sum())
    return m


def create_poi_model(N):
    m = gurobi.Model()
    x = m.add_variables(range(N), range(N))
    y = m.add_variables(range(N), range(N))
    for i in range(N):
        for j in range(N):
            m.add_linear_constraint(x[i, j] - y[i, j], poi.Geq, i)
            m.add_linear_constraint(x[i, j] + y[i, j], poi.Geq, 0)
    expr = poi.ExprBuilder()
    poi.quicksum_(expr, x, lambda x: 2 * x)
    poi.quicksum_(expr, y)
    m.set_objective(expr)
    return m


N = 800

t0 = time.time()
model = create_model(N)
t1 = time.time()
print('Time to create linopy model:', t1 - t0)

t0 = time.time()
raw_model = model.to_gurobipy()
t1 = time.time()
print('Time to convert to gurobipy model:', t1 - t0)

t0 = time.time()
raw_model = to_poi(model)
t1 = time.time()
print('Time to convert to poi model:', t1 - t0)

t0 = time.time()
raw_model = create_poi_model(N)
t1 = time.time()
print('Time to create poi model:', t1 - t0)