text:
import time
import numpy as np
from pycqed.analysis import analysis_toolbox as a_tools
import pycqed.analysis_v2.base_analysis as ba
# import dataprep for tomography module
# import tomography module
# using the data prep module of analysis V2
# from pycqed.analysis_v2 import tomography_dataprep as dataprep
from pycqed.analysis import measurement_analysis as ma
try:
    import qutip as qt
except ImportError as e:
    pass
    # logging.warning('Could not import qutip, tomo code will not work')

def reshape_block(shots_data, segments_per_block=16, block_size=4092, mode='truncate'):
    """
    inputs: shots_data 1D array of dimension N
    organizes data in blocks of dimension block_size.
    num of blocks is N/block_size
    """
    N = len(shots_data)
    # Data dimension needs to be an integer multiple of block_size
    assert(N%block_size==0)
    num_blocks = N//block_size
    full_segments = block_size//segments_per_block
    orfan_segments = block_size % segments_per_block
    missing_segments = segments_per_block - orfan_segments
#     print(N,num_blocks,full_segments,orfan_segments,missing_segments)
    reshaped_data = shots_data.reshape((num_blocks,block_size))
    if mode.lower()=='truncate':
        truncate_idx = full_segments*segments_per_block
        return reshaped_data[:,:truncate_idx]
    elif mode.lower()=='padd':
        padd_dim = (full_segments+1)*segments_per_block
        return_block = np.nan*np.ones((num_blocks,padd_dim))
        return_block[:,:block_size] = reshaped_data
        return return_block
    else:
        raise ValueError('Mode not understood. Needs to be truncate or padd')

def all_repetitions(shots_data,segments_per_block=16):
    flat_dim = shots_data.shape[0]*shots_data.shape[1]
    # Data dimension needs to divide the segments_per_block
    assert(flat_dim%segments_per_block==0)
    num_blocks = flat_dim // segments_per_block
    block_data = shots_data.reshape((num_blocks,segments_per_block))
    return block_data

def get_segments_average(shots_data, segments_per_block=16, block_size=4092, mode='truncate', average=True):
    reshaped_data = reshape_block(shots_data=shots_data,
                                      segments_per_block=segments_per_block,
                                      block_size=block_size,
                                      mode=mode)
    all_reps = all_repetitions(shots_data=reshaped_data,
                                       segments_per_block=segments_per_block)
    if average:
        return np.mean(all_reps,axis=0)
    else:
        return all_reps

class ExpectationValueCalculation_mmt_computer:

    def __init__(self, auto=True, label='', timestamp=None,
                 fig_format='png',
                 q0_label='q0',
                 q1_label='q1', close_fig=True, **kw):
        self.label = label
        self.timestamp = timestamp
        self.fig_format = fig_format
        # q0 == D2
        self.q0_label = q0_label
        # q1 == A
        self.q1_label = q1_label
        self.n_states = 2 ** 2
        self.ma_obj = ma.MeasurementAnalysis(auto=False, label=label,
                                             timestamp=timestamp)
        self.ma_obj.get_naming_and_values()
        # self.get_naming_and_values()
        # hard coded number of segments for a 2 qubit state tomography
        # constraint imposed by UHFLI
        self.nr_segments = 16
        # self.exp_name = os.path.split(self.folder)[-1][7:]

        avg_h1 = self.ma_obj.measured_values[0]
        avg_h2 = self.ma_obj.measured_values[1]
        avg_h12 = self.ma_obj.measured_values[2]
        h1_00 = np.mean(avg_h1[36:36+7])
        h1_01 = np.mean(avg_h1[43:43+7])
        h1_10 = np.mean(avg_h1[50:50+7])
        h1_11 = np.mean(avg_h1[57:])

        h2_00 = np.mean(avg_h2[36:36+7])
        h2_01 = np.mean(avg_h2[43:43+7])
        h2_10 = np.mean(avg_h2[50:50+7])
        h2_11 = np.mean(avg_h2[57:])

        h12_00 = np.mean(avg_h12[36:36+7])
        h12_01 = np.mean(avg_h12[43:43+7])
        h12_10 = np.mean(avg_h12[50:50+7])
        h12_11 = np.mean(avg_h12[57:])

        measurement_channel_1 = np.array([avg_h1[0], avg_h1[1], avg_h1[7],
                                          avg_h1[8], avg_h1[14], avg_h1[21],
                                          avg_h1[28], avg_h1[35]])
        measurement_channel_2 = np.array([avg_h2[0], avg_h2[1], avg_h2[7], avg_h2[8],
                                          avg_h2[14], avg_h2[21], avg_h2[28], avg_h2[35]])
        measurement_channel_3 = np.array([avg_h12[0], avg_h12[1], avg_h12[7],
                                          avg_h12[8], avg_h12[14], avg_h12[21],
                                          avg_h12[28],avg_h12[35]])
        self.measurements_tomo = np.array([measurement_channel_1,
                                           measurement_channel_2,
                                           measurement_channel_3]).flatten()
        # print(self.measurements_tomo)
        # print(len(self.measurements_tomo))

        # 108 x 1
        # get the calibration points by averaging over the five measurements
        # taken knowing the initial state we put in
        self.measurements_cal = np.array(
            [h1_00, h1_01, h1_10, h1_11,
             h2_00, h2_01, h2_10, h2_11,
             h12_00, h12_01, h12_10, h12_11])

    def _calibrate_betas(self):
        """
        calculates betas from calibration points for the initial measurement
        operator

        Betas are ordered by B0 -> II B1 -> IZ etc(binary counting)
        <0|Z|0> = 1, <1|Z|1> = -1

        Keyword arguments:
        measurements_cal --- array(2 ** n_qubits) should be ordered
            correctly (00, 01, 10, 11) for 2 qubits
        """
        cal_matrix = np.zeros((self.n_states, self.n_states))
        # get the coefficient matrix for the betas
        for i in range(self.n_states):
            for j in range(self.n_states):
                # perform bitwise AND and count the resulting 1s
                cal_matrix[i, j] = (-1)**(bin((i & j)).count("1"))
        # invert solve the simple system of equations
        # print(cal_matrix)
        # print(np.linalg.inv(cal_matrix))
        betas = np.zeros(12)
        # print(self.measurements_cal[0:4])
        betas[0:4] = np.dot(np.linalg.inv(cal_matrix),
                            self.measurements_cal[0:4])
        self.betas_up = betas[0:4]
        betas[4:8] = np.dot(np.linalg.inv(cal_matrix),
                            self.measurements_cal[4:8])
        self.betas_p = betas[4:8]
        betas[8:] = np.dot(np.linalg.inv(cal_matrix),
                           self.measurements_cal[8:12])
        self.betas_pp = betas[8:]
        return betas

    def assemble_M_matrix_single_block(self, beta_array):
        M_matrix_single_block_row_1 = np.array([beta_array[0], beta_array[1],
                                                beta_array[2], beta_array[3],
                                                0, 0, 0, 0, 0, 0])
        M_matrix_single_block_row_2 = np.array([beta_array[0],
                                                -1*beta_array[1],
                                                beta_array[2],
                                                -1*beta_array[3],
                                                0, 0, 0, 0, 0, 0])
        M_matrix_single_block_row_3 = np.array([beta_array[0],
                                                beta_array[1],
                                                -1*beta_array[2],
                                                -1*beta_array[3],
                                                0, 0, 0, 0, 0, 0])
        M_matrix_single_block_row_4 = np.array([beta_array[0],
                                                -1*beta_array[1],
                                                -1*beta_array[2],
                                                beta_array[3],
                                                0, 0, 0, 0, 0, 0])
        M_matrix_single_block_row_5 = np.array([beta_array[0],
                                                0, 0, 0, -beta_array[1],
                                                -beta_array[2],
                                                beta_array[3], 0, 0, 0])
        M_matrix_single_block_row_6 = np.array([beta_array[0], 0, 0, 0,
                                                beta_array[1],
                                                beta_array[2],
                                                beta_array[3],
                                                0, 0, 0])
        M_matrix_single_block_row_7 = np.array([beta_array[0], 0, 0,
                                                0, 0, 0, 0, beta_array[1],
                                                beta_array[2],
                                                beta_array[3]])
        M_matrix_single_block_row_8 = np.array([beta_array[0], 0, 0, 0, 0,
                                                0, 0, -beta_array[1],
                                                -beta_array[2],
                                                beta_array[3]])
        M_matrix_single_block = np.vstack((M_matrix_single_block_row_1,
                                           M_matrix_single_block_row_2,
                                           M_matrix_single_block_row_3,
                                           M_matrix_single_block_row_4,
                                           M_matrix_single_block_row_5,
                                           M_matrix_single_block_row_6,
                                           M_matrix_single_block_row_7,
                                           M_matrix_single_block_row_8))
        M_matrix_single_block = M_matrix_single_block.reshape(8, 10)
        return M_matrix_single_block

    def assemble_M_matrix(self):
        Block1 = self.assemble_M_matrix_single_block(self.betas_up)
        Block2 = self.assemble_M_matrix_single_block(self.betas_p)
        Block3 = self.assemble_M_matrix_single_block(self.betas_pp)
        self.M_matrix = np.vstack((Block1, Block2, Block3)).reshape(24, 10)
        return self.M_matrix

    def invert_M_matrix(self):
        self.inverse_matrix = np.linalg.pinv(self.M_matrix)
        return self.inverse_matrix

    def execute_error_signalling(self, ev):
        II = (ev[0] - ev[3])/(1 - ev[3])
        IZ = (ev[1] - ev[2])/(1 - ev[3])
        ZI = (ev[1] - ev[2])/(1 - ev[3])
        ZZ = (ev[3] - ev[0])/(1 - ev[3])
        XX = (ev[4] + ev[5])/(1 - ev[3])
        YY = (ev[4] + ev[5])/(1 - ev[3])
        ev_error_signalling = np.array([II, IZ, ZI, ZZ, XX, YY])
        return ev_error_signalling

    def execute_expectation_value_calculation(self):
        # assemble matrix that connects RO with terms
        self._calibrate_betas()
        self.assemble_M_matrix()
        self.invert_M_matrix()
        # use it to get terms back from RO
        rescaled_measurements_tomo = self.measurements_tomo
        self.expect_values = np.dot(self.inverse_matrix,
                                    rescaled_measurements_tomo)
        expect_values_VQE = np.array([self.expect_values[0],
                                      self.expect_values[1],
                                      self.expect_values[2],
                                      self.expect_values[3],
                                      self.expect_values[6],
                                      self.expect_values[9]])
        return expect_values_VQE

    def execute_expectation_value_calculation_traceone(self):
        # assemble matrix that connects RO with terms
        self._calibrate_betas()
        self.assemble_M_matrix()
        self.inverse_matrix = np.linalg.pinv(self.M_matrix[:, 1:])
        # use it to get terms back from RO
        rescaled_measurements_tomo = self.measurements_tomo
        self.expect_values = np.dot(self.inverse_matrix,
                                    rescaled_measurements_tomo)
        expect_values_VQE = np.array([1,
                                      self.expect_values[0],
                                      self.expect_values[1],
                                      self.expect_values[2],
                                      self.expect_values[5],
                                      self.expect_values[8]])
        return expect_values_VQE

    def execute_expectation_value_calculation_T1signaling(self):
        # assemble matrix that connects RO with terms
        self._calibrate_betas()
        self.assemble_M_matrix()
        self.inverse_matrix = np.linalg.pinv(self.M_matrix[:, 1:])
        # use it to get terms back from RO
        rescaled_measurements_tomo = self.measurements_tomo
        self.expect_values = np.dot(self.inverse_matrix,
                                    rescaled_measurements_tomo)
        expect_values_VQE = np.array([1,
                                      self.expect_values[0],
                                      self.expect_values[1],
                                      self.expect_values[2],
                                      self.expect_values[5],
                                      self.expect_values[8]])
        expect_values_VQE = self.execute_error_signalling(expect_values_VQE)
        return expect_values_VQE

