%%capture
# Time Access and Conversionsi
import time

# Unix Pattern Extensions
import glob

# System-specific Parameters and Functions
import sys

# Operating System Interfaces
import os

# Plotting 
import matplotlib.pyplot as plt
%pylab inline --no-import-all

# Scientific Computing
import numpy as np

# Notebook Displays
from IPython.display import display, HTML, Image

# Importing the grmpy package by editing the PYTHONPATH 
sys.path.insert(0, '../grmpy')

# Package
import grmpy as gp

# Hidden function
from tests._auxiliary import random_init
from tools.user.processing import _add_auxiliary

"""
Explanation: Private Package
Packages are a way to structure Python namespaces and help to orgainze our code in logical units. This allows makes our code more comprehensible and extensible.  
Roadmap

Structuring the Package
Exploring Performance
Profiling
Quality Assurance

Exploring the Package
This is a first look at the structure of the grmpy package.
grmpy/
  tests/               # This directory collects all functions
    __init__.py        # related to the testing of the package.
    _auxiliary.py      #

  tools/
    economics/         # All functions related to the economics
      __init__.py      # of our version of the generalized Roy
      simulation.py    # model.

    optimization/      # All functions related to the maximum
      __init__.py      # likelihood estimation of the model.
      estimation.py

    user/              # Some functions that allow the user to 
      __init__.py      # to interact with the toolbox.
      inspection.py
      processing.py

  __init__.py

Let us take a look at around the package in PyCharm in more detail. For example, we will study how to actually distribute a Python package across several files and folders. Throughout, we will encounter some very useful feature of our IDE that allow to navigate our Python package.
Before we get started using the package, let us import some basic libraries.
End of explanation
"""


# Process initialization file
init_dict = gp.process('material/msc/init.ini')

# Simulate synthetic sample
gp.simulate(init_dict)

# Estimate model
rslt = gp.estimate(init_dict)

# Write results
gp.inspect(rslt, init_dict)

# Inspect the results
%cat results.grm.txt

"""
Explanation: Now that we have an idea about the structure of the package, let us make some simple calls to the public functions.
End of explanation
"""


display(HTML('material/images/grmpy.html')) 

"""
Explanation: From the perspective of the user, actually nothing changed compared to our implementation of a simple module. How do we get access to the relevant set of functions?
End of explanation
"""


%%capture 
# Magic function to suppress output

# Ensure recomputability
np.random.seed(7892)

# Set the number of tests to run
NUM_TESTS = 10

# Run repeated tests
for _ in range(NUM_TESTS):

    # Generate random request
    init_dict = random_init()

    # Ensure same starting value. If we choose the random 
    # starting values instead, the points of evaluation 
    # differ for the slow and fast implementations.
    init_dict['ESTIMATION']['start'] = 'init'
    init_dict['ESTIMATION']['maxiter'] = 0

    # Simulate sample
    gp.simulate(init_dict)
    
    # Initialize result container
    rslt = dict()
    
    # Estimate generalized Roy model
    for version in ['slow', 'fast', 'object']:

        init_dict['ESTIMATION']['version'] = version

        rslt[version] = gp.estimate(init_dict)['fval']
    
    # Assert equality of results
    np.testing.assert_allclose(rslt['slow'], rslt['fast'])
    np.testing.assert_allclose(rslt['slow'], rslt['object'])
    
    # Cleanup
    os.remove(init_dict['BASICS']['source'])

"""
Explanation: Exploring Performance
There are several implementations of the likelihood function available: slow, fast, and object-oriented. We can select the implementation to use in the initialization file. Let us make sure that the result is the same regardless of the implementation.
End of explanation
"""


%%capture 
# Magic function to suppress output

# Set grid for varying number of agents
agent_grid = range(1000, 10000, 500)

# Generate random request
init_dict = random_init()

# Request single function evaluation
init_dict['ESTIMATION']['maxiter'] = 0

# Simulate sample
rslt = dict()

