diff --git a/Homework4/README.md b/Homework4/README.md index 2ea7b2f..e0379ce 100644 --- a/Homework4/README.md +++ b/Homework4/README.md @@ -1,60 +1,65 @@ # SP4E - Homework 4 ---- ## General Info This file provides a brief documentation and information related to the fourth (last) Homework of the course "Scientific Programming for Engineers", fall 2019. This homework is done by **O. Ashtari** and **A. Sieber**. Last update: 15.01.2020 ---- ## Project Description +The first part of the Homework aims at using the `Pybind11` library in order to create Python binding for the C++ particles code. The ultimate goal of this part consists in allowing the user to run this code through a Python interface. +The second part of the Homework focuses on the planet trajectory simulation module of the particles code. There, a Python optimization routine is written which main purpose is to modify the initial planets states (one planet at a time) in order to minimize the error of their computed trajectories when compared to reference trajectories. ---- ## Executable Files ---- ## Optimization routine The optimization routine `error_minimization.py` is intended to minimize the error on Mercury trajectory (or any other planet of the solar system) by scaling the initial velocity of the planet of interest. The minimization is performed by computing the error between the computed planet trajectory and a reference planet trajectory. This routine calls `main.py`, therefore, in order to run it, the user must first make sure she/he can launch the particles code through the python interface. Once the requirement is satisfied the `error_minimization.py` routine can be called as follow: ``` -$ python3 directory_comp directory_ref planet input_file scale nb_steps freq +$ python3 directory_comp directory_ref planet input_file scale nb_steps freq ``` where: * `directory_comp`, the path to the directory containing the computed trajectory data. * `directory_ref`, the path to the directory containing the reference trajectory data. * `planet`, the planet of interest. * `input_file`, the file containing the initial state of the different planets. * `scale`, the initial scaling factor (initial guess of the minimization). * `nb_steps`, the number of time steps of the simulation within the particles code. * `freq`, the dumping frequency at which the particles code writes outputs. -As an example a minimization on Mercury trajectory over 365 day with a dumping frequency of 1 day could be run as such: +As an example a minimization on Mercury trajectory over 365 days with a dumping frequency of 1 day could be run as such: ``` -$ python3 dumps trajectories mercury init.csv 1.0 365 1 +$ python3 dumps trajectories mercury init.csv 1.0 365 1 ``` -The `error_minimization.py` routine the prints the scaling factorminimizing the error, the value of this error and the amount of minimization iterations needed to reach an optimum. It moreover plots the evolution of the error versus the scaling factor. +The `error_minimization.py` routine then prints the scaling factor minimizing the error, the value of this error and the amount of minimization iterations needed to reach an optimum. It moreover plots the evolution of the error versus the scaling factor as shown in the link below: +![Optimization on Mercury initial velocity](/diffusion/9482/browse/master/Homework4/optimization_example/minimization_evolution.png) + +The above graph shows that the optimized scaling is evaluated at 0.3989. Comparing this result to the exact scaling of 0.4 (found by dividing the reference initial velocity by the one stored in `init.csv`) leads to a relative error of less than 0.3% in the optimization procedure. ---- ## Comment on how createSimulation function is overloaded To comment on this function overload, let's first discuss what is the role of `createComputes` in the whole code. In the `ParticlesFactoryInterface` class, `createComputes` is a Real to voind function which is to be defined later, if needed. In the constructor of `MaterialPointsFactory`, `createComputes` is defined to be the default function (i.e. to be `createDefaultComputes`.) In the default function, `ComputeTemperature` is added to the system evolution object. Similarly in the `PlanetsFactory` where `createComputes` is defined to be `createDefaultComputes` in which `verlet` is added to the system evolution object. -As a summary, by constructing an object of type `MaterialPointsFactory` or `PlanetsFactory`, `createComputes` takes a default definition. It is then used inside the first definition of `createSimulation` method whose arguments are file name and time step size. The whole idea behind overloading `createSimulation` and using the templated functor is to give the flexibility to the python user to tailor compute cobjects (i.e. `ComputeTemperature` or `verlet` etc.) based on his/her needs more easily. The reason is that when a C++ code is wrapped to be used through python, the C++ side is supposed to not be touched anymore. However, default settings for material point, for instance, is hard-coded inside its factory. Therefore, to change default values, the python user should manipulate the C++ files. Or assume, like the previous homework, one wants to get length, conductivity, etc. from command line arguments. Or one wants to develop the code such that these variables are read from a file. Using the overloaded `createSimulation` with very little insight into how the C++ code is designed, one can do all these by just writing few lines in python, with neither changing any C++ file nor editing the binding file. +As a summary, by constructing an object of type `MaterialPointsFactory` or `PlanetsFactory`, `createComputes` takes a default definition. It is then used inside the first definition of `createSimulation` method whose arguments are file name and time step size. The whole idea behind overloading `createSimulation` and using the templated functor is to give the flexibility to the python user to tailor compute cobjects (i.e. `ComputeTemperature` or `verlet` etc.) based on his/her needs more easily. The reason is that when a C++ code is wrapped to be used through python, the C++ side is supposed to not be touched anymore. However, default settings for material point, for instance, is hard-coded inside its factory. Therefore, to change default values, the python user should manipulate the C++ files. Or assume, like the previous homework, one wants to get length, conductivity, etc. from command line arguments. Or one wants to develop the code such that these variables are read from a file. Using the overloaded `createSimulation` with very little insight into how the C++ code is designed, one can do all these by just writing few lines in python, with neither changing any C++ file nor editing the binding file. diff --git a/Homework4/optimization_example/minimization_evolution.png b/Homework4/optimization_example/minimization_evolution.png new file mode 100644 index 0000000..e7d8932 Binary files /dev/null and b/Homework4/optimization_example/minimization_evolution.png differ diff --git a/Homework4/src/error_minimization.py b/Homework4/src/error_minimization.py index 9b2f7e5..a624933 100644 --- a/Homework4/src/error_minimization.py +++ b/Homework4/src/error_minimization.py @@ -1,135 +1,136 @@ import numpy as np import argparse import matplotlib.pyplot as plt import glob from scipy import optimize import main # Python function ################# # Reads the position of a given planet over time and stores it in an array def readPosition(planet_name, directory): file_list = sorted(glob.glob(directory + '/*.csv')) m = len(file_list) n = 3 positions = np.zeros([m,n]) i = 0 for file_path in file_list: data = np.genfromtxt(file_path, delimiter = ' ', dtype=np.str) loc = np.where(data == planet_name) line = loc[0] positions[i, 0] = float(data[line, 0]) positions[i, 1] = float(data[line, 1]) positions[i, 2] = float(data[line, 2]) i = i + 1 return positions # Computes the l2-norm error between the reference and computed positions of a given planet def computeError(positions, positions_ref): e = 0.0 for i in range(len(positions)): diff = (positions_ref[i, 0] - positions[i, 0])**2 + (positions_ref[i, 1] - positions[i, 1])**2 + (positions_ref[i, 2] - positions[i, 2])**2 e = e + diff error = np.sqrt(e) return error # Modifies the input file by scaling the velocity components of a given planet def generateInput(scale, planet_name, input_filename, output_filename): data = np.genfromtxt(input_filename, delimiter = ' ', dtype=np.str) loc = np.where(data == planet_name) line = loc[0] data[line, 3] = float(data[line, 3])*scale data[line, 4] = float(data[line, 4])*scale data[line, 5] = float(data[line, 5])*scale np.savetxt(output_filename, data, delimiter=' ', fmt="%s") return 0 # Calls main.py to run the particle code based on a specified input file def launchParticles(input_file, nb_steps, freq): particle_type = 'planet' timestep = 1 main.main(nb_steps, freq, input_file, particle_type, timestep) return 0 # Runs the particle code and computes the l2-norm error based on a specified input file def runAndComputeError(scale, *args): planet_name = args[0] input_file = args[1] nb_steps = args[2] freq = args[3] generateInput(scale, planet_name, input_file, 'scaled_input.csv') launchParticles('scaled_input.csv', nb_steps, freq) positions = readPosition(planet_name, dir_comp) positions_ref = readPosition(planet_name, dir_ref) error = computeError(positions, positions_ref) return error # Arguments parser ################## parser = argparse.ArgumentParser(description='Run particles code and computes simulation error') parser.add_argument('directory_comp', type=str, help='path to the directory of the computed data') parser.add_argument('directory_ref', type=str, help='path to the directory of the reference data') parser.add_argument('planet', type=str, help='Name of the planet to be investigated') parser.add_argument('input_file', type=str, help='Name of the inputs file') parser.add_argument('scale', type=float, help='Initial velocity scaling factor') parser.add_argument('nb_steps', type=int, help='Number of time steps') parser.add_argument('freq', type=int, help='Outputs dumping frequency') args = parser.parse_args() dir_comp = args.directory_comp dir_ref = args.directory_ref planet_name = args.planet input_file = args.input_file scale = args.scale nb_steps = args.nb_steps freq = args.freq # Minimization ############## # Stores scaling factor and error value after each minimization iteration def reporter(scale_r): global scale_int global error_int scale_int = np.vstack([scale_int, scale_r]) error_r = runAndComputeError(scale_r, planet_name, input_file, nb_steps, freq) error_int = np.vstack([error_int, error_r]) return 0 scale_0 = scale scale_int = scale_0 error_int = runAndComputeError(scale_0, planet_name, input_file, nb_steps, freq) solution = optimize.fmin(runAndComputeError, scale_0, args = (planet_name, input_file, nb_steps, freq), callback = reporter, xtol = 0.0001, ftol = 0.0001, maxiter = 100, full_output = True) +print('Optimized scaling factor for %s trajectory:' % planet_name, *scale_int[-1]) # Post-Processing ################# fig = plt.figure(1) ax = fig.add_subplot(111) ax.plot(scale_int, error_int, 'r--o') # Axis labels ax.set_xlabel('scaling', fontsize = 14) ax.set_ylabel('error: L-2 norm', fontsize = 14) # Title of the window ax.set_title('%s trajectory error minimization' % planet_name, fontsize = 14) plt.savefig('minimization_evolution.png') plt.show() \ No newline at end of file