PS9

gistfile1.txt

# 6.00 Problem Set 9

import numpy
import random
import pylab
from ps8b_precompiled_27 import *

#
# PROBLEM 1
#        
def simulationTwoDrugsDelayedTreatment(numViruses, maxPop, maxBirthProb, clearProb,
                               resistances, mutProb, numTrials):
    """
    Runs simulations and make histograms. 
    *Note: I used this for problem 1 AND 2. Problem 2 suggested that the only variable passed into 
    simulationTwoDrugsDelayedTreatment would be numTrials. As you can see above, I didn't do that.                    

    Runs numTrials simulations to show the relationship between administration
    of multiple drugs and patient outcome.

    Histograms of final total virus populations are displayed for lag times of
    300, 150, 75, 0 timesteps between adding drugs (followed by an additional
    150 timesteps of simulation).                     

    numTrials: number of simulation runs to execute (an integer)
    """
    totalViruses = None
    resistantViruses = None
    totalList = []
    lessThan50Count = 0
    percentCured = 0.0

    """
    Run the simulation for 150 time steps before administering guttagonol to the patient. 
    Then run the simulation for 300, 150, 75, and 0 time steps before administering a second drug, grimpex, to the patient.
    Finally, run the simulation for an additional 150 time steps.
    Change the below values to simulate the above.
    """
    stepsBeforeDrug1 = 150
    stepsBeforeDrug2 = 0
    finalSimulation = 150
    totalNumSteps = stepsBeforeDrug1 + stepsBeforeDrug2 + finalSimulation
    print str(totalNumSteps)
    
    for i in xrange(0, numTrials): #Trials loop
        (total, resistant,) = runDrugSimulation(numViruses, maxPop, maxBirthProb, clearProb,
                               resistances, mutProb, stepsBeforeDrug1, stepsBeforeDrug2, totalNumSteps)


        #print total[-1]
        if total[-1] < 50:
            lessThan50Count += 1

        totalList.append(total[-1])
        
        if totalViruses == None:
            totalViruses = total
            resistantViruses = resistant
        else:
            for j in xrange(0, len(total)):
                totalViruses[j] += total[j]
                resistantViruses[j] += resistant[j]

    #% Cured
    percentCured = float(lessThan50Count)/float(numTrials)
    print percentCured

    numBins = 50
    pylab.hist(totalList, bins = numBins, label = str(percentCured))
    pylab.title('Treated Patients Cured - D1 App: ' + str(stepsBeforeDrug1)
                + ' D2 App: ' + str(stepsBeforeDrug2) + ' Total Steps: ' + str(totalNumSteps))
    xmin, xmax = pylab.xlim()
    ymin, ymax = pylab.ylim()
    pylab.xlabel('Viruses')
    pylab.ylabel('Trials')
    pylab.legend(loc='best')
    pylab.figure(1)
    pylab.savefig('Histogram: ' + str(stepsBeforeDrug1) + ' ' + str(stepsBeforeDrug2) + ' ' + str(totalNumSteps))
    pylab.show()

def runDrugSimulation(numViruses, maxPop, maxBirthProb, clearProb, resistances,
                      mutProb, stepsBeforeDrug1, stepsBeforeDrug2, totalNumSteps):
    """ Helper function for doing one actual simulation run with drug applied """
    assert stepsBeforeDrug1 <= totalNumSteps
    assert stepsBeforeDrug2 <= totalNumSteps
    viruses = []
    for i in xrange(0, numViruses):
        viruses.append(ResistantVirus(maxBirthProb, clearProb, resistances, mutProb))
 
    patient = TreatedPatient(viruses, maxPop)
    numVirusesEachStep = []
    numResistantVirusesEachStep = []
    for i in xrange(0, totalNumSteps):
        if i == stepsBeforeDrug1:
            patient.addPrescription('guttogonal')
        numVirusesEachStep.append(patient.update())
        numResistantVirusesEachStep.append(patient.getResistPop(['guttogonal']))
        if i == stepsBeforeDrug2:
            patient.addPrescription('grimpex')
        numVirusesEachStep.append(patient.update())
        numResistantVirusesEachStep.append(patient.getResistPop(['grimpex']))

    #print numVirusesEachStep[-1]
    
    return (numVirusesEachStep, numResistantVirusesEachStep)

#Testing START
simulationTwoDrugsDelayedTreatment(100, 1000, 0.1, 0.05, {'guttogonal': False, 'grimpex': False}, 0.005, 200)
#Testing END