diff --git a/hw4-pybind/README.md b/hw4-pybind/README.md index ca62769b..d44ff355 100644 --- a/hw4-pybind/README.md +++ b/hw4-pybind/README.md @@ -1,23 +1,30 @@ # Homework 4 -### Exercise 1: +Remember to include the folders `trajectories`, `pybind11` and `googletest`. +Also make a folder `build` which contains a folder named `dumps`. -1.2. : The role of function overloading "createSimulation" to associate one method or another one depending of the type or number of arguments that takes the method. +### Exercise 1: +1.2: The role of function overloading `createSimulation` to associate one method or another one depending of the type or number of arguments that takes the method. ### Exercise 2: -2.2. : To ensure that references are correctly managed to python bindings we used shared pointers : std::shared_ptr<>. Using this, the memory of these pointers will be DEALLOCATED when the python's reference goes to 0! +2.2: To ensure that references are correctly managed to python bindings we used shared pointers: `std::shared_ptr<>`. Using this, the memory of these pointers will be DEALLOCATED when the python's reference goes to 0! ### Exercise 4: -Can run the example using the command `./particles 365 1 ../init.csv planet 1` from the `build` folder. From C++ the gravitational constant G must be changed accordingly. +Can run the example using the command +`./particles 365 1 ../init.csv planet 1` +or +`python3 main.py 365 1 ../init.csv planet 1` +from the `build` folder. +The gravitational constant G has been accordingly changed in the c++ code to give same results as python. + +Note that we assume 5 integers in the filenames (e.g. `step-00000.csv`), i.e., not 4 as in the trajectory folder. +We can easily rename multiple files using the command `rename 's/step-0/step-00/g' *` in order to get same naming convention. The trajectories can be viewed using Paraview, and we observe trajectories as expected. Except for Mercury. ### Exercise 5 through 7: See the file `python_functions.py`. Run the file from the `build` folder with `python3 ../python_functions.py`. -Note that we assume 5 integers in the filenames (e.g. `step-00000.csv`), i.e., not 4 as in the trajectory folder. -We can easily rename multiple files using the command `rename 's/step-0/step-00/g' *` in order to get same naming convention. - -We find that the scaling factor for Mercury's velocity is approximately 0.4. +We find that the scaling factor for Mercury's velocity is approximately 0.39. diff --git a/hw4-pybind/python_functions.py b/hw4-pybind/python_functions.py index cb5d994d..319aa081 100644 --- a/hw4-pybind/python_functions.py +++ b/hw4-pybind/python_functions.py @@ -1,128 +1,128 @@ import numpy as np import subprocess import os from scipy import optimize import matplotlib.pyplot as plt # reads the trajectory of a given planet and makes a numpy array out of it # returns: 365 × 3 matrix, the columns being the three components of the planet position. def readPositions(planet_name, directory): out = np.zeros((365, 3)) # initialize out matrix for s in range(0, 365): # for each time step filename = directory + "/step-" + str(s).zfill(5) + ".csv" #print(filename) with open(filename) as file: for line in file: line = line.split() for word in line: #print(word) if word == planet_name: out[s, 0] = line[0] out[s, 1] = line[1] out[s, 2] = line[2] break; return out # computes L2 norm of the error in position def computeError(positions, positions_ref): return np.sqrt( np.sum( (positions-positions_ref)**2) ) # scales velocity of given planet in given input file def generateInput(scale, planet_name, input_filename, output_filename): v = np.zeros(3) # find and scale input velocity content = [] with open(input_filename) as file: for line in file: content.append(line) line = line.split() for word in line: if word == planet_name: v[0] = scale * float(line[3]) v[1] = scale * float(line[4]) v[2] = scale * float(line[5]) break; # break for loop of word search if planet found. Go to next line. # create new output file content new_content = [] for line in content: line = line.split() if planet_name in line: new_content.append( line[0] + " " + line[1] + " " + line[2] + " " + str(v[0]) + " " + str(v[1]) + " " + str(v[2]) + " " + " ".join(line[6:]) + "\n" ) else: new_content.append( " ".join(line) + "\n" ) # write content to file with open(output_filename, 'w') as new_file: new_file.writelines(new_content) # launches the particle code on an provided input def launchParticles(input, nb_steps, freq): - ### Calling c++ executable: - res = subprocess.call(['./particles', str(nb_steps), str(freq), input, "planet", "1"]) - ### Calling python function: - #res = subprocess.call(['python', "main.py", str(nb_steps), str(freq), input, "planet", "1"]) - + ### Alt 1: Calling c++ executable: + #res = subprocess.call(['./particles', str(nb_steps), str(freq), input, "planet", "1"]) + ### Alt 2: Calling python function: + res = subprocess.call(['python', "main.py", str(nb_steps), str(freq), input, "planet", "1"]) +# Comment to above code: Do not really need subprocess to call a python function from python as we could simply "import" it # gives the error for a given scaling velocity factor of a given planet def runAndComputeError(scale, planet_name, input, nb_steps, freq): scaled_input = input[:-4] + "_scaled.csv" #print(scaled_input) generateInput(scale, planet_name, input, scaled_input) launchParticles(scaled_input, nb_steps, freq) positions = readPositions(planet_name, "dumps") positions_ref = readPositions(planet_name, "../trajectories") error = computeError(positions, positions_ref) #print("runAndComputeError: The L2 error for planet ", planet, " is ", error) return error ################################### ####### Testing functions ######### ################################### ### Note: ### relative file paths are given here such that the script can be run from the build folder! ######### Exercise 5 - Question 1 ############ planet = "mercury" # positions = readPositions(planet, "dumps") # positions_ref = readPositions(planet, "../trajectories") ######### Exercise 5 - Question 2 ############# # error = computeError(positions, positions_ref) # print("The L2 error for planet ", planet, " is ", error) ######### Exercise 6 - Question 1 ############# #generateInput(2.0, planet, "../init.csv", "../init_scaled.csv") ######### Exercise 6 - Question 2 ############# #launchParticles("../init.csv", 365, 1) ######### Exercise 6 - Question 3 ############## runAndComputeError(1.0, planet, "../init.csv", 365, 1) ######### Exercise 7 - Question 1 ############## init_guess = 1.0 scaling_factors = [init_guess] errors = [] def callback(xk): scaling_factors.append(xk[0]) approximate_scaling = optimize.fmin(runAndComputeError, 1.0, args=(planet, "../init.csv", 365, 1), xtol=0.0001, ftol=0.0001, maxiter=100000, callback=callback) print("The approximate scaling of velocity is ", approximate_scaling[0]) print("scale \t error") for scale in scaling_factors: e = runAndComputeError(scale, planet, "../init.csv", 365, 1) print(scale, "\t", e) errors.append(e) plt.figure() plt.plot(scaling_factors, errors, 'k-*') plt.xlabel("Scaling factor") plt.ylabel("Error") plt.savefig("error.png") plt.show() plt.close()