Duality Theorum Solver

duality.py

from copy import deepcopy

import networkx as nx
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

G = nx.DiGraph()

class Literal(str):
    def __new__(cls, literal):
        return str.__new__(cls, literal)

    def __init__(self, literal):
        self.literals = [literal]

    def __invert__(self):
        return Not(self)

    def __and__(self, other):
        return And(self, other)

    def __or__(self, other):
        return Or(self, other)

    def __copy__(self):
        return Literal(self.literals[0])

    def __deepcopy__(self, memodict={}):
        return Literal(deepcopy(self.literals[0], memodict))

class And(Literal):
    def __new__(cls, *literals):
        literals = literals
        assert all(type(xi) in (Literal, Not) for xi in literals), f"This package only supports single nesting: {[type(l) for l in literals]}"
        this = f"({str.join(' and ', literals)})"
        return str.__new__(cls, this)

    def __init__(self, *literals):
        self.literals = literals

    def __copy__(self):
        return And(*self.literals)

    def __deepcopy__(self, memodict={}):
        return And(*[deepcopy(l, memodict) for l in self.literals])

class Or(Literal):
    def __new__(cls, *literals):
        literals = literals
        assert all(type(xi) in (Literal, Not) for xi in literals), f"This package only supports single nesting: {[type(l) for l in literals]}"
        this = f"({str.join(' or ', literals)})"
        return str.__new__(cls, this)

    def __init__(self, *literals):
        self.literals = literals

    def __copy__(self):
        return Or(*self.literals)

    def __deepcopy__(self, memodict={}):
        return Or(*[deepcopy(l, memodict) for l in self.literals])

class Not(Literal):
    def __new__(cls, literal):
        if isinstance(literal, Not):
            return literal.literals[0]
        elif isinstance(literal, And):
            return Or(*[Not(xi) for xi in literal.literals])
        elif isinstance(literal, Or):
            return And(*[Not(xi) for xi in literal.literals])
        else:
            literals = [literal]
            this = f'(Not {literals[0]})'
            return str.__new__(cls, this)

    def __init__(self, literal):
        self.literals = [literal]

    def __copy__(self):
        return Not(self.literals[0])

    def __deepcopy__(self, memodict={}):
        return Not(deepcopy(self.literals[0], memodict))

def test_not():
    assert str(Not(Not(Not('a'))))=='(Not a)'
    assert str(Not(Not('a')))=='a'
    assert str(Not('a'))=='(Not a)'
    assert str(Not(Or(Literal('p'), Literal('q'))))==str(And(Not('p'), Not('q')))
    assert str(Not(And(Literal('p'), Literal('q'))))==str(Or(Not('p'), Not('q')))
    assert type(Not(Not(Not(Literal('a'))))) is Not
    assert type(Not(Not(Literal('a')))) is Literal

test_not()


PrimalFeasible = Literal('Primal Feasible')
PrimalBounded = Literal('Primal Bounded')
DualFeasible = Literal('Dual Feasible')
DualBounded = Literal('Dual Bounded')
WeakDuality = Literal('cTx <= bTy')
StrongDuality = Literal('cTx==bTy')
PrimalOptimal = Literal('Primal Optimal')
DualOptimal = Literal('Dual Optimal')
# ManualCalculation1 = Literal('ManualCalculation1')
# ManualCalculation2 = Literal('ManualCalculation2')
# ManualCalculation3 = Literal('ManualCalculation3')
# ManualCalculation4 = Literal('ManualCalculation4')

PrimalInfeasible = ~PrimalFeasible
DualInfeasible = ~DualFeasible
PrimalUnbounded = ~PrimalBounded
DualUnbounded = ~DualBounded

# Unbounded LP
G.add_edge(DualUnbounded, PrimalInfeasible, source='Unbounded LP')
G.add_edge(PrimalUnbounded, DualInfeasible, source='Unbounded LP')

# Check Unbounded
G.add_edge(DualInfeasible, PrimalInfeasible | PrimalUnbounded, source='Check Unbounded')

# Weak Duality
G.add_edge(PrimalFeasible & DualFeasible, WeakDuality, source='Weak Duality')

# Matching Values Lecture
G.add_edge(And(PrimalFeasible, DualFeasible, StrongDuality), PrimalOptimal, source='Matching Values Lecture')

# Primal Optimal iff Dual Optimal
G.add_edge(PrimalOptimal, DualOptimal, source='Strong Duality')
G.add_edge(DualOptimal, PrimalOptimal, source='Strong Duality')

# Feasible and bounded iff optimal
G.add_edge(PrimalFeasible & PrimalBounded, PrimalOptimal, source='Strong Duality')
G.add_edge(PrimalOptimal, PrimalFeasible & PrimalBounded, source='Strong Duality')
G.add_edge(DualFeasible & DualBounded, DualOptimal, source='Strong Duality')
G.add_edge(DualOptimal, DualFeasible & DualBounded, source='Strong Duality')

# My own findings
# G.add_edge(ManualCalculation1, PrimalBounded, source='manual')
# G.add_edge(ManualCalculation2, PrimalUnbounded, source='manual')
# G.add_edge(ManualCalculation3, DualBounded, source='manual')
# G.add_edge(ManualCalculation4, DualUnbounded, source='manual')

