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