block_submodular.py

from gurobipy import *
from itertools import *
from functools import partial, reduce
import random
import sys
import math
from operator import mul


### Set and partition operations

# return all pairs
def pairs(V):
    return combinations(V,2)

# generate a iterable of sets
def powerset(iterable):
  xs = list(iterable)
  # note we return an iterator rather than a list
  return map(set,chain.from_iterable( combinations(xs,n) for n in range(len(xs)+1) ))


def algorithm_u(ns,m):
    def visit(n, a):
        ps = [[] for i in range(m)]
        for j in range(n):
            ps[a[j + 1]].append(ns[j])
        return ps

    def f(mu, nu, sigma, n, a):
        if mu == 2:
            yield visit(n, a)
        else:
            for v in f(mu - 1, nu - 1, (mu + sigma) % 2, n, a):
                yield v
        if nu == mu + 1:
            a[mu] = mu - 1
            yield visit(n, a)
            while a[nu] > 0:
                a[nu] = a[nu] - 1
                yield visit(n, a)
        elif nu > mu + 1:
            if (mu + sigma) % 2 == 1:
                a[nu - 1] = mu - 1
            else:
                a[mu] = mu - 1
            if (a[nu] + sigma) % 2 == 1:
                for v in b(mu, nu - 1, 0, n, a):
                    yield v
            else:
                for v in f(mu, nu - 1, 0, n, a):
                    yield v
            while a[nu] > 0:
                a[nu] = a[nu] - 1
                if (a[nu] + sigma) % 2 == 1:
                    for v in b(mu, nu - 1, 0, n, a):
                        yield v
                else:
                    for v in f(mu, nu - 1, 0, n, a):
                        yield v

    def b(mu, nu, sigma, n, a):
        if nu == mu + 1:
            while a[nu] < mu - 1:
                yield visit(n, a)
                a[nu] = a[nu] + 1
            yield visit(n, a)
            a[mu] = 0
        elif nu > mu + 1:
            if (a[nu] + sigma) % 2 == 1:
                for v in f(mu, nu - 1, 0, n, a):
                    yield v
            else:
                for v in b(mu, nu - 1, 0, n, a):
                    yield v
            while a[nu] < mu - 1:
                a[nu] = a[nu] + 1
                if (a[nu] + sigma) % 2 == 1:
                    for v in f(mu, nu - 1, 0, n, a):
                        yield v
                else:
                    for v in b(mu, nu - 1, 0, n, a):
                        yield v
            if (mu + sigma) % 2 == 1:
                a[nu - 1] = 0
            else:
                a[mu] = 0
        if mu == 2:
            yield visit(n, a)
        else:
            for v in b(mu - 1, nu - 1, (mu + sigma) % 2, n, a):
                yield v

    n = len(ns)
    a = [0] * (n + 1)
    for j in range(1, m + 1):
        a[n - m + j] = j - 1
    return f(m, n, 0, n, a)

def allPartitions(n,k):
    return [y for y in [frozenset(list(map(frozenset,x))) for x in algorithm_u(range(n),k)]]

def to_cycles(perm):
    pi = {i: perm[i] for i in range(len(perm))}
    cycles = []

    while pi:
        elem0 = next(iter(pi)) # arbitrary starting element
        this_elem = pi[elem0]
        next_item = pi[this_elem]

        cycle = []
        while True:
            cycle.append(this_elem)
            del pi[this_elem]
            this_elem = next_item
            if next_item in pi:
                next_item = pi[next_item]
            else:
                break
        if len(cycle)>1:
            cycles.append([c+1 for c in cycle])

    return cycles

# return all integer vectors of length k that sums to n
def vector_sums_to_n(k, n):
    if k == 1:
        yield (n,)
    else:
        for value in range(n + 1):
            for permutation in vector_sums_to_n(k - 1, n - value):
                yield (value,) + permutation

### Model

# given a partition, sum the variables
def partition_to_expr(partition, vars):
    return quicksum([vars[X] for X in partition])

