Created
April 14, 2014 17:40
-
-
Save jasdumas/10668438 to your computer and use it in GitHub Desktop.
MITx - 6.02x
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#problem 2-1 | |
class SimpleVirus(object): | |
""" | |
Representation of a simple virus (does not model drug effects/resistance). | |
""" | |
def __init__(self, maxBirthProb, clearProb): | |
""" | |
Initialize a SimpleVirus instance, saves all parameters as attributes | |
of the instance. | |
maxBirthProb: Maximum reproduction probability (a float between 0-1) | |
clearProb: Maximum clearance probability (a float between 0-1). | |
""" | |
self.maxBirthProb = maxBirthProb | |
self.clearProb = clearProb | |
def getMaxBirthProb(self): | |
""" | |
Returns the max birth probability. | |
""" | |
return self.maxBirthProb | |
def getClearProb(self): | |
""" | |
Returns the clear probability. | |
""" | |
return self.clearProb | |
def doesClear(self): | |
""" | |
Stochastically determines whether this virus particle is cleared from the | |
patient's body at a time step. | |
returns: True with probability self.getClearProb and otherwise returns | |
False. | |
""" | |
if random.random() < self.getClearProb(): | |
return True | |
else: | |
return False | |
def reproduce(self, popDensity): | |
""" | |
Stochastically determines whether this virus particle reproduces at a | |
time step. Called by the update() method in the Patient and | |
TreatedPatient classes. The virus particle reproduces with probability | |
self.getMaxBirthProb * (1 - popDensity). | |
If this virus particle reproduces, then reproduce() creates and returns | |
the instance of the offspring SimpleVirus (which has the same | |
maxBirthProb and clearProb values as its parent). | |
popDensity: the population density (a float), defined as the current | |
virus population divided by the maximum population. | |
returns: a new instance of the SimpleVirus class representing the | |
offspring of this virus particle. The child should have the same | |
maxBirthProb and clearProb values as this virus. Raises a | |
NoChildException if this virus particle does not reproduce. | |
""" | |
if random.random() < self.getMaxBirthProb() * (1 - popDensity): | |
return SimpleVirus(self.maxBirthProb, self.clearProb) | |
else: | |
raise NoChildException ("This child does not reproduce") | |
class Patient(object): | |
""" | |
Representation of a simplified patient. The patient does not take any drugs | |
and his/her virus populations have no drug resistance. | |
""" | |
def __init__(self, viruses, maxPop): | |
""" | |
Initialization function, saves the viruses and maxPop parameters as | |
attributes. | |
viruses: the list representing the virus population (a list of | |
SimpleVirus instances) | |
maxPop: the maximum virus population for this patient (an integer) | |
""" | |
self.viruses = viruses | |
self.maxPop = maxPop | |
def getViruses(self): | |
""" | |
Returns the viruses in this Patient. | |
""" | |
return self.viruses | |
def getMaxPop(self): | |
""" | |
Returns the max population. | |
""" | |
return self.maxPop | |
def getTotalPop(self): | |
""" | |
Gets the size of the current total virus population. | |
returns: The total virus population (an integer) | |
""" | |
return len(self.viruses) | |
def update(self): | |
""" | |
Update the state of the virus population in this patient for a single | |
time step. update() should execute the following steps in this order: | |
- Determine whether each virus particle survives and updates the list | |
of virus particles accordingly. | |
- The current population density is calculated. This population density | |
value is used until the next call to update() | |
- Based on this value of population density, determine whether each | |
virus particle should reproduce and add offspring virus particles to | |
the list of viruses in this patient. | |
returns: The total virus population at the end of the update (an | |
integer) | |
""" | |
children = [] | |
# remove cleared viruses from list | |
for virus in self.getViruses(): | |
if virus.doesClear(): | |
self.viruses.remove(virus) | |
# update population density | |
popDens = float(self.getTotalPop()) / float(self.maxPop) | |
# determine which viruses reproduce | |
for virus in self.getViruses(): | |
try: | |
children.append(virus.reproduce(popDens)) | |
except NoChildException: pass | |
for c in children: | |
if self.getTotalPop() < self.maxPop: | |
self.viruses.append(c) | |
return self.getTotalPop() | |
#problem 3-1 | |
# Enter your definition for simulationWithoutDrug in this box | |
def simulationWithoutDrug(numViruses, maxPop, maxBirthProb, clearProb, | |
numTrials): | |
""" | |
Run the simulation and plot the graph for problem 3 (no drugs are used, | |
viruses do not have any drug resistance). | |
For each of numTrials trial, instantiates a patient, runs a simulation | |
for 300 timesteps, and plots the average virus population size as a | |
function of time. | |
numViruses: number of SimpleVirus to create for patient (an integer) | |
maxPop: maximum virus population for patient (an integer) | |
maxBirthProb: Maximum reproduction probability (a float between 0-1) | |
clearProb: Maximum clearance probability (a float between 0-1) | |
numTrials: number of simulation runs to execute (an integer) | |
""" | |
viruses = [] | |
result = [] | |
# instantiate viruses | |
for v in range(numViruses): | |
viruses.append(SimpleVirus(maxBirthProb, clearProb)) | |
# run simulation on patient | |
for t in range(numTrials): | |
patient = Patient(viruses, maxPop) | |
for s in range(300): | |
result.append(patient.update()) | |
# average results | |
finalResult = [] | |
for r in result: | |
finalResult.append(float(r) / float(numTrials)) | |
pylab.plot(finalResult) | |
pylab.title('SimpleVirus simulation') | |
pylab.xlabel('Time Steps') | |
pylab.ylabel('Average Virus Population') | |
pylab.legend('Average Virus Population') | |
pylab.show() | |
#problem 4-1 | |
# Enter your definition for the ResistantVirus class in this box. | |
# You'll enter your code for TreatedPatient on the next page. | |
class ResistantVirus(SimpleVirus): | |
""" | |
Representation of a virus which can have drug resistance. | |
""" | |
def __init__(self, maxBirthProb, clearProb, resistances, mutProb): | |
""" | |
Initialize a ResistantVirus instance, saves all parameters as attributes | |
of the instance. | |
maxBirthProb: Maximum reproduction probability (a float between 0-1) | |
clearProb: Maximum clearance probability (a float between 0-1). | |
resistances: A dictionary of drug names (strings) mapping to the state | |
of this virus particle's resistance (either True or False) to each drug. | |
e.g. {'guttagonol':False, 'srinol':False}, means that this virus | |
particle is resistant to neither guttagonol nor srinol. | |
mutProb: Mutation probability for this virus particle (a float). This is | |
the probability of the offspring acquiring or losing resistance to a drug. | |
""" | |
SimpleVirus.__init__(self, maxBirthProb, clearProb) | |
self.resistances = resistances | |
self.mutProb = mutProb | |
def getResistances(self): | |
""" | |
Returns the resistances for this virus. | |
""" | |
return self.resistances | |
def getMutProb(self): | |
""" | |
Returns the mutation probability for this virus. | |
""" | |
return self.mutProb | |
def isResistantTo(self, drug): | |
""" | |
Get the state of this virus particle's resistance to a drug. This method | |
is called by getResistPop() in TreatedPatient to determine how many virus | |
particles have resistance to a drug. | |
drug: The drug (a string) | |
returns: True if this virus instance is resistant to the drug, False | |
otherwise. | |
""" | |
if drug in self.getResistances(): | |
if self.getResistances()[drug] == True: return True | |
else: return False | |
else: | |
raise KeyError ('Not in resistances (neither true nor false)') | |
def reproduce(self, popDensity, activeDrugs): | |
""" | |
Stochastically determines whether this virus particle reproduces at a | |
time step. Called by the update() method in the TreatedPatient class. | |
A virus particle will only reproduce if it is resistant to ALL the drugs | |
in the activeDrugs list. For example, if there are 2 drugs in the | |
activeDrugs list, and the virus particle is resistant to 1 or no drugs, | |
then it will NOT reproduce. | |
Hence, if the virus is resistant to all drugs | |
in activeDrugs, then the virus reproduces with probability: | |
self.getMaxBirthProb * (1 - popDensity). | |
If this virus particle reproduces, then reproduce() creates and returns | |
the instance of the offspring ResistantVirus (which has the same | |
maxBirthProb and clearProb values as its parent). The offspring virus | |
will have the same maxBirthProb, clearProb, and mutProb as the parent. | |
For each drug resistance trait of the virus (i.e. each key of | |
self.resistances), the offspring has probability 1-mutProb of | |
inheriting that resistance trait from the parent, and probability | |
mutProb of switching that resistance trait in the offspring. | |
For example, if a virus particle is resistant to guttagonol but not | |
srinol, and self.mutProb is 0.1, then there is a 10% chance that | |
that the offspring will lose resistance to guttagonol and a 90% | |
chance that the offspring will be resistant to guttagonol. | |
There is also a 10% chance that the offspring will gain resistance to | |
srinol and a 90% chance that the offspring will not be resistant to | |
srinol. | |
popDensity: the population density (a float), defined as the current | |
virus population divided by the maximum population | |
activeDrugs: a list of the drug names acting on this virus particle | |
(a list of strings). | |
returns: a new instance of the ResistantVirus class representing the | |
offspring of this virus particle. The child should have the same | |
maxBirthProb and clearProb values as this virus. Raises a | |
NoChildException if this virus particle does not reproduce. | |
""" | |
# check if resistant to ALL active drugs | |
isResistant = True | |
for drug in activeDrugs: | |
if not self.isResistantTo(drug): | |
isResistant = False | |
# if resistant, may reproduce | |
if isResistant: | |
if (self.getMaxBirthProb() * (1 - popDensity)) > random.random(): # if passes random chance to reproduce | |
childResistances = self.getResistances().copy() | |
for drug in childResistances: | |
if self.mutProb > random.random(): # mutProb that resistance will switch | |
childResistances[drug] = not (self.isResistantTo(drug)) | |
return ResistantVirus(self.maxBirthProb, self.clearProb, childResistances, self.mutProb) | |
else: # not reproduce due to chance | |
raise NoChildException ('This child does not reproduce') | |
else: # not resistant to all drugs | |
raise NoChildException ('This child does not reproduce') | |
#problem 4-2 | |
# Enter your definitions for the ResistantVirus and TreatedPatient classes in this box. | |
class ResistantVirus(SimpleVirus): | |
""" | |
Representation of a virus which can have drug resistance. | |
""" | |
def __init__(self, maxBirthProb, clearProb, resistances, mutProb): | |
""" | |
Initialize a ResistantVirus instance, saves all parameters as attributes | |
of the instance. | |
maxBirthProb: Maximum reproduction probability (a float between 0-1) | |
clearProb: Maximum clearance probability (a float between 0-1). | |
resistances: A dictionary of drug names (strings) mapping to the state | |
of this virus particle's resistance (either True or False) to each drug. | |
e.g. {'guttagonol':False, 'srinol':False}, means that this virus | |
particle is resistant to neither guttagonol nor srinol. | |
mutProb: Mutation probability for this virus particle (a float). This is | |
the probability of the offspring acquiring or losing resistance to a drug. | |
""" | |
SimpleVirus.__init__(self, maxBirthProb, clearProb) | |
self.resistances = resistances | |
self.mutProb = mutProb | |
def getResistances(self): | |
""" | |
Returns the resistances for this virus. | |
""" | |
return self.resistances | |
def getMutProb(self): | |
""" | |
Returns the mutation probability for this virus. | |
""" | |
return self.mutProb | |
def isResistantTo(self, drug): | |
""" | |
Get the state of this virus particle's resistance to a drug. This method | |
is called by getResistPop() in TreatedPatient to determine how many virus | |
particles have resistance to a drug. | |
drug: The drug (a string) | |
returns: True if this virus instance is resistant to the drug, False | |
otherwise. | |
""" | |
if drug in self.getResistances(): | |
if self.getResistances()[drug] == True: return True | |
else: return False | |
else: | |
raise KeyError ('Not in resistances (neither true nor false)') | |
def reproduce(self, popDensity, activeDrugs): | |
""" | |
Stochastically determines whether this virus particle reproduces at a | |
time step. Called by the update() method in the TreatedPatient class. | |
A virus particle will only reproduce if it is resistant to ALL the drugs | |
in the activeDrugs list. For example, if there are 2 drugs in the | |
activeDrugs list, and the virus particle is resistant to 1 or no drugs, | |
then it will NOT reproduce. | |
Hence, if the virus is resistant to all drugs | |
in activeDrugs, then the virus reproduces with probability: | |
self.getMaxBirthProb * (1 - popDensity). | |
If this virus particle reproduces, then reproduce() creates and returns | |
the instance of the offspring ResistantVirus (which has the same | |
maxBirthProb and clearProb values as its parent). The offspring virus | |
will have the same maxBirthProb, clearProb, and mutProb as the parent. | |
For each drug resistance trait of the virus (i.e. each key of | |
self.resistances), the offspring has probability 1-mutProb of | |
inheriting that resistance trait from the parent, and probability | |
mutProb of switching that resistance trait in the offspring. | |
For example, if a virus particle is resistant to guttagonol but not | |
srinol, and self.mutProb is 0.1, then there is a 10% chance that | |
that the offspring will lose resistance to guttagonol and a 90% | |
chance that the offspring will be resistant to guttagonol. | |
There is also a 10% chance that the offspring will gain resistance to | |
srinol and a 90% chance that the offspring will not be resistant to | |
srinol. | |
popDensity: the population density (a float), defined as the current | |
virus population divided by the maximum population | |
activeDrugs: a list of the drug names acting on this virus particle | |
(a list of strings). | |
returns: a new instance of the ResistantVirus class representing the | |
offspring of this virus particle. The child should have the same | |
maxBirthProb and clearProb values as this virus. Raises a | |
NoChildException if this virus particle does not reproduce. | |
""" | |
# check if resistant to ALL active drugs | |
isResistant = True | |
for drug in activeDrugs: | |
if not self.isResistantTo(drug): | |
isResistant = False | |
# if resistant, may reproduce | |
if isResistant: | |
if (self.getMaxBirthProb() * (1 - popDensity)) > random.random(): # if passes random chance to reproduce | |
childResistances = self.getResistances().copy() | |
for drug in childResistances: | |
if self.mutProb > random.random(): # mutProb that resistance will switch | |
childResistances[drug] = not (self.isResistantTo(drug)) | |
return ResistantVirus(self.maxBirthProb, self.clearProb, childResistances, self.mutProb) | |
else: # not reproduce due to chance | |
raise NoChildException ('This child does not reproduce') | |
else: # not resistant to all drugs | |
raise NoChildException ('This child does not reproduce') | |
class TreatedPatient(Patient): | |
""" | |
Representation of a patient. The patient is able to take drugs and his/her | |
virus population can acquire resistance to the drugs he/she takes. | |
""" | |
def __init__(self, viruses, maxPop): | |
""" | |
Initialization function, saves the viruses and maxPop parameters as | |
attributes. Also initializes the list of drugs being administered | |
(which should initially include no drugs). | |
viruses: The list representing the virus population (a list of | |
virus instances) | |
maxPop: The maximum virus population for this patient (an integer) | |
""" | |
Patient.__init__(self, viruses, maxPop) | |
self.drugList = [] | |
def addPrescription(self, newDrug): | |
""" | |
Administer a drug to this patient. After a prescription is added, the | |
drug acts on the virus population for all subsequent time steps. If the | |
newDrug is already prescribed to this patient, the method has no effect. | |
newDrug: The name of the drug to administer to the patient (a string). | |
postcondition: The list of drugs being administered to a patient is updated | |
""" | |
if newDrug not in self.drugList: | |
self.drugList.append(newDrug) | |
def getPrescriptions(self): | |
""" | |
Returns the drugs that are being administered to this patient. | |
returns: The list of drug names (strings) being administered to this | |
patient. | |
""" | |
return self.drugList | |
def getResistPop(self, drugResist): | |
""" | |
Get the population of virus particles resistant to the drugs listed in | |
drugResist. | |
drugResist: Which drug resistances to include in the population (a list | |
of strings - e.g. ['guttagonol'] or ['guttagonol', 'srinol']) | |
returns: The population of viruses (an integer) with resistances to all | |
drugs in the drugResist list. | |
""" | |
resistantViruses = 0 | |
for drug in drugResist: | |
for virus in self.viruses: | |
if virus.isResistantTo(drug): | |
resistantViruses += 1 | |
return resistantViruses | |
def update(self): | |
""" | |
Update the state of the virus population in this patient for a single | |
time step. update() should execute these actions in order: | |
- Determine whether each virus particle survives and update the list of | |
virus particles accordingly | |
- The current population density is calculated. This population density | |
value is used until the next call to update(). | |
- Based on this value of population density, determine whether each | |
virus particle should reproduce and add offspring virus particles to | |
the list of viruses in this patient. | |
The list of drugs being administered should be accounted for in the | |
determination of whether each virus particle reproduces. | |
returns: The total virus population at the end of the update (an | |
integer) | |
""" | |
children = [] | |
# remove cleared viruses from list | |
for virus in self.getViruses(): | |
if virus.doesClear(): | |
self.viruses.remove(virus) | |
# update population density | |
popDens = float(self.getTotalPop()) / float(self.maxPop) | |
# determine which viruses reproduce | |
for virus in self.getViruses(): | |
try: | |
children.append(virus.reproduce(popDens, self.getPrescriptions())) | |
except NoChildException: pass | |
for c in children: | |
if self.getTotalPop() < self.maxPop: | |
self.viruses.append(c) | |
return self.getTotalPop() | |
#problem 5 | |
def simulationWithDrug(numViruses, maxPop, maxBirthProb, clearProb, resistances, mutProb, numTrials): | |
def runTrialDrug(elapsedTimeSteps, numViruses, maxPop, maxBirthProb, clearProb, resistances, mutProb, tempResistPopList): | |
viruses = []; | |
for i in range(numViruses): | |
virus = ResistantVirus(maxBirthProb, clearProb, resistances, mutProb); | |
viruses.append(virus); | |
patient = TreatedPatient(viruses, maxPop); | |
virusLevelsThisTrial = []; | |
ResistPopListThisTrial = []; | |
ResistPopListThisTrial.append(patient.getResistPop(patient.getPrescriptions())); | |
virusLevelsThisTrial.append(patient.getTotalPop()); | |
for i in range(elapsedTimeSteps): | |
if i == 150: | |
patient.addPrescription('guttagonol'); | |
virusLevelsThisTrial.append(patient.update()); | |
ResistPopListThisTrial.append(patient.getResistPop(['guttagonol'])); | |
for i in range(len(ResistPopListThisTrial)): | |
tempResistPopList.append(ResistPopListThisTrial[i]); | |
return virusLevelsThisTrial; | |
accumulatedVirusLevels = []; | |
virusLevelResults = []; | |
resistPopList = []; | |
for trial in range(numTrials): | |
tempResistPopList = []; | |
virusLevelResults = runTrialDrug(300, numViruses, maxPop, maxBirthProb, clearProb, resistances.copy(), mutProb, tempResistPopList); | |
if trial == 0: | |
accumulatedVirusLevels = virusLevelResults; | |
resistPopList = tempResistPopList[:]; | |
else: | |
for i in range(len(virusLevelResults)): | |
accumulatedVirusLevels[i] += virusLevelResults[i]; | |
for i in range(len(tempResistPopList)): | |
resistPopList[i] += tempResistPopList[i]; | |
accumulatedVirusLevels.remove(accumulatedVirusLevels[0]); | |
resistPopList.remove(resistPopList[0]); | |
for i in range(len(accumulatedVirusLevels)): | |
accumulatedVirusLevels[i] /= float(numTrials); | |
for i in range(len(resistPopList)): | |
resistPopList[i] /= float(numTrials); | |
pylab.plot(range(0, len(accumulatedVirusLevels)), accumulatedVirusLevels, label = "Total") | |
pylab.plot(range(0, len(resistPopList)), resistPopList, | |
label = "ResistantVirus") | |
pylab.title("ResistantVirus simulation") | |
pylab.xlabel("time step") | |
pylab.ylabel("# viruses") | |
pylab.legend(loc = "best") | |
pylab.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment