Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62747431
error_minimization.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, May 15, 08:07
Size
4 KB
Mime Type
text/x-python
Expires
Fri, May 17, 08:07 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17683385
Attached To
R9482 SP4E_Homework_Ashtari_Sieber
error_minimization.py
View Options
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
)
# 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
()
Event Timeline
Log In to Comment