Windows Process Activation Service Could Not Be Started

I began today by starting a different post about publishing an app to localhost. Things soon went sideways. That post will have to wait while I peel back this onion:

When I open up IIS, I can see that the Default Site and the DefaultAppPool are stopped. No problem, let’s start them.

You knew it wasn’t going to be that simple, right? From an admin command prompt:

C:\WINDOWS\system32>net start was
The Windows Process Activation Service service is starting.
The Windows Process Activation Service service could not be started.

A system error has occurred.

System error 2 has occurred.

The system cannot find the file specified.

C:\WINDOWS\system32>net start was

The Windows Process Activation Service service is starting.

The Windows Process Activation Service service could not be started.

A system error has occurred.

System error 2 has occurred.

The system cannot find the file specified.

Now, WAS and IIS are both installed, so I follow the recommendation from iis.net to go ahead and reinstall them, but save yourself the restart. That didn’t do it.

Windows Features I’m not eager to Regedit based on a StackOverflow answer, but it’s simply removing a file and Scott Hanselman’s foray into PerfMon lead me to several errors, but they did not get me closer to a resolution. I’ve reinstalled WAS once, so I deleted NanoSetup.

C:\WINDOWS\system32>net start was
The Windows Process Activation Service service is starting.
The Windows Process Activation Service service was started successfully.

C:\WINDOWS\system32>net start was

The Windows Process Activation Service service is starting.

The Windows Process Activation Service service was started successfully.

Now the App Pool can be started, but we’re missing one last step to turn on the DefaultSite.

One more Admin Command Prompt later and we’re there.

C:\WINDOWS\system32>net start w3svc
The World Wide Web Publishing Service service is starting.
The World Wide Web Publishing Service service was started successfully.

C:\WINDOWS\system32>net start w3svc

The World Wide Web Publishing Service service is starting.

The World Wide Web Publishing Service service was started successfully.

Now, we can navigate to localhost in the browser.

Retrieving Items from a List in a Model

Let’s look at how you would get the List<Book> Books from the Author model. We’ll use Author’s Detail View and then extract that to a Partial View. A common error is to get a null pointer exception here:

@foreach (var item in Model.Books)

1	@foreach (var item in Model.Books)

Object reference not set to an instance of an object.

Description: An unhandled exception occurred during the execution of the current web request. Please review the stack trace for more information about the error and where it originated in the code.

Exception Details: System.NullReferenceException: Object reference not set to an instance of an object.

This is fixed by adding Virtual to the Books property. Now, when we go to the Author’s Detail page, the Books property is instantiated and the list can be displayed.

using System.Collections.Generic;

namespace BookStore.Models
{
    public class Author
    {
        public int AuthorId { get; set; }
        public string FirstName { get; set; }
        public string LastName { get; set; }
        public virtual List<Book> Books { get; set; }
    }
}

using System.Collections.Generic;

namespace BookStore.Models

{

public class Author

{

public int AuthorId { get; set; }

public string FirstName { get; set; }

public string LastName { get; set; }

public virtual List<Book> Books { get; set; }

}

namespace BookStore.Models
{
    public class Book
    {
        public int BookId { get; set; }
        public int AuthorId { get; set; }
        public Author Author { get; set; }
        public string Title { get; set; }
    }
}

namespace BookStore.Models

{

public class Book

{

public int BookId { get; set; }

public int AuthorId { get; set; }

public Author Author { get; set; }

public string Title { get; set; }

}

PM> Add-Migration Virtual
Scaffolding migration 'Virtual'.
PM> Update-Database
Applying explicit migration: 201811301319400_Virtual.
Running Seed method.

PM> Add-Migration Virtual

Scaffolding migration 'Virtual'.

PM> Update-Database

Applying explicit migration: 201811301319400_Virtual.

Running Seed method.

@model BookStore.Models.Author
@{
    ViewBag.Title = "Details";
}
<h2>Details</h2>
<div>
    <h4>Author</h4>
    <hr />
    <dl class="dl-horizontal">
        <dt>
            @Html.DisplayNameFor(model => model.FirstName)
        </dt>
        <dd>
            @Html.DisplayFor(model => model.FirstName)
        </dd>
        <dt>
            @Html.DisplayNameFor(model => model.LastName)
        </dt>
        <dd>
            @Html.DisplayFor(model => model.LastName)
        </dd>
        <dt>
            @foreach (var item in Model.Books)
            {
                @Html.DisplayFor(modelItem => item.Title)
            }
        </dt>
    </dl>
</div>

@model BookStore.Models.Author

ViewBag.Title = "Details";

}

<h2>Details</h2>

<div>

<h4>Author</h4>

<hr />

<dt>

@Html.DisplayNameFor(model => model.FirstName)

</dt>

<dd>

@Html.DisplayFor(model => model.FirstName)

</dd>

<dt>

@Html.DisplayNameFor(model => model.LastName)

</dt>

<dd>

@Html.DisplayFor(model => model.LastName)

</dd>

<dt>

@foreach (var item in Model.Books)

{

@Html.DisplayFor(modelItem => item.Title)

}

</dt>

</dl>

</div>

Using the Code First Approach with Entity Framework

Using Entity Framework, we can take a Code First Approach developing the data models and allow EF to translate the model classes into database tables. The properties of classes become the attributes of the tables. Defining an ID for a class will result in a primary key and referencing that ID from another class will create a foreign key relationship. The first step is to import EF to the application.

using System.Collections.Generic;

namespace BookStore.Models
{
    public class Author
    {
        public int AuthorId { get; set; }
        public string FirstName { get; set; }
        public string LastName { get; set; }
        public List<Book> Books { get; set; }
    }
}

using System.Collections.Generic;

namespace BookStore.Models

{

public class Author

{

public int AuthorId { get; set; }

public string FirstName { get; set; }

public string LastName { get; set; }

public List<Book> Books { get; set; }

}

namespace BookStore.Models
{
    public class Book
    {
        public int BookId { get; set; }
        public int AuthorId { get; set; }
        public Author Author { get; set; }
        public string Title { get; set; }
    }
}

namespace BookStore.Models

{

public class Book

{

public int BookId { get; set; }

public int AuthorId { get; set; }

public Author Author { get; set; }

public string Title { get; set; }

}

using System.Data.Entity;

namespace BookStore.Models
{
    public class BookStoreDb : DbContext
    {
        public BookStoreDb()
        {

        }
        public DbSet<Author> Authors { get; set; }
        public DbSet<Book> Books { get; set; }
    }
}

using System.Data.Entity;

namespace BookStore.Models

{

public class BookStoreDb : DbContext

{

public BookStoreDb()

{

}

public DbSet<Author> Authors { get; set; }

public DbSet<Book> Books { get; set; }

}

It will be useful to create Controllers and Views from the Models. With a DbContext Model in place: run build, then enable migrations from the PM Console, add an initial migration and run update-database. At this point, EF creates the database in LocalDb. See the initial create file in the Migrations folder for the operations run on update. Repeat with add-migration [name] whenever the models change and update-database.

PM> Enable-Migrations -ContextTypeName BookStoreDb
Checking if the context targets an existing database...
Code First Migrations enabled for project BookStore.
PM> Add-Migration InitialCreate
Scaffolding migration 'InitialCreate'.
The Designer Code for this migration file includes a snapshot of your current Code First model. 
This snapshot is used to calculate the changes to your model when you scaffold the next migration. 
If you make additional changes to your model that you want to include in this migration, 
then you can re-scaffold it by running 'Add-Migration InitialCreate' again.
PM> Update-Database
Specify the '-Verbose' flag to view the SQL statements being applied to the target database.
Applying explicit migrations: [201811280425589_InitialCreate].
Applying explicit migration: 201811280425589_InitialCreate.
Running Seed method.
PM>

PM> Enable-Migrations -ContextTypeName BookStoreDb

Checking if the context targets an existing database...

Code First Migrations enabled for project BookStore.

PM> Add-Migration InitialCreate

Scaffolding migration 'InitialCreate'.

The Designer Code for this migration file includes a snapshot of your current Code First model.

This snapshot is used to calculate the changes to your model when you scaffold the next migration.

If you make additional changes to your model that you want to include in this migration,

then you can re-scaffold it by running 'Add-Migration InitialCreate' again.

PM> Update-Database

Specify the '-Verbose' flag to view the SQL statements being applied to the target database.

Applying explicit migrations: [201811280425589_InitialCreate].

Applying explicit migration: 201811280425589_InitialCreate.

Running Seed method.

PM>

Further control can be gained by defining the connection string in the Web.config file and modify the Configuration.cs. For example, seeding the database with preliminary data. The connection string is also the place to attach the DB if you want to see it in the App_Data folder. It will appear after building and executing the app.

<?xml version="1.0" encoding="utf-8"?>
<configuration>
...
  <connectionStrings>
    <add name="BookStoreDb"
         connectionString="Data Source=(LocalDb)\MSSQLLocalDb; AttachDBFileName=|DataDirectory|\BookStore.mdf; Initial Catalog=BookStore; Integrated Security=SSPI"
         providerName="System.Data.SqlClient"/>
  </connectionStrings>
...
</configuration>

<?xml version="1.0" encoding="utf-8"?>

...

<add name="BookStoreDb"

connectionString="Data Source=(LocalDb)\MSSQLLocalDb; AttachDBFileName=|DataDirectory|\BookStore.mdf; Initial Catalog=BookStore; Integrated Security=SSPI"

providerName="System.Data.SqlClient"/>

</connectionStrings>

...

</configuration>

Storing Values in App.config

Storing values in the App.config file will allow others to update values without needing to access the source code and rebuilding the project. For example, when delivering a solution to another team which may need to change variables such as url, username or emails, then building the project with this in mind and retrieving values from App.config will allow them to easily make those changes.

Add references and using statements to System.Configuration
Add a key value pair to App.config
Get the value with ConfigurationManager.AppSettings
Edit the key value pair from the XML file MyApp.exe.config

On editing the config file, there will be no need to recompile the executable, but the values will be updated.

<?xml version="1.0" encoding="utf-8" ?>
<configuration>
  <appSettings>
    <add key="Username" value="Hoopy"/>
    <add key="Environment" value="Testing"/>
  </appSettings>
    <startup> 
        <supportedRuntime version="v4.0" sku=".NETFramework,Version=v4.6.1" />
    </startup>
</configuration>

<?xml version="1.0" encoding="utf-8" ?>

</appSettings>

</startup>

</configuration>

using System;
using System.Configuration;

namespace StoreValues
{
    class Program
    {
        static void Main(string[] args)
        {
            string user = ConfigurationManager.AppSettings["Username"];
            string env = ConfigurationManager.AppSettings["Environment"];
            Console.WriteLine("Hello, " + user + ".");
            Console.WriteLine("Welcome to the " + env + " Environment.");
            Console.ReadLine();
        }
    }
}

using System;

using System.Configuration;

namespace StoreValues

{

class Program

{

static void Main(string[] args)

{

string user = ConfigurationManager.AppSettings["Username"];

string env = ConfigurationManager.AppSettings["Environment"];

Console.WriteLine("Hello, " + user + ".");

Console.WriteLine("Welcome to the " + env + " Environment.");

Console.ReadLine();

}

Hello, HoopyTest.
Welcome to the Testing Environment.
>

Hello, HoopyTest.

Welcome to the Testing Environment.

<?xml version="1.0" encoding="utf-8" ?>
<configuration>
  <appSettings>
    <add key="Username" value="HoopyProd"/>
    <add key="Environment" value="Production"/>
  </appSettings>
    <startup> 
        <supportedRuntime version="v4.0" sku=".NETFramework,Version=v4.6.1" />
    </startup>
</configuration>

<?xml version="1.0" encoding="utf-8" ?>

</appSettings>

</startup>

</configuration>

Hello, HoopyProd.
Welcome to the Production Environment.
>

Hello, HoopyProd.

Welcome to the Production Environment.

Sending Email from .NET using Gmail SMTP

I automate lots of jobs to run daily. Mostly just getting Task Scheduler to run a bat file or python script and I can just check from my desktop if something failed and re-run it. However, if I am in .NET and something is going to run from another server that I won’t necessarily be able to access whenever I want, receiving an email that a job succeeded or failed would be nice. I am hoping for a more permanent solution, but in pinch, a gmail account will work.

using System;
using System.Net;
using System.Net.Mail;

namespace SmtpWithGmail
{
    class Program
    {
        static void Main(string[] args)
        {
            // This may be deprecated but it works for now
            // docs.microsoft.com/en-us/dotnet/api/system.net.mail.smtpclient
            MailMessage msg = new MailMessage();
            msg.Body = "Sleep tight. All went well.";
            msg.Subject = "Scheduled job succeeded";
            msg.From = new MailAddress("[mySender]@gmail.com");
            msg.To.Add(new MailAddress("[myReceiver]@gmail.com"));
            SmtpClient smtp = new SmtpClient();
            smtp.Host = "smtp.gmail.com";
            smtp.Port = 587;
            smtp.Credentials = new NetworkCredential("[mySender]@gmail.com"
                , "[myPass]");
            smtp.EnableSsl = true;
            smtp.Send(msg);
        }
    }
}

using System;

using System.Net;

using System.Net.Mail;

namespace SmtpWithGmail

{

class Program

{

static void Main(string[] args)

{

// This may be deprecated but it works for now

// docs.microsoft.com/en-us/dotnet/api/system.net.mail.smtpclient

MailMessage msg = new MailMessage();

msg.Body = "Sleep tight. All went well.";

msg.Subject = "Scheduled job succeeded";

msg.From = new MailAddress("[mySender]@gmail.com");

msg.To.Add(new MailAddress("[myReceiver]@gmail.com"));

SmtpClient smtp = new SmtpClient();

smtp.Host = "smtp.gmail.com";

smtp.Port = 587;

smtp.Credentials = new NetworkCredential("[mySender]@gmail.com"

, "[myPass]");

smtp.EnableSsl = true;

smtp.Send(msg);

}

Connect to Oracle Data Source from .NET – A Console App

In order to connect to an Oracle database from .NET, use the Oracle.ManagedDataAccess package, otherwise called ODP.NET. It can be downloaded with NuGet by searching “oracle.”

You can get started easily with defining a connection string. Later, you can put this into App.Config with different Data Sources for Testing and Production environments. Don’t forget the using statement.

using System;
using Oracle.ManagedDataAccess.Client;

namespace OracleConnect
{
    class Program
    {
        static void Main(string[] args)
        {
            string conString = "User Id=HR;Password=HR;" +
                "Data Source=localhost:1521/xe";

            OracleConnection con = new OracleConnection(conString);
            con.Open();
            OracleCommand cmd = con.CreateCommand();
            cmd.CommandText = "Select first_name, last_name from employees";
            OracleDataReader reader = cmd.ExecuteReader();
            while (reader.Read())
            {
                Console.WriteLine(reader.GetString(0) + "\t" 
				+ reader.GetString(1));
            }
            con.Dispose();
            Console.ReadLine();
        }
    }
}

using System;

using Oracle.ManagedDataAccess.Client;

namespace OracleConnect

{

class Program

{

static void Main(string[] args)

{

string conString = "User Id=HR;Password=HR;" +

"Data Source=localhost:1521/xe";

OracleConnection con = new OracleConnection(conString);

con.Open();

OracleCommand cmd = con.CreateCommand();

cmd.CommandText = "Select first_name, last_name from employees";

OracleDataReader reader = cmd.ExecuteReader();

while (reader.Read())

{

Console.WriteLine(reader.GetString(0) + "\t"

+ reader.GetString(1));

}

con.Dispose();

Console.ReadLine();

}

...
Nandita Sarchand
Ismael  Sciarra
John    Seo
Sarath  Sewall
Lindsey Smith
William Smith
Stephen Stiles
Martha  Sullivan
Patrick Sully
>_

...

Nandita Sarchand

Ismael Sciarra

John Seo

Sarath Sewall

Lindsey Smith

William Smith

Stephen Stiles

Martha Sullivan

Patrick Sully

Python: Problem Set 11

Problem 11 is a weighting digraph problem, finding your way across campus given building numbers and sidewalks connecting them. In order to tackle this type of problem, I did use some outside tutorials on graph problems and depth-first search from TopCoder and the Python.org websites.

Part one involves finding the shortest path using a brute force method. One key to solving this is to keep from going into an infinite loop with cycles through the same nodes over and over. My chosen method is to look for the next node in the path I’ve already taken. This way [path] does double duty as my result and keeps track of nodes visited.

I banged my head on the dynamic programming portion of the problem long enough that I have to wrap this up without taking into account path weights. First, all hops are measured equally. Think number of airport connections made on a trip, but ignores miles flown. Second, I simply could not keep track of the current shortest distance as a break/continue condition without passing it as an argument when I recursed. There may not be anything wrong with this per se, but it is outside the problem specifications, so I feel I should disclose.

Thanks to MIT for providing all the course materials. It was fun.

# **PART OF** campus-map.txt
# start, end, dist, distOutdoors
32 36 70 0
36 32 70 0
32 57 30 0
57 32 30 0
32 76 80 50

# **PART OF** campus-map.txt

# start, end, dist, distOutdoors

32 36 70 0

36 32 70 0

32 57 30 0

57 32 30 0

32 76 80 50

# ps11.py
import time
class Node(object):
    def __init__(self, name):
        self.name = str(name)
    def getName(self):
        return self.name
    def __str__(self):
        return self.name
    def __eq__(self, other):
        return self.name == other.name
    def __repr__(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 str(self.src) + '->' + str(self.dest)

class Wedge(Edge):
    def __init__(self, *a, **k):
        if 'tot' in k:
            self.tot = k['tot']
        if 'out' in k:
            self.out = k['out']
        Edge.__init__(self, *a)
    def getTotal(self):
        return self.tot
    def getOutside(self):
        return self.out

class Digraph(object):
    def __init__(self):
        self.nodes = set([])
        self.edges = {}
    def addNode(self, node):
        unique = True
        for i in self.nodes:
            if i == node:
                unique = False
        if unique:
            self.nodes.add(node)
            self.edges[node] = []
    def addEdge(self, edge):
        self.edges[edge.src].append(edge)
    def childrenOf(self, node):
        return self.edges[node]
    def getNode(self, node):
        for j in self.nodes:
            if j == node:
                return j
    def __str__(self):
        res = ''
        for k in self.edges:
            for d in self.edges[k]:
                res = res + str(k) + '->' + str(d) + '\n'
        return res[:-1]        
        
def load_map(mapFilename):

    campus = Digraph()
    print "Loading map from file..."
    f = open(mapFilename, 'r')
    for line in f:
        data = line.split()
        # if nodes on campus, copy them
        if campus.getNode(Node(data[0])):
            a = campus.getNode(Node(data[0]))
        else:
            a = Node(data[0])
        if campus.getNode(Node(data[1])):
            b = campus.getNode(Node(data[1]))
        else:
            b = Node(data[1])
        c = Wedge(a, b, tot=data[2], out=data[3])
        campus.addNode(a)
        campus.addNode(b)
        # since nodes are not duplicates (new instances), edge.src won't keyError edges{}
        campus.addEdge(c)
        
    return campus

def find_path(graph, start, end, path=[]):
    start = Node(start)
    end = Node(end)        
    # because I want to copy the nodes, not make new instances
    for i in graph.nodes:
        if i == start:
            start = i
        if i == end:
            end = i
    path = path + [start]            
    if start == end:
        return path
    shortest = None
    for wedge in graph.childrenOf(start):
        if wedge.getDestination() not in path:
            newpath = find_path(graph, wedge.getDestination(), end, path=path)
            if newpath:
                if not shortest or len(newpath) < len(shortest):
                    shortest = newpath
    return shortest

graph = load_map("mit_map.txt")
begin = time.time()
print "find_path(graph, '32', '56')"
print find_path(graph, '32', '56')
finish = time.time()
print "seconds: ", finish - begin

# had to add a shortest arg to make work, else scope and recursion headache
def dp_find_path(graph, start, end, path=[], shortest=[0]*999):
    start = Node(start)
    end = Node(end)        
    # because I want to copy the nodes, not make new instances
    for i in graph.nodes:
        if i == start:
            start = i
        if i == end:
            end = i
    path = path + [start]
    if start == end:
        return path
    for wedge in graph.childrenOf(start):
        if len(path) > len(shortest):
            continue
        if wedge.getDestination() not in path:
            newpath = dp_find_path(graph, wedge.getDestination(), end, path, shortest)
            if newpath:
                if not shortest or len(newpath) < len(shortest):
                    shortest = newpath
    return shortest
#graph = load_map("campus-map.txt")
begin = time.time()
print "dp_find_path(graph, '32', '56')"
print dp_find_path(graph, '32', '56')
finish = time.time()
print "seconds: ", finish - begin

begin = time.time()
print "dp_find_path(graph, '2', '9')"
print dp_find_path(graph, '2', '9')
finish = time.time()
print "seconds: ", finish - begin

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

# ps11.py

import time

class Node(object):

def __init__(self, name):

self.name = str(name)

def getName(self):

return self.name

def __str__(self):

return self.name

def __eq__(self, other):

return self.name == other.name

def __repr__(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 str(self.src) + '->' + str(self.dest)

class Wedge(Edge):

def __init__(self, *a, **k):

if 'tot' in k:

self.tot = k['tot']

if 'out' in k:

self.out = k['out']

Edge.__init__(self, *a)

def getTotal(self):

return self.tot

def getOutside(self):

return self.out

class Digraph(object):

def __init__(self):

self.nodes = set([])

self.edges = {}

def addNode(self, node):

unique = True

for i in self.nodes:

if i == node:

unique = False

if unique:

self.nodes.add(node)

self.edges[node] = []

def addEdge(self, edge):

self.edges[edge.src].append(edge)

def childrenOf(self, node):

return self.edges[node]

def getNode(self, node):

for j in self.nodes:

if j == node:

return j

def __str__(self):

res = ''

for k in self.edges:

for d in self.edges[k]:

res = res + str(k) + '->' + str(d) + '\n'

return res[:-1]

def load_map(mapFilename):

campus = Digraph()

print "Loading map from file..."

f = open(mapFilename, 'r')

for line in f:

data = line.split()

# if nodes on campus, copy them

if campus.getNode(Node(data[0])):

a = campus.getNode(Node(data[0]))

else:

a = Node(data[0])

if campus.getNode(Node(data[1])):

b = campus.getNode(Node(data[1]))

else:

b = Node(data[1])

c = Wedge(a, b, tot=data[2], out=data[3])

campus.addNode(a)

campus.addNode(b)

# since nodes are not duplicates (new instances), edge.src won't keyError edges{}

campus.addEdge(c)

return campus

def find_path(graph, start, end, path=[]):

start = Node(start)

end = Node(end)

# because I want to copy the nodes, not make new instances

for i in graph.nodes:

if i == start:

start = i

if i == end:

end = i

path = path + [start]

if start == end:

return path

shortest = None

for wedge in graph.childrenOf(start):

if wedge.getDestination() not in path:

newpath = find_path(graph, wedge.getDestination(), end, path=path)

if newpath:

if not shortest or len(newpath) < len(shortest):

shortest = newpath

return shortest

graph = load_map("mit_map.txt")

begin = time.time()

print "find_path(graph, '32', '56')"

print find_path(graph, '32', '56')

finish = time.time()

print "seconds: ", finish - begin

# had to add a shortest arg to make work, else scope and recursion headache

def dp_find_path(graph, start, end, path=[], shortest=[0]*999):

start = Node(start)

end = Node(end)

# because I want to copy the nodes, not make new instances

for i in graph.nodes:

if i == start:

start = i

if i == end:

end = i

path = path + [start]

if start == end:

return path

for wedge in graph.childrenOf(start):

if len(path) > len(shortest):

continue

if wedge.getDestination() not in path:

newpath = dp_find_path(graph, wedge.getDestination(), end, path, shortest)

if newpath:

if not shortest or len(newpath) < len(shortest):

shortest = newpath

return shortest

#graph = load_map("campus-map.txt")

begin = time.time()

print "dp_find_path(graph, '32', '56')"

print dp_find_path(graph, '32', '56')

finish = time.time()

print "seconds: ", finish - begin

begin = time.time()

print "dp_find_path(graph, '2', '9')"

print dp_find_path(graph, '2', '9')

finish = time.time()

print "seconds: ", finish - begin

Python: Problem Set 10

Problem Set 10 uses K-Means clustering to group counties into similar groups based on 14 demographic features such as income, housing, education, farm acres, age, population, population change and a several others.

The first task is to measure the sum of the squared error (the distance between the points in a cluster and the center of the cluster).

Next, we look at the effects of using increasing numbers of initial clusters (K). If the number of clusters were equal to the the number of points in the dataset, then SSE would be zero and we would be finished. This, however, does not tell us much (if anything) about the relationship between points. So choosing a good K number is important.

In the county dataset, we see that for increasing numbers of clusters (k-vals), the SSE is nearly always on a decreasing trend.

The final part of the problem is to weight some features so as to cluster tighter on the poverty rate without using that feature in the measurement. Weighting all features except poverty ends in about a 44,000 SSE at 125 clusters. But, weighting for income and high school graduation rate was a noticeable decrease, about 31,000.

# Problem Set 10
# Name: Casey

#Code shared across examples
import pylab, random, string, copy, math

# creates a point class to hold points, including names, attributes and methods
class Point(object):
    def __init__(self, name, originalAttrs, normalizedAttrs = None):
        """normalizedAttrs and originalAttrs are both arrays"""
        self.name = name
        self.unNormalized = originalAttrs
        self.attrs = normalizedAttrs
    def dimensionality(self):
        return len(self.attrs)
    def getAttrs(self):
        return self.attrs
    def getOriginalAttrs(self):
        return self.unNormalized
    def distance(self, other):
        #Euclidean distance metric
        difference = self.attrs - other.attrs
        return sum(difference * difference) ** 0.5
    def getName(self):
        return self.name
    def toStr(self):
        return self.name + str(self.attrs)
    def __str__(self):
        return self.name

# county is a type of point ( a subclass ) just to handle weighing features
class County(Point):
    weights = pylab.array([1.0] * 14)
    # Override Point.distance to use County.weights to decide the
    # significance of each dimension
    # when do custom weights get passed in?
    def distance(self, other):
        difference = self.getAttrs() - other.getAttrs()
        # sq root the sum of weights times differences in attributes gives dist
        # is this the Minkowski metric?
        return sum(County.weights * difference * difference) ** 0.5

class Cluster(object):
    def __init__(self, points, pointType):
        self.points = points
        self.pointType = pointType
        # using func below
        self.centroid = self.computeCentroid()
    def getCentroid(self):
        return self.centroid
    def computeCentroid(self):
        # dimensionality from superclass
        dim = self.points[0].dimensionality()
        # init totVals to 0.0
        totVals = pylab.array([0.0]*dim)
        for p in self.points:
            totVals += p.getAttrs()
        meanPoint = self.pointType('mean',
                                   totVals/float(len(self.points)),
                                   totVals/float(len(self.points)))
        return meanPoint
    def update(self, points):
        oldCentroid = self.centroid
        self.points = points
        if len(points) > 0:
            self.centroid = self.computeCentroid()
            return oldCentroid.distance(self.centroid)
        else:
            return 0.0
    def getPoints(self):
        return self.points
    def contains(self, name):
        for p in self.points:
            if p.getName() == name:
                return True
        return False
    def toStr(self):
        result = ''
        for p in self.points:
            result = result + p.toStr() + ', '
        return result[:-2]
    def __str__(self):
        result = ''
        for p in self.points:
            result = result + str(p) + ', '
        return result[:-2]

def kmeans(points, k, cutoff, pointType, minIters = 3, maxIters = 100,
toPrint = False):
    """ Returns (Cluster list, max dist of any point to its cluster) """
    #Uses random initial centroids
    initialCentroids = random.sample(points,k)
    clusters = []
    for p in initialCentroids:
        clusters.append(Cluster([p], pointType))
    numIters = 0
    biggestChange = cutoff
    while (biggestChange >= cutoff and numIters < maxIters) or numIters < minIters:
        print "Starting iteration " + str(numIters)
        newClusters = []
        for c in clusters:
            newClusters.append([])
        for p in points:
            smallestDistance = p.distance(clusters[0].getCentroid())
            index = 0
            for i in range(len(clusters)):
                distance = p.distance(clusters[i].getCentroid())
                if distance < smallestDistance:
                    smallestDistance = distance
                    index = i
            newClusters[index].append(p)
        biggestChange = 0.0
        for i in range(len(clusters)):
            change = clusters[i].update(newClusters[i])
            biggestChange = max(biggestChange, change)
        numIters += 1
        if toPrint:
            print 'Iteration count =', numIters
    maxDist = 0.0
    for c in clusters:
        for p in c.getPoints():
            if p.distance(c.getCentroid()) > maxDist:
                maxDist = p.distance(c.getCentroid())
    print 'Total Number of iterations =', numIters, 'Max Diameter =', maxDist
    print biggestChange
    return clusters, maxDist

#US Counties example
def readCountyData(fName, numEntries = 14):
    dataFile = open(fName, 'r')
    dataList = []
    nameList = []
    maxVals = pylab.array([0.0]*numEntries)
    #Build unnormalized feature vector
    for line in dataFile:
        if len(line) == 0 or line[0] == '#':
            continue
        dataLine = string.split(line)
        name = dataLine[0] + dataLine[1]
        features = []
        #Build vector with numEntries features
        for f in dataLine[2:]:
            try:
                f = float(f)
                features.append(f)
                if f > maxVals[len(features)-1]:
                    maxVals[len(features)-1] = f
            except ValueError:
                name = name + f
        if len(features) != numEntries:
            continue
        dataList.append(features)
        nameList.append(name)
    return nameList, dataList, maxVals

def buildCountyPoints(fName):
    """
    Given an input filename, reads County values from the file and returns
    them all in a list.
    """
    nameList, featureList, maxVals = readCountyData(fName)
    points = []
    for i in range(len(nameList)):
        originalAttrs = pylab.array(featureList[i])
        normalizedAttrs = originalAttrs/pylab.array(maxVals)
        points.append(County(nameList[i], originalAttrs, normalizedAttrs))
    return points

def randomPartition(l, p):
    """
    Splits the input list into two partitions, where each element of l is
    in the first partition with probability p and the second one with
    probability (1.0 - p).

    l: The list to split
    p: The probability that an element of l will be in the first partition

    Returns: a tuple of lists, containing the elements of the first and
    second partitions.
    """
    l1 = []
    l2 = []
    for x in l:
        if random.random() < p:
            l1.append(x)
        else:
            l2.append(x)
    return (l1,l2)

def getAveIncome(cluster):
    """
    Given a Cluster object, finds the average income field over the members
    of that cluster.

    cluster: the Cluster object to check

    Returns: a float representing the computed average income value
    """
    tot = 0.0
    numElems = 0
    for c in cluster.getPoints():
        tot += c.getOriginalAttrs()[1]
    return float(tot) / len(cluster.getPoints())

def getAvePoverty(cluster):

    tot = 0.0
    numElems = 0
    for c in cluster.getPoints():
        tot += c.getOriginalAttrs()[2]
    return float(tot) / len(cluster.getPoints())

def test(points, k = 200, cutoff = 0.1):
    """
    A sample function to show you how to do a simple kmeans run and graph
    the results.
    """
    # but this does not use any weights, no?
    incomes = []
    print ''
    clusters, maxSmallest = kmeans(points, k, cutoff, County)

    for i in range(len(clusters)):
        if len(clusters[i].points) == 0: continue
        incomes.append(getAveIncome(clusters[i]))

    pylab.hist(incomes)
    pylab.xlabel('Ave. Income')
    pylab.ylabel('Number of Clusters')
    pylab.show()

# create a list of points for all counties
points = buildCountyPoints('counties.txt')
#random.seed(123)
# create a sample using only 10% of the points from counties
testPoints = random.sample(points, len(points)/10)
# test(testPoints)
# print testPoints


def graphRemovedErr(points, kvals = [25, 50, 75, 100, 125, 150], cutoff = 0.1):
    """
    Should produce graphs of the error in training and holdout point sets, and
    the ratio of the error of the points, after clustering for the given
    values of k.
    For details see Problem 1.
    """

    # Your Code Here
    # 1.1 partition training and holdout sets
    training, holdout = randomPartition(points, 0.8)
    # 1.2 sum the distance of all points to their cluster's centroid
    t_list = []
    h_list = []
    ratios = []
    for k in kvals:
        # find sum of squared dist of each point in training set
        clusters, maxSmallest = kmeans(training, k, cutoff, County)
        dists = []
        for c in clusters:
            for p in c.getPoints():
                dists.append(p.distance(c.getCentroid())**2)
        t_tot = sum(dists)
        t_list.append(t_tot)

    # find sum of squared dist of each point in training set
        clusters, maxSmallest = kmeans(holdout, k, cutoff, County)
        dists = []
        for c in clusters:
            for p in c.getPoints():
                dists.append(p.distance(c.getCentroid())**2)
        h_tot = sum(dists)
        h_list.append(h_tot)
        ratios.append(h_tot/t_tot)
    # 2.1 graph errors against k
    pylab.plot(h_list)
    pylab.plot(t_list)
    pylab.xlabel('k vals')
    pylab.xticks(range(len(kvals)), kvals)
    pylab.ylabel('Distance Error')
    pylab.show()
    # 2.2 graph ratio of error for holdout set and training set in separate
    plot
    pylab.plot(ratios)
    pylab.xlabel('k vals')
    pylab.xticks(range(len(kvals)), kvals)
    pylab.ylabel('holdout over training ratio')
    pylab.show()
    return t_list, h_list # a tuple of two lists
# testPoints is too small for this @10% of Counties.txt w/ holdout of 20%
#print graphRemovedErr(points)

def kmeansAndYou(points, k = 50, cutoff = 0.1):
    for i in range(3):
        clusters, maxSmallest = kmeans(points, k, cutoff, County)
        for c in clusters:
            for p in c.getPoints():
                if p.getName() == "AZCoconino":
                    print c, "\n\n"
#print kmeansAndYou(points)

def graphPredictionErr(points, dimension, kvals = [25, 50, 75, 100, 125], cutoff = 0.1):
    """
    Given input points and a dimension to predict, should cluster on the
    appropriate values of k and graph the error in the resulting predictions,
    as described in Problem 3.
    """

	# TODO
    # weight features
    for p in points:
        p.attrs = dimension * p.getAttrs()
    hold_dist = []
    for k in kvals:
       # find sum of squared dist of each cluster
        clusters, maxSmallest = kmeans(points, k, cutoff, County)
        poverty = []
        for c in clusters:
            avePoverty = getAvePoverty(c)
            for p in c.getPoints():
                ptPoverty = p.getOriginalAttrs()[2]
                poverty.append((ptPoverty - avePoverty)**2)
        tot = sum(poverty)
        hold_dist.append(tot)
    return hold_dist, kvals
    
def plotAll(points, otherWeights, zeroWeighted, fullWeighted):
    featureBlind, kvals = graphPredictionErr(points, zeroWeighted)
    featurePredict, kvals = graphPredictionErr(points, otherWeights)
    
    pylab.plot(featurePredict, label="HS grad & Income")
    pylab.plot(featureBlind, label="All except Poverty %")
    
    pylab.xlabel("k vals")
    pylab.xticks(range(len(kvals)), kvals)    
    pylab.ylabel("Sum Squared Error from Poverty")
    pylab.legend()
    pylab.show()

'''
(('HomeVal', '1'), ('Income', '0'), ('Poverty', '1'),
                    ('Population', '1'), ('Pop Change', '1'),
                    ('Prcnt 65+', '1'), ('Below 18', '1'),
                    ('Prcnt Female', '1'), ('Prcent HS Grad', '1'),
                    ('Prcnt College', '0'), ('Unemployed', '0'),
                    ('Prcent Below 18', '1'), ('Life Expect', '1'),
                    ('Farm Acres', '1'))
'''
otherWeights = [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
zeroWeighted = [1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
fullWeighted = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
plotAll(points, otherWeights, zeroWeighted, fullWeighted)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

# Problem Set 10

# Name: Casey

#Code shared across examples

import pylab, random, string, copy, math

# creates a point class to hold points, including names, attributes and methods

class Point(object):

def __init__(self, name, originalAttrs, normalizedAttrs = None):

"""normalizedAttrs and originalAttrs are both arrays"""

self.name = name

self.unNormalized = originalAttrs

self.attrs = normalizedAttrs

def dimensionality(self):

return len(self.attrs)

def getAttrs(self):

return self.attrs

def getOriginalAttrs(self):

return self.unNormalized

def distance(self, other):

#Euclidean distance metric

difference = self.attrs - other.attrs

return sum(difference * difference) ** 0.5

def getName(self):

return self.name

def toStr(self):

return self.name + str(self.attrs)

def __str__(self):

return self.name

# county is a type of point ( a subclass ) just to handle weighing features

class County(Point):

weights = pylab.array([1.0] * 14)

# Override Point.distance to use County.weights to decide the

# significance of each dimension

# when do custom weights get passed in?

def distance(self, other):

difference = self.getAttrs() - other.getAttrs()

# sq root the sum of weights times differences in attributes gives dist

# is this the Minkowski metric?

return sum(County.weights * difference * difference) ** 0.5

class Cluster(object):

def __init__(self, points, pointType):

self.points = points

self.pointType = pointType

# using func below

self.centroid = self.computeCentroid()

def getCentroid(self):

return self.centroid

def computeCentroid(self):

# dimensionality from superclass

dim = self.points[0].dimensionality()

# init totVals to 0.0

totVals = pylab.array([0.0]*dim)

for p in self.points:

totVals += p.getAttrs()

meanPoint = self.pointType('mean',

totVals/float(len(self.points)),

totVals/float(len(self.points)))

return meanPoint

def update(self, points):

oldCentroid = self.centroid

self.points = points

if len(points) > 0:

self.centroid = self.computeCentroid()

return oldCentroid.distance(self.centroid)

else:

return 0.0

def getPoints(self):

return self.points

def contains(self, name):

for p in self.points:

if p.getName() == name:

return True

return False

def toStr(self):

result = ''

for p in self.points:

result = result + p.toStr() + ', '

return result[:-2]

def __str__(self):

result = ''

for p in self.points:

result = result + str(p) + ', '

return result[:-2]

def kmeans(points, k, cutoff, pointType, minIters = 3, maxIters = 100,

toPrint = False):

""" Returns (Cluster list, max dist of any point to its cluster) """

#Uses random initial centroids

initialCentroids = random.sample(points,k)

clusters = []

for p in initialCentroids:

clusters.append(Cluster([p], pointType))

numIters = 0

biggestChange = cutoff

while (biggestChange >= cutoff and numIters < maxIters) or numIters < minIters:

print "Starting iteration " + str(numIters)

newClusters = []

for c in clusters:

newClusters.append([])

for p in points:

smallestDistance = p.distance(clusters[0].getCentroid())

index = 0

for i in range(len(clusters)):

distance = p.distance(clusters[i].getCentroid())

if distance < smallestDistance:

smallestDistance = distance

index = i

newClusters[index].append(p)

biggestChange = 0.0

for i in range(len(clusters)):

change = clusters[i].update(newClusters[i])

biggestChange = max(biggestChange, change)

numIters += 1

if toPrint:

print 'Iteration count =', numIters

maxDist = 0.0

for c in clusters:

for p in c.getPoints():

if p.distance(c.getCentroid()) > maxDist:

maxDist = p.distance(c.getCentroid())

print 'Total Number of iterations =', numIters, 'Max Diameter =', maxDist

print biggestChange

return clusters, maxDist

#US Counties example

def readCountyData(fName, numEntries = 14):

dataFile = open(fName, 'r')

dataList = []

nameList = []

maxVals = pylab.array([0.0]*numEntries)

#Build unnormalized feature vector

for line in dataFile:

if len(line) == 0 or line[0] == '#':

continue

dataLine = string.split(line)

name = dataLine[0] + dataLine[1]

features = []

#Build vector with numEntries features

for f in dataLine[2:]:

try:

f = float(f)

features.append(f)

if f > maxVals[len(features)-1]:

maxVals[len(features)-1] = f

except ValueError:

name = name + f

if len(features) != numEntries:

continue

dataList.append(features)

nameList.append(name)

return nameList, dataList, maxVals

def buildCountyPoints(fName):

"""

Given an input filename, reads County values from the file and returns

them all in a list.

"""

nameList, featureList, maxVals = readCountyData(fName)

points = []

for i in range(len(nameList)):

originalAttrs = pylab.array(featureList[i])

normalizedAttrs = originalAttrs/pylab.array(maxVals)

points.append(County(nameList[i], originalAttrs, normalizedAttrs))

return points

def randomPartition(l, p):

"""

Splits the input list into two partitions, where each element of l is

in the first partition with probability p and the second one with

probability (1.0 - p).

l: The list to split

p: The probability that an element of l will be in the first partition

Returns: a tuple of lists, containing the elements of the first and

second partitions.

"""

l1 = []

l2 = []

for x in l:

if random.random() < p:

l1.append(x)

else:

l2.append(x)

return (l1,l2)

def getAveIncome(cluster):

"""

Given a Cluster object, finds the average income field over the members

of that cluster.

cluster: the Cluster object to check

Returns: a float representing the computed average income value

"""

tot = 0.0

numElems = 0

for c in cluster.getPoints():

tot += c.getOriginalAttrs()[1]

return float(tot) / len(cluster.getPoints())

def getAvePoverty(cluster):

tot = 0.0

numElems = 0

for c in cluster.getPoints():

tot += c.getOriginalAttrs()[2]

return float(tot) / len(cluster.getPoints())

def test(points, k = 200, cutoff = 0.1):

"""

A sample function to show you how to do a simple kmeans run and graph

the results.

"""

# but this does not use any weights, no?

incomes = []

print ''

clusters, maxSmallest = kmeans(points, k, cutoff, County)

for i in range(len(clusters)):

if len(clusters[i].points) == 0: continue

incomes.append(getAveIncome(clusters[i]))

pylab.hist(incomes)

pylab.xlabel('Ave. Income')

pylab.ylabel('Number of Clusters')

pylab.show()

# create a list of points for all counties

points = buildCountyPoints('counties.txt')

#random.seed(123)

# create a sample using only 10% of the points from counties

testPoints = random.sample(points, len(points)/10)

# test(testPoints)

# print testPoints

def graphRemovedErr(points, kvals = [25, 50, 75, 100, 125, 150], cutoff = 0.1):

"""

Should produce graphs of the error in training and holdout point sets, and

the ratio of the error of the points, after clustering for the given

values of k.

For details see Problem 1.

"""

# Your Code Here

# 1.1 partition training and holdout sets

training, holdout = randomPartition(points, 0.8)

# 1.2 sum the distance of all points to their cluster's centroid

t_list = []

h_list = []

ratios = []

for k in kvals:

# find sum of squared dist of each point in training set

clusters, maxSmallest = kmeans(training, k, cutoff, County)

dists = []

for c in clusters:

for p in c.getPoints():

dists.append(p.distance(c.getCentroid())**2)

t_tot = sum(dists)

t_list.append(t_tot)

# find sum of squared dist of each point in training set

clusters, maxSmallest = kmeans(holdout, k, cutoff, County)

dists = []

for c in clusters:

for p in c.getPoints():

dists.append(p.distance(c.getCentroid())**2)

h_tot = sum(dists)

h_list.append(h_tot)

ratios.append(h_tot/t_tot)

# 2.1 graph errors against k

pylab.plot(h_list)

pylab.plot(t_list)

pylab.xlabel('k vals')

pylab.xticks(range(len(kvals)), kvals)

pylab.ylabel('Distance Error')

pylab.show()

# 2.2 graph ratio of error for holdout set and training set in separate

plot

pylab.plot(ratios)

pylab.xlabel('k vals')

pylab.xticks(range(len(kvals)), kvals)

pylab.ylabel('holdout over training ratio')

pylab.show()

return t_list, h_list # a tuple of two lists

# testPoints is too small for this @10% of Counties.txt w/ holdout of 20%

#print graphRemovedErr(points)

def kmeansAndYou(points, k = 50, cutoff = 0.1):

for i in range(3):

clusters, maxSmallest = kmeans(points, k, cutoff, County)

for c in clusters:

for p in c.getPoints():

if p.getName() == "AZCoconino":

print c, "\n\n"

#print kmeansAndYou(points)

def graphPredictionErr(points, dimension, kvals = [25, 50, 75, 100, 125], cutoff = 0.1):

"""

Given input points and a dimension to predict, should cluster on the

appropriate values of k and graph the error in the resulting predictions,

as described in Problem 3.

"""

# TODO

# weight features

for p in points:

p.attrs = dimension * p.getAttrs()

hold_dist = []

for k in kvals:

# find sum of squared dist of each cluster

clusters, maxSmallest = kmeans(points, k, cutoff, County)

poverty = []

for c in clusters:

avePoverty = getAvePoverty(c)

for p in c.getPoints():

ptPoverty = p.getOriginalAttrs()[2]

poverty.append((ptPoverty - avePoverty)**2)

tot = sum(poverty)

hold_dist.append(tot)

return hold_dist, kvals

def plotAll(points, otherWeights, zeroWeighted, fullWeighted):

featureBlind, kvals = graphPredictionErr(points, zeroWeighted)

featurePredict, kvals = graphPredictionErr(points, otherWeights)

pylab.plot(featurePredict, label="HS grad & Income")

pylab.plot(featureBlind, label="All except Poverty %")

pylab.xlabel("k vals")

pylab.xticks(range(len(kvals)), kvals)

pylab.ylabel("Sum Squared Error from Poverty")

pylab.legend()

pylab.show()

'''

(('HomeVal', '1'), ('Income', '0'), ('Poverty', '1'),

('Population', '1'), ('Pop Change', '1'),

('Prcnt 65+', '1'), ('Below 18', '1'),

('Prcnt Female', '1'), ('Prcent HS Grad', '1'),

('Prcnt College', '0'), ('Unemployed', '0'),

('Prcent Below 18', '1'), ('Life Expect', '1'),

('Farm Acres', '1'))

'''

otherWeights = [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]

zeroWeighted = [1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

fullWeighted = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

plotAll(points, otherWeights, zeroWeighted, fullWeighted)

Python: Problem Set 9

The topic this week is optimization and clustering. This problem is less intensive than last week, but it does demonstrate the trade offs between finding an optimal solution (time intensive) and an adequate solution to an optimization problem.

The point here is to optimize a student’s schedule given two factors: workload and value as defined in the given dictionary. E.g. {‘6.00′: (10, 1)}, the class 6.00 is worth ten times its estimated workload. Naturally, there must be a constraint to our optimization, in this case the student’s desired maximum workload e.g. 12 might return {‘6.00′: (10, 1), ‘6.01’: (5, 4), ‘6.02’: (5, 6)}. With a workload of 11, no additional courses could be added to the schedule without exceeding the maximum.

The initial optimizations went quickly, creating three possible sorting functions: work, value or val/work ratio. The brute force best possible function was more tricky. I first used itertools.permutation method to calculate all possible solutions, in order to then loop through and pick one. This died on the short_subject_list.txt with just 20 courses. My function here works well enough on eight items, but I will have to check the staff solution to see how we are expected to carry out such a calculation on 20 items.

EDIT: OCW has a solution for PS9, but it is for a completely different assignment. But, I realized I need a combination of all possible schedules, not a permutation of each. This makes a sample size of twenty courses manageable.

# 6 .00 Problem Set 9
#
# Intelligent Course Advisor
#
#
import itertools

SUBJECT_FILENAME = "subjects.txt"
SHORT_SUBJECT_FILENAME = "shortened_subjects.txt"
VALUE, WORK = 0, 1

#
# Problem 1: Building A Subject Dictionary
#
def loadSubjects(filename):
    """
    Returns a dictionary mapping subject name to (value, work), where the name
    is a string and the value and work are integers. The subject information is
    read from the file named by the string filename. Each line of the file
    contains a string of the form "name,value,work".

    returns: dictionary mapping subject name to (value, work)
    """

    # TODO: Instead of printing each line, modify the above to parse the name,
    # value, and work of each subject and create a dictionary mapping the name
    # to the (value, work).
    inputFile = open(filename)
    catalog = {}
    for line in inputFile:
        line = line.split(",")
        line[2] = line[2].strip()
        catalog[line[0]] = (int(line[1]), int(line[2]))
        
    return catalog

#loadSubjects(SHORT_SUBJECT_FILENAME)

def printSubjects(subjects):
    """
    Prints a string containing name, value, and work of each subject in
    the dictionary of subjects and total value and work of all subjects
    """
    totalVal, totalWork = 0,0
    if len(subjects) == 0:
        return 'Empty SubjectList'
    res = 'Course\tValue\tWork\n======\t====\t=====\n'
    subNames = subjects.keys()
    subNames.sort()
    for s in subNames:
        val = subjects[s][VALUE]
        work = subjects[s][WORK]
        res = res + s + '\t' + str(val) + '\t' + str(work) + '\n'
        totalVal += val
        totalWork += work
    res = res + '\nTotal Value:\t' + str(totalVal) +'\n'
    res = res + 'Total Work:\t' + str(totalWork) + '\n'
    print res

#
# Problem 2: Subject Selection By Greedy Optimization
#

def cmpValue(subInfo1, subInfo2):
    """
    Returns True if value in (value, work) tuple subInfo1 is GREATER than
    value in (value, work) tuple in subInfo2
    """
    # TODO...
    return subInfo1[VALUE] > subInfo2[VALUE]

def cmpWork(subInfo1, subInfo2):
    """
    Returns True if work in (value, work) tuple subInfo1 is LESS than than work
    in (value, work) tuple in subInfo2
    """
    # TODO...
    return subInfo1[WORK] < subInfo2[WORK]

def cmpRatio(subInfo1, subInfo2):
    """
    Returns True if value/work in (value, work) tuple subInfo1 is 
    GREATER than value/work in (value, work) tuple in subInfo2
    """
    # TODO...
    return (float(subInfo1[VALUE])/subInfo1[WORK]) > (float(subInfo2[VALUE])/subInfo2[WORK])

def greedyAdvisor(subjects, maxWork, comparator):
    """
    Returns a dictionary mapping subject name to (value, work) which includes
    subjects selected by the algorithm, such that the total work of subjects in
    the dictionary is not greater than maxWork.  The subjects are chosen using
    a greedy algorithm.  The subjects dictionary should not be mutated.

    subjects: dictionary mapping subject name to (value, work)
    maxWork: int >= 0
    comparator: function taking two tuples and returning a bool
    returns: dictionary mapping subject name to (value, work)
    """
    # TODO...
    courses = {}
    work = 0
    choices = []    
    for i in subjects.items():
        choices.append(i)
    #print choices
    # bubble sort!
    # assume sorted
    isSorted = False
    while not isSorted:
        isSorted = True
        for i in range(0, len(choices)-1):
            temp = 0
            if comparator(choices[i][1],choices[i+1][1]):
                isSorted = False
                temp = choices[i]
                choices[i] = choices[i+1]
                choices[i+1] = temp
                break
            else:
                continue
    choices = list(reversed(choices))
    for k,v in choices:
        if work + v[WORK] <= maxWork:
            courses[k] = v
            work += v[WORK]
    return courses
    
    
#subjects = loadSubjects(SHORT_SUBJECT_FILENAME)
#print greedyAdvisor(subjects, 5, cmpRatio)


#
# Problem 3: Subject Selection By Brute Force
#
def bruteForceAdvisor(subjects, maxWork):
    """
    Returns a dictionary mapping subject name to (value, work), which
    represents the globally optimal selection of subjects using a brute force
    algorithm.

    subjects: dictionary mapping subject name to (value, work)
    maxWork: int >= 0
    returns: dictionary mapping subject name to (value, work)
    """
    # TODO...
    # build a list of all possible solutions, pick best
    subjects = subjects.items()
    options = []
    best = []
    for i in xrange(1, len(subjects)+1):
        options.append(itertools.combinations(subjects, i))
    for schedules in options:
        for courses in schedules:
            work = 0
            value = 0
            bestValue = 0
            for i in courses:
                work += i[1][WORK]
            if work <= maxWork:
                for i in courses:
                    value += i[1][VALUE]
                for i in best:
                    bestValue += i[1][VALUE]
                if value > bestValue:
                    best = courses
    return best

subjects = loadSubjects(SHORT_SUBJECT_FILENAME)
print bruteForceAdvisor(subjects, 7)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

# 6 .00 Problem Set 9

# Intelligent Course Advisor

import itertools

SUBJECT_FILENAME = "subjects.txt"

SHORT_SUBJECT_FILENAME = "shortened_subjects.txt"

VALUE, WORK = 0, 1

# Problem 1: Building A Subject Dictionary

def loadSubjects(filename):

"""

Returns a dictionary mapping subject name to (value, work), where the name

is a string and the value and work are integers. The subject information is

read from the file named by the string filename. Each line of the file

contains a string of the form "name,value,work".

returns: dictionary mapping subject name to (value, work)

"""

# TODO: Instead of printing each line, modify the above to parse the name,

# value, and work of each subject and create a dictionary mapping the name

# to the (value, work).

inputFile = open(filename)

catalog = {}

for line in inputFile:

line = line.split(",")

line[2] = line[2].strip()

catalog[line[0]] = (int(line[1]), int(line[2]))

return catalog

#loadSubjects(SHORT_SUBJECT_FILENAME)

def printSubjects(subjects):

"""

Prints a string containing name, value, and work of each subject in

the dictionary of subjects and total value and work of all subjects

"""

totalVal, totalWork = 0,0

if len(subjects) == 0:

return 'Empty SubjectList'

res = 'Course\tValue\tWork\n======\t====\t=====\n'

subNames = subjects.keys()

subNames.sort()

for s in subNames:

val = subjects[s][VALUE]

work = subjects[s][WORK]

res = res + s + '\t' + str(val) + '\t' + str(work) + '\n'

totalVal += val

totalWork += work

res = res + '\nTotal Value:\t' + str(totalVal) +'\n'

res = res + 'Total Work:\t' + str(totalWork) + '\n'

print res

# Problem 2: Subject Selection By Greedy Optimization

def cmpValue(subInfo1, subInfo2):

"""

Returns True if value in (value, work) tuple subInfo1 is GREATER than

value in (value, work) tuple in subInfo2

"""

# TODO...

return subInfo1[VALUE] > subInfo2[VALUE]

def cmpWork(subInfo1, subInfo2):

"""

Returns True if work in (value, work) tuple subInfo1 is LESS than than work

in (value, work) tuple in subInfo2

"""

# TODO...

return subInfo1[WORK] < subInfo2[WORK]

def cmpRatio(subInfo1, subInfo2):

"""

Returns True if value/work in (value, work) tuple subInfo1 is

GREATER than value/work in (value, work) tuple in subInfo2

"""

# TODO...

return (float(subInfo1[VALUE])/subInfo1[WORK]) > (float(subInfo2[VALUE])/subInfo2[WORK])

def greedyAdvisor(subjects, maxWork, comparator):

"""

Returns a dictionary mapping subject name to (value, work) which includes

subjects selected by the algorithm, such that the total work of subjects in

the dictionary is not greater than maxWork. The subjects are chosen using

a greedy algorithm. The subjects dictionary should not be mutated.

subjects: dictionary mapping subject name to (value, work)

maxWork: int >= 0

comparator: function taking two tuples and returning a bool

returns: dictionary mapping subject name to (value, work)

"""

# TODO...

courses = {}

work = 0

choices = []

for i in subjects.items():

choices.append(i)

#print choices

# bubble sort!

# assume sorted

isSorted = False

while not isSorted:

isSorted = True

for i in range(0, len(choices)-1):

temp = 0

if comparator(choices[i][1],choices[i+1][1]):

isSorted = False

temp = choices[i]

choices[i] = choices[i+1]

choices[i+1] = temp

break

else:

continue

choices = list(reversed(choices))

for k,v in choices:

if work + v[WORK] <= maxWork:

courses[k] = v

work += v[WORK]

return courses

#subjects = loadSubjects(SHORT_SUBJECT_FILENAME)

#print greedyAdvisor(subjects, 5, cmpRatio)

# Problem 3: Subject Selection By Brute Force

def bruteForceAdvisor(subjects, maxWork):

"""

Returns a dictionary mapping subject name to (value, work), which

represents the globally optimal selection of subjects using a brute force

algorithm.

subjects: dictionary mapping subject name to (value, work)

maxWork: int >= 0

returns: dictionary mapping subject name to (value, work)

"""

# TODO...

# build a list of all possible solutions, pick best

subjects = subjects.items()

options = []

best = []

for i in xrange(1, len(subjects)+1):

options.append(itertools.combinations(subjects, i))

for schedules in options:

for courses in schedules:

work = 0

value = 0

bestValue = 0

for i in courses:

work += i[1][WORK]

if work <= maxWork:

for i in courses:

value += i[1][VALUE]

for i in best:

bestValue += i[1][VALUE]

if value > bestValue:

best = courses

return best

subjects = loadSubjects(SHORT_SUBJECT_FILENAME)

print bruteForceAdvisor(subjects, 7)

# shortened_subjects.txt
6.00,10,1
6.01,5,4
6.02,5,6
6.03,2,9
6.04,1,2
6.05,1,18
6.06,5,19
6.07,2,10
6.08,1,10
6.09,3,7
6.10,8,18
6.11,6,8
6.12,6,3
6.13,9,16
6.14,10,8
6.15,10,6
6.16,6,9
6.17,9,3
6.18,10,4
6.19,8,19

# shortened_subjects.txt

6.00,10,1

6.01,5,4

6.02,5,6

6.03,2,9

6.04,1,2

6.05,1,18

6.06,5,19

6.07,2,10

6.08,1,10

6.09,3,7

6.10,8,18

6.11,6,8

6.12,6,3

6.13,9,16

6.14,10,8

6.15,10,6

6.16,6,9

6.17,9,3

6.18,10,4

6.19,8,19

Python: Problem Set 8

Problem set 8 built on the virus and patient model created in ps7. The main things to look at are comparisons of patients with different medication regimens. In part 2, we look at the results of patients prescribed ‘guttagonal’ with a control group that receives no medication.

In part five, we see that there are very similar looking graphs for two groups of patients: group A receives one prescription followed by a delay a second prescription, while the second group, B, receives two anti-virals simultaneously. What the graph does not show is that all 100 patients in group A are virus free by the end of the trial, but group B has 20-25 patients with single digit virus counts remaining.

# 6.00 Problem Set 8
#
# Name: Casey
# Collaborators:
# Time: ~8hr

import numpy
import random
import pylab
from ps7 import *

#
# PROBLEM 1
#
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, 'grimpex',False}, means that this virus
        particle is resistant to neither guttagonol nor grimpex.

        mutProb: Mutation probability for this virus particle (a float). This is
        the probability of the offspring acquiring or losing resistance to a drug.        

        """

        # TODO
        assert 0 < maxBirthProb < 1 and 0 < clearProb < 1 and 0 < mutProb < 1
        assert isinstance(resistances, dict)
        
        self.maxBirthProb = maxBirthProb
        self.clearProb = clearProb
        self.resistances = resistances
        self.mutProb = mutProb

    def isResistantTo(self, drug):

        """
        Get the state of this virus particle's resistance to a drug. This method
        is called by getResistPop() in Patient 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.
        """

        # TODO
        # resistances: A dictionary of drug names (strings)
        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 Patient class.

        If the virus particle is not resistant to any drug in activeDrugs,
        then it does not reproduce. Otherwise, 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 ResistantVirus (which has the same
        maxBirthProb and clearProb values as its 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
        grimpex, 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
        grimpex and a 90% chance that the offspring will not be resistant to
        grimpex.

        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.         
        """
        # TODO
        # if not resistant to a drug in activeDrugs: don't reproduce
        # else: reproduce with birthProb and clearProb

        resistant = False
        newResistances = self.resistances
        # if no active drugs, go to random.random()
        if len(activeDrugs) > 0:
            for i in activeDrugs:
                if self.isResistantTo(i):
                    resistant = True
            if resistant:
                for k,v in self.resistances.items():
                    if random.random() < self.mutProb:
                        # reverse resistance in newResistances
                        newResistances[k] = not v
            else:
                return NoChildException()        

        if random.random() < self.maxBirthProb * (1 - popDensity):
            return ResistantVirus(self.maxBirthProb, self.clearProb, newResistances, self.mutProb)
        else:
            return NoChildException()
        
#a = ResistantVirus(.99, .25, {'guttagonal': True, 'grimpex': True}, .99)
#b = a.reproduce(.1, ['grimpex'])
#print b.resistances

            

class Patient(SimplePatient):

    """
    Representation of a patient. The patient is able to take drugs and his/her
    virus population can acquire resistance to the drugs he/she takes.
    """

    def __init__(self, viruses, maxPop):
        """
        Initialization function, saves the viruses and maxPop parameters as
        attributes. Also initializes the list of drugs being administered
        (which should initially include no drugs).               

        viruses: the list representing the virus population (a list of
        SimpleVirus instances)
        
        maxPop: the  maximum virus population for this patient (an integer)
        """
        # TODO
        assert isinstance(viruses, list)
        assert isinstance(maxPop, int)
        self.viruses = viruses
        self.maxPop = maxPop
    

    def addPrescription(self, newDrug):

        """
        Administer a drug to this patient. After a prescription is added, the 
        drug acts on the virus population for all subsequent time steps. If the
        newDrug is already prescribed to this patient, the method has no effect.

        newDrug: The name of the drug to administer to the patient (a string).

        postcondition: list of drugs being administered to a patient is updated
        """
        # TODO
        # should not allow one drug being added to the list multiple times
        # if patient has no drug list: make one 
        # if newDrug on list: don't add
        # else: add drug
        assert isinstance(newDrug, str)
        
        try:
            if newDrug in self.drugs:
                return
            self.drugs.append(newDrug) 
        except AttributeError:
            self.drugs = []      
            self.drugs.append(newDrug)


    def getPrescriptions(self):

        """
        Returns the drugs that are being administered to this patient.
        returns: The list of drug names (strings) being administered to this
        patient.
        """

        # TODO
        try:
            return self.drugs
        except AttributeError:
            return []
        

    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', 'grimpex'])

        returns: the population of viruses (an integer) with resistances to all
        drugs in the drugResist list.
        """
        # TODO
        assert isinstance(drugResist, list)
        
        count = 0
        for i in drugResist:
            count += self.viruses.count(i)
        return count


    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().
        - Determine whether each virus particle should reproduce and add
          offspring virus particles to the list of viruses in this patient. 
          The listof 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)
        """
        # TODO
        survive_list = []
        for i in range(0, len(self.viruses)):
            if not self.viruses[i].doesClear():
                survive_list.append(self.viruses[i])
# debug        print "survive list: ", len(survive_list)
        self.viruses = survive_list
        # update the pop density
        now_density  =  self.getTotalPop() / float(self.maxPop)
        # for each virus particle, reproduce or not
        new_viruses = []
        for i in range(0, len(self.viruses)):
            virus = self.viruses[i]
            a = virus.reproduce(now_density, self.getPrescriptions())
            if isinstance(a, ResistantVirus):
                new_viruses.append(a)
# debug        print "new viruses: ", len(new_viruses)        
        self.viruses = self.viruses + new_viruses
        # return total virus pop as an int
        return self.getTotalPop()


#
# PROBLEM 2
#

def simulationWithDrug():

    """

    Runs simulations and plots graphs for problem 4.
    Instantiates a patient, runs a simulation for 150 timesteps, adds
    guttagonol, and runs the simulation for an additional 150 timesteps.
    total virus population vs. time and guttagonol-resistant virus population
    vs. time are plotted
    """
    # TODO
    a = ResistantVirus(0.1, 0.05, {'guttagonal': False}, 0.005)
    viruses = [a]*100
    trialPatient = Patient(viruses, 1000)
    controlPatient = Patient(viruses, 1000)
    guttagonal = []
    control = []
    for i in range(0, 150):
        guttagonal.append(trialPatient.update())
        control.append(controlPatient.update())
    trialPatient.addPrescription('guttagonal')
    for i in range(0, 150):
        guttagonal.append(trialPatient.update()) 
        control.append(controlPatient.update())
    pylab.plot(guttagonal, 'r-', label='guttagonal prescribed')
    pylab.plot(control, 'b-', label='without medication')
    pylab.xlabel("Time")
    pylab.ylabel("Virus count")
    pylab.legend(loc='best')
    pylab.title("Virus when drugs applied at time 150")    
    pylab.show()

#simulationWithDrug()
'''
    virus *was* dying off before we get to timestep 150
    was not resistant to empty activeDrugs
'''
#
# PROBLEM 3
#        

def simulationDelayedTreatment(stepsDelayed, numTrials):

    """
    Runs simulations and make histograms for problem 5.
    Runs multiple simulations to show the relationship between delayed treatment
    and patient outcome.
    Histograms of final total virus populations are displayed for delays of 300,
    150, 75, 0 timesteps (followed by an additional 150 timesteps of
    simulation).    
    """

    # TODO
    # run 30 trials for each scenario
    # displays histograms of final virus counts as four colored bars
    virusPops = []
    a = ResistantVirus(0.1, 0.05, {'guttagonal': True}, 0.005)
    viruses = [a]*100
    for k in range(0, numTrials):
        # reset patients for each trial
        myPatients = []
        for i in range(0, len(stepsDelayed)):
            a = Patient(viruses, 1000)
            myPatients.append(a)
        # run through all steps with delayed prescriptions
        # capture final viral counts
        statsStep = []
        counter = 0
        for i in range(0, max(stepsDelayed)+150):
            if counter in stepsDelayed:
                myPatients[stepsDelayed.index(counter)].addPrescription('guttagonal')
            for j in myPatients:
                j.update()
            counter += 1
        # one final result for each patient each trial
        for z in myPatients:
            statsStep.append(z.getTotalPop())
        print statsStep
        virusPops.append(statsStep)
    # for hist to work, list should convert to numpy.array()
    virusPops = numpy.array(virusPops)
    '''
        README
        100 trials shows 389 patients free of virus and 11 with 1 virus left
        needs legend
        and axis labels
    '''
    stepsDelayed = [str(i) for i in stepsDelayed] # convert ints in list to strs
    pylab.hist(virusPops, 8, histtype='bar', color=['blue', 'red', 'green', 'purple'],
    label=stepsDelayed)
    pylab.xlabel('Virus cells remaining')
    pylab.ylabel('Patients')
    pylab.title('Outcomes for delayed treatment')
    pylab.legend() # labels won't appear without it
    pylab.show()
        
#simulationDelayedTreatment([0,75,150,300], 100)

#
# PROBLEM 4
#

def simulationTwoDrugsDelayedTreatment(stepsDelayed, numTrials):

    """
    Runs simulations and make histograms for problem 6.
    Runs multiple simulations to show the relationship between administration
    of multiple drugs and patient outcome.
   
    Histograms of final total virus populations are displayed for lag times of
    150, 75, 0 timesteps between adding drugs (followed by an additional 150
    timesteps of simulation).
    """

    # TODO
    # run 30 trials for each scenario
    # plot 
    virusPops = []
    a = ResistantVirus(0.1, 0.05, {'guttagonal': False, 'grimpex': False}, 0.005)
    viruses = [a]*100
    for k in range(0, numTrials):
        # reset patients for each trial
        myPatients = []
        for i in range(0, len(stepsDelayed)):
            a = Patient(viruses, 1000)
            myPatients.append(a)

        # run an initial 150 steps prior to guttagonal
        for i in range(0, 150):
            for j in myPatients:
                j.update()
        for j in myPatients:
            j.addPrescription('guttagonal')                        

        # apply grimpex at delays defined in stepsDelayed
        counter = 0
        for i in range(0, max(stepsDelayed)+150):
            if counter in stepsDelayed:
                myPatients[stepsDelayed.index(counter)].addPrescription('grimpex')
            for j in myPatients:
                j.update()
            counter += 1

        # one final result for each patient each trial
        statsStep = []
        for z in myPatients:
            statsStep.append(z.getTotalPop())
        print statsStep
        virusPops.append(statsStep)
    # for hist to work, list should convert to numpy.array()
    virusPops = numpy.array(virusPops)
    stepsDelayed = [str(i) for i in stepsDelayed] # convert ints in list to strs
    pylab.hist(virusPops, 8, histtype='bar', color=['blue', 'red', 'green', 'purple'],
    label=stepsDelayed)
    pylab.xlabel('Virus cells remaining')
    pylab.ylabel('Patients')
    pylab.title('Outcomes for delayed treatment w/ 2 drugs')
    pylab.legend() # labels won't appear without it
    pylab.show()

#simulationTwoDrugsDelayedTreatment([0,75,150,300], 10)

#
# PROBLEM 5
#    

def simulationTwoDrugsVirusPopulations():

    """

    Run simulations and plot graphs examining the relationship between
    administration of multiple drugs and patient outcome.
    Plots of total and drug-resistant viruses vs. time are made for a
    simulation with a 300 time step delay between administering the 2 drugs and
    a simulations for which drugs are administered simultaneously.        

    """
    #TODO
    # go with stacked graphs
    # scatter plot for multiple trials
    # draw a best fit line
    # groupA 150, guttagonal, 300, grimpex, 150
    # groupB 150, guttagonal + grimpex, 150
    
    a = ResistantVirus(0.1, 0.05, {'guttagonal':False, 'grimpex':False}, 0.005)
    viruses = [a]*100
    patientsA = []
    patientsB = []
    groupA = []
    groupB = []
    for i in range(0,100):
        patientsA.append(Patient(viruses, 1000))
        patientsB.append(Patient(viruses, 1000))        
    for i in range(0,150):
        statsA = []
        statsB = []
        for j in patientsA:
            statsA.append(j.update())
        for j in patientsB:
            statsB.append(j.update())
        if i % 5 == 0:
            groupA.append(statsA)
            groupB.append(statsB)
    # prescribe guttagonal
    for j in patientsA:
        j.addPrescription('guttagonal')
    for k in patientsB:
        k.addPrescription('guttagonal')
        k.addPrescription('grimpex')
    for i in range(0,300):
        stats = []
        for j in patientsA:
            stats.append(j.update())
        if i % 5 == 0:
            groupA.append(stats)
    # prescribe grimpex
    for m in patientsA:
        m.addPrescription('grimpex')
    for i in range(0,150):
        statsA = []
        statsB = []
        for j in patientsA:
            statsA.append(j.update())
        for j in patientsB:
            statsB.append(j.update())
        if i % 5 == 0:
            groupA.append(statsA)
            groupB.append(statsB)
    # no virus left in simultaneous group
    # 25 patients had > 0 (e.g. 1-2) virus count remaining in delayed group
    print 100 - groupA[-1].count(0)
    print 100 - groupB[-1].count(0)
    # averaged line
    axisA = [sum(i)/len(i) for i in groupA]
    pylab.subplot(2, 1, 1)
    pylab.xlim([20,60])
    pylab.ylabel("Virus count")
    pylab.title('Two stage prescriptions')
    pylab.plot(groupA, 'bo')
    pylab.plot(axisA, 'r', linewidth=4)
    # stacked graphs
    axisB = [sum(i)/len(i) for i in groupB]
    pylab.subplot(2, 1, 2)
    pylab.xlim([20,60])
    pylab.title('Simultaneous prescriptions')
    pylab.xlabel("Time")
    pylab.plot(groupB, 'go')
    pylab.plot(axisB, 'r', linewidth=4)    
    pylab.show()
'''
    cannot tell a difference between one drug and two
    so printed the final numbers at line 501
'''    
#simulationTwoDrugsVirusPopulations()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

# 6.00 Problem Set 8

# Name: Casey

# Collaborators:

# Time: ~8hr

import numpy

import random

import pylab

from ps7 import *

# PROBLEM 1

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, 'grimpex',False}, means that this virus

particle is resistant to neither guttagonol nor grimpex.

mutProb: Mutation probability for this virus particle (a float). This is

the probability of the offspring acquiring or losing resistance to a drug.

"""

# TODO

assert 0 < maxBirthProb < 1 and 0 < clearProb < 1 and 0 < mutProb < 1

assert isinstance(resistances, dict)

self.maxBirthProb = maxBirthProb

self.clearProb = clearProb

self.resistances = resistances

self.mutProb = mutProb

def isResistantTo(self, drug):

"""

Get the state of this virus particle's resistance to a drug. This method

is called by getResistPop() in Patient 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.

"""

# TODO

# resistances: A dictionary of drug names (strings)

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 Patient class.

If the virus particle is not resistant to any drug in activeDrugs,

then it does not reproduce. Otherwise, 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 ResistantVirus (which has the same

maxBirthProb and clearProb values as its 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

grimpex, 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

grimpex and a 90% chance that the offspring will not be resistant to

grimpex.

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.

"""

# TODO

# if not resistant to a drug in activeDrugs: don't reproduce

# else: reproduce with birthProb and clearProb

resistant = False

newResistances = self.resistances

# if no active drugs, go to random.random()

if len(activeDrugs) > 0:

for i in activeDrugs:

if self.isResistantTo(i):

resistant = True

if resistant:

for k,v in self.resistances.items():

if random.random() < self.mutProb:

# reverse resistance in newResistances

newResistances[k] = not v

else:

return NoChildException()

if random.random() < self.maxBirthProb * (1 - popDensity):

return ResistantVirus(self.maxBirthProb, self.clearProb, newResistances, self.mutProb)

else:

return NoChildException()

#a = ResistantVirus(.99, .25, {'guttagonal': True, 'grimpex': True}, .99)

#b = a.reproduce(.1, ['grimpex'])

#print b.resistances

class Patient(SimplePatient):

"""

Representation of a patient. The patient is able to take drugs and his/her

virus population can acquire resistance to the drugs he/she takes.

"""

def __init__(self, viruses, maxPop):

"""

Initialization function, saves the viruses and maxPop parameters as

attributes. Also initializes the list of drugs being administered

(which should initially include no drugs).

viruses: the list representing the virus population (a list of

SimpleVirus instances)

maxPop: the maximum virus population for this patient (an integer)

"""

# TODO

assert isinstance(viruses, list)

assert isinstance(maxPop, int)

self.viruses = viruses

self.maxPop = maxPop

def addPrescription(self, newDrug):

"""

Administer a drug to this patient. After a prescription is added, the

drug acts on the virus population for all subsequent time steps. If the

newDrug is already prescribed to this patient, the method has no effect.

newDrug: The name of the drug to administer to the patient (a string).

postcondition: list of drugs being administered to a patient is updated

"""

# TODO

# should not allow one drug being added to the list multiple times

# if patient has no drug list: make one

# if newDrug on list: don't add

# else: add drug

assert isinstance(newDrug, str)

try:

if newDrug in self.drugs:

return

self.drugs.append(newDrug)

except AttributeError:

self.drugs = []

self.drugs.append(newDrug)

def getPrescriptions(self):

"""

Returns the drugs that are being administered to this patient.

returns: The list of drug names (strings) being administered to this

patient.

"""

# TODO

try:

return self.drugs

except AttributeError:

return []

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', 'grimpex'])

returns: the population of viruses (an integer) with resistances to all

drugs in the drugResist list.

"""

# TODO

assert isinstance(drugResist, list)

count = 0

for i in drugResist:

count += self.viruses.count(i)

return count

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().

- Determine whether each virus particle should reproduce and add

offspring virus particles to the list of viruses in this patient.

The listof 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)

"""

# TODO

survive_list = []

for i in range(0, len(self.viruses)):

if not self.viruses[i].doesClear():

survive_list.append(self.viruses[i])

# debug print "survive list: ", len(survive_list)

self.viruses = survive_list

# update the pop density

now_density = self.getTotalPop() / float(self.maxPop)

# for each virus particle, reproduce or not

new_viruses = []

for i in range(0, len(self.viruses)):

virus = self.viruses[i]

a = virus.reproduce(now_density, self.getPrescriptions())

if isinstance(a, ResistantVirus):

new_viruses.append(a)

# debug print "new viruses: ", len(new_viruses)

self.viruses = self.viruses + new_viruses

# return total virus pop as an int

return self.getTotalPop()

# PROBLEM 2

def simulationWithDrug():

"""

Runs simulations and plots graphs for problem 4.

Instantiates a patient, runs a simulation for 150 timesteps, adds

guttagonol, and runs the simulation for an additional 150 timesteps.

total virus population vs. time and guttagonol-resistant virus population

vs. time are plotted

"""

# TODO

a = ResistantVirus(0.1, 0.05, {'guttagonal': False}, 0.005)

viruses = [a]*100

trialPatient = Patient(viruses, 1000)

controlPatient = Patient(viruses, 1000)

guttagonal = []

control = []

for i in range(0, 150):

guttagonal.append(trialPatient.update())

control.append(controlPatient.update())

trialPatient.addPrescription('guttagonal')

for i in range(0, 150):

guttagonal.append(trialPatient.update())

control.append(controlPatient.update())

pylab.plot(guttagonal, 'r-', label='guttagonal prescribed')

pylab.plot(control, 'b-', label='without medication')

pylab.xlabel("Time")

pylab.ylabel("Virus count")

pylab.legend(loc='best')

pylab.title("Virus when drugs applied at time 150")

pylab.show()

#simulationWithDrug()

'''

virus *was* dying off before we get to timestep 150

was not resistant to empty activeDrugs

'''

# PROBLEM 3

def simulationDelayedTreatment(stepsDelayed, numTrials):

"""

Runs simulations and make histograms for problem 5.

Runs multiple simulations to show the relationship between delayed treatment

and patient outcome.

Histograms of final total virus populations are displayed for delays of 300,

150, 75, 0 timesteps (followed by an additional 150 timesteps of

simulation).

"""

# TODO

# run 30 trials for each scenario

# displays histograms of final virus counts as four colored bars

virusPops = []

a = ResistantVirus(0.1, 0.05, {'guttagonal': True}, 0.005)

viruses = [a]*100

for k in range(0, numTrials):

# reset patients for each trial

myPatients = []

for i in range(0, len(stepsDelayed)):

a = Patient(viruses, 1000)

myPatients.append(a)

# run through all steps with delayed prescriptions

# capture final viral counts

statsStep = []

counter = 0

for i in range(0, max(stepsDelayed)+150):

if counter in stepsDelayed:

myPatients[stepsDelayed.index(counter)].addPrescription('guttagonal')

for j in myPatients:

j.update()

counter += 1

# one final result for each patient each trial

for z in myPatients:

statsStep.append(z.getTotalPop())

print statsStep

virusPops.append(statsStep)

# for hist to work, list should convert to numpy.array()

virusPops = numpy.array(virusPops)

'''

README

100 trials shows 389 patients free of virus and 11 with 1 virus left

needs legend

and axis labels

'''

stepsDelayed = [str(i) for i in stepsDelayed] # convert ints in list to strs

pylab.hist(virusPops, 8, histtype='bar', color=['blue', 'red', 'green', 'purple'],

label=stepsDelayed)

pylab.xlabel('Virus cells remaining')

pylab.ylabel('Patients')

pylab.title('Outcomes for delayed treatment')

pylab.legend() # labels won't appear without it

pylab.show()

#simulationDelayedTreatment([0,75,150,300], 100)

# PROBLEM 4

def simulationTwoDrugsDelayedTreatment(stepsDelayed, numTrials):

"""

Runs simulations and make histograms for problem 6.

Runs multiple simulations to show the relationship between administration

of multiple drugs and patient outcome.

Histograms of final total virus populations are displayed for lag times of

150, 75, 0 timesteps between adding drugs (followed by an additional 150

timesteps of simulation).

"""

# TODO

# run 30 trials for each scenario

# plot

virusPops = []

a = ResistantVirus(0.1, 0.05, {'guttagonal': False, 'grimpex': False}, 0.005)

viruses = [a]*100

for k in range(0, numTrials):

# reset patients for each trial

myPatients = []

for i in range(0, len(stepsDelayed)):

a = Patient(viruses, 1000)

myPatients.append(a)

# run an initial 150 steps prior to guttagonal

for i in range(0, 150):

for j in myPatients:

j.update()

for j in myPatients:

j.addPrescription('guttagonal')

# apply grimpex at delays defined in stepsDelayed

counter = 0

for i in range(0, max(stepsDelayed)+150):

if counter in stepsDelayed:

myPatients[stepsDelayed.index(counter)].addPrescription('grimpex')

for j in myPatients:

j.update()

counter += 1

# one final result for each patient each trial

statsStep = []

for z in myPatients:

statsStep.append(z.getTotalPop())

print statsStep

virusPops.append(statsStep)

# for hist to work, list should convert to numpy.array()

virusPops = numpy.array(virusPops)

stepsDelayed = [str(i) for i in stepsDelayed] # convert ints in list to strs

pylab.hist(virusPops, 8, histtype='bar', color=['blue', 'red', 'green', 'purple'],

label=stepsDelayed)

pylab.xlabel('Virus cells remaining')

pylab.ylabel('Patients')

pylab.title('Outcomes for delayed treatment w/ 2 drugs')

pylab.legend() # labels won't appear without it

pylab.show()

#simulationTwoDrugsDelayedTreatment([0,75,150,300], 10)

# PROBLEM 5

def simulationTwoDrugsVirusPopulations():

"""

Run simulations and plot graphs examining the relationship between

administration of multiple drugs and patient outcome.

Plots of total and drug-resistant viruses vs. time are made for a

simulation with a 300 time step delay between administering the 2 drugs and

a simulations for which drugs are administered simultaneously.

"""

#TODO

# go with stacked graphs

# scatter plot for multiple trials

# draw a best fit line

# groupA 150, guttagonal, 300, grimpex, 150

# groupB 150, guttagonal + grimpex, 150

a = ResistantVirus(0.1, 0.05, {'guttagonal':False, 'grimpex':False}, 0.005)

viruses = [a]*100

patientsA = []

patientsB = []

groupA = []

groupB = []

for i in range(0,100):

patientsA.append(Patient(viruses, 1000))

patientsB.append(Patient(viruses, 1000))

for i in range(0,150):

statsA = []

statsB = []

for j in patientsA:

statsA.append(j.update())

for j in patientsB:

statsB.append(j.update())

if i % 5 == 0:

groupA.append(statsA)

groupB.append(statsB)

# prescribe guttagonal

for j in patientsA:

j.addPrescription('guttagonal')

for k in patientsB:

k.addPrescription('guttagonal')

k.addPrescription('grimpex')

for i in range(0,300):

stats = []

for j in patientsA:

stats.append(j.update())

if i % 5 == 0:

groupA.append(stats)

# prescribe grimpex

for m in patientsA:

m.addPrescription('grimpex')

for i in range(0,150):

statsA = []

statsB = []

for j in patientsA:

statsA.append(j.update())

for j in patientsB:

statsB.append(j.update())

if i % 5 == 0:

groupA.append(statsA)

groupB.append(statsB)

# no virus left in simultaneous group

# 25 patients had > 0 (e.g. 1-2) virus count remaining in delayed group

print 100 - groupA[-1].count(0)

print 100 - groupB[-1].count(0)

# averaged line

axisA = [sum(i)/len(i) for i in groupA]

pylab.subplot(2, 1, 1)

pylab.xlim([20,60])

pylab.ylabel("Virus count")

pylab.title('Two stage prescriptions')

pylab.plot(groupA, 'bo')

pylab.plot(axisA, 'r', linewidth=4)

# stacked graphs

axisB = [sum(i)/len(i) for i in groupB]

pylab.subplot(2, 1, 2)

pylab.xlim([20,60])

pylab.title('Simultaneous prescriptions')

pylab.xlabel("Time")

pylab.plot(groupB, 'go')

pylab.plot(axisB, 'r', linewidth=4)

pylab.show()

'''

cannot tell a difference between one drug and two

so printed the final numbers at line 501

'''

#simulationTwoDrugsVirusPopulations()