for version in ['slow', 'fast', 'object']:
    
    # Initialzie container for results
    rslt[version] = []
    
    # Update initialization file
    init_dict['ESTIMATION']['version'] = version
    
    for agents in agent_grid:
        
        # Update initialization file
        init_dict['BASICS']['agents'] = agents
        
        # Simulate sample
        gp.simulate(init_dict)
        
        # Measure execution time of for estimation
        start_time = time.time()

        gp.estimate(init_dict)

        rslt[version] += [time.time() - start_time]

"""
Explanation: Let's see how the execution time differes as we increase the number of observations.
End of explanation
"""


# Initialize canvas
ax = plt.figure(figsize=(12,8)).add_subplot(111)

# Plot execution times by implementations 
ax.plot(agent_grid, rslt['slow'], label='Slow')
ax.plot(agent_grid, rslt['fast'], label='Fast')
ax.plot(agent_grid, rslt['object'], label='Object')

# Set axis labels
ax.set_xlabel('Simulated Agents', fontsize=20)
ax.set_ylabel('Execution Time (Seconds)', fontsize=20)

# Change background color
ax.set_axis_bgcolor('white')

# Set up legend
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.10),
    fancybox=False, frameon=False, shadow=False, ncol=3,
    fontsize=20)
     
# Remove first element on y-axis
ax.yaxis.get_major_ticks()[0].set_visible(False)

"""
Explanation: Now we can easily visualize the difference in performance.
End of explanation
"""


# Import package for profiling 
import cProfile

# Process initialization file
init_dict = gp.process('material/msc/init.ini')

# Modify estimation request
init_dict['ESTIMATION']['maxiter'] = 0
init_dict['BASICS']['agents'] = 10000

# Simulate synthetic sample
gp.simulate(init_dict)

# Profiling the fast and slow implementation
for version in ['fast', 'slow']:
    
    # Select implementation
    init_dict['ESTIMATION']['version'] = version

    # Run profiler and print to screen
    cProfile.run('gp.estimate(init_dict)')

    # Run profiler and print to file
    cProfile.run('gp.estimate(init_dict)', 'grmpy_' + version + '.grm.prof')


"""
Explanation: What implementation choices are driving these differences in performance?
Profiling
Profiling tools allow to determine the performance of your computer program and identify the frequency and duration of function calls. Often just slight modifications make lead to more rapid executions and less memory use. Profiling our package is straightforward:
End of explanation
"""


Image(filename='material/images/runsnake_slow.png')

"""
Explanation: More details about profiling Python code is available here.  I also use RunSnakeRun to visualize the information provided by the profiler. PyCharm just added an integreated Python profiler, which is probably worth exploring.
End of explanation
"""


display(HTML('material/images/criterion_slow.html')) 

"""
Explanation: With this information at hand, let us revist our baseline implementation for the likelihood function.
End of explanation
"""


display(HTML('material/images/criterion_fast.html')) 

"""
Explanation: Now we turn to the faster alternative.
End of explanation
"""


Image(filename='material/images/runsnake_fast.png')

"""
Explanation: How about now?
End of explanation
"""


%%capture 
# Magic function to suppress output

# Ensure recomputability 
np.random.seed(7516)

# Set the number of tests to run
NUM_TESTS = 10

# Run repeated tests
for _ in range(NUM_TESTS):

    # Generate random request
    init_dict = random_init()
    
    # Let us restrict attention to either zero or just one 
    # step of the estimator.
    init_dict['ESTIMATION']['maxiter'] = np.random.choice([0, 1])
    
    # Simulate sample
    gp.simulate(init_dict)
    
    # Estimate model
    rslt = gp.estimate(init_dict)

    # Write results
    gp.inspect(rslt, init_dict)

"""
Explanation: Quality Assurance
Is our code working for all the different version of the model we want to be able to estimate? Often, parts of the code are only visited in special cases which results in errors remaining undetected for quite a while. Let us just hit the package with some random requests and see if it all runs. Later, of course, we will focus on whether it is actually working correctly.
End of explanation
"""


# Ensure recomputability 
np.random.seed(7523216)