class ExpectationValueCalculation:

    def __init__(self, auto=True, label='', timestamp=None,
                 fig_format='png',
                 q0_label='q0',
                 q1_label='q1', close_fig=True, **kw):
        self.label = label
        self.timestamp = timestamp
        self.fig_format = fig_format
        # q0 == D2
        self.q0_label = q0_label
        # q1 == A
        self.q1_label = q1_label
        self.n_states = 2 ** 2
        self.ma_obj = ma.MeasurementAnalysis(auto=False, label=label,
                                             timestamp=timestamp)
        self.ma_obj.get_naming_and_values()
        # self.get_naming_and_values()
        # hard coded number of segments for a 2 qubit state tomography
        # constraint imposed by UHFLI
        # self.nr_segments = 16
        # self.exp_name = os.path.split(self.folder)[-1][7:]

        avg_h1 = self.ma_obj.measured_values[0]
        avg_h2 = self.ma_obj.measured_values[1]
        avg_h12 = self.ma_obj.measured_values[2]

        #this should be implemented with a flag
        #but
        h1_00 = np.mean(avg_h1[36:36+7])
        h1_01 = np.mean(avg_h1[43:43+7])
        h1_10 = np.mean(avg_h1[50:50+7])
        h1_11 = np.mean(avg_h1[57:])

        h2_00 = np.mean(avg_h2[36:36+7])
        h2_01 = np.mean(avg_h2[43:43+7])
        h2_10 = np.mean(avg_h2[50:50+7])
        h2_11 = np.mean(avg_h2[57:])

        h12_00 = np.mean(avg_h12[36:36+7])
        h12_01 = np.mean(avg_h12[43:43+7])
        h12_10 = np.mean(avg_h12[50:50+7])
        h12_11 = np.mean(avg_h12[57:])

        mean_h1 = (h1_00+h1_10+h1_01+h1_11)/4
        mean_h2 = (h2_00+h2_01+h2_10+h2_11)/4
        mean_h12 = (h12_00+h12_11+h12_01+h12_10)/4

        #subtract beta 0 from all measurements
        #rescale them
        avg_h1 -= mean_h1
        avg_h2 -= mean_h2
        avg_h12 -= mean_h12

        scale_h1 = (h1_00+h1_10-h1_01-h1_11)/4
        scale_h2 = (h2_00+h2_01-h2_10-h2_11)/4
        scale_h12 = (h12_00+h12_11-h12_01-h12_10)/4

        avg_h1 = (avg_h1)/scale_h1
        avg_h2 = (avg_h2)/scale_h2
        avg_h12 = (avg_h12)/scale_h12
        #The averages have been redefined so redefine the cal terms
        h1_00 = np.mean(avg_h1[36:36+7])
        h1_01 = np.mean(avg_h1[43:43+7])
        h1_10 = np.mean(avg_h1[50:50+7])
        h1_11 = np.mean(avg_h1[57:])

        h2_00 = np.mean(avg_h2[36:36+7])
        h2_01 = np.mean(avg_h2[43:43+7])
        h2_10 = np.mean(avg_h2[50:50+7])
        h2_11 = np.mean(avg_h2[57:])

        h12_00 = np.mean(avg_h12[36:36+7])
        h12_01 = np.mean(avg_h12[43:43+7])
        h12_10 = np.mean(avg_h12[50:50+7])
        h12_11 = np.mean(avg_h12[57:])

        measurement_channel_1 = np.array([avg_h1[0],avg_h1[1],avg_h1[7],avg_h1[8],avg_h1[14],avg_h1[21],avg_h1[28],avg_h1[35]])
        measurement_channel_2 = np.array([avg_h2[0],avg_h2[1],avg_h2[7],avg_h2[8],avg_h2[14],avg_h2[21],avg_h2[28],avg_h2[35]])
        measurement_channel_3 = np.array([avg_h12[0],avg_h12[1],avg_h12[7],avg_h12[8],avg_h12[14],avg_h12[21],avg_h12[28],avg_h12[35]])
        self.measurements_tomo = np.array([measurement_channel_1,measurement_channel_2,measurement_channel_3]).flatten()

        # print(self.measurements_tomo)
        # print(len(self.measurements_tomo))

        # 108 x 1
        # get the calibration points by averaging over the five measurements
        # taken knowing the initial state we put in
        self.measurements_cal=np.array([h1_00, h1_01, h1_10, h1_11, h2_00, h2_01, h2_10, h2_11, h12_00, h12_01, h12_10, h12_11])

    def _calibrate_betas(self):
        """
        calculates betas from calibration points for the initial measurement
        operator

        Betas are ordered by B0 -> II B1 -> IZ etc(binary counting)
        <0|Z|0> = 1, <1|Z|1> = -1

        Keyword arguments:
        measurements_cal --- array(2 ** n_qubits) should be ordered
            correctly (00, 01, 10, 11) for 2 qubits
        """
        cal_matrix = np.zeros((self.n_states, self.n_states))
        # get the coefficient matrix for the betas
        for i in range(self.n_states):
            for j in range(self.n_states):
                # perform bitwise AND and count the resulting 1s
                cal_matrix[i, j] = (-1)**(bin((i & j)).count("1"))
        # invert solve the simple system of equations
        # print(cal_matrix)
        # print(np.linalg.inv(cal_matrix))
        self.betas = np.zeros(12)
        # print(self.measurements_cal[0:4])
        self.betas[0:4] = np.dot(np.linalg.inv(cal_matrix),
                            self.measurements_cal[0:4])
        self.betas_up = self.betas[0:4]
        self.betas[4:8] = np.dot(np.linalg.inv(cal_matrix),
                            self.measurements_cal[4:8])
        self.betas_p = self.betas[4:8]
        self.betas[8:] = np.dot(np.linalg.inv(cal_matrix),
                           self.measurements_cal[8:12])
        self.betas_pp = self.betas[8:]
        return self.betas

    def assemble_M_matrix_single_block(self, beta_array):
        M_matrix_single_block_row_1 = np.array([beta_array[1],
                                                beta_array[2], beta_array[3],
                                                0, 0, 0, 0, 0, 0])
        M_matrix_single_block_row_2 = np.array([-1*beta_array[1],
                                                beta_array[2],
                                                -1*beta_array[3],
                                                0, 0, 0, 0, 0, 0])
        M_matrix_single_block_row_3 = np.array([beta_array[1],
                                                -1*beta_array[2],
                                                -1*beta_array[3],
                                                0, 0, 0, 0, 0, 0])
        M_matrix_single_block_row_4 = np.array([-1*beta_array[1],
                                                -1*beta_array[2],
                                                beta_array[3],
                                                0, 0, 0, 0, 0, 0])
        M_matrix_single_block_row_5 = np.array([0, 0, 0, -beta_array[1],
                                                -beta_array[2],
                                                beta_array[3], 0, 0, 0])
        M_matrix_single_block_row_6 = np.array([0, 0, 0,
                                                beta_array[1],
                                                beta_array[2],
                                                beta_array[3],
                                                0, 0, 0])
        M_matrix_single_block_row_7 = np.array([0, 0,
                                                0, 0, 0, 0, beta_array[1],
                                                beta_array[2],
                                                beta_array[3]])
        M_matrix_single_block_row_8 = np.array([0, 0, 0, 0,
                                                0, 0, -beta_array[1],
                                                -beta_array[2],
                                                beta_array[3]])
        M_matrix_single_block = np.vstack((M_matrix_single_block_row_1,
                                           M_matrix_single_block_row_2,
                                           M_matrix_single_block_row_3,
                                           M_matrix_single_block_row_4,
                                           M_matrix_single_block_row_5,
                                           M_matrix_single_block_row_6,
                                           M_matrix_single_block_row_7,
                                           M_matrix_single_block_row_8))
        M_matrix_single_block = M_matrix_single_block.reshape(8, 9)
        return M_matrix_single_block

    def assemble_M_matrix(self):
        Block1 = self.assemble_M_matrix_single_block(self.betas_up)
        Block2 = self.assemble_M_matrix_single_block(self.betas_p)
        Block3 = self.assemble_M_matrix_single_block(self.betas_pp)
        self.M_matrix = np.vstack((Block1, Block2, Block3)).reshape(24, 9)

        return self.M_matrix

    def invert_M_matrix(self):
        self.inverse_matrix = np.linalg.pinv(self.M_matrix)
        return self.inverse_matrix

    def execute_error_signalling(self, ev):
        II = (ev[0] - ev[3])/(1 - ev[3])
        IZ = (ev[1] - ev[2])/(1 - ev[3])
        ZI = (ev[1] - ev[2])/(1 - ev[3])
        ZZ = (ev[3] - ev[0])/(1 - ev[3])
        XX = (ev[4] + ev[5])/(1 - ev[3])
        YY = (ev[4] + ev[5])/(1 - ev[3])
        ev_error_signalling = np.array([II, IZ, ZI, ZZ, XX, YY])
        return ev_error_signalling

    def execute_expectation_value_calculation_traceone(self):
        # assemble matrix that connects RO with terms
        self._calibrate_betas()
        self.assemble_M_matrix()
        self.inverse_matrix = np.linalg.pinv(self.M_matrix)

        # use it to get terms back from RO
        rescaled_measurements_tomo = self.measurements_tomo
        self.expect_values = np.dot(self.inverse_matrix,
                                    rescaled_measurements_tomo)
        expect_values_VQE = np.array([1,
                                      self.expect_values[0],
                                      self.expect_values[1],
                                      self.expect_values[2],
                                      self.expect_values[5],
                                      self.expect_values[8]])
        return expect_values_VQE

    def execute_expectation_value_calculation_T1signaling(self):
        # assemble matrix that connects RO with terms
        self._calibrate_betas()
        self.assemble_M_matrix()
        self.inverse_matrix = np.linalg.pinv(self.M_matrix[:, 1:])
        # use it to get terms back from RO
        rescaled_measurements_tomo = self.measurements_tomo
        self.expect_values = np.dot(self.inverse_matrix,
                                    rescaled_measurements_tomo)
        expect_values_VQE = np.array([1,
                                      self.expect_values[0],
                                      self.expect_values[1],
                                      self.expect_values[2],
                                      self.expect_values[5],
                                      self.expect_values[8]])
        expect_values_VQE = self.execute_error_signalling(expect_values_VQE)
        return expect_values_VQE


text:
#!/usr/bin/python
import matplotlib
matplotlib.use('Agg')
from scipy import *
from pylab import *
from numpy import *
from matplotlib.collections import LineCollection
import glob
import os
from os.path import join as pjoin

def comparacio(a,b):

	(Sepa,numa) = a.split('_')
	(pa,exta) = numa.split('.')

	(Sepb,numb) = b.split('_')
	(pb,extb) = numb.split('.')
        value=float(pa)-float (pb)
        return int(value/abs(value))

#Kv = os.getenv('K_Conn')
alpha=float(os.getenv('alpha'))
beta=float(os.getenv('beta'))
nNodes = int(os.getenv('Nodes')) # Column 0 also counts

Trans=2500
path = os.getenv('P_Dir')

font = {'family' : 'monospace',
        'weight' : 'bold',
        'size'   : 12}
rc('xtick', labelsize=11)
rc('ytick', labelsize=11)

rc('font', **font)

infiles = sorted(glob.glob( '%s/Data_*.dat' %path), cmp=comparacio)
total=zeros([401])
for infile in infiles:

	#Load data sets

	s = loadtxt(infile, unpack=True) # Voxels 1,2,3

	#Compute the auto-correlation and plot it.

	for i in range(1,nNodes+1):
		b=acorr(s[i][Trans:]-s[i][Trans:].mean(),normed=True,usevlines=False,maxlags=200,ls='-',marker='None')
		fout1=open('%s/A.dat' %(path),'w')
		fout=open('%s/A_%s.dat' %(path,i),'w')
		#xlim([-200,200])
		#savefig(pjoin(path,'Autocorrelation_%i.pdf' %i))
		close()

		for j in range(len(b[1])):
			total[j]=total[j]+b[1][j]
			print >>fout, str(b[1][j]).strip('[]')
			print >>fout1, total[j]/(nNodes)

	fout1.close()
	fout.close()
	plot(total/(nNodes))
	savefig(pjoin(path,'Autocorr_mean_alpha_%s_beta_%s.pdf' %(alpha,beta)))
	close()

ext: py

reason: from numpy,from scipy


text:
import os
import subprocess
import findspark
from pyspark import SparkContext, SparkConf
import numpy as np
import pandas as pd
from numpy import linalg as LA
from statistics import mean
import ipyleaflet
from ipyleaflet import (
    Map,
    Marker,
    TileLayer, ImageOverlay,
    Polyline, Polygon, Rectangle, Circle, CircleMarker,
    GeoJSON,
    DrawControl
)
from matplotlib.colors import rgb2hex
from matplotlib.pyplot import xticks, yticks, figure
import pylab as plt

from ipyleaflet import Heatmap, WidgetControl, FullScreenControl
from ipywidgets import IntSlider, jslink

class leaflet_old:
    """
    Plots circles on a map (one per station) whose size is proportional to the number of feature measurements, and whose color is equal to their aggregated value (min, max, avg)

    :param featureStr: the feature you'd like to visualize given as a string, eg. 'SNWD'
    :param aggregateType: 'avg', 'min', 'max'
    """

    def __init__(self, sqlctxt, featureStr):
        self.feat = featureStr
        self.pdfMaster = None
        self.sqlctxt = sqlctxt

        self.maxLong = None
        self.minLong = None
        self.maxLat = None
        self.minLat = None

        self.minAggVal = None
        self.maxAggVal = None

        self.cmap = plt.get_cmap('jet_r')
        self.m = None

    def add(self, dataframe):
        """
        Adds (or concatenates) a dataframe to the instance's master dataframe. The dataframe to be added must have atleast
        the following columns: station, latitude, longitude, featureStr, value
        """

        self.sqlctxt.registerDataFrameAsTable(dataframe, "temp")
        query = f"""
        SELECT Station, latitude, longitude, COUNT(Measurement), MEAN(Values)
        FROM temp
        WHERE Measurement=='SNWD'
        GROUP BY Station, latitude, longitude
        """
        tempdf = self.sqlctxt.sql(query).toPandas()

        if (self.pdfMaster is None):
            # set as new pdfMaster
            self.pdfMaster = tempdf

        else:
            # append to existing dfMaster
            self.pdfMaster = pd.concat([self.pdfMaster, tempdf])

    def plot_all(self):
        # update internals
        self.maxLong = self.pdfMaster['longitude'].max()
        self.minLong = self.pdfMaster['longitude'].min()
        self.maxLat = self.pdfMaster['latitude'].max()
        self.minLat = self.pdfMaster['latitude'].min()
        self.minAggVal = self.pdfMaster['avg(Values)'].min()
        self.maxAggVal = self.pdfMaster['avg(Values)'].max()

        # update center of map
        self.center = [(self.minLat + self.maxLat)/2, (self.minLong + self.maxLong)/2]
        self.zoom = 6
        self.m = Map(default_tiles=TileLayer(opacity=1.0), center=self.center, zoom=self.zoom)

        # loop over all the points in given dataframe, adding them to self.map
        circles = []
        for index,row in self.pdfMaster.iterrows():
            _lat=row['latitude']
            _long=row['longitude']
            _count=row['count(Measurement)']
            _coef=row['avg(Values)']
#             pdb.set_trace()
            # taking sqrt of count so that the  area of the circle corresponds to the count
            c = Circle(location=(_lat,_long), radius=int(1200*np.sqrt(_count+0.0)), weight=1,
                    color='#AAA', opacity=0.8, fill_opacity=0.4,
                    fill_color=self.get_color(_coef))
            circles.append(c)
            self.m.add_layer(c)
        self.m

    def color_legend(self):
        self.cfig = figure(figsize=[10,1])
        ax = plt.subplot(111)
#         vals = self.cmap(np.arange(0,1,.005))[:,:3]
        vals = self.cmap(np.arange(0,1,.005))[:,:3]
        vals3 = np.stack([vals]*10)
        vals3.shape
        ax.imshow(vals3)
        midpoint = 200. * -self.minAggVal/(self.maxAggVal - self.minAggVal)
        xticks((0,midpoint,200),["%4.1f"%v for v in (self.minAggVal,0.,self.maxAggVal)])
        yticks(());

    def get_color(self, val):
        x = (val - self.minAggVal)/(self.maxAggVal - self.minAggVal)
        return(rgb2hex(self.cmap(x)[:3]))

class leaflet_eig:
    """
    Plots circles on a map (one per station) whose size is proportional to the number of feature measurements, and whose color is equal to their aggregated value (min, max, avg)

    :param featureStr: the feature you'd like to visualize given as a string, eg. 'SNWD'
    :param aggregateType: 'avg', 'min', 'max'
    """

    def __init__(self, sqlctxt, featureStr):
        self.feat = featureStr
        self.pdfMaster = None
        self.sqlctxt = sqlctxt

        self.maxLong = None
        self.minLong = None
        self.maxLat = None
        self.minLat = None

        self.minAggVal = None
        self.maxAggVal = None

        self.cmap = plt.get_cmap('jet_r')
        self.m = None

    def add(self, dataframe):
        """
        Adds (or concatenates) a dataframe to the instance's master dataframe. The dataframe to be added must have atleast
        the following columns: station, latitude, longitude, featureStr, value
        """

        self.sqlctxt.registerDataFrameAsTable(dataframe, "temp")
        query = f"""
        SELECT Station, latitude, longitude, COUNT(Year) NumMeas, MEAN(Values)
        FROM temp
        GROUP BY Station, latitude, longitude
        """
        tempdf = self.sqlctxt.sql(query).toPandas()

        if (self.pdfMaster is None):
            # set as new pdfMaster
            self.pdfMaster = tempdf

        else:
            # append to existing dfMaster
            self.pdfMaster = pd.concat([self.pdfMaster, tempdf])

    def plot_all(self):
        # update internals
        self.maxLong = self.pdfMaster['longitude'].max()
        self.minLong = self.pdfMaster['longitude'].min()
        self.maxLat = self.pdfMaster['latitude'].max()
        self.minLat = self.pdfMaster['latitude'].min()
        self.minAggVal = self.pdfMaster['avg(Values)'].min()
        self.maxAggVal = self.pdfMaster['avg(Values)'].max()

        # update center of map
        self.center = [(self.minLat + self.maxLat)/2, (self.minLong + self.maxLong)/2]
        self.zoom = 6
        self.m = Map(default_tiles=TileLayer(opacity=1.0), center=self.center, zoom=self.zoom)

        # loop over all the points in given dataframe, adding them to self.map
        circles = []
        for index,row in self.pdfMaster.iterrows():
            _lat=row['latitude']
            _long=row['longitude']
            _count=row['NumMeas']
            _coef=row['avg(Values)']
#             pdb.set_trace()
            # taking sqrt of count so that the  area of the circle corresponds to the count
            c = Circle(location=(_lat,_long), radius=int(2000*np.sqrt(_count+0.0)), weight=1,
                    color='#AAA', opacity=0.8, fill_opacity=0.4,
                    fill_color=self.get_color(_coef))
            circles.append(c)
            self.m.add_layer(c)
        self.m.add_control(FullScreenControl())

    def color_legend(self):
        self.cfig = figure(figsize=[10,1])
        ax = plt.subplot(111)
#         vals = self.cmap(np.arange(0,1,.005))[:,:3]
        vals = self.cmap(np.arange(0,1,.005))[:,:3]
        vals3 = np.stack([vals]*10)
        vals3.shape
        ax.imshow(vals3)
        midpoint = 200. * -self.minAggVal/(self.maxAggVal - self.minAggVal)
        xticks((0,midpoint,200),["%4.1f"%v for v in (self.minAggVal,0.,self.maxAggVal)])
        yticks(());

    def get_color(self, val):
        x = (val - self.minAggVal)/(self.maxAggVal - self.minAggVal)
        return(rgb2hex(self.cmap(x)[:3]))

class leaflet:
    """
    Plots circles on a map (one per station) whose size is proportional to the number of feature measurements, and whose color is equal to their aggregated value (min, max, avg)

    :param featureStr: the feature you'd like to visualize given as a string, eg. 'SNWD'
    :param aggregateType: 'avg', 'min', 'max'
    """

    def __init__(self, sqlctxt, featureStr):
        self.feat = featureStr
        self.pdfMaster = None
        self.sqlctxt = sqlctxt

        self.maxLong = None
        self.minLong = None
        self.maxLat = None
        self.minLat = None

        self.minAggVal = None
        self.maxAggVal = None

        self.cmap = plt.get_cmap('jet')
        self.m = None

    def add(self, dataframe):
        """
        Adds (or concatenates) a dataframe to the instance's master dataframe. The dataframe to be added must have atleast
        the following columns: station, latitude, longitude, featureStr, value
        """

        self.sqlctxt.registerDataFrameAsTable(dataframe, "temp")
        query = f"""
        SELECT Station, latitude, longitude, NumMeas, MEAN(Values)
        FROM temp
        GROUP BY Station, latitude, longitude, NumMeas
        """
        tempdf = self.sqlctxt.sql(query).toPandas()

        if (self.pdfMaster is None):
            # set as new pdfMaster
            self.pdfMaster = tempdf

        else:
            # append to existing dfMaster
            self.pdfMaster = pd.concat([self.pdfMaster, tempdf])

    def plot_all(self):
        # update internals
        self.maxLong = self.pdfMaster['longitude'].max()
        self.minLong = self.pdfMaster['longitude'].min()
        self.maxLat = self.pdfMaster['latitude'].max()
        self.minLat = self.pdfMaster['latitude'].min()
        self.minAggVal = self.pdfMaster['avg(Values)'].min()
        self.maxAggVal = self.pdfMaster['avg(Values)'].max()

        # update center of map
        self.center = [(self.minLat + self.maxLat)/2, (self.minLong + self.maxLong)/2]
        self.zoom = 6
        self.m = Map(default_tiles=TileLayer(opacity=1.0), center=self.center, zoom=self.zoom)

        # loop over all the points in given dataframe, adding them to self.map
        circles = []
        for index,row in self.pdfMaster.iterrows():
            _lat=row['latitude']
            _long=row['longitude']
            _count=row['NumMeas']
            _coef=row['avg(Values)']
            # taking sqrt of count so that the  area of the circle corresponds to the count
            c = Circle(location=(_lat,_long), radius=int(1200*(_count+0.0)), weight=1,
                    color='#AAA', opacity=0.8, fill_opacity=0.4,
                    fill_color=self.get_color(_coef))
            circles.append(c)
            self.m.add_layer(c)
        self.m

    def heatmap(self):
        # update internals
        self.maxLong = self.pdfMaster['longitude'].max()
        self.minLong = self.pdfMaster['longitude'].min()
        self.maxLat = self.pdfMaster['latitude'].max()
        self.minLat = self.pdfMaster['latitude'].min()
        self.minAggVal = self.pdfMaster['avg(Values)'].min()
        self.maxAggVal = self.pdfMaster['avg(Values)'].max()

        # update center of map
        self.center = [(self.minLat + self.maxLat)/2, (self.minLong + self.maxLong)/2]
        self.zoom = 4
        self.m = Map(default_tiles=TileLayer(opacity=1.0), center=self.center, zoom=self.zoom)

        #
        locations = [[self.pdfMaster.iloc[i]['latitude'], self.pdfMaster.iloc[i]['longitude'], (self.pdfMaster.iloc[i]['avg(Values)'])*1000/(self.maxAggVal)] for i in range(len(self.pdfMaster)) if self.pdfMaster.iloc[i]['avg(Values)'] >= 0]
        hm_hot = Heatmap(locations=locations,
                     radius=10,
                     blur=5,
                     max=.8,
                     max_zoom=100,
                     gradient={0.0: 'yellow', .5: 'orange', 1.0: 'red'}
                    )
        self.m.add_layer(hm_hot)

        #
        locations = [[self.pdfMaster.iloc[i]['latitude'], self.pdfMaster.iloc[i]['longitude'], (abs(self.pdfMaster.iloc[i]['avg(Values)']))*1000/(abs(self.minAggVal))] for i in range(len(self.pdfMaster)) if self.pdfMaster.iloc[i]['avg(Values)'] <= 0]
        hm_cold = Heatmap(locations=locations,
                     radius=10,
                     blur=5,
                     max=.8,
                     max_zoom=100,
                     gradient={0.0: 'lime', 0.5: 'cyan', 1.0: 'blue'}
                    )
        self.m.add_layer(hm_cold)

        radius_slider = IntSlider(description='Saturation', min=1, max=25, value=12)
        jslink((radius_slider, 'value'), (hm_hot, 'radius'))
        jslink((radius_slider, 'value'), (hm_cold, 'radius'))
        radius_control = WidgetControl(widget=radius_slider, position='topright')
        self.m.add_control(radius_control)

        self.m.add_control(FullScreenControl())

    def color_legend(self):
        self.cfig = figure(figsize=[10,1])
        ax = plt.subplot(111)