# for a given partition X, add inequality of the form 
# expression(Y) - expression(X) >=0 for all X in dominated
def addMinPartitionInequalities(m, min_partition, partitions, vars):
    X_expr = partition_to_expr(min_partition, vars)
    for Y in partitions:
        if frozenset(Y) != frozenset(min_partition):
            Y_expr = partition_to_expr(Y, vars)

            expr = Y_expr - X_expr
            z = m.addConstr(expr >= 0)
            m.update()
            #print("expr>=0", expr)
    m.update()

# add all submodular inequalities over groundset U
def addSubmodularInequalities(m, U, var):
    # submodular inequalities
    allset = map(frozenset,powerset(U))
    for X in allset:
        for (a,b) in pairs(U-X):
            sa = frozenset([a])
            sb = frozenset([b])
            sab = frozenset([a,b])
            m.addConstr(var[X|sa]+var[X|sb]-var[X|sab]-var[X]>=0)

# test if the "expr >= constant" is a redundant constraint of m
def test_redundant(m, expr, constant):
    z = m.addConstr(expr >= constant-1)
    m.setObjective(expr, GRB.MINIMIZE)
    m.update()
    m.optimize()
    obj = m.getObjective()
    result = obj.getValue()
    m.remove(z)
    m.update()
    if(result<constant):
        return False
    return True

# name a particular variable
def name(x):
    return "e"+str(x)

# initialize the model, and create a variable for each subset of V
def initalize_model(V):
    m = Model("yay")
    m.setParam("LogToConsole",0)
    m.setParam("DualReductions",0)
    m.setParam("InfUnbdInfo",1)
    m.params.tuneResults = 5
    m.params.FeasibilityTol = 1e-9
    m.params.OptimalityTol = 1e-9
    # for set, add an variable
    variables = dict([])
    for X in powerset(V):
        variables[frozenset(X)] = m.addVar(name=name(X))
        m.update()
    return m, variables

# this is used to cut down certain results
# def symmetric_permutations():
#     Xs = frozenset([frozenset([0,1,2]),frozenset([3,4,5])])
#     Ys = frozenset([frozenset([0,3]),frozenset([1,4]),frozenset([2,5])])


def apply_permutation_to_partition(permutation, partition):
    return frozenset([frozenset([permutation[x] for x in X]) for X in partition])
def symmetric_permutations(k):
    n = k*(k-1)

    X = frozenset([frozenset(range(i*k,i*k+k)) for i in range(k-1)])
    Y = frozenset([frozenset([j*k+i for j in range(k-1)]) for i in range(k)])
    result = []
    for permutation in itertools.permutations(range(n)):
        if apply_permutation_to_partition(permutation, X)==X and apply_permutation_to_partition(permutation, Y)==Y:
            result.append(permutation)
    return result

#Z = symmetric_permutations(3)
#for X in Z:
#    print(to_cycles(list(X)))
#print(len(symmetric_permutations(3)))



# find all size constraint k partitions on n vertices
def all_size_constraint_k_partitions(n, sizes, Z):
    sizes = sorted(sizes)
    k = len(sizes)
    results = []
    for X in allPartitions(n,k):
        X_sizes = sorted([sum([Z[x] for x in S]) for S in X])
        if compare(X_sizes,sizes):
            results.append(X)
    return results

# We test if X_1 or X_2 is contained in some subset of Y
def partition_noncrossing(X, Y):
    for Xp in X:
        for Yp in Y:
            if Xp <= Yp:
                return True
    return False

def partition_nice(X, Y, interesting_2_partitions):
    if (frozenset([0,1,2,3]) in Y or 
       frozenset([0,1,2,4]) in Y or
       frozenset([0,1,2,5]) in Y or
       frozenset([3,4,5,0]) in Y or
       frozenset([3,4,5,1]) in Y or
       frozenset([3,4,5,2]) in Y):
        count = 0
        V = frozenset([0,1,2,3,4,5])
        for Z in Y:
            #print([V-Z,Z])
            if frozenset([V-Z,Z]) in interesting_2_partitions:
                #print("hey")
                count+=1
        if count >= 2:
            return True
    return False 


