matrixchemistryNTOPmolecules.py

# -*- coding: utf-8 -*-
import numpy
#import hashlib; import base64
import binascii
global urllib

global DirectoryAC


def writeoptions(directory): #unused
    """a function that writes a file options.py with all the value of the options.
    The options are stored as python variables, and it is enough to execute the options.py program to get those options"""
    options=""
    print "writing options"
    global urllib
    global DirectoryAC
    options+="AC=\"matrixchemistry\"%s"%os.linesep
    options+="MainToPrimary=%s%s"%(MainToPrimary,os.linesep)
    options+="MainToSecondary=%s%s"%(MainToSecondary,os.linesep)
    options+="Result2Main=%s%s"%(Result2Main,os.linesep)
    options+="IllegalMolecules=%s%s"%(IllegalMolecules,os.linesep)
    #    m = hashlib.md5().update(options)
    #    urllib=base64.urlsafe_b64encode(m.digest())
    urllib="matrixchemistry_%s"%(binascii.crc32(options))
    libdir="%s%s%s%s%s"%(directory,os.sep,os.pardir,os.sep,urllib)
    #directoryAC=libdir

    if not os.path.isdir(libdir):
            os.mkdir(libdir)
    open("%s%sACType.py"%(libdir,os.sep),'w').write(options)

    options+="urllib='%s'%s"%(urllib,os.linesep)
    options+="mutation_rate_bit_flip=%s%s"%(mutation_rate_bit_flip,os.linesep)
    options+="starting_length=%s%s"%(starting_length,os.linesep)
    options+="starting_cell_size=%s%s"%(starting_cell_size,os.linesep)
    options+="LengthofExperiments=%s%s"%(LengthofExperiments,os.linesep)
    options+="TimeBetweenSavings=%s%s"%(TimeBetweenSavings,os.linesep)
    open("%s%smatrixchemistryoptions.py"%(directory,os.sep),'w').write(options)



class MoleculeRactionsDB(object):
    """a class the handles the Database of molecular reactions. Every time a reaction is tried it will call the function check.
    if the reaction is present in the DB, the result is returned. If not it is calculated, added to the DB and returned
    """
    def __repr__(self):  return 'MoleculeRactionsDB('+`self.DB`+')'
    def __str__(self):   return self.DB.__str__()
    #def __init__(self,DB={}): self.DB=DB
    def __init__(self,directory=".",DB={}):
        self.DB=shelve.open("%s%sMoleculeRactions.shelve"%(directory,os.sep),writeback=True)
        for k in DB.keys():
            self.DB[str(k)]=DB[k]
    def check(self,m0,m1):
        s01=`m0`+'vs'+`m1`
        try:
            return self.DB[s01]
        except KeyError:
            s10=`m1`+'vs'+`m0`
            rs01=basicreaction(m0,m1)
            WriteReaction(m0,m1,rs01) #sotres the reaction in the REA file
            self.DB[s01]=rs01  #slightly faster by taking away a hash call
            return rs01
##        except DBRunRecoveryError:
##        except :
##            s10=`m1`+'vs'+`m0`
##            rs01=basicreaction(m0,m1)
##            self.DB[s01]=rs01  #slightly faster by taking away a hash call
##            return rs01

def InitialiseReaFile(Molecules):
    """a function that initialisese the rea file"""
    global DirectoryAC, urllib
    Output=""
    Output+="# reactions rulesFrom: %s\n"%urllib
    Output+="# (generated by Speroni's AC. If you use this please cite my papers :-), especially the ones where I present it )\n"
    Output+="\n"
    Output+="\n"
    Output+="\n"
    Output+="# Number of Components: \n"
    Output+="%s\n"%len(Molecules)
    Output+="\n"
    Output+="# Components:\n"
    for m in Molecules:
        Output+="%s\n"%m
    Output+="\n"
    Output+="# Number of Reactions: \n"
    Output+="#to be calculated by hand\n"
    Output+="\n"
    Output+="\n"
    Output+="# Reactions:\n"
    open("%s%sACReactions.rea"%(DirectoryAC,os.sep),'a').write(Output)

def WriteReaction(m1,m0,rs01):
    """a function that writes a reaction to the REA file"""
    global DirectoryAC
    if rs01:
        ReactionLine= "1 %s  1 %s -> 1 %s &\n"%(m1,m0,rs01)
        open("%s%sACReactions.rea"%(DirectoryAC,os.sep),'a').write(ReactionLine)
	

class MoleculeMutationsDB(object):
    """
    a class the handles the Database of molecular mutations: given two molecules what is the Levensteinh distance between them.
    Once we need to show the distance between two molecules, we call the check function.
    If the distance is present in the DB, the result is returned. If not it is calculated, added to the DB and returned
    """
    def __repr__(self):  return 'MoleculeMutationsDB('+`self.DB`+')'
    def __str__(self):   return self.DB.__str__()
#    def __init__(self,DB={}): self.DB=DB
    def __init__(self,directory=".",DB={}):
        self.DB=shelve.open("%s%sMoleculeMutations.shelve"%(directory,os.sep),writeback=True)
        for k in DB.keys():
            self.DB[str(k)]=DB[k]

    def check(self,m0,m1):
        s01=`m0`+'vs'+`m1`
        try:
            return self.DB[s01]
        except KeyError:
            s10=`m1`+'vs'+`m0`
            rs01=LevDistance(moleculelibrary.check(m0).mainlist,moleculelibrary.check(m1).mainlist)
            self.DB[s01]=rs01
            return rs01


#mutationsDB=MoleculeMutationsDB()


def basicreaction(m0,m1):
    """A molecule is generated by the product of the primary matrix per the secondary matrix of the second."""
    matrix1=moleculelibrary.check(m0).primarymatrix
    matrix2=moleculelibrary.check(m1).secondarymatrix
    resultingmatrix=matrix1 * matrix2
    mainlist=resultmatrix2mainlist(resultingmatrix)
    value,size=mainlist2valuesize(mainlist)
    if tuple(mainlist) in IllegalMolecules:
#        print "I just tried to calculate %s and %s but the result is illegal"%(m0,m1)
        return None
    name=valuesize2name(value,size)
    return name

def react(molecule0, molecule1):
    """picks two molecules, checks if they can react
    finds out the resulting molecule and returns it"""
    return reactionDB.check(molecule0,molecule1) 

def ReturnAllMolecules():
    size=starting_length
    print NRowsFolded
    startingvalue=int('0'*size,2)
    endvalue=int('1'*size,2)
    mol=[]
    for v in range(startingvalue,endvalue+1):
        structuremolecule=tuple(valuesize2mainlist(v,size))
        if structuremolecule not in IllegalMolecules:
            namemolecule=valuesize2name(v,size)
            print "molecule %s, name = %s, structure = %s"%(v,namemolecule,structuremolecule)
            mol.append(namemolecule)
    return mol

def testrightidentity(molecule,context):
    for m in context:
        if reactionDB.check(m,molecule)!=m:
            return False
    return True

def testleftidentity(molecule,context):
    for m in context:
        if reactionDB.check(molecule,m)!=m:
            return False
    return True

def testrightzero(molecule,context):
    for m in context:
        if reactionDB.check(m,molecule)!=molecule:
            return False
    return True

def testleftzero(molecule,context):
    for m in context:
        if reactionDB.check(molecule,m)!=molecule:
            return False
    return True


    

    
###########################################
######## Molecule #########################
###########################################


class Molecule(object):
    """the class Molecule holds each molecule. Each molecule is present in various formats:
    As a list, as a string, as a couple value, size, or as a single Long number (name).
    Various functions are used to convert from one representation to the other.
    Then all the representations are stored in the dictionary "moleculelibrary".
    The molecule has a primary structure and a folded (aka secondary) structure.
    There are various possible foldings. AT the moment the default one is organised by taking one number every two in the primary string."""
    resultmatrix=primarymatrix=secondarymatrix=numpy.matrix([0])
    mainlist=[]      #the primary structure of the molecule, i.e. the list of bases//needed for (insertions and deletion mutations)       #the secondary structure of the molecule, i.e. the structure once it is folded
    primarystring=''    #the primary structure coded as a string for ease of certain operations
    name=value=0L
    size=0
        
    def __init__(self,ValueOrName=0L,SizeIfValue=0):
        """
        the initializer of the MatrixMolecule class.
        If you pass no valuse it will make a default molecule which with the present code is the molecule '100001010'.
        by just uncommenting two lines it can be the molecule '1'*starting length, or a random molecule again of length starting_length
        if you pass a single value it consider that value to be the 'name' of the molecule.
        If two values are passed they are considered to the the Value and the size.
        all the various representations of the string are calculated.
        """
        if not SizeIfValue: #only one value has been passed
            if not ValueOrName:#no value has been passed, I make it up
                self.size=starting_length
