#!/usr/bin/env python
# coding: utf-8

# # Advanced Examples
# ## Making a Monte Carlo parameter study
# In this example, we will make a Monte Carlo study. This is a case where the python library shows it's advantage. 
# 
# We will make a Monte Carlo study on the position of patella tendon insertion and origin in the simplified knee model used in the first tutorial. Thus, we need some macros that change two `sRel` variables in the model. In this case, we choose the values from a truncated normal distribution, but any statistical distribution could have been used. 
# 

# In[1]:


from scipy.stats import distributions
# Truncated normal between +/- 2SD.
patella_tendon_insertion = distributions.truncnorm(-2,2,[0.02, 0.12, 0], [0.01,0.01,0.01]) 
patella_tendon_origin = distributions.truncnorm(-2,2,[0.0,-0.03, 0], [0.01,0.01,0.01]) 


# In[2]:


from anypytools import macro_commands as mc

macro = [
    mc.SetValue_random('Main.MyModel.Tibia.Patella2.sRel', patella_tendon_insertion),
    mc.SetValue_random('Main.MyModel.Patella.Lig.sRel', patella_tendon_origin)
]

macro 


# The default representation of the macro is 50% percentile values from the distribution (i.e. the mean value). To generate something random we need the `AnyMacro` helper class. 

# In[3]:


from anypytools import AnyMacro

mg = AnyMacro(macro)
monte_carlo_macros = mg.create_macros_MonteCarlo(100)


# The first generated macro just has the default values.

# In[4]:


monte_carlo_macros[0]


# The next two macros have random offsets. 

# In[5]:


monte_carlo_macros[1:4]


# Now let us expand the macro to also load and run the model. 

# In[6]:


macro = [
    mc.Load('Knee.any'),
    mc.SetValue_random('Main.MyModel.Tibia.Patella2.sRel', patella_tendon_insertion ) ,
    mc.SetValue_random('Main.MyModel.Patella.Lig.sRel', patella_tendon_origin ) ,
    mc.OperationRun('Main.MyStudy.InverseDynamics'),
    mc.Dump('Main.MyStudy.Output.Abscissa.t'),
    mc.Dump('Main.MyStudy.Output.MaxMuscleActivity')
]
mg = AnyMacro(macro, seed=1)
monte_carlo_macros = mg.create_macros_MonteCarlo(100)


# ### Running the Monte Carlo macro

# The macro is passed to the AnyPyProcess object which excutes the macros

# In[7]:


from anypytools import AnyPyProcess 

app = AnyPyProcess()
monte_carlo_results = app.start_macro( monte_carlo_macros )


# The output object (`monte_carlo_result`) is a list-like object where each element is a dictionary with the output of the corresponding simulation.
# 

# In[8]:


print('Length :', len(monte_carlo_results) )
print('Data keys in first element: ', list(monte_carlo_results[0].keys()))


# ### Filtering results with Errors
# If any model errors occured in some of the simulations, then the output may be missing. That can be a problem for futher processing. So we may want to remove those simulations. 
# 
# That is easily done. Results from simuations with errors will contain a  an 'ERROR' key. We can use that to filter the results.

# In[9]:


monte_carlo_results[:] = [
    output 
    for output in monte_carlo_results 
    if 'ERROR' not in output
]


# ### Extracting data from the resutls
# The object looks and behaves like a list, but it can do more things than a standard Python list. If we use it as a dictionary and pass the variable names, it will return that variable concatenated over all runs. 
# 
# 

# In[10]:


monte_carlo_results['Main.MyStudy.Output.MaxMuscleActivity']


# In fact, we don't even have to use the full name of the variable. As long as the string uniquely defines the variable in the data set.

# In[11]:


max_muscle_activity = monte_carlo_results['MaxMus']


# ### Plotting the results
# Finally we can plot the result of the Monte Carlo study

# In[12]:


get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt

time = monte_carlo_results['Abscissa.t'][0]
plt.fill_between(time, max_muscle_activity.min(0), max_muscle_activity.max(0),alpha=0.4  )
for trace in max_muscle_activity:
    plt.plot(time, trace,'b', lw=0.2 )
# Plot result with the mean of the inputs ( stored in the first run)
plt.plot(time, max_muscle_activity[0],'r--', lw = 2, ) 


# ## Latin Hypercube sampling

# Monte Carlo studies are not very efficient when investigating the effect of many parameters. It quickly becomes necessary to run the model thousands of times. Not very convenient if the AnyBody model takes a long time to run. 
# 
# Another approach is to use Latin Hypercube sampling. From Wikipedia
# > Latin hypercube sampling (LHS) is a statistical method for generating a sample of plausible collections of parameter values from a multidimensional distribution. The sampling method is often used to construct computer experiments.
# 
# Using LHS we can generate a sample that better spans the whole multidimensional space. Thus, fever model evaluations are necessary. See [pyDOE](http://pythonhosted.org/pyDOE/randomized.html) for examples (and explanation of the `criterion`/`iterations` parameters which can be parsed to `create_macros_LHS()`).
# 
# To following uses LHS to do the same as in the previous example:

# In[13]:


patella_tendon_insertion = distributions.norm([0.02, 0.12, 0], [0.01,0.01,0.01]) 
patella_tendon_origin = distributions.norm([0.0,-0.03, 0], [0.01,0.01,0.01]) 

macro = [ 
    mc.Load('Knee.any'),
    mc.SetValue_random('Main.MyModel.Tibia.Patella2.sRel', patella_tendon_insertion ) ,
    mc.SetValue_random('Main.MyModel.Patella.Lig.sRel', patella_tendon_origin ) ,
    mc.OperationRun('Main.MyStudy.InverseDynamics'),
    mc.Dump('Main.MyStudy.Output.Abscissa.t'),
    mc.Dump('Main.MyStudy.Output.MaxMuscleActivity')
]
mg = AnyMacro(macro, seed=1)
LHS_macros = mg.create_macros_LHS(25)


# In[14]:


from anypytools import AnyPyProcess 

app = AnyPyProcess()
lhs_results = app.start_macro(LHS_macros)


# In[15]:


get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt

time = lhs_results['Abscissa.t'][0]
plt.fill_between(time, lhs_results['MaxMus'].min(0), lhs_results['MaxMus'].max(0),alpha=0.4  )
for trace in lhs_results['MaxMus']:
    plt.plot(time, trace,'b', lw=0.2 )
# Plot the mean value that was stored in the first results
plt.plot(time, lhs_results['MaxMus'].mean(0),'r--', lw = 2, )