Revised Monte Carlo Tutorial Page (#1511)

* Adding new monte carlo tutorial images

* Delete Trick-QP-Random.png

* Delete Trick-QP-Calculated.png

* Delete Trick-QP.png

* Delete Trick-QP-File.png

* Delete Trick-DRE.png

* Delete Trick-DP.png

* Added Trick-DP image

* Broad Tutorial Modifications

Lots of the content was changed. The navigation and headers are not correct yet.

* Fixed table of contents

* Fixed toc again
This commit is contained in:
Adam M Garsha 2023-06-07 11:10:59 -05:00 committed by GitHub
parent 506bac5de7
commit e102dd7ef2
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
10 changed files with 344 additions and 454 deletions

View File

@ -1,454 +1,344 @@
| [Home](/trick) → [Tutorial Home](Tutorial) → Monte Carlo | | [Home](/trick) → [Tutorial Home](Tutorial) → Monte Carlo |
|--------------------------------------------------------| |--------------------------------------------------------|
Monte Carlo is an advanced simulation capability provided by Trick that allows users to repeatedly run copies of a simulation with different input values. Users can vary the input space of a simulation via input file, random value generation, or by calculating values from previous Monte Carlo runs in a process called optimization. This tutorial will show you how to develop simulations that can take advantage of this capability. # Monte Carlo
**For a thorough explanation of Monte Carlo and its features, read the [Monte Carlo User Guide](/trick/documentation/simulation_capabilities/UserGuide-Monte-Carlo).** **Contents**
## Core Simulation * [What is Monte Carlo?](#what-is-monte-carlo)
The core simulation this tutorial is built around is a slight modification of the Analytical Cannon simulation created earlier in the tutorial. **The primary modification made was to change the *init_angle* attribute from radians to degrees.** * [Example Task](#example-task)
* [Input Files](#input-files)
### cannon.h - [Listing - **angle_value_list**](#listing-value-list)
This header file contains the cannon structure we are going to use to create our simulation object and function prototypes we will use to specify our Trick jobs. - [Listing - **input.py**](#listing-input_1)
* [Random Input Generation](#random-input-generation)
```C - [Listing - **input.py**](#listing-input-2)
/* * [Optimization](#optimization)
PURPOSE: Cannon structure and simulation functions. - [Listing - **optimization.h**](#listing-optimization-h)
*/ - [Listing - **optimization.c**](#listing-optimization-c)
#ifndef CANNON_H - [Listing - **input.py**](#listing-input-3)
#define CANNON_H - [Listing - **S_Define**](#listing-s-define)
#include "trick_utils/comm/include/tc.h"
#include "trick_utils/comm/include/tc_proto.h" ***
typedef struct <a id=what-is-monte-carlo></a>
{ ## What is Monte Carlo?
double vel0[2]; // *i m Initial velocity of the cannonball.
double pos0[2]; // *i m Initial position of the cannonball. Monte Carlo is an advanced simulation capability provided by Trick that allows users to repeatedly run copies of a simulation with different input values. Users can vary the input space of a simulation via input file, random value generation, or by calculating values from previous Monte Carlo runs in a process called optimization. This tutorial will show you how to modify the cannon_numeric simulation to take advantage of this capability.
double init_speed; // *i m/s Initial barrel speed.
double init_angle; // *i ° Launch angle of cannon. **For a thorough explanation of Monte Carlo and its features, read the [Monte Carlo User Guide](/trick/documentation/simulation_capabilities/UserGuide-Monte-Carlo).**
double acc[2]; // m/s2 xy-acceleration <a id=example-task></a>
double vel[2]; // m/s xy-velocity ## Example Task
double pos[2]; // m xy-position **What would be the optimal launch angle required to ensure our cannonball travels the furthest distance?** Let us assume that we have no conception of physics or trigonometry and that we don't already know the answer.
double time; // s Model time. <p align="center">
double timeRate; // -- Model time per sim time. <img src="images/OptimalLaunchAngle.png" width=550px/>
</p>
int impact; // -- Has impact occured?
double impactTime; // s Time of impact. <a id=input-files></a>
## Input Files
} CANNON; Input files allow you to specify the exact values you want on a particular simulation run. Input files are the most precise implementation, but they require more effort to setup and modify later down the road. Input files can contain multiple (tab or space) delimited columns filled with numerical information.
#ifdef __cplusplus ### Value List
extern "C" { Create the following text file in your simulation directory with the name **angle\_value\_list**:
#endif
<a id=listing-value-list></a>
// Function prototypes. **Listing - angle_value_list**
int cannon_default_data(CANNON*);
int cannon_init(CANNON*); ```
int cannon_analytic(CANNON*); 0.1
0.2
#ifdef __cplusplus 0.3
} 0.4
#endif 0.5
#endif // #ifndef CANNON_H 0.6
``` 0.7
0.8
### cannon.c 0.9
This .c file contains the logic used to advance our simulation as well as the default data and initial data we wish to with our cannon object. 1
```C 1.1
/* 1.2
PURPOSE: Cannon model simulation jobs and logic. 1.3
*/ 1.4
1.5
#include <math.h>
#include "../include/cannon.h" ```
#include "trick/exec_proto.h" This text file will be used to assign the cannon's initial angle. Remember that this angle is in radians.
// Default data job. ### Updating our Input File
int cannon_default_data( CANNON* C ) The simulation input file must be adjusted to run the simulation using Monte Carlo techniques.
{
C->acc[0] = 0.0; ```
C->acc[1] = -9.81; % cd $HOME/trick_sims/SIM_cannon_analytic/RUN_test
% vi input.py
C->init_angle = 15; ```
C->init_speed = 50.0;
The simulation must be told to enable Monte Carlo, recognize a Monte Carlo variable, and use the angle_value_list file created above. To accomplish this, change the input file to include the following:
C->pos0[0] = 0.0;
C->pos0[1] = 0.0; <a id=listing-input_1></a>
**Listing - input.py**
C->time = 0.0;
```python
C->impact = 0;
C->impactTime = 0.0; exec(open("Modified_data/cannon.dr").read())
return 0; # Enable Monte Carlo.
} trick.mc_set_enabled(1)
// Initialization job. # Sets the number of runs to perform to 15. Trick will not exceed the number of values in an input file.
int cannon_init( CANNON* C) trick.mc_set_num_runs(15)
{
C->vel0[0] = C->init_speed*cos(C->init_angle * .0174533); # Create and add a new Monte Carlo File variable to the simulation.
C->vel0[1] = C->init_speed*sin(C->init_angle * .0174533); mcvar_launch_angle = trick.MonteVarFile("dyn.cannon.init_angle", "angle_value_list", 1, "rad")
trick.mc_add_variable(mcvar_launch_angle)
C->vel[0] = C->vel0[0];
C->vel[1] = C->vel0[1]; # Stop Monte Carlo runs after 25 seconds of simulation time
trick.stop(25)
C->impactTime = 0.0;
C->impact = 0.0; # Stop Monte Carlo runs if they take longer than 1 second of real time
trick.mc_set_timeout(1)
return 0; ```
}
After the file has been adjusted, save it and run the simulation.
// Model simulation logic.
int cannon_analytic( CANNON* C ) ```
{ % cd ..
C->acc[0] = 0.00; % ./S_main*.exe RUN_test/input.py
C->acc[1] = -9.81; ```
C->vel[0] = C->vel0[0] + C->acc[0] * C->time; The terminal will display a fairly verbose Monte Carlo process detailing each run. After completion, information about each run will be put into a MONTE_RUN_test directory.
C->vel[1] = C->vel0[1] + C->acc[1] * C->time;
In order to complete our task of finding the optimal launch angle, it will be necessary to plot each run to compare the distance achieved by each cannon. Open up the Trick Data Product application with the following command.
C->pos[0] = C->pos0[0] + (C->vel0[0] + (0.5) * C->acc[0] * C->time) * C->time;
C->pos[1] = C->pos0[1] + (C->vel0[1] + (0.5) * C->acc[1] * C->time) * C->time; ```
if (C->pos[1] < 0.0) trick-dp &
{ ```
C->impactTime = (- C->vel0[1] - sqrt( C->vel0[1] * C->vel0[1] - 2 * C->pos0[1]))/C->acc[1];
C->pos[0] = C->impactTime * C->vel0[0]; <p align="center">
C->pos[1] = 0.0; <img src="images/Trick-DP.png" width=750px/>
C->vel[0] = 0.0; </p>
C->vel[1] = 0.0;
if (!C->impact) Right click the MONTE_RUN_test directory and select **Add run(s)**. Then open quick plot by clicking the blue lightning bolt. Expand the dyn.cannon.pos[0-1] variable in the left pane and create a curve with pos[1] as the Y axis and pos[0] as the X axis. Finally, click the **comparison plot** button in the actions menu.
{
C->impact = 1; <p align="center">
fprintf(stderr, "\n\nIMPACT: ModelTime = %.9f, pos = %.9f\n\n", C->impactTime, C->pos[0]); <img src="images/MONTE_list_plot.png" width=750px/>
} </p>
}
The various curves show the trajectories of each cannon run. It may be necessary to hide the legend if all the run names cover up the plot.
C->time += 0.01 ;
return 0 ; <a id=random-input-generation></a>
} ## Random Input Generation
``` Random Input Generation provides users with the ability to statistically generate input values along a Gaussian or Poisson distribution. Random generation is less precise than an input file, but it is more extensible and much easier to modify. Modify the input file again to use a gaussian distribution to generate launch angles.
### Simulation Definition <a id=listing-input-2></a>
The simulation definition file simply sets up our simulation objects and configures their respective Trick jobs. **Listing - input.py**
```C++
/* ```python
PURPOSE: Monte tutorial simulation definition. exec(open("Modified_data/cannon.dr").read())
LIBRARY DEPENDENCIES: ((cannon/src/cannon.c))
*/ # Enable Monte Carlo.
trick.mc_set_enabled(1)
#include "sim_objects/default_trick_sys.sm"
##include "cannon/include/cannon.h" # Run 100 randomly generated variables.
trick.mc_set_num_runs(100)
class CannonSimObject : public Trick::SimObject
{ # Create a Monte Carlo Random variable.
public: mcvar_launch_angle = trick.MonteVarRandom("dyn.cannon.init_angle", trick.MonteVarRandom.GAUSSIAN, "rad")
CANNON cannon;
# Set the random number generator seed.
CannonSimObject() mcvar_launch_angle.set_seed(1)
{
("default_data") cannon_default_data(&cannon); # Set the standard deviation for this bellcurve.
("initialization") cannon_init(&cannon); mcvar_launch_angle.set_sigma(30)
(0.01, "scheduled") cannon_analytic(&cannon);
} # Set the center of the bellcurve.
}; mcvar_launch_angle.set_mu(3.141592 / 4) # PI/4
CannonSimObject cso; # Set the maximum and minimum values to be generated.
``` mcvar_launch_angle.set_max(3.141592 / 2) # PI/2
mcvar_launch_angle.set_min(3.141592 / 12) #PI/12
### Data Recording
Visualizing the differences between the various Monte Carlo methods requires us to first establish a data recording file. In order to do this, the simulation must first be compiled; run the following command in your shell and compile the simulation above. # The min and max are absolute values, not relative to mu.
mcvar_launch_angle.set_min_is_relative(False)
``` mcvar_launch_angle.set_max_is_relative(False)
trick-CP
``` # Add the variable.
trick.mc_add_variable(mcvar_launch_angle)
After successfully compiling the simulation, open the Trick Data Recording Editor with the following
shell command: # Stop Monte Carlo runs after 25 seconds of simulation time
trick.stop(25)
```
trick-dre & # Stop Monte Carlo runs if they take longer than 1 second of real time
``` trick.mc_set_timeout(1)
```
<p align="center">
<img src=images/Trick-DRE.png width="750"/> Run the script and plot the curves the same way as before. You will end up with something similar to this:
</p>
<p align="center">
#### Steps <img src="images/MONTE_gauss_plot.png" width=750px/>
01. In the left pane will be *cso*, the CannonSimObject we created. Expand this tree and double click on the pos[2] variable. Ensure the variable appears in the "Selected Variables" section at the bottom. </p>
01. At the top, change DR Cycle from 0.1 to 0.01.
01. Save the data recording file in the simulation directory. <a id=optimization></a>
## Optimization
### Test Run Optimization is the process of evaluating the previous run's data for use in the next run. In essence, you are optimizing each subsequent run and closing in on a specific value; in this instance, we are closing in on the optimal launch angle.
Create a sub-directory called *RUN_test* in your simulation directory. In this new directory create an input file named *test.py*. This input file executes the data recording file you saved above and stops the simulation after 10 seconds of simulation time.
<a id=listing-optimization-h></a>
```python **Listing - optimization.h**
exec(open("monte_cannon.dr").read())
trick.stop(10) ### optimization.h
``` We need to create two new Trick jobs to house our optimization logic. Create this file in your include directory.
After running the input file, open up the Trick Data Product application with the following command. ```C
/******************************* TRICK HEADER ****************************
``` PURPOSE: Function prototypes for Monte Carlo optimization.
trick-dp & *************************************************************************/
```
#ifndef OPTIMIZATION_H
<p align="center"> #define OPTIMIZATION_H
<img src="images/Trick-DP.png" width=750px/>
</p> #include "../include/cannon.h"
Find the RUN_test in the Runs Tree panel and add it to your run selections. Click the blue lightning bolt to open Trick Quick Plot. #ifdef __cplusplus
extern "C" {
#endif
<p align="center">
<img src="images/Trick-QP.png" width=750px/> int cannon_slave_post(CANNON *);
</p> int cannon_master_post(CANNON *);
Expand the cso.cannon.pos[0-1] variable in the left pane and create a curve with pos[1] as the Y axis and pos[0] as the X axis. Once your DP Content looks like the above image, click the comparison plot button in the actions menu. #ifdef __cplusplus
}
<p align="center"> #endif
<img src="images/RUN_test_plot.png" width=750px/> #endif
</p> ```
In each of the bellow sections, you can repeat these steps to visualize the differences between Monte Carlo implementations. ### optimization.c
What we are doing in these two functions is sending the slave's cannon structure from after the run has completed back to the master. The master then analyzes the data and sends the new run information to the slave. This cycles over and over again until we hit the number of runs specified in our input script. Create this file in your src directory.
## The Task
**What would be the optimal launch angle required to ensure our cannonball travels the furthest distance?** Let us assume that we have no conception of physics or trigonometry and that we don't already know the answer. <a id=listing-optimization-c></a>
**Listing - optimization.c**
<p align="center">
<img src="images/OptimalLaunchAngle.png" width=550px/> ```C
</p> /******************************* TRICK HEADER ****************************
PURPOSE: Monte Carlo optimization functions.
## Input Files *************************************************************************/
Input files allow you to specify the exact values you want on a particular simulation run. Input files are the most precise implementation, but they require more effort to setup and modify later down the road. Input files can contain multiple (tab or space) delimited columns filled with numerical information.
#include "../include/optimization.h"
### Value List #include "../include/cannon.h"
Create the following text file in your simulation directory with the name **angle\_value\_list**: #include "sim_services/MonteCarlo/include/montecarlo_c_intf.h"
``` int cannon_slave_post(CANNON *C)
5 {
10 mc_write((char*) C, sizeof(CANNON));
15 return 0;
20 }
25
30 int cannon_master_post(CANNON *C)
35 {
40 CANNON run_cannon;
45 static double previous_distance = 0;
50 static double increment = 0.2; // Remember radians
55
60 // Get the run's cannon back from slave.
65 mc_read((char*) &run_cannon, sizeof(CANNON));
70
75 // Optimization logic.
80 if(run_cannon.pos[0] < previous_distance)
85 {
// Cut the increment in half and reverse the direction.
``` increment /= 2;
increment *= -1;
### Script }
Create a new directory called RUN_file and place the following python script in it with the name **file.py**:
C->init_angle += increment;
```python previous_distance = run_cannon.pos[0];
# -*- coding: UTF-8 -*- return 0;
}
exec(open("monte_cannon.dr").read()) ```
# Enable Monte Carlo. ### Modifications to input.py
trick.mc_set_enabled(1)
<a id=listing-input-3></a>
# Sets the number of runs to perform to 20. Trick will not exceed the number of values in an input file. **Listing - input.py**
trick.mc_set_num_runs(20)
```python
# Create and add a new Monte Carlo File variable to the simulation. exec(open("Modified_data/cannon.dr").read())
mcvar_launch_angle = trick.MonteVarFile("cso.cannon.init_angle", "angle_value_list", 1, "°")
trick_mc.mc.add_variable(mcvar_launch_angle) # Enable Monte Carlo.
trick.mc_set_enabled(1)
# Stop the simulation run after 15 seconds of simulation time.
trick.stop(15) # Run 25 optimizations.
# The more runs, the more precise your variable will end up assuming you wrote your optimization logic correctly.
``` trick.mc_set_num_runs(25)
Run the simulation with file.py as your input file and plot the runs like you did earlier. You will see something similar to the following: # Create a calculated variable and add it to Monte Carlo.
mcvar_launch_angle = trick.MonteVarCalculated("dyn.cannon.init_angle", "rad")
<p align="center"> trick.mc_add_variable(mcvar_launch_angle)
<img src="images/Trick-QP-File.png" alt="Trick-QP-File"/>
</p> # Stop Monte Carlo runs after 25 seconds of simulation time
trick.stop(25)
## Random Input Generation
Random Input Generation provides users with the ability to statistically generate input values along a Gaussian or Poisson distribution. Random generation is less precise than an input file, but it is more extensible and much easier to modify. # Stop Monte Carlo runs if they take longer than 1 second of real time
trick.mc_set_timeout(1)
### Script ```
```python
# -*- coding: UTF-8 -*- ### Modifications to S_Define
exec(open("data/monte_cannon.dr").read()) The last thing that we need to do is modify our simulation definition file and add the two new Trick jobs. As you can see, we have added a new library dependency, a new ## inclusion, and two new constructor jobs.
# Enable Monte Carlo. <a id=listing-s-define></a>
trick.mc_set_enabled(1) **Listing - S_Define**
# Run 100 randomly generated variables. ```C++
trick.mc_set_num_runs(100)
/************************TRICK HEADER*************************
# Create a Monte Carlo Random variable. PURPOSE:
mcvar_launch_angle = trick.MonteVarRandom("cso.cannon.init_angle", trick.MonteVarRandom.GAUSSIAN, "°") (S_define file for SIM_cannon_numeric)
LIBRARY DEPENDENCIES:
# Set the random number generator seed. (
mcvar_launch_angle.set_seed(1) (cannon/src/cannon_init.c)
(cannon/src/cannon_numeric.c)
# Set the standard deviation for this bellcurve. (cannon/src/cannon_shutdown.c)
mcvar_launch_angle.set_sigma(30) (cannon/src/optimization.c)
)
# Set the center of the bellcurve. *************************************************************/
mcvar_launch_angle.set_mu(45)
#include "sim_objects/default_trick_sys.sm"
# Set the maximum and minimum values to be generated. ##include "cannon/include/cannon_numeric.h"
mcvar_launch_angle.set_max(90) ##include "cannon/include/optimization.h"
mcvar_launch_angle.set_min(0)
class CannonSimObject : public Trick::SimObject {
# The min and max are absolute values, not relative to mu.
mcvar_launch_angle.set_min_is_relative(False) public:
mcvar_launch_angle.set_max_is_relative(False) CANNON cannon;
# Add the variable. CannonSimObject() {
trick_mc.mc.add_variable(mcvar_launch_angle) ("default_data") cannon_default_data( &cannon ) ;
("initialization") cannon_init( &cannon ) ;
# Stop the run after 15 simulation seconds. ("derivative") cannon_deriv( &cannon ) ;
trick.stop(15) ("integration") trick_ret= cannon_integ( & cannon ) ;
``` ("dynamic_event") cannon_impact( &cannon ) ;
("monte_slave_post") cannon_slave_post( &cannon ) ;
Run the script and plot the curves. You will end up with something similar to this: ("monte_master_post") cannon_master_post( &cannon ) ;
("shutdown") cannon_shutdown( &cannon ) ;
<p align="center"> }
<img src="images/Trick-QP-Random.png" alt="Trick-QP-Random.png"/> } ;
</p>
CannonSimObject dyn ;
## Optimization IntegLoop dyn_integloop (0.01) dyn ;
Optimization is the process of evaluating the previous run's data for use in the next run. In essence, you are optimizing each subsequent run and closing in on a specific value; in this instance, we are closing in on the optimal launch angle. void create_connections() {
dyn_integloop.getIntegrator(Runge_Kutta_4, 4);
### optimization.h }
We need to create two new Trick jobs to house our optimization logic. Create this file in your include directory. ```
```C Run the script and plot the curves.
/*
PURPOSE: Function prototypes for Monte Carlo optimization. ```
*/ % trick-CP
...
#ifndef OPTIMIZATION_H % ./S_main*.exe RUN_test/input.py
#define OPTIMIZATION_H ```
#include "../include/cannon.h" Once the simulation is complete, the x-y position plot should look like this:
#ifdef __cplusplus <p align="center">
extern "C" { <img src="images/MONTE_calculated_plot.png" alt="Trick-QP-Calculated"/>
#endif </p>
int cannon_slave_post(CANNON *);
int cannon_master_post(CANNON *);
#ifdef __cplusplus
}
#endif
#endif
```
### optimization.c
What we are doing in these two functions is sending the slave's cannon structure from after the run has completed back to the master. The master then analyzes the data and sends the new run information to the slave. This cycles over and over again until we hit the number of runs specified in our input script. Create this file in your src directory.
```C
/*
PURPOSE: Monte Carlo optimization functions.
*/
#include "../include/optimization.h"
#include "../include/cannon.h"
#include "sim_services/MonteCarlo/include/montecarlo_c_intf.h"
int cannon_slave_post(CANNON *C)
{
mc_write((char*) C, sizeof(CANNON));
return 0;
}
int cannon_master_post(CANNON *C)
{
CANNON run_cannon;
static double previous_distance = 0;
static double increment = 10;
// Get the run's cannon back from slave.
mc_read((char*) &run_cannon, sizeof(CANNON));
// Optimization logic.
if(run_cannon.pos[0] < previous_distance)
{
// Cut the increment in half and reverse the direction.
increment /= 2;
increment *= -1;
}
C->init_angle += increment;
previous_distance = run_cannon.pos[0];
return 0;
}
```
### Script
```python
# -*- coding: UTF-8 -*-
exec(open("data/monte_cannon.dr").read())
# Enable Monte Carlo.
trick.mc_set_enabled(1)
# Run 50 optimizations.
# The more runs, the more precise your variable will end up assuming you wrote your logic correctly.
trick.mc_set_num_runs(50)
# Create a calculated variable and add it to Monte Carlo.
mcvar_launch_angle = trick.MonteVarCalculated("cso.cannon.init_angle", "°")
trick_mc.mc.add_variable(mcvar_launch_angle)
# Stop the run after 15 seconds.
trick.stop(15)
```
### Simulation Definition
The last thing that we need to do is modify our simulation definition file and add the two new Trick jobs. As you can see, we have added a new library dependency, a new ## inclusion, and two new constructor jobs.
```C++
/*
PURPOSE: Monte tutorial simulation definition.
LIBRARY DEPENDENCIES: ((models/src/cannon.c) (models/src/optimization.c))
*/
#include "sim_objects/default_trick_sys.sm"
##include "models/include/cannon.h"
##include "models/include/optimization.h"
class CannonSimObject : public Trick::SimObject
{
public:
// Cannon object we are wrapping with CannonSimObject.
CANNON cannon;
// Constructor.
CannonSimObject()
{
("default_data") cannon_default_data(&cannon);
("initialization") cannon_init(&cannon);
(0.01, "scheduled") cannon_analytic(&cannon);
("monte_slave_post") cannon_slave_post(&cannon);
("monte_master_post") cannon_master_post(&cannon);
}
} ;
CannonSimObject cso;
```
Run the script and plot the curves. You will see this:
<p align="center">
<img src="images/Trick-QP-Calculated.png" alt="Trick-QP-Calculated"/>
</p>

Binary file not shown.

After

Width:  |  Height:  |  Size: 63 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 209 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 45 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 60 KiB

After

Width:  |  Height:  |  Size: 50 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 64 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 197 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 244 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 853 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 55 KiB