#                self.size=9
                #self.value=int('1'*starting_length,2)
                #self.value=int('100001010',2)
#                self.value=int('123456789',10)
                self.value=random.getrandbits(self.size)
                self.name=int(self.value+2**self.size)
            else:#only one value has been passed: name
                assert ValueOrName>3
                self.name=int(ValueOrName)
                self.size=int(math.floor(math.log(self.name,2)))
                self.value=self.name-2**self.size
        else:#two values has been passed: value, size
            assert ValueOrName<2**SizeIfValue
            self.value=ValueOrName
            self.size=SizeIfValue
            self.name=int(self.value+2**self.size)

        self.mainlist=valuesize2mainlist(self.value,self.size)
#        self.primarystring=''.join(self.primarylist)

        self.primarymatrix=mainlist2primarymatrix(self.mainlist)
        self.secondarymatrix=mainlist2secondarymatrix(self.mainlist)

        self.primarystring=self.primarymatrix.tolist()
        self.secondarystring=self.secondarymatrix.tolist()

    def __repr__(self):  return 'Molecule('+`self.name`+')'
    def __str__(self):   return "['"+self.primarymatrix+"','"+self.secondarymatrix+"']"

    def reproduce(self):
        """
        Makes a copy of the molecule, with mutations. The mutations depends on the constants in the options file that define
        the probability to have an added base, to lose a base, and to have a base change from one to another.
        """
        daughtervalue=self.value
        mutations=0L
        for t in range(0,self.size):
            if random.random()<mutation_rate_bit_flip:
                    mutations=mutations+2**t
        daughtervalue=daughtervalue^mutations
        return daughtervalue,self.size
    
def mainlist2primarymatrix(mainlist):
    " calculates and returns the primary list from the primary matrix"
    tofold=[mainlist[t] for t in MainToPrimary]
    #return numpy.matrix(tofold).reshape(3,3)
	#    return numpy.matrix(tofold).reshape(2,2)
    return numpy.matrix(tofold).reshape(NRowsFolded,NColumnsFolded)

def mainlist2secondarymatrix(mainlist):
    " calculates and returns the primary list from the primary matrix"
#    print mainlist
#    print MainToSecondary
    tofold=[mainlist[t] for t in MainToSecondary]
#    return numpy.matrix(tofold).reshape(2,2)
    return numpy.matrix(tofold).reshape(NRowsFolded,NColumnsFolded)
#    return numpy.matrix(tofold).reshape(3,3)

def valuesize2mainlist(value,size):
    "A function that taking the value and size of a molecule returns the primary structure"
    mainstructurelist=[0]*size
    offset=size-1
    while value:
        mainstructurelist[offset]=value&1
        value=value>>1
        offset-=1
    return mainstructurelist
def mainlist2valuesize(mainlist):
    " calculates and returns the value and the size from the primary structure"
    temp=[str(t) for t in mainlist]
    primarystring=''.join(temp)
    value=int(primarystring,2)
    size=len(mainlist)
    return value,size
def name2valuesize(name):
    "calculates and returns the value and the size from the name"
    size=int(math.floor(math.log(name,2)))
    return int(name-2**size),int(size)
def valuesize2name(value,size):
    "calculates and returns the name from the value and size"
    return int(value+2**size)


def resultmatrix2mainlist(resmatrix):
    "calculates and returns the value and size"
    #aslist=resmatrix.reshape(1,9).tolist()[0]
    aslist=resmatrix.reshape(1,starting_length).tolist()[0]	
    mainlist=[Result2Main[int(x)]   for x in aslist]
    return mainlist