# Generate random request
init_dict = random_init()

# Let us make sure to have a setup that is very favourable to
# the performance of our estimator for now. In my experience,
# small unobserved variability in agent choices and outcomes
# and a large agent count does the trick. Of course, later 
# we will investigate the performance of your estimator
# for more challenging tasks.
for key_ in ['COST', 'TREATED', 'UNTREATED']:
    init_dict[key_]['sd'] = 0.05

init_dict['BASICS']['agents'] = 10000

# We need to ensure that the random request actually entails
# a serious estimation run.
init_dict['ESTIMATION']['maxiter'] = 1000
init_dict['ESTIMATION']['start'] = 'auto'
init_dict['ESTIMATION']['version'] = 'fast'
init_dict['ESTIMATION']['optimizer'] = 'bfgs'

# We need to update the auxiliary entries in the initialization
# dictionary. Otherwise the `init_values' do not entry does
# not correspond to the true values.
init_dict = _add_auxiliary(init_dict)

# Simulate synthetic sample
gp.simulate(init_dict)

# Estimate model
rslt = gp.estimate(init_dict)

# Check for convergence, which is very different from our 
# notion of a SUCCESSFUL estimation run. 
assert (rslt['success'] is True)

# Write results
gp.inspect(rslt, init_dict)

# Inspect the results
%cat results.grm.txt

"""
Explanation: Now, let us finally have a look whether the program actually recovers the true structural parameters.
End of explanation
"""


# Load true and estimated parameters
true_values = rslt['AUX']['init_values']
est_values = rslt['AUX']['x_internal']
start_values = init_dict['AUX']['start_values']
    
# Auxiliary objects
num_paras = len(true_values)

# Formatting
fmt = '{0:10.2f}{1:10.2f}{2:10.2f}{3:10.2f}'

# Print both parameters
print '     Start  Estimate     Truth    Difference \n'
for i in range(num_paras):
            
    start, est, true = start_values[i], est_values[i], true_values[i]

    diff = est - true
    
    print fmt.format(start, est, true, diff)


"""
Explanation: How does this compare to some of the structural parameters we used for the simulation of the dataset?
End of explanation
"""


def rmse(rslt):
    """ Calculate the root-mean squared error.
    """
    # Antibugging
    assert (isinstance(rslt, dict))

    # Distribute information
    x_internal = rslt['AUX']['x_internal']
    init_values = rslt['AUX']['init_values']

    # Calculate statistic
    rslt = ((x_internal - init_values) ** 2).mean()

    # Antibugging
    assert (np.isfinite(rslt))
    assert (rslt > 0.0)

    # Finishing
    return rslt

%%capture 
# Magic function to suppress output

# Ensure recomputability
np.random.seed(456)

# Set grid for varying degree of noise
noise_grid = range(20)

init_dict = random_init()

# Modify estimation request ensuring fast executation, 
# large sample size, and very favourable starting values.
init_dict['BASICS']['agents'] = 1000
init_dict['ESTIMATION']['version'] = 'fast'
init_dict['ESTIMATION']['maxiter'] = 10000
init_dict['ESTIMATION']['start'] = 'auto'
        
# Simulate synthetic sample
rslt = dict()

for algorithm in ['bfgs', 'nm']:

    # Initialize containers
    rslt[algorithm] = []

    # Ensure same simulated setup
    np.random.seed(123)

    for i in noise_grid:

        # Increase noise in observed sample
        for key_ in ['COST', 'TREATED', 'UNTREATED']:
            init_dict[key_]['sd'] = 0.05 + i*0.50

        # We need to update the auxiliary entries in the initialization
        # dictionary. Otherwise the `init_values' do not entry does
        # not correspond to the true values.
        init_dict = _add_auxiliary(init_dict)

        # Simulate dataset
        gp.simulate(init_dict)

        # Select estimation setup
        init_dict['ESTIMATION']['algorithm'] = algorithm

        # Calculate performance statistic
        stat = rmse(gp.estimate(init_dict)) 

        # Collect results
        rslt[algorithm] += [stat]

