text:
'''
Collections of several functions that used to data pre-processing.
'''

__all__=["detect_operation","min_max_scale","min_max_recover","filter_na_outlier","scale_by_name","detect_defrost","nan_mean","create_unit_input","create_input_values","get_scale_const_template","get_last_outputs","get_wetbulb"]
import numpy as np
import pandas as pd
import re
import pickle
#from metpy.calc import dewpoint_from_relative_humidity,wet_bulb_temperature
#from metpy.units import units
import psypy.psySI as SI

def psypy_wb(temp,rh,p_atm=101325):
    '''Calculate Wet-bulb temperature from dry-bulb temperature and relative humidity via psypy package.
    Args:
        temp (numeric): temperature in celsius [C].
        rh (numeric): relative humidity in 0-1 scale [-].
        p_atm (numeric) : atmospheric pressure [Pa]. Default is 101325 Pa.
    Raises: return np.nan when wet-bulb is too low.
    Returns:
        wb (numeric): wet-bulb temperature in celsius.
    '''
    temp_k=temp+273.15 # Kelvin
    try:
        S=SI.state("DBT",temp_k,"RH",rh,p_atm)
        wb=S[5]-273.15 # Kelvin to Celsius
    except:
        wb=np.nan # return np.nan when wet-bulb is too low.
    return wb

def get_wetbulb(temp,rh,p_atm=101325):
    '''Vectorization of wet-bulb calculation (psypy_wb function) while skipping nan.
    Args:
        temp (np.array [1d]): temperature in celsius [C].
        rh (np.array [1d]): relative humidity in 0-1 scale [-].
        p_atm (np.array [1d] or numeric) : atmospheric pressure [Pa]. Default is 101325 Pa.
    Raises:
    Returns:
        wb (np.array [1d]): wet-bulb temperature in np.array
    '''
    nan_index=np.any(np.isnan(np.stack([temp,rh],axis=1)),axis=1) # find nan data in row
    base_=np.zeros([temp.shape[0]]) # original data dimension
    base_[nan_index]=np.nan # put nan for nan values.
    n_data=temp[~nan_index].shape[0] # non-nan data
    vec_psypy_wb=np.vectorize(psypy_wb) # vectorize psypy_wb function
    wb=vec_psypy_wb(temp=temp[~nan_index],rh=rh[~nan_index],p_atm=p_atm) # calculate wet-bulb for non-nan data
    base_[~nan_index]=wb # put calculated wet-bulb of non-nan data for the original data dimension
    return(base_)

def min_max_scale(x,x_min,x_max,lower=0.,upper=1.):
    '''
    min_max_scale (see rescaling in https://en.wikipedia.org/wiki/Feature_scaling).
    The data in [x_min,x_max] is linearly rescaled to [lower,upper].
    # it can hanldes np.nan.
    Args:
        x (1d np.array): data that needs to be rescaled.
        x_min (numeric): minimum value of x.
        x_max (numeric): maximum value of x.
        lower (numeric): minimum value of rescaled x. default is 0.
        upper (numeric): maximum value of rescaled x. default is 1.
    Raises:
    Returns:
        rescaled x (numeric): rescaled x from [x_min,x_max] to [lower,upper]
    '''
    if type(x)!=np.ndarray:
        raise ValueError("Put data as numpy ndarray.format.")
    x=np.array(x)
    return (lower+((x-x_min)*(upper-lower))/(x_max-x_min))

def min_max_recover(x,x_min,x_max,lower=0.,upper=1.):
    '''
    reverse of min_max_scale function. The rescaled x is recovered to original scale.
    The data in [lower,upper] is recovered to [x_min,x_max] range.
    # it can hanldes np.nan.
    Args:
        x (1d np.array): data that needs to be rescaled.
        x_min (numeric): minimum value of x.
        x_max (numeric): maximum value of x.
        lower (numeric): minimum value of rescaled x. default is 0.
        upper (numeric): maximum value of rescaled x. default is 1.
    Raises:
    Returns:
        original scale x (numeric): rescaled x is recovered to [lower,upper] range from [x_min,x_max].
    '''
    if type(x)!=np.ndarray:
        raise ValueError("Put data as numpy ndarray.format.")
    x=np.array(x)
    return ((x-lower)*(x_max-x_min)/(upper-lower)+x_min)

def scale_by_name(df,scale_const,var_name):
    '''
    Apply min_max_scale for df['var_name']
    Args:
        df (pandas dataframe): data frame that contains various variables that needs to be scaled.
        scale_const (dictionary): dictionary that contains x_min, x_max, lower, and upper for var_name variable in df.
        var_name (str): a variable that needs to be scaled. One of the columns in df.
    Raises:
    Returns:
        scaled df['var_name] (numeric): scaled df['var_name] as a numpy array.
    '''
    if np.isin(var_name,df.columns).item():
        pass
    else:
        raise ValueError(f"{var_name} does not exist in the input dataframe.")

    x_max=scale_const[f'{var_name}_max']
    x_min=scale_const[f'{var_name}_min']
    lower=scale_const[f'{var_name}_lower']
    upper=scale_const[f'{var_name}_upper']
    x=df[var_name].to_numpy()
    x_scale=min_max_scale(x,x_min,x_max,lower,upper)
    #df[var_name]=df_scale.apply(lambda x :min_max_scale(x=x['sp_cool'],x_min=t_min,x_max=t_max,lower=-1.,upper=1.),axis=1)
    return x_scale

def get_scale_const_template():
    '''
    The tempelate for scale_const. Feel free to revise for your case.

    Args:
    Raises:
    Returns:
        scale_const (dictionary): dictionary that contains min, max, lower, upper information of T_out, T_in, wb_in, net, heat, cool,df,aux.
    '''
    scale_const={"T_out_max":40.,"T_out_min":-20.,"T_out_upper":1.,"T_out_lower":-1.,\
                "T_in_max":30.,"T_in_min":10.,"T_in_upper":1.,"T_in_lower":-1.,\
                "wb_in_max":30.,"wb_in_min":0.,"wb_in_upper":1.,"wb_in_lower":-1.,\
             "net_max":10000.0,"net_min":0.0,"heat_max":2000,"cool_max":2000,"df_max":4000,"aux_max":6000}
    # what is 0 in scaled T_out?
    #-1/3=-1+((0+20)*(1+1))/(40+20)
    return scale_const

def get_last_outputs(output_path):
    # get the latest input_values and prior_values file path given the output_path
    xx=output_path.glob("*")
    xxx=[x for x in xx if (x.is_file()) and (re.search('prior_values',x.__str__())) ]
    last_prior_values_path=xxx[-1]
    xx=output_path.glob("*")
    xxx2=[x for x in xx if (x.is_file()) and (re.search('input_values',x.__str__())) ]
    last_input_values_path=xxx2[-1]
    print(f'last_prior_values_path is {last_prior_values_path}')
    print(f'last_input_values_path is {last_input_values_path}')
    with open(last_input_values_path.__str__(),"rb") as handle:
        last_input_values=pickle.load(handle)
    with open(last_prior_values_path.__str__(),"rb") as handle:
        last_prior_values=pickle.load(handle)

    return last_input_values, last_prior_values,last_input_values_path,last_prior_values_path

############################################################################33
def detect_operation(operation,operation_state):
    # operation is vector of operation_state (numpy object vector with string data)
    # operation_state is hc system operation state \in {heat1,cool1,aux1,heat1_aux1,idle}
    # To detect heatpump heating, it should be ['heat1','heat1_aux1']
    if operation.dtype=="O":
        pass
    else:
        # change dtype to Object to have str and nan together
        operation=operation.astype("O")
    # check if any 'nan' and put np.nan for the Object type vector
    if np.any(operation==str(np.nan)):
        operation[operation==str(np.nan)]=np.nan

    if type(operation_state)==np.ndarray:
        if operation_state.dtype=="O":
            pass
        else:
            operation_state=operation_state.astype("O")
    else:
        operation_state=np.array([operation_state],dtype=object).flatten()

    i_vec=np.zeros_like(operation).astype("float32")
    # put nan for missing inputs
    i_vec[pd.isna(operation)]=np.nan

    i_vec[(np.isin(operation,operation_state)&(~pd.isna(operation)))]=1.0
    return i_vec

# ,T_out=None,df_cutpoint=0.0
def detect_defrost(i_vec,T_out,df_cutpoint=0):
    i_heat=i_vec.copy()
    i_df=i_vec.copy()
    i_heat[(T_out<=df_cutpoint) ]=0.0
    i_df[(T_out>df_cutpoint) ]=0.0
    return i_heat, i_df

def filter_na_outlier(x,x_max,x_min):
    x[x>(x_max)]=np.nan
    x[x<(x_min)]=np.nan
    return x

def nan_mean(vec):
    if np.all(np.isnan(vec)):
        return(np.nan)
    else:
        vec[np.isnan(vec)]=0.

        return(np.mean(vec))

def create_unit_input(thermostat_data,meter_data,scale_const,time_interval=15,training=False):

    # checking timestamp is pd.Timestamp object.
    if pd.core.dtypes.common.is_datetime_or_timedelta_dtype(thermostat_data['timestamp']):
        pass
    else:
        thermostat_data['timestamp']=pd.to_datetime(thermostat_data['timestamp'])
    if pd.core.dtypes.common.is_datetime_or_timedelta_dtype(meter_data['timestamp']):
        pass
    else:
        meter_data['timestamp']=pd.to_datetime(meter_data['timestamp'])

    start_date=thermostat_data['timestamp'][0].strftime("%Y-%m-%d")
    end_date=thermostat_data.tail(1)['timestamp'].reset_index(drop=True)[0].strftime("%Y-%m-%d")

    # create regular time inerval data
    # 5min time grid
    num_time_grid=((pd.Timestamp(end_date)-pd.Timestamp(start_date)+pd.Timedelta("1day")).total_seconds())/300
    time_grid=pd.date_range(pd.Timestamp(start_date), periods=num_time_grid, freq='5min')
    thermostat_data_base=pd.DataFrame(data={"timestamp":time_grid})
    thermostat_data=pd.merge(thermostat_data_base, thermostat_data, how='left', on=['timestamp'])

    thermostat_time_delta=((thermostat_data['timestamp'][1]-thermostat_data['timestamp'][0]).total_seconds())/60 # in minutes
    if (thermostat_time_delta!=5):
        raise ValueError("thermostat_data time interval should be 5 minute interval.")

    # time grid for meter_data
    meter_time_delta=((meter_data['timestamp'][1]-meter_data['timestamp'][0]).total_seconds())/60 # in minutes
    num_time_grid_meter=((pd.Timestamp(end_date)-pd.Timestamp(start_date)+pd.Timedelta("1day")).total_seconds())/(meter_time_delta*60) #
    time_grid_meter=pd.date_range(pd.Timestamp(start_date), periods=num_time_grid_meter, freq=f'{meter_time_delta}min')
    meter_data_base=pd.DataFrame(data={"timestamp":time_grid_meter})
    meter_data=pd.merge(meter_data_base, meter_data, how='left', on=['timestamp'])

    # make meter_data to 5 minute interval data
    if (meter_time_delta>5) and (meter_time_delta%5==0):
        meter_data=meter_data.iloc[np.repeat(np.arange(len(meter_data)), int(meter_time_delta/5))].reset_index(drop=True).copy()
        meter_data['timestamp']=time_grid
    elif (meter_time_delta==5):
        pass
    else:
        raise ValueError("meter_data time interval should be 5 minute interval or larger interval one with multiple of 5.")

    wb_in=get_wetbulb(temp=thermostat_data['T_in'].to_numpy(),rh=thermostat_data['rh_in'].to_numpy())
    wb_out=get_wetbulb(temp=thermostat_data['T_out'].to_numpy(),rh=thermostat_data['rh_out'].to_numpy())
    thermostat_data['wb_in']=wb_in
    thermostat_data['wb_out']=wb_out
    T_out=thermostat_data['T_out'].to_numpy()

    i_heat1_all=detect_operation(operation=thermostat_data['operation'],operation_state=['heat1','heat1_aux1'])
    i_heat2_all=detect_operation(operation=thermostat_data['operation'],operation_state=['heat2','heat2_aux1'])
    i_cool1=detect_operation(operation=thermostat_data['operation'],operation_state=['cool1'])
    i_cool2=detect_operation(operation=thermostat_data['operation'],operation_state=['cool2'])
    i_aux1=detect_operation(operation=thermostat_data['operation'],operation_state=['aux1','heat1_aux1','heat2_aux1'])
    i_fan1=detect_operation(operation=thermostat_data['operation'],operation_state=['fan1'])

    unit_data=pd.merge(thermostat_data, meter_data, how='left', on=['timestamp', 'unitcode'])

    net_max2=scale_const['net_max']*2.
    heat_max2=scale_const['heat_max']*2.
    cool_max2=scale_const['cool_max']*2.
    aux_max2=scale_const['aux_max']*2.
    df_max2=scale_const['df_max']*2.

    if type(unit_data['timestamp'][0])==str:
        unit_data['timestamp']=pd.to_datetime(unit_data['timestamp'].str.extract(r'(\d\d\d\d-\d\d-\d\d \d\d:\d\d:\d\d)')[0])

    unit_input=pd.DataFrame(data={"timestamp":unit_data['timestamp'].reset_index(drop=True),
                        "T_out":unit_data['T_out'].to_numpy(),
                        "T_in":unit_data['T_in'].to_numpy(),
                        "wb_in":unit_data['wb_in'].to_numpy(),
                        "net":unit_data['net'].to_numpy(),
                        "i_heat1_all":i_heat1_all,
                        "i_heat2_all":i_heat2_all,
                        "i_cool1":i_cool1,
                        "i_cool2":i_cool2,
                        "i_aux1":i_aux1,
                        "i_fan1":i_fan1,
                        "operation":unit_data['operation']
                        })

    unit_input['net']=filter_na_outlier(x=unit_input['net'].to_numpy(),x_max=net_max2,x_min=0.)
    if training:
        # training only no heatpump ahu data
        pass
    else:
        unit_input['ahu']=unit_data['ahu'].to_numpy()
        unit_input['heatpump']=unit_data['heatpump'].to_numpy()
        unit_input['misc']=unit_data['net'].to_numpy()-unit_data['ahu'].to_numpy()-unit_data['heatpump'].to_numpy()
        ahu=unit_input['ahu'].to_numpy()
        unit_input['ahu']=filter_na_outlier(x=ahu,x_max=aux_max2,x_min=0.)
        heatpump=unit_input['heatpump'].to_numpy()
        unit_input['heatpump']=filter_na_outlier(x=heatpump,x_max=heat_max2,x_min=0.)
        #unit_input['misc']=fill_na_outlier(x=unit_input['misc'].to_numpy(),x_max=misc_max2,x_min=0.)

    if time_interval==5:
        output=unit_input.copy()
        #pass
    else:
        #unit_input_15=unit_input.copy()
        unit_input['timestamp']=unit_input['timestamp'].dt.floor(f'{time_interval}min')
        output=unit_input.groupby('timestamp').agg(nan_mean).reset_index().copy()
        #output=unit_input.groupby('timestamp').mean().reset_index().copy()
    return output

def create_input_values(unit_input,scale_const=None,training=True):

    if training:
        print("Drop NA data to create training data.")
        nan_index=np.any(unit_input.isnull().to_numpy(),axis=1)
        timestamp=unit_input['timestamp']
        unit_input=unit_input.dropna()
    else:
        timestamp=unit_input['timestamp']
        nan_index=np.repeat(True,unit_input.shape[0])

    time_delta=(unit_input['timestamp'][1]-unit_input['timestamp'][0])/np.timedelta64(1,'m')
    net=unit_input['net'].to_numpy()

    T_out=unit_input['T_out'].to_numpy()
    T_in=unit_input['T_in'].to_numpy()
    wb_in=unit_input['wb_in'].to_numpy()
    i_heat1_all=unit_input['i_heat1_all'].to_numpy()
    i_heat2_all=unit_input['i_heat2_all'].to_numpy()
    i_cool1=unit_input['i_cool1'].to_numpy()
    i_cool2=unit_input['i_cool2'].to_numpy()

    i_aux1=unit_input['i_aux1'].to_numpy()
    i_fan1=unit_input['i_fan1'].to_numpy()
    i_hc=i_heat1_all+i_heat2_all+i_cool1+i_cool2+i_aux1+i_fan1 # hc signal
    #i_hc=i_heat1_all+i_cool1+i_aux1+i_fan1 # hc signal

    if scale_const is None:
        scale_const=get_scale_const_template()
        print(f"Template scale_const is used. The values are {scale_const}. Please specify scale_const as a function input by calling `get_scale_const_template()` function.")
    else:
        print(f'scale const is {scale_const}.')

    s_T_out=min_max_scale(T_out,x_min=scale_const['T_out_min'],x_max=scale_const['T_out_max'],
                  lower=scale_const['T_out_lower'],upper=scale_const['T_out_upper'])

    s_T_in=min_max_scale(T_in,x_min=scale_const['T_in_min'],x_max=scale_const['T_in_max'],
                  lower=scale_const['T_in_lower'],upper=scale_const['T_in_upper'])

    s_wb_in=min_max_scale(wb_in,x_min=scale_const['wb_in_min'],x_max=scale_const['wb_in_max'],
                  lower=scale_const['wb_in_lower'],upper=scale_const['wb_in_upper'])

    s_net=net/scale_const['net_max'] # 0-1 scale. Watt

    input_values={"net":net,
                  "s_net":s_net,
                  "T_out":T_out,
                  "s_T_out":s_T_out,
                  "T_in":T_in,
                  "s_T_in":s_T_in,
                  "wb_in":wb_in,
                  "s_wb_in":s_wb_in,
                  "i_heat1_all":i_heat1_all,
                  "i_heat2_all":i_heat2_all,
                  "i_cool1":i_cool1,
                  "i_cool2":i_cool2,
                  "i_aux1":i_aux1,
                  "i_fan1":i_fan1,
                  "i_hc":i_hc,
                  "nan_index":nan_index,
                  "scale_const":scale_const,
                  "timestamp":timestamp,
                  "time_delta":time_delta
                  }
    # in case we have validation data
    if ('ahu' in unit_input.columns) and ('heatpump' in unit_input.columns):
        hc=unit_input['heatpump'].to_numpy()+unit_input['ahu'].to_numpy()
        input_values.update({"hc":hc})

    return input_values

max_stars_repo_name: ecosang/rhced

max_issues_repo_name: ecosang/rhced

max_forks_repo_name: ecosang/rhced


text:
import tensorflow as tf
from tensorflow.keras import datasets,Sequential,layers,optimizers,metrics
import os
import numpy as np
from tensorflow import keras
tf.random.set_seed(11)
np.random.seed(111)

total_words=10000
max_len=80
batchsz=128
embedding_len=100
(x_train,y_train),(x_test,y_test)=datasets.imdb.load_data(num_words=total_words)
x_train=tf.keras.preprocessing.sequence.pad_sequences(x_train,maxlen=max_len)
x_test=tf.keras.preprocessing.sequence.pad_sequences(x_test,maxlen=max_len)
db_train=tf.data.Dataset.from_tensor_slices((x_train,y_train))
db_train=db_train.shuffle(1000).batch(batchsz,drop_remainder=True)
db_test=tf.data.Dataset.from_tensor_slices((x_test,y_test))
db_test=db_test.batch(batchsz,drop_remainder=True)
print('x_train:',x_train.shape,tf.reduce_min(y_train),tf.reduce_max(y_train))
print('x_test:',x_test.shape)

class Myrnn(keras.Model):
    def __init__(self,units):
        super(Myrnn,self).__init__()
        self.state0=[tf.zeros([batchsz,units]),tf.zeros([batchsz,units])]
        self.state1=[tf.zeros([batchsz,units]),tf.zeros([batchsz,units])]
        self.embedding=layers.Embedding(total_words,embedding_len,
                                        input_length=max_len)
        self.rnn_cell0=layers.LSTMCell(units,dropout=0.2)
        self.rnn_cell1=layers.LSTMCell(units,dropout=0.2)
        self.fc=layers.Dense(1)

    def call(self,inputs,training=None):
        x=inputs
        x=self.embedding(x)
        state0=self.state0
        state1=self.state1
        for word in tf.unstack(x,axis=1):
            out0,state0=self.rnn_cell0(word,state0,training)
            out1,state1=self.rnn_cell1(out0,state1,training)

        x=self.fc(out1)
        prob=tf.sigmoid(x)
        return prob

def main():
    units=64
    epochs=4
    model=Myrnn(units)
    # model.build(input_shape=[None,units])
    # model.summary()
    model.compile(optimizer=optimizers.Adam(lr=1e-3),
                  loss=tf.losses.BinaryCrossentropy(),
                  metrics=['accuracy'])
    model.fit(db_train,epochs=epochs,validation_data=db_test,validation_freq=1)
    model.summary()

    model.evaluate(db_test)

if __name__ == '__main__':
    main()

ext: py

reason: import numpy


text:
from datetime import datetime
from math import sin, cos, sqrt, fabs, atan2
from gps_time import GPSTime, utc_to_gpst
import numpy as np
from constants import SPEED_OF_LIGHT, SECS_IN_MIN, SECS_IN_HR, SECS_IN_DAY, EARTH_ROTATION_RATE, EARTH_GM
from helpers import get_constellation

def read4(f, rinex_ver):
  line = f.readline()[:-1]
  if rinex_ver == 2:
    line = ' ' + line  # Shift 1 char to the right
  line = line.replace('D', 'E')  # Handle bizarro float format
  return float(line[4:23]), float(line[23:42]), float(line[42:61]), float(line[61:80])

def convert_ublox_ephem(ublox_ephem):
  ephem = {}
  if ublox_ephem.gpsWeek < 1024:
    week = ublox_ephem.gpsWeek + 1024
  else:
    week = ublox_ephem.gpsWeek
  ephem['toe'] = GPSTime(week, ublox_ephem.toe)
  ephem['toc'] = GPSTime(week, ublox_ephem.toc)
  ephem['af0'] = ublox_ephem.af0
  ephem['af1'] = ublox_ephem.af1
  ephem['af2'] = ublox_ephem.af2
  ephem['tgd'] = ublox_ephem.tgd

  ephem['sqrta'] = np.sqrt(ublox_ephem.a)
  ephem['dn'] = ublox_ephem.deltaN
  ephem['m0'] = ublox_ephem.m0

  ephem['ecc'] = ublox_ephem.ecc
  ephem['w'] = ublox_ephem.omega
  ephem['cus'] = ublox_ephem.cus
  ephem['cuc'] = ublox_ephem.cuc
  ephem['crc'] = ublox_ephem.crc
  ephem['crs'] = ublox_ephem.crs
  ephem['cic'] = ublox_ephem.cic
  ephem['cis'] = ublox_ephem.cis

  ephem['inc'] = ublox_ephem.i0
  ephem['inc_dot'] = ublox_ephem.iDot
  ephem['omegadot'] = ublox_ephem.omegaDot
  ephem['omega0'] = ublox_ephem.omega0

  return ephem

class EphemerisType(object):
  # TODO this isn't properly supported
  NAV = 0
  FINAL_ORBIT = 1
  RAPID_ORBIT = 2
  ULTRA_RAPID_ORBIT = 3
  QCOM_POLY = 4

class Ephemeris(object):
  def valid(self, time):
    return abs(time - self.epoch) <= self.max_time_diff

class GLONASSEphemeris(Ephemeris):
  def __init__(self, data, epoch, healthy=True):
    self.prn = data['prn']
    self.epoch = epoch
    self.healthy = healthy
    self.data = data
    self.max_time_diff = 25*SECS_IN_MIN
    self.type = EphemerisType.NAV
    self.channel = data['freq_num']

  def get_sat_info(self, time):
    if not self.healthy:
      return None
    # see the russian doc for this:
    # http://gauss.gge.unb.ca/GLONASS.ICD.pdf

    eph = self.data
    # TODO should handle leap seconds better
    toc_gps_time = utc_to_gpst(eph['toc'])
    tdiff = time - toc_gps_time

    # Clock correction (except for general relativity which is applied later)
    clock_err = eph['min_tauN'] + tdiff * (eph['GammaN'])
    clock_rate_err = eph['GammaN']

    def glonass_diff_eq(state, acc):
      J2 = 1.0826257e-3
      mu = 3.9860044e14
      omega = 7.292115e-5
      ae = 6378136.0
      r = np.sqrt(state[0]**2 + state[1]**2 + state[2]**2)
      ders = np.zeros(6)
      if r**2 < 0:
        return ders
      a = 1.5 * J2 * mu * (ae**2)/ (r**5)
      b = 5 * (state[2]**2) / (r**2)
      c = -mu/(r**3) - a*(1-b)

      ders[0:3] = state[3:6]
      ders[3] = (c + omega**2)*state[0] + 2*omega*state[4] + acc[0]
      ders[4] = (c + omega**2)*state[1] - 2*omega*state[3] + acc[1]
      ders[5] = (c - 2*a)*state[2] + acc[2]
      return ders

    init_state = np.empty(6)
    init_state[0] = eph['x']
    init_state[1] = eph['y']
    init_state[2] = eph['z']
    init_state[3] = eph['x_vel']
    init_state[4] = eph['y_vel']
    init_state[5] = eph['z_vel']
    init_state = 1000*init_state
    acc = 1000*np.array([eph['x_acc'], eph['y_acc'], eph['z_acc']])
    state = init_state
    tstep = 90
    if tdiff < 0:
      tt = -tstep
    elif tdiff > 0:
      tt = tstep
    while abs(tdiff) > 1e-9:
      if abs(tdiff) < tstep:
        tt = tdiff
      k1 = glonass_diff_eq(state, acc)
      k2 = glonass_diff_eq(state + k1*tt/2, -acc)
      k3 = glonass_diff_eq(state + k2*tt/2, -acc)
      k4 = glonass_diff_eq(state + k3*tt, -acc)
      state += (k1 + 2*k2 + 2*k3 + k4)*tt/6.0
      tdiff -= tt

    pos = state[0:3]
    vel = state[3:6]
    return pos, vel, clock_err, clock_rate_err