#         vals = self.cmap(np.arange(0,1,.005))[:,:3]
        vals = self.cmap(np.arange(0,1,.005))[:,:3]
        vals3 = np.stack([vals]*10)
        vals3.shape
        ax.imshow(vals3)
        midpoint = 200. * -self.minAggVal/(self.maxAggVal - self.minAggVal)
        xticks((0,midpoint,200),["%4.1f"%v for v in (self.minAggVal,0.,self.maxAggVal)])
        yticks(());

    def get_color(self, val):
        x = (val - self.minAggVal)/(self.maxAggVal - self.minAggVal)
        return(rgb2hex(self.cmap(x)[:3]))

def decadeMeasurementDelta(featureStr, states, data_dir, sqlContext):
    deltaDFs = []
    # create the sparksql context
    for s in states:
        parquet = s + '.parquet'
        parquet_path = data_dir + '/' + parquet
        df = sqlContext.read.parquet(parquet_path)
        sqlContext.registerDataFrameAsTable(df,f'table_{s}')

        # 70s

        ###
        Query = f"""
        SELECT Station, Measurement, Values, longitude, latitude, Year
        FROM table_{s}
        WHERE Measurement=={featureStr} and (Year >= 1970 and Year < 1980)
        """
        query70s = sqlContext.sql(Query)
        rdd70s = query70s.rdd.map(lambda x: replaceSNWD(x, 'Values'))
        dfs70 = sqlContext.createDataFrame(rdd70s)
        sqlContext.registerDataFrameAsTable(dfs70, f'table_{s}_70s')

        ###
        Query = f"""
        SELECT Station, COUNT(Measurement) as NumMeas70s, MEAN(Values) as AvgVal70s, longitude, latitude
        FROM table_{s}_70s
        GROUP BY Station, Measurement, longitude, latitude
        ORDER BY Station ASC
        """
        groupbystation70s = sqlContext.sql(Query)
        sqlContext.registerDataFrameAsTable(groupbystation70s, "meas70s")
        ###

        # 00s

        ###
        Query = f"""
        SELECT Station, Measurement, Values, longitude, latitude, Year
        FROM table_{s}
        WHERE Measurement=={featureStr} and (Year >= 2000 and Year < 2010)
        """
        query00s = sqlContext.sql(Query)
        rdd00s = query00s.rdd.map(lambda x: replaceSNWD(x, 'Values'))
        dfs00 = sqlContext.createDataFrame(rdd00s)
        sqlContext.registerDataFrameAsTable(dfs00, f'table_{s}_00s')

        ###
        Query = f"""
        SELECT Station, COUNT(Measurement) as NumMeas00s, MEAN(Values) as AvgVal00s, longitude, latitude
        FROM table_{s}_00s
        GROUP BY Station, Measurement, longitude, latitude
        ORDER BY Station ASC
        """
        groupbystation00s = sqlContext.sql(Query)
        sqlContext.registerDataFrameAsTable(groupbystation00s, "meas00s")
        ###

        ### Combined

        ###
        Query = f"""
        SELECT meas70s.Station, meas70s.longitude, meas70s.latitude, AvgVal70s, AvgVal00s, NumMeas70s, NumMeas00s
        FROM meas70s
        INNER JOIN meas00s on meas70s.Station = meas00s.Station
        ORDER BY Station ASC
        """
        combined = sqlContext.sql(Query)
        rddcombined = combined.rdd.map(lambda x: computeDelta(x, 'AvgVal70s', 'AvgVal00s'))
        deltaDFs.append(sqlContext.createDataFrame(rddcombined))
        print(f"Done adding {s}")

    return(deltaDFs)

def decadeMeasurementDelta_states(featureStr, states, data_dir, sqlContext):
    deltaDFs = []
    # create the sparksql context
    for s in states:
        parquet = s + '.parquet'
        parquet_path = data_dir + '/' + parquet
        df = sqlContext.read.parquet(parquet_path)
        sqlContext.registerDataFrameAsTable(df,f'table_{s}')

        # 70s

        ###
        Query = f"""
        SELECT state, Station, Measurement, Values, longitude, latitude, Year
        FROM table_{s}
        WHERE Measurement=={featureStr} and (Year >= 1970 and Year < 1980)
        """
        query70s = sqlContext.sql(Query)
        rdd70s = query70s.rdd.map(lambda x: replaceSNWD(x, 'Values'))
        dfs70 = sqlContext.createDataFrame(rdd70s)
        sqlContext.registerDataFrameAsTable(dfs70, f'table_{s}_70s')

        ###
        Query = f"""
        SELECT state, Station, COUNT(Measurement) as NumMeas70s, MEAN(Values) as AvgVal70s, longitude, latitude
        FROM table_{s}_70s
        GROUP BY state, Station, Measurement, longitude, latitude
        ORDER BY Station ASC
        """
        groupbystation70s = sqlContext.sql(Query)
        sqlContext.registerDataFrameAsTable(groupbystation70s, "meas70s")
        ###

        # 00s

        ###
        Query = f"""
        SELECT state, Station, Measurement, Values, longitude, latitude, Year
        FROM table_{s}
        WHERE Measurement=={featureStr} and (Year >= 2000 and Year < 2010)
        """
        query00s = sqlContext.sql(Query)
        rdd00s = query00s.rdd.map(lambda x: replaceSNWD(x, 'Values'))
        dfs00 = sqlContext.createDataFrame(rdd00s)
        sqlContext.registerDataFrameAsTable(dfs00, f'table_{s}_00s')

        ###
        Query = f"""
        SELECT state, Station, COUNT(Measurement) as NumMeas00s, MEAN(Values) as AvgVal00s, longitude, latitude
        FROM table_{s}_00s
        GROUP BY state, Station, Measurement, longitude, latitude
        ORDER BY Station ASC
        """
        groupbystation00s = sqlContext.sql(Query)
        sqlContext.registerDataFrameAsTable(groupbystation00s, "meas00s")
        ###

        ### Combined

        ###
        Query = f"""
        SELECT meas70s.state, meas70s.Station, meas70s.longitude, meas70s.latitude, AvgVal70s, AvgVal00s, NumMeas70s, NumMeas00s
        FROM meas70s
        INNER JOIN meas00s on meas70s.Station = meas00s.Station
        ORDER BY Station ASC
        """
        combined = sqlContext.sql(Query)
        rddcombined = combined.rdd.map(lambda x: computeDelta(x, 'AvgVal70s', 'AvgVal00s'))
        deltaDFs.append(sqlContext.createDataFrame(rddcombined))
        print(f"Done adding {s}")

    return(deltaDFs)

def computeDelta(x, key70s, key00s):
    d = x.asDict()
    if (key70s in d) and (key00s in d):
        delta = d[key00s] - d[key70s]
        d['Values'] = delta
        d['NumMeas'] = d['NumMeas70s']+d['NumMeas00s']
    d.pop('NumMeas70s')
    d.pop('NumMeas00s')
    d.pop(key70s)
    d.pop(key00s)
    return(d)

def replaceSNWD(x, key):
    d = x.asDict()
    if key in d:
        odd_indices = list(x.Values)[1::2]
        winter_odd_indices = odd_indices[:90]
        d[key] = float(mean(winter_odd_indices))
    return(d)

def remove0sAndAverage(x, key):
    d = x.asDict()
    if key in d:
        vals = list(x.Values)
        vals = [v for v in vals if v > 0]
        if len(vals)==0:
            vals = [0]
        d[key] = float(mean(vals))
    return(d)

def getStateData(state, data_dir, tarname, parquet):

    shell = f'''
    rm -rf {data_dir}/{tarname}
    curl https://mas-dse-open.s3.amazonaws.com/Weather/by_state/{tarname} > {data_dir}/{tarname}
    ls -lh {data_dir}/{tarname}
    pwd
    cd {data_dir}
    tar -xzf {tarname}
    du ./{parquet}
    '''

    out = subprocess.getoutput(shell)
    print(out)

def getStations():
    data_dir = '../DataHW3'
    tarname = 'Weather_Stations.tgz'
    parquet = 'stations.parquet'
    shell = f'''
    curl https://mas-dse-open.s3.amazonaws.com/Weather/{tarname} > {data_dir}/{tarname}
    cd {data_dir}
    tar -xzf {tarname}
    du ./{parquet}
    '''
    out = subprocess.getoutput(shell)
    print(out)

def create_sc(pyFiles):
    sc_conf = SparkConf()
    sc_conf.setAppName("Weather_PCA")
    sc_conf.set('spark.executor.memory', '3g')
    sc_conf.set('spark.executor.cores', '1')
    sc_conf.set('spark.cores.max', '4')
    sc_conf.set('spark.default.parallelism','10')
    sc_conf.set('spark.logConf', True)
    print(sc_conf.getAll())

    sc = SparkContext(conf=sc_conf,pyFiles=pyFiles)

    return sc

def outerProduct(X):
    """Computer outer product and indicate which locations in matrix are undefined"""
    O=np.outer(X,X)
    N=1-np.isnan(O)
    return (O,N)

def sumWithNan(M1,M2):
    """Add two pairs of (matrix,count)"""
    (X1,N1)=M1
    (X2,N2)=M2
    N=N1+N2
    X=np.nansum(np.dstack((X1,X2)),axis=2)
    return (X,N)

### HW: Replace the RHS of the expressions in this function (They need to depend on S and N.
def HW_func(S,N):
    E=      np.ones([365]) # E is the sum of the vectors
    NE=     np.ones([365]) # NE is the number of not-nan antries for each coordinate of the vectors
    Mean=   np.ones([365]) # Mean is the Mean vector (ignoring nans)
    O=      np.ones([365,365]) # O is the sum of the outer products
    NO=     np.ones([365,365]) # NO is the number of non-nans in the outer product.
    return  E,NE,Mean,O,NO

def computeCov(RDDin):
    """computeCov recieves as input an RDD of np arrays, all of the same length,
    and computes the covariance matrix for that set of vectors"""
    RDD=RDDin.map(lambda v:np.array(np.insert(v,0,1),dtype=np.float64)) # insert a 1 at the beginning of each vector so that the same
                                           #calculation also yields the mean vector
    OuterRDD=RDD.map(outerProduct)   # separating the map and the reduce does not matter because of Spark uses lazy execution.
    (S,N)=OuterRDD.reduce(sumWithNan)

    E,NE,Mean,O,NO=HW_func(S,N)

    Cov=O/NO - np.outer(Mean,Mean)
    # Output also the diagnal which is the variance for each day
    Var=np.array([Cov[i,i] for i in range(Cov.shape[0])])
    return {'E':E,'NE':NE,'O':O,'NO':NO,'Cov':Cov,'Mean':Mean,'Var':Var}

def computeStatistics(sqlContext,df):
    """Compute all of the statistics for a given dataframe
    Input: sqlContext: to perform SQL queries
            df: dataframe with the fields
            Station(string), Measurement(string), Year(integer), Values (byteArray with 365 float16 numbers)
    returns: STAT, a dictionary of dictionaries. First key is measurement,
             second keys described in computeStats.STAT_Descriptions
    """

    sqlContext.registerDataFrameAsTable(df,'weather')
    STAT={}  # dictionary storing the statistics for each measurement
    measurements=['TMAX', 'SNOW', 'SNWD', 'TMIN', 'PRCP', 'TOBS']

    for meas in measurements:
        t=time()
        Query="SELECT * FROM weather\n\tWHERE measurement = '%s'"%(meas)
        mdf = sqlContext.sql(Query)
        print(meas,': shape of mdf is ',mdf.count())

        data=mdf.rdd.map(lambda row: unpackArray(row['Values'],np.float16))

        #Compute basic statistics
        STAT[meas]=computeOverAllDist(data)   # Compute the statistics

        # compute covariance matrix
        OUT=computeCov(data)

        #find PCA decomposition
        eigval,eigvec=LA.eig(OUT['Cov'])

        # collect all of the statistics in STAT[meas]
        STAT[meas]['eigval']=eigval
        STAT[meas]['eigvec']=eigvec
        STAT[meas].update(OUT)

        print('time for',meas,'is',time()-t)

    return STAT

# Compute the overall distribution of values and the distribution of the number of nan per year
def find_percentiles(SortedVals,percentile):
    L=int(len(SortedVals)/percentile)
    return SortedVals[L],SortedVals[-L]

def computeOverAllDist(rdd0):
    UnDef=np.array(rdd0.map(lambda row:sum(np.isnan(row))).sample(False,0.01).collect())
    flat=rdd0.flatMap(lambda v:list(v)).filter(lambda x: not np.isnan(x)).cache()
    count,S1,S2=flat.map(lambda x: np.float64([1,x,x**2]))\
                  .reduce(lambda x,y: x+y)
    mean=S1/count
    std=np.sqrt(S2/count-mean**2)
    Vals=flat.sample(False,0.0001).collect()
    SortedVals=np.array(sorted(Vals))
    low100,high100=find_percentiles(SortedVals,100)
    low1000,high1000=find_percentiles(SortedVals,1000)
    return {'UnDef':UnDef,\
          'mean':mean,\
          'std':std,\
          'SortedVals':SortedVals,\
          'low100':low100,\
          'high100':high100,\
          'low1000':low100,\
          'high1000':high1000
          }

# description of data returned by computeOverAllDist
STAT_Descriptions=[
('SortedVals', 'Sample of values', 'vector whose length varies between measurements'),
 ('UnDef', 'sample of number of undefs per row', 'vector whose length varies between measurements'),
 ('mean', 'mean value', ()),
 ('std', 'std', ()),
 ('low100', 'bottom 1%', ()),
 ('high100', 'top 1%', ()),
 ('low1000', 'bottom 0.1%', ()),
 ('high1000', 'top 0.1%', ()),
 ('E', 'Sum of values per day', (365,)),
 ('NE', 'count of values per day', (365,)),
 ('Mean', 'E/NE', (365,)),
 ('O', 'Sum of outer products', (365, 365)),
 ('NO', 'counts for outer products', (365, 365)),
 ('Cov', 'O/NO', (365, 365)),
 ('Var', 'The variance per day = diagonal of Cov', (365,)),
 ('eigval', 'PCA eigen-values', (365,)),
 ('eigvec', 'PCA eigen-vectors', (365, 365))
]

reason: import numpy,from numpy


text:
# prepare-data.py
#
# Generate consistent word order, word-image co-occurrence matrix,
# and randomly project image features.
# Input: data files in the format specified in the readme
#        extracted image features
# Output: sparse image-word co-occurrence matrix
#         JSON dict from word to column in the above matrix
#         randomly-projected image features

import argparse
import numpy as np
import scipy.sparse
import os
import json
from collections import defaultdict

def parse_args():
    parser = argparse.ArgumentParser()

    parser.add_argument(
        '--output_dir',
        type=str)

    parser.add_argument(
        '--vocab',
        type=str)

    parser.add_argument(
        '--doc2images',
        type=str)

    parser.add_argument(
        '--doc2words',
        type=str)

    parser.add_argument(
        '--image2row',
        type=str)

    parser.add_argument(
        '--image_features',
        type=str)

    parser.add_argument(
        '--sparse_mask',
        default=False,
        action='store_true')

    return parser.parse_args()

# randomly project image embeddings
# from a 1664-dimensional-space down to 256 dimensions
def random_projection(image_features):

    reduced_dims = 256
    new_fn = '{}_random-projection-{}.npy'.format(image_features[:-4], reduced_dims)

    print('Loading image embeddings.')
    image_embeddings = np.load(image_features, mmap_mode='r')
    print('Image embedding shape:', image_embeddings.shape)
    num_images = image_embeddings.shape[0]
    num_dims = image_embeddings.shape[1]

    print('Getting random basis vectors.')
    # V random 256-dimensional directions
    # from gaussian with 0 mean and 1 variance
    random_basis = np.zeros((num_dims, reduced_dims))
    for i in range(num_dims):
        random_basis[i,:] = np.random.randn(reduced_dims)
    print('Random basis shape (old dims x new dims): ', random_basis.shape)

    # project the documents into the lower-dimensional space
    print('Projecting.')
    embeddings_projected = np.matmul(image_embeddings, random_basis)

    print('Normalizing rows.')
    normalizer = 1.0 / np.sqrt(np.sum(embeddings_projected ** 2, axis=1))
    embeddings_projected *= normalizer[:, np.newaxis]
    print('Saving low-rank projection!')
    np.save(new_fn, embeddings_projected)

def main():

    args = parse_args()

    save_dir = args.output_dir
    if not os.path.exists(save_dir):
        os.makedirs(save_dir)

    with open(args.vocab, 'r') as f:
        vocab = set(json.load(f))
    print('Size of vocabulary: %d' % len(vocab))

    with open(args.doc2words, 'r') as f:
        doc_to_words = json.load(f)
        # remove a bunch of non-alphanumeric characters and empties
        doc_to_words = {doc: [w.replace('‘', '')
                               .replace('’', '')
                               .replace('“', '')
                               .replace('”', '')
                               .replace('…', '')
                              for w in words if w.replace('‘', '')
                                                 .replace('’', '')
                                                 .replace('“', '')
                                                 .replace('”', '')
                                                 .replace('…', '') != '']
                              for doc, words in doc_to_words.items()}
        # filter out words not in the vocabulary
        doc_to_words = {doc: [w for w in words if w in vocab]
                             for doc, words in doc_to_words.items()}
    num_tokens = sum([len(ws) for ws in doc_to_words.values()])

    with open(args.image2row, 'r') as f:
        image_name_to_row = json.load(f)

    num_images = len(image_name_to_row)
    # we need to map rows to images as well
    image_row_to_name = {row: name for name, row in image_name_to_row.items()}
    image_names_in_order = [image_row_to_name[i] for i in range(num_images)]

    print('Number of documents: {}'.format(len(doc_to_words)))
    print('Number of tokens: {}'.format(num_tokens))
    print('Number of images: {}'.format(num_images))

    with open(args.doc2images, 'r') as f:
        doc_to_images = json.load(f)
        # filter out the few images that don't have extracted features
        doc_to_images = {doc: [im for im in ims if im in image_name_to_row] for doc, ims in doc_to_images.items()}

    # filter out docs with 0 images
    docs_to_remove = []
    for doc, images in doc_to_images.items():
        if len(images) == 0:
            docs_to_remove.append(doc)
    print("Docs with no images: %d" % len(docs_to_remove))
    assert len(docs_to_remove) == 0

    unique_words = sorted(list(vocab))
    num_clusters = len(unique_words)
    word_to_index = {word: i for i, word in enumerate(unique_words)}
    print('Number of unique words / clusters = %d' % num_clusters)

    # get the words each image co-occurs with (these are fixed from documents)
    # run once! and save
    print('Getting word counts for each image as a dict.')
    image_name_to_word_counts = defaultdict(lambda: defaultdict(int))
    for doc, images in doc_to_images.items():
        this_doc_to_word_counts = defaultdict(int)
        # get counts for this doc
        for word in doc_to_words[doc]:
            this_doc_to_word_counts[word] = this_doc_to_word_counts[word] + 1
        # update all co-occurring images by merging this doc with existing counts
        for image in images:
            current = image_name_to_word_counts[image]
            for word, count in this_doc_to_word_counts.items():
                image_name_to_word_counts[image][word] = current[word] + count

    print('Constructed image membership mask:')
    image_membership_mask = np.zeros((num_images, num_clusters))
    for num, (image, word_counts) in enumerate(image_name_to_word_counts.items()):
        if num % 50000 == 0:
            print('    {}/{} images done so far.'.format(num, num_images))
        image_row = image_name_to_row[image]
        for word, count in word_counts.items():
            word_col = word_to_index[word]
            image_membership_mask[image_row][word_col] = 1

    # save everything!
    with open('{}/word-to-col.json'.format(save_dir), 'w') as f:
        json.dump(word_to_index, f, indent=4)

    # sparse matrix has a much smaller size
    # but saving the dense enables memory mapping later for much faster load times
    if args.sparse_mask:
        print('Saving sparse image membership mask.')
        sparse_mask = scipy.sparse.coo_matrix(image_membership_mask)
        scipy.sparse.save_npz('{}/image-membership-mask_sparse.npz'.format(save_dir), sparse_mask)
    else:
        print('Saving dense image membership mask.')
        np.save('{}/image-membership-mask.npy'.format(save_dir), image_membership_mask)

    # finally randomly project the image features
    random_projection(args.image_features)

if __name__ == '__main__':
    main()

reason: import numpy,import scipy


text:
"""
Test the observation process load and write
"""
import os
import shutil
import numpy as np
import flopy
from nose.tools import raises

try:
    import pymake
except:
    print('could not import pymake')

cpth = os.path.join('temp', 't041')
# delete the directory if it exists
if os.path.isdir(cpth):
    shutil.rmtree(cpth)

exe_name = 'mf2005'
v = flopy.which(exe_name)

run = True
if v is None:
    run = False

def test_hob_simple():
    """
    test041 create and run a simple MODFLOW-2005 OBS example
    """
    pth = os.path.join(cpth, 'simple')
    modelname = 'hob_simple'
    nlay, nrow, ncol = 1, 11, 11
    shape3d = (nlay, nrow, ncol)
    shape2d = (nrow, ncol)
    ib = np.ones(shape3d, dtype=np.int)
    ib[0, 0, 0] = -1
    m = flopy.modflow.Modflow(modelname=modelname, model_ws=pth,
                              verbose=False, exe_name=exe_name, )
    dis = flopy.modflow.ModflowDis(m, nlay=1, nrow=11, ncol=11, nper=2,
                                   perlen=[1, 1])

    bas = flopy.modflow.ModflowBas(m, ibound=ib, strt=10.)
    lpf = flopy.modflow.ModflowLpf(m)
    pcg = flopy.modflow.ModflowPcg(m)
    obs = flopy.modflow.HeadObservation(m, layer=0, row=5, column=5,
                                        time_series_data=[[1., 54.4],
                                                          [2., 55.2]])
    hob = flopy.modflow.ModflowHob(m, iuhobsv=51, hobdry=-9999.,
                                   obs_data=[obs])

    # Write the model input files
    m.write_input()

    # run the modflow-2005 model
    if run:
        success, buff = m.run_model(silent=False)
        assert success, 'could not run simple MODFLOW-2005 model'

    return

def test_obs_load_and_write():
    """
    test041 load and write of MODFLOW-2005 OBS example problem
    """
    pth = os.path.join('..', 'examples', 'data', 'mf2005_obs')
    opth = os.path.join(cpth, 'tc1-true', 'orig')
    # delete the directory if it exists
    if os.path.isdir(opth):
        shutil.rmtree(opth)
    os.makedirs(opth)
    # copy the original files
    files = os.listdir(pth)
    for file in files:
        src = os.path.join(pth, file)
        dst = os.path.join(opth, file)
        shutil.copyfile(src, dst)

    # load the modflow model
    mf = flopy.modflow.Modflow.load('tc1-true.nam', verbose=True,
                                    model_ws=opth, exe_name=exe_name)

    # run the modflow-2005 model
    if run:
        success, buff = mf.run_model(silent=False)
        assert success, 'could not run original MODFLOW-2005 model'

        try:
            iu = mf.hob.iuhobsv
            fpth = mf.get_output(unit=iu)
            pth0 = os.path.join(opth, fpth)
            obs0 = np.genfromtxt(pth0, skip_header=1)
        except:
            raise ValueError('could not load original HOB output file')

    npth = os.path.join(cpth, 'tc1-true', 'new')
    mf.change_model_ws(new_pth=npth, reset_external=True)

    # write the lgr model in to the new path
    mf.write_input()

    # run the modflow-2005 model
    if run:
        success, buff = mf.run_model(silent=False)
        assert success, 'could not run new MODFLOW-2005 model'

        # compare parent results
        try:
            pth1 = os.path.join(npth, fpth)
            obs1 = np.genfromtxt(pth1, skip_header=1)

            msg = 'new simulated heads are not approximately equal'
            assert np.allclose(obs0[:, 0], obs1[:, 0], atol=1e-4), msg

            msg = 'new observed heads are not approximately equal'
            assert np.allclose(obs0[:, 1], obs1[:, 1], atol=1e-4), msg
        except:
            raise ValueError('could not load new HOB output file')