# takes all possible configurations, and a function takes a configuration, and 
# generates all interesting 2 and 3 partitions
def block_noncrossing_3_partition(configurations, all_interesting_2_partitions, all_interesting_3_partitions):
    n = 6
    V = range(n)
    V = set(V)

    # all partition of the blocks
    all_3_partitions = allPartitions(n,3)
    all_2_partitions = allPartitions(n,2)

    #print(all_2_partitions)

    X1 = frozenset([0,1,2])
    X2 = frozenset([3,4,5])

    Y1 = frozenset([0,3])
    Y2 = frozenset([1,4])
    Y3 = frozenset([2,5])

    min_interesting_2_partition = X = [X1,X2]
    min_interesting_3_partition = Y = [Y1,Y2,Y3]

    for Z in configurations:
        m, variables = initalize_model(V)

        interesting_2_partitions = all_interesting_2_partitions(Z)
        interesting_3_partitions = all_interesting_3_partitions(Z)

        # nice_interesting_2_partitions

        # super_interesting_2_partitions = []

        # for TT in interesting_2_partitions:
        #     zz = sorted([len(A) for A in TT])
        #     if zz == [1,5]:
        #         super_interesting_2_partitions.append(TT)

        # super_interesting_3_partitions = []

        # for Yprime in interesting_3_partitions:
        #     if partition_nice(X, Yprime, interesting_2_partitions):
        #         super_interesting_3_partitions.append(Yprime)

        addMinPartitionInequalities(m, X, interesting_2_partitions, variables)
        addMinPartitionInequalities(m, Y, interesting_3_partitions, variables)
        addSubmodularInequalities(m, V, variables)

        # presolve to eliminiate redundent constraints
        m.presolve()

        success = False
        # check if any interesting 3 partition Y' we have f(Y')<=f(Y)
        for Yprime in interesting_3_partitions:
            #if partition_nice(X, Yprime, interesting_2_partitions):
            if partition_noncrossing(X, Yprime): 
                # generate expression expr = 
                expr = partition_to_expr(Y,variables) - partition_to_expr(Yprime,variables)
                if test_redundant(m,expr,0):
                    success = True
                    break

        if not success:
            return False, Z
            break
    return True, None

# In order to use the program, all we need is to make sure
SP = symmetric_permutations(3)
def canonicalization(X):
    # consider all symmetries
    # that is, all permutations that maintains X[0],X[1],X[2] is the same
    z = sorted([[X[pi[i]] for i in range(n)] for pi in SP])
    return tuple(z[-1])

n = 6

sizes_2 = [4,4]
sizes_3 = [1,2,2]

sizes_specific_3 = [5,5,5]
vsize = 19

#print(SP)

# print("lol")
# INTERESTING = [1,1,2,1,1,1]
# p = (3,4,5,0,1,2)
# print([INTERESTING[p[i]] for i in range(n)])

# print("huh?")
# print(sorted([[INTERESTING[pi[i]] for i in range(n)] for pi in SP]))

# check if a component wise >= b
def compare(a,b):
    a = sorted(a)
    b = sorted(b)
    for i in range(len(a)):
        if a[i]<b[i]:
            return False
    return True

def configurations():
    # generate all configurations
    result = []
    for Y in vector_sums_to_n(n, vsize):
        X = list(range(len(Y)))
        for i in range(len(Y)):
            X[i] = min(Y[i],max(sizes_specific_3+sizes_3+sizes_2))
        X = tuple(X)
        if not compare([X[0]+X[1]+X[2], X[3]+X[4]+X[5]],sizes_2):
            continue
        if not compare([X[0]+X[3], X[1]+X[4], X[2]+X[5]],sizes_specific_3):
            continue
        #print(X,Y)
        if canonicalization(X)==X:
            result.append(X)
    return list(set(result))

interesting_2_partitions = partial(all_size_constraint_k_partitions,n,sizes_2)
interesting_3_partitions = partial(all_size_constraint_k_partitions,n,sizes_3)

#print(configurations())
print(len(configurations()))
print(block_noncrossing_3_partition(configurations(),interesting_2_partitions, interesting_3_partitions))