class PolyEphemeris(Ephemeris):
  def __init__(self, prn, data, epoch, healthy=True, eph_type=None, tgd=0):
    self.prn = prn
    self.epoch = epoch
    self.healthy = healthy
    self.data = data
    self.tgd = tgd
    self.max_time_diff = SECS_IN_HR
    self.type = eph_type

  def get_sat_info(self, time):
    if not self.healthy:
      return None
    dt = time - self.data['t0']
    deg = self.data['deg']
    deg_t = self.data['deg_t']
    sat_pos = np.array([sum([(dt**p)*self.data['x'][deg-p] for p in xrange(deg+1)]),
                        sum([(dt**p)*self.data['y'][deg-p] for p in xrange(deg+1)]),
                        sum([(dt**p)*self.data['z'][deg-p] for p in xrange(deg+1)])])
    sat_vel = np.array([sum([p*(dt**(p-1))*self.data['x'][deg-p] for p in xrange(1,deg+1)]),
                        sum([p*(dt**(p-1))*self.data['y'][deg-p] for p in xrange(1,deg+1)]),
                        sum([p*(dt**(p-1))*self.data['z'][deg-p] for p in xrange(1,deg+1)])])
    time_err = sum([(dt**p)*self.data['clock'][deg_t-p] for p in xrange(deg_t+1)])
    time_err_rate = sum([p*(dt**(p-1))*self.data['clock'][deg_t-p] for p in xrange(1,deg_t+1)])
    time_err_with_rel = time_err - 2*np.inner(sat_pos, sat_vel)/SPEED_OF_LIGHT**2
    return sat_pos, sat_vel, time_err_with_rel, time_err_rate

class GPSEphemeris(Ephemeris):
  def __init__(self, data, epoch, healthy=True):
    self.prn = 'G%02i' % data['prn']
    self.epoch = epoch
    self.healthy = healthy
    self.data = data
    self.max_time_diff = 2*SECS_IN_HR
    self.max_time_diff_tgd = SECS_IN_DAY
    self.type = EphemerisType.NAV

  def get_tgd(self):
    return self.data['tgd']

  def get_sat_info(self, time):
    if not self.healthy:
      return None
    eph = self.data
    tdiff = time - eph['toe']  # Time of clock
    clock_err = eph['af0'] + tdiff * (eph['af1'] + tdiff * eph['af2'])
    clock_rate_err = eph['af1'] + 2 * tdiff * eph['af2']

    # Orbit propagation
    tdiff = time - eph['toe']  # Time of ephemeris (might be different from time of clock)

    # Calculate position per IS-GPS-200D p 97 Table 20-IV
    a = eph['sqrta'] * eph['sqrta']  # [m] Semi-major axis
    ma_dot = sqrt(EARTH_GM / (a * a * a)) + eph['dn']  # [rad/sec] Corrected mean motion
    ma = eph['m0'] + ma_dot * tdiff  # [rad] Corrected mean anomaly

    # Iteratively solve for the Eccentric Anomaly (from Keith Alter and David Johnston)
    ea = ma  # Starting value for E

    ea_old = 2222
    while fabs(ea - ea_old) > 1.0E-14:
      ea_old = ea
      tempd1 = 1.0 - eph['ecc'] * cos(ea_old)
      ea = ea + (ma - ea_old + eph['ecc'] * sin(ea_old)) / tempd1
    ea_dot = ma_dot / tempd1

    # Relativistic correction term
    einstein = -4.442807633E-10 * eph['ecc'] * eph['sqrta'] * sin(ea)

    # Begin calc for True Anomaly and Argument of Latitude
    tempd2 = sqrt(1.0 - eph['ecc'] * eph['ecc'])
    # [rad] Argument of Latitude = True Anomaly + Argument of Perigee
    al = atan2(tempd2 * sin(ea), cos(ea) - eph['ecc']) + eph['w']
    al_dot = tempd2 * ea_dot / tempd1

    # Calculate corrected argument of latitude based on position
    cal = al + eph['cus'] * sin(2.0 * al) + eph['cuc'] * cos(2.0 * al)
    cal_dot = al_dot * (1.0 + 2.0 * (eph['cus'] * cos(2.0 * al) -
                                     eph['cuc'] * sin(2.0 * al)))

    # Calculate corrected radius based on argument of latitude
    r = a * tempd1 + eph['crc'] * cos(2.0 * al) + eph['crs'] * sin(2.0 * al)
    r_dot = (a * eph['ecc'] * sin(ea) * ea_dot +
             2.0 * al_dot * (eph['crs'] * cos(2.0 * al) -
                             eph['crc'] * sin(2.0 * al)))

    # Calculate inclination based on argument of latitude
    inc = (eph['inc'] + eph['inc_dot'] * tdiff +
           eph['cic'] * cos(2.0 * al) +
           eph['cis'] * sin(2.0 * al))
    inc_dot = (eph['inc_dot'] +
               2.0 * al_dot * (eph['cis'] * cos(2.0 * al) -
                               eph['cic'] * sin(2.0 * al)))

    # Calculate position and velocity in orbital plane
    x = r * cos(cal)
    y = r * sin(cal)
    x_dot = r_dot * cos(cal) - y * cal_dot
    y_dot = r_dot * sin(cal) + x * cal_dot

    # Corrected longitude of ascending node
    om_dot = eph['omegadot'] - EARTH_ROTATION_RATE
    om = eph['omega0'] + tdiff * om_dot - EARTH_ROTATION_RATE * eph['toe'].tow

    # Compute the satellite's position in Earth-Centered Earth-Fixed coordiates
    pos = np.empty(3)
    pos[0] = x * cos(om) - y * cos(inc) * sin(om)
    pos[1] = x * sin(om) + y * cos(inc) * cos(om)
    pos[2] = y * sin(inc)

    tempd3 = y_dot * cos(inc) - y * sin(inc) * inc_dot

    # Compute the satellite's velocity in Earth-Centered Earth-Fixed coordiates
    vel = np.empty(3)
    vel[0] = -om_dot * pos[1] + x_dot * cos(om) - tempd3 * sin(om)
    vel[1] = om_dot * pos[0] + x_dot * sin(om) + tempd3 * cos(om)
    vel[2] = y * cos(inc) * inc_dot + y_dot * sin(inc)

    clock_err += einstein

    return pos, vel, clock_err, clock_rate_err

