Page MenuHomec4science
Files F61170967

error_minimization.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 23:53

error_minimization.py
View Options

import numpy as np
import argparse
import matplotlib.pyplot as plt
import glob
from scipy import optimize
import main
# Functions
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
def computeError(positions, positions_ref):
e = 0.0
print(len(positions))
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
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
def launchParticles(input_file, nb_steps, freq):
particle_type = 'planet'
timestep = 1
main.main(nb_steps, freq, input_file, particle_type, timestep)
return 0
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)
print(error)
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
#positions = readPosition(planet_name, dir_comp)
#positions_ref = readPosition(planet_name, dir_ref)
#errors = computeError(positions, positions_ref)
#print(errors)
def reporter(scale_r):
global scale_int
scale_int = np.vstack([scale_int, scale_r])
return 0
scale_0 = scale
scale_int = scale_0
solution = optimize.fmin(runAndComputeError, scale_0, args = (planet_name, input_file, nb_steps, freq), callback = reporter, xtol = 0.0001, ftol = 0.0001, maxiter = 10)
print(scale_int)

Event Timeline

c4science · Help