{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Advanced Examples" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Making a Monte Carlo parameter study\n", "In this example, we will make a Monte Carlo study. This is a case where the python library shows it's advantage. \n", "\n", "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. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from scipy.stats import distributions\n", "# Truncated normal between +/- 2SD.\n", "patella_tendon_insertion = distributions.truncnorm(-2,2,[0.02, 0.12, 0], [0.01,0.01,0.01]) \n", "patella_tendon_origin = distributions.truncnorm(-2,2,[0.0,-0.03, 0], [0.01,0.01,0.01]) " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[classoperation Main.MyModel.Tibia.Patella2.sRel \"Set Value\" --value=\"{0.02,0.12,2.78291642467e-18}\",\n", " classoperation Main.MyModel.Patella.Lig.sRel \"Set Value\" --value=\"{2.78291642467e-18,-0.03,2.78291642467e-18}\"]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from anypytools import macro_commands as mc\n", "\n", "macro = [\n", " mc.SetValue_random('Main.MyModel.Tibia.Patella2.sRel', patella_tendon_insertion),\n", " mc.SetValue_random('Main.MyModel.Patella.Lig.sRel', patella_tendon_origin)\n", "]\n", "\n", "macro " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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. \n", "\n", "The `AnyMacro` helper class can generate multiple macros from a single macro." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from anypytools import AnyMacro\n", "\n", "mg = AnyMacro(macro)\n", "monte_carlo_macros = mg.create_macros_MonteCarlo(100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first generated macro just has the default values." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['classoperation Main.MyModel.Tibia.Patella2.sRel \"Set Value\" --value=\"{0.02,0.12,2.78291642467e-18}\"',\n", " 'classoperation Main.MyModel.Patella.Lig.sRel \"Set Value\" --value=\"{2.78291642467e-18,-0.03,2.78291642467e-18}\"']" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "monte_carlo_macros[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next two macros have random offsets. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[['classoperation Main.MyModel.Tibia.Patella2.sRel \"Set Value\" --value=\"{0.0144267872987,0.109982579448,-0.01776142425}\"',\n", " 'classoperation Main.MyModel.Patella.Lig.sRel \"Set Value\" --value=\"{0.00182301737147,-0.0385242952422,0.0183692935384}\"'],\n", " ['classoperation Main.MyModel.Tibia.Patella2.sRel \"Set Value\" --value=\"{0.0242051783405,0.103809162986,0.00739209161022}\"',\n", " 'classoperation Main.MyModel.Patella.Lig.sRel \"Set Value\" --value=\"{-0.00636881243957,-0.0180115562062,0.0118074789971}\"'],\n", " ['classoperation Main.MyModel.Tibia.Patella2.sRel \"Set Value\" --value=\"{0.0249230579594,0.115461116587,-0.00698883214514}\"',\n", " 'classoperation Main.MyModel.Patella.Lig.sRel \"Set Value\" --value=\"{0.00483548485266,-0.014494203176,0.00744105538977}\"']]" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "monte_carlo_macros[1:4]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let us expand the macro to also load and run the model. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "macro = [\n", " mc.Load('Knee.any'),\n", " mc.SetValue_random('Main.MyModel.Tibia.Patella2.sRel', patella_tendon_insertion ) ,\n", " mc.SetValue_random('Main.MyModel.Patella.Lig.sRel', patella_tendon_origin ) ,\n", " mc.RunOperation('Main.MyStudy.InverseDynamics'),\n", " mc.Export('Main.MyStudy.Output.Abscissa.t', \"time\"),\n", " mc.Export('Main.MyStudy.Output.MaxMuscleActivity', \"MaxMuscleActivity\")\n", "]\n", "mg = AnyMacro(macro, seed=1)\n", "monte_carlo_macros = mg.create_macros_MonteCarlo(100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Running the Monte Carlo macro" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The macro is passed to the AnyPyProcess object which excutes the macros" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "1cb74f95431d49bcab3f5a4bc0ddcf49", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Completed: \u001b[1;36m100\u001b[0m\n" ] }, { "data": { "text/html": [ "
\n" ], "text/plain": [] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from anypytools import AnyPyProcess \n", "\n", "app = AnyPyProcess()\n", "monte_carlo_results = app.start_macro( monte_carlo_macros )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The output object (`monte_carlo_result`) is a list-like object where each element is a dictionary with the output of the corresponding simulation.\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Length : 100\n", "Data keys in first element: ['time', 'MaxMuscleActivity', 'task_macro_hash', 'task_id', 'task_work_dir', 'task_name', 'task_processtime', 'task_macro', 'task_logfile']\n" ] } ], "source": [ "print('Length :', len(monte_carlo_results) )\n", "print('Data keys in first element: ', list(monte_carlo_results[0].keys()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Filtering results with Errors\n", "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. \n", "\n", "That is easily done. Results from simuations with errors will contain a an 'ERROR' key. We can use that to filter the results." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "monte_carlo_results[:] = [\n", " output \n", " for output in monte_carlo_results \n", " if 'ERROR' not in output\n", "]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Extracting data from the resutls\n", "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. \n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.00890538, 0.00927552, 0.00986515, ..., 0.00986515, 0.00927552,\n", " 0.00890538],\n", " [0.00761831, 0.00793615, 0.00844398, ..., 0.00844398, 0.00793614,\n", " 0.00761831],\n", " [0.00965061, 0.01004897, 0.0106812 , ..., 0.0106812 , 0.01004897,\n", " 0.00965062],\n", " ...,\n", " [0.01007754, 0.0104962 , 0.01116247, ..., 0.01116247, 0.01049619,\n", " 0.01007754],\n", " [0.01026166, 0.01068757, 0.01136452, ..., 0.01136452, 0.01068757,\n", " 0.01026167],\n", " [0.00981378, 0.01021894, 0.01086207, ..., 0.01086207, 0.01021893,\n", " 0.00981379]], shape=(100, 100))" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "monte_carlo_results['MaxMuscleActivity']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the results\n", "Finally we can plot the result of the Monte Carlo study" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[Processing tasks ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100/100 0:01:07 0:00:00\n\n", "text/plain": "Processing tasks \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m 100/100 \u001b[33m0:01:07\u001b[0m \u001b[36m0:00:00\u001b[0m\n" }, "metadata": {}, "output_type": "display_data" } ] } }, "5b2de447f3ec43a29cfed31968dbbdb3": { "model_module": "@jupyter-widgets/output", "model_module_version": "1.0.0", "model_name": "OutputModel", "state": { "layout": "IPY_MODEL_f5ee03c5e3a7486bb81231019293ee28", "outputs": [ { "data": { "text/html": "
Processing tasks ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 25/25 0:00:22 0:00:00\n\n", "text/plain": "Processing tasks \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m 25/25 \u001b[33m0:00:22\u001b[0m \u001b[36m0:00:00\u001b[0m\n" }, "metadata": {}, "output_type": "display_data" } ] } }, "f5ee03c5e3a7486bb81231019293ee28": { "model_module": "@jupyter-widgets/base", "model_module_version": "2.0.0", "model_name": "LayoutModel", "state": {} } }, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }