Fixed timing offsets are added multiple times in estimateTimingAndVelocity function. #45
Description
I'm working on a fix for a problem in wmpl/Trajectory/Trajectory.py which causes fixed timing offsets to be added multiple times.
A recent bugfix added these lines which allowed the outputted offsets to be shown correctly:
# Reset fixed times to 0, as the timing offsets have already been applied
for station in self.fixed_time_offsets:
self.fixed_time_offsets[station] = 0.0
However, this does not address the underlying problem, as the reset to zero does not happen early enough in the code sequence. The offsets are first applied in the infillTrajectory() function:
# Add a fixed offset to time data if given
if str(station_id) in self.fixed_time_offsets:
time_data += self.fixed_time_offsets[str(station_id)]
It would be simple to reset the offsets at this point in the code, but that creates the problem that we can lo longer access the fixed offset data for output purposes. My proposed fix is to go around required places in the code and make a copy set of the original offset data, to allow it to be accessed in the output stages. Starting with the above changed to:
# Add a fixed offset to time data if given
if str(station_id) in self.fixed_time_offsets:
time_data += self.fixed_time_offsets[str(station_id)]
# Make a copy of the fixed offset values for output display purposes.
self.fixed_time_offsets_copy[str(station_id)] = self.fixed_time_offsets[str(station_id)]
# set the original value to zero as it has already been applied.
self.fixed_time_offsets[str(station_id)] = 0
Early in the def init stage of the class Trajectory(object), we have to also initialise the self.fixed_time_offsets_copy array and give it identical values to the self.fixed_time_offsets array:
# Estimating the difference in timing between stations, and the initial velocity if this flag is True
self.fixed_time_offsets = {}
self.fixed_time_offsets_copy = {}
if isinstance(estimate_timing_vel, str):
# If a list of fixed timing offsets was given, parse it into a dictionary
for entry in estimate_timing_vel.split(','):
station, offset = entry.split(":")
self.fixed_time_offsets[station] = float(offset)
self.fixed_time_offsets_copy[station] = float(offset)
print("Fixed timing given:", self.fixed_time_offsets)
self.estimate_timing_vel = False
elif isinstance(estimate_timing_vel, bool):
self.estimate_timing_vel = estimate_timing_vel
else:
self.estimate_timing_vel = True
# Extract the fixed times from the fixed time offsets
self.fixed_times = fixed_times
if isinstance(estimate_timing_vel, bool) and isinstance(self.fixed_times, str):
if estimate_timing_vel:
self.fixed_time_offsets = {}
for entry in self.fixed_times.split(','):
station, offset = entry.split(":")
self.fixed_time_offsets[station] = float(offset)
# Make a copy of the offsets for display output use.
self.fixed_time_offsets_copy[station] = self.fixed_time_offsets[station]
Next in estimateTimingAndVelocity, timing offset dictionaries are created, so we have to copy those:
# Run timing offset estimation if it needs to be done
if estimate_timing_vel:
# Make a dictionary of station IDs and the time offset estimation status
# - if the time is fixed, a number is given
# - if the time is to be estimated, True is set
self.stations_time_dict = collections.OrderedDict()
self.stations_time_dict_copy = collections.OrderedDict()
station_list = [str(obs.station_id) for obs in observations]
for obs in observations:
if str(obs.station_id) in self.fixed_time_offsets:
self.stations_time_dict[str(obs.station_id)] = self.fixed_time_offsets[str(obs.station_id)]
self.stations_time_dict_copy[str(obs.station_id)] = self.fixed_time_offsets_copy[str(obs.station_id)]
else:
self.stations_time_dict[str(obs.station_id)] = True
self.stations_time_dict_copy[str(obs.station_id)] = True
# If no fixed times are given, set the station with the longest track as the reference station
# (time difference = 0)
if len(self.fixed_time_offsets) == 0:
obs_points = [obs.kmeas for obs in observations]
ref_index = obs_points.index(max(obs_points))
self.stations_time_dict[station_list[ref_index]] = 0
self.stations_time_dict_copy[station_list[ref_index]] = 0
# Generate an initial guess for stations which have no fixed time
p0 = np.zeros(len([val for val in self.stations_time_dict.values() if val is True]))
Now we are set to use the copy array to display fixed offsets in the output where these occur:
# If the minimization was successful, apply the time corrections
if estimate_velocity:
stat_count = 0
for i, obs in enumerate(observations):
# Check the station timing dictionary to see if the station is fixed
stat_status = self.stations_time_dict[str(obs.station_id)]
# Make a copy of the time offset from the original values stored before they were zeroed.
stat_status_copy = self.stations_time_dict_copy[str(obs.station_id)]
# If the station has a fixed time offset, read it
if not isinstance(stat_status, bool):
t_diff = stat_status
t_diff_copy = stat_status_copy
if self.verbose:
print('STATION ' + str(obs.station_id) + ' TIME OFFSET = ' + str(t_diff_copy) + ' s (fixed offset applied)')
# Add the final time difference of the site to the list
time_diffs[i] = t_diff_copy
# Otherwise read the estimated offset
else:
t_diff = timing_mini.x[stat_count]
stat_count += 1
if self.verbose:
print('STATION ' + str(obs.station_id) + ' TIME OFFSET = ' + str(t_diff) + ' s')
# Add the final time difference of the site to the list
time_diffs[i] = t_diff
# Skip NaN and inf time offsets
if np.isnan(t_diff) or np.isinf(t_diff):
continue
# Apply the time shift to original time data
obs.time_data = obs.time_data + t_diff
# Apply the time shift to the excluded time
if obs.excluded_time is not None:
obs.excluded_time = [ex_time + t_diff for ex_time in obs.excluded_time]
# Add the final time difference of the site to the list
# time_diffs[i] = t_diff_copy
Note that these last two lines are commented out as earler the addition of time_diffs[i] is divided into two cases where they are either the solver's own iterated values or the preset fixed offsets.
One final issue is that currently one cannot enter a table of offsets for every station. This produces an error, because the array used to create offset bounds for the optimizer has zero length if the offset array has no floating values. To get round this, the case can be trapped with the addition of 3 lines:
# Generate an initial guess for stations which have no fixed time
p0 = np.zeros(len([val for val in self.stations_time_dict.values() if val is True]))
if self.verbose:
print('Initial function evaluation:', timingResiduals(p0, observations,
self.stations_time_dict,
weights=weights))
# Set bounds for timing to +/- given maximum time offset
bounds = []
# Trap the case where all stations are supplied with a time offset value.
if len([val for val in self.stations_time_dict.values() if val is True]) == 0:
p0 = np.zeros(1)
for i in range(len(p0)):
bounds.append([-self.max_toffset, self.max_toffset])
One reason one may wish to manually enter offsets for all cameras is if some camera timings were very inaccurate, and there was some other timing information to allow them to be anchored accurately to UT. This could arise if an external video was read as a camera station, with a large timing error, yet the solver decided to make it the reference camera with zero offset.
It is essential when analysing long duration fireballs to be able to experiment with fixed offsets which may be subtly or less subtly varied from values found by the solver. This is to allow accurate manual adjustment of timing offsets in circumstances where the solver is unable to fit the lag curve perfectly. This happens particularly where not all the data sets overlap. Sometimes also if there is very low drag in the case of an Earthgrazer, the solver will confuse camera timing error with deceleration and correct these when the absolute camera timings are known to be more accurate.