def parse_sp3_orbits(file_names, SUPPORTED_CONSTELLATIONS):
  ephems = []
  data = {}
  for file_name in file_names:
    f = open(file_name, 'r')
    while True:
      line = f.readline()[:-1]
      if not line:
        break
      # epoch header
      if line[0:2] == '* ':
        year = int(line[3:7])
        month = int(line[8:10])
        day = int(line[11:13])
        hour = int(line[14:16])
        minute = int(line[17:19])
        second = int(float(line[20:31]))
        epoch = GPSTime.from_datetime(datetime(year, month, day, hour, minute, second))
      # pos line
      elif line[0] == 'P':
        prn = line[1:4].replace(' ','0')
        if get_constellation(prn) not in SUPPORTED_CONSTELLATIONS:
          continue
        if prn not in data:
          data[prn] = []
        #TODO this is a crappy way to deal with overlapping ultra rapid
        if len(data[prn]) < 1 or epoch - data[prn][-1][0] > 0:
          parsed = [epoch,
                    1e3*float(line[4:18]),
                    1e3*float(line[18:32]),
                    1e3*float(line[32:46]),
                    1e-6*float(line[46:60])]
          if (np.array(parsed[1:]) != 0).all():
            data[prn].append(parsed)
    f.close()
  deg = 16
  deg_t = 1
  for prn in data:
    # TODO Handle this properly
    # Currently don't even bother with satellites that have unhealthy times
    if (np.array(data[prn])[:,4] > .99).any():
      continue
    for i in xrange(len(data[prn]) - deg):
      times, x, y, z, clock = [],[],[],[],[]
      epoch = data[prn][i + deg//2][0]
      for j in xrange(deg + 1):
        times.append(data[prn][i + j][0] - epoch)
        x.append(data[prn][i + j][1])
        y.append(data[prn][i + j][2])
        z.append(data[prn][i + j][3])
        clock.append(data[prn][i + j][4])
      if (np.diff(times) != 900).any():
        continue
      poly_data = {}
      poly_data['t0'] = epoch
      poly_data['x'] = np.polyfit(times, x, deg)
      poly_data['y'] = np.polyfit(times, y, deg)
      poly_data['z'] = np.polyfit(times, z, deg)
      poly_data['clock'] = [(data[prn][i + deg//2 + 1][4] - data[prn][i + deg//2 - 1][4])/1800, data[prn][i + deg//2][4]]
      poly_data['deg'] = deg
      poly_data['deg_t'] = deg_t
      ephems.append(PolyEphemeris(prn, poly_data, epoch, healthy=True, eph_type=EphemerisType.RAPID_ORBIT))
  return ephems

def parse_rinex_nav_msg_gps(file_name):
  ephems = []
  got_header = False
  rinex_ver = None
  #ion_alpha = None
  #ion_beta = None
  f = open(file_name, 'r')
  while True:
    line = f.readline()[:-1]
    if not line:
      break
    if not got_header:
      if rinex_ver is None:
        if line[60:80] != "RINEX VERSION / TYPE":
          raise RuntimeError("Doesn't appear to be a RINEX file")
        rinex_ver = int(float(line[0:9]))
        if line[20] != "N":
          raise RuntimeError("Doesn't appear to be a Navigation Message file")
      #if line[60:69] == "ION ALPHA":
      #  line = line.replace('D', 'E')  # Handle bizarro float format
      #  ion_alpha= [float(line[3:14]), float(line[15:26]), float(line[27:38]), float(line[39:50])]
      #if line[60:68] == "ION BETA":
      #  line = line.replace('D', 'E')  # Handle bizarro float format
      #  ion_beta= [float(line[3:14]), float(line[15:26]), float(line[27:38]), float(line[39:50])]
      if line[60:73] == "END OF HEADER":
        #ion = ion_alpha + ion_beta
        got_header = True
      continue
    if rinex_ver == 3:
      if line[0] != 'G':
        continue
    if rinex_ver == 3:
      prn = int(line[1:3])
      epoch = GPSTime.from_datetime(datetime.strptime(line[4:23], "%y %m %d %H %M %S"))
    elif rinex_ver == 2:
      prn = int(line[0:2])
      epoch = GPSTime.from_datetime(datetime.strptime(line[3:20], "%y %m %d %H %M %S"))
      line = ' ' + line  # Shift 1 char to the right

    line = line.replace('D', 'E')  # Handle bizarro float format
    e = {'epoch': epoch, 'prn': prn}
    e['toc'] = epoch
    e['af0'] = float(line[23:42])
    e['af1'] = float(line[42:61])
    e['af2'] = float(line[61:80])

    e['iode'], e['crs'], e['dn'], e['m0'] = read4(f, rinex_ver)
    e['cuc'], e['ecc'], e['cus'], e['sqrta'] = read4(f, rinex_ver)
    toe_tow, e['cic'], e['omega0'], e['cis'] = read4(f, rinex_ver)
    e['inc'], e['crc'], e['w'], e['omegadot'] = read4(f, rinex_ver)
    e['inc_dot'], e['l2_codes'], toe_week, e['l2_pflag'] = read4(f, rinex_ver)
    e['sv_accuracy'], e['health'], e['tgd'], e['iodc'] = read4(f, rinex_ver)
    f.readline()  # Discard last row

    e['toe'] = GPSTime(toe_week, toe_tow)
    e['healthy'] = (e['health'] == 0.0)

    ephems.append(GPSEphemeris(e, epoch))
  f.close()
  return ephems

def parse_rinex_nav_msg_glonass(file_name):
  ephems = []
  f = open(file_name, 'r')
  got_header = False
  rinex_ver = None
  while True:
    line = f.readline()[:-1]
    if not line:
      break
    if not got_header:
      if rinex_ver is None:
        if line[60:80] != "RINEX VERSION / TYPE":
          raise RuntimeError("Doesn't appear to be a RINEX file")
        rinex_ver = int(float(line[0:9]))
        if line[20] != "G":
          raise RuntimeError("Doesn't appear to be a Navigation Message file")
      if line[60:73] == "END OF HEADER":
        got_header = True
      continue
    if rinex_ver == 3:
      prn = line[:3]
      epoch = GPSTime.from_datetime(datetime.strptime(line[4:23], "%y %m %d %H %M %S"))
    elif rinex_ver == 2:
      prn = 'R%02i' % int(line[0:2])
      epoch = GPSTime.from_datetime(datetime.strptime(line[3:20], "%y %m %d %H %M %S"))
      line = ' ' + line  # Shift 1 char to the right

    line = line.replace('D', 'E')  # Handle bizarro float format
    e = {'epoch': epoch, 'prn': prn}
    e['toc'] = epoch
    e['min_tauN'] = float(line[23:42])
    e['GammaN'] = float(line[42:61])
    e['tk'] = float(line[61:80])

    e['x'], e['x_vel'], e['x_acc'], e['health'] = read4(f, rinex_ver)
    e['y'], e['y_vel'], e['y_acc'], e['freq_num'] = read4(f, rinex_ver)
    e['z'], e['z_vel'], e['z_acc'], e['age'] = read4(f, rinex_ver)

    e['healthy'] = (e['health'] == 0.0)

    ephems.append(GLONASSEphemeris(e, epoch))
  f.close()
  return ephems

'''
def parse_ublox_ephems(ublox_ephems):
  ephems = []
  for ublox_ephem in ublox_ephems:
    svId = ublox_ephem.ubloxGnss.ephemeris.svId
    data = convert_ublox_ephem(ublox_ephem.ubloxGnss.ephemeris)
    epoch = data['toe']
    ephems.append(GPSEphemeris(svId, data, epoch))
  return ephems

def parse_qcom_ephems(qcom_polys, current_week):
  ephems = []
  for qcom_poly in qcom_polys:
    svId = qcom_poly.qcomGnss.drSvPoly.svId
    data = qcom_poly.qcomGnss.drSvPoly
    t0 = data.t0
    # fix glonass time
    if get_constellation(svId) == 'GLONASS':
      # TODO should handle leap seconds better
      epoch = GPSTime(current_week, (t0 + 3*SECS_IN_WEEK) % (SECS_IN_WEEK) + 18)
    else:
      epoch = GPSTime(current_week, t0)
    poly_data = {}
    poly_data['t0'] = epoch
    poly_data['x'] = [data.xyzN[2], data.xyzN[1], data.xyzN[0], data.xyz0[0]]
    poly_data['y'] = [data.xyzN[5], data.xyzN[4], data.xyzN[3], data.xyz0[1]]
    poly_data['z'] = [data.xyzN[8], data.xyzN[7], data.xyzN[6], data.xyz0[2]]
    poly_data['clock'] = [1e-3*data.other[3], 1e-3*data.other[2], 1e-3*data.other[1], 1e-3*data.other[0]]
    poly_data['deg'] = 3
    poly_data['deg_t'] = 3
    ephems.append(PolyEphemeris(svId, poly_data, epoch, eph_type=EphemerisType.QCOM_POLY))
  return ephems
'''


text:
import multiprocessing as mp
import collections
import itertools
import numpy as np
import matplotlib.pyplot as plt
from functools import partial
import qecsim
from qecsim import app
from qecsim.models.generic import PhaseFlipErrorModel
from qecsim.models.planar import PlanarCode, PlanarMPSDecoder

# set models
codes = [PlanarCode(*size) for size in [(3, 3), (5, 5)]]
error_model = PhaseFlipErrorModel()
decoder = PlanarMPSDecoder()
# set physical error probabilities
error_probability_min, error_probability_max = 0, 0.4
error_probabilities = np.linspace(error_probability_min, error_probability_max, 5)
# set max_runs for each probability
max_runs = 10

# print run parameters
print('Codes:', [code.label for code in codes])
print('Error model:', error_model.label)
print('Decoder:', decoder.label)
print('Error probabilities:', error_probabilities)
print('Maximum runs:', max_runs)

def parallel_step_p(code, error_model, decoder, max_runs,error_probability):
    return app.run(code, error_model, decoder, error_probability, max_runs=100)

pL_list =np.zeros((len(codes),len(error_probabilities)))
std_list=np.zeros((len(codes),len(error_probabilities)))

for code_index,code in enumerate(codes):
    p=mp.Pool()
    func=partial(parallel_step_p,code, error_model, decoder, max_runs)
    pL_list[code_index]=p.map(func, error_probabilities)
    p.close()
    p.join()

max_stars_repo_name: dua-arpit/qecsim

max_issues_repo_name: dua-arpit/qecsim

max_forks_repo_name: dua-arpit/qecsim


text:
from abc import ABC, abstractmethod

import numpy as np
import pandas as pd

from pvrpm.core.enums import ConfigKeys as ck
from pvrpm.core.case import SamCase
from pvrpm.core.utils import sample, get_higher_components

class Monitor(ABC):
    """
    Abstract class for monitoring types
    """

    def __init__(
        self,
        level: str,
        comp_level_df: pd.DataFrame,
        case: SamCase,
    ):
        """
        Initalizes a monitoring instance

        Args:
            level (str): The component level this monitoring is apart of
            comp_level_df (:obj:`pd.DataFrame`): The component level dataframe containing the simulation data
            case (:obj:`SamCase`): The SAM case for this simulation
        """
        super().__init__()
        self.level = level
        self.df = comp_level_df
        self.case = case
        self.initialize_components()

    @abstractmethod
    def initialize_components(self):
        """
        Initalizes monitoring data for all components to be tracked during simulation for this monitor type

        Note:
            Updates the underlying dataframes in place
        """
        pass

    @abstractmethod
    def reinitialize_components(self, df: pd.DataFrame) -> pd.DataFrame:
        """
        Reinitalize components in a dataframe similiar to the inital initalization. Used for when repairs or other things may occur

        Args:
            df (:obj:`pd.DataFrame`): The dataframe containing the components to reinitalize

        Returns:
            :obj:`pd.DataFrame`: The reinitalized components
        """
        pass

    @abstractmethod
    def update(self, day: int):
        """
        Perform a monitoring update for one day in the simulation:

        Modifies time to detection for components

        Args:
            day (int): Current day in the simulation

        Note:
            Updates the underlying dataframes in place
        """
        pass

class LevelMonitor(Monitor):
    """
    Defines monitoring at a component level
    """

    def initialize_components(self):
        component_info = self.case.config[self.level]
        df = self.df
        monitor_modes = list(component_info.get(ck.MONITORING, {}).keys())
        failure_modes = list(component_info.get(ck.FAILURE, {}).keys())

        # monitoring times, these will be added to the repair time for each component
        # basically, the time until each failure is detected

        if len(monitor_modes) == 1:
            # same monitor mode for every failure
            monitor = component_info[ck.MONITORING][monitor_modes[0]]
            monitor = sample(monitor[ck.DIST], monitor[ck.PARAM], component_info[ck.NUM_COMPONENT])
            df["monitor_times"] = monitor
            df["time_to_detection"] = df["monitor_times"].copy()
        elif len(monitor_modes) > 1:
            failure_ind = [failure_modes.index(m) for m in df["failure_type"]]
            modes = [monitor_modes[i] for i in failure_ind]
            monitor = np.array(
                [
                    sample(component_info[ck.MONITORING][m][ck.DIST], component_info[ck.MONITORING][m][ck.PARAM], 1)[0]
                    for m in modes
                ]
            )
            df["monitor_times"] = monitor
            df["time_to_detection"] = df["monitor_times"].copy()

    def reinitialize_components(self, df: pd.DataFrame) -> pd.DataFrame:
        component_info = self.case.config[self.level]
        monitor_modes = list(component_info.get(ck.MONITORING, {}).keys())
        failure_modes = list(component_info.get(ck.FAILURE, {}).keys())

        if len(monitor_modes) == 1:
            # same monitor mode for every failure
            monitor = component_info[ck.MONITORING][monitor_modes[0]]
            df["monitor_times"] = sample(monitor[ck.DIST], monitor[ck.PARAM], len(df))
            df["time_to_detection"] = df["monitor_times"].copy()
        elif len(monitor_modes) > 1:
            failure_ind = [failure_modes.index(m) for m in df["failure_type"]]
            modes = [monitor_modes[i] for i in failure_ind]
            monitor = np.array(
                [
                    sample(component_info[ck.MONITORING][m][ck.DIST], component_info[ck.MONITORING][m][ck.PARAM], 1)[0]
                    for m in modes
                ]
            )
            df["monitor_times"] = monitor
            df["time_to_detection"] = df["monitor_times"].copy()

        return df

    def update(self, day: int):
        df = self.df

        mask = (df["state"] == 0) & (df["time_to_detection"] > 1) & (df["indep_monitor"] == False)
        df.loc[mask, "time_to_detection"] -= 1

class CrossLevelMonitor(Monitor):
    """
    Defines monitoring across different levels
    """

    def __init__(
        self,
        level: str,
        comp_level_df: pd.DataFrame,
        case: SamCase,
    ):
        """
        Initalizes a cross level monitoring instance

        Args:
            level (str): The component level this monitoring is apart of
            comp_level_df (:obj:`pd.DataFrame`): The component level dataframe containing the simulation data
            case (:obj:`SamCase`): The SAM case for this simulation
        """
        super().__init__(level, comp_level_df, case)
        if self.case.config[level][ck.COMP_MONITOR].get(ck.FAIL_PER_THRESH, None) is not None:
            top_level = self.case.config[level][ck.COMP_MONITOR][ck.LEVELS]
            bins = np.zeros(len(self.df))
            num_per_top = get_higher_components(top_level, level, self.case)
            indicies = np.floor(self.df.index / num_per_top)
            self.monitor_bins = {
                "top_level": top_level,
                "indicies": indicies,
                "bins": bins,
                "num_per_top": num_per_top,
            }

    def update_bins(self, df: pd.DataFrame):
        indicies, counts, _ = get_higher_components(
            self.monitor_bins["top_level"],
            self.level,
            self.case,
            start_level_df=df,
        )
        self.monitor_bins["bins"][indicies] = counts

    def initialize_components(self):
        component_info = self.case.config[self.level]
        df = self.df

        monitor = component_info[ck.COMP_MONITOR]
        df["monitor_times"] = sample(monitor[ck.DIST], monitor[ck.PARAM], component_info[ck.NUM_COMPONENT])
        df["time_to_detection"] = df["monitor_times"].copy()

    def reinitialize_components(self, df: pd.DataFrame) -> pd.DataFrame:
        component_info = self.case.config[self.level]
        monitor = component_info[ck.COMP_MONITOR]

        df["monitor_times"] = sample(monitor[ck.DIST], monitor[ck.PARAM], len(df))
        df["time_to_detection"] = df["monitor_times"].copy()

        return df

    def update(self, day: int):
        df = self.df

        # fraction failed is number of failed components, no matter if its being repaired, etc
        mask = df["state"] == 0
        failed_df = df.loc[mask]
        num_failed = len(failed_df)
        mask = mask & (df["time_to_detection"] > 1) & (df["indep_monitor"] == False)

        conf = self.case.config[self.level][ck.COMP_MONITOR]
        global_threshold = False
        if conf.get(ck.FAIL_THRESH, None) is not None:
            frac_failed = num_failed / len(df)
            # only decrement time to detection once failure threshold is met
            if frac_failed > conf[ck.FAIL_THRESH]:
                # TODO: compound the failures... later
                df.loc[mask, "time_to_detection"] -= 1
                global_threshold = True
            elif conf.get(ck.FAIL_PER_THRESH, None) is None:
                df.loc[mask, "monitor_times"] += 1

        if conf.get(ck.FAIL_PER_THRESH, None) is not None and not global_threshold:
            # TODO: idk how to make this more efficent, if i can reduce this down to just a single iloc or loc, and try to elimate the call to get_higher_components i can bump up run time, cause right now it doubles runtime per realization
            # first, need to get the proper components of this level under the top monitor level
            self.update_bins(failed_df)

            data = self.monitor_bins
            undetected_failed = df.loc[mask].copy()
            indicies, counts, total = get_higher_components(
                data["top_level"],
                self.level,
                self.case,
                start_level_df=undetected_failed,
            )
            # calculate failure threshold for each top level monitoring component for the monitoried components
            subtract = []
            add = []
            prev_cnt = 0
            # here, cnt is the number of failed components per top level component, no matter if its being repaired, etc
            for ind, cnt in zip(indicies, counts):
                # total currently failed regardless of detected, being repaired, etc
                frac_failed = data["bins"][ind] / total
                if frac_failed > conf[ck.FAIL_PER_THRESH]:
                    subtract += np.arange(prev_cnt, prev_cnt + cnt).tolist()
                else:
                    add += np.arange(prev_cnt, prev_cnt + cnt).tolist()
                prev_cnt += cnt

            if subtract:
                undetected_failed.iloc[subtract, undetected_failed.columns.get_loc("time_to_detection")] -= 1
            if add:
                undetected_failed.iloc[add, undetected_failed.columns.get_loc("monitor_times")] += 1
            df.loc[mask] = undetected_failed

class IndepMonitor(Monitor):
    """
    Defines independent monitoring that is component level independent
    """

    def __init__(
        self,
        case: SamCase,
        comps: dict,
        costs: dict,
        dc_availability: np.array,
    ):
        """
        Initalizes an independent monitoring instance

        Args:
            case (:obj:`SamCase`): The SAM case for this simulation
            comps (dict): The dataframes for each component level that tracks information during the simulation
            costs (dict): The dictionary for costs accrued for each component level
            dc_availability (:obj:`np.array`): The dc availability for every day in the simulation
        """
        super().__init__(None, None, case)
        self.comps = comps
        self.costs = costs
        self.labor_rate = self.case.config[ck.LABOR_RATE]
        self.dc_availability = dc_availability
        # keep track of static monitoring for every level
        index = []
        data = []
        for name, monitor_config in self.case.config.get(ck.INDEP_MONITOR, {}).items():
            index.append(name)
            if ck.INTERVAL in monitor_config:
                data.append(monitor_config[ck.INTERVAL])
            else:
                data.append(np.nan)

        if index:
            self.indep_monitoring = pd.Series(index=index, data=data)
        else:
            raise AttributeError("No static monitoring defined for any component level!")

    def initialize_components(self):
        pass

    def update_labor_rate(self, new: float):
        """
        Updates the labor rate for this repair

        Args:
            new (float): The new labor rate
        """
        self.labor_rate = new

    def reinitialize_components(self, df: pd.DataFrame) -> pd.DataFrame:
        df["monitor_times"] = None
        df["time_to_detection"] = None
        df["indep_monitor"] = False

        return df

    def update(self, day: int):
        # need to update time to detections for components with a distribution and need to decrement time to detection
        # this overrides other monitoring methods, if they are defined, using the indep_monitor column
        for name, monitor_config in self.case.config.get(ck.INDEP_MONITOR, {}).items():
            if ck.DIST in monitor_config:  # if no dist, time to detection is None or 0
                for level in monitor_config[ck.LEVELS]:
                    if not self.case.config[level][ck.CAN_FAIL]:
                        continue

                    df = self.comps[level]
                    mask = (df["state"] == 0) & (df["time_to_detection"] >= 1) & df["indep_monitor"]
                    df.loc[mask, "time_to_detection"] -= 1

        self.indep_monitoring -= 1
        mask = self.indep_monitoring < 1
        current_monitors = self.indep_monitoring.loc[mask].index
        for monitor in current_monitors:
            # independent monitoring method commences, instantly detecting all failed components in their specified levels
            # which have not been detected by other monitoring methods yet.
            config = self.case.config[ck.INDEP_MONITOR][monitor]
            for level in config[ck.LEVELS]:
                if not self.case.config[level][ck.CAN_FAIL]:
                    continue
                df = self.comps[level]
                mask = (
                    (df["state"] == 0)
                    & ((df["time_to_detection"] >= 1) | (df["time_to_detection"].isna()))
                    & (df["indep_monitor"] == False)
                )
                monitor_comps = df.loc[mask].copy()
                if monitor_comps["time_to_detection"].isna().sum() == 0:
                    monitor_comps["monitor_times"] -= monitor_comps["time_to_detection"]

                if ck.DIST in config:
                    monitor_comps["time_to_detection"] = sample(config[ck.DIST], config[ck.PARAM], len(monitor_comps))
                else:
                    monitor_comps["time_to_detection"] = 0

                monitor_comps["indep_monitor"] = True

                df.loc[mask] = monitor_comps
                # add cost to total costs (split evenly among all the levels for this day)
                self.costs[level][day] += (config[ck.COST] + self.labor_rate * config[ck.LABOR]) / len(
                    config[ck.LEVELS]
                )

            # reset static monitoring time
            self.indep_monitoring[monitor] = self.case.config[ck.INDEP_MONITOR][monitor][ck.INTERVAL]

        # threshold monitoring calculations
        for name, monitor_config in self.case.config.get(ck.INDEP_MONITOR, {}).items():
            # ignore threshold if the interval was reached for this monitor
            if name in current_monitors or (
                ck.FAIL_THRESH not in monitor_config and ck.FAIL_PER_THRESH not in monitor_config
            ):
                continue

            threshold_met = False
            if ck.FAIL_THRESH in monitor_config:
                # calculate dc availability and determine if indep monitoring is needed
                # this does use the previous day since this is calculated at the end of each day, but i figured
                # that wouldn't matter in the end and saves computing it again here
                # subtracting 1 also works since the inital dc availability is calculated for day 0
                # and day starts at 1
                if (1 - (self.dc_availability[day - 1])) > monitor_config[ck.FAIL_THRESH]:
                    threshold_met = True

            if not threshold_met and ck.FAIL_PER_THRESH in monitor_config:
                for level in monitor_config[ck.LEVELS]:
                    if not self.case.config[level][ck.CAN_FAIL]:
                        continue

                    df = self.comps[level]
                    mask = df["state"] == 0

                    if (mask.sum() / len(df)) > monitor_config[ck.FAIL_PER_THRESH]:
                        threshold_met = True
                        break

            if not threshold_met:
                continue

            # threshold met, run independent monitor
            # for levels with interval and threshold, you can reverse engineer the day by taking current day, adding the time to indep monitoring to it, then subtracting the monitor time for the components. then you can take the difference from current day and fail day to get the monitor time
            for level in monitor_config[ck.LEVELS]:
                if not self.case.config[level][ck.CAN_FAIL]:
                    continue

                df = self.comps[level]
                mask = (
                    (df["state"] == 0)
                    & ((df["time_to_detection"] >= 1) | (df["time_to_detection"].isna()))
                    & (df["indep_monitor"] == False)
                )
                monitor_comps = df.loc[mask].copy()

                if len(monitor_comps) == 0:
                    continue

                if monitor_comps["time_to_detection"].isna().sum() != 0:
                    if ck.INTERVAL in monitor_config:
                        monitor_comps["monitor_times"] = day - (
                            (day + self.indep_monitoring[name]) - (monitor_comps["monitor_times"])
                        )
                    else:
                        monitor_comps["monitor_times"] = day - monitor_comps["monitor_times"]
                else:
                    monitor_comps["monitor_times"] -= monitor_comps["time_to_detection"]

                if ck.DIST in monitor_config:
                    monitor_comps["time_to_detection"] = sample(
                        monitor_config[ck.DIST], monitor_config[ck.PARAM], len(monitor_comps)
                    )
                else:
                    monitor_comps["time_to_detection"] = 0

                monitor_comps["indep_monitor"] = True

                df.loc[mask] = monitor_comps
                # add cost to total costs (split evenly among all the levels for this day)
                self.costs[level][day] += (monitor_config[ck.COST] + self.labor_rate * monitor_config[ck.LABOR]) / len(
                    monitor_config[ck.LEVELS]
                )

max_stars_repo_name: FSEC-Photovoltaics/pvrpm-lcoe

max_issues_repo_name: FSEC-Photovoltaics/pvrpm-lcoe

max_forks_repo_name: FSEC-Photovoltaics/pvrpm-lcoe


text:
# -*- coding: utf-8 -*-
'''
Generate test data from a model
-------------------------------
Code to generate synthetic test data from one of the models.
'''

from __future__ import print_function
import numpy as np
import argparse
import sys
from astropy.time import Time
import json

import lal
import lalsimulation as lalsim

from em_pe.models import model_dict
from em_pe.utils import calc_mej, calc_vej

parser = argparse.ArgumentParser(description='Generate synthetic data for PE tests')
parser.add_argument('--m', help='Name of model to use')
parser.add_argument('--out', help='Directory to output data files to')
parser.add_argument('--p', action='append', nargs=2, help='Parameter value. Must set each parameter value for the model being used, as a pair (e.g. "--p dist 40")')
parser.add_argument('--tmin', type=float, help='Minimum time (in days)')
parser.add_argument('--tmax', type=float, help='Maximum time (in days)')
parser.add_argument('--n', type=int, help='Number of data points per band')
parser.add_argument('--err', type=float, help='Error std. dev.')
parser.add_argument('--orientation', nargs=2, help='Orientation profile and angle')
parser.add_argument('--t0', type=float, default=0, help='Start time for event')
parser.add_argument('--time-format', default='gps', help='Time format for t0 (gps or mjd)')
parser.add_argument('--json', action='store_true', help='Export data in JSON format')
parser.add_argument('--bns-params', action='store_true', help='Calculate EM parameters using BNS parameters')
parser.add_argument('--sigma', type=float, default=0.0, help='Extra error estimate to be fit as a parameter')
args = parser.parse_args()

def convert_time(t0):
    t = Time(t0, format='gps')
    return t.mjd

delta_t = args.t0

if args.time_format == 'gps':
    delta_t = convert_time(delta_t)

### create a dictionary mapping parameter names to values
params = dict(args.p)
for p in params:
    params[p] = float(params[p])

### convert BNS parameters to EM parameters, if needed
### borrowed from EOSManager.py
def lambda_from_m(m):
    eos = lalsim.SimNeutronStarEOSByName("AP4")
    eos_fam = lalsim.CreateSimNeutronStarFamily(eos)
    if m<10**15:
        m=m*lal.MSUN_SI
    k2=lalsim.SimNeutronStarLoveNumberK2(m, eos_fam)
    r=lalsim.SimNeutronStarRadius(m, eos_fam)
    m=m*lal.G_SI/lal.C_SI**2
    lam=2./(3*lal.G_SI)*k2*r**5
    dimensionless_lam=lal.G_SI*lam*(1/m)**5
    return dimensionless_lam

if args.bns_params:
    lambda1 = lambda_from_m(params['m1'])
    lambda2 = lambda_from_m(params['m2'])
    params['mej'] = calc_mej(params['m1'], lambda1, params['m2'], lambda2)
    params['vej'] = calc_vej(params['m1'], lambda1, params['m2'], lambda2)

### initialize the model
if args.orientation is not None:
    profile = args.orientation[0]
    angle = float(args.orientation[1])
    model = model_dict['oriented'](args.m, profile)
    params['angle'] = angle
else:
    model = model_dict[args.m]()

t_bounds = [max(0.01, args.tmin), args.tmax]

### set model parameters
model.set_params(params, t_bounds)

### generate times
#tdays = np.logspace(np.log10(args.tmin), np.log10(args.tmax), args.n)
tdays = np.exp(np.random.uniform(np.log(args.tmin), np.log(args.tmax), args.n))
tdays = np.sort(tdays)

data_dict = {}

### generate the data
for band in model.bands:
    m, _ = model.evaluate(tdays, band)
    if 'dist' in params:
        dist = params['dist']
        m += 5*(np.log10(dist*1e6) - 1)
    data = np.empty((4, args.n))
    data[0] = tdays + delta_t
    data[2] = m + np.random.normal(loc=0.0, scale=args.err, size=args.n) + np.random.normal(loc=0.0, scale=args.sigma, size=args.n)# generate errors
    data[3] = np.ones(args.n) * args.err
    data_dict[band] = data

### if not json format, just save as internal format
if not args.json:
    for band in data_dict:
        data = data_dict[band]
        filename = args.out + band + '.txt'
        np.savetxt(filename, data.T)

### json
else:
    photometry = []
    json_dict = {
        'SyntheticEvent': {
            'photometry': photometry
        }
    }
    for band in data_dict:
        data = data_dict[band]
        for i in range(args.n):
            [t, _, m, m_err] = data[:,i]
            photometry.append({
                'band':band,
                'time':t,
                'magnitude':m,
                'e_magnitude':m_err
            })
    with open(args.out + 'SyntheticEvent.json', 'w') as f:
        json.dump(json_dict, f, indent=4)

### save true values

truths = []
for param in model.param_names:
    truths.append(params[param])
header = " ".join(model.param_names)
if 'dist' in params:
    truths.append(params['dist'])
    header += " dist"
filename = args.out + 'test_truths.txt'
np.savetxt(filename, truths, header=header)

reason: import numpy,from astropy


text:
import numpy as np
import pandas as pd
import os
import functions
import scipy.io
import matplotlib.pyplot as plt

# LOAD DATA----------------------------------------------------------------------------------
sport_folder = 'PoleVault/'
dataDir = "dataset_dtw/" + sport_folder

finalPaths = pd.read_csv('Paths/PoleVault/path_polevault.csv')
finalPaths = np.array(finalPaths)

# SET ARRAYS---------------------------------------------------------------------------------
names = []
fc7 = []
delete_videos = [34, 35, 39, 54, 75, 87, 96, 97, 112, 113, 136, 137, 141, 145]
hm = []

'''
i = 0
for file in sorted(os.listdir(dataDir)):
    if file[0] != '.':

        mats = []
        mats = (scipy.io.loadmat(dataDir + file))
        tmp_fc7 = mats['fc7']

        if len(tmp_fc7.shape) > 2:
            names.append(os.path.basename(file.title()).lower())
            if tmp_fc7.shape[0] > 1:
                tmp_fc7 = tmp_fc7[0]

            tmp_fc7 = np.squeeze(tmp_fc7)

            fc7.append(tmp_fc7)

            tmp_hm = np.squeeze(mats['heatmaps'])
            hm.append(tmp_hm)
        else:
            delete_videos.append(i)

        i = i+1

print(delete_videos)
'''

# LOAD PICKLES ANNOTATIONS AND PLOT ----------------------------------------------------------

pickles, g = functions.open_pickles_annots(sport_folder, delete_videos)
print(pickles[1])

milestone_polevault = ['pole_down', 'over_bar', 'end_jump']

real_pole_down = []
stretched_pole_down = []

real_over_bar = []
stretched_over_bar = []

real_end_jump = []
stretched_end_jump = []

print(finalPaths.shape[0])

for i in range(finalPaths.shape[0]):
    real_pole_down.append(pickles[i][milestone_polevault[0]][1])

check = 0
for i in range(finalPaths.shape[0]):
    for k in range(len(finalPaths[i])):
        if check == 0:
            if finalPaths[i][k] == real_pole_down[i]:
                index = k
                check = 1
    check = 0
    stretched_pole_down.append(index)

for i in range(finalPaths.shape[0]):
    real_over_bar.append(pickles[i][milestone_polevault[1]][1])

check = 0
for i in range(finalPaths.shape[0]):
    for k in range(len(finalPaths[i])):
        if check == 0:
            if finalPaths[i][k] == real_over_bar[i]:
                index = k
                check = 1
    check = 0
    stretched_over_bar.append(index)

for i in range(finalPaths.shape[0]):
    real_end_jump.append(pickles[i][milestone_polevault[2]][1])

check = 0
for i in range(finalPaths.shape[0]):
    for k in range(len(finalPaths[i])):
        if check == 0:
            if finalPaths[i][k] == real_end_jump[i]:
                index = k
                check = 1
    check = 0
    stretched_end_jump.append(index)

# NORMALIZE ----------------------------------------------------------------------

stretched_pole_down = [x / len(finalPaths[i]) for x in stretched_pole_down]
stretched_over_bar = [x / len(finalPaths[i]) for x in stretched_over_bar]
stretched_end_jump = [x / len(finalPaths[i]) for x in stretched_end_jump]

real_pole_down_norm = []
for i in range(len(real_pole_down)):
    real_pole_down_norm.append(real_pole_down[i] / max(finalPaths[i]))

real_over_bar_norm = []
for i in range(len(real_pole_down)):
    real_over_bar_norm.append(real_over_bar[i] / max(finalPaths[i]))

real_end_jump_norm = []
for i in range(len(real_pole_down)):
    real_end_jump_norm.append(real_end_jump[i] / max(finalPaths[i]))

# VARIANCE -----------------------------------------------------------------------
print('-----pole down variance:')
print(np.std(real_pole_down_norm))
print(np.std(stretched_pole_down))
print('-----over bar variance:')
print(np.std(real_over_bar_norm))
print(np.std(stretched_over_bar))
print('-----end_jump variance:')
print(np.std(real_end_jump_norm))
print(np.std(stretched_end_jump))

# PLOTS -------------------------------------------------------------------------
X = np.arange(len(real_pole_down))

plt.figure()
plt.title('pole_down')
plt.scatter(X, sorted(real_pole_down_norm))
plt.scatter(X, sorted(stretched_pole_down))
plt.show()

plt.figure()
plt.title('over_bar')
plt.scatter(X, sorted(real_over_bar_norm))
plt.scatter(X, sorted(stretched_over_bar))
plt.show()

plt.figure()
plt.title('end_jump')
plt.scatter(X, sorted(real_end_jump_norm))
plt.scatter(X, sorted(stretched_end_jump))
plt.show()
#
# plt.figure()
# plt.hist(real_pole_down_norm, bins=20)
# plt.title('real_pole_down_norm')
# plt.show()
# plt.figure()
# plt.hist(stretched_pole_down, bins=20)
# plt.title('stretched_pole_down_norm')
# plt.show()
#
# plt.figure()
# plt.hist(real_over_bar_norm, bins=20)
# plt.title('real_over_bar_norm')
# plt.show()
# plt.figure()
# plt.hist(stretched_over_bar, bins=20)
# plt.title('stretched_over_bar_norm')
# plt.show()
#
# plt.figure()
# plt.hist(real_end_jump_norm, bins=20)
# plt.title('real_end_jump_norm')
# plt.show()
# plt.figure()
# plt.hist(stretched_end_jump, bins=20)
# plt.title('stretched_end_jump_norm')
# plt.show()

#functions.global_show_frames(struct, names, samples, finalAnnots, sport_folder)

reason: import numpy,import scipy


text:
# Demo usage of persisting task results
# Each task result will be stored in a serialized form to disk

import numpy as np
import prefect
from prefect import task, Flow
from prefect import context as prectx
from prefect.engine.results import LocalResult

import cloudpickle

# NOTE: Only works if env variable is set!
# export PREFECT__FLOWS__CHECKPOINTING=true

#@task(result=LocalResult(dir='./localtask'))
@task
def sum(a, b):
    logger = prectx.get("logger")
    c = a + b
    logger.info(f"Result: {c}")
    return c

def example():
    with Flow('Example 002', result=LocalResult(dir='./checkpoint_results')) as flow:
        sum(*np.random.randint(1000, size=2))
        sum(*np.random.randint(1000, size=2))
        sum(*np.random.randint(1000, size=2))
        sum(*np.random.randint(1000, size=2))
        sum(*np.random.randint(1000, size=2))

    state = flow.run()
    # Weired thing here, one needs to extract the storage locations from the state object
    # So make sure to persist the state object for later analysis.
    task_result_locations = [state.result[t]._result.location for t in flow.get_tasks()]

    for idx, loc in enumerate(task_result_locations):
        with open(loc, 'rb') as f:
            result = cloudpickle.load(f)
            print(f"Task[{idx}] return: {result}")

if __name__ == '__main__':
    example()


text:
def one_way_ANOVA(df):
    # 列と行の長さを取得
    f1 = len(df.columns)
    n = len(df.index)

    # 因子の効果を求めた後、因子変動S1を求める
    f_mean = df.mean().mean()
    f1_mean = df.mean()
    f1_effect = f1_mean - f_mean
    S1 = ((f1_effect**2) * n).sum()

    # 誤差変動Seを求める
    df_error = df.copy() - f1_mean
    Se = (df_error**2).sum().sum()

    # 総変動Stを求める
    St = ((df-f_mean)**2).sum().sum()

    # 自由度dfを求める
    df1 = f1 - 1
    dfe = f1*(n-1)
    dft = df1 + dfe

    # 不偏分散Vを求める
    V1 = S1 / df1
    Ve = Se / dfe

    # F値とp値を求める
    F = V1 / Ve
    p = 1 - st.f.cdf(F, dfn=df1, dfd=dfe)

    # 分散分析表を作成する
    df_S = pd.Series([S1, Se, St])
    df_df = pd.Series([df1, dfe, dft])
    df_V = pd.Series([V1, Ve])
    df_F = pd.Series([F])
    df_p = pd.DataFrame([p], columns=['p'])
    df_p['sign'] = df_p['p'].apply(lambda x : '**' if x < 0.01 else '*' if x < 0.05 else '')
    df_ANOVA = pd.concat([df_S, df_df, df_V, df_F, df_p], axis=1).set_axis(['S','df','V','F','p','sign'], axis=1).set_axis(['Columns','Error', 'Total']).rename_axis('ANOVA_table', axis=1).fillna('')

    # 因子の効果をデータフレームにまとめる
    df_effect = pd.DataFrame(f1_effect).set_axis(['mean'], axis=1).rename_axis('Effect', axis=1)
    df_effect['mean(95%CL)'] = df_effect['mean'].map(lambda x : st.t.interval(0.95, dfe, loc=x, scale=np.sqrt(Ve/n)))

    # 各水準毎の取得データの予測値をデータフレームにまとめる
    df_prediction = pd.DataFrame(df_effect['mean'].map(lambda x : st.t.interval(0.95, dfe, loc=x+f_mean, scale=np.sqrt(Ve/n + Ve)))).T.set_axis(['value']).rename_axis('Prediction(95%CL)', axis=1)

    # 結果を出力する
    return df_ANOVA, df_effect.T, df_prediction

# テスト
import numpy as np
import pandas as pd
import scipy.stats as st

df_1 = pd.DataFrame([[8,10,12], [4,8,12]]).set_axis(['column0', 'column1', 'column2'], axis=1)
print("入力：")
display(df_1)

df_ANOVA, df_effect, df_prediction = one_way_ANOVA(df_1)
print("出力：")
display(df_ANOVA, df_effect, df_prediction)


text:
from os.path import join as opjoin

import numpy as np

from metrics import FROCMetricsCalculator, draw_several, save_csv, MetricsCalculator
from roc_eval import SimpleTest, PatientTest, run_test

def label2model(wa, cv_stage, epoch):
    prefix = 'with_augmented' if wa else 'baseline'
    return {
        'path2model': '{}_128_cv_{}/detector_{}.ckpt'.format(prefix, cv_stage, str(epoch).zfill(3))
    }

def str_label(wa, cv_stage, epoch):
    prefix = 'with_augmented' if wa else 'baseline'
    return '{}_cv_{}_epoch_{}'.format(prefix, cv_stage, str(epoch).zfill(3))

class TestWrapper:
    def __init__(self, wa, cv_stage, epoch, mode='froc', dtype='crop'):
        self.wa = wa
        self.cv_stage = cv_stage
        self.epoch = epoch
        self.label = str_label(wa, cv_stage, epoch)
        self.mode = mode
        self.dtype = dtype

    def run(self, r_rand=0.5, stage=0):
        path2model = label2model(self.wa, self.cv_stage, self.epoch)['path2model']
        if self.dtype == 'crop':
            self.test = SimpleTest(paths2model=[path2model], r_rand=r_rand, stage=stage)
        else:
            self.test = PatientTest(paths2model=[path2model], r_rand=r_rand, stage=stage)

    def eval_metrics(self, left=-1, thr_number=10):
        result = run_test(self.test, mode=self.mode, left=left, thr_number=thr_number)
        self.mc = FROCMetricsCalculator(result, label=self.label) if self.mode == 'froc' else MetricsCalculator(result)

    def load(self, subfolder=None):
        self.mc = FROCMetricsCalculator(label=self.label)
        subfolder = '' if subfolder is None else subfolder
        load_path = opjoin(subfolder, self.label)
        self.mc.load(load_path + '.npy')

    def save(self, subfolder):
        self.mc.save(subfolder)

class DoubleTestWrapper:
    def __init__(self, cv_stage, baseline_epoch, augmented_epoch, all_tta=False):
        self.cv_stage = cv_stage
        self.baseline_epoch = baseline_epoch
        self.augmented_epoch = augmented_epoch
        self.baseline_label = str_label(False, cv_stage, baseline_epoch)
        self.augmented_label = str_label(True, cv_stage, augmented_epoch)
        self.all_tta = all_tta

    def run(self, r_rand=0.5, stage=0):
        path2baseline = label2model(False, self.cv_stage, self.baseline_epoch)['path2model']
        path2augmented = label2model(True, self.cv_stage, self.augmented_epoch)['path2model']
        self.test = SimpleTest(paths2model=[path2baseline, path2augmented], r_rand=r_rand, stage=stage, all_tta=self.all_tta)

    def eval_metrics(self):
        baseline_result = run_test(self.test, mode='froc', left=-1, thr_number=10, net_number=0)
        augmented_result = run_test(self.test, mode='froc', left=-1, thr_number=10, net_number=1)
        self.baseline_mc = FROCMetricsCalculator(baseline_result, label=self.baseline_label)
        self.augmented_mc = FROCMetricsCalculator(augmented_result, label=self.augmented_label)

    def load(self, subfolder=None):
        self.baseline_mc = FROCMetricsCalculator(label=self.baseline_label)
        subfolder = '' if subfolder is None else subfolder
        baseline_load_path = opjoin(subfolder, self.baseline_label)
        self.baseline_mc.load(baseline_load_path + '.npy')

        self.augmented_mc = FROCMetricsCalculator(label=self.augmented_label)
        aug_load_path = opjoin(subfolder, self.augmented_label)
        self.augmented_mc.load(aug_load_path + '.npy')

    def save(self, subfolder):
        self.baseline_mc.save(subfolder)
        self.augmented_mc.save(subfolder)

class TripleTestWrapper(DoubleTestWrapper):
    def add_baseline(self, path2model):
        self.path2second_baseline = path2model

    def run(self, r_rand=0.5, stage=0):
        path2baseline = label2model(False, self.cv_stage, self.baseline_epoch)['path2model']
        path2augmented = label2model(True, self.cv_stage, self.augmented_epoch)['path2model']
        self.test = SimpleTest(paths2model=[path2baseline, path2augmented, self.path2second_baseline], r_rand=r_rand, stage=stage)

    def eval_metrics(self):
        baseline_result = run_test(self.test, mode='froc', left=-1, thr_number=10, net_number=0)
        augmented_result = run_test(self.test, mode='froc', left=-1, thr_number=10, net_number=1)
        self.baseline_mc = FROCMetricsCalculator(baseline_result, label=self.baseline_label)
        self.augmented_mc = FROCMetricsCalculator(augmented_result, label=self.augmented_label)
        second_baseline_result = run_test(self.test, mode='froc', left=-1, thr_number=10, net_number=2)
        second_baseline_mc = FROCMetricsCalculator(second_baseline_result, label='label')
        self.baseline_mc.roc_result = np.add(self.baseline_mc.roc_result, second_baseline_mc.roc_result)

class CrossValidation:
    def __init__(self, subfolder, all_tta=False):
        self.stage2tw = {}
        self.subfolder = subfolder
        self.all_tta = all_tta

    def init_stage(self, number, baseline_epoch, wa_epoch):
        tw = DoubleTestWrapper(number, baseline_epoch, wa_epoch, self.all_tta)
        self.stage2tw[number] = tw

    def init_triple_stage(self, number, baseline_epoch, wa_epoch,  path2model):
        tw = TripleTestWrapper(number, baseline_epoch, wa_epoch, self.all_tta)
        tw.add_baseline(path2model)
        self.stage2tw[number] = tw

    def run_stage(self, number, r_rand=0.5):
        self.stage2tw[number].run(r_rand, number)

    def eval_stage_metrics(self, number):
        self.stage2tw[number].eval_metrics()

    def save_stage(self, number):
        self.stage2tw[number].save(self.subfolder)

    def load_stage(self, number):
        self.stage2tw[number].load(self.subfolder)

    def run_stages(self):
        for i in range(5):
            self.run_stage(i)

    def save_stages(self):
        for i in range(5):
            self.save_stage(i)

    def load_stages(self):
        for i in range(5):
            self.load_stage(i)

    def draw(self, stage_n):
        tw = self.stage2tw[stage_n]
        draw_several([tw.baseline_mc, tw.augmented_mc])

    def save_csv(self, filename):
        mcs = []
        for tw in self.stage2tw.values():
            mcs.append(tw.baseline_mc)
            mcs.append(tw.augmented_mc)
        save_csv(mcs, filename)

    def rename_labels(self):
        for i in range(5):
            self.stage2tw[i].baseline_mc.label = 'baseline_{}'.format(i + 1)
            self.stage2tw[i].augmented_mc.label = 'augmented_{}'.format(i + 1)


text:
"""
cluster parcels based on diffusion
"""

import numpy
import igraph
import os
import scipy.stats
import sklearn.cluster
import nibabel.gifti.giftiio
from consensus_cluster import pair_consensus
from run_shell_cmd import run_shell_cmd

basedir='/Users/poldrack/Dropbox/data/selftracking/DTI'
data=numpy.loadtxt('/Users/poldrack/Dropbox/data/selftracking/DTI/tracksumm_distcorr.txt')
data=data[:620,:620]

method='infomap'

size_thresh=8
densities=numpy.arange(0.00005,0.03,0.005)
assignments=numpy.zeros((620,len(densities)))

for dens in range(len(densities)):
    edge_density=densities[dens]
    cc=data.copy()
    cc[numpy.tril_indices(620)]=0
    edlabel='_%0.2f'%edge_density
    pct=scipy.stats.scoreatpercentile(cc[numpy.triu_indices(620,1)],100.0 - 100.0*edge_density)

    cc[cc<pct]=0
    cc[cc>0]=1
    gg=igraph.Graph.Adjacency(cc.tolist(),mode='UNDIRECTED')
    modclust_infomap=gg.community_infomap() #edge_weights='weight')
    print modclust_infomap.modularity
    assgn=numpy.array([i+1 for i in modclust_infomap.membership])  # zero treated as unclustered
    # exclude small modules by setting to zero
    for c in numpy.unique(assgn):
        if numpy.sum(assgn==c)<size_thresh:
            assgn[assgn==c]=0
    print 'density %0.3f (%0.3f): %d clusters'%(densities[dens],pct,len(numpy.unique(assgn)))
    assignments[:,dens]=assgn

numpy.savetxt('assignments.txt',assignments)

#cons=pair_consensus(assignments.T.astype('int'))
run_shell_cmd('/Applications/MATLAB_R2015a.app/bin/matlab -nosplash -nodesktop < Consensus_infomap_simple.m')

# modules start at zero
# relabel all single-node modules as zero and ignore them

labels=numpy.loadtxt('assignments_minsize5_regularized_consensus.txt')

SHUFFLE=False
if SHUFFLE:
    numpy.random.shuffle(labels)
# load parcel info
lines=open('/Users/poldrack/Dropbox/data/selftracking/rsfmri/module_assignments.txt').readlines()
modules=numpy.array([float(l.strip().split()[3]) for l in lines])
unique_modules=numpy.append(numpy.unique(modules)[1:],100)  # drop the zeros
module_hist_cutoffs=unique_modules-0.01

# get task distribution for each module

network_names=['1 Default','2 Second Visual','3 Frontal-Parietal','4.5 First Visual (V1+)','5 First Dorsal Attention','6 Second Dorsal Attention','7 Ventral Attention/Language','8 Salience','9 Cingulo-opercular','10 Somatomotor','11.5 Frontal-Parietal Other','15 Parietal Episodic Retrieval','16 Parieto-Occipital']

for l in numpy.unique(labels):
    if l < 1:
        continue
    members=labels==l
    networks=modules[members]
    networks=networks[networks>0]
    h=numpy.histogram(networks,bins=module_hist_cutoffs)[0]
    maxnet=numpy.where(h==numpy.max(h))[0][0]
    print 'Module',l,numpy.sum(members),'members'
    meanzstat=numpy.mean(data[:,members],1)
    idx=numpy.argsort(meanzstat)[::-1]
    #print 'Top contrasts:'
    #for i in idx[:5]:
    #    print copenames[i]
    hg=scipy.stats.hypergeom.sf(h[maxnet]-1,620,numpy.sum(modules==unique_modules[maxnet]),numpy.sum(members))
    if numpy.isnan(hg):
        asdf
    print 'Top network match: %s (pct match=%0.2f, dice = %0.3f, p < %0.3f)'%(network_names[maxnet],
        h[maxnet]/float(numpy.sum(members)),
        2*numpy.sum(modules==unique_modules[maxnet] * members)/float(numpy.sum(modules==unique_modules[maxnet]) + numpy.sum(members)),
        hg)

    print ''

# make gii files showing all modules

import nibabel.gifti.giftiio

lh=nibabel.gifti.giftiio.read('/Users/poldrack/Dropbox/data/selftracking/parcellation/all_selected_L_parcel_renumbered.func.gii')
rh=nibabel.gifti.giftiio.read('/Users/poldrack/Dropbox/data/selftracking/parcellation/all_selected_R_parcel_renumbered.func.gii')

lh_labels=lh.darrays[0].data.copy()
rh_labels=rh.darrays[0].data.copy()
labeled_regions_lh=numpy.unique(lh_labels)
labeled_regions_rh=numpy.unique(rh_labels)

lhimg=nibabel.gifti.GiftiImage()

rhimg=nibabel.gifti.GiftiImage()

labelvals=numpy.unique(labels)
modules_lh=labels[:310]
modules_rh=labels[310:]

for i in labelvals:
    if i < 1:
        continue
    darray=numpy.zeros(lh.darrays[0].data.shape,dtype=numpy.float32)
    matches=numpy.where(modules_lh==i)[0]+1
    for m in matches:
        darray[lh_labels==m]=1

    lhimg.add_gifti_data_array(nibabel.gifti.GiftiDataArray.from_array(darray,intent=lh.darrays[0].intent,datatype=lh.darrays[0].datatype,ordering='F',meta={'AnatomicalStructurePrimary':'CortexLeft','Name':'module%02d'%i}))

    darray=numpy.zeros(rh.darrays[0].data.shape,dtype=numpy.float32)
    matches=numpy.where(modules_rh==i)[0]+311
    for m in matches:
        darray[rh_labels==m]=1

    rhimg.add_gifti_data_array(nibabel.gifti.GiftiDataArray.from_array(darray,intent=rh.darrays[0].intent,datatype=rh.darrays[0].datatype,ordering='F',meta={'AnatomicalStructurePrimary':'CortexRight','Name':'module%02d'%i}))

if not SHUFFLE:
    nibabel.gifti.giftiio.write(lhimg,os.path.join(basedir,'lh_%s_consensus.func.gii'%method))
    nibabel.gifti.giftiio.write(rhimg,os.path.join(basedir,'rh_%s_consensus.func.gii'%method))

max_stars_repo_name: poldrack/myconnectome

max_issues_repo_name: poldrack/myconnectome

max_forks_repo_name: poldrack/myconnectome


text:
import numpy as np

class DefaultOrder:
    def __init__(self, order):
        self.order_map = {}
        for idx, elem in enumerate(order):
            self.order_map[elem] = idx
    def __call__(self, elem):
        if elem not in self.order_map:
            self.order_map[elem] = len(self.order_map)
        return self.order_map[elem]

def find_metric_name(automl):
    try:
        statistics = automl.sprint_statistics()
        lines = statistics.split("\n")
        my_line = list(filter(lambda x: x.upper().find("METRIC") > -1, lines))[0]
        metric = my_line.split(":")[1].strip().upper()
        return metric
    except Exception as e:
        return "METRIC"

def find_ensemble_weights(all_params, models_with_weights):
    weights = np.zeros(len(all_params))
    try:
      for idx, param in enumerate(all_params):
          for weight, pipeline in models_with_weights:
              weighted_param = pipeline.config.get_dictionary()
              if param == weighted_param:
                  weights[idx] = weight
                  break
    except Exception:
        pass
    return weights

def import_autosklearn(automl, source='auto-sklearn'):
    cv_results = automl.cv_results_
    weights = find_ensemble_weights(automl.cv_results_['params'], automl.get_models_with_weights())
    node_order = DefaultOrder(['balancing', 'data_preprocessing', 'feature_preprocessor','classifier', 'regressor'])
    n_models = len(cv_results['mean_test_score'])
    pipelines = []
    metric_name = find_metric_name(automl)
    for i in range(n_models):
        # Computing auxiliary pipeline structure
        struct = {}
        for key in cv_results['params'][i]:
            split = key.split(":")
            if split[0] == 'balancing' and split[1] == 'strategy':
                if cv_results['params'][i][key] != 'none':
                    struct['balancing'] = {'class_balancing': {'strategy': cv_results['params'][i][key]}}
                continue
            if split[1] == '__choice__':
                continue
            module_type = split[0]
            if module_type not in struct:
                struct[module_type] = {}
            module_name = split[1]
            if module_name not in struct[module_type]:
                struct[module_type][module_name] = {}
            module_param = '_'.join(filter(lambda x: x != '__choice__', split[2:]))
            struct[module_type][module_name][module_param] = cv_results['params'][i][key]
        ordered_types = struct.keys()
        ordered_types = sorted(ordered_types, key=lambda x: node_order(x))

        # Creating pipeline graph
        pipeline = {
            'inputs': [{'name': 'input dataset'}],
            'steps': [],
            'scores': [{
                'metric': {'metric': metric_name, 'params': {'pos_label': '1'}},
                'normalized': cv_results['mean_test_score'][i],
                'value': cv_results['mean_test_score'][i],
            },
            {
                'metric': {'metric': 'ENSEMBLE WEIGHT', 'params': {'pos_label': '2'}},
                'normalized': weights[i],
                'value': weights[i]
            },{
                'metric': {'metric': 'MEAN FIT TIME', 'params': {'pos_label': '3'}},
                'normalized': cv_results['mean_fit_time'][i],
                'value': cv_results['mean_fit_time'][i],
            }],
            'pipeline_source': {'name': source},
            'pipeline_digest': '{}'.format(i),
        }
        prev_list = ['inputs.0']
        cur_step_idx = 0
        for module_type in ordered_types:
            new_prev_list = []
            modules = struct[module_type]
            for module in modules:
                step_ref = 'steps.{}.produce'.format(cur_step_idx)
                cur_step_idx += 1
                new_prev_list.append(step_ref)
                step = {
                    'primitive': {'python_path':'auto_sklearn.primitives.{}.{}'.format(module_type, module), 'name': module},
                    'arguments':{},
                    'outputs': [{'id': 'produce'}],
                    'reference': {'type': 'CONTAINER', 'data': step_ref},
                    'hyperparams': {}
                }
                for param in struct[module_type][module]:
                    step['hyperparams'][param] = {
                        'type': 'VALUE',
                        'data': struct[module_type][module][param]
                    }
                for idx, prev in enumerate(prev_list):
                    cur_argument_idx = 'input{}'.format(idx)
                    step['arguments'][cur_argument_idx] = {
                        'data': prev
                    }
                pipeline['steps'].append(step)
            prev_list = new_prev_list
        pipeline['outputs'] = []
        for prev in prev_list:
            pipeline['outputs'].append({'data': prev})
        pipelines.append(pipeline)
    return pipelines

def import_autosklearn_test_data(automl, test_X, test_y, metric_eval, source='auto-sklearn'):
    node_order = DefaultOrder(['balancing', 'data_preprocessing', 'feature_preprocessor','classifier', 'regressor'])
    models_with_weights = automl.get_models_with_weights()
    weights = []
    models = []
    for weight, model in models_with_weights:
        weights.append(weight)
        models.append(model)
    n_models = len(models)
    pipelines = []
    for i in range(n_models):
        # Computing auxiliary pipeline structure
        config = models[i].config.get_dictionary()
        struct = {}
        for key in config:
            split = key.split(":")
            if split[0] == 'balancing' and split[1] == 'strategy':
                if config[key] != 'none':
                    struct['balancing'] = {'class_balancing': {'strategy': config[key]}}
                continue
            if split[1] == '__choice__':
                continue
            module_type = split[0]
            if module_type not in struct:
                struct[module_type] = {}
            module_name = split[1]
            if module_name not in struct[module_type]:
                struct[module_type][module_name] = {}
            module_param = '_'.join(filter(lambda x: x != '__choice__', split[2:]))
            struct[module_type][module_name][module_param] = config[key]
        ordered_types = struct.keys()
        ordered_types = sorted(ordered_types, key=lambda x: node_order(x))

        # Computing score
        pred_y = models[i].predict(test_X)
        score = metric_eval(test_y, pred_y)

        # Creating pipeline graph
        pipeline = {
            'inputs': [{'name': 'input dataset'}],
            'steps': [],
            'scores': [{
                'metric': {'metric': str(metric_eval).capitalize(), 'params': {'pos_label': '1'}},
                'normalized': score,
                'value': score,
            },
            {
                'metric': {'metric': 'Ensemble Weight', 'params': {'pos_label': '2'}},
                'normalized': weights[i],
                'value': weights[i]
            }],
            'pipeline_source': {'name': source},
            'pipeline_digest': '{}'.format(i),
        }
        prev_list = ['inputs.0']
        cur_step_idx = 0
        for module_type in ordered_types:
            new_prev_list = []
            modules = struct[module_type]
            for module in modules:
                step_ref = 'steps.{}.produce'.format(cur_step_idx)
                cur_step_idx += 1
                new_prev_list.append(step_ref)
                step = {
                    'primitive': {'python_path':'auto_sklearn.primitives.{}.{}'.format(module_type, module), 'name': module},
                    'arguments':{},
                    'outputs': [{'id': 'produce'}],
                    'reference': {'type': 'CONTAINER', 'data': step_ref},
                    'hyperparams': {}
                }
                for param in struct[module_type][module]:
                    step['hyperparams'][param] = {
                        'type': 'VALUE',
                        'data': struct[module_type][module][param]
                    }
                for idx, prev in enumerate(prev_list):
                    cur_argument_idx = 'input{}'.format(idx)
                    step['arguments'][cur_argument_idx] = {
                        'data': prev
                    }
                pipeline['steps'].append(step)
            prev_list = new_prev_list
        pipeline['outputs'] = []
        for prev in prev_list:
            pipeline['outputs'].append({'data': prev})
        pipelines.append(pipeline)
    return pipelines

max_stars_repo_name: VIDA-NYU/PipelineVis

max_issues_repo_name: VIDA-NYU/PipelineVis

max_forks_repo_name: VIDA-NYU/PipelineVis


text:
import pytest

from skimpy.nullspace import left_integer_nullspace

import numpy as np
from scipy.sparse import random
from scipy import stats

class ThisCustomRandomState(np.random.RandomState):

    def randint(self, k):
        i = np.random.randint(k)
        return i

    def choice(self, mn, size, replace):
        return np.random.choice(mn, size=size, replace=replace)

def test_left_nullspace():

    rs = ThisCustomRandomState()
    rvs = stats.poisson(2, loc=10).rvs
    S = random(5,6, density=0.1, random_state=rs, data_rvs=rvs, dtype=np.int)

    # print(S.todense())

    ns = left_integer_nullspace(S.todense())

    # print(ns)
    null = ns @ S.todense()

    # print(null)
    # assert(np.any(null) == False)
    assert(null.max() < 1e-12)

reason: import numpy,from scipy


text:
"""Test visserver."""

import os
import tempfile

import numpy as np
import pytest

import pyabc
import pyabc.visserver.server as server

db_path = "sqlite:///" + tempfile.mkstemp(suffix='.db')[1]

def setup_module():
    """Run an analysis and create a database.

    Called once at the beginning.
    """

    def model(p):
        return {
            'ss0': p['p0'] + 0.1 * np.random.uniform(),
            'ss1': p['p1'] + 0.1 * np.random.uniform(),
        }

    p_true = {'p0': 3, 'p1': 4}
    observation = {'ss0': p_true['p0'], 'ss1': p_true['p1']}
    limits = {'p0': (0, 5), 'p1': (1, 8)}
    prior = pyabc.Distribution(
        **{
            key: pyabc.RV(
                'uniform', limits[key][0], limits[key][1] - limits[key][0]
            )
            for key in p_true.keys()
        }
    )
    distance = pyabc.PNormDistance(p=2)

    abc = pyabc.ABCSMC(model, prior, distance, population_size=50)
    abc.new(db_path, observation)
    abc.run(minimum_epsilon=0.1, max_nr_populations=4)

def teardown_module():
    """Remove history.

    Called once after all tests.
    """
    os.remove(db_path[len("sqlite:///") :])

@pytest.fixture
def client():
    """A fake server client."""
    history = pyabc.History(db_path)
    server.app.config["HISTORY"] = history
    with server.app.test_client() as client:
        yield client

def test_visserver_basic(client):
    """Test that the visserver does at least something."""
    ok_requests = [
        '/',
        '/abc',
        '/abc/1',
        '/abc/1/model/0/t/3',
        '/info',
    ]
    for request in ok_requests:
        rv = client.get(request)
        assert rv.status == '200 OK'

    internal_server_error_requests = [
        '/abc/2',
    ]
    for request in internal_server_error_requests:
        rv = client.get(request)
        assert rv.status == '500 INTERNAL SERVER ERROR'

    not_found_requests = [
        '/britney',
    ]
    for request in not_found_requests:
        rv = client.get(request)
        assert rv.status == '404 NOT FOUND'

max_stars_repo_name: Gabriel-p/pyABC

max_issues_repo_name: Gabriel-p/pyABC

max_forks_repo_name: Gabriel-p/pyABC


text:
import os
import numpy as np
from soundfile import write

def transpose(samples):
    return np.array(samples).transpose(1, 0)

data_sin = np.sin(np.array(range(0, 360, 15)) * np.pi / 180)
data_sins = transpose([data_sin] * 2)

def path(filename):
    return os.path.join('fuzzing_audio_files', filename)

def write_float_sin():
    write(path('sin_441_float.wav'), data_sin, 44100, subtype='FLOAT')
    write(path('sins_441_float.wav'), data_sins, 44100, subtype='FLOAT')
    write(path('sin_96_float.wav'), data_sin, 96000, subtype='FLOAT')
    write(path('sins_96_float.wav'), data_sins, 96000, subtype='FLOAT')

def write_PCM16_sin():
    write(path('sin_441_pcm16.wav'), data_sin, 44100, subtype='PCM_16')
    write(path('sins_441_pcm16.wav'), data_sins, 44100, subtype='PCM_16')
    write(path('sin_96_pcm16.wav'), data_sin, 96000, subtype='PCM_16')
    write(path('sins_96_pcm16.wav'), data_sins, 96000, subtype='PCM_16')

def write_PCM24_sin():
    write(path('sin_441_pcm24.wav'), data_sin, 44100, subtype='PCM_24')
    write(path('sins_441_pcm24.wav'), data_sins, 44100, subtype='PCM_24')
    write(path('sin_96_pcm24.wav'), data_sin, 96000, subtype='PCM_24')
    write(path('sins_96_pcm24.wav'), data_sins, 96000, subtype='PCM_24')

def write_long_sin():
    write(path('sin_long_PCM_16.wav'), np.tile(data_sin, 65536), 44100, subtype='PCM_16')

def write_multichannel_sins():
    write(path('sins_7ch.wav'), transpose([data_sin]*7), 44100, subtype='PCM_32')

def main():
    write_float_sin()
    write_PCM16_sin()
    write_PCM24_sin()
    write_long_sin()
    write_multichannel_sins()

if __name__ == "__main__":
    main()

max_stars_repo_name: RCJacH/monolizer

max_issues_repo_name: RCJacH/monolizer

max_forks_repo_name: RCJacH/monolizer


text:
import numpy as np

def txt2arr(path: str, delimiter: str = ' ') -> np.array:
    with open(path, 'r', encoding='utf-8') as f:
        mat = f.read()
    row_list = mat.splitlines()
    data_list = [[float(i)
                  for i in row.strip().split(delimiter)]
                 for row in row_list]
    return np.array(data_list)

def RMSE(x: np.array, y: np.array) -> float:
    return (np.sum((x - y)**2))**0.5 / len(x)

if __name__ == '__main__':
    mat1 = txt2arr('mat-A-2048.txt')
    mat2 = txt2arr('mat-B-2048.txt')
    my_prod = txt2arr('out2048.txt')
    print(np.dot(mat1, mat2))  # to check if the answer is generally correct
    print('RMSE = ', RMSE(np.dot(mat1, mat2), my_prod))

max_stars_repo_path: Project 2/src/check.py

max_issues_repo_path: Project 2/src/check.py

max_forks_repo_path: Project 2/src/check.py


text:
import gym
import numpy as np
import pybullet as p

from igibson.external.pybullet_tools.utils import joints_from_names
from igibson.robots.robot_locomotor import LocomotorRobot

class JR2_Kinova(LocomotorRobot):
    """
    JR2 Kinova robot
    Reference: https://cvgl.stanford.edu/projects/jackrabbot/
    Uses joint velocity control
    """

    def __init__(self, config):
        self.config = config
        self.wheel_velocity = config.get("wheel_velocity", 0.3)
        self.wheel_dim = 2
        self.arm_velocity = config.get("arm_velocity", 1.0)
        self.arm_dim = 5

        LocomotorRobot.__init__(
            self,
            "jr2_urdf/jr2_kinova.urdf",
            action_dim=self.wheel_dim + self.arm_dim,
            scale=config.get("robot_scale", 1.0),
            is_discrete=config.get("is_discrete", False),
            control="velocity",
            self_collision=True,
        )

    def set_up_continuous_action_space(self):
        """
        Set up continuous action space
        """
        self.action_high = np.array([self.wheel_velocity] * self.wheel_dim + [self.arm_velocity] * self.arm_dim)
        self.action_low = -self.action_high
        self.action_space = gym.spaces.Box(shape=(self.wheel_dim + self.arm_dim,), low=-1.0, high=1.0, dtype=np.float32)

    def set_up_discrete_action_space(self):
        """
        Set up discrete action space
        """
        assert False, "JR2_Kinova does not support discrete actions"

    def get_end_effector_position(self):
        """
        Get end-effector position
        """
        return self.parts["m1n6s200_end_effector"].get_position()

    def robot_specific_reset(self):
        """
        JR2 Kinova robot specific reset.
        Initialize JR's arm to about the same height at its neck, facing forward
        """
        super(JR2_Kinova, self).robot_specific_reset()
        self.ordered_joints[2].reset_joint_state(-np.pi / 2.0, 0.0)
        self.ordered_joints[3].reset_joint_state(np.pi / 2.0, 0.0)
        self.ordered_joints[4].reset_joint_state(np.pi / 2.0, 0.0)
        self.ordered_joints[5].reset_joint_state(np.pi / 2.0, 0.0)
        self.ordered_joints[6].reset_joint_state(0.0, 0.0)

    def load(self):
        """
        Load the robot into pybullet. Filter out unnecessary self collision
        due to modeling imperfection in the URDF
        """
        ids = super(JR2_Kinova, self).load()
        robot_id = self.robot_ids[0]

        disable_collision_names = [
            ["base_chassis_joint", "pan_joint"],
            ["base_chassis_joint", "tilt_joint"],
            ["base_chassis_joint", "camera_joint"],
            ["jr2_fixed_body_joint", "pan_joint"],
            ["jr2_fixed_body_joint", "tilt_joint"],
            ["jr2_fixed_body_joint", "camera_joint"],
        ]
        for names in disable_collision_names:
            link_a, link_b = joints_from_names(robot_id, names)
            p.setCollisionFilterPair(robot_id, robot_id, link_a, link_b, 0)

        return ids

max_stars_repo_name: suresh-guttikonda/iGibson

max_issues_repo_name: suresh-guttikonda/iGibson

max_forks_repo_name: suresh-guttikonda/iGibson