Page MenuHomec4science
Files F61126111

python_functions.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 17:06

python_functions.py
View Options

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"])
# 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()

Event Timeline

c4science · Help