def test_obs_create_and_write():
    """
    test041 create and write of MODFLOW-2005 OBS example problem
    """
    pth = os.path.join('..', 'examples', 'data', 'mf2005_obs')
    opth = os.path.join(cpth, 'create', 'orig')
    # delete the directory if it exists
    if os.path.isdir(opth):
        shutil.rmtree(opth)
    os.makedirs(opth)
    # copy the original files
    files = os.listdir(pth)
    for file in files:
        src = os.path.join(pth, file)
        dst = os.path.join(opth, file)
        shutil.copyfile(src, dst)

    # load the modflow model
    mf = flopy.modflow.Modflow.load('tc1-true.nam', verbose=True,
                                    model_ws=opth, exe_name=exe_name,
                                    forgive=False)
    # remove the existing hob package
    iuhob = mf.hob.unit_number[0]
    mf.remove_package('HOB')

    # create a new hob object
    obs_data = []

    # observation location 1
    tsd = [[1., 1.], [87163., 2.], [348649., 3.],
           [871621., 4.], [24439070., 5.], [24439072., 6.]]
    names = ['o1.1', 'o1.2', 'o1.3', 'o1.4', 'o1.5', 'o1.6']
    obs_data.append(flopy.modflow.HeadObservation(mf, layer=0, row=2, column=0,
                                                  time_series_data=tsd,
                                                  names=names, obsname='o1'))
    # observation location 2
    tsd = [[0., 126.938], [87163., 126.904], [871621., 126.382],
           [871718.5943, 115.357], [871893.7713, 112.782]]
    names = ['o2.1', 'o2.2', 'o2.3', 'o2.4', 'o2.5']
    obs_data.append(flopy.modflow.HeadObservation(mf, layer=0, row=3, column=3,
                                                  time_series_data=tsd,
                                                  names=names, obsname='o2'))
    hob = flopy.modflow.ModflowHob(mf, iuhobsv=51, obs_data=obs_data,
                                   unitnumber=iuhob)
    # write the hob file
    hob.write_file()

    # run the modflow-2005 model
    if run:
        success, buff = mf.run_model(silent=False)
        assert success, 'could not run original MODFLOW-2005 model'

        try:
            iu = mf.hob.iuhobsv
            fpth = mf.get_output(unit=iu)
            pth0 = os.path.join(opth, fpth)
            obs0 = np.genfromtxt(pth0, skip_header=1)
        except:
            raise ValueError('could not load original HOB output file')

    npth = os.path.join(cpth, 'create', 'new')
    mf.change_model_ws(new_pth=npth, reset_external=True)

    # write the model at the new path
    mf.write_input()

    # run the modflow-2005 model
    if run:
        success, buff = mf.run_model(silent=False)
        assert success, 'could not run new MODFLOW-2005 model'

        # compare parent results
        try:
            pth1 = os.path.join(npth, fpth)
            obs1 = np.genfromtxt(pth1, skip_header=1)

            msg = 'new simulated heads are not approximately equal'
            assert np.allclose(obs0[:, 0], obs1[:, 0], atol=1e-4), msg

            msg = 'new observed heads are not approximately equal'
            assert np.allclose(obs0[:, 1], obs1[:, 1], atol=1e-4), msg
        except:
            raise ValueError('could not load new HOB output file')

def test_filenames():
    """
    test041 load and run a simple MODFLOW-2005 OBS example with specified
    filenames
    """
    print('test041 load and run a simple MODFLOW-2005 OBS example with'
          ' specified filenames')
    pth = os.path.join(cpth, 'simple')
    modelname = 'hob_simple'
    pkglst = ['dis', 'bas6', 'pcg', 'lpf']
    m = flopy.modflow.Modflow.load(modelname + '.nam', model_ws=pth,
                                   check=False, load_only=pkglst,
                                   verbose=False, exe_name=exe_name,
                                   forgive=False)

    obs = flopy.modflow.HeadObservation(m, layer=0, row=5, column=5,
                                        time_series_data=[[1., 54.4],
                                                          [2., 55.2]])
    f_in = modelname + '_custom_fname.hob'
    f_out = modelname + '_custom_fname.hob.out'
    filenames = [f_in, f_out]
    hob = flopy.modflow.ModflowHob(m, iuhobsv=51, hobdry=-9999.,
                                   obs_data=[obs], options=['NOPRINT'],
                                   filenames=filenames)

    # Write the model input files
    m.write_input()

    s = 'output filename ({}) does ' \
        'not match specified name'.format(m.get_output(unit=51))
    assert m.get_output(unit=51) == f_out, s
    s = 'specified HOB input file not found'
    assert os.path.isfile(os.path.join(pth, f_in)), s

    # run the modflow-2005 model
    if run:
        success, buff = m.run_model(silent=False)
        assert success, 'could not run simple MODFLOW-2005 model'

    return

def test_multilayerhob_pr():
    """
    test041 test multilayer obs PR == 1 criteria with problematic PRs
    """
    ml = flopy.modflow.Modflow()
    dis = flopy.modflow.ModflowDis(ml, nlay=3, nrow=1, ncol=1, nper=1,
                                   perlen=[1])
    flopy.modflow.HeadObservation(ml,layer=-3,row=0,column=0,
                                  time_series_data=[[1.0, 0]],
                                  mlay={0:0.19, 1:0.69, 2:0.12})
    return

@raises(ValueError)
def test_multilayerhob_prfail():
    """
    test041 failure of multilayer obs PR == 1 criteria
    """
    ml = flopy.modflow.Modflow()
    dis = flopy.modflow.ModflowDis(ml, nlay=3, nrow=1, ncol=1, nper=1,
                                   perlen=[1])
    flopy.modflow.HeadObservation(ml,layer=-3,row=0,column=0,
                                  time_series_data=[[1.0, 0]],
                                  mlay={0:0.50, 1:0.50, 2:0.01})
    return

def test_flwob_load():
    """
    test041 create, write, and load ModflowFlwob package.
    """
    # load the modflow model
    opth = os.path.join(cpth, 'tc1-true', 'orig')
    m = flopy.modflow.Modflow.load('tc1-true.nam', verbose=True,
                                   model_ws=opth, exe_name=exe_name)

    npth = os.path.join(cpth, 'tc1-true', 'flwob')
    m.change_model_ws(new_pth=npth, reset_external=True)

    # write the lgr model in to the new path
    m.write_input()

    # add DRN package
    spd = {0: [[0, 5, 5, .5, 8e6],
               [0, 8, 8, .7, 8e6]]}
    drn = flopy.modflow.ModflowDrn(m, 53, stress_period_data=spd)

    # flow observation

    # Lists of length nqfb
    nqobfb = [1, 1]
    nqclfb = [1, 1]

    # Lists of length nqtfb
    obsnam = ['drob_1', 'drob_2']
    irefsp = [0, 0]
    toffset = [0, 0]
    flwobs = [-5.678, -6.874]

    # Lists of length (nqfb, nqclfb)
    layer = [[0], [0]]
    row = [[5], [8]]
    column = [[5], [8]]
    factor = [[1.], [1.]]

    drob = flopy.modflow.ModflowFlwob(m,
                                      nqfb=len(nqclfb),
                                      nqcfb=np.sum(nqclfb),
                                      nqtfb=np.sum(nqobfb),
                                      nqobfb=nqobfb,
                                      nqclfb=nqclfb,
                                      obsnam=obsnam,
                                      irefsp=irefsp,
                                      toffset=toffset,
                                      flwobs=flwobs,
                                      layer=layer,
                                      row=row,
                                      column=column,
                                      factor=factor,
                                      flowtype='drn',
                                      options=['NOPRINT'])
    # Write the model input files
    m.write_input(check=False)

    # Load the DROB package and compare it to the original
    pkglst = ['drob']
    m = flopy.modflow.Modflow.load('tc1-true.nam', model_ws=npth,
                                   check=False, load_only=pkglst,
                                   verbose=False, exe_name=exe_name,
                                   forgive=False)

    # check variables were read properly
    s = 'nqfb loaded from {} read incorrectly'.format(m.drob.fn_path)
    assert(drob.nqfb == m.drob.nqfb), s
    s = 'nqcfb loaded from {} read incorrectly'.format(m.drob.fn_path)
    assert (drob.nqcfb == m.drob.nqcfb), s
    s = 'nqtfb loaded from {} read incorrectly'.format(m.drob.fn_path)
    assert (drob.nqtfb == m.drob.nqtfb), s
    s = 'obsnam loaded from {} read incorrectly'.format(m.drob.fn_path)
    assert (list([n for n in drob.obsnam]) ==
            list([n for n in m.drob.obsnam])), s
    s = 'flwobs loaded from {} read incorrectly'.format(m.drob.fn_path)
    assert np.array_equal(drob.flwobs, m.drob.flwobs), s
    s = 'layer loaded from {} read incorrectly'.format(m.drob.fn_path)
    assert np.array_equal(drob.layer, m.drob.layer), s
    s = 'row loaded from {} read incorrectly'.format(m.drob.fn_path)
    assert np.array_equal(drob.row, m.drob.row), s
    s = 'column loaded from {} read incorrectly'.format(m.drob.fn_path)
    assert np.array_equal(drob.column, m.drob.column), s

    return

if __name__ == '__main__':
    test_hob_simple()
    test_obs_create_and_write()
    test_obs_load_and_write()
    test_filenames()
    test_flwob_load()

max_stars_repo_name: geowat/flopy

max_issues_repo_name: geowat/flopy

max_forks_repo_name: geowat/flopy

reason: import numpy


text:
#!/usr/bin/env python

import rospy
from std_msgs.msg import Int32
from geometry_msgs.msg import PoseStamped, TwistStamped
from styx_msgs.msg import Lane, Waypoint
import numpy as np
from scipy.spatial import KDTree
import math

#EN: at waypotin 6392 area, lag occurs, and then after stop light vehcile stops updating position...
# dropping rospy refresh rate down to 30 might be a fix

LOOKAHEAD_WPS = 200 # Number of waypoints we will publish. You can change this number

class WaypointUpdater(object):
    def __init__(self):
        rospy.init_node('waypoint_updater')

        rospy.Subscriber('/current_pose', PoseStamped, self.pose_cb)
        rospy.Subscriber('/base_waypoints', Lane, self.waypoints_cb)
        rospy.Subscriber('/current_velocity',TwistStamped, self.velocity_cb)

        rospy.Subscriber('/traffic_waypoint', Int32, self.traffic_cb)

        self.final_waypoints_pub = rospy.Publisher('final_waypoints', Lane, queue_size=1)
        self.current_waypoint_idx = rospy.Publisher('current_waypoint_idx', Int32, queue_size=1)

        # TODO: Add other member variables you need below
        self.decel_limit = rospy.get_param('~decel_limit', -5)
        self.pose = None
        self.waypoints = None
        self.waypoints_2d = None
        self.waypoint_tree = None
        self.light_idx = None
        self.target_speed = 25.755836 # in mph
        self.current_vel = None

        self.loop()

    def loop(self):
        rate = rospy.Rate(30)
        while not rospy.is_shutdown():
            rate.sleep()
            lane = Lane()
            if self.waypoints and self.pose and self.current_vel:
                closet_waypoint_index = self.get_closest_waypoint_idx()
                wp_length = len(self.waypoints.waypoints)
                looked_ahead_index = closet_waypoint_index+LOOKAHEAD_WPS
                wrap_around = looked_ahead_index - wp_length

                # generate default target speed profile:
                lane.waypoints = self.waypoints.waypoints[closet_waypoint_index:looked_ahead_index]

                if wrap_around > 0:
                    lane.waypoints = np.append(lane.waypoints, self.waypoints.waypoints[0:wrap_around])
                # set default vehicle speed:
                for waypoint in lane.waypoints:
                    waypoint.twist.twist.linear.x = self.target_speed

                # Check for light status:
                if self.light_idx != -1 and self.light_idx >= closet_waypoint_index and (self.light_idx <=looked_ahead_index or (wrap_around>0 and  self.light_idx <= wrap_around)):
                    # number of stopping points:
                    wp2stop = self.light_idx - closet_waypoint_index
                    # check the current speed and see how many waypoints are needed to come to 0 speed...
                    minWaypoints2stop = abs(self.current_vel/0.8*self.decel_limit)
                    if minWaypoints2stop > wp2stop:
                        # this will attempt to break hard, but possibly fail.
                        for waypoint in lane.waypoints[0:self.light_idx - closet_waypoint_index]:
                            waypoint.twist.twist.linear.x = 0
                    else:
                        # set last 5 waypoints to 0:
                        end_idx = self.light_idx - closet_waypoint_index
                        start_idx = end_idx-5
                        if start_idx < 0:
                            start_idx = 0
                        for waypoint in lane.waypoints[start_idx:end_idx]:
                            waypoint.twist.twist.linear.x = 0
                        # smooth deceleration up to last 5 waypoints:
                        if start_idx > 0:
                            for i in range(1, start_idx):
                                lane.waypoints[i].twist.twist.linear.x = lane.waypoints[i-1].twist.twist.linear.x - self.decel_limit

                self.final_waypoints_pub.publish(lane)
                self.current_waypoint_idx.publish(closet_waypoint_index)

    def get_closest_waypoint_idx(self):
        # short circuit out of here:
        if self.waypoint_tree is None:
            return -1
        x = self.pose.pose.position.x
        y = self.pose.pose.position.y
        closest_idx = self.waypoint_tree.query([x,y],1)[1]

        #check if closest is ahead of behind the vehicle:
        closest_coord = self.waypoints_2d[closest_idx]
        prev_coord = self.waypoints_2d[closest_idx-1]

        #hyperplane:
        cl_vect = np.array(closest_coord)
        prev_vect = np.array(prev_coord)
        pos_vect = np.array([x,y])

        val = np.dot(cl_vect - prev_vect, pos_vect - cl_vect)

        if val > 0:
            closest_idx = (closest_idx + 1) % len(self.waypoints_2d)

        return closest_idx

    def velocity_cb(self, msg):
        self.current_vel = msg.twist.linear.x

    def pose_cb(self, msg):
        self.pose = msg

    def waypoints_cb(self, waypoints):
        self.waypoints = waypoints
        if not self.waypoints_2d:
                self.waypoints_2d = [[waypoint.pose.pose.position.x, waypoint.pose.pose.position.y] for waypoint in waypoints.waypoints]
                self.waypoint_tree = KDTree(self.waypoints_2d)

    def traffic_cb(self, msg):
        self.light_idx = msg.data

    def obstacle_cb(self, msg):
        pass

    def get_waypoint_velocity(self, waypoint):
        return waypoint.twist.twist.linear.x

    def set_waypoint_velocity(self, waypoints, waypoint, velocity):
        waypoints[waypoint].twist.twist.linear.x = velocity

    def distance(self, waypoints, wp1, wp2):
        dist = 0
        dl = lambda a, b: math.sqrt((a.x-b.x)**2 + (a.y-b.y)**2  + (a.z-b.z)**2)
        for i in range(wp1, wp2+1):
            dist += dl(waypoints[wp1].pose.pose.position, waypoints[i].pose.pose.position)
            wp1 = i
        return dist

if __name__ == '__main__':
    try:
        WaypointUpdater()
    except rospy.ROSInterruptException:
        rospy.logerr('Could not start waypoint updater node.')

max_stars_repo_name: evgenyslab/CarND-Capstone

max_issues_repo_name: evgenyslab/CarND-Capstone

max_forks_repo_name: evgenyslab/CarND-Capstone

reason: import numpy,from scipy


text:
import numpy as np

def rmse(y, t):
	"""
	Returns the Root Mean Square Error of predicted values y
	from target values t.
	Args:
		y: numpy array of size (samples, features).
		t: numpy array of size (samples, features).
	Returns:
		float, the RMSE.
	"""
	N = y.shape[0]
	error_norms = np.linalg.norm(y - t, axis=1);
	return np.sqrt(np.dot(error_norms, error_norms)/N)


text:
from numpy import linspace

# version
VERSION = '0.3.0'

# constants
IDEAL_GAS_CONSTANT_KCAL = 1.987204118E-3
TEMPERATURE_CELSIUS = 25.

# calculated constants
ZERO_KELVIN = 273.15
TEMPERATURE_KELVIN = ZERO_KELVIN + TEMPERATURE_CELSIUS
RT = IDEAL_GAS_CONSTANT_KCAL * TEMPERATURE_KELVIN

# error
FITTING_PENALTY = 1e10

# display properties
FONT_SIZE = 18
MARKER_SIZE = 8
LABEL_SIZE = 18
LINE_WIDTH = 3

# error options
CONFIDENCE_INTERVAL = 95

# default range to simulate fit function
XSIM = linspace(0.,10.,100)

if __name__ == "__main__":
	pass

max_stars_repo_name: quantumjot/PyFolding

max_issues_repo_name: quantumjot/PyFolding

max_forks_repo_name: quantumjot/PyFolding

reason: from numpy


text:
#-------------------------------------------------------------------------
# KJUN
from numpy import *
import gzip, numpy, numpy as np, scipy.linalg as sla, os, operator, traceback, sys, ipdb
import numpy.random as ra, ipdb, pickle as pickle, time, numpy.linalg as la, scipy.stats as st;
import itertools          # For Multiprobe
import copy
from datetime import datetime
from collections import OrderedDict
#import iprod
#from distcolors import get_distinguishable_colors

np.set_printoptions(precision=4);

################################################################################
# Pickle
################################################################################
def LoadPickle(fName):
    """ load a pickle file. Assumes that it has one dictionary object that points to
 many other variables."""
    if type(fName) == str:
        try:
            fp = open(fName, 'rb')
        except:
            print("Couldn't open %s" % (fName))
            traceback.print_exc(file=sys.stderr)
            ipdb.set_trace();
    else:
        fp = fName
    try:
        ind = pickle.load(fp)
        fp.close()
        return ind
    except:
        print("Couldn't read the pickle file", fName)
        traceback.print_exc(file=sys.stderr)
        ipdb.set_trace();

def SavePickle(filename, var, protocol=2):
    try:
        with open(filename, 'wb') as f:
            pickle.dump(var, f, protocol=protocol)
        statinfo = os.stat(filename,)
        if statinfo:
            print("Wrote out", statinfo.st_size, "bytes to", \
                filename)
    except:
        print("Couldn't pickle the file", filename)
        traceback.print_exc(file=sys.stderr)
        ipdb.set_trace();

def LoadPickleGzip(fName):
    if type(fName) == str:
        try:
            fp = gzip.GzipFile(fName, 'rb')
        except:
            print("Couldn't open %s" % (fName))
            return None
    else:
        fp = fName
    try:
        ind = pickle.load(fp)
        fp.close()
        return ind
    except:
        print("Couldn't read pickle file", fName)
        traceback.print_exc(file=sys.stderr)

def SavePickleGzip(filename, var, protocol=2):
    try:
        with gzip.GzipFile(filename, 'wb') as f:
            pickle.dump(var, f, protocol=protocol)
        statinfo = os.stat(filename,)
        if statinfo:
            print("Wrote out", statinfo.st_size, "bytes to", \
                filename)
    except:
        print("Couldn't pickle index to file", filename)
        traceback.print_exc(file=sys.stderr)

# saves variables in 'varList' to filename, and it searches variables from given 'dic'
# typically, if you want to save variables from current python context
# set 'locals()' as parameter 'dic'
def savePickleFromDic(varList, fileName, dic):
    varDic = {};
    for k in varList:
        varDic[k] = dic[k];
    f = open(fileName, 'wb');
    pickle.dump(varDic, f);
    f.close();

################################################################################
# other utility function
################################################################################

def SaveToDict(dic, varList):
    varDic = {};
    for k in varList:
        varDic[k] = dic[k];
    return varDic;

from types import SimpleNamespace
def SaveToSimpleNamespace(dic, varList):
    ret = SimpleNamespace()
    for k in varList:
        varDic[k] = dic[k];
    return varDic;

    import inspect
    frame = inspect.currentframe()
    try:
        loc = frame.f_back.f_locals
        print(expr, '= ', end=' ')
        if (bPretty):
            pprint(eval(expr, globals(), loc));
        else:
            print((eval(expr, globals(), loc)));
    finally:
        del frame

#    locals(), ['opts', 'res'])

def importVarsFromDict(srcDict, destDict):
    for (k,v) in srcDict.items():
        assert k not in destDict;
        destDict[k] = v;
    pass

def tic():
    """
    equivalent to Matlab's tic. It start measuring time.
    returns handle of the time start point.
    """
    global gStartTime
    gStartTime = datetime.utcnow();
    return gStartTime

def toc(prev=None):
    """
    get a timestamp in seconds. Time interval is from previous call of tic() to current call of toc().
    You can optionally specify the handle of the time ending point.
    """
    if prev==None: prev = gStartTime;
    return (datetime.utcnow() - prev).total_seconds();

def printExpr(expr, bPretty=True):
    """ Print the local variables in the caller's frame."""
    from pprint import pprint
    import inspect
    frame = inspect.currentframe()
    try:
        loc = frame.f_back.f_locals
        print(expr, '= ', end=' ')
        if (bPretty):
            pprint(eval(expr, globals(), loc));
        else:
            print((eval(expr, globals(), loc)));
    finally:
        del frame

def DictInvert(aDict):
    """\\
    Inverts a dictionary; {(key,value)} is turned into {(value, [key1, key2, ...])}.
    For example, {(1,'a'), (2,'b'), (3,'a')} is turned into {('a',[1,3]), ('b',2)}
    """
    invDict = {};
    for k,v in aDict.items():
        invDict[v] = invDict.get(v, []);
        invDict[v].append(k);

    return invDict;

################################################################################
# numpy
################################################################################

def calcListStat(aList):
    res = Struct()
    res.nItems = len(aList);
    res.minVal = min(aList);
    res.maxVal = max(aList);
    res.mean = np.mean(aList);
    res.median = np.median(aList);
    sortedList = sorted(aList)

    idx = max([int(round(float(len(aList)) * 0.95)) - 1, 0]);
    res.upper5Perc = sortedList[idx];
    idx = max([int(round(float(len(aList)) * 0.05)) - 1, 0]);
    res.lower5Perc = sortedList[idx];
    return res;

def listFindAll(searchList, elem):
    return [[i for i, x in enumerate(searchList) if x == e] for e in elem]

def fillNFromKFairly(N,K,randStream=ra):
    """
    In [28]: res = fillNFromKFairly(10,3)
    In [29]: res
    Out[29]: array([1, 2, 0, 1, 2, 0, 1, 2, 0, 1])
    """
    perm = ra.permutation(K);
    res = np.zeros(N);
    res = np.remainder(np.arange(N), K);
    res = perm[res];
    return res;

def ListOf2dArrayTo3d(mat):
    tmp = np.zeros((len(mat), mat[0].shape[0], mat[0].shape[1]));
    for i in range(len(mat)): tmp[i,:,:] = mat[i]
    return tmp;

def argkmax(npary,k):
    return npary.argsort()[-k:][::-1]

def get_time_now_kwang():
    return time.strftime('%Y%m%d%a-%H%M%S');

def mahalanobis_norm_sq(x, M):
    assert(x.ndim == 1);
    return np.dot(x, np.dot(M,x));

def mahalanobis_norm_sq_batch_old(X, M):
    """ X is N by d, M is d by d """
    return (np.dot(X, M) * X).sum(1)

def mahalanobis_norm_sq_batch(X, M):
    """ X is N by d, M is d by d """
    return np.einsum('...i,...i', np.dot(X, M), X);

def qoful_construct_query(q_vector, q_mat):
    assert(q_vector.ndim == 1);
    d = len(q_vector);
    q = np.zeros( (d + d**2,1) );
    q[:d,0] = q_vector;
    q[d:,0] = q_mat.ravel();
    return q;

def chooseTopIdxList(ary, topRatio):
    """
      choose topRatio largest member (returns their indexes)
    """
    sidx = np.argsort(ary);
    topIdx = int(np.round(len(ary)*(1-topRatio)));
    return sidx[topIdx:]

def chooseInitPoint(expectedRewardAry, rand_stream=ra):
    """
      choose the initial point by (1) draw rewards (2) choose one randomly
    """
    n = len(expectedRewardAry);
    rewardAry = rand_stream.rand(n) < expectedRewardAry;
    if np.sum(rewardAry) == 0:
      idxAry = np.arange(n);  # when no positive, just choose anything.
    else:
      idxAry = np.where(rewardAry)[0];
    return rand_stream.choice(idxAry);

