diff --git a/onepass.py b/onepass.py index 71f9e5b..047de1a 100644 --- a/onepass.py +++ b/onepass.py @@ -1,163 +1,157 @@ import numpy as np import timeit import gurobipy as gp from gurobipy import GRB def addvariables(Z): upperbounds=[] I=[] for M in database_indices: m=len(data['database_ncharges'][M]) I=I+[(i,j,M) for i in range(m) for j in range(n)] upperbounds.append(int(n/m)) x=Z.addVars(I, vtype=GRB.BINARY) # dummy variables y associated to number of times molecule M is picked y=Z.addVars(database_indices, vtype="I", lb=0, ub=upperbounds) print("Variables added.") return x,y,I def addconstraints(Z,x,y): - # injection into [n] - Z.addConstrs(x.sum('*',j,'*') <= 1 for j in range(n)) - # use at least 80% of target - Z.addConstr(x.sum() >= 0.8*n) + # bijection into [n] + Z.addConstrs(x.sum('*',j,'*') == 1 for j in range(n)) # additional constraints: take whole molecule or leave it out, and matching charges for M in database_indices: CM=data['database_ncharges'][M] m=len(CM) - # the number of indices of M used (counting multiple use) is at least 55% of its size - Z.addConstr(x.sum('*','*',M) >= 0.55*m*y[M]) + # the number of indices of M used (counting multiple use) is at least mol_percent of its size + Z.addConstr(x.sum('*','*',M) >= mol_percent*m*y[M]) # each index in M is used at most y[M] times Z.addConstrs(x.sum(i,'*',M) <= y[M] for i in range(m)) # ignore incompatible charges and atoms of charge 1 for i in range(m): for j in range(n): if(CM[i] != CT[j]): - #if(CM[i] != CT[j] or CM[i]==1): Z.addConstr(x[i,j,M]==0) - print("Constraints added.") return 0 def setobjective(Z,x,y): expr=gp.QuadExpr() print("Constructing objective function... ") key=0 for M in database_indices: key=key+1 Mol=data['database_CMs'][M] m=len(Mol) expr.addTerms(-m, y[M]) for i in range(m): for j in range(m): for k in range(n): for l in range(k,n): expr.add(x[i,k,M] * x[j,l,M], np.abs(T[k,l]-Mol[i,j])**2) print(key, " / ", size_database) Z.setObjective(expr, GRB.MINIMIZE) print("Objective function set.") return 0 # prints mappings of positions (indices+1) of each molecule (before preprocess) to positions inside target (before preprocess, but the hydrogens are at the end anyway) def print_sols(Z, x, y): SolCount=Z.SolCount for solnb in range(SolCount): print() print("--------------------------------") Z.setParam("SolutionNumber",solnb) print("Solution number", solnb+1, ", objective value", Z.PoolObjVal) for M in database_indices: amount_picked=int(np.rint(y[M].Xn)) if amount_picked != 0: m=len(data['database_ncharges'][M]) U=np.zeros((m,amount_picked)) # constructing U for i in range(m): k=0 for j in range(n): - if x[i,j,M].Xn==1 and sum(U[:,k]!=0) < 0.55*m: + if x[i,j,M].Xn==1 and sum(U[:,k]!=0) < mol_percent*m: U[i,k]=j+1 k=k+1 # reading U for k in range(amount_picked): if np.any(U[:,k] != 0): if k==0: print("Molecule", data['database_labels'][M], "has been picked", amount_picked, "time(s) ( size", len(data['database_ncharges'][M]), ", used", sum([x[i,j,M].Xn for i in range(m) for j in range(n)]), ")") print(k+1, end=": ") for i in range(m): if U[i,k]!=0: print(oldindex(i,M)+1, "->", U[i,k], end=", ") print("used", sum(U[:,k]!=0)) # converts new index (with hydrogens removed) to old index. -def oldindex(i, M, target): +def oldindex(i, M): k=0 notones=0 CM=olddata['database_ncharges'][M] while notones