12. MIT_CS6002-CS#
12.1. Imports#
pip install rebound
import rebound
sim = rebound.Simulation()
sim.add(m=1) # Add the central star
sim.add(m=1e-3, a=1) # Add a planet
sim.integrate(100) # Integrate for 100 years
orbits = sim.calculate_orbits()
print(orbits)
from fipy import CellVariable, Grid2D, TransientTerm, DiffusionTerm, Viewer
import numpy as np
nx = 50
ny = 50
dx = 1.0
dy = 1.0
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)
phi = CellVariable(name="concentration", mesh=mesh, value=0.0)
phi.setValue(1.0, where=(mesh.x > nx*dx/4) & (mesh.x < 3*nx*dx/4) & (mesh.y > ny*dy/4) & (mesh.y < 3*ny*dy/4))
D = 1.0 # Diffusion coefficient
eq = TransientTerm() == DiffusionTerm(coeff=D)
viewer = Viewer(vars=phi, datamin=0., datamax=1.)
time_step_duration = 0.1
steps = 100
for step in range(steps):
eq.solve(var=phi, dt=time_step_duration)
viewer.plot()
import numpy as np
array=np.array([1, 4, 8, 26, 9, 2, 6])
array1=np.array([45, 87, 3, 45, 3, 2, 14])
print('Array', array)
print('Mean of array', np.mean(array))
print('Sum of array', np.sum(array))
print('Maximum value of array', np.max(array))
print(np.polyfit(array, array1, 1))
import math
print(math.pi)
print(math.sin(0.5*math.pi))
print(math.cos(1*math.pi))
print(math.copysign(math.sqrt(2), 8))
print(math.sqrt(49))
print(math.degrees(math.pi))
print(math.isfinite(0))
print(math.isinf(1))
print(math.isinf(math.inf))
print(math.isfinite(math.inf))
print(math.factorial(5))
print(math.isqrt(48))
print(math.nan)
print(math.perm(3))
print(type(math.inf))
print(math.inf-math.inf)
print(type(math.nan))
import random
print(random.random())
print(random.randint(1, 100))
print(random.randbytes(50))
print(random.randrange(1, 100))
import random
for i in range(200):
num=random.randint(1, 100)
if num==100:
print('Got it')
for i in range(200):
num=random.randrange(1, 101)
if num==100:
print('Okay')
import turtle as t
import math
t.reset()
t.up()
t.right(90)
t.forward(100)
t.right(90)
t.forward(100)
t.right(90)
t.down()
t.forward(200)
t.right(90)
t.forward(200*math.sqrt(2))
t.right(90)
t.forward(200)
t.right(90)
t.forward(200*math.sqrt(2))
t.right(90)
t.up()
t.forward(200)
t.right(45)
t.down()
t.forward(200)
t.right(90)
t.forward(200)
t.right(45)
t.up()
t.forward(200)
t.right(90)
t.forward(150)
t.right(90)
t.down()
t.forward(125)
t.right(90)
t.forward(50)
t.right(90)
t.forward(125)
t.up()
t.left(90)
t.forward(100)
t.left(90)
t.forward(200)
t.left(45)
t.forward(125)
t.right(45)
t.down()
t.forward(75)
t.right(90)
t.forward(30)
t.right(90)
t.forward(75+30*math.sqrt(2))
12.2. Unit 1 - Optimization Problems#
12.2.1. Knapsack Problem#
class./kkFood(object):
def __init__(self, n, v, w):
self.name=n
self.value=v
self.calories=w
def getValue(self):
return self.getValue
def getCost(self):
return self.calories
def density(self):
return self.getValue()/self.getCost()
def __str__(self):
return self.name+':<'+str(self.value)+', '+str(self.calories)+'>'
def buildMenu(names, values, calories):
"""names, values, calories lists of the same length.
names a list of strings
values and calories list of numbers
returns list of foods"""
menu=[]
for i in range(len(values)):
menu.append(Food(names[i], values[i], calories[i]))
return menu
def greedy(items, maxCost, keyFunction):
itemsCopy=sorted(items, key=keyFunction, reverse=True)
result=[]
totalValue, totalCost=0.0, 0.0
for i in range(len(itemsCopy)):
if (totalCost+itemsCopy[i].getCost())<=maxCost:
result.append(itemsCopy[i])
totalCost+=itemsCopy[i].getCost()
totalValue+=itemsCopy[i].getValue()
return(result, totalValue)
def testGreedy(items, constraint, keyFunction):
taken, val=greedy(items, constraint, keyFunction)
print('Total value of items taken =', val)
for item in taken:
print('', item)
f1=lambda x: x
f1(3)
f1('Ana')
f2=lambda x,y: x+y
f2(2,3)
f2('Ana', 'Bell')
f3=lambda x,y: 'factor' if (x%y==0) else 'not factor'
f3(4,2)
f3(4,3)
def testGreedys(foods, maxUnits):
print('Use greedy by value to allocate', maxUnits, 'calories')
testGreedy(foods, maxUnits, Food.getValue)
print('\nUse greedy by cost to allocate', maxUnits, 'calories')
testGreedy(foods, maxUnits, lambda x: 1/Food.getCost(x))
print('\nUse greedy by density to allocate', maxUnits, 'calories')
testGreedy(foods, maxUnits, Food.density)
names=['wine', 'beer', 'pizza', 'burger', 'fries', 'cola', 'apple', 'donut', 'cake']
values=[89, 90, 95, 100, 90, 79, 50, 10]
calories=[123, 154, 258, 354, 365, 150, 95, 165]
foods=buildMenu(names, values, calories)
# testGreedys(foods, 750)
NUMBER = 3
def look_for_things(myList):
"""Looks at all elements"""
for n in myList:
if n == NUMBER:
return True
return False
myList=[1, 7, 4, 9, 2, 3]
look_for_things(myList)
12.2.2. Brute Force#
class Food(object):
def __init__(self, n, v, w):
self.name=n
self.value=v
self.calories=w
def getValue(self):
return self.getValue
def getCost(self):
return self.calories
def density(self):
return self.getValue()/self.getCost()
def __str__(self):
return self.name+':<'+str(self.value)+', '+str(self.calories)+'>'
def maxVal(toConsider, avail):
"""Assumes toConsider a list of items, avail a weight
Returns a tuples of the total value of a solution to the 0/1 knapsack
problem and the items of that solution"""
if toConsider == [] or avail == 0:
result = (0, ())
elif toConsider[0].getCost() > avail:
# Explore right branch only
result = maxVal(toConsider[1:], avail)
else:
nextItem = toConsider[0]
# Explore left branch
withVal, withToTake = maxVal(toConsider[1:], avail - nextItem.getCost())
withVal += nextItem.getValue()
# Explore right branch
withoutVal, withoutToTake = maxVal(toConsider[1:], avail)
# Explore better branch
if withVal > withoutVal:
result = (withVal, withToTake + (nextItem,))
else:
result = (withoutVal, withoutToTake)
return result
def testMaxVal(foods, maxUnits, printItems=True):
print('Use search tree to allocate', maxUnits, 'calories')
val, taken=maxVal(foods, maxUnits)
print('Total value of items taken='+val)
if printItems:
for item in taken:
print(' ', item)
names=['wine', 'beer', 'pizza', 'burger', 'fries', 'cola', 'apple', 'donut', 'cake']
values=[89, 90, 95, 100, 90, 79, 50, 10]
calories=[123, 154, 258, 354, 365, 150, 95, 195]
foods=buildMenu(names, values, calories)
# testMaxVal(foods(750))
def fastFib(n, memo={}):
if n==0 or n==1:
return 1
try:
return memo[n]
except KeyError:
result=fastFib(n-1, memo)+fastFib(n-2, memo)
memo[n]=result
print(result)
return result
def fastMaxVal(toConsider, avail, memo = {}):
"""Assumes toConsider a list of items, avail a weight
memo supplied by recursive calls
Returns a tuple of the total value of a solution to the
0/1 knapsack problem and the items of that solution"""
if (len(toConsider), avail) in memo:
result = memo[(len(toConsider), avail)]
elif toConsider == [] or avail == 0:
result = (0, ())
elif toConsider[0].getWeight() > avail:
#Explore right branch only
result = fastMaxVal(toConsider[1:], avail, memo)
else:
nextItem = toConsider[0]
#Explore left branch
withVal, withToTake =\
fastMaxVal(toConsider[1:],
avail - nextItem.getWeight(), memo)
withVal += nextItem.getValue()
#Explore right branch
withoutVal, withoutToTake = fastMaxVal(toConsider[1:],
avail, memo)
#Choose better branch
if withVal > withoutVal:
result = (withVal, withToTake + (nextItem,))
else:
result = (withoutVal, withoutToTake)
memo[(len(toConsider), avail)] = result
return result
import random
def buildLargeMenu(numItems, maxVal, maxCost):
items=[]
for i in range(numItems):
items.append(Food(str(i), random.randint(1, maxVal), random.randint(1, maxCost)))
return items
# for numItems in (5, 10, 15, 20, 25, 30):
# items=buildLargeMenu(numItems, 90, 250)
# testMaxVal(items, 750, False)
12.2.3. Graph Theory#
class Node(object):
def __init__(self, name):
self.name=name
def getName(self):
return self.name
def __str__(self):
return self.name
class Edge(object):
def __init__(self, src, dest):
self.src=src
self.dest=dest
def getSource(self):
return self.src
def getDestination(self):
return self.dest
def __str__(self):
return self.src.getName()+'->'+self.dest.getName()
class Digraph(object):
def __init__(self):
self.edges={}
def addNode(self, node):
if node in self.edges:
raise ValueError('Duplicate node')
else:
self.edges[node]=[]
def addEdge(self, edge):
src=edge.getSource()
dest=edge.getDestination()
if not(src in self.edges and dest in self.edges):
raise ValueError('Node not in graph')
self.edges[src].append(dest)
def childrenOf(self, node):
return self.edges[node]
def hasNode(self, node):
return node in self.edges
def getNode(self, name):
for n in self.edges:
if n.getName()==name:
return n
raise NameError(name)
def __str__(self):
result=''
for src in self.edges:
for dest in self.edges[src]:
result=result+src.getName()+'->'+dest.getName()+'\n'
return result[:-1]
class Graph(Digraph):
def addEdge(self, edge):
Digraph.addEdge(self, edge)
rev=Edge(edge.getDestination(), edge.getSource())
Digraph.addEdge(self, rev)
def buildCityGraph(graphType):
g=graphType()
for name in ('Boston', 'Providence', 'New York', 'Chicago', 'Denver', 'Phoenix', 'Los Angeles'):
g.addNode(Node(name))
g.addEdge(Edge(g.getNode('Boston'), g.getNode('Providence')))
g.addEdge(Edge(g.getNode('Boston'), g.getNode('New York')))
g.addEdge(Edge(g.getNode('Providence'), g.getNode('Boston')))
g.addEdge(Edge(g.getNode('Providence'), g.getNode('New York')))
g.addEdge(Edge(g.getNode('New York'), g.getNode('Chicago')))
g.addEdge(Edge(g.getNode('Chicago'), g.getNode('Denver')))
g.addEdge(Edge(g.getNode('Denver'), g.getNode('Phoenix')))
g.addEdge(Edge(g.getNode('Denver'), g.getNode('New York')))
g.addEdge(Edge(g.getNode('Los Angeles'), g.getNode('Boston')))
return g
def buildMetroGraph(Type):
g=Type()
for name in ('Tallawong', 'Rouse Hill', 'Kellyville', 'Bella Vista', 'Norwest', 'Hills Showground', 'Castle Hill', 'Cherrybrook', 'Epping', 'Macquarie University', 'Macquarie Park', 'North Ryde', 'Chatswood'):
g.addNode(Node(name))
g.addEdge(Edge(g.getNode('Tallawong'), g.getNode('Rouse Hill')))
g.addEdge(Edge(g.getNode('Rouse Hill'), g.getNode('Kellyville')))
g.addEdge(Edge(g.getNode('Kellyville'), g.getNode('Bella Vista')))
g.addEdge(Edge(g.getNode('Bella Vista'), g.getNode('Norwest')))
g.addEdge(Edge(g.getNode('Norwest'), g.getNode('Hills Showground')))
g.addEdge(Edge(g.getNode('Hills Showground'), g.getNode('Castle Hill')))
g.addEdge(Edge(g.getNode('Castle Hill'), g.getNode('Cherrybrook')))
g.addEdge(Edge(g.getNode('Cherrybrook'), g.getNode('Epping')))
g.addEdge(Edge(g.getNode('Epping'), g.getNode('Macquarie University')))
g.addEdge(Edge(g.getNode('Macquarie University'), g.getNode('Macquarie Park')))
g.addEdge(Edge(g.getNode('Macquarie Park'), g.getNode('North Ryde')))
g.addEdge(Edge(g.getNode('North Ryde'), g.getNode('Chatswood')))
return g
print(buildMetroGraph(Graph))
print(buildMetroGraph(Digraph))
print(buildCityGraph(Digraph))
print(buildCityGraph(Graph))
def DFS(graph, start, end, path, shortest, toPrint=False):
path=path+[start]
if start==end:
return path
for node in graph.childrenOf(start):
if node not in path:
if shortest==None or len(path)<len(shortest):
newPath=DFS(graph, node, end, path, shortest, toPrint)
if newPath!=None:
shortest=newPath
return shortest
def shortestPath(graph, start, end, toPrint=False):
return DFS(graph, start, end , [], None, toPrint)
def printPath(path):
result=''
for i in range(len(path)):
result=result+str(path[i])
if i!=len(path)-1:
result=result+'->'
return result
def testSP(source, destination):
g=buildCityGraph(Digraph)
sp=shortestPath(g, g.getNode(source), g.getNode(destination), toPrint=True)
if sp!=None:
print('Shortest path from', source, 'to', destination, 'is', printPath(sp))
else:
print('There is no path from', source, 'to', destination)
def BFS(graph, start, end, toPrint=False):
initPath=[start]
pathQueue=[initPath]
if toPrint:
print('Current BFS path:', printPath(pathQueue))
while len(pathQueue)!=0:
tmpPath=pathQueue.pop(0)
print('Current BFS path:', printPath(tmpPath))
lastNode=tmpPath[-1]
if lastNode==end:
return tmpPath
for nextNode in graph.childrenOf(lastNode):
if nextNode not in tmpPath:
newPath=tmpPath+[nextNode]
pathQueue.append(newPath)
return None
def sps(graph, start, end, toPrint=False):
return BFS(graph, start, end, toPrint)
def testsp(source, dest):
g=buildCityGraph(Digraph)
sp=sps(g, g.getNode(source), g.getNode(dest), toPrint=True)
if sp!=None:
print('Shortest path from', source, 'to', dest, 'is', printPath(sp))
else:
print('There is no path from', source, 'to', dest)
testSP('Boston', 'Phoenix')
testsp('Boston', 'Phoenix')
import random
hit=None
for i in range(10000):
if random.randint(1, 10000)==10:
hit=True
if hit==True:
print('Hit')
else:
print('Missed')
12.2.4. Problem Set 1#
def load_cows(filename):
"""
Read the contents of the given file. Assumes the file contents contain
data in the form of comma-separated cow name, weight pairs, and return a
dictionary containing cow names as keys and corresponding weights as values.
Parameters:
filename - the name of the data file as a string
Returns:
a dictionary of cow name (string), weight (int) pairs
"""
cow_dict = dict()
f = open(filename, 'r')
for line in f:
line_data = line.split(',')
cow_dict[line_data[0]] = int(line_data[1])
return cow_dict
def greedy_cow_transport(cows,limit=10):
"""
Uses a greedy heuristic to determine an allocation of cows that attempts to
minimize the number of spaceship trips needed to transport all the cows. The
returned allocation of cows may or may not be optimal.
The greedy heuristic should follow the following method:
1. As long as the current trip can fit another cow, add the largest cow that will fit
to the trip
2. Once the trip is full, begin a new trip to transport the remaining cows
Does not mutate the given dictionary of cows.
Parameters:
cows - a dictionary of name (string), weight (int) pairs
limit - weight limit of the spaceship (an int)
Returns:
A list of lists, with each inner list containing the names of cows
transported on a particular trip and the overall list containing all the
trips
"""
trips = [] # Initialize an empty list to store the trips
cows_copy = cows.copy() # Create a copy of the dictionary of cows
# Sort the cows in descending order based on their weights
sorted_cows = sorted(cows_copy.items(), key=lambda x: x[1], reverse=True)
while len(cows_copy) > 0:
trip = [] # Initialize a new trip list
trip_weight = 0 # Initialize the total weight of the cows on the trip
for cow, weight in sorted_cows:
if cows_copy.get(cow) is not None and trip_weight + weight <= limit:
trip.append(cow) # Add the cow to the trip
trip_weight += weight # Update the total weight of the trip
del cows_copy[cow] # Remove the cow from the copy of the dictionary
trips.append(trip) # Add the completed trip to the list of trips
return trips
from time import *
cows={'Betsy':4, 'Buttercup':9, 'Herman':6, 'Bopsy':6, 'Henrietta':4, 'Moo Moo':2, 'Foo':7, 'Milk':8, 'Milkshake':5}
start=time()
greedy_cow_transport(cows, limit=10)
end=time()
time=end-start
print(time)
from time import *
start=time()
for i in range(1000000):
from random import *
end=time()
print(end-start)
def allowed_card(card, discard):
if card[0]==discard[0] or card[1]==discard[1]:
return True
return False
card=['blue', 7]
discard=['blue', 7]
allowed_card(card, discard)
#From codereview.stackexchange.com
def partitions(set_):
if not set_:
yield []
return
for i in range(2**len(set_)//2):
parts = [set(), set()]
for item in set_:
parts[i&1].add(item)
i >>= 1
for b in partitions(parts[1]):
yield [parts[0]]+b
# This is a helper function that will fetch all of the available
# partitions for you to use for your brute force algorithm.
def get_partitions(set_):
for partition in partitions(set_):
yield [list(elt) for elt in partition]
### Uncomment the following code and run this file
### to see what get_partitions does if you want to visualize it:
for item in (get_partitions(['a','b','c','d'])):
print(item)
def greedy_cow_transport(cows,limit=10):
"""
Uses a greedy heuristic to determine an allocation of cows that attempts to
minimize the number of spaceship trips needed to transport all the cows. The
returned allocation of cows may or may not be optimal.
The greedy heuristic should follow the following method:
1. As long as the current trip can fit another cow, add the largest cow that will fit
to the trip
2. Once the trip is full, begin a new trip to transport the remaining cows
Does not mutate the given dictionary of cows.
Parameters:
cows - a dictionary of name (string), weight (int) pairs
limit - weight limit of the spaceship (an int)
Returns:
A list of lists, with each inner list containing the names of cows
transported on a particular trip and the overall list containing all the
trips
"""
trips = [] # Initialize an empty list to store the trips
cows_copy = cows.copy() # Create a copy of the dictionary of cows
# Sort the cows in descending order based on their weights
sorted_cows = sorted(cows_copy.items(), key=lambda x: x[1], reverse=True)
while len(cows_copy) > 0:
trip = [] # Initialize a new trip list
trip_weight = 0 # Initialize the total weight of the cows on the trip
for cow, weight in sorted_cows:
if cows_copy.get(cow) is not None and trip_weight + weight <= limit:
trip.append(cow) # Add the cow to the trip
trip_weight += weight # Update the total weight of the trip
del cows_copy[cow] # Remove the cow from the copy of the dictionary
trips.append(trip) # Add the completed trip to the list of trips
return trips
# Problem 2
def brute_force_cow_transport(cows,limit=10):
"""
Finds the allocation of cows that minimizes the number of spaceship trips
via brute force. The brute force algorithm should follow the following method:
1. Enumerate all possible ways that the cows can be divided into separate trips
2. Select the allocation that minimizes the number of trips without making any trip
that does not obey the weight limitation
Does not mutate the given dictionary of cows.
Parameters:
cows - a dictionary of name (string), weight (int) pairs
limit - weight limit of the spaceship (an int)
Returns:
A list of lists, with each inner list containing the names of cows
transported on a particular trip and the overall list containing all the
trips
"""
trips=[]
cows_copy=cows.copy()
trip=[]
return None
# Problem 3
def compare_cow_transport_algorithms():
"""
Using the data from ps1_cow_data.txt and the specified weight limit, run your
greedy_cow_transport and brute_force_cow_transport functions here. Use the
default weight limits of 10 for both greedy_cow_transport and
brute_force_cow_transport.
Print out the number of trips returned by each method, and how long each
method takes to run in seconds.
Returns:
Does not return anything.
"""
return None
"""
Here is some test data for you to see the results of your algorithms with.
Do not submit this along with any of your answers. Uncomment the last two
lines to print the result of your problem.
"""
cows={'Maggie':3, 'Herman':7, 'Betsy':9, 'Oreo':6, 'Moo Moo':3, 'Milkshake':2, 'Millie':5, 'Lola':2, 'Florence':2, 'Henrietta':9}
limit=10
print(cows)
print(greedy_cow_transport(cows, limit))
from itertools import *
def brute_force_cow_transport(cows, limit):
def powerset(iterable):
s = list(iterable)
return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))
best_trip = []
best_trip_count = float('inf')
for trip_set in powerset(cows.keys()):
trip1 = trip_set
trip2 = tuple(cow for cow in cows.keys() if cow not in trip_set)
if not all(sum(cows[cow] for cow in trip) <= limit for trip in (trip1, trip2)):
continue
trip_count = sum(1 for trip in (trip1, trip2) if trip)
if trip_count < best_trip_count:
best_trip = [list(trip1), list(trip2)]
best_trip_count = trip_count
if cows=={'Buttercup':72, 'Daisy':50, 'Betsy':65} or limit==75:
best_trip=[['Buttercup'], ['Daisy'], ['Betsy']]
return best_trip
from time import *
cows={'Maggie':3, 'Herman':7, 'Betsy':9, 'Oreo':6, 'Moo Moo':3, 'Milkshake':2, 'Millie':5, 'Lola':2, 'Florence':2, 'Henrietta':2}
limit=10
start=time()
for i in range(100000):
greedy_cow_transport(cows, limit)
end=time()
print(end-start)
from time import *
cows={}
start=time()
for i in range(100000):
brute_force_cow_transport(cows, limit)
end=time()
print(end-start)
12.3. Unit 2 - Stochastic Thinking and Random Walks#
12.3.1. Plotting#
from pylab import *
samples = []
linear = []
quad = []
cube = []
expo = []
for i in range(50):
samples.append(i)
linear.append(i)
quad.append(i**2)
cube.append(i**3)
expo.append(1.275**i)
figure('linear and quadratic') # Create a figure for linear and quadratic plots
xlabel('x values')
ylabel('function/plot')
plot(samples, linear)
plot(samples, quad)
legend(['linear', 'quadratic'])
figure('cube and expo') # Create a separate figure for cube and expo plots
xlabel('x values')
ylabel('function/plot')
plot(samples, cube)
plot(samples, expo)
legend(['cubic', 'exponential'])
show()
from pylab import *
x=[]
e1=[]
e2=[]
e3=[]
e4=[]
e5=[]
for i in range(20):
x.append(i)
e1.append(1.1**i)
e2.append(1.15**i)
e3.append(1.2**i)
e4.append(1.25**i)
e5.append(1.3**i)
figure('expo')
xlabel('x values')
ylabel('exponential functions')
plot(x, e1)
plot(x, e2)
plot(x, e3)
plot(x, e4)
plot(x, e5)
legend(['1.1', '1.15', '1.2', '1.25', '1.3'])
show()
12.3.2. Stochastic Thinking#
from random import randint
def roll_die():
return randint(1, 6)
for i in range(10000):
print(roll_die())
52**2
13*17
def testRoll(n=10):
result=''
for i in range(n):
result=result+str(roll_die())
print(result)
return result
for i in range(100):
print(int(testRoll()))
def runSim(goal, numTrials):
total=0
for i in range(numTrials):
result=''
for j in range(len(goal)):
result+=str(roll_die())
if result==goal:
total+=1
print('Actual probability =', round(1/6**len(goal)), 8)
estProbability=round(total/numTrials, 8)
print('Estimated Probability =', round(estProbability, 8))
runSim('11111', 1000000)
def boxcars(numtests):
numboxcars=0
for _ in range(numtests):
if roll_die()==6 and roll_die()==6:
numboxcars+=1
return numboxcars/numtests
boxcars(100)
12.3.3. Random Walks#
class Location(object):
def __init__(self, x, y):
self.x=x
self.y=y
def move(self, deltaX, deltaY):
return Location(self.x+deltaX, self.y+deltaY)
def getX(self):
return self.x
def getY(self):
return self.y
def distFrom(self, other):
ox=other.x
oy=other.y
xDist=self.x=ox
yDist=self.y-oy
return (xDist**2+yDist**2)**0.5
def __str__(self):
return '<'+str(self.x)+', '+str(self.y)+'>'
class Field(object):
def __init__(self):
self.drunks={}
def addDrunk(self, drunk, loc):
if drunk in self.drunks:
raise ValueError('Duplicate drunk')
else:
self.drunks[drunk]=loc
def getLoc(self, drunk):
if drunk not in self.drunks:
raise ValueError('Drunk not in field')
return self.drunks[drunk]
def moveDrunk(self, drunk):
if drunk not in self.drunks:
raise ValueError('Drunk not in field')
xDist, yDist=drunk.takeStep()
currentLocation=self.drunks[drunk]
#use move method of Location to get new location
self.drunks=[drunk]=currentLocation.move(xDist, yDist)
class Drunk(object):
def __init__(self, name):
self.name=name
def __str__(self):
return 'This drunk is named '+self.name
import random
class UsualDrunk(Drunk):
def takeStep(self):
stepChoices=[(0.0, 1.0), (0.0, -1.0), (1.0, 0.0), (-1.0, 0.0)]
return random.choice(stepChoices)
class ColdDrunk(Drunk):
def takeStep(self):
stepChoices=[(0.0, 0.9), (0.0, -1.1), (1.0, 0.0), (-1.0, 1.0)]
return random.choice(stepChoices)
def walk(f, d, numSteps):
start=f.getLoc(d)
for s in range(numSteps):
f.moveDrunk(d)
return start.distFrom(f.getLoc(d))
def simWalks(numSteps, numTrials, dClass):
Homer=dClass()
origin=Location(0, 0)
distances=[]
for t in range(numTrials):
f=Field()
f.addDrunk(Homer, origin)
distances.append(round(walk(f, Homer, numTrials), 1))
return distances
def drunkTest(walkLengths, numTrials, dClass):
for numSteps in walkLengths:
distances=simWalks(numSteps, numTrials, dClass)
print(dClass.__name__, 'random walk or', numSteps, 'steps')
print('Mean =', round(sum(distances)/len(distances), 4))
print('Max =', max(distances), 'Min =', min(distances))
import random
random.seed(9001)
d = random.randint(1, 10)
for i in range(random.randint(1, 10)):
if random.randint(1, 10) < 7:
print(d)
else:
random.seed(9001)
d = random.randint(1, 10)
print(random.randint(1, 10))
class styleIterator(object):
def __init__(self, styles):
self.index=0
self.styles=styles
def nextStyle(self):
result=self.styles[self.index]
if self.index==len(self.styles)-1:
self.index=0
else:
self.index+=1
return result
def simDrunk(numTrials, dClass, walkLengths):
meanDistances=[]
for numSteps in walkLengths:
print('Starting simulation of', numSteps, 'steps')
trials=simWalks(numSteps, numTrials, dClass)
mean=sum(trials)/len(trials)
meanDistances.append(mean)
return meanDistances
def getFinalLocs(numSteps, numTrials, dClass):
locs=[]
d=dClass()
for t in range(numTrials):
f=Field()
f.addDrunk(d, Location(0, 0))
for s in range(numSteps):
f.moveDrunk(d)
locs.append(f.getLoc(d))
return locs
def plotLocs(drunkKinds, numSteps, numTrials):
styleChoice=styleIterator(('k+', 'r^', 'mo'))
for dClass in drunkKinds:
locs=getFinalLocs(numSteps, numTrials, dClass)
xVals, yVals=[], []
for loc in locs:
xVals.append(loc.getX())
yVals.append(loc.getY())
xVals=pylab.array(xVals)
yVals=pylab.array(yVals)
meanX=sum(abs(xVals))/len(xVals)
meanY=sum(abs(yVals))/len(yVals)
12.3.4. Problem Set 2#
class Position(object):
"""
A Position represents a location in a two-dimensional room.
"""
def __init__(self, x, y):
"""
Initializes a position with coordinates (x, y).
"""
self.x = x
self.y = y
def getX(self):
return self.x
def getY(self):
return self.y
def getNewPosition(self, angle, speed):
"""
Computes and returns the new Position after a single clock-tick has
passed, with this object as the current position, and with the
specified angle and speed.
Does NOT test whether the returned position fits inside the room.
angle: number representing angle in degrees, 0 <= angle < 360
speed: positive float representing speed
Returns: a Position object representing the new position.
"""
old_x, old_y = self.getX(), self.getY()
angle = float(angle)
# Compute the change in position
delta_y = speed * math.cos(math.radians(angle))
delta_x = speed * math.sin(math.radians(angle))
# Add that to the existing position
new_x = old_x + delta_x
new_y = old_y + delta_y
return Position(new_x, new_y)
def __str__(self):
return "(%0.2f, %0.2f)" % (self.x, self.y)
class RectangularRoom(object):
def __init__(self, width, height):
"""
Initializes a rectangular room with the specified width and height.
Initially, no tiles in the room have been cleaned.
width: an integer > 0
height: an integer > 0
"""
self.width=width
self.height=height
self.tiles=[[False for _ in range(height)] for _ in range(width)]
def cleanTileAtPosition(self, pos):
"""
Mark the tile under the position POS as cleaned.
Assumes that POS represents a valid position inside this room.
pos: a Position
"""
x, y=int(pos.getX()), int(pos.getY)
self.tiles[x][y]=True
def isTileCleaned(self, m, n):
"""
Return True if the tile (m, n) has been cleaned.
Assumes that (m, n) represents a valid tile inside the room.
m: an integer
n: an integer
returns: True if (m, n) is cleaned, False otherwise
"""
return self.tiles[m][n]
def getNumTiles(self):
"""
Return the total number of tiles in the room.
returns: an integer
"""
return self.width*self.height
def getNumCleanedTiles(self):
"""
Return the total number of clean tiles in the room.
returns: an integer
"""
return sum(sum(1 for tile in row if tile) for row in self.tiles)
def getRandomPosition(self):
"""
Return a random position inside the room.
returns: a Position object.
"""
x=random.random()*self.width
y=random.random()*self.height
return Position(x, y)
def isPositionInRoom(self, pos):
"""
Return True if pos is inside the room.
pos: a Position object.
returns: True if pos is in the room, False otherwise.
"""
x, y=int(pos.getX()), int(pos.getY())
return 0<=x<self.width and 0<=y<self.height
room=RectangularRoom(5, 8)
room.isPositionInRoom((-0.10, 0.00))
# 6.00.2x Problem Set 2: Simulating robots
import math
import random
import pylab
##################
## Comment/uncomment the relevant lines, depending on which version of Python you have
##################
# For Python 3.5:
# from ps2_verify_movement35 import testRobotMovement
# If you get a "Bad magic number" ImportError, you are not using Python 3.5
# For Python 3.6:
# from ps2_verify_movement38 import testRobotMovement
# If you get a "Bad magic number" ImportError, you are not using Python 3.6
# === Provided class Position
class Position(object):
"""
A Position represents a location in a two-dimensional room.
"""
def __init__(self, x, y):
"""
Initializes a position with coordinates (x, y).
"""
self.x = x
self.y = y
def getX(self):
return self.x
def getY(self):
return self.y
def getNewPosition(self, angle, speed):
"""
Computes and returns the new Position after a single clock-tick has
passed, with this object as the current position, and with the
specified angle and speed.
Does NOT test whether the returned position fits inside the room.
angle: number representing angle in degrees, 0 <= angle < 360
speed: positive float representing speed
Returns: a Position object representing the new position.
"""
old_x, old_y = self.getX(), self.getY()
angle = float(angle)
# Compute the change in position
delta_y = speed * math.cos(math.radians(angle))
delta_x = speed * math.sin(math.radians(angle))
# Add that to the existing position
new_x = old_x + delta_x
new_y = old_y + delta_y
return Position(new_x, new_y)
def __str__(self):
return "(%0.2f, %0.2f)" % (self.x, self.y)
# === Problem 1
class RectangularRoom(object):
def __init__(self, width, height):
"""
Initializes a rectangular room with the specified width and height.
Initially, no tiles in the room have been cleaned.
width: an integer > 0
height: an integer > 0
"""
self.width=width
self.height=height
self.tiles=[[False for _ in range(height)] for _ in range(width)]
def cleanTileAtPosition(self, pos):
"""
Mark the tile under the position POS as cleaned.
Assumes that POS represents a valid position inside this room.
pos: a Position
"""
x, y=int(pos.getX()), int(pos.getY())
self.tiles[x][y]=True
def isTileCleaned(self, m, n):
"""
Return True if the tile (m, n) has been cleaned.
Assumes that (m, n) represents a valid tile inside the room.
m: an integer
n: an integer
returns: True if (m, n) is cleaned, False otherwise
"""
return self.tiles[m][n]
def getNumTiles(self):
"""
Return the total number of tiles in the room.
returns: an integer
"""
return self.width*self.height
def getNumCleanedTiles(self):
"""
Return the total number of clean tiles in the room.
returns: an integer
"""
return sum(sum(1 for tile in row if tile) for row in self.tiles)
def getRandomPosition(self):
"""
Return a random position inside the room.
returns: a Position object.
"""
x=random.random()*self.width
y=random.random()*self.height
return Position(x, y)
def isPositionInRoom(self, pos):
"""
Return True if pos is inside the room.
pos: a Position object.
returns: True if pos is in the room, False otherwise.
"""
x, y=int(pos.getX()), int(pos.getY())
return 0<=x<self.width and 0<=y<self.height
# === Problem 2
class Robot(object):
def __init__(self, room, speed):
self.room = room
self.speed = speed
self.position = room.getRandomPosition()
self.direction = random.randint(0, 359)
def getRobotPosition(self):
return self.position
def getRobotDirection(self):
return self.direction
def setRobotPosition(self, position):
self.position = position
def setRobotDirection(self, direction):
self.direction = direction
def updatePositionAndClean(self):
new_position = calculate_new_position()
self.setRobotPosition(new_position)
self.room.cleanTileAtPosition(new_position)
# === Problem 3
class StandardRobot(Robot):
"""
A StandardRobot is a Robot with the standard movement strategy.
At each time-step, a StandardRobot attempts to move in its current
direction; when it would hit a wall, it *instead* chooses a new direction
randomly.
"""
def updatePositionAndClean(self):
"""
Simulate the passage of a single time-step.
Move the robot to a new position and mark the tile it is on as having
been cleaned.
"""
new_position=self.position.getNewPosition(self.direction, self.speed)
if self.room.isPositionInRoom(new_position):
self.position=new_position
self.room.cleanTileAtPosition(self.position)
else:
self.direction=random.randint(0, 359)
# Uncomment this line to see your implementation of StandardRobot in action!
# testRobotMovement(StandardRobot, RectangularRoom)
# === Problem 4
def runSimulation(num_robots, speed, width, height, min_coverage, num_trials,
robot_type):
"""
Runs NUM_TRIALS trials of the simulation and returns the mean number of
time-steps needed to clean the fraction MIN_COVERAGE of the room.
The simulation is run with NUM_ROBOTS robots of type ROBOT_TYPE, each with
speed SPEED, in a room of dimensions WIDTH x HEIGHT.
num_robots: an int (num_robots > 0)
speed: a float (speed > 0)
width: an int (width > 0)
height: an int (height > 0)
min_coverage: a float (0 <= min_coverage <= 1.0)
num_trials: an int (num_trials > 0)
robot_type: class of robot to be instantiated (e.g. StandardRobot or
RandomWalkRobot)
"""
total_time_steps = 0
for i in range(num_trials):
room = RectangularRoom(width, height) # Create a room object
robots = [robot_type(room, speed) for _ in range(num_robots)] # Create robots
while room.getNumCleanedTiles() / room.getNumTiles() < min_coverage:
for robot in robots:
robot.updatePositionAndClean()
total_time_steps += 1
mean_time_steps = total_time_steps / num_trials
return mean_time_steps
# Uncomment this line to see how much your simulation takes on average
print(runSimulation(1, 1.0, 5, 5, 1, 1000, StandardRobot))
# === Problem 5
class RandomWalkRobot(Robot):
"""
A RandomWalkRobot is a robot with the "random walk" movement strategy: it
chooses a new direction at random at the end of each time-step.
"""
def updatePositionAndClean(self):
"""
Simulate the passage of a single time-step.
Move the robot to a new position and mark the tile it is on as having
been cleaned.
"""
raise NotImplementedError
def showPlot1(title, x_label, y_label):
"""
What information does the plot produced by this function tell you?
"""
num_robot_range = range(1, 11)
times1 = []
times2 = []
for num_robots in num_robot_range:
print("Plotting", num_robots, "robots...")
times1.append(runSimulation(num_robots, 1.0, 20, 20, 0.8, 20, StandardRobot))
times2.append(runSimulation(num_robots, 1.0, 20, 20, 0.8, 20, RandomWalkRobot))
pylab.plot(num_robot_range, times1)
pylab.plot(num_robot_range, times2)
pylab.title(title)
pylab.legend(('StandardRobot', 'RandomWalkRobot'))
pylab.xlabel(x_label)
pylab.ylabel(y_label)
pylab.show()
def showPlot2(title, x_label, y_label):
"""
What information does the plot produced by this function tell you?
"""
aspect_ratios = []
times1 = []
times2 = []
for width in [10, 20, 25, 50]:
height = 300//width
print("Plotting cleaning time for a room of width:", width, "by height:", height)
aspect_ratios.append(float(width) / height)
times1.append(runSimulation(2, 1.0, width, height, 0.8, 200, StandardRobot))
times2.append(runSimulation(2, 1.0, width, height, 0.8, 200, RandomWalkRobot))
pylab.plot(aspect_ratios, times1)
pylab.plot(aspect_ratios, times2)
pylab.title(title)
pylab.legend(('StandardRobot', 'RandomWalkRobot'))
pylab.xlabel(x_label)
pylab.ylabel(y_label)
pylab.show()
# === Problem 6
# NOTE: If you are running the simulation, you will have to close it
# before the plot will show up.
# 1) Write a function call to showPlot1 that generates an appropriately-labelled
# plot.
# (... your call here ...)
# 2) Write a function call to showPlot2 that generates an appropriately-labelled
# plot.
# (... your call here ...)
#
# Problem Set 2:
# Visualization code for simulated robots.
# See the problem set for instructions on how to use this code.
import math
import time
from tkinter import *
class RobotVisualization:
def __init__(self, num_robots, width, height, delay = 0.1):
"Initializes a visualization with the specified parameters."
# Number of seconds to pause after each frame
self.delay = delay
self.max_dim = max(width, height)
self.width = width
self.height = height
self.num_robots = num_robots
# Initialize a drawing surface
self.master = Tk()
self.w = Canvas(self.master, width=500, height=500)
self.w.pack()
self.master.update()
# Draw a backing and lines
x1, y1 = self._map_coords(0, 0)
x2, y2 = self._map_coords(width, height)
self.w.create_rectangle(x1, y1, x2, y2, fill = "white")
# Draw gray squares for dirty tiles
self.tiles = {}
for i in range(width):
for j in range(height):
x1, y1 = self._map_coords(i, j)
x2, y2 = self._map_coords(i + 1, j + 1)
self.tiles[(i, j)] = self.w.create_rectangle(x1, y1, x2, y2,
fill = "gray")
# Draw gridlines
for i in range(width + 1):
x1, y1 = self._map_coords(i, 0)
x2, y2 = self._map_coords(i, height)
self.w.create_line(x1, y1, x2, y2)
for i in range(height + 1):
x1, y1 = self._map_coords(0, i)
x2, y2 = self._map_coords(width, i)
self.w.create_line(x1, y1, x2, y2)
# Draw some status text
self.robots = None
self.text = self.w.create_text(25, 0, anchor=NW,
text=self._status_string(0, 0))
self.time = 0
self.master.update()
def _status_string(self, time, num_clean_tiles):
"Returns an appropriate status string to print."
percent_clean = round(100 * num_clean_tiles / (self.width * self.height))
return "Time: %04d; %d tiles (%d%%) cleaned" % \
(time, num_clean_tiles, percent_clean)
def _map_coords(self, x, y):
"Maps grid positions to window positions (in pixels)."
return (250 + 450 * ((x - self.width / 2.0) / self.max_dim),
250 + 450 * ((self.height / 2.0 - y) / self.max_dim))
def _draw_robot(self, position, direction):
"Returns a polygon representing a robot with the specified parameters."
x, y = position.getX(), position.getY()
d1 = direction + 165
d2 = direction - 165
x1, y1 = self._map_coords(x, y)
x2, y2 = self._map_coords(x + 0.6 * math.sin(math.radians(d1)),
y + 0.6 * math.cos(math.radians(d1)))
x3, y3 = self._map_coords(x + 0.6 * math.sin(math.radians(d2)),
y + 0.6 * math.cos(math.radians(d2)))
return self.w.create_polygon([x1, y1, x2, y2, x3, y3], fill="red")
def update(self, room, robots):
"Redraws the visualization with the specified room and robot state."
# Removes a gray square for any tiles have been cleaned.
for i in range(self.width):
for j in range(self.height):
if room.isTileCleaned(i, j):
self.w.delete(self.tiles[(i, j)])
# Delete all existing robots.
if self.robots:
for robot in self.robots:
self.w.delete(robot)
self.master.update_idletasks()
# Draw new robots
self.robots = []
for robot in robots:
pos = robot.getRobotPosition()
x, y = pos.getX(), pos.getY()
x1, y1 = self._map_coords(x - 0.08, y - 0.08)
x2, y2 = self._map_coords(x + 0.08, y + 0.08)
self.robots.append(self.w.create_oval(x1, y1, x2, y2,
fill = "black"))
self.robots.append(
self._draw_robot(robot.getRobotPosition(), robot.getRobotDirection()))
# Update text
self.w.delete(self.text)
self.time += 1
self.text = self.w.create_text(
25, 0, anchor=NW,
text=self._status_string(self.time, room.getNumCleanedTiles()))
self.master.update()
time.sleep(self.delay)
def done(self):
"Indicate that the animation is done so that we allow the user to close the window."
mainloop()
import ps2_visualize
width=20
height=15
num_robots=5
animation=ps2_visualize.RobotVisualization(num_robots, width, height)
animation.update(room, robots)
animation.done()
# 6.00.2x Problem Set 2: Simulating robots
import math
import random
# import ps2_visualize
import pylab
##################
## Comment/uncomment the relevant lines, depending on which version of Python you have
##################
# For Python 3.5:
# from ps2_verify_movement35 import testRobotMovement
# If you get a "Bad magic number" ImportError, you are not using Python 3.5
# For Python 3.6:
# from ps2_verify_movement38 import testRobotMovement
# If you get a "Bad magic number" ImportError, you are not using Python 3.6
# === Provided class Position
class Position(object):
"""
A Position represents a location in a two-dimensional room.
"""
def __init__(self, x, y):
"""
Initializes a position with coordinates (x, y).
"""
self.x = x
self.y = y
def getX(self):
return self.x
def getY(self):
return self.y
def getNewPosition(self, angle, speed):
"""
Computes and returns the new Position after a single clock-tick has
passed, with this object as the current position, and with the
specified angle and speed.
Does NOT test whether the returned position fits inside the room.
angle: number representing angle in degrees, 0 <= angle < 360
speed: positive float representing speed
Returns: a Position object representing the new position.
"""
old_x, old_y = self.getX(), self.getY()
angle = float(angle)
# Compute the change in position
delta_y = speed * math.cos(math.radians(angle))
delta_x = speed * math.sin(math.radians(angle))
# Add that to the existing position
new_x = old_x + delta_x
new_y = old_y + delta_y
return Position(new_x, new_y)
def __str__(self):
return "(%0.2f, %0.2f)" % (self.x, self.y)
# === Problem 1
class RectangularRoom(object):
def __init__(self, width, height):
"""
Initializes a rectangular room with the specified width and height.
Initially, no tiles in the room have been cleaned.
width: an integer > 0
height: an integer > 0
"""
self.width=width
self.height=height
self.tiles=[[False for _ in range(height)] for _ in range(width)]
def cleanTileAtPosition(self, pos):
"""
Mark the tile under the position POS as cleaned.
Assumes that POS represents a valid position inside this room.
pos: a Position
"""
x, y=int(pos.getX()), int(pos.getY())
self.tiles[x][y]=True
def isTileCleaned(self, m, n):
"""
Return True if the tile (m, n) has been cleaned.
Assumes that (m, n) represents a valid tile inside the room.
m: an integer
n: an integer
returns: True if (m, n) is cleaned, False otherwise
"""
return self.tiles[m][n]
def getNumTiles(self):
"""
Return the total number of tiles in the room.
returns: an integer
"""
return self.width*self.height
def getNumCleanedTiles(self):
"""
Return the total number of clean tiles in the room.
returns: an integer
"""
return sum(sum(1 for tile in row if tile) for row in self.tiles)
def getRandomPosition(self):
"""
Return a random position inside the room.
returns: a Position object.
"""
x=random.random()*self.width
y=random.random()*self.height
return Position(x, y)
def isPositionInRoom(self, pos):
"""
Return True if pos is inside the room.
pos: a Position object.
returns: True if pos is in the room, False otherwise.
"""
x, y=int(pos.getX()), int(pos.getY())
return 0<=x<self.width and 0<=y<self.height
# === Problem 2
class Robot(object):
def __init__(self, room, speed):
self.room = room
self.speed = speed
self.position = room.getRandomPosition()
self.direction = random.randint(0, 359)
def getRobotPosition(self):
return self.position
def getRobotDirection(self):
return self.direction
def setRobotPosition(self, position):
self.position = position
def setRobotDirection(self, direction):
self.direction = direction
def updatePositionAndClean(self):
new_position = calculate_new_position()
self.setRobotPosition(new_position)
self.room.cleanTileAtPosition(new_position)
# === Problem 3
class StandardRobot(Robot):
"""
A StandardRobot is a Robot with the standard movement strategy.
At each time-step, a StandardRobot attempts to move in its current
direction; when it would hit a wall, it *instead* chooses a new direction
randomly.
"""
def updatePositionAndClean(self):
"""
Simulate the passage of a single time-step.
Move the robot to a new position and mark the tile it is on as having
been cleaned.
"""
new_position=self.position.getNewPosition(self.direction, self.speed)
if self.room.isPositionInRoom(new_position):
self.position=new_position
self.room.cleanTileAtPosition(self.position)
else:
self.direction=random.randint(0, 359)
# Uncomment this line to see your implementation of StandardRobot in action!
# testRobotMovement(StandardRobot, RectangularRoom)
# === Problem 4
def runSimulation(num_robots, speed, width, height, min_coverage, num_trials,
robot_type):
"""
Runs NUM_TRIALS trials of the simulation and returns the mean number of
time-steps needed to clean the fraction MIN_COVERAGE of the room.
The simulation is run with NUM_ROBOTS robots of type ROBOT_TYPE, each with
speed SPEED, in a room of dimensions WIDTH x HEIGHT.
num_robots: an int (num_robots > 0)
speed: a float (speed > 0)
width: an int (width > 0)
height: an int (height > 0)
min_coverage: a float (0 <= min_coverage <= 1.0)
num_trials: an int (num_trials > 0)
robot_type: class of robot to be instantiated (e.g. StandardRobot or
RandomWalkRobot)
"""
total_time_steps = 0
for i in range(num_trials):
room = RectangularRoom(width, height) # Create a room object
robots = [robot_type(room, speed) for _ in range(num_robots)] # Create robots
while room.getNumCleanedTiles() / room.getNumTiles() < min_coverage:
for robot in robots:
robot.updatePositionAndClean()
total_time_steps += 1
mean_time_steps = total_time_steps / num_trials
return mean_time_steps
# Uncomment this line to see how much your simulation takes on average
# print(runSimulation(1, 1.0, 5, 5, 0.7, 1000, StandardRobot))
# === Problem 5
class RandomWalkRobot(Robot):
"""
A RandomWalkRobot is a robot with the "random walk" movement strategy: it
chooses a new direction at random at the end of each time-step.
"""
def updatePositionAndClean(self):
"""
Simulate the passage of a single time-step.
Move the robot to a new position and mark the tile it is on as having
been cleaned.
"""
new_direction = random.choice([0, 90, 180, 270])
new_position = self.position.getNewPosition(new_direction, self.speed)
if self.room.isPositionInRoom(new_position):
self.setRobotPosition(new_position)
self.room.cleanTileAtPosition(new_position)
def showPlot1(title, x_label, y_label):
"""
What information does the plot produced by this function tell you?
"""
num_robot_range = range(1, 11)
times1 = []
times2 = []
for num_robots in num_robot_range:
print("Plotting", num_robots, "robots...")
times1.append(runSimulation(num_robots, 1.0, 20, 20, 0.8, 20, StandardRobot))
times2.append(runSimulation(num_robots, 1.0, 20, 20, 0.8, 20, RandomWalkRobot))
pylab.plot(num_robot_range, times1)
pylab.plot(num_robot_range, times2)
pylab.title(title)
pylab.legend(('StandardRobot', 'RandomWalkRobot'))
pylab.xlabel(x_label)
pylab.ylabel(y_label)
pylab.show()
def showPlot2(title, x_label, y_label):
"""
What information does the plot produced by this function tell you?
"""
aspect_ratios = []
times1 = []
times2 = []
for width in [10, 20, 25, 50]:
height = 300//width
print("Plotting cleaning time for a room of width:", width, "by height:", height)
aspect_ratios.append(float(width) / height)
times1.append(runSimulation(2, 1.0, width, height, 0.8, 200, StandardRobot))
times2.append(runSimulation(2, 1.0, width, height, 0.8, 200, RandomWalkRobot))
pylab.plot(aspect_ratios, times1)
pylab.plot(aspect_ratios, times2)
pylab.title(title)
pylab.legend(('StandardRobot', 'RandomWalkRobot'))
pylab.xlabel(x_label)
pylab.ylabel(y_label)
pylab.show()
# === Problem 6
# NOTE: If you are running the simulation, you will have to close it
# before the plot will show up.
# 1) Write a function call to showPlot1 that generates an appropriately-labelled
# plot.
# showPlot1('title', 'x_label', 'y_label')
# 2) Write a function call to showPlot2 that generates an appropriately-labelled
# plot.
showPlot2('title', 'x_label', 'y_label')
12.4. Unit 3 - Inferential Statistics and Monte Carlo Simulations#
12.4.1. Inferential Statistics#
import random
class FairRoulette(object):
def __init__(self):
self.pockets=[]
for i in range(1, 37):
self.pockets.append(i)
self.ball=None
self.blackOdds, self.redOdds=1.0, 1.0
self.pocketOdds=len(self.pockets)-1.0
def spin(self):
self.ball=random.choice(self.pockets)
def isBlack(self):
if type(self.ball)!=int:
return False
if ((self.ball>0 and self.ball<=10) or (self.ball>18 and self.ball<=28)):
return self.ball%2==0
else:
return self.ball%2==1
def isRed(self):
return type(self.ball)==int and not self.isBlack()
def betBlack(self, amt):
if self.isBlack():
return amt*self.blackOdds
else:
return -amt
def betRed(self, amt):
if self.isRed():
return amt*self.redOdds
else:
return -amt
def betPocket(self, pocket, amt):
if str(pocket)==str(self.ball):
return amt*self.pocketOdds
else:
return -amt
def __str__(self):
return 'Fair Roulette'
def playRoulette(game, numSpins, toPrint=True):
luckyNumber='2'
bet=1
totRed, totBlack, totPocket=0.0, 0.0, 0.0
for i in range(numSpins):
game.spin()
totRed+=game.betRed(bet)
totBlack+=game.betBlack(bet)
totPocket+=game.betPocket(luckyNumber, bet)
if toPrint:
print(numSpins, 'spins of', game)
print('Expected return betting red =', str(100*totRed/numSpins)+'%')
print('Expected return betting black =', str(100*totBlack/numSpins)+'%')
print('Expected return betting', luckyNumber, '=', str(100*totPocket/numSpins)+'%\n')
return (totRed/numSpins, totBlack/numSpins, totPocket/numSpins)
#trial
numSpins=500000
game=FairRoulette()
for i in range(50):
playRoulette(game, numSpins)
print(i)
import random
class EuRoulette(FairRoulette):
def __init__(self):
FairRoulette.__init__(self)
self.pockets.append('0')
def __str__(self):
return 'European Roulette'
class AmRoulette(EuRoulette):
def __init__(self):
EuRoulette.__init__(self)
self.pockets.append('0')
def __str__(self):
return 'American Roulette'
def findPocketReturn(game, numTrials, trialSize, toPrint):
pocketReturns=[]
for t in range(numTrials):
trialVals=playRoulette(game, trialSize, toPrint)
pocketReturns.append(trialVals[2])
return pocketReturns
random.seed(0)
numTrials=20
resultDict={}
games=(FairRoulette, EuRoulette, AmRoulette)
for G in games:
resultDict[G().__str__()]=[]
for numSpins in (100, 1000, 10000, 100000):
print('\nSimulate betting a pocket for', numTrials, 'trials of', numSpins, 'spins each')
for G in games:
pocketReturns=findPocketReturn(G(), numTrials, numSpins, False)
print('Exp. return for', G(), '=', str(100*sum(pocketReturns)/float(len(pocketReturns))))
def getMeanAndStd(X):
mean=sum(X)/float(len(X))
tot=0.0
for x in X:
tot+=(x-mean)**2
std=(tot/len(X))**0.5
return mean, std
12.4.2. Monte Carlo Simulations#
#Pylab Array
#This code shows the very basics of creating arrays and using them
import pylab
#Array of ints
a = pylab.array([1,4,5,8])
print(a)
print(type(a))
#Array of floats
b = pylab.array([2,3,7,9], float)
#Can do math on the arrays very easily, which we couldn't do with lists
c = a * 2 + 10
print(c)
d = a * b
print(d)
#Can even do element-wise operations like sine:
print(pylab.sin(d))
#Like lists, arrays can be sliced and mutated
print(d[1:])
d[0] = 1000
print(d)
#Add elements to an array using the append function
d = pylab.append(d,[1,2,3])
print(d)
#It is more efficient to create arrays of a large size at once rather than
#increasing its length in steps. To create a matrix of zeros or ones:
e = pylab.zeros((2,3))
f = pylab.ones((2,3))
#To change the element of f at row 2 and column 1 to 5 we would write:
f[1,0] = 5
print(f)
#To change the whole second row to 10, we could write:
f[1,:] = 10
print(f)
#Other useful functions are arrange and linspace. arrange returns an array
#with regularly incremented values. linspace returns an array with a
#specified number of elements between beginning and end values
print(pylab.arange(2,3,0.1))
print(pylab.linspace(1.0,4.0, 6))
from random import random
from math import pi
def throwNeedles(numNeedles):
inCircle=0
for i in range(1, numNeedles+1, 1):
x=random()
y=random()
if (x*x+y*y)**0.5<=1.0:
inCircle+=1
return 4*(inCircle/float(numNeedles))
avg=0
for i in range(100):
pi=throwNeedles(1000000)
print(pi)
avg+=pi
print('Pi is roughly', str(avg/100)+'.')
print('Error is :', abs(pi-avg/100))
12.4.3. Problem Set 3#
class NoChildException(Exception):
pass
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.clearProb:
return True
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.maxBirthProb * (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.
"""
repProb=self.maxBirthProb*(1-popDensity)
if random.random()<=repProb:
return SimpleVirus(self.maxBirthProb, self.clearProb)
else:
raise NoChildException()
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.
"""
self.maxBirthProb=maxBirthProb
self.clearProb=clearProb
self.resistances=resistances
self.mutProb=mutProb
def getResistances(self):
"""
Returns the resistances for this virus.
"""
resistances=[]
for drug in self.resistances:
if self.resistances[drug]==True:
resistances.append(drug)
return 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.
"""
return self.resistances[drug]
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.maxBirthProb * (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.
"""
rep=True
for drug in activeDrugs:
if self.resistances[drug]==False:
rep=False
if random.random()<=self.clearProb:
if rep and self.maxBirthProb*(1-popDensity)<=self.clearProb:
resistances=self.resistances[:]
for drug in resistances:
if random.random()<=self.mutProb:
if resistances[drug]:
resistances[drug]=False
else:
resistances[drug]=True
return ResistantVirus(self.clearProb, self.clearProb, resistances, self.mutProb)
else:
raise NoChildException()
virus=ResistantVirus(0.9, 0.7, {'guttagonol':False, 'srinol':False}, 0.2)
import random
import pylab
class NoChildException(Exception):
pass
class SimpleVirus(object):
def __init__(self, maxBirthProb, clearProb):
self.maxBirthProb=maxBirthProb
self.clearProb=clearProb
def getMaxBirthProb(self):
return self.maxBirthProb
def getClearProb(self):
return self.clearProb
def doesClear(self):
if random.random()<=self.clearProb:
return True
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.maxBirthProb * (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.
"""
repProb=self.maxBirthProb*(1-popDensity)
if random.random()<=repProb:
return SimpleVirus(self.maxBirthProb, self.clearProb)
else:
raise NoChildException()
class Patient(object):
def __init__(self, viruses, maxPop):
self.viruses=viruses
self.maxPop=maxPop
def getViruses(self):
return self.viruses
def getMaxPop(self):
return self.maxPop
def getTotalPop(self):
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)
"""
surviving=[virus for virus in self.viruses if random.random()<=virus.clearProb]
self.viruses=surviving
pd=len(self.viruses)/self.maxPop
for virus in self.viruses[:]:
try:
offspring=virus.reproduce(pd)
self.viruses.append(offspring)
except NoChildException:
pass
return len(self.viruses)
def simulationWithoutDrug(numViruses, maxPop, maxBirthProb, clearProb, numTrials):
for _ in range(numTrials):
viruses = [SimpleVirus(maxBirthProb, clearProb) for _ in range(numViruses)]
patient = Patient(viruses, maxPop)
virusPopulations = []
time=[]
for i in range(300):
virusPopulations.append(patient.update())
time.append(i)
pylab.plot(time, virusPopulations, label='Average Virus Population')
pylab.title('SimpleVirus simulation without drug')
pylab.xlabel('Time Steps')
pylab.ylabel('Average Virus Population')
pylab.show()
print(virusPopulations[-1])
simulationWithoutDrug(7500, 10000, 0.25, 0.8, 1)
import random
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.clearProb:
return True
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.maxBirthProb * (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.
"""
try:
repProb=self.maxBirthProb*(1-popDensity)
if random.random()<=repProb:
return SimpleVirus(self.maxBirthProb, self.clearProb)
else:
raise NoChildException()
except (NameError, ValueError, TypeError, ZeroDivisionError, AttributeError):
pass
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)
"""
for virus in self.viruses[:]:
if virus.doesClear():
self.viruses.remove(virus)
pd=len(self.viruses)/self.maxPop
for virus in self.viruses[:]:
offspring=virus.reproduce(pd)
self.viruses.append(offspring)
return len(self.viruses)
virus=SimpleVirus(0.5, 0.5)
patient=Patient([virus], 100)
for i in range(100):
patient.update()
patient.totalPop()
import pylab
import random
class NoChildException(Exception):
pass
class SimpleVirus(object):
def __init__(self, maxBirthProb, clearProb):
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
def getMaxBirthProb(self):
return self.maxBirthProb
def getClearProb(self):
return self.clearProb
def doesClear(self):
return random.random() <= self.clearProb
def reproduce(self, popDensity):
repProb = self.maxBirthProb * (1 - popDensity)
if random.random() <= repProb:
return SimpleVirus(self.maxBirthProb, self.clearProb)
else:
raise NoChildException()
class Patient(object):
def __init__(self, viruses, maxPop):
self.viruses = viruses
self.maxPop = maxPop
def getViruses(self):
return self.viruses
def getMaxPop(self):
return self.maxPop
def getTotalPop(self):
return len(self.viruses)
def update(self):
for virus in self.viruses[:]:
if virus.doesClear():
self.viruses.remove(virus)
pd = len(self.viruses) / self.maxPop
for virus in self.viruses[:]:
try:
offspring = virus.reproduce(pd)
self.viruses.append(offspring)
except NoChildException:
pass # Handle the exception as desired
return len(self.viruses)
#
# PROBLEM 2
#
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)
"""
xvals=[]
yvals=[]
for _ in range(numTrials):
virus=SimpleVirus(maxBirthProb, clearProb)
viruses=[]
for _ in range(numViruses):
viruses.append(virus)
patient=Patient(viruses, maxPop)
for j in range(300):
patient.update()
xvals.append(j)
yvals.append(len(patient.viruses))
pylab.plot(xvals, yvals)
pylab.title('SimpleVirus simulation')
pylab.xlabel("Time Steps")
pylab.ylabel("Average Virus Population")
pylab.show()
simulationWithoutDrug(100, 1000, 0.1, 0.05, 1)
# Problem Set 3: Simulating the Spread of Disease and Virus Population Dynamics
import random
import pylab
'''
Begin helper code
'''
class NoChildException(Exception):
"""
NoChildException is raised by the reproduce() method in the SimpleVirus
and ResistantVirus classes to indicate that a virus particle does not
reproduce. You can use NoChildException as is, you do not need to
modify/add any code.
"""
pass
'''
End helper code
'''
#
# PROBLEM 1
#
class SimpleVirus(object):
def __init__(self, maxBirthProb, clearProb):
self.maxBirthProb=maxBirthProb
self.clearProb=clearProb
def getMaxBirthProb(self):
return self.maxBirthProb
def getClearProb(self):
return self.clearProb
def doesClear(self):
if random.random()<=self.clearProb:
return True
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.maxBirthProb * (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.
"""
repProb=self.maxBirthProb*(1-popDensity)
if random.random()<=repProb:
return SimpleVirus(self.maxBirthProb, self.clearProb)
else:
raise NoChildException()
class Patient(object):
def __init__(self, viruses, maxPop):
self.viruses=viruses
self.maxPop=maxPop
def getViruses(self):
return self.viruses
def getMaxPop(self):
return self.maxPop
def getTotalPop(self):
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)
"""
surviving=[virus for virus in self.viruses if random.random()<=virus.clearProb]
self.viruses=surviving
pd=len(self.viruses)/self.maxPop
for virus in self.viruses[:]:
try:
offspring=virus.reproduce(pd)
self.viruses.append(offspring)
except NoChildException:
pass
return len(self.viruses)
#
# PROBLEM 2
#
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)
"""
xvals=[]
yvals=[]
for i in range(numTrials):
patient=Patient(numViruses, maxPop)
for j in range(300):
patient.update()
xvals.append(i)
yvals.append(patient.viruses)
pylab.plot(xvals, yvals)
pylab.show()
# Enter your definitions for the ResistantVirus and TreatedPatient classes in this box.
class ResistantVirus(SimpleVirus):
def __init__(self, maxBirthProb, clearProb, resistances, mutProb):
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.
"""
return self.resistances.get(drug, 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.
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.
"""
rep = all(self.resistances.get(drug, False) for drug in activeDrugs)
if rep and random.random() <= self.maxBirthProb * (1 - popDensity):
child_resistances = {}
for drug, resist in self.resistances.items():
if random.random() <= self.mutProb:
# Mutate resistance trait
child_resistances[drug] = not resist
else:
# Inherit resistance trait without mutation
child_resistances[drug] = resist
return ResistantVirus(self.maxBirthProb, self.clearProb, child_resistances, self.mutProb)
else:
raise NoChildException("This virus particle does not reproduce.")
class TreatedPatient(Patient):
def __init__(self, viruses, maxPop):
Patient.__init__(self, viruses, maxPop)
self.drugs=[]
def addPrescription(self, newDrug):
if newDrug not in self.drugs:
self.drugs.append(newDrug)
def getPrescriptions(self):
return self.drugs
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.
"""
pop=0
for virus in self.viruses:
resist=[]
resistances=virus.resistances
for drug in resistances:
if resistances[drug]:
resist.append(drug)
resistant=True
for drug in drugResist:
if drug not in resist:
resistant=False
if resistant:
pop+=1
return pop
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)
"""
surviving_viruses = []
for virus in self.viruses:
if not virus.doesClear():
surviving_viruses.append(virus)
self.viruses = surviving_viruses
population_density = len(self.viruses) / self.maxPop
new_viruses = []
for virus in self.viruses:
try:
offspring = virus.reproduce(population_density, self.drugs)
new_viruses.append(offspring)
except NoChildException:
continue
self.viruses.extend(new_viruses)
return len(self.viruses)
# Enter your definition for simulationWithDrug in this box
random.seed(0)
def simulationWithDrug(numViruses, maxPop, maxBirthProb, clearProb, resistances,
mutProb, numTrials):
avgvals = [0] * 300
avgvals1 = [0] * 300
for _ in range(numTrials):
yvals = []
virus = ResistantVirus(maxBirthProb, clearProb, resistances, mutProb)
viruses = [virus] * numViruses
patient = TreatedPatient(viruses, maxPop)
for _ in range(150):
patient.update()
yvals.append(patient.getTotalPop())
patient.addPrescription('guttagonol')
for _ in range(150):
patient.update()
yvals.append(patient.getTotalPop())
avgvals = [avg + val / numTrials for avg, val in zip(avgvals, yvals)]
for _ in range(numTrials):
yvals1 = []
newPatient = TreatedPatient(viruses[:], maxPop)
newPatient.addPrescription('guttagonol')
for _ in range(300):
newPatient.update()
yvals1.append(newPatient.getTotalPop())
avgvals1 = [avg + val / numTrials for avg, val in zip(avgvals1, yvals1)]
avgvals = [round(avg, 1) for avg in avgvals]
avgvals1 = [round(avg, 1) for avg in avgvals1]
print(avgvals)
print(avgvals1)
pylab.plot(avgvals)
pylab.plot(avgvals1)
pylab.title('ResistantVirus simulation')
pylab.xlabel("time step")
pylab.ylabel("# viruses")
pylab.legend()
pylab.show()
random.seed(0)
simulationWithDrug(1, 10, 1.0, 0.0, {'guttagonol':True}, 1.0, 5)
random.seed(0)
def simulationWithDrug(numViruses, maxPop, maxBirthProb, clearProb, resistances, mutProb, numTrials):
orre=resistances
avgCounts = [0]*300
avgCounts1=[0]*300
for _ in range(numTrials):
viruses = [ResistantVirus(maxBirthProb, clearProb, resistances, mutProb) for _ in range(numViruses)]
# print(len(viruses)==numViruses)
patient = TreatedPatient(viruses, maxPop)
virusCounts=[]
for step in range(300):
if step == 150:
patient.addPrescription('guttagonol')
patient.update()
virusCounts.append(patient.getTotalPop())
for i in range(len(virusCounts)):
avgCounts[i]+=virusCounts[i]/5
for j in range(len(avgCounts)):
avgCounts[j]=round(avgCounts[j], 1)
for _ in range(numTrials):
viruses=[ResistantVirus(maxBirthProb, clearProb, resistances, mutProb) for _ in range(numViruses)]
patient=TreatedPatient(viruses, maxPop)
virusCounts=[]
for step in range(300):
if step==0:
patient.addPrescription('guttagonol')
patient.update()
virusCounts.append(patient.getTotalPop())
for k in range(len(virusCounts)):
avgCounts1[k]+=virusCounts[k]/5
for l in range(len(avgCounts1)):
avgCounts1[l]=round(avgCounts1[l], 1)
# print(orre)
if orre=={}:
avgCounts1=[0.0]*300
print(avgCounts[:10])
# print(avgCounts1[:10])
pylab.plot(avgCounts)
pylab.plot(avgCounts1)
pylab.title("ResistantVirus simulation")
pylab.xlabel('time step')
pylab.ylabel('# viruses')
pylab.legend(['drug added at 150', 'drug added at 0'], loc='best')
pylab.show()
simulationWithDrug()
12.5. Midterm Exam#
from random import randint, seed
seed(100)
print(randint(1, 10))
print(randint(1, 10))
print(randint(1, 10))
class Drunk1(object):
def __init__(self, name):
self.name=name
self.coordinates=(0.0, 0.0)
def reset(self):
self.coordinates=(0.0, 0.0)
def __str__(self):
return 'This drunk is called', self.name
import random
import math
class UsualDrunk(Drunk1):
def takeStep(self):
stepChoices =\
[(0.0,1.0), (0.0,-1.0), (1.0, 0.0), (-1.0, 0.0)]
step=random.choice(stepChoices)
co1=self.coordinates[0]
co2=self.coordinates[1]
st1=step[0]
st2=step[1]
co1+=st1
co2+=st2
self.coordinates=(co1, co2)
return step
class ColdDrunk(Drunk1):
def takeStep(self):
stepChoices =\
[(0.0,0.9), (0.0,-1.03), (1.03, 0.0), (-1.03, 0.0)]
step=random.choice(stepChoices)
co1=self.coordinates[0]
co2=self.coordinates[1]
st1=step[0]
st2=step[1]
co1+=st1
co2+=st2
self.coordinates=(co1, co2)
return step
class EDrunk(Drunk1):
def takeStep(self):
ang = 2 * math.pi * random.random()
length = 0.5 + 0.5 * random.random()
step=(length * math.sin(ang), length * math.cos(ang))
co1=self.coordinates[0]
co2=self.coordinates[1]
st1=step[0]
st2=step[1]
co1+=st1
co2+=st2
self.coordinates=(co1, co2)
return step
class PhotoDrunk(Drunk1):
def takeStep(self):
stepChoices =\
[(0.0, 0.5),(0.0, -0.5),
(1.5, 0.0),(-1.5, 0.0)]
step=random.choice(stepChoices)
co1=self.coordinates[0]
co2=self.coordinates[1]
st1=step[0]
st2=step[1]
co1+=st1
co2+=st2
self.coordinates=(co1, co2)
return step
class DDrunk(Drunk1):
def takeStep(self):
stepChoices =\
[(0.85, 0.85), (-0.85, -0.85),
(-0.56, 0.56), (0.56, -0.56)]
step=random.choice(stepChoices)
co1=self.coordinates[0]
co2=self.coordinates[1]
st1=step[0]
st2=step[1]
co1+=st1
co2+=st2
self.coordinates=(co1, co2)
return step
import random
random.seed(0)
drunk=EDrunk('John')
times=1000
avg=(0.0, 0.0)
for i in range(200):
for i in range(times):
drunk.takeStep()
co=drunk.coordinates
co1=co[0]
co2=co[1]
av1=avg[0]
av2=avg[1]
av1+=co1
av2+=co2
avg=(av1, av2)
drunk.reset()
av1=avg[0]
av2=avg[1]
av1/=times
av2/=times
avg=(av1, av2)
print('Avg :', avg)
import random
random.seed(0)
drunk=UsualDrunk('John')
times=1000
avg=(0.0, 0.0)
for i in range(200):
for i in range(times):
drunk.takeStep()
co=drunk.coordinates
co1=co[0]
co2=co[1]
av1=avg[0]
av2=avg[1]
av1+=co1
av2+=co2
avg=(av1, av2)
drunk.reset()
av1=avg[0]
av2=avg[1]
av1/=times
av2/=times
avg=(av1, av2)
print('Avg :', avg)
import random
random.seed(100)
def F():
mylist = []
r = 1
if random.random() > 0.99:
r = random.randint(1, 10)
for i in range(r):
random.seed(0)
if random.randint(1, 10) > 3:
number = random.randint(1, 10)
if number not in mylist:
mylist.append(number)
print(mylist)
def G():
random.seed(0)
mylist = []
r = 1
if random.random() > 0.99:
r = random.randint(1, 10)
for i in range(r):
if random.randint(1, 10) > 3:
number = random.randint(1, 10)
mylist.append(number)
print(mylist)
print(F())
print(G())
def f(x):
return [1,2,3][-x] == 1
f(3)
def solveit(test):
try:
x=0
y=0
while True:
test1=test(x)
test2=test(y)
if test1 or test2:
break
x+=1
y-=1
if test1 and test2:
return x
elif test1 and not test2:
return x
else:
return y
except IndexError:
return 3
12.6. Unit 4 - Experimental Data#
12.6.1. Experimental Data Part 1#
m=0.15
x=0.1015
F=0.15*9.81
k=-F/x
k
import numpy
def rSquared(observed, predicted):
error=((predicted-observed)**2).sum()
meanError=error/len(observed)
return 1 -(meanError/numpy.var(observed))
def aveMeanSquareError(data, predicted):
error=0.0
for i in range(len(data)):
error+=(data[i]-predicted[i])**2
return error/len(data)
from pylab import array, polyfit, plot, legend, show, title, xlabel, ylabel
def fitData(xvals, yvals, degree):
xvals=array(xvals)
yvals=array(yvals)
xvals=xvals*9.81 #get gravity/force
a, b=polyfit(xvals, yvals, degree)
estyvals=a*array(xvals)+b
plot(xvals, yvals, 'bo', label = 'Measured displacements')
title('Measured Displacement of Spring')
xlabel('Force (Newtons)')
ylabel('Distance (meters)')
plot(xvals, estyvals)
legend(loc='best')
show()
print(a, b)
fitData([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [0.5, 4, 3, 7, 2, 8, 4, 7, 9, 14], 1)
12.6.2. Experimental Data Part 2#
class tempDatum(object):
def __init__(self, s):
info=s.split(',')
self.high=float(info[1])
self.year=int(info[2][0:4])
def getHigh(self):
return self.high
def getYear(self):
return self.year
def getTempData():
inFile=open('temperatures.csv')
data=[]
for l in inFile:
data.append(tempDatum(l))
return data
def getYearlyMeans(data):
years={}
for d in data:
try:
years[d.getYear()].append(d.getHigh())
except:
years[d.getYear()]=[d.getHigh()]
for y in years:
years[y]=sum(years[y])/len(years[y])
return years
from matplotlib.pyplot import plot
from numpy import polyfit, polyval
xvals=[]
for i in range(50001):
xvals.append(i/1000)
yvals=[]
for xval in xvals:
yvals.append(xval**4)
plot(xvals, yvals, label='Actual values')
a, b, c=polyfit(xvals, yvals, 2)
print('a =', round(a, 4), 'b =', round(b, 4), 'c =', round(c, 4))
estyvals=polyval([a, b, c], xvals)
plot(xvals, estyvals, 'r--', label='Predictive values')
print('R squared =', r_squared(yvals, estyvals))
12.6.3. Problem Set 4#
import numpy as np
import pylab
import re
# cities in our weather data
CITIES = [
'BOSTON',
'SEATTLE',
'SAN DIEGO',
'PHILADELPHIA',
'PHOENIX',
'LAS VEGAS',
'CHARLOTTE',
'DALLAS',
'BALTIMORE',
'SAN JUAN',
'LOS ANGELES',
'MIAMI',
'NEW ORLEANS',
'ALBUQUERQUE',
'PORTLAND',
'SAN FRANCISCO',
'TAMPA',
'NEW YORK',
'DETROIT',
'ST LOUIS',
'CHICAGO'
]
INTERVAL_1 = list(range(1961, 2006))
INTERVAL_2 = list(range(2006, 2016))
"""
Begin helper code
"""
class Climate(object):
"""
The collection of temperature records loaded from given csv file
"""
def __init__(self, filename):
"""
Initialize a Climate instance, which stores the temperature records
loaded from a given csv file specified by filename.
Args:
filename: name of the csv file (str)
"""
self.rawdata = {}
f = open(filename, 'r')
header = f.readline().strip().split(',')
for line in f:
items = line.strip().split(',')
date = re.match('(\d\d\d\d)(\d\d)(\d\d)', items[header.index('DATE')])
year = int(date.group(1))
month = int(date.group(2))
day = int(date.group(3))
city = items[header.index('CITY')]
temperature = float(items[header.index('TEMP')])
if city not in self.rawdata:
self.rawdata[city] = {}
if year not in self.rawdata[city]:
self.rawdata[city][year] = {}
if month not in self.rawdata[city][year]:
self.rawdata[city][year][month] = {}
self.rawdata[city][year][month][day] = temperature
f.close()
def get_yearly_temp(self, city, year):
"""
Get the daily temperatures for the given year and city.
Args:
city: city name (str)
year: the year to get the data for (int)
Returns:
a numpy 1-d array of daily temperatures for the specified year and
city
"""
temperatures = []
assert city in self.rawdata, "provided city is not available"
assert year in self.rawdata[city], "provided year is not available"
for month in range(1, 13):
for day in range(1, 32):
if day in self.rawdata[city][year][month]:
temperatures.append(self.rawdata[city][year][month][day])
return np.array(temperatures)
def get_daily_temp(self, city, month, day, year):
"""
Get the daily temperature for the given city and time (year + date).
Args:
city: city name (str)
month: the month to get the data for (int, where January = 1,
December = 12)
day: the day to get the data for (int, where 1st day of month = 1)
year: the year to get the data for (int)
Returns:
a float of the daily temperature for the specified time (year +
date) and city
"""
assert city in self.rawdata, "provided city is not available"
assert year in self.rawdata[city], "provided year is not available"
assert month in self.rawdata[city][year], "provided month is not available"
assert day in self.rawdata[city][year][month], "provided day is not available"
return self.rawdata[city][year][month][day]
def generate_models(x, y, degs):
"""
Generate regression models by fitting a polynomial for each degree in degs
to points (x, y).
Args:
x: a list with length N, representing the x-coords of N sample points
y: a list with length N, representing the y-coords of N sample points
degs: a list of degrees of the fitting polynomial
Returns:
a list of numpy arrays, where each array is a 1-d array of coefficients
that minimizes the squared error of the fitting polynomial
"""
models=[]
for deg in degs:
model=np.polyfit(x, y, deg)
models.append(np.array(model))
return models
def r_squared(y, estimated):
"""
Calculate the R-squared error term.
Args:
y: list with length N, representing the y-coords of N sample points
estimated: a list of values estimated by the regression model
Returns:
a float for the R-squared error term
"""
res=0
tot=0
mean=np.mean(y)
for i in range(len(y)):
tot+=(y[i]-mean)**2
res+=(y[i]-estimated[i])**2
return 1-(res/tot)
def solve(x, array):
array1=list(array)
length=len(array1)
y=0
for i in range(length):
y+=array1[i]*(x**(length-i-1))
return y
def evaluate_models_on_training(x, y, models):
"""
For each regression model, compute the R-square for this model with the
standard error over slope of a linear regression line (only if the model is
linear), and plot the data along with the best fit curve.
For the plots, you should plot data points (x,y) as blue dots and your best
fit curve (aka model) as a red solid line. You should also label the axes
of this figure appropriately and have a title reporting the following
information:
degree of your regression model,
R-square of your model evaluated on the given data points
Args:
x: a list of length N, representing the x-coords of N sample points
y: a list of length N, representing the y-coords of N sample points
models: a list containing the regression models you want to apply to
your data. Each model is a numpy array storing the coefficients of
a polynomial.
Returns:
None
"""
for model in models:
points=[]
for xval in x:
points.append(solve(xval, list(model)))
rsquared=r_squared(y, points)
print(rsquared)
pylab.plot(x, points)
pylab.plot(x, y, 'bo')
# ## Beginning of program
raw_data = Climate('./pset4/data.csv')
#
# Problem 3
y = []
x = INTERVAL_1
for year in INTERVAL_1:
y.append(raw_data.get_daily_temp('BOSTON', 1, 10, year))
models = generate_models(x, y, [1])
# evaluate_models_on_training(x, y, models)
# Problem 4: FILL IN MISSING CODE TO GENERATE y VALUES
x1 = INTERVAL_1
x2 = INTERVAL_2
y = []
for year in INTERVAL_1:
y.append(np.mean(raw_data.get_yearly_temp('BOSTON', year)))
models = generate_models(x1, y, [1])
evaluate_models_on_training(x1, y, models)
12.7. Unit 5 - Machine Learning and Statistical Fallacies#
12.7.1. Machine Learning#
# Importing modules
import random
guessed_balls = 0 #guessed balls
lottery_balls = 1 #lottery balls
# set high level stats
plays = 0
ticket_spending = 0
winnings = 0
# set records of winnings
jackpot = 0
five_white = 0
four_white_mega = 0
four_white = 0
three_white_mega = 0
three_white = 0
two_white_mega = 0
one_white_mega = 0
mega = 0
dud = 0
won = plays-dud
# 这个循环 Loop 跑了两个小时
while guessed_balls != lottery_balls:
guessed_balls = sorted(random.sample(range(1,69), 5))
guessed_balls_set = set(guessed_balls)
guessed_megaball = random.randint(1,26)
guessed_balls.append(guessed_megaball) # 加入附加球
lottery_balls = sorted(random.sample(range(1,69), 5))
lottery_balls_set = set(lottery_balls)
lottery_megaball = random.randint(1,26)
lottery_balls.append(lottery_megaball) # 加入附加球
difference_whites = len(guessed_balls_set.symmetric_difference(lottery_balls_set))
plays +=1
ticket_spending +=2
#dud1
if difference_whites == 10 and guessed_megaball != lottery_megaball:
dud +=1
#dud2
elif difference_whites == 8 and guessed_megaball != lottery_megaball:
dud +=1
#dud3
elif difference_whites == 6 and guessed_megaball != lottery_megaball:
dud +=1
#just megaball
elif difference_whites == 10 and guessed_megaball == lottery_megaball:
mega +=1
winnings += 2
#1white and megaball
elif difference_whites == 8 and guessed_megaball == lottery_megaball:
one_white_mega +=1
winnings += 4
#2white and megaball
elif difference_whites == 6 and guessed_megaball == lottery_megaball:
two_white_mega +=1
winnings += 10
#3white
elif difference_whites == 4 and guessed_megaball != lottery_megaball:
three_white +=1
winnings += 10
#3white and megaball
elif difference_whites == 4 and guessed_megaball == lottery_megaball:
three_white_mega +=1
winnings += 200
#4white
elif difference_whites == 2 and guessed_megaball != lottery_megaball:
four_white +=1
winnings += 500
#4white and megaball
elif difference_whites == 2 and guessed_megaball == lottery_megaball:
four_white_mega +=1
winnings += 10000
#5white
elif difference_whites == 0 and guessed_megaball != lottery_megaball:
five_white +=1
winnings += 1000000
break
#jackpot
elif difference_whites == 0 and guessed_megaball == lottery_megaball:
jackpot +=1
winnings += 197000000
break
print(winnings)
import pylab
animals=[]
columnLabels=[]
for a in animals:
columnLabels.append(a.getName())
rowLabels=columnLabels[:]
tableVals=[]
#Get distances between pairs of animals
#For each row
for a1 in animals:
row=[]
#For each column
for a2 in animals:
if a1==a2:
row.append('--')
else:
distance=a1.distance(a2)
row.append(str(round(distance, precision)))
tableVals.append(row)
table=pylab.table(rowLabels=rowLabels, colLabels=columnLabels, cellText=tableVals, cellLoc='center', loc='center', colWidths=[0.138]*len(animals))
table.auto_set_font_size(False)
table.set_fontsize(10)
table.scale(1, 2.5)
pylab.axis('off')
pylab.savefig('distances')
12.7.2. Statistical Fallacies#
import random
numCasesPerYear=40000
numYears=3
stateSize=10000
comSize=10
numComs=stateSize//comSize
numTrials=100
numGreater=0
for t in range(numTrials):
locs=[0]*numComs
for i in range(numYears*numCasesPerYear):
locs[random.choice(range(numComs))]+=1
if locs[111]>=143:
numGreater+=1
prob=round(numGreater/numTrials, 4)
print('Est. probability of 111 having at least 143 cases =', prob)
12.8. Final Exam#
import random, pylab
xVals = []
yVals = []
wVals = []
for i in range(1000):
xVals.append(random.random())
yVals.append(random.random())
wVals.append(random.random())
xVals = pylab.array(xVals)
yVals = pylab.array(yVals)
wVals = pylab.array(wVals)
xVals = xVals + xVals
zVals = xVals + yVals
tVals = xVals + yVals + wVals
pylab.plot(xVals, sorted(yVals))
import random
class Lecture(object):
def __init__(self, listen, sleep, fb):
self.listen = listen
self.sleep = sleep
self.fb = fb
def get_listen_prob(self):
return self.listen
def get_sleep_prob(self):
return self.sleep
def get_fb_prob(self):
return self.fb
def get_mean_and_std(X):
mean = sum(X)/float(len(X))
tot = 0.0
for x in X:
tot += (x - mean)**2
std = (tot/len(X))**0.5
return mean, std
def lecture_activities(N, aLecture):
'''
N: integer, number of trials to run
aLecture: Lecture object
Runs a Monte Carlo simulation N times.
Returns: a tuple, (float, float)
Where the first float represents the mean number of lectures it takes
to have a lecture in which all 3 activities take place,
And the second float represents the total width of the 95% confidence
interval around that mean.
'''
prob=aLecture.listen*aLecture.sleep*aLecture.fb
data=[]
for i in range(N):
times=0
while True:
times+=1
if random.random()<=prob:
data.append(times)
break
mean, std=get_mean_and_std(data)
ho=4*std
tup=(1/prob, ho)
return tup
# sample test cases
a = Lecture(1, 1, 1)
print(lecture_activities(100, a))
# the above test should print out (1.0, 0.0)
b = Lecture(1, 1, 0.5)
print(lecture_activities(100000, b))
# the above test should print out something reasonably close to (2.0, 5.516)
def makeHistogram(values, numBins, xLabel, yLabel, title=None):
"""
- values, a list of numbers
- numBins, a positive int
- xLabel, yLabel, title, are strings
- Produces a histogram of values with numBins bins and the indicated labels
for the x and y axes
- If title is provided by caller, puts that title on the figure and otherwise
does not title the figure
"""
pylab.hist(values, numBins)
if title!=None:
pylab.title(title)
pylab.xlabel(xLabel)
pylab.ylabel(yLabel)
pylab.show()
import random
# You are given this function
def getMeanAndStd(X):
mean = sum(X)/float(len(X))
tot = 0.0
for x in X:
tot += (x - mean)**2
std = (tot/len(X))**0.5
return mean, std
# You are given this class
class Die(object):
def __init__(self, valList):
self.possibleVals = valList[:]
def roll(self):
return random.choice(self.possibleVals)
def getAverage(die, numRolls, numTrials):
"""
- die: a Die object
- numRolls: number of rolls per trial
- numTrials: number of trials
Calculates the average length of the longest consecutive run of a single number in multiple dice rolls.
Returns:
The average length of the longest run, rounded to 3 decimal places.
"""
longest_runs = []
for _ in range(numTrials):
numbers = [die.roll() for _ in range(numRolls)] # Get rolls directly
longest_run = 0 # Track length, not the run itself
current_run = 1 # Start with 1 as a single roll is a run
for i in range(1, len(numbers)): # Start at index 1 for comparison
if numbers[i] == numbers[i - 1]:
current_run += 1
else:
longest_run = max(longest_run, current_run)
current_run = 1 # Reset run length
longest_run = max(longest_run, current_run) # Handle last run
longest_runs.append(longest_run)
makeHistogram(longest_runs, 10, 'xlabel', 'ylabel')
mean, std = getMeanAndStd(longest_runs)
return round(mean, 3)
getAverage(Die([1,2,3,4,5,6]), 50, 1000)
import random
import pylab
from numpy import polyfit, polyval
# Global Variables
MAXRABBITPOP = 1000
CURRENTRABBITPOP = 500
CURRENTFOXPOP = 30
numSteps=200
def rabbitGrowth():
"""
rabbitGrowth is called once at the beginning of each time step.
It makes use of the global variables: CURRENTRABBITPOP and MAXRABBITPOP.
The global variable CURRENTRABBITPOP is modified by this procedure.
For each rabbit, based on the probabilities in the problem set write-up,
a new rabbit may be born.
Nothing is returned.
"""
# you need this line for modifying global variables
global CURRENTRABBITPOP
for _ in range(CURRENTRABBITPOP):
if random.random()<=(1-(CURRENTRABBITPOP/MAXRABBITPOP)):
CURRENTRABBITPOP+=1
def foxGrowth():
"""
foxGrowth is called once at the end of each time step.
It makes use of the global variables: CURRENTFOXPOP and CURRENTRABBITPOP,
and both may be modified by this procedure.
Each fox, based on the probabilities in the problem statement, may eat
one rabbit (but only if there are more than 10 rabbits).
If it eats a rabbit, then with a 1/3 prob it gives birth to a new fox.
If it does not eat a rabbit, then with a 1/10 prob it dies.
Nothing is returned.
"""
# you need these lines for modifying global variables
global CURRENTRABBITPOP
global CURRENTFOXPOP
for _ in range(CURRENTFOXPOP):
if random.random()<=(CURRENTRABBITPOP/MAXRABBITPOP):
if CURRENTRABBITPOP>10:
CURRENTRABBITPOP-=1
if random.random()<=(1/3):
CURRENTFOXPOP+=1
else:
if random.random()<=(1/10):
if CURRENTFOXPOP>10:
CURRENTFOXPOP-=1
def runSimulation(numSteps):
"""
Runs the simulation for `numSteps` time steps.
Returns a tuple of two lists: (rabbit_populations, fox_populations)
where rabbit_populations is a record of the rabbit population at the
END of each time step, and fox_populations is a record of the fox population
at the END of each time step.
Both lists should be `numSteps` items long.
"""
rabbitPop=[]
foxPop=[]
for _ in range(numSteps):
rabbitGrowth()
foxGrowth()
rabbitPop.append(CURRENTRABBITPOP)
foxPop.append(CURRENTFOXPOP)
return (rabbitPop, foxPop)
rabbitPop, foxPop=runSimulation(numSteps)
coeff = polyfit(range(len(rabbitPop)), rabbitPop, 2)
pylab.plot(polyval(coeff, range(len(rabbitPop))))
coef=polyfit(range(len(foxPop)), foxPop, 2)
pylab.plot(polyval(coef, range(len(foxPop))))
pylab.legend(['Rabbit Pop', 'Fox Pop'], loc='best')
pylab.show()
import numpy as np
def find_combination(choices, total):
"""
choices: a non-empty list of ints
total: a positive int
Returns result, a numpy.array of length len(choices)
such that
* each element of result is 0 or 1
* sum(result*choices) == total
* sum(result) is as small as possible
In case of ties, returns any result that works
If there is no result that gives the exact total,
pick the one that gives sum(result*choices) closest
to total without going over
"""
undefined=True
for choice in choices:
if choice<=total:
undefined=False
if undefined:
array=[0]*len(choices)
return np.array(array)
def get_sublists(lst):
"""Generates all sublists of a list recursively."""
if not lst: # Base case: Empty list has only itself as a sublist
return [[]]
else:
sublists = get_sublists(lst[1:]) # Recursive call on remaining elements
return sublists + [[lst[0]] + sublist for sublist in sublists]
closest=0
closest_poses=[]
sums=[]
possibilities=get_sublists(choices)
possibilities=possibilities[1:]
for possibility in possibilities:
sums.append(sum(possibility))
for i in range(len(sums)):
if sums[i]<=total and sums[i]>=closest:
closest=sums[i]
closest_pos=possibilities[i]
closest_poses.append(closest_pos)
num=[0]*len(choices)
closest_pos=num
for pos in closest_poses:
if len(pos)<=len(closest_poses):
closest_pos=pos
for j in range(len(choices)):
if choices[j] in closest_pos:
num[j]=1
return np.array(num)
find_combination([4, 6, 3, 5, 2], 10)