# Collect all literals
n = list(G.nodes())
print(n)

# Include all nots as contradictions
for n in list(G.nodes()):
    G.add_edge(n, Not(n), source='contradiction')

# Basic Logic

# And implies each
for n in list(G.nodes()):
    if type(n) is And:
        for a in n.literals:
            G.add_edge(n, a, source='and each')

# Include all nots as contradictions (again)
for n in list(G.nodes()):
    G.add_edge(n, Not(n), source='contradiction')


# Power graph
# def directed_power(G, k, source=None):
#     """ Networkx doesn't have a power graph implementation for directed graphs for some reason. """
#     g = G
#     G = G.copy()
#     assert len(G.edges()) == len(g.edges())
#     assert G is not g
#     new_edges = False
#     for i in range(k):
#         for n in list(G.nodes()):
#             for m1 in list(G.successors(n)):
#                 for m2 in list(G.successors(m1)):
#                     if not G.has_edge(n, m2):
#                         G.add_edge(n, m2, source=source)
#                         print(f"New Edge {n} -> {m2}")
#                         new_edges = True

#     if new_edges:
#         assert len(G.edges()) != len(g.edges())

#     return G

def directed_power(G, k, source='calculated'):
    """ Networkx doesn't have a power graph implementation for directed graphs for some reason. """
    g = G
    G = G.copy()
    assert len(G.edges()) == len(g.edges())
    assert G is not g
    new_edges = False
    for i in range(k):
        for n in list(G.nodes()):
            for m1 in list(G.successors(n)):
                if 'contradiction' in G[n][m1]['source']:
                    m1 = Not(m1)
                for m2 in list(G.successors(m1)):
                    c = 'contradiction' in G[m1][m2]['source']
                    if not G.has_edge(n, m2):
                        G.add_edge(n, m2, source=source if not c else 'contradiction '+source)
                        print(f"New Edge {n} -> {m2}")
                        new_edges = True

    if new_edges:
        assert len(G.edges()) != len(g.edges())

    return G

def test_directed_power():
    G2 = directed_power(G, 1)
    assert (len(G2.edges()) > len(G.edges()))
    # print((set(G2.edges()) - set(G.edges())))

test_directed_power()

g = G # A previous version
G = directed_power(G, 100)

# Combine my custom "Needs calculation" nodes
# ManualCalculation = Literal('ManualCalculation')
# G.add_node(ManualCalculation)
# for n in [ManualCalculation1, ManualCalculation2, ManualCalculation3, ManualCalculation4]:
#     for m in list(G.successors(n)):
#         G.add_edge(ManualCalculation, m)
#     for m in list(G.successors(Not(n))):
#         G.add_edge(ManualCalculation, m)
#     G.remove_node(n)
#     G.remove_node(Not(n))

# Include all nots as contradictions (again)
for n in list(G.nodes()):
    G.add_edge(n, Not(n), source='contradiction')

# I want identity edges (has to be last to avoid infinite loops)
for n in list(G.nodes()):
    G.add_edge(n, n, source="identity")

# Plot the changes
fig, (ax1, ax2) = plt.subplots(2, figsize=(10,10))
nx.draw_circular(g, with_labels=True, ax=ax1)
ax1.set_title("Before")
nx.draw_circular(G, with_labels=True, ax=ax2)
ax2.set_title("After")
plt.savefig('duality_before_after.pdf')

# Plot the result alone
fig, ax = plt.subplots(1, figsize=(10,10))
nx.draw_circular(G, with_labels=True, ax=ax)
plt.savefig('duality.pdf')

# Get the numpy array
def sortkey(self):
    if type(self) is Not:
        print("here", str(self), Not(self))
        # assert Not(self) != ManualCalculation4
        return str(Not(self))
    print("here2", str(self))
    return str(self)

nodelist = list(sorted(list(G.nodes()), key=sortkey))
outdf = pd.DataFrame([['' for i in range(len(nodelist))] for i in range(len(nodelist))],
                     columns=nodelist, index=nodelist)
for i, n in enumerate(nodelist):
    for j, m in enumerate(nodelist):
        if G.has_edge(n, m):
            if 'source' in G[n][m]:
                source = G[n][m]['source']
            else:
                source = 'calculated'
            outdf.iloc[i,j] = source

print(outdf)

# Sort by index of identity value
rowsort = []
for idx, row in outdf.iterrows():
    i = np.where(row=="identity")[0][0]
    rowsort.append(i)

outdf = outdf.iloc[rowsort, rowsort]

# Remove identity rows
drop_indexes = []
for idx, row in outdf.iterrows():
    if (row!='').sum() == 2:
        drop_indexes.append(idx)

# Remove identity columns
drop_columns = []
for column in outdf:
    col = outdf[column]
    if (col!='').sum() == 2:
        drop_columns.append(column)

outdf.drop(drop_indexes, axis='index', inplace=True)
outdf.drop(drop_columns, axis='columns', inplace=True)



print(outdf)
outdf.to_csv('duality.csv')