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<p_\textrm{max}$, 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 $<p=0, s=0, l=18>$, then finding the optimal solution for $<p=0,s=1,l=18>$ involves redoing many of the calculations for $<p=1,s=*,l=*>$. 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[0]
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;
}
```
