# Using MILP to place a set of primers along a genome

Define variables $$p_i,u_i\in\{0,1\}^G$$, for $$i=1,\ldots,8$$ and $$G=30000$$.

Let $$v$$ be a constant vector also in $$\{0,1\}^G$$, with approximately 25% of its entries equal to $$1$$ (randomly located).

Let $$c$$ be the vector $$c=[1,2,3,\ldots,G]$$.

Let $$A$$ be the matrix with $$1$$'s along its diagonal and $$-1$$'s along its subdiagonal (with zeros everywhere else). Note that I'm pretty sure $$A$$ is totally unimodular.

Let $$lb_u,ub_u,lb_p,ub_p$$ be reasonably chosen nonnegative integers, e.g.

$$lb_u=18$$ $$ub_u=60$$ $$lb_p=18$$ $$ub_p=24$$

I would like to solve the following BILP,

\begin{align} \underset{p}{\text{min}}&\;\;\;v^T\sum_i p_i\\\\ \text{subject to}&\;\;\;\sum_ip_i\leq{\bf 1}\\\\ &\;\;\;Ap_i\leq u_i\\\\ &\;\;\;{\bf 1}^Tu_i=1\\\\ &\;\;\;lb_u\leq c^T(u_{i+1}-u_i)\leq ub_u\\\\ &\;\;\;lb_p\leq{\bf 1}^Tp_i\leq ub_p \end{align}

To give the reader some context, I have a set of eight primers, $$p_1,\ldots,p_8$$ that I need to place along a genome such that their overlap with a mutation indicator vector $$v$$ is as small as possible.

The first constraint enforces primer non-overlap. The second and third constraints enforce primer contiguity. The fourth constraint dictates that adjacent primers must have a minimum and maximum distance between their respective starting positions of $$lb_u$$ and $$ub_u$$. The fifth constraint requires that the length of each primer must be between $$lb_p$$ and $$ub_p$$.

The size of the genome $$G$$ is approximately 30,000, which means that this is an integer programming problem with 480,000 binary variables.

Empirical evidence suggests that I can obtain a high quality relaxation by performing the standard LP relaxation on just the $$p_i$$ variables. This requires introducing the additional nonnegativity constraints $$p_i\geq0$$. Note that the $$p_i\leq1$$ constraints are unnecessary since they are already covered by the non-overlap constraint. A feasible solution may then be recovered by setting the starting positions of the primers to be the location of the $$1$$ in primer $$i$$'s corresponding $$u_i$$ vector.

This relaxation reduces the number of binary variables to 240,000 and leads to a 10-100x speed up. Unfortunately this is still far too large a problem for GLPK to handle. Can anyone suggest an approach/re-formulation/other software that would make this problem tractable?

Thanks!

• Hi, if I understand correctly you are using GLPK to solve the mixed integer program resulting from relaxing $p$'s integrality? That's a first issue I see there - other solvers like CBC (github.com/coin-or/Cbc) should be used. Besides this, CLP (github.com/coin-or/Clp) tends to outperform GLPK for linear programming. I have some advances regarding this question but code isn't working well yet and needs more debugging. May I ask, is this from a studies, work or research? (looked for publications on genetic marker placement and wasn't able to find one that mentions the same optimization model) Jun 29, 2020 at 4:54
• @dhasson, correct, it's a MILP with $p$ relaxed. Could you elaborate on what makes CBC more performant than GLPK in this context? This is a research question that I confronted as part of my work, as far as I could tell, there wasn't any existing research that took this approach to primer design.
– Set
Jun 29, 2020 at 17:47
• A couple benchmarks are available, for example this one by Hans Mittelmann. In general CLP and CBC perform better than GLPK for linear programs and mixed integer programs, respectively. Can you use the NEOS Server? (see terms of use) In that case you'd be able to run the model on high performance optimization solvers (only for academic non-commercial research). Besides this, how big can $lb_u, lb_p, ub_u, ub_p$ be in real world scenarios (ie.for $G=30000$)? Still smaller than 100? Jun 29, 2020 at 22:38
• @dhasson Strictly speaking it's work related, although at this point it's unlikely that I would actually purchase a primer design based on the results; computing a solution would mostly be about scratching an intellectual itch. As per your question about the bounds, yes, they're all less than 100 in most real world scenarios that I can imagine.
– Set
Jun 29, 2020 at 22:46