def colormolecule1(molname):
    saturation=0.5
    if reactionDB.check(molname,molname)==molname:
        hue=0.0
    else:
        hue=0.66666
    lightness=0.9
    col=colorsys.hsv_to_rgb(hue, saturation, lightness)
    return col

def colormolecule2(molname):
    mol=moleculelibrary.check(molname)
    red=mol.mainlist[0]*128+mol.mainlist[1]*64+mol.mainlist[1]*32
    green=mol.mainlist[3]*128+mol.mainlist[4]*64+mol.mainlist[5]*32
    blue=mol.mainlist[6]*128+mol.mainlist[7]*64+mol.mainlist[8]*32
    
    return (red/256.0,green/256.0,blue/256.0)

def ColorLayout():
    while 1:
        yield colormolecule1
        yield colormolecule2


def writereact(t,h,option,molecules,numbermolecules):
    op=option%6
    if not op:   return ""
    if op==1 :
        res=reactionDB.check(t,h)
        if res:  return "%i\t"%res
        else  :  return "\t"
    if op==2 :
        res=reactionDB.check(t,h)
        if t==h:
            if   res==t: return "*\t"
            elif res   : return "%i\t"%res
            else       : return "\t"
        else:
            if   res==t: return ">\t"
            elif res==h: return "^\t"
            elif res   : return "%i\t"%res
            else       : return "\t"
    if op==3 :
        res=reactionDB.check(t,h)
        if t==h:
            if   res==t:        return "*\t"
            elif res   :        return "%i\t"%res
            else       :        return "\t"
        else:
            ser=reactionDB.check(h,t)
            if   res==t:
                if   ser==t:    return "*->\t"
                elif ser==h:    return "==>\t"
                elif ser:       return "->?\t"
                else :          return " ?\t"
            elif   res==h:
                if   ser==t:    return "==^\t"
                elif ser==h:    return "*^\t"
                elif ser:       return "^?\t"
                else :          return "^\t"
            elif res:
                if   ser==t:    return "%i^\t"%res
                elif ser==h:    return "%i>\t"%res
                elif ser:
                    if ser==res:return "*%i\t"%res
                    else:       return "%i?\t"%res
                else :          return "%i\t"%res
            else:
                if ser==t:      return " ^\t"%res
                elif ser==h:    return " >\t"%res
                elif ser:       return " ?\t"%res
                else:           return "\t"
    if op==4 :
        res=reactionDB.check(t,h)
        try:                        ind=molecules.index(res)+1
        except ValueError:          ind=0
        if t==h:
            if   res==t:                    return "*\t"
            elif res:
                if ind==0:                  return "NEW\t"
                elif ind<=numbermolecules:   return "# %i\t"%ind
                else:                       return "%i\t"%res
            else:                           return "\t"
        else:
            if   res==t:                    return ">\t"
            elif res==h:                    return "^\t"
            elif res:
                if ind==0:                  return "NEW\t"
                elif ind<=numbermolecules:   return "# %i\t"%ind
                else:                       return "%i\t"%res
            else:                           return "\t"
    if op==5 :
        res=reactionDB.check(t,h)
        try:                        ind=molecules.index(res)+1
        except ValueError:          ind=0
        if t==h:
            if res:
                if ind==0:                  return "NEW\t"
                elif ind<=numbermolecules:   return "# %i\t"%ind
                else:                       return "%i\t"%res
            else:                           return "\t"
        else:
            if res:
                if ind==0:                  return "NEW\t"
                elif ind<=numbermolecules:   return "# %i\t"%ind
                else:                       return "%i\t"%res
            else:                           return "\t"



       
## moved to aclib making it general
##def GeneratedMolecules(molecules):
##    """
##    Returns the list of all the molecules generated by a set of molecules, with their relative multiplicity
##    """
##    results=mydict()
##    for m in molecules:
##        for n in molecules:
##            res=reactionDB.check(m,n)
##            if res!= None:
##                try:
##                    results[res]+=1
##                except KeyError:
##                    results[res]=1
##    return results
##
##