Commit 31ace2bd authored by Matthias Schroeder's avatar Matthias Schroeder

minor changes and some speed optimizations

parent 9eae25d0
#!/usr/bin/python3
#
import numpy as np
import os
import sys
from time import gmtime, strftime
import numpy as np
MAX_MASS = 2000
LINES_MOL = 2000
......@@ -10,15 +12,15 @@ OUTDIR = 'masses/'
# description of elements
elements = (
{ "formel":'C' , "weight": 12 , "min":1, "max":167, "val": 4, "spec": 0 },
{ "formel":'C13' , "weight": 13.003354835 , "min":0, "max": 3, "val": 4, "spec": 0 },
{ "formel":'C' , "weight": 12 , "min":1, "max":170, "val": 4, "spec": 0 },
{ "formel":'C13' , "weight": 13.003354835 , "min":0, "max": 3, "val": 4, "spec": 2 },
{ "formel":'H' , "weight": 1.0078250 , "min":1, "max":800, "val": 1, "spec": 0 },
{ "formel":'O' , "weight": 15.99491461956, "min":0, "max": 90, "val": 2, "spec": 0 },
{ "formel":'O18' , "weight": 17.9991610 , "min":0, "max": 1, "val": 2, "spec": 0 },
{ "formel":'N' , "weight": 14.0030740048 , "min":0, "max": 50, "val": 5, "spec": 0 },
{ "formel":'N15' , "weight": 15.000108899 , "min":0, "max": 1, "val": 5, "spec": 0 },
{ "formel":'S' , "weight": 31.972072 , "min":0, "max": 30, "val": 6, "spec": 0 },
{ "formel":'S34' , "weight": 33.96786690 , "min":0, "max": 1, "val": 6, "spec": 0 },
{ "formel":'O' , "weight": 15.99491461956, "min":0, "max": 91, "val": 2, "spec": 0 },
{ "formel":'O18' , "weight": 17.9991610 , "min":0, "max": 1, "val": 2, "spec": 2 },
{ "formel":'N' , "weight": 14.0030740048 , "min":0, "max": 51, "val": 5, "spec": 0 },
{ "formel":'N15' , "weight": 15.000108899 , "min":0, "max": 1, "val": 5, "spec": 2 },
{ "formel":'S' , "weight": 31.972072 , "min":0, "max": 31, "val": 6, "spec": 0 },
{ "formel":'S34' , "weight": 33.96786690 , "min":0, "max": 1, "val": 6, "spec": 2 },
{ "formel":'P' , "weight": 30.973763 , "min":0, "max": 20, "val": 5, "spec": 0 },
{ "formel":'Cl' , "weight": 34.9688527 , "min":0, "max": 2, "val": 1, "spec": 1 },
{ "formel":'Cui' , "weight": 62.92959772 , "min":0, "max": 2, "val": 1, "spec": 1 },
......@@ -45,19 +47,24 @@ for i in range(len(elements)):
idx[elements[i]['formel']] = i
for i in range(len(elements)):
if elements[i]['max']>10:
elements[i]['max'] = 10
if elements[i]['spec'] == 1:
# if elements[i]['max']>10:
# elements[i]['max'] = 10
if elements[i]['spec'] > 0:
elements[i]['max'] = 0
elements[idx['P']]['max'] = 0
elements[idx['S34']]['max'] = 0
elements[idx['C13']]['max'] = 0
elements[idx['P']]['min'] = 0 # replace by
elements[idx['P']]['max'] = elements[idx['P']]['min']
# prepare some arrays
moltype = np.dtype([('strct', np.int16, len(elements)), ('weight', np.float64)])
moltype = np.dtype([('weight', np.float64), ('name', np.str_, 512)])
mols= np.ndarray(shape=(MAX_MASS,LINES_MOL), dtype=moltype)
posi= [0] * MAX_MASS
def info(s):
print("%s %s" % (strftime("%H:%M:%S", gmtime()), s))
sys.stdout.flush()
def check_strct(s,w):
""" Check for relations between element, total mass, bla bla,
return True for seemingly valid element,
......@@ -68,9 +75,9 @@ def check_strct(s,w):
# Valence check
val = 0
specs = 0
for i in range(LE):
val += s[i] * elements[i]['val']
if elements[i]['spec'] == 1:
for (i,e) in enumerate(elements):
val += s[i] * e['val']
if e['spec'] == 1:
specs += s[i]
if val % 2 == 1:
......@@ -80,24 +87,27 @@ def check_strct(s,w):
return False
totC= (s[idx['C']] + s[idx['C13']])
totS= (s[idx['S']] + s[idx['S34']])
totO= (s[idx['O']] + s[idx['O18']])
totN= (s[idx['N']] + s[idx['N15']])
db = 1 + 0.5 * (totC * 2 - s[idx['H']] - s[idx['Cl']] + s[idx['N']] + s[idx['P']])
db = 1 + 0.5 * (totC * 2 - s[idx['H']] - s[idx['Cl']] + totN + s[idx['P']])
if db<0:
return False
if (s[idx['O']] + s[idx['O18']]) / totC >2.1:
if totO / totC >2.1:
return False
if s[idx['H']] / totC > 8:
return False
if s[idx['S']] > totC:
if totS > totC:
return False
if s[idx['P']] > totC:
return False
if s[idx['N']] > 1.25 * totC:
if totN > 1.25 * totC:
return False
return True
......@@ -106,16 +116,16 @@ def check_strct(s,w):
def format_strct(s):
""" Formats molecule in a nice way """
m=''
for i in range(LE):
if s[0][i]>0:
m += elements[i]['formel'] + '_' + str(s[0][i]) + ' '
return "%16.10f\t%s" % (s[1], m)
for (i,e) in enumerate(elements):
if s[i]>0:
m += e['formel'] + '_' + str(s[i]) + ' '
return m
def calc_weight(s):
""" return total weight of molekule descripte by s """
w=0
for i in range(LE):
w= w + s[i] * elements[i]['weight']
for (i,e) in enumerate(elements):
w= w + s[i] * e['weight']
return w
def flush_mass(w):
......@@ -124,12 +134,10 @@ def flush_mass(w):
in case there is something """
if posi[w]>0:
filename= FILEPAT % w
print("Filename %s, cnt=%d" % (filename,posi[w]))
return
f = open(os.path.join(OUTDIR, filename), 'a')
for l in range(posi[w]):
f.write(format_strct(mols[w][l]) + "\n")
f.close()
# print("Filename %s, cnt=%d" % (filename,posi[w]))
s = [ "%16.10f\t%s\n" % (mols[w][l][0], mols[w][l][1]) for l in range(posi[w]) ]
with open(os.path.join(OUTDIR, filename), 'a') as f:
f.write("".join(s))
posi[w]=0
......@@ -138,26 +146,29 @@ if __name__ == "__main__":
strct= [ e['min'] for e in elements ]
# start the work
info("Started")
while True:
w= calc_weight(strct)
if check_strct(strct, w):
# we have a valid molecule
widx= int(w)
mols[widx][posi[widx]]= (strct[:], w)
mols[widx][posi[widx]]= (w, format_strct(strct))
posi[widx]+= 1
if posi[widx]>=LINES_MOL:
flush_mass(widx)
# build next molecule
inc_done = False
for i in range(LE):
if strct[i]>=elements[i]['max']:
for (i,e) in enumerate(elements):
if strct[i]>=e['max']:
# hit top, -> fall down and continue to increment next position
strct[i]=elements[i]['min']
strct[i]=e['min']
else:
# increment next position and done
strct[i]= strct[i]+1
inc_done =True
if i>6:
info("increment %s to %d" % (e['formel'], strct[i]))
break
# no increment - so we had had reached the last element
if not inc_done:
......@@ -166,3 +177,5 @@ if __name__ == "__main__":
# flush the rest ....
for l in range(MAX_MASS):
flush_mass(l)
info("Done")
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment