From 788879501b7f4da72bc23e92bba4b1032a07071e Mon Sep 17 00:00:00 2001
From: "oceanpdoshi@g.ucla.edu" <oceanpdoshi@.ucla.edu>
Date: Sun, 19 May 2019 23:12:23 -0700
Subject: [PATCH] In progress debugging

---
 ...perative_Localization_Journal2_sim_data.py | 110 ++++++++++++++++++
 CoLo-AT/comparison_plot.py                    |  55 +++++++++
 CoLo-AT/dataset_manager/data_generator.py     |  34 ++++--
 .../simulated_dataset_manager.py              |  32 +++--
 4 files changed, 209 insertions(+), 22 deletions(-)
 create mode 100644 CoLo-AT/Cooperative_Localization_Journal2_sim_data.py
 create mode 100644 CoLo-AT/comparison_plot.py

diff --git a/CoLo-AT/Cooperative_Localization_Journal2_sim_data.py b/CoLo-AT/Cooperative_Localization_Journal2_sim_data.py
new file mode 100644
index 0000000..23cdf34
--- /dev/null
+++ b/CoLo-AT/Cooperative_Localization_Journal2_sim_data.py
@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May 16 15:41:51 2019
+
+@author: william
+"""
+
+import os, sys
+import getpass
+from math import pi
+sys.path.append(os.path.join(os.path.dirname(__file__), "."))
+from dataset_manager.realworld_dataset_manager import RW_Dataset_Manager
+from dataset_manager.simulated_dataset_manager import SimulatedDataSetManager
+from simulation_process.sim_manager import  SimulationManager
+from robots.robot_system import RobotSystem
+from simulation_process.state_recorder import StatesRecorder
+from data_analysis.data_analyzer import Analyzer
+from data_analysis.realtime_plot import animate_plot
+
+import comparison_plot
+
+# load algorithms 
+sys.path.append(os.path.join(os.path.dirname(__file__), "localization_algos"))
+from centralized_ekf import Centralized_EKF # works
+from ekf_ls_bda import EKF_LS_BDA
+from ekf_ls_ci import EKF_LS_CI
+from ekf_gs_ci2 import EKF_GS_CI2
+from gs_ci_bound import GS_CI_Bound
+from ekf_gs_sci2 import EKF_GS_SCI2
+
+landmarks = {10: [0,0] }
+
+robot_labels = [1,2,3,4,5]
+boundary = 200 
+duration = 30
+testing_dataset = SimulatedDataSetManager('testing', boundary, landmarks, duration, robot_labels, 
+velocity_noise=0.0125, angular_velocity_noise=0.0, measurement_range_noise=0.1, bearing_noise=0.0349, robot_fov=2*pi)
+start_time, starting_states, dataset_data, time_arr = testing_dataset.simulate_dataset('random', test=False)
+
+analyzer = Analyzer('analyzer', robot_labels)
+
+#############################################################################
+
+loc_algo = EKF_GS_CI2('algo')
+robot = RobotSystem('robot', robot_labels, loc_algo, distr_sys = True)
+sim = SimulationManager('sim gs_ci')
+state_recorder_gs_ci = StatesRecorder('GS_CI', robot_labels)
+
+sim.sim_process_native(robot_labels, testing_dataset, robot, state_recorder_gs_ci, simple_plot = False)
+loc_err_per_run, trace_per_run, t_arr = analyzer.calculate_loc_err_and_trace_state_variance_per_run(state_recorder_gs_ci, plot_graphs = False)
+
+##############################################################################
+
+
+loc_algo = GS_CI_Bound('algo')
+robot = RobotSystem('robot', robot_labels, loc_algo, distr_sys = True)
+
+sim = SimulationManager('sim gs_ci_bound')
+state_recorder_bound = StatesRecorder('GS_CI_bound', robot_labels, state_var_only = True)
+sim.sim_process_native(robot_labels, testing_dataset, robot, state_recorder_bound, simple_plot = False)
+loc_err_per_run, trace_per_run, t_arr = analyzer.calculate_loc_err_and_trace_state_variance_per_run(state_recorder_bound, plot_graphs = False)
+
+
+##############################################################################
+
+
+loc_algo = Centralized_EKF('algo')
+robot = RobotSystem('robot', robot_labels, loc_algo, distr_sys = False)
+
+sim = SimulationManager('sim cen_ekf')
+state_recorder_LS_cen = StatesRecorder('LS_cen', robot_labels)
+sim.sim_process_native(robot_labels, testing_dataset, robot, state_recorder_LS_cen, simple_plot = False)
+loc_err_per_run, trace_per_run, t_arr = analyzer.calculate_loc_err_and_trace_state_variance_per_run(state_recorder_LS_cen, plot_graphs = False)
+
+##############################################################################
+
+
+loc_algo = EKF_LS_BDA('algo')
+robot = RobotSystem('robot', robot_labels, loc_algo, distr_sys = False)
+
+sim = SimulationManager('sim ls_bda')
+state_recorder_LS_BDA = StatesRecorder('LS_BDA', robot_labels)
+sim.sim_process_native(robot_labels, testing_dataset, robot, state_recorder_LS_BDA, simple_plot = False)
+loc_err_per_run, trace_per_run, t_arr = analyzer.calculate_loc_err_and_trace_state_variance_per_run(state_recorder_LS_BDA, plot_graphs = False)
+
+
+##############################################################################
+
+
+loc_algo = EKF_LS_CI('algo')
+robot = RobotSystem('robot', robot_labels, loc_algo, distr_sys = False)
+
+sim = SimulationManager('sim ls_ci')
+state_recorder_LS_CI = StatesRecorder('LS_CI', robot_labels)
+sim.sim_process_native(robot_labels, testing_dataset, robot, state_recorder_LS_CI, simple_plot = False)
+loc_err_per_run, trace_per_run, t_arr = analyzer.calculate_loc_err_and_trace_state_variance_per_run(state_recorder_LS_CI, plot_graphs = False)
+
+##############################################################################
+
+
+loc_algo = EKF_GS_SCI2('algo')
+robot = RobotSystem('robot', robot_labels, loc_algo, distr_sys = True)
+
+sim = SimulationManager('sim gs_sci')
+state_recorder_GS_SCI = StatesRecorder('GS_SCI', robot_labels)
+sim.sim_process_native(robot_labels, testing_dataset, robot, state_recorder_GS_SCI, simple_plot = False)
+loc_err_per_run, trace_per_run, t_arr = analyzer.calculate_loc_err_and_trace_state_variance_per_run(state_recorder_GS_SCI, plot_graphs = False)
+
+comparison_plot.robot1_algo_comparison_plot([state_recorder_gs_ci, state_recorder_LS_cen, state_recorder_LS_BDA, state_recorder_LS_CI, state_recorder_GS_SCI, state_recorder_bound], only_trace = ['GS_CI_bound'])
\ No newline at end of file
diff --git a/CoLo-AT/comparison_plot.py b/CoLo-AT/comparison_plot.py
new file mode 100644
index 0000000..6262f43
--- /dev/null
+++ b/CoLo-AT/comparison_plot.py
@@ -0,0 +1,55 @@
+from matplotlib import pyplot as plt
+import IPython
+
+# [time_stamps, loc_err_per_run, data_recorder.name]
+# [time_stamps, trace_per_run, data_recorder.name]
+def generate_comparison_plot(arr_loc_err, arr_trace, only_trace):
+    fig = plt.figure()
+    fig1 = fig.add_subplot(211)
+    fig2 = fig.add_subplot(212)
+    for loc_err in arr_loc_err:
+        fig1.plot(loc_err[0], loc_err[1], label = loc_err[2])
+    
+    for trace in arr_trace:
+        if (only_trace != None) and trace[2] in only_trace:
+            fig2.plot(trace[0], trace[1], '--' ,label = trace[2])
+        else:
+            fig2.plot(trace[0], trace[1], label = trace[2])
+
+    fig1.set_title('Location error', fontsize = 24)
+    fig1.set_xlabel('Time[s]', fontsize = 22)
+    fig1.set_ylabel('RMS[m]', fontsize = 22) 
+    fig1.set_ylim(0, 1)
+    fig1.legend(loc=1, fontsize = 16)
+    fig1.tick_params(labelsize=18)
+
+    fig2.set_title('Trace of state covariance', fontsize = 24)
+    fig2.set_xlabel('Time [s]', fontsize = 22)
+    fig2.set_ylabel('Sigma_s [m^2]', fontsize = 22)
+    fig2.set_ylim(0, 0.1)
+    fig2.legend(loc=1, fontsize = 20)
+    fig2.tick_params(labelsize=18)
+
+    # plt.subplots_adjust(hspace=.6)
+    plt.show()
+    return
+
+def robot1_algo_comparison_plot(arr_data_recorder, only_trace=None):
+    print("************Algorithm Comparison***************")
+    arr_loc_err = []
+    arr_trace = []
+
+    for data_recorder in arr_data_recorder:
+        r1_label = data_recorder.get_dataset_labels()[0]
+
+        loc_err_per_run = data_recorder.loc_err_arr[r1_label]
+        trace_per_run = data_recorder.trace_sigma_s_arr[r1_label]
+        time_stamps = data_recorder.get_time_arr(r1_label)
+
+        if only_trace == None or data_recorder.name not in only_trace:
+            arr_loc_err.append([time_stamps, loc_err_per_run, data_recorder.name])
+        arr_trace.append([time_stamps, trace_per_run, data_recorder.name] )
+    print('Plotting Comparison Graphs')
+    generate_comparison_plot(arr_loc_err, arr_trace, only_trace)
+    return arr_loc_err, arr_trace
+
diff --git a/CoLo-AT/dataset_manager/data_generator.py b/CoLo-AT/dataset_manager/data_generator.py
index 207210c..545703b 100644
--- a/CoLo-AT/dataset_manager/data_generator.py
+++ b/CoLo-AT/dataset_manager/data_generator.py
@@ -1,7 +1,9 @@
 import numpy as np
-from math import pi, sqrt, atan2, hypot, sin, cos
+from math import pi, sqrt, atan2, hypot, sin, cos, pow, trunc
 import matplotlib.pyplot as plt
 
+import IPython
+
 # TODO - add collision/boundary checking for circle/random data
 class DataGenerator():
     def __init__(self, boundary, landmarks, duration, robot_labels, starting_states, start_time, delta_t,
@@ -35,10 +37,12 @@ class DataGenerator():
             if landmark_id < 5:
                 raise Exception("Invalid landmark ID: landmark IDs must be bigger than 5.")
 
+        # TODO - user specified sample rate
         self.time_arr = {'odometry': [[] for robot in robot_labels], 'measurement': [[] for robot in robot_labels], 'groundtruth': [[] for robot in robot_labels]}
         for data_type in self.time_arr:
             for i in range(self.num_robots):
-                self.time_arr[data_type][i] = np.arange(self.start_time, self.end_time)
+                arr = np.arange(self.start_time, self.end_time, self.delta_t)
+                self.time_arr[data_type][i] = arr
 
         self.groundtruth_data = [ [] for robot in range(self.num_robots)]
         self.odometry_data = [[] for robot in range(self.num_robots)]
@@ -92,7 +96,7 @@ class DataGenerator():
         path_circumference = 2*pi*path_radius
 
         radian_step = 2*pi/self.duration
-        velocity_step = path_circumference/self.duration
+        velocity_step = path_circumference/self.duration # May want to randomize velocity w/normal distribution
         distance_step = velocity_step*self.delta_t
 
         # Generate ground truth data
@@ -141,15 +145,15 @@ class DataGenerator():
                 curr_y = self.groundtruth_data[i][time_idx-1]['y_pos']
                 curr_orientation = self.groundtruth_data[i][time_idx-1]['orientation']
 
-                # rotate some random number of radians between pi/4 and -pi/4
+                # rotate some random number of radians between pi/3 and -pi/3
                 next_orientation = curr_orientation + np.random.uniform(-pi/3, pi/3)
                 next_orientation = self.converge_to_angle_range(next_orientation)
-                rand_dist = np.random.uniform(0, 0.5) # units of distance change are in 
+                rand_dist = 1 + np.random.normal(0, sqrt(0.25)) # units of distance change are in meters, if delta t is one second... v = m/s
                 next_x = curr_x + rand_dist*cos(next_orientation)
                 next_y = curr_y + rand_dist*sin(next_orientation)
                 
                 self.groundtruth_data[i].append({'time' : self.time_arr['groundtruth'][i][time_idx], 'x_pos': next_x, 'y_pos': next_y, 'orientation': next_orientation})
-       
+
         self.generate_odometry_data()
         self.generate_measurement_data()
 
@@ -278,6 +282,10 @@ class DataGenerator():
 
                         meas_diff = actual_meas - generated_measurement
                         bearing_diff = actual_bearing - generated_bearing
+                        
+                        # double check for edge cases where both are approximately pi (one may be (-), the other may be positive
+                        if (abs(bearing_diff) > 1):
+                            print(actual_bearing, generated_bearing)
 
                         meas_test_arr[i].append({'meas_diff': meas_diff, 'bearing_diff': bearing_diff, 'time': curr_time})
                     else:
@@ -295,6 +303,10 @@ class DataGenerator():
                         meas_diff = actual_meas - generated_measurement
                         bearing_diff = actual_bearing - generated_bearing
 
+                        # double check for edge cases where both are approximately pi (one may be (-), the other may be positive
+                        if (abs(bearing_diff) > 1):
+                            print(actual_bearing, generated_bearing)
+
                         meas_test_arr[i].append({'meas_diff': meas_diff, 'bearing_diff': bearing_diff, 'time': curr_time})
 
         fig, ax_arr = plt.subplots(6)
@@ -330,10 +342,10 @@ class DataGenerator():
             ax_arr[5].scatter(meas_time_arr, bearing_diff_arr)
 
         fig.legend(labels, prop={'size': 12})
-        ax_arr[0].set_ylim([-10*self.velocity_noise, 10*self.velocity_noise])
-        ax_arr[1].set_ylim([-10*self.angular_velocity_noise, 10*self.angular_velocity_noise])
-        ax_arr[4].set_ylim([-10*self.measurement_range_noise, 10*self.measurement_range_noise])
-        ax_arr[5].set_ylim([-10*self.bearing_noise, 10*self.bearing_noise])
+        # ax_arr[0].set_ylim([-10*self.velocity_noise, 10*self.velocity_noise])
+        # ax_arr[1].set_ylim([-10*self.angular_velocity_noise, 10*self.angular_velocity_noise])
+        # ax_arr[4].set_ylim([-10*self.measurement_range_noise, 10*self.measurement_range_noise])
+        # ax_arr[5].set_ylim([-10*self.bearing_noise, 10*self.bearing_noise])
         
         plt.show()
         plt.close('all')
@@ -364,6 +376,8 @@ class DataGenerator():
     def within_vision(self, robot_fov, robot_loc, landmark_loc, robot_orientation=0):    
         # between pi and -pi with respect to LOCAL FRAME
         bearing = atan2((landmark_loc[1]-robot_loc[1]),(landmark_loc[0]-robot_loc[0])) - robot_orientation
+        if (robot_fov == 2*pi):
+            return True
         return abs(bearing) <= robot_fov/2
 
     # get angle between -pi and pi
diff --git a/CoLo-AT/dataset_manager/simulated_dataset_manager.py b/CoLo-AT/dataset_manager/simulated_dataset_manager.py
index e51247b..6ac1925 100644
--- a/CoLo-AT/dataset_manager/simulated_dataset_manager.py
+++ b/CoLo-AT/dataset_manager/simulated_dataset_manager.py
@@ -13,6 +13,13 @@ sys.path.insert(0,parentdir)
 
 from dataset_manager.data_generator import DataGenerator
 
+def find_nearest_time_idx(l, value):
+    if len(l) != 0:
+        array = np.asarray(l)
+        idx = (np.abs(array-value)).argmin()
+        return idx
+    else:
+        return None
 
 # TODO 
 # communication
@@ -36,8 +43,8 @@ class SimulatedDataSetManager:
         self.curr_request = 'odometry'
         self.curr_robot_idx = 0
         
-        # time step of the simulation (secs)
-        self.delta_t = 1.0 
+        # time step of the simulation (secs) TODO - user specified
+        self.delta_t = 1.0
 
         # landmarks = {landmark ID: [x,y]} (m)
         self.landmark_map = landmarks
@@ -47,7 +54,7 @@ class SimulatedDataSetManager:
 
         # starting_states =	{robot label: [time, x_pos, y_pos, orientation], ... }
         # TODO - create a generate starting states function (potentially randomized to fit boundary)
-        self.starting_states = {label : [self.start_time, -boundary + float(i+1)*boundary/self.num_robots, 0, 0.0] for i,label in enumerate(self.robot_labels)}
+        self.starting_states = {label : [self.start_time, i, 0, 0.0] for i,label in enumerate(self.robot_labels)}
 
         # groundtruth & odometry are of fixed length=duration
         # measurement_data is of variable length
@@ -71,7 +78,6 @@ class SimulatedDataSetManager:
             return self.start_time, self.starting_states, self.dataset_data, self.time_arr
         elif(path == 'random'):
             self.time_arr, self.groundtruth_data, self.odometry_data, self.measurement_data, self.dataset_data = self.generator.generate_random_data(test)
-            self.measurement_data[0] = [] # TODO - to test simulated communication, trace of state variance should increase unless a communication update
             return self.start_time, self.starting_states, self.dataset_data, self.time_arr
         else:
             raise Exception("Unsupported robot path specified.")
@@ -162,16 +168,17 @@ class SimulatedDataSetManager:
         valid_dataline = current_time <= self.end_time
 
         # Select message data
+        # NOTE - use difference rather than inequality because of float difficulties
         if valid_dataline:
             if req_type == 'odometry':
-                odo_time_idx = np.where(self.time_arr[req_type][robot_idx] == current_time)[0][0] 
+                odo_time_idx = find_nearest_time_idx(self.time_arr[req_type][robot_idx], current_time)
                 message['data'] = self.odometry_data[robot_idx][odo_time_idx]
             # if no measurement data available (ex: no landmarks were in range) returns
             # {'time': -1, 'subject_ID': None, 'measurment_range': None,'bearing': None}
             elif req_type == 'measurement':
                 message['data'] = self.find_measurement(current_time, robot_idx)
 
-        gt_time_idx = np.where(self.time_arr['groundtruth'][robot_idx] == current_time)[0][0]
+        gt_time_idx = find_nearest_time_idx(self.time_arr['groundtruth'][robot_idx], current_time)
         message['groundtruth'] = self.groundtruth_data[robot_idx][gt_time_idx]
 
         return valid_dataline, message, req_type, robot_idx
@@ -181,18 +188,19 @@ class SimulatedDataSetManager:
     #           {'time': current_time, 'subject_ID': None, 'measurment_range': None,'bearing': None} if no measurement available
     def find_measurement(self, current_time, robot_idx):
         
-        # {'idx': idx, 'distance': measurment_range}
         available_measurements = []
+        measurement_indices = []
         
         for data_idx, measurement in enumerate(self.measurement_data[robot_idx]):
-            if (measurement['time'] == current_time):
-                available_measurements.append({'idx': data_idx, 'distance': measurement['measurment_range']})
+            if (abs(measurement['time'] - current_time) < 0.001):
+                available_measurements.append(measurement['measurment_range'])
+                measurement_indices.append(data_idx)
 
         # Select the "best"=shortest distance measurement
         if (len(available_measurements) != 0):
-            sorted_measurements = sorted(available_measurements, key=lambda k: k['distance'])
-            min_distance_idx = sorted_measurements[0]['idx']
-            return self.measurement_data[robot_idx][min_distance_idx]
+            least_index = np.argmin(available_measurements)
+            idx = measurement_indices[least_index]
+            return self.measurement_data[robot_idx][idx]
 
         # No measurement data available (ex: no landmarks in vision at given time)
         return {'time': current_time, 'subject_ID': None, 'measurment_range': None,'bearing': None}
\ No newline at end of file
-- 
GitLab