2

I will be grateful to anyone who can help me write some Python code to enumerate the 21×2×3 arrays, indexed with i, j and k, which are two-thirds filled with 0's and one-third filled with the values 'Ava', 'Bob', 'Joe', 'Mia', 'Sam', 'Tom', 'Zoe' in such a way that:

  1. fixed the index i you have exactly two empty 2-tuples and one 2-tuple with different non-zero values;

  2. fixed the index k you have exactly fourteen empty 2-tuples and seven 2-tuple with different non-zero values;

  3. fixed the indexes j and k you have a 21-tuple with fourteen zero values and exactly one occurrence of each of the non-zero values, respecting the following constraints:

    a) "Ava" can appear only in a row with index 0, 1, 4, 6, 10, 11, 13, 14, 15, 19 or 20;

    b) "Bob" can appear only in a row with index 2, 3, 5, 7, 8, 9, 12, 16, 17 or 18;

    c) "Joe" can appear only in a row with index 2, 4, 5, 7, 8, 10, 14, 15, 18 or 20;

    d) "Mia" can appear only in a row with index 0, 1, 3, 6, 9, 12, 13, 16, 17 or 19;

    e) "Sam" can appear only in a row with index 1, 2, 7, 9, 15, 17 or 20;

    f) "Tom" can appear only in a row with index 0, 3, 8, 10, 12, 16 or 19;

    g) "Zoe" can appear only in a row with index 4, 5, 6, 11, 13, 14 or 18.

As a result I would like to obtain something like this:

[ 0   0       [Tom Mia      [ 0   0
  0   0        Ava Sam        0   0
  0   0        Sam Bob        0   0
  0   0        Bob Tom        0   0
  0   0         0   0        Joe Zoe
  0   0        Joe Zoe        0   0
  0   0         0   0        Zoe Ava
 Joe Sam        0   0         0   0
  0   0         0   0        Tom Bob
  0   0         0   0        Mia Sam
 Tom Ava        0   0         0   0
 Ava Zoe        0   0         0   0
 Bob Mia        0   0         0   0
  0   0        Mia Ava        0   0
  0   0        Zoe Joe        0   0
 Sam Joe        0   0         0   0
  0   0         0   0        Bob Tom
  0   0         0   0        Sam Mia
 Zoe Bob        0   0         0   0
 Mia Tom        0   0         0   0
  0   0 ]       0   0 ]      Ava Joe]

Rows represent school classes, columns represent school terms (there are 2 of them), tubes represent class days (there are 3 of them: Monday, Wednesday and Friday). So the first horizontal slice of the above solution means that class 1A has lesson only on Wednesday, in the the first term with teacher Tom and in the second term with teacher Mia. (Teachers can only work in some classes and not in others.)

Thanks in advance!


Update n. 1

As a starting point, I tried to attack the following toy problem:
Enumerate all arrays with a given number of rows and 3 columns which are two-thirds filled with "0" and one-third filled with "1" in such a way that summing the values in each row you always get 1 and summing the values in each column you always get rows / 3.
Finally, after struggling a bit, I think I managed to get a solution with the following code, that I kindly ask you to correct or improve. (I have set rows = 6 because the number of permutations of the obvious solution is 6!/(2!*2!*2!) = 90, whereas setting rows = 21 I would have got 21!/(7!*7!*7!) = 399,072,960 solutions.)

from ortools.sat.python import cp_model

# Create model.
model = cp_model.CpModel()

# Create variables.
rows = 6
columns = 3
x = []
for i in range(rows):
  x.append([model.NewBoolVar(f'x[{i}][{j}]') for j in range(columns)])

# Add constraints.
for i in range(rows):
  model.Add(sum(x[i]) == 1)

# Uncomment the following four lines of code if you want to solve the slightly more general problem that asks to enumerate
# all boolean arrays, with a given number of rows and columns, filled in such a way that summing the values in each
# row you always get 1 and summing the values in each column you always get no more than the ceiling of (rows / columns).
# if rows % columns != 0:
#   for j in range(columns):
#     model.Add(sum(x[i][j] for i in range(rows)) <= rows // columns + 1)
# else:
  for j in range(columns):
    model.Add(sum(x[i][j] for i in range(rows)) == rows // columns)


class MyPrintedSolution():

    def __init__(self, sol, sol_number):
        self.sol = sol
        self.sol_number = sol_number

    def PrintReadableTable(self):
        print(f'Solution {self.sol_number}, printed in readable form:')
        counter = 0
        for v in self.sol:
          if counter % columns != columns-1:
            print(v, end = ' ')
          else:
            print(v)
          counter += 1
        print()

    def PrintRawSolution(self):
        print(f'Solution {self.sol_number}, printed in raw form:')
        counter = 0
        for v in self.sol:
          print(f'{v}', end = '')
          counter += 1
        print('\n')


class VarArraySolutionPrinter(cp_model.CpSolverSolutionCallback):

    def __init__(self, variables, limit):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__variables = variables
        self.__solution_count = 0
        self.__solution_limit = limit

    def solution_count(self):
        return self.__solution_count

    def on_solution_callback(self):
        self.__solution_count += 1
        solution = [self.Value(v) for v in self.__variables]
        myprint = MyPrintedSolution(solution, self.__solution_count)
        myprint.PrintReadableTable()
        # myprint.PrintRawSolution()
        if self.__solution_count >= self.__solution_limit:
          print(f'Stop search after {self.__solution_limit} solutions')
          self.StopSearch()


# Create solver and solve model.
solver = cp_model.CpSolver()
# solver.parameters.num_workers = 16  # Solver works better with more workers. (At least 8, 16 if enough cores.)
# solver.parameters.log_search_progress = True
solver.parameters.enumerate_all_solutions = True
# solver.parameters.max_time_in_seconds = 10.0
solution_limit = 100000
solution_printer = VarArraySolutionPrinter([x[i][j] for i in range(rows) for j in range(columns)], solution_limit)
solver.Solve(model, solution_printer)

Update n. 2

Following @Christopher Hamkins' initial roadmap and subsequent precious suggestions, I think I finally got what I wanted, using the following code (although I am of course always open to corrections or further suggestions).

from ortools.sat.python import cp_model

# Create model.
model = cp_model.CpModel()

# Create variables.
classes = 21 # indexed with "i", but one could as well have chosen "c"
terms = 2 # indexed with "j", but one could as well have chosen "t"
days = 3 # indexed with "k", but one could as well have chosen "d"
persons = 8 # indexed with "p"
persons_names = [' 0 ', 'Ava', 'Bob', 'Joe', 'Mia', 'Sam', 'Tom', 'Zoe']
classes_names = ['1A', '1B', '1C', '1D', '1E', '1F', '1G', '2A', '2B', '2C', '2D', '2E', '2F', '2G', '3A', '3B', '3C', '3D', '3E', '3F', '3G']
classes_p = [[] for _ in range(persons)]
classes_p[0] = list(range(classes))
classes_p[1] = [0, 1, 4, 6, 10, 11, 13, 14, 15, 19, 20] # list of classes in which person 1 can work
classes_p[2] = [2, 3, 5, 7, 8, 9, 12, 16, 17, 18] # list of classes in which person 2 can work
classes_p[3] = [2, 4, 5, 7, 8, 10, 14, 15, 18, 20] # list of classes in which person 3 can work
classes_p[4] = [0, 1, 3, 6, 9, 12, 13, 16, 17, 19] # list of classes in which person 4 can work
classes_p[5] = [1, 2, 7, 9, 15, 17, 20] # list of classes in which person 5 can work
classes_p[6] = [0, 3, 8, 10, 12, 16, 19] # list of classes in which person 6 can work
classes_p[7] = [4, 5, 6, 11, 13, 14, 18] # list of classes in which person 7 can work
x = {}
for i in range(classes):
  for j in range(terms):
    for k in range(days):
      for p in range(persons):
        x[i, j, k, p] = model.NewBoolVar(f'x[{i}, {j}, {k}, {p}]')

# Add constraints.
"""
For all i, j, k constrain the sum of x[i, j, k, p] over p in the range of people to be equal to 1,
so exactly nobody or one person is selected at a given slot.
"""
for i in range(classes):
  for j in range(terms):
    for k in range(days):
      model.Add(sum(x[i, j, k, p] for p in range(persons)) == 1)

"""
For all i constrain the sum of x[i, j, k, p] over all j, k, p in their respective ranges (except p = 0)
to be exactly equal to 2, so exactly two people are in a given row.
"""
for i in range(classes):
  model.Add(sum(x[i, j, k, p] for j in range(terms) for k in range(days) for p in range(1, persons)) == 2)

"""
For all i, k, and for p = 0, add the implications
x[i, 0, k, 0] == x[i, 1, k, 0]
"""
for i in range(classes):
 for k in range(days):
    model.Add(x[i, 0, k, 0] == x[i, 1, k, 0])

"""
For all i, p (except p = 0), constrain the sum of x[i, j, k, p] over all j and k
to be at most 1.
"""
for i in range(classes):
  for p in range(1, persons):
    model.Add(sum(x[i, j, k, p] for j in range(terms) for k in range(days)) <= 1)
    # for k in range(days): # Equivalent alternative to the previous line of code
    #   model.AddBoolOr([x[i, 0, k, p].Not(), x[i, 1, k, p].Not])

"""
For all j, k constrain the sum of x[i, j, k, p] over all i, p in their respective ranges (except p = 0)
to be exactly equal to 7, so exactly seven people are in a given column.
"""
for j in range(terms):
  for k in range(days):
    model.Add(sum(x[i, j, k, p] for i in range(classes) for p in range(1, persons)) == 7)

"""
For all j, k, p (except p = 0) constrain the sum of x[i, j, k, p] over all i
to be exactly equal to 1, so each person appears exactly once in the column.
"""
for j in range(terms):
  for k in range(days):
    for p in range(1, persons):
      model.Add(sum(x[i, j, k, p] for i in range(classes)) == 1)

"""
For all j and k, constrain x[i, j, k, p] == 0 for the row i in which each person p can't appear.
"""
for p in range(persons):
  for i in enumerate(set(range(classes)) - set(classes_p[p])):
    for j in range(terms):
      for k in range(days):
        model.Add(x[i[1], j, k, p] == 0)


class MyPrintedSolution():

    def __init__(self, sol, sol_number):
        self.sol = sol
        self.sol_number = sol_number

    def PrintReadableTable1(self):
        print(f'Solution {self.sol_number}, printed in first readable form:')
        print('    |      Mon      |      Wed      |      Fri     ')
        print(' Cl | Term1   Term2 | Term1   Term2 | Term1   Term2') 
        print('----------------------------------------------------', end='')
        q = [_ for _ in range(8)] + [_ for _ in range(24, 32)] + [_ for _ in range(8, 16)] + [_ for _ in range(32, 40)] + [_ for _ in range(16, 24)] + [_ for _ in range(40, 48)]
        r = []
        for i in range(21):
            r += [n+48*i for n in q]
        shuffled_sol = [self.sol[m] for m in tuple(r)]
        counter = 0
        for w in shuffled_sol:
          if (counter % (persons * days * terms)) == 0:
            print('\n ', classes_names[counter // (terms * days * persons)], sep='', end=' |')
          if w:
            print('  ', persons_names[counter % persons], sep='', end='   ')
          counter += 1
        print('\n')

    def PrintReadableTable2(self):
        print(f'Solution {self.sol_number}, printed in second readable form:')
        print(' Cl |      Term1             Term2       ')
        print(' Cl | Mon   Wed   Fri   Mon   Wed   Fri  ')
        print('----------------------------------------', end = '')
        counter = 0
        for v in self.sol:
          if (counter % (persons * days * terms)) == 0:
            print('\n ', classes_names[counter // (terms * days * persons)], sep = '', end = ' |')
          if v:
            print(' ', persons_names[counter % persons], sep = '', end = '  ')
          counter += 1
        print('\n')

    def PrintRawSolution(self):
        print(f'Solution {self.sol_number}, printed in raw form:')
        counter = 0
        for v in self.sol:
          print(f'{v}', end = '')
          counter += 1
        print('\n')


class VarArraySolutionPrinter(cp_model.CpSolverSolutionCallback):

    def __init__(self, variables, limit):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__variables = variables
        self.__solution_count = 0
        self.__solution_limit = limit

    def solution_count(self):
        return self.__solution_count

    def on_solution_callback(self):
        self.__solution_count += 1
        solution = [self.Value(v) for v in self.__variables]
        myprint = MyPrintedSolution(solution, self.__solution_count)
        myprint.PrintReadableTable1()
        # myprint.PrintReadableTable2()
        # myprint.PrintRawSolution()
        if self.__solution_count >= self.__solution_limit:
          print(f'Stop search after {self.__solution_limit} solutions')
          self.StopSearch()


# Create solver and solve model.
solver = cp_model.CpSolver()
# solver.parameters.num_workers = 16  # Solver works better with more workers. (At least 8, 16 if enough cores.)
# solver.parameters.log_search_progress = True
solver.parameters.enumerate_all_solutions = True
# solver.parameters.max_time_in_seconds = 10.0
solution_limit = 20
solution_printer = VarArraySolutionPrinter([x[i, j, k, p] for i in range(classes) for j in range(terms) for k in range(days) for p in range(persons)], solution_limit)
status = solver.Solve(model, solution_printer)

Update n. 3

@AirSquid proposed a solution using PuLP which is to me almost as valuable as the one using CP-SAT. It provides only one solution at a time, but (it has other advantages and) one can always get around this by adding some ad hoc further constraints, for example to see a different solution with a certain person in a specific position.

S.R.
  • 23
  • 4
  • Can you show what you tried? If you already know about Pulp, your easier example is rather trivial to *model* because it just consists of sums over some matrix-view of variables. Somewhat like the linear-assignment problem, but a bit more general. What's not working out with this approach? | The full problem (if i interpret it correctly) can be tackled in the same way plus some complications. I think i would prefer or-tools cp-sat solver here to Pulp. – sascha Dec 26 '22 at 17:21
  • @sascha I've just discovered the existence of PuLP via Google, but I know almost nothing about it or linear programming (and very little about Python, too). I tried this code: `import pulp` `vars = pulp.LpVariable.matrix('vars', (range(21), range(3)), 0, 1, pulp.LpInteger)` `prob = pulp.LpProblem('array_problem', pulp.LpMinimize)` `for i in range(21):` `prob += sum(vars[i]) == 1` `for j in range(3):` `prob += sum(vars[:,j]) == 7` `prob += sum(vars)` `prob.solve()` `result = [[vars[i][j].value() for j in range(3)] for i in range(21)]` `print(result)` but I got a compilation error. – S.R. Dec 26 '22 at 20:13
  • I'm afraid, you need to get more familiar with python and those discrete-optimization tools first. The task at hand will ask for some *technical* stuff (due to some aggregations / hierarchies needed when modelling it). Some remarks: 1) Putting code into comments isn't helpful 2) The general problem is an integer-program, not an LP (pulp handles both). 3) You style of indexing into the matrix looks wrong to me but i don't use pulp. Check what pulp expects when indexing one dim of a matrix and compare with your `sum(vars[i]) == 1` 4) Python isn't compiled. You probably see a runtime-error. – sascha Dec 26 '22 at 20:15
  • I'm wondering where `pulp.LpVariable.matrix` is coming from. That's not supported according to the [docs](https://coin-or.github.io/pulp/technical/pulp.html#pulp.LpVariable). – sascha Dec 26 '22 at 20:21
  • You mentioned index i, j and k for the array but later in the text you refer to z. Did you mean k? This problem would be a great match for Google's OR-Tools CP-Sat solver. – Christopher Hamkins Dec 27 '22 at 15:01
  • @ChristopherHamkins Yes, I meant "k", thanks for pointing out the error! And thank you also for your suggestion, in agreement with sascha's one. (If you could suggest me at least some more detailed roadmap to follow it would be amazing, but I am already grateful to both you and sascha.) – S.R. Dec 27 '22 at 15:41
  • I think your reformulation to get 7 people in a column is correct. Do you get solutions when you remove the constraints about certain people not appearing in certain rows? – Christopher Hamkins Dec 30 '22 at 15:46
  • I think I made a mistake with the formulation of the implications x[i, 0, k, p] implies x[i, 1, k, p].Not and x[i, 1, k, p] implies x[i, 0, k, p].Not This ought to do the trick, however: ``model.AddBoolOr([x[i, 0, k, p].Not(), x[i, 1, k, p].Not])`` (just this line instead of the current two). If it works I'll edit my answer. – Christopher Hamkins Dec 30 '22 at 16:13
  • The constraint "For all k, p (except p = 0) constrain the sum of x[i, j, k, p] over all i and j to be exactly equal to 1, so each person appears exactly once in the column." should also read "For all j, k, and p (except p = 0) constrain the sum of x[i, j, k, p] over all i to be exactly equal to 1, so each person appears exactly once in the column." I had not considered that a person appears once on the left and once on the right. – Christopher Hamkins Dec 30 '22 at 17:08
  • I edited my answer to reflect the last three comments. – Christopher Hamkins Jan 02 '23 at 15:51
  • Thanks again, @Christopher Hamkins, your help has been invaluable to me! I am sorry it took me so long to reply, but I was sucked into a vortex of Python and CP-SAT concepts I didn't know enough about, and I had to delve deeper to understand a little what was going on under the hood and solve some problems encountered along the way. – S.R. Jan 04 '23 at 23:46
  • Regarding your suggestions, here I would like to point out only that: 1) the mistake with logical implications was an oversight on my part (I know very well that 'a → b' is equivalent '¬a ∨ b'), sorry! 2) I think a probably simpler alternative is the following: "For all i and p (except p = 0), constrain the sum of x[i, j, k, p] for all j and k to be to be at most 1" (not "exactly equal to 1", as you stated in your original answer). – S.R. Jan 04 '23 at 23:46

2 Answers2

1

Your "toy" problem is definitely going in the right direction.

For your actual problem, try making a 21×2×3x8 array x indexed with i, j, k and p (for person) of BoolVar's. The last index represents the person, it will need 0 to represent "nobody" and for the rest Ava = 1, Bob = 2, etc., so its max value will be one more than the number of people. If the variable X[i,j,k,p] is true (1) it means that the given person p is present at the index i, j, k. If X[i,j,k,0] is true, it means a 0 = nobody is present at i, j, k.

For all i, j, k, constrain the sum of x[i, j, k, p] for p in the range of people to be equal to 1, so exactly nobody or one person is selected at a given slot.

For point 1: fixed the index i you have exactly two empty 2-tuples and one 2-tuple with different non-zero values:

For all i constrain the sum of x[i, j, k, p] for all j, k, p in their respective ranges (except p = 0) to be exactly equal to 2, so exactly two people are in a given row.

For all i, k, and for p = 0, add the implications

x[i, 0, k, 0] == x[i, 1, k, 0]

This will ensure that if one of the pair is 0, so is the other.

For all i, k and p except p = 0, add the implications

x[i, 0, k, p] implies x[i, 1, k, p].Not and

x[i, 1, k, p] implies x[i, 0, k, p].Not

(Actually one of these alone should be sufficient)

You can directly add an implication with the AddImplication(self, a, b) method, or you can realize that "a implies b" means the same thing as "b or not a" and add the implication with the AddBoolOr method. For the first implication, with x[i, 0, k, p] as a, and x[i, 1, k, p].Not as b, therefore adding:

AddBoolOr([x[i, 0, k, p].Not(), x[i, 1, k, p].Not])

Note that both variables are negated with Not in the expression.

Since the other implication assigns x[i, 1, k, p] as a, and x[i, 0, k, p].Not as b, the resulting expression is exactly the same

AddBoolOr([x[i, 0, k, p].Not(), x[i, 1, k, p].Not])

so you only need to add it once.

This will ensure a tuple will consist of two different people.

Alternative formulation of the last part:

For all i and p except p = 0, constraint the sum of x[i, j, k, p] for all j and k to be exactly equal to 1.

For point 2: fixed the index k you have exactly fourteen empty 2-tuples and seven 2-tuple with different non-zero values;

For all j and k constrain the sum of x[i, j, k, p] for all i and p (except p=0) in their respective ranges to be exactly equal to 7, so exactly seven people are in a given column.

For all j, k, and p (except p = 0) constrain the sum of x[i, j, k, p] over all i to be exactly equal to 1, so each person appears exactly once in the column (that is, once for each value of the indices j and k, for some value of i).

For point 3:

For all j and k, Constrain x[i, j, k, p] == 0 for the row i in which each person p can't appear.

Let us know how it works.

Christopher Hamkins
  • 1,442
  • 9
  • 18
1

You're taking a pretty big swing if you are new to the trifecta of python, linear programming, and pulp, but the problem you describe is very doable...perhaps the below will get you started. It is a smaller example that should work just fine for the data you have, I just didn't type it all in.

A couple notes:

  • The below is a linear program. It is "naturally integer" as coded, preventing the need to restrict the domain of the variables to integers, so it is much easier to solve. (A topic for you to research, perhaps).
  • You could certainly code this up as a constraint problem as well, I'm just not as familiar. You could also code this up like a matrix as you are doing with i, j, k, but most frameworks allow more readable names for the sets.
  • The teaching day M/W/F is arbitrary and not linked to anything else in the problem, so you can (externally to the problem), just pick 1/3 of the assignments per day from the solution for each course & term.
  • The transition from the verbiage to the constraint formulation is most of the magic in linear programming and you'd be well suited with an introductory text if you continue along!

Code:

# teacher assignment

import pulp
from itertools import chain

# some data...
teach_days = {'M', 'W', 'F'}
terms = {'Spring', 'Fall'}
courses = {'Math 101', 'English 203', 'Physics 201'}
legal_asmts = { 'Bob': {'Math 101', 'Physics 201'},
                'Ann': {'Math 101', 'English 203'},
                'Tim': {'English 203'},
                'Joe': {'Physics 201'}}

# quick sanity check
assert courses == set.union(*chain(legal_asmts.values())), 'course mismatch'

# set up the problem

prob = pulp.LpProblem('teacher_assignment', pulp.LpMaximize)

# make a 3-tuple index of the term, class, teacher
idx = [(term, course, teacher) for term in terms for course in courses for teacher in legal_asmts.keys()]

assn = pulp.LpVariable.dicts('assign', idx, cat=pulp.LpContinuous, lowBound=0)

# OBJECTIVE:  teach as many courses as possible within constraints...
prob += pulp.lpSum(assn)

# CONSTRAINTS
# teach each class no more than once per term
for term in terms:
    for course in courses:
        prob += pulp.lpSum(assn[term, course, teacher] for teacher in legal_asmts.keys()) <= 1

# each teacher no more than 1 course per term
for term in terms:
    for teacher in legal_asmts.keys():
        prob += pulp.lpSum(assn[term, course, teacher] for course in courses) <= 1

# each teacher can only teach within legal assmts, and if legal, only teach it once
for teacher in legal_asmts.keys():
    for course in courses:
        if course in legal_asmts.get(teacher):
            prob += pulp.lpSum(assn[term, course, teacher] for term in terms) <= 1
        else:  # it is not legal assignment
            prob += pulp.lpSum(assn[term, course, teacher] for term in terms) <= 0


prob.solve()

#print(prob)

# Inspect results...
for i in idx:
    if assn[i].varValue:  # will be true if value is non-zero
        print(i, assn[i].varValue)

Output:

Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Presolve 16 (-10) rows, 12 (-12) columns and 32 (-40) elements
Perturbing problem by 0.001% of 1 - largest nonzero change 0.00010234913 ( 0.010234913%) - largest zero change 0
0  Obj -0 Dual inf 11.99913 (12)
10  Obj 5.9995988
Optimal - objective value 6
After Postsolve, objective 6, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 6 - 10 iterations time 0.002, Presolve 0.00
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

('Spring', 'Math 101', 'Bob') 1.0
('Spring', 'Physics 201', 'Joe') 1.0
('Spring', 'English 203', 'Ann') 1.0
('Fall', 'Math 101', 'Ann') 1.0
('Fall', 'Physics 201', 'Bob') 1.0
('Fall', 'English 203', 'Tim') 1.0

EDIT: corrected form of problem

Misunderstood part of the problem statement. The below is fixed. Needed to introduce a binary indicator variable for the day assignment per form and had some fun with tabulate.

Using an LP has the advantage that (with the included obj statement) it will do the best possible within the constraints to teach as much as possible, even if there is a teacher shortage, where a CP will not. A CP on the other hand can enumerate all the combos that satisfy the constraints, the LP cannot.

Code

# teacher assignment

import pulp
from tabulate import tabulate

# some data...
teach_days = {'M', 'W', 'F'}
terms = {'Spring', 'Fall'}
forms = list(range(20))
teach_capable = {   "Ava" : [ 0, 1, 4, 6, 10, 11, 13, 14, 15, 19, 20],
                    "Bob" : [ 2, 3, 5, 7, 8, 9, 12, 16, 17, 18],
                    "Joe" : [ 2, 4, 5, 7, 8, 10, 14, 15, 18, 20],
                    "Mia" : [ 0, 1, 3, 6, 9, 12, 13, 16, 17, 19],
                    "Sam" : [ 1, 2, 7, 9, 15, 17, 20],
                    "Tom" : [ 0, 3, 8, 10, 12, 16, 19],
                    "Zoe" : [ 4, 5, 6, 11, 13, 14, 18],}


# set up the problem

prob = pulp.LpProblem('teacher_assignment', pulp.LpMaximize)

# make a 4-tuple index of the day, term, class, teacher
idx = [(day, term, form, teacher) 
        for day in teach_days
        for term in terms
        for form in forms 
        for teacher in teach_capable.keys()]

# variables
assn = pulp.LpVariable.dicts('assign', idx, cat=pulp.LpContinuous, lowBound=0)
form_day = pulp.LpVariable.dicts('form-day', 
    [(form, day) for form in forms for day in teach_days],
    cat=pulp.LpBinary)   # inidicator for which day the form uses

# OBJECTIVE:  teach as many courses as possible within constraints...
prob += pulp.lpSum(assn)

# CONSTRAINTS
# 1.  Teach each form on no more than 1 day
for form in forms:
    prob += pulp.lpSum(form_day[form, day] for day in teach_days) <= 1  # limit to 1 day per form
for form in forms:
    for day in teach_days:
        for term in terms:
            # no more than 1 assignment, if this day is the designated "form-day"
            prob += pulp.lpSum(assn[day, term, form, teacher] for teacher in teach_capable.keys()) \
                <= form_day[form, day]


# 2.  Each teacher can only teach within legal assmts, and limit them to teaching that form once
for teacher in teach_capable.keys():
    for form in forms:
        if form in teach_capable.get(teacher):
            prob += pulp.lpSum(assn[day, term, form, teacher] for day in teach_days for term in terms) <= 1
        else:  # it is not legal assignment
            prob += pulp.lpSum(assn[day, term, form, teacher] for day in teach_days for term in terms) <= 0

# 3.  Each teacher can only teach on once per day per term
for teacher in teach_capable.keys():
    for term in terms:
        for day in teach_days:
            prob += pulp.lpSum(assn[day, term, form, teacher] for form in forms) <= 1

prob.solve()
print("Status = %s" % pulp.LpStatus[prob.status])

#print(prob)

# gather results...
selections = []
for i in idx:
    if assn[i].varValue:  # will be true if value is non-zero
        selections.append(i)
        #print(i, assn[i].varValue)


# Let's try to make some rows for tabulate...  hacky but fun
def row_index(label):
    """return the form, column index, and name"""
    col = 1
    if 'W' in label: col += 2
    elif 'F' in label: col += 4
    if 'Fall' in label: col += 1
    return label[2], col, label[-1]

headers = ['Form', 'Mon-1', 'Mon-2', 'Wed-1', 'Wed-2', 'Fri-1', 'Fri-2']
row_data = [[f,'','','','','',''] for f in forms]
for selection in selections:
    form, col, name = row_index(selection)
    row_data[form][col] = name

print(tabulate(row_data, headers=headers, tablefmt='grid'))

Output:

Status = Optimal
+--------+---------+---------+---------+---------+---------+---------+
|   Form | Mon-1   | Mon-2   | Wed-1   | Wed-2   | Fri-1   | Fri-2   |
+========+=========+=========+=========+=========+=========+=========+
|      0 |         |         | Ava     | Tom     |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|      1 | Mia     | Sam     |         |         |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|      2 |         |         |         |         | Sam     | Joe     |
+--------+---------+---------+---------+---------+---------+---------+
|      3 |         |         |         |         | Bob     | Mia     |
+--------+---------+---------+---------+---------+---------+---------+
|      4 | Ava     | Zoe     |         |         |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|      5 | Zoe     | Joe     |         |         |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|      6 |         |         | Mia     | Zoe     |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|      7 |         |         | Joe     | Sam     |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|      8 | Bob     | Tom     |         |         |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|      9 |         |         | Sam     | Mia     |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|     10 |         |         |         |         | Ava     | Tom     |
+--------+---------+---------+---------+---------+---------+---------+
|     11 |         |         | Zoe     | Ava     |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|     12 |         |         | Tom     | Bob     |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|     13 |         |         |         |         | Zoe     | Ava     |
+--------+---------+---------+---------+---------+---------+---------+
|     14 |         |         |         |         | Joe     | Zoe     |
+--------+---------+---------+---------+---------+---------+---------+
|     15 | Joe     | Ava     |         |         |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|     16 |         |         |         |         | Tom     | Bob     |
+--------+---------+---------+---------+---------+---------+---------+
|     17 | Sam     | Bob     |         |         |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|     18 |         |         | Bob     | Joe     |         |         |
+--------+---------+---------+---------+---------+---------+---------+
|     19 | Tom     | Mia     |         |         |         |         |
+--------+---------+---------+---------+---------+---------+---------+
[Finished in 165ms]
AirSquid
  • 10,214
  • 2
  • 7
  • 31
  • Thanks, AirSquid, but I think you have solved a different problem (probably because I am a non-native speaker and I wasn't clear enough). With 'class' here I mean a group of students, not a course, as it is common in some non-US school systems. Maybe we can call it a 'form' to avoid confusion. Anyway, I tried to adapt your code, but I think I'm stuck when I have to add the constraint that each form (group of students) has lessons on only one day of the week for all the year (the same day for both terms), with exactly one teacher per term, one for the first term and another for the second term. – S.R. Jan 04 '23 at 23:47
  • Got it. It isn't too daunting of a modification, and it would look pretty similar to the solution you are working with now in the CP setup.... and it looks like the issue is solved... LMK if it is needed. – AirSquid Jan 05 '23 at 18:58
  • Although the problem is now solved in the CP setting, I'd love to see solutions in other settings as well: it's always good to have many strings to one's bow! I admit that in the next few days I won't have much time to do any further reading, but I'd really appreciate if, at least, you could explain to me why the following naive constrain formulation doesn't seem to work: `for f in forms:` `for d in days:` `prob += (pulp.lpSum(assn[d, term, f, t] for term in terms for t in legal_asmts.keys()) == 0 or pulp.lpSum(assn[d, term, f, t] for t in terms for t in legal_asmts.keys()) == 2)`. Thx again! – S.R. Jan 05 '23 at 20:31
  • See my edit above with revisions per your clarification on the problem statement. Regarding your "constraint" ... 'or' statements are not legal inside a constraint. I used a new binary variable to handle the day assignments.... see code! ;) – AirSquid Jan 05 '23 at 23:39
  • Thank you very much, @AirSquid, you have been very kind to me! I love having a working solution in the PuLP setting as well, the code you provided is crystal-clear and your comments are insightful (not even to mention that you let me discover _tabulate_ as a good byproduct). At the moment I'm not allowed to upvote your answer, but if I ever get the privilege I will certainly do it right away, because to me it's as valuable as Christopher Hamkins' one. (I just point out two typos, if you want to correct them: `range(20)` should be everywhere `range(21)`, and "inidicator" of course "indicator".) – S.R. Jan 06 '23 at 12:42
  • All good. Interesting problem. And tabulate is kinda neat. ;) And yes, you are correct on the edits! – AirSquid Jan 06 '23 at 15:12