############################################################################### # Copyright (c) 2011 - 2012 Trimble Navigation Ltd # $Id: NoPiSG_Process.py,v 1.3 2012/03/08 04:38:46 shankar Exp $ ############################################################################### # # NoPiSG_Process.py # # Module to load a single day's diffs/acq_diffs_combo file(s), compute # statistics for each combo file and write the statistics to a file # ############################################################################### import sys import os from NoPiSG_Utils import add_utils_dir_to_path add_utils_dir_to_path() from NoPiUT_Common_Meas import * import shlex # Import matplotlib to set the matplotlib backend. In particular, must be done # before importing pylab import matplotlib matplotlib.use('agg') from pylab import * from numpy import * from scipy import stats IDX_INVALID = -1 ############################################################################### # Class definition containing measurement statistics ############################################################################### class cl_stats : def __init__( self ) : self.sv_id = 0 self.tot_epochs = 0 self.mean = 0.0 self.median = 0.0 self.std = 0.0 self.madn = 0.0 self.skewness = 0.0 self.three_sigma = 0.0 self.rob_three_sigma = 0.0 ############################################################################### # Class definition containing normative values ############################################################################### class cl_norm_vals : def __init__( self ) : self.norm_mean = [] self.norm_std_dev = [] ############################################################################### # Wrapper function to load a diffs_combo file to memory. Actual fread # operations performed within load_combo_file() in # NoPi_Common/NoPi_Utils.py ############################################################################### def load_combo( diffs_summ, fname ) : global din global flags [ din, flags ] = load_combo_file( diffs_summ, fname ) ############################################################################### # Parse a file containing normative ( population ) mean and std. deviation # values. Store values and corresponding combo name to global variables ############################################################################### def parse_norm_vals( norm_vals_file ) : global norm_vals global combo_names parse_ok = False # Use three separate objects to store normative values for each meas type norm_vals = [ cl_norm_vals(), cl_norm_vals(), cl_norm_vals() ] combo_names = [ ] # If the user specified normative values file exists, parse it. Else # return False if os.access( norm_vals_file, os.R_OK ) : fin = open( norm_vals_file, 'rt' ) norm_values = fin.readlines() fin.close() for value in norm_values : tmp_val = shlex.split( value, True ) len_tmp_val = len( tmp_val ) if ( len_tmp_val > 0 ) : # Not a comment line # Normative values line should have 10 columns if ( len_tmp_val == 10 ) : tmp_combo_name = '%s %s %s %s' % \ ( tmp_val[ len_tmp_val-4 ], \ tmp_val[ len_tmp_val-3 ], \ tmp_val[ len_tmp_val-2 ], \ tmp_val[ len_tmp_val-1 ] ) combo_names.append( tmp_combo_name ) # Store normative values for each meas type for mtype in range( 3 ) : norm_vals[ mtype ].norm_mean.append( float( tmp_val[ mtype*2 ] ) ) norm_vals[ mtype ].norm_std_dev.append( float( tmp_val[ mtype*2 + 1 ] ) ) parse_ok = True else : print 'Normative values file format mismatch' print 'Please refer to the README for the expected file format' print 'Verify: \n %s\'s column format' % norm_vals_file print 'Row with column format mismatch is: \n %s' % value # Flag a warning if the parsed file does not have valid normative values. if ( not( parse_ok ) ) : print 'The specified normative values file does not contain valid normative values' print 'Using sample mean/std dev for all valid combos' # Flag a warning if the user specified norm vals file cannot be located else : print 'Normative values file cannot be found' print 'Using sample mean/std dev. for all valid combos' return( parse_ok ) ############################################################################### # Process NoPi results for a given combo file and measurement type. Compute # statistics over the entire dataset comprising of all satellites observed ############################################################################### def process_data( diffs_summ, combo, meas_type, data_path, data_type, norm_idx ) : data_col = get_data_column( diffs_summ, meas_type ) data_scale = get_data_scale( meas_type ) din[:, data_col] = din[:, data_col] * data_scale; sv = 0 compute_stats( sv, din[ :, data_col ], flags, combo, meas_type, data_path, data_type, norm_idx ) ############################################################################### # Compute a set of statistics to evaluate receiver performance. The statistics # are written to a file for subsequent processing of further analysis in the # event of a receiver anomaly ############################################################################### def compute_stats( sv, din_sv, flags_sv, combo, meas_type, data_path, data_type, norm_idx ) : global idx_vld three_sigma_count = 0 obs_stats = cl_stats() flags_msk = get_flags_mask( meas_type ) obs_stats.sv = sv idx_vld = find( flags_sv & flags_msk ) # Compute statistics only if there are at least 10 valid observations. # Identical threshold used in NoPi_Display/NoPiDI_Process.py as well. if ( len( idx_vld ) > 10 ) : obs_stats.tot_epochs = len( idx_vld ) # Mean and Median obs_stats.mean = mean( din_sv[ idx_vld ] ) obs_stats.median = median( din_sv[ idx_vld ] ) # Standard Deviation and MADN ( Robust Std ) obs_stats.std = std( din_sv[ idx_vld ] ) # MADN obs_stats.madn = compute_madn( din_sv[ idx_vld ], obs_stats.median ) # Skewness obs_stats.skewness = stats.skew( din_sv[ idx_vld ], bias='False' ) # Store % of epochs exceeding the 3-sigma rule. Determine count based on # normative and robust first and second moments # Use normative values if available. Else use sample mean/std deviation # values if ( norm_idx == IDX_INVALID ) : tmp_mean = obs_stats.mean tmp_std = obs_stats.std else : tmp_mean = norm_vals[ meas_type ].norm_mean[ norm_idx ] tmp_std = norm_vals[ meas_type ].norm_std_dev[ norm_idx ] obs_stats.three_sigma = compute_three_sigma_perc( din_sv[ idx_vld ], tmp_mean, tmp_std ) # Equivalent robust 3-sigma percentage obs_stats.rob_three_sigma = compute_three_sigma_perc( din_sv[ idx_vld ], obs_stats.median, obs_stats.madn ) # Write observation statistics to stats_combo file write_results( obs_stats, combo, meas_type, data_path, data_type, norm_idx ) return() ############################################################################### # MADN: Median Absolute Deviation - Normalized. 0.6745 is necessary to # normalize the standard deviation of the distribution to a N(0,1) # Ref: Robust Statistics: Theory and Methods, Ricardo Maronna et.al. Chapter 1 ############################################################################### def compute_madn( din_sv, obs_median ) : dev_from_median = fabs( array( din_sv ) - obs_median ) madn = median( array( dev_from_median ) )/0.6745 return( madn ) ############################################################################### # Compute percentage of data exceeding 3-sigma Expects observation data, # expected value of random variable and measure of spread ############################################################################### def compute_three_sigma_perc( din_sv, exp_val, spread ) : normalized_din_sv = fabs(true_divide( ( din_sv - exp_val ), spread ) ) count_three_sigma = len( find( normalized_din_sv > 3 ) ) three_sigma_perc = true_divide( count_three_sigma, len( din_sv ) )*100 return( three_sigma_perc ) ############################################################################### # Write computed stats to file ############################################################################### def write_results( obs_stats, combo, meas_type, data_path, data_type, norm_idx ) : if ( data_type == 'DATA' ) : fname = data_path + 'stats_combo_' + str( int( combo ) ) + '.mtb' elif ( data_type == 'ACQ_DATA' ) : fname = data_path + 'acq_stats_combo_' + str( int( combo ) ) + '.mtb' else : print 'Unsupported stats data type' return() # If the required file exists, append to the file. If the data is being # reprocessed, any existing file will be overwritten if ( os.path.exists( fname ) and meas_type != 0 ) : fout = open( fname, 'a' ) else : fout = open( fname, 'wt' ) # Enable Bit Flag if using sample mean/std deviation instead of # normative values to compute three sigma percentage if ( norm_idx == IDX_INVALID ) : bit_flag = 1 else : bit_flag = 0 fout.write('%2d %d %d %d %d %6d %3.3f %3.3f %3.3f %3.3f %2.2f %4.3f %4.3f\n'% ( obs_stats.sv, meas_type, 0, # Elevation Mask 0, # CNo Mask bit_flag, # Bit Flag: # b0 = TRUE = Use sample mean and sample std # deviation to compute three sigma percentage # FALSE = Use normative values of mean # and std deviation to compute three sigma # percentage # b1:b31 = unused obs_stats.tot_epochs, obs_stats.mean, obs_stats.median, obs_stats.std, obs_stats.madn, obs_stats.skewness, obs_stats.three_sigma, obs_stats.rob_three_sigma ) ) fout.close() return() ############################################################################### # Return index of combo in a user specified normative value file ############################################################################### def get_norm_val_idx( cur_combo_name, parse_ok ) : combo_idx = IDX_INVALID # Only look for combo number corresponding to cur_combo_name if # normative values were read from a user specified file if parse_ok : for idx, name in enumerate( combo_names ) : if( cur_combo_name == name ) : combo_idx = idx break # Return -1 as the default index indicating the requested combo # cannot be found. This case happen either because: # 1. The combo's normative values are not specified in the provided normative # values file or # 2. The specified file does not exists. return( combo_idx )