Page MenuHomec4science
Files F62839492

onepass.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 16, 00:30

onepass.py
View Options

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[targetname+"_amons_ncharges"][M])
CM=data[targetname+"_amons_ncharges"][M]
I=I+[(i,j,M,G) for G in range(maxduplicates) for i in range(m) for j in range(n) if CM[i] == CT[j]] # if condition excludes j; i always takes all m values
x=Z.addVars(I, vtype=GRB.BINARY)
print("Variables added.")
return x,I
def addconstraints(Z,x,I):
# bijection into [n]
Z.addConstrs(x.sum('*',j,'*', '*') == 1 for j in range(n))
# each i of each group is used at most once
for M in database_indices:
CM=data[targetname+"_amons_ncharges"][M]
m=len(CM)
Z.addConstrs(x.sum(i,'*',M,G) <= 1 for i in range(m) for G in range(maxduplicates))
print("Constraints added.")
return 0
# objective value should then be square rooted in the end (doesn't change optimality)
def setobjective(Z,x,I):
print("Constructing objective function... ")
key=0
expr=gp.QuadExpr()
for k in range(n):
for l in range(n):
expr += T[k,l]**2
for M in database_indices:
key=key+1
for G in range(maxduplicates):
Mol=data[targetname+"_amons_CMs"][M]
m=len(Mol)
for (i,k) in [v[:2] for v in I if v[2:]==(M,G)]:
for (j,l) in [v[:2] for v in I if v[2:]==(M,G)]:
expr += (Mol[i,j]**2 - 2*T[k,l]*Mol[i,j])*x[i,k,M,G]*x[j,l,M,G]
print(key, " / ", size_database)
Z.setObjective(expr, GRB.MINIMIZE)
print("Objective function set.")
return 0
# prints mappings of positions (indices+1) of each molecule to positions inside target
def print_sols(Z, x, I):
SolCount=Z.SolCount
print("Target has size", n)
for solnb in range(SolCount):
print()
print("--------------------------------")
Z.setParam("SolutionNumber",solnb)
print("Solution number", solnb+1, ", objective value (square rooted)", np.sqrt(Z.PoolObjVal))
for M in database_indices:
groups=[]
for G in range(maxduplicates):
if sum([x[v].Xn for v in I if v[2:]==(M,G)]) >= 1:
groups.append(G)
amount_picked=len(groups)
for k in range(amount_picked):
G=groups[k]
m=len(data[targetname+"_amons_ncharges"][M])
label=data[targetname+"_amons_labels"][M]
if k==0:
print("Molecule", label, "has been picked", amount_picked, "time(s) ( size", m, ", used", sum([x[v].Xn for v in I if v[2]==M]), ")")
print(k+1, end=": ")
for (i,j) in [v[:2] for v in I if v[2:]==(M,G) and x[v].Xn==1]:
print(i+1, "->", j+1, end=", ")
print()
def main():
# construction of the model
start=timeit.default_timer()
Z = gp.Model()
Z.setParam('OutputFlag',1)
x,I=addvariables(Z)
addconstraints(Z,x,I)
setobjective(Z,x,I)
stop=timeit.default_timer()
print("Model setup: ", stop-start, "s")
# optimization
# trying parameters from Z.tune() and others.
Z.setParam("DegenMoves", 2)
Z.setParam("ZeroHalfCuts", 2)
Z.setParam("BQPCuts", 2)
Z.setParam("RLTCuts", 2)
# time (seconds) and gap (%) limit, (max) number of solutions kept in memory
#Z.setParam("MIPGap", 0)
Z.setParam("TimeLimit", timelimit)
Z.setParam("PoolSolutions", numbersolutions)
print("------------")
print("Optimization")
print("------------")
Z.optimize()
print("Optimization runtime: ", Z.RunTime, "s")
if(Z.status == 3):
print("Model was proven to be infeasible.")
return 1
print_sols(Z,x,I)
return 0
# modifiable global settings
target_index=0 # 0, 1, or 2
maxduplicates=3 # number of possible copies of each molecule of the database
timelimit=180
numbersolutions=5
# global constants
data=np.load("amons_data.npz", allow_pickle=True)
targetdata=np.load("datafilter.npz", allow_pickle=True)
CT=targetdata['target_ncharges'][target_index]
T=targetdata['target_CMs'][target_index]
n=len(CT)
targetname=["qm9", "vitc", "vitd"][target_index]
size_database=len(data[targetname+"_amons_labels"])
database_indices=range(size_database)
main()

Event Timeline

c4science · Help