%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from datetime import datetime, timedelta
from collections import defaultdict

full_life_table = pd.read_csv('data/life_table.csv')
full_life_table.set_index('Age',inplace=True)

full_life_table[80:90]

# grab the values for the "Male" column
prob_death = full_life_table.M.values

# Create a new life table
l = len(prob_death)
life_table = np.ones((l,l))
for i in range(0,l,1):
    for j in range(i,l,1):
        life_table[i][j]=np.prod(1 - prob_death[i:j+1])

# Load current apostle ages and seniority
apostle_data = pd.read_csv('data/apostles.csv')
apostle_data

class Apostle:
    def __init__(self,name,birth,twelve,ordained,seniority):
        self.name=name
        self.birth=datetime.strptime(birth,'%m/%d/%Y')
        self.twelve=datetime.strptime(twelve,'%m/%d/%Y')
        self.ordained=datetime.strptime(ordained,'%m/%d/%Y')
        self.seniority=seniority
    
    def __str__(self):
        return self.name
    
    def age(self,year):
        """Return the apostle's age in a particular year."""
        return year-self.birth.year
    
    def current_age(self):
        """Return the apostle's age in the current year."""
        current_year = datetime.now().year
        return self.age(current_year)
    
    def most_probable_death_year(self,life_table):
        """Return the apostle's most likely year of death given
        a particular life table."""
        current_age = self.current_age()
        age_of_death = np.argmax(life_table[current_age]<.5)
        return self.birth.year + age_of_death
    
    def most_probable_life_state(self,year,life_table):
        """Return True if, given the most probable year of death,
        the apostle is alive in any part of a particular year
        given a particular life table."""
        return year <= self.most_probable_death_year(life_table)
    
    def simulate_death_year(self,life_table):
        """Return a death year randomly drawn from the distribution
        of likely death years given a particular life table."""
        death_year = np.argmin(life_table[self.current_age()] > np.random.random(len(life_table[0])))
        return self.birth.year + death_year
    
    def simulate_life_state(self,year,life_table):
        """Return True if the apostle is alive given a year of death
        drawn from the distribution of likely death years given
        by a particular life table."""
        life_state = self.simulate_death_year(life_table)
        return year <= self.simulate_death_year(life_table)
        

# Create apostle objects and print their ages
apostles = []
for i,row in apostle_data.iterrows():
    apostle = Apostle(row.Name,row.Birth,row.Twelve,row.Ordained,row.Seniority)
    apostles.append(apostle)
    print "{}. {} is {} years old".format(apostle.seniority,apostle.name,apostle.current_age())
    

# define some helper functions for calculating who is prophet in a particular year
# using two different methods

def most_probable_prophet_in_year(apostles,year,life_table):
    """Return the apostle (given a list of apostles) most likely
    to be prophet in a particular year"""
    apostles_alive = [apostle for apostle in apostles
                      if apostle.most_probable_life_state(year,life_table)]
    
    if len(apostles_alive) == 0:
        return None
    
    apostle_index = np.argmin([apostle.seniority for apostle in apostles_alive])
    return apostles_alive[apostle_index]

def simulate_prophet_in_year(apostles,year,life_table):
    """Return the apostle (given a list of apostles) who is prophet
    in a particular year after simulating each apostle's life state
    in the given year."""
    apostles_alive = [apostle for apostle in apostles
                      if apostle.simulate_life_state(year,life_table)]
    if len(apostles_alive) == 0:
        return None
    
    apostle_index = np.argmin([apostle.seniority for apostle in apostles_alive])
    return apostles_alive[apostle_index]

# Plot a histogram of each apostle's likely death years
for apostle in apostles:
    death_year_dist = []
    for i in range(10000):
        death_year_dist.append(apostle.simulate_death_year(life_table))

    plt.hist(death_year_dist,bins=range(2014,2040))
    plt.title("Histogram of year of death of {}".format(apostle.name))
    plt.show()

# Given our model, who is most likely to be prophet in the year 2020?
print most_probable_prophet_in_year(apostles,2020,life_table)

# Given our model, who is most likely to be prophet in every year from now until 2040?
for year in range(2014,2040):
    print "{}: {}".format(year, most_probable_prophet_in_year(apostles,year,life_table))

def run_simulation(n_trials,apostles,life_table):
    for year in range(2014,2040):
        prophets = defaultdict(list)
        for i in range(n_trials):
            prophet = simulate_prophet_in_year(apostles,year,life_table)
            if prophet is not None:
                prophets[prophet.name].append(1)
            else:
                prophets["other"].append(1)
        probabilities = []
        for key,count in prophets.items():
            probabilities.append((key,float(len(count))/n_trials))
        probabilities = sorted(probabilities,key=lambda x: x[1],reverse=True)
        print year
        for name,probability in probabilities:
            print " {} ({})".format(name,probability)
        print

run_simulation(10000,apostles,life_table)

default = 8
apostle_age_adjustments = {
    "Thomas S. Monson":   default - 2,
    "Boyd K. Packer":     default - 2,
    "L. Tom Perry":       default + 1,
    "Russell M. Nelson":  default + 2,
    "Dallin H. Oaks":     default + 1,
    "M. Russell Ballard": default,
    "Richard G. Scott":   default - 2,
    "Robert D. Hales":    default + 2,
    "Jeffrey R. Holland": default,
    "Henry B. Eyring":    default,
    "Dieter F. Uchtdorf": default,
    "David A. Bednar":    default,
    "Quentin L. Cook":    default,
    "D. Todd Christofferson": default,
    "Neil L. Andersen":   default}

adj_apostles = []
for apostle in apostles:
    age_adj = apostle_age_adjustments[apostle.name]
    adj_apostles.append(Apostle(apostle.name,
                                (apostle.birth + timedelta(days=age_adj*365)).strftime("%m/%d/%Y"),
                                apostle.twelve.strftime("%m/%d/%Y"),
                                apostle.ordained.strftime("%m/%d/%Y"),
                                apostle.seniority))

run_simulation(10000,adj_apostles,life_table)