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
Out[2]:
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]
Out[4]:
The next two macros have random offsets.
In [5]:
monte_carlo_macros[1:4]
Out[5]:
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)
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()))
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
]
In [10]:
monte_carlo_results['Main.MyStudy.Output.MaxMuscleActivity']
Out[10]:
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']
In [12]:
%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, )
Out[12]:
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 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]:
%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, )
Out[15]: