Noli, thanks for the confirmation.

I have attached my initial attempt (schedule.mod). As described in that
file:

/***
aim: to allocate patients to linear accelerators for treatment.
patients may require linacs with certain capabilities.
patients may prefer to avoid certain times of the day.
each patient needs to be scheduled on exactly one linac (with sufficient
capabilities)
a linac can only be used by up to one person at a time.
each linac only operates during a certain window during the day.
where not all patients' time preferences can be met, preference should be
given to those with more flexibility
***/


This runs, but I am not convinced I have specified everything correctly. In
particular, I notice that the 'badness' value is nonzero for patients even
where they have not provided a blackout window. ie: I believe no s.t.
constraint requires them to go above zero.
Could it be that I am seeing a nonzero "badness" because my "patient" param
is a number, rather than a reference to a specific member of a set?


For completeness, I have also attached a simple python file I whipped up
which wraps the schedule and allows me to dynamically inject data. (the
existing python libraries I saw are bitrotting, so I am wrapping the
commandline tool at least for now). If I increase the number of patients
(line 150 of parse.py) significantly, glpsol takes a long time and
ultimately thinks it cannot find a solution, even though I believe it
should be possible. My suspicion is that the "badness" calculation above is
close to the root cause.

For anyone who takes pity on me and has the time to help, I would be
interested in knowing whether my general approach is sound. I was not sure
how to best represent the fact that each patient needs exactly one session
on one machine: my approach of having a time on each machine (x) coupled
with a binary flag to indicate whether this is valid (lp) feels like a hack.

Thanks for taking the time to read this.

Nick.


On 28 March 2016 at 11:59, Noli Sicad <[email protected]> wrote:

> Hi Nick,
>
> Yes, this is the right place to ask questions about GLPK and MathProg.
>
> Noli
>
> On 3/28/16, Nick Farrell <[email protected]> wrote:
> > Is this the best place to ask for mathprog-specific queries?
> >
> > I am using glpsol to solve a problem, and feel like I am missing
> something.
> > I've looked through the examples and the wiki pages but haven't really
> > found what I want - which might simply be an introduction to mathprog. I
> > suspect I am misusing params where I should be using a set, resulting in
> a
> > more complex solution space.
>

Attachment: schedule.mod
Description: audio/mod

#!/usr/bin/env python3

from collections import defaultdict
import re
from pprint import pprint
import tempfile
from subprocess import check_call


reRecord = re.compile(r'''^\s*(\d+) ([^\[]+)\[([^\]]+)\][^\d]+([\d\.]+)''',
    re.MULTILINE)


class Glpsol:
    outfile = None

    def __init__(self, modfile):
        self.modfile = modfile
        self.sequence = []

    def add(self, data):
        self.sequence.append(data)

    def solve(self):
        datafile = self.render_data()
        self.outfile = "/tmp/o"
        check_call(["glpsol",
                    "--math", self.modfile,
                    "--data", datafile,
                    "--tmlim", "25",
                    "-o", self.outfile])

    def parse(self, variables_of_interest):
        result = defaultdict(dict)
        with open(self.outfile, "rt") as f:
            for match in reRecord.finditer(f.read()):
                # print(match.groups())
                var_name = match.group(2)
                parameter = match.group(3)
                value = match.group(4)
                if var_name in variables_of_interest:
                    result[var_name][parameter] = value
        return result

    def render_data(self):
        f = tempfile.NamedTemporaryFile(delete=False, mode="w+t")
        for s in self.sequence:
            s.render(f)
            f.write("\n")
        f.write("end;\n")
        f.close()
        return f.name


class Data:
    def render(self, destination):
        raise NotImplementedError("whoops")


class Param1D(Data):
    def __init__(self, param, mapping):
        self.param = param
        self.mapping = mapping

    def render(self, destination):
        destination.write("param {} :=\n".format(self.param))
        for k, v in self.mapping.items():
            destination.write("\t[{}]\t{}\n".format(k, v))
        destination.write(";")