def chooseInitPoint_v2(expectedRewardAry, nTry=10, rand_stream=ra):
    """
      choose the initial point by (1) draw rewards 10 times (2) choose one randomly
    """
    n = len(expectedRewardAry);
    sumRewardAry = np.zeros(n);
    for iTry in range(nTry):
        rewardAry = rand_stream.rand(n) < expectedRewardAry;
        sumRewardAry += rewardAry;

    if np.sum(rewardAry) == 0:
        idxAry = np.arange(n);  # when no positive, just choose anything.
    else:
        idxAry = np.where(np.max(sumRewardAry) == sumRewardAry)[0]
    return rand_stream.choice(idxAry);

def calcRmse(yhat, y):
    assert yhat.shape == y.shape;
    return np.sqrt(np.mean((yhat - y)**2));

def cmapGen(cmapName='copper', minVal=0.0, maxVal=1.0):
    import matplotlib.pyplot as plt;
    import matplotlib.colors as colors
    import matplotlib.cm as cmx
    myCmap = cm = plt.get_cmap(cmapName);
    cNorm  = colors.Normalize(vmin=minVal, vmax=maxVal)
    scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=myCmap);
    return scalarMap;

def cmapGetColor(scalarMap, v):
    return scalarMap.to_rgba(v);

def range_ensure_endpoint(a, b, step):
    """
      In [15]: range_ensure_endpoint(0,10,2)
      Out[15]: [0, 2, 4, 6, 8, 9]

      In [16]: range_ensure_endpoint(0,9,2)
      Out[16]: [0, 2, 4, 6, 8]

      In [18]: range(0,9,2)
      Out[18]: [0, 2, 4, 6, 8]

      In [19]: range(0,10,2)
      Out[19]: [0, 2, 4, 6, 8]
    """
    aList = list(range(a,b,step));
    if (aList[-1] != b-1):
        aList += [b-1];
    return aList;

def translateIndex(ia, ib):
    """ #- ia and ib can be string lists.
        ia = np.array([2,4,5,1,0,3]);
        ib = np.array([5,3,4,2,0,1]);
        b2a,a2b = translateIndex(ia,ib);
        assert np.all(ia == ib[b2a])
    """
    sa = np.argsort(ia);
    sb = np.argsort(ib);
    inv_sa = np.zeros_like( sa );
    inv_sa[sa] = np.arange(len(ia));

    b2a = sb[inv_sa];
    a2b = np.zeros_like( sa );

    a2b[b2a] = np.arange(len(ia));
    return b2a, a2b

def gen_seeds(seed, n_seeds):
    ra.seed(seed);
    return ra.randint(0,np.iinfo(np.uint32).max+1, n_seeds);

def kjunSeed(baseSeed, i):
    """
        generates i-th seeds
        Due to the way it generates the seed, do not use i that is too large..
    """
    assert i <= 100000
    rs = ra.RandomState(baseSeed);
    randVals = rs.randint(np.iinfo(np.uint32).max+1, size=i+1);
    return randVals[-1];

def kjunSeedList(baseSeed, n):
    """
        generates n seeds
        Due to the way it generates the seed, do not use i that is too large..
    """
    assert n <= 100000
    rs = ra.RandomState(baseSeed);
    randVals = rs.randint(np.iinfo(np.uint32).max+1, size=n);
    return randVals;

def nans(*args):
    ary = np.zeros(*args);
    ary.fill(np.nan);
    return ary;

def dstack_product(x, y):
    return np.dstack(numpy.meshgrid(x, y)).reshape(-1, 2)

################################################################################
# Quadratic optimization
################################################################################

def projectOntoL1Ball(v, b):
    """
    PROJECTONTOL1BALL Projects point onto L1 ball of specified radius.

    w = ProjectOntoL1Ball(v, b) returns the vector w which is the solution
    to the following constrained minimization problem:

     min   ||w - v||_2
     s.t.  ||w||_1 <= b.

    That is, performs Euclidean projection of v to the 1-norm ball of radius
    b.

    Author: John Duchi (jduchi@cs.berkeley.edu)
    """
    assert b >= 0, 'Radius of L1 ball is negative';
    assert v.dtype == np.float64;
    if (np.linalg.norm(v,1) < b):
        return v;
    u = np.abs(v);    # this makes a copy
    u[::-1].sort();   # reverse sorting (in-place sort)
    sv = np.cumsum(u);
    rho = np.where(u > ((sv-b) / np.arange(1,len(u)+1,dtype=np.float64)))[0][-1];
    th = np.maximum(0.0, (sv[rho] - b) / (rho+1));
    w = np.sign(v) * np.maximum(abs(v) - th, 0.0);
    return w;

class QuadOptimData(object):
    def __init__(self, dim):
        self.dim = dim;
        self.A = np.zeros((dim,dim));
        self.b = np.zeros(dim);

    def set_A(self, A):
        self.A[:,:] = A;

    def set_b(self, b):
        self.b[:] = b;

def objQuad(th, data, opt):
    f = None;
    g = None;
    v = dot(data.A,th);
    if (opt == 1 or opt == 3):
        f = .5*dot(th,v) + dot(th,data.b); #data.x) - dot(data.thHat,v);
    if (opt == 2 or opt == 3):
        g = v + data.b;
    return f,g;

def minFuncQuadL1Options():
    ret = dict();
    ret['debug'] = False;
    ret['maxIter'] = 400;
    ret['maxLineSearch'] = 10;
    ret['tolX'] = 1e-7;
    ret['tolObj'] = 1e-7;
    ret['alpha0'] = .01; # 2.0/maxEigVal(data.A) is a good heuristic.
    ret['line_c'] = 1e-3; # line_* is for line search
    ret['line_tau'] = .1;
    ret['line_tau0'] = 1.5; # factor to multiply before line search

    return ret;

def minFuncQuadL1(qoData, maxL1Norm, th0, opt):
    myObj = lambda th, opt: objQuad(th, qoData, opt);

    debug = opt['debug'];
    th = th0;
    if (debug):
        nLineSearchAry = [];
        objValAry = [];
        alphaAry = [];
    maxLineSearch = opt['maxLineSearch'];
    tolX = opt['tolX'];
    tolObj = opt['tolObj'];
    alphaOld = opt['alpha0'];
    line_c = opt['line_c'];
    line_tau = opt['line_tau'];
    line_tau0 = opt['line_tau0'];
    fCnt = 0;
    gCnt = 0;
    bConvergedFVal = False;
    bConvergedX = False;
    #tic();
    f, trash = myObj(th,1);
    for k in range(opt['maxIter']):
        trash, g = myObj(th,3);
        p = -g;
#        p = -np.dot(qoData.invA, g);
        gCnt += 1;

        #- perform backtracking
        alpha = line_tau0 * alphaOld;
        thCur = th;
        for cnt in range(1,maxLineSearch+1):
#         cnt = 1;
#         while cnt <= maxLineSearch:
            thNew = projectOntoL1Ball(thCur + alpha*p, maxL1Norm);
            fNew, trash = myObj(thNew, 1);
            fCnt += 1;
            if (fNew < f - alpha*line_c):
                break;
            elif la.norm(thNew - thCur) / (1+la.norm(thCur)) < tolX:
                if (fNew > f):
                    thNew = thCur;
                    fNew = f;
                break;
            alpha *= line_tau;
#            cnt += 1;

        if (debug):
            nLineSearchAry.append(cnt);
            objValAry.append(f);
            alphaAry.append(alpha);

        if (abs(f - fNew) / (1+abs(f)) <= tolObj):
            bConvergedFVal = True;
            break;
        elif la.norm(thNew - th) / (1 + la.norm(th)) < tolX:
            bConvergedX = True;
            break;

        th = thNew;
        f = fNew;
        alphaOld = alpha;
    #- thNew: the solution

    #printExpr('toc()');
    iterCnt = k + 1;

    info = {'iterCnt':iterCnt, 'fCnt':fCnt, 'gCnt':gCnt, 'bConverged': bConvergedFVal or bConvergedX};

    debugDict = None;
    if (debug):
        info['debugDict'] = {'nLineSearchAry':nLineSearchAry, 'objValAry':objValAry, 'alphaAry':alphaAry };
        return thNew, fNew, info;
    else:
        return thNew, fNew, info;

################################################################################
#- for error bound / ml evaluation
################################################################################

def confidenceFactor(N, alpha=0.05):
    return st.t.isf(alpha/2.0, N-1)

def getDeviation(ary, alpha=0.05):
    N = len(ary);
    return confidenceFactor(N, alpha=alpha) * (ary.std(ddof=1) / np.sqrt(N));

def getDeviationMat(mat, alpha=0.05):
    """
    ary: (nMethod) by (nTry)
    """
    nMethod, nTry = mat.shape;
    mat = np.array(mat);
    std = mat.std(1, ddof=1);
    return confidenceFactor(nTry, alpha=alpha) * std / np.sqrt(nTry);

def getErrorBar(ary, alpha=0.05):
    ary = np.array(ary);
    me = ary.mean();
    st = ary.std(ddof=1);
    dev = st / np.sqrt(len(ary)) * confidenceFactor(len(ary), alpha);
    return me, dev

def getErrorBarMat(mat, alpha=0.05):
    me = np.mean(np.array(mat),1);
    return me, getDeviationMat(mat, alpha);

def evalTestSignificance(aAry, bAry, alpha=0.05):
    diffAry = aAry - bAry;

    n = len(aAry);
    threshold = st.t.isf(alpha/2.0, n-1);
    me = diffAry.mean();
    std = diffAry.std(ddof=1);
    x = diffAry - me;

    dev = np.abs(me / (std / np.sqrt(n)));
    if (dev == np.nan):
        tf = false;
    else:
        tf = (dev > threshold);

    return tf, dev, threshold;

def evalTestSignificanceMat(scoreMat, alpha=0.05, maxOrMin='max'):
    assert maxOrMin in ['min', 'max'];
    nMethod,nTry = scoreMat.shape;

    me = np.mean(scoreMat,1);
    # st = np.std(scoreMat,1);

    significanceMat = np.full( (nMethod,nMethod), False, dtype=bool);
    # devAry = evalCalcDeviation(scoreMat, alpha=alpha)';

    for i in range(nMethod-1):
        for j in range(i+1,nMethod):
            tf, dev, thres = evalTestSignificance(scoreMat[i,:], scoreMat[j,:], alpha=alpha);
            significanceMat[i,j] = tf;
            significanceMat[j,i] = tf;

    if (maxOrMin == 'max'):
        bestIdx = me.argmax()
    else:
        bestIdx = me.argmin()

    indistinctBestAry = np.where(significanceMat[bestIdx,:] == False)[0];

    resDict = {};
    resDict['mean'] = me;
    resDict['alpha'] = alpha;
    resDict['significanceMat'] = significanceMat;
    return bestIdx, indistinctBestAry, resDict;

def nullspace(A, atol=1e-13, rtol=0):
    """Compute an approximate basis for the nullspace of A.

    The algorithm used by this function is based on the singular value
    decomposition of `A`.

    Parameters
    ----------
    A : ndarray
        A should be at most 2-D.  A 1-D array with length k will be treated
        as a 2-D with shape (1, k)
    atol : float
        The absolute tolerance for a zero singular value.  Singular values
        smaller than `atol` are considered to be zero.
    rtol : float
        The relative tolerance.  Singular values less than rtol*smax are
        considered to be zero, where smax is the largest singular value.

    If both `atol` and `rtol` are positive, the combined tolerance is the
    maximum of the two; that is::
        tol = max(atol, rtol * smax)
    Singular values smaller than `tol` are considered to be zero.

    Return value
    ------------
    ns : ndarray
        If `A` is an array with shape (m, k), then `ns` will be an array
        with shape (k, n), where n is the estimated dimension of the
        nullspace of `A`.  The columns of `ns` are a basis for the
        nullspace; each element in numpy.dot(A, ns) will be approximately
        zero.
    """

    A = np.atleast_2d(A)
    u, s, vh = la.svd(A)
    tol = max(atol, rtol * s[0])
    nnz = (s >= tol).sum()
    ns = vh[nnz:].conj().T
    return ns

#------------- stdout redirector
from contextlib import contextmanager

@contextmanager
def stdout_redirector(stream):
    old_stdout = sys.stdout
    sys.stdout = stream
    try:
        yield
    finally:
        sys.stdout = old_stdout

max_stars_repo_name: jajajang/efalb

max_issues_repo_name: jajajang/efalb

max_forks_repo_name: jajajang/efalb