# Formulating the problem

This MILP method you've constructed is pretty cool! But it's not the way I would choose to solve this problem. Rather, I would use dynamic programming.

To do so, recognize that we're essentially sliding a series of non-overlapping sum-windows across the dataset, subject to some constraints on the location of those windows.

Next, recognize that we can specify a state for the problem. For any given primer $$p, starting location for that primer $$s$$, and length of that primer $$l$$, there is an optimal solution.

Next, recognize that if I calculate the optimal solution for $$$$, then finding the optimal solution for $$$$ involves redoing many of the calculations for $$$$. That is, the problem has overlapping subproblems and optimal substructure. Recognizing states we've seen before and avoiding recalculating is therefore key to good performance.

Next, note that calculating how many mutated bases the primers intersect with naively takes $$O(p)$$ time per array position per primer length, for $$O(Gp^2)$$ time total. However, for any given primer length we can use a sliding window to calculate this in $$O(Gp)$$ time. Additionally, we can cache that calculation so we only perform it once per primer length. (We could do fancy things to avoid the $$p$$ factor, but that's unnecessary.)

Next, note that the problem as you've constructed it provides no benefit for choosing longer primer lengths. The solution is either to run the problem multiple times for each primer length being considered, or to prefer solutions that use longer primer lengths. I choose the latter option here. In your original formulation you could model this by modifying the objective function, like so: $$\min_p v^T\sum_i p_i-\left(\frac{1}{8 ub_p+1}\right) {\bf 1}^T \sum_i p_i$$ That is, we take the total number of "hot bits" in the primers, divide by the upper bound of hot bits, and subtract one. The effect is that for solutions with equal mutation coverage, the one with longer primers is preferred; however, the additional gain of long primers will never be sufficient to out-weight less mutation coverage.

The nice thing about this formulation is that we can efficiently solve the problem to optimality: on the randomly generated datasets I use below, I observed optimal objective values of approximately 8 with primer lengths of, e.g., 24 18 24 19 20 24 23 23.

Below, I describe both a Python and a C++ solution.

# A Python Solution

The Python solution takes 5.9 minutes (354s) and 1.9GB of RAM using the pypy3 interpreter (which is usually much faster than the standard python3 interpreter).

#!/usr/bin/env python3

from collections import deque
from functools import lru_cache
import copy
import random

def sliding_window_sum(a, size):
assert size>0
out     = []
the_sum = 0
q       = deque()
for i in a:
if len(q)==size:
the_sum -= q
q.popleft()
q.append(i)
the_sum += i
if len(q)==size:
out.append(the_sum)
return out

class Scoreifier:
def __init__(
self,
v,               #Array of mutations
lb_u:int   = 18, #Lower bound on inter-primer spacing
ub_u:int   = 60, #Upper bound on inter-primer spacing
lb_p:int   = 18, #Lower bound on primer length
ub_p:int   = 24, #Upper bound on primer length
pcount:int = 8   #Number of primers
):
#Problem attributes
self.v      = v
self.lb_u   = lb_u
self.ub_u   = ub_u
self.lb_p   = lb_p
self.ub_p   = ub_p
self.pcount = pcount
#Cache some handy information for later (pulls a factor len(p) out of the
#time complexity). Code is simplified at low cost of additional space by
#calculating subarray sums we won't use.
self.sub_sums = [[]] + [sliding_window_sum(v, i) for i in range(1, ub_p+1)]

@staticmethod
def _get_best(current_best, ret):
if current_best is None:
current_best = copy.deepcopy(ret)
elif ret["score"]<current_best["score"]:
current_best = copy.deepcopy(ret)
elif ret["score"]==current_best["score"] and ret["cum_len"]>current_best["cum_len"]:
current_best = copy.deepcopy(ret)

return current_best

@lru_cache(maxsize=None)
def _find_best_helper(
self,
p,        #Primer we're currently considering
start,    #Starting position for this primer
plen      #Length of this primer
):
#Don't consider primer location-length combinations that put us outside the
#dataset
if start>=len(self.sub_sums[plen]):
return {
"score":     float('inf'),
"cum_len":   -float('inf'),
"lengths":   [],
"positions": []
}
elif p==self.pcount-1:
return {
"score":     self.sub_sums[plen][start],
"cum_len":   plen,
"lengths":   [plen],
"positions": [start]
}

#Otherwise, find the best arrangement starting from the current location
current_best = None
for next_start in range(start+self.lb_u, start+self.ub_u+1):
for next_plen in range(self.lb_p, self.ub_p+1):
ret = self._find_best_helper(p=p+1, start=next_start, plen=next_plen)
current_best = self._get_best(current_best, ret)

current_best["score"]   += self.sub_sums[plen][start]
current_best["cum_len"] += plen
current_best["lengths"].append(plen)
current_best["positions"].append(start)

return current_best

def find_best(self):
#Consider all possible starting locations
current_best = None
for start in range(len(v)):
print(f"Start: {start}")
for plen in range(self.lb_p, self.ub_p+1):
ret = self._find_best_helper(p=0, start=start, plen=plen)
current_best = self._get_best(current_best, ret)

return current_best

G = 30_000
v = random.choices(population=[0,1], weights=[0.75, 0.25], k=G)

ret = Scoreifier(v=v).find_best()
print(ret)


# A C++ Solution

The C++ solution takes 56s on my machine using 295MB of RAM. With some care, it could be parallelized for faster performance. Better memory management would also give better performance.

#include <boost/container_hash/extensions.hpp>

#include <cassert>
#include <cstdlib>
#include <deque>
#include <iostream>
#include <vector>
#include <utility>
#include <unordered_map>

typedef std::vector<int> ivec;

struct Score {
double score   = std::numeric_limits<double>::infinity();
double cum_len = -std::numeric_limits<double>::infinity();
ivec lengths;
ivec positions;
bool operator<(const Score &o) const {
if(score<o.score)
return true;
else if(score==o.score && cum_len>o.cum_len)
return true;
else
return false;
}
};

typedef std::tuple<int,int,int> find_best_arg_type;

struct FBAThash {
std::size_t operator()(const find_best_arg_type &key) const {
return boost::hash_value(key);
}
};

using FBATmap = std::unordered_map<find_best_arg_type, Score, FBAThash>;

template<class T>
std::vector<T> sliding_window_sum(const std::vector<T> &v, const int size){
assert(size>0);
std::vector<T> out;
T the_sum = 0;
std::deque<T> q;
for(const auto &x: v){
if(q.size()==size){
the_sum -= q.front();
q.pop_front();
}
q.push_back(x);
the_sum += x;
if(q.size()==size)
out.push_back(the_sum);
}
return out;
}

class Scoreifier {
public:
ivec v;
const int lb_u;
const int ub_u;
const int lb_p;
const int ub_p;
const int pcount;

Scoreifier(const ivec &v, int lb_u, int ub_u, int lb_p, int ub_p, int pcount):
v(v), lb_u(lb_u), ub_u(ub_u), lb_p(lb_p), ub_p(ub_p), pcount(pcount)
{
//Cache some handy information for later (pulls a factor len(p) out of the
//time complexity). Code is simplified at low cost of additional space by
//calculating subarray sums we won't use.
sub_sums.emplace_back(); //Empty array for 0
for(int i=1;i<ub_p+1;i++)
sub_sums.push_back(sliding_window_sum(v, i));
}

Score find_best(){
//Consider all possible starting locations
Score current_best;
for(int start=0;start<v.size();start++){
std::cout<<"Start: "<<start<<"\n";
for(int plen=lb_p;plen<ub_p+1;plen++)
current_best = std::min(current_best,find_best_helper(0, start, plen));
}
return current_best;
}

private:
FBATmap visited;

std::vector<ivec> sub_sums;

Score find_best_helper(
const int p,     //Primer we're currently considering
const int start, //Starting position for this primer
const int plen   //Length of this primer
){
//Don't repeat if we've already solved this problem
const auto key = find_best_arg_type(p,start,plen);
if(visited.count(key)!=0)
return visited.at(key);

//Don't consider primer location-length combinations that put us outside the
//dataset
if(start>=sub_sums.at(plen).size())
return {};
else if(p==pcount-1)
return {(double)sub_sums.at(plen).at(start), (double)plen, {plen}, {start}};

//Otherwise, find the best arrangement starting from the current location
Score current_best;
for(int next_start=start+lb_u; next_start<start+ub_u+1; next_start++)
for(int next_plen=lb_p; next_plen<ub_p+1; next_plen++)
current_best = std::min(current_best, find_best_helper(p+1, next_start, next_plen));

current_best.score   += sub_sums[plen][start];
current_best.cum_len += plen;
current_best.lengths.push_back(plen);
current_best.positions.push_back(start);

visited[key] = current_best;

return current_best;
}
};

int main(){
const int G=30'000;

ivec v;
for(int i=0;i<G;i++){
v.push_back(rand()%100<25);
}

const auto sc = Scoreifier(v, 18, 60, 18, 24, 8).find_best();

std::cout<<"best_score      = "<<sc.score<<std::endl;
std::cout<<"best_cum_length = "<<sc.cum_len<<std::endl;

std::cout<<"best_lengths    = ";
for(const auto &x: sc.lengths)
std::cout<<x<<" ";
std::cout<<std::endl;

std::cout<<"best_positions  = ";
for(const auto &x: sc.positions)
std::cout<<x<<" ";
std::cout<<std::endl;

return 0;
}
$$$$
`
• I agree with you, was thinking another approach like a greedy construction heuristic or DP should suffice for this case, as primer lengths (and length range) are quite small comparing to the 30.000, added to the fact that ~75% bits of v are 0. By the way, nice solution. Jul 1, 2020 at 23:29

Let's consider the following model where $$p$$ integrality has been relaxed and call it P:

\begin{align} \underset{p}{\text{min}}&\;\;\;v^T\sum_i p_i\\ \text{subject to}&\;\;\;\sum_ip_i\leq{\bf 1} &(1)\\ &\;\;\;Ap_i - u_i \leq 0 &(2)\\ &\;\;\;{\bf 1}^Tu_i=1 &(3)\\ &\;\;\;lb_u\leq c^T(u_{i+1}-u_i)\leq ub_u &(4)\\ &\;\;\;lb_p\leq{\bf 1}^Tp_i\leq ub_p &(5)\\ &\;\;\;p\geq 0 &(6)\\ &\;\;\;u\in \{ 0,1 \} &(7) \end{align}

As stated in the comments, a first step is changing GLPK for another solver like CBC for performance reasons. Besides this I see a couple ways to accelerate solution time for this quite large MIP:

1. Add branching priorities to the integer variables such that the structure of the solution, the primer location patterns, are defined early in the branch-and-bound (B&B) procedure. For example, setting higher priority on branching the $$u_i$$ for the subset of indices $$I' = \{ i \in 1 \ldots G : |i-j| \leq ub_p \text{ for some j where v_j = 1} \}$$, does B&B performance improve?
2. Apply Lagrangian relaxation to the model: Note that constraints (2) are the only ones linking $$u$$ and $$p$$. Taking them out yields a simpler model that can be solved pretty fast (I made some tests locally and in many cases this has optimal value $$0$$) and can be decomposed into a linear program for $$p$$ and an integer program for $$u$$ solved separately.
3. Reformulation: It may be possible to state a different model where the variables represent primer combinations. This path sounds interesting as the primer lengths are rather small. Even though the resulting model will be large, a column generation or branch-and-price approach can be used.
4. Cutting plane generation: Relax model P, omitting the difficult constraint class (2). Iteratively solve the relaxation, find a subset of unfulfilled constraints from (2), and add that subset of constraints to the model, solving a sequence of several related sub-problems.
5. Find valid inequalities: If you can derive inequalities that can strengthen the formulation, turning its linear relaxation more tight but without discarding feasible solutions, it could also help with solution time. This is, inequalities that all feasible solutions $$(u, p)$$ must satisfy.

I will elaborate on some of them:

## Lagrangian relaxation

In this sense, you can introduce penalizations $$\mu \geq 0$$ to the linking constraint, adding a term $$\sum_{i} \mu_i^T (Ap_i - u_i)$$ to the objective function. The resulting problem can be decomposed into two smaller subproblems that depend on $$\mu$$ values:

A linear program for $$p$$ \begin{align} (LR^1_\mu) \underset{p} {\text{min}}&\;\;\;v^T\sum_i p_i + \sum_{i} \mu_i^T Ap_i\\ \text{subject to}&\;\;\;\sum_ip_i\leq{\bf 1} &(1)\\ &\;\;\;lb_p\leq{\bf 1}^Tp_i\leq ub_p &(5)\\ &\;\;\;p\geq 0 &(6) \end{align} and an integer program for $$u$$ \begin{align} (LR^2_\mu) \underset{u} {\text{max}}&\;\;\; \sum_{i} \mu_i^T u_i\\ \text{subject to}&\;\;\;{\bf 1}^Tu_i=1 &(3)\\ &\;\;\;lb_u\leq c^T(u_{i+1}-u_i)\leq ub_u &(4)\\ &\;\;\;u\in \{ 0,1 \} &(7) \end{align}

this kind of approach is known as Lagrangian decomposition. To solve it, the subgradient method can be used, and you'll end up with both a lower bound and an upper bound for the optimal value of P. I made tests on my computer and this works substantially faster than solving P. If you can take advantage of a subproblems' structure to devise a method to solve it without calling an LP or MIP solver, you could achieve an impressive performance. You'll also need to design a heuristic to approximate the Lagrangian relaxation's solution back to a solution of the original problem, that is make $$Ap \leq u$$ feasible.

## Column Generation formulation

In this case, solutions are modelled in a completely different way: Let's consider $$\mathcal{P}$$ the set of all possible primers, i.e. all sequences of $$1$$'s of length between $$lb_p$$ and $$ub_p$$. You want to select 8 (not necessarily distinct) elements of $$\mathcal{C}$$.

These decision variables can work: $$x_{pi} = 1$$ if pattern $$p \in \mathcal{P}$$ is put beginning at $$i$$, and $$0$$ otherwise. We know $$p$$'s length from beforehand, let's call it $$L_p$$. The constraints to take into account:

• Primer non-overlap - no pair of primers can use the same coordinate. For example: $$x_{pi} = 1 \wedge x_{qj} = 1 \Rightarrow$$ both their coordinates must not overlap with the other one's (use $$L_p$$ and $$L_q$$ to make this explicit)

• Primer contiguity - this is achieved by construction of the $$p \in \mathcal{P}$$.

• Respect minimum and maximum distance between primer positions - we must take into acount this, can be modeled as a logical constraint $$x_{pi} = 1 \Rightarrow x_{qj} = 0 \, \forall q \in \mathcal{C}, \, \forall j$$ outside the range $$[i + lb, i + ub]$$.

• The pattern positions must respect total length $$G$$, that is, $$x_{pi} = 1$$ => $$i \leq G - L_p$$.

• Must locate 8 primers: $$\sum_{p,i} = 8$$.

Depending on the modelling interface, can use these logical constraints as is or need the standard ways to turn them into inequality constraints. Some point might be missing on this last formulation, but hopefully you get the idea. For example, wiill need to build the sum of primers from these new variables for the new objective function.

This results in a large model, with a lot of variables and columns. The good thing is that from your comment, primers tend to be small (with respect to $$G$$). Moreover, column generation approaches work with a subset of columns (and variables). It's just a different paradigm for solving large MIPs, you can read more on it and on branch-and-price. You'll need to build a starting feasible model, e.g. choose a small subset of different patterns and locations that guarantee the initial relaxed model is feasible. Tipically that step isn't difficult.

What I would do: try a combination of 1. with 2. and/or 4. If the resulting performance isn't enough, or you're feeling creative, try to find valid inequalities to add. Alternatively, try column generation.

After that, a feasible solution for the original integer program can then be recovered as stated by @Thoth in the question, by setting the starting positions of the primers to be the location of the $$1$$ in primer $$i$$'s corresponding $$u_i$$ vector.