class Set(Data):
    def __init__(self, name, values):
        self.name = name
        self.values = sorted(values)

    def render(self, destination):
        destination.write("set {} :=".format(self.name))
        for v in self.values:
            destination.write("\t{}".format(v))
        destination.write(";")


class Param0D(Data):
    def __init__(self, name, value):
        self.name = name
        self.value = value

    def render(self, destination):
        destination.write("param {} := {};\n".format(self.name, self.value))


class Param2D(Data):
    def __init__(self, name, rows, columns, default=None):
        self.default = default
        self.name = name
        self.rows = sorted(rows)
        self.columns = sorted(columns)
        self.values = dict()  # key is (r, c)

    def insert(self, row, column, value):
        assert row in self.rows
        assert column in self.columns
        self.values[row, column] = value

    def render(self, destination):
        destination.write("param {} : {} :=\n".format(self.name,
            "\t".join(self.columns)))
        for r in self.rows:
            destination.write("\t\t{}".format(r))
            for c in self.columns:
                if self.default is None:
                    destination.write("\t{}".format(self.values[r, c]))
                else:
                    destination.write("\t{}".format(
                        self.values.get((r, c), self.default)))
            destination.write("\n")
        destination.write(";")


class MultiParam(Data):
    def __init__(self, name, columns, data=None):
        self.name = name
        assert isinstance(columns, list)
        self.columns = columns
        self.data = []
        if data is not None:
            for d in data:
                assert len(d) == len(columns)
                self.data.append(d)

    def append(self, datum):
        assert len(datum) == len(self.columns)
        self.data.append(datum)

    def render(self, destination):
        destination.write("param : {} : {} :=\n".format(self.name,
            "\t".join(self.columns)))
        for i, datum in enumerate(self.data):
            destination.write("\t\t{}".format(i))
            for d in datum:
                destination.write("\t{}".format(d))

            destination.write("\n")
        destination.write(";")


if __name__ == '__main__':
    LINACS = ["linac{}".format(i) for i in range(2)]
    PATIENTS = ["MRN{:03d}".format(i) for i in range(10)]
    patient_index = range(1, len(PATIENTS) + 1)
    CAPABILITIES = {"IMRT", "VMAT", "MRI"}
    g = Glpsol("schedule.mod")
    g.add(Set(name="C", values=CAPABILITIES))
    g.add(Param0D(name="nPatients", value=len(PATIENTS)))
    g.add(Set(name="L", values=LINACS))

    c = Param2D("capabilities", LINACS, CAPABILITIES, default=0)
    for i, l in enumerate(LINACS):
        c.insert(l, "IMRT", 1)
        if i % 2 == 0:
            c.insert(l, "VMAT", 1)
        if i % 4 == 0:
            c.insert(l, "MRI", 1)
    g.add(c)

    r = Param2D("requirements", patient_index, CAPABILITIES, default=0)
    for p in patient_index:
        if p % 2 == 0:
            r.insert(p, "VMAT", 1)
        else:
            r.insert(p, "IMRT", 1)
    g.add(r)

    g.add(Param1D("firstTime", {l: 480 + 120 * i
                                for i, l in enumerate(LINACS)}))
    g.add(Param1D("lastTime", {l: 1200 for l in LINACS}))
    times = MultiParam("TIMES", ["patient", "startWindow", "endWindow"])
    for p in patient_index:
        if p % 2 == 0:
            times.append((p, 5 * 60, 12 * 60))  # not available in the morning
        if p % 3 == 0:
            times.append((p, 12 * 60, 22 * 60))  # not available after noon
        if p % 5 == 0:
            times.append((p, 8 * 60, 18 * 60))  # not available during work
    g.add(times)
    g.add(Param1D("duration", {p: 20 for p in patient_index}))

    g.solve()
    parsed = g.parse(["x"])
    pprint(parsed)
_______________________________________________
Help-glpk mailing list
[email protected]
https://lists.gnu.org/mailman/listinfo/help-glpk

Reply via email to