"""
Explanation: Let us repeat this exercise, but use the Nelder Mead algorithm instead of BFGS. How does performance differ? What about using different starting values?
Let's investigate the performance of our estimator as we increase the unobserved variability in the choice and outcome equation. We will judge quality of the estimations by the Root-Mean-Squared Error (RMSE).
End of explanation
"""


# Initialize canvas
ax = plt.figure(figsize=(12,8)).add_subplot(111)

# Plot execution times by implementations 
ax.plot(noise_grid, rslt['bfgs'], label='BFGS')
ax.plot(noise_grid, rslt['nm'], label='Nelder-Mead')

# Set axis labels
ax.set_xlabel('Unobserved Variability', fontsize=20)
ax.set_ylabel('Root-Mean-Squared Error', fontsize=20)

# Change background color
ax.set_axis_bgcolor('white')

# Set up legend
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.10),
    fancybox=False, frameon=False, shadow=False, ncol=2,
    fontsize=20)
     
# Remove first element on y-axis
ax.yaxis.get_major_ticks()[0].set_visible(False)

"""
Explanation: Let us visualize the performance.
End of explanation
"""


%%capture 
# Magic function to suppress output

# Ensure recomputability
np.random.seed(456)

# Set grid for varying degree of noise
noise_grid = range(20)

init_dict = random_init()

init_dict['BASICS']['agents'] = 1000
init_dict['ESTIMATION']['version'] = 'fast'
init_dict['ESTIMATION']['maxiter'] = 10000
init_dict['ESTIMATION']['start'] = 'auto'
        
# Simulate synthetic sample
rslt = dict()

for algorithm in ['bfgs', 'nm']:

    # Initialize containers
    rslt[algorithm] = []

    # Ensure same simulated setup
    np.random.seed(123)

    for i in noise_grid:

        # Increase noise in observed sample
        for key_ in ['COST', 'TREATED', 'UNTREATED']:
            init_dict[key_]['sd'] = 0.05 + i*0.50

        # Simulate dataset
        gp.simulate(init_dict)

        # Select estimation setup
        init_dict['ESTIMATION']['algorithm'] = algorithm

        # Measure execution time of for estimation
        start_time = time.time()

        gp.estimate(init_dict)

        # Collect results
        rslt[algorithm] += [time.time() - start_time]

"""
Explanation: In this very favourable setting (correctly specified model, small unobserved variability in agent outcomes and choices, large sample), both optimzers perform quite well. But what about the total time for an estimation run?
End of explanation
"""


# Initialize canvas
ax = plt.figure(figsize=(12,8)).add_subplot(111)

# Plot execution times by implementations 
ax.plot(noise_grid, rslt['bfgs'], label='BFGS')
ax.plot(noise_grid, rslt['nm'], label='Nelder-Mead')

# Set axis labels
ax.set_xlabel('Unobserved Variability', fontsize=20)
ax.set_ylabel('Execution Time (Seconds)', fontsize=20)

# Change background color
ax.set_axis_bgcolor('white')

# Set up legend
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.10),
    fancybox=False, frameon=False, shadow=False, ncol=2,
    fontsize=20)
     
# Remove first element on y-axis
ax.yaxis.get_major_ticks()[0].set_visible(False)

"""
Explanation: So, how different are they and does relative performance differ as we increase the unobserved variability?
End of explanation
"""


# Create list of all files generated by the grmpy package
files = glob.glob('*.grm.*')

# Remove files
for file_ in files:
    
    os.remove(file_)

"""
Explanation: The relative performance of the BFGS algorithm deteriorates as we increase the unobserved variability. The noise has a negative influence on the numerical approximation of the derivatives of our criterion function. The thus provide less reliable information for the optimizer's next step.
Cleanup
This concludes this lecture and all that remains done it to clean up.
End of explanation
"""


import urllib; from IPython.core.display import HTML
HTML(urllib.urlopen('http://bit.ly/1OKmNHN').read())

"""
Explanation: Formatting
End of explanation
"""


