text:
# Generic/Built-in
import datetime
import math
import os
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import pvlib
from dataclasses_json import dataclass_json
from typing import Optional

from dataclasses import dataclass
from functools import lru_cache
from hisim.simulationparameters import SimulationParameters
# Owned
from hisim import component as cp
from hisim import loadtypes as lt
from hisim import utils
from hisim import log
from hisim.components.weather import Weather

__authors__ = "Vitor Hugo Bellotto Zago"
__copyright__ = "Copyright 2021, the House Infrastructure Project"
__credits__ = ["Noah Pflugradt"]
__license__ = "MIT"
__version__ = "0.1"
__maintainer__ = "Vitor Hugo Bellotto Zago"
__email__ = "vitor.zago@rwth-aachen.de"
__status__ = "development"

"""
The functions cited in this module are at some degree based on the tsib project:

[tsib-kotzur]:
Kotzur, Leander, Detlef Stolten, and Hermann-Josef Wagner. Future grid load of the residential building sector. No. RWTH-2018-231872. Lehrstuhl für Brennstoffzellen (FZ Jülich), 2019.
ID: http://hdl.handle.net/2128/21115
    http://nbn-resolving.org/resolver?verb=redirect&identifier=urn:nbn:de:0001-2019020614

The implementation of the tsib project can be found under the following repository:
https://github.com/FZJ-IEK3-VSA/tsib
"""

temp_model = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS["sapm"]["open_rack_glass_glass"]

@lru_cache(maxsize=16)
def simPhotovoltaicFast(
    dni_extra=None,
    DNI=None,
    DHI=None,
    GHI=None,
    azimuth=None,
    apparent_zenith=None,
    temperature=None,
    wind_speed=None,
    surface_azimuth : float = 180,
    surface_tilt : float = 30 ):

    """
    Simulates a defined PV array with the Sandia PV Array Performance Model.
    The implementation is done in accordance with following tutorial:
    https://github.com/pvlib/pvlib-python/blob/master/docs/tutorials/tmy_to_power.ipynb

    Parameters
    ----------
    surface_tilt: int or float, optional (default:30)
        Tilt angle of of the array in degree.
    surface_azimuth: int or float, optional (default:180)
        Azimuth angle of of the array in degree. 180 degree means south,
        90 degree east and 270 west.
    losses: float, optional (default: 0.1)
        Losses due to soiling, mismatch, diode connections, dc wiring etc.
    Returns
    --------
    """

    poa_irrad = pvlib.irradiance.get_total_irradiance( surface_tilt,
                                                        surface_azimuth,
                                                        apparent_zenith,
                                                        azimuth,
                                                        DNI,
                                                        GHI,
                                                        DHI,
                                                        dni_extra )

    pvtemps = pvlib.temperature.sapm_cell(poa_irrad["poa_global"], temperature, wind_speed, **temp_model)

    pv_dc = pvlib.pvsystem.pvwatts_dc( poa_irrad[ "poa_global" ],
                                      temp_cell = pvtemps,
                                      pdc0 = 1,
                                      gamma_pdc = -0.002,
                                      temp_ref = 25.0 )
    if math.isnan(pv_dc):
        pv_dc = 0
    return pv_dc

def simPhotovoltaicSimple(
    dni_extra=None,
    DNI=None,
    DHI=None,
    GHI=None,
    azimuth=None,
    apparent_zenith=None,
    temperature=None,
    wind_speed=None,
    surface_tilt=30,
    surface_azimuth=180,
    albedo=0.2):
    """
    Simulates a defined PV array with the Sandia PV Array Performance Model.
    The implementation is done in accordance with following tutorial:
    https://github.com/pvlib/pvlib-python/blob/master/docs/tutorials/tmy_to_power.ipynb

    Based on the tsib project @[tsib-kotzur] (Check header)

    Parameters
    ----------
    tmy_data: pandas.DataFrame(), required
        Weatherfile in the format of a tmy file.
    surface_tilt: int or float, optional (default:30)
        Tilt angle of of the array in degree.
    surface_azimuth: int or float, optional (default:180)
        Azimuth angle of of the array in degree. 180 degree means south,
        90 degree east and 270 west.
    albedo: float, optional (default: 0.2)
        Reflection coefficient of the surrounding area.
    losses: float, optional (default: 0.1)
        Losses due to soiling, mismatch, diode connections, dc wiring etc.
    load_module_data: Boolean, optional (default: False)
        If True the module data base is loaded from the Sandia Website.
        Otherwise it is loaded from this relative path
            '\\profiles\\PV-Modules\\sandia_modules.csv'.
    module_name: str, optional (default:'Hanwha_HSL60P6_PA_4_250T__2013_')
        Module name. The string must be existens in Sandia Module database.
    integrateInverter: bool, optional (default: True)
        If an inverter shall be added to the simulation, providing the photovoltaic output after the inverter.
    inverter_name: str, optional (default: 'ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_')
        Type of inverter.

    Returns
    --------
    """
    # automatic pd time series in future pvlib version
    # calculate airmass
    airmass = pvlib.atmosphere.get_relative_airmass(apparent_zenith)
    # use perez model to calculate the plane of array diffuse sky radiation
    poa_sky_diffuse = pvlib.irradiance.perez(
        surface_tilt,
        surface_azimuth,
        DHI,
        np.float64(DNI),
        dni_extra,
        apparent_zenith,
        azimuth,
        airmass,
    )
    # calculate ground diffuse with specified albedo
    poa_ground_diffuse = pvlib.irradiance.get_ground_diffuse(
        surface_tilt, GHI, albedo=albedo
    )
    # calculate angle of incidence
    aoi = pvlib.irradiance.aoi(surface_tilt, surface_azimuth, apparent_zenith, azimuth)
    # calculate plane of array irradiance
    poa_irrad = pvlib.irradiance.poa_components(aoi, np.float64(DNI), poa_sky_diffuse, poa_ground_diffuse)
    # calculate pv cell and module temperature
    temp_model = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS["sapm"]["open_rack_glass_glass"]
    pvtemps = pvlib.temperature.sapm_cell(poa_irrad["poa_global"], temperature, wind_speed, **temp_model)

    pv_dc = pvlib.pvsystem.pvwatts_dc(poa_irrad["poa_global"], temp_cell=pvtemps, pdc0=1, gamma_pdc=-0.002,
                                      temp_ref=25.0)
    if math.isnan(pv_dc):
        pv_dc = 0
    return pv_dc

@dataclass_json
@dataclass
class PVSystemConfig:
    parameter_string: str
    time: int
    location: str
    module_name:str
    integrate_inverter: bool
    inverter_name:str
    power: float

    def __init__(self,
                 my_simulation_parameters: SimulationParameters,
                 time:int,
                 location:str,
                 power:float,
                 module_name:str,
                 integrate_inverter:bool,
                 inverter_name:str ):
        self.parameter_string = my_simulation_parameters.get_unique_key()
        self.time = time
        self.location = location
        self.module_name = module_name
        self.integrate_inverter = integrate_inverter
        self.inverter_name = inverter_name
        self.power = power

class PVSystem(cp.Component):
    """
    Parameters:
    -----------------------------------------------------
    time:
        simulation timeline
    location: Location
        object Location with temperature and solar data
    power: float
        Power in kWp to be provided by the PV System

    Returns:
    -----------------------------------------------------
    pass
    """
    # Inputs
    TemperatureOutside = "TemperatureOutside"
    DirectNormalIrradiance = "DirectNormalIrradiance"
    DirectNormalIrradianceExtra = "DirectNormalIrradianceExtra"
    DiffuseHorizontalIrradiance = "DiffuseHorizontalIrradiance"
    GlobalHorizontalIrradiance = "GlobalHorizontalIrradiance"
    Azimuth = "Azimuth"
    ApparentZenith = "ApparentZenith"
    WindSpeed = "WindSpeed"

    # Outputs
    ElectricityOutput = "ElectricityOutput"

    #Forecasts
    PV_Forecast_24h = "PV_Forecast_24h"

    # Similar components to connect to:
    # 1. Weather
    @utils.measure_execution_time
    def __init__(self,
                 my_simulation_parameters: SimulationParameters,
                 my_simulation_repository : Optional[ cp.SimRepository ] = None,
                 time : int = 2019,
                 location : str = "Aachen",
                 power : float = 10E3,
                 load_module_data : bool = False,
                 module_name : str = "Hanwha_HSL60P6_PA_4_250T__2013_",
                 integrateInverter : bool = True,
                 inverter_name : str = "ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_" ):

        super().__init__( "PVSystem", my_simulation_parameters = my_simulation_parameters )
        self.pvconfig = PVSystemConfig(my_simulation_parameters=my_simulation_parameters,
                                       location=location, power = power, module_name=module_name,
                                       integrate_inverter=integrateInverter, inverter_name=inverter_name,
                                       time=time)
        self.build( load_module_data, my_simulation_repository )

        self.t_outC : cp.ComponentInput = self.add_input(self.ComponentName,
                                                        self.TemperatureOutside,
                                                        lt.LoadTypes.Temperature,
                                                        lt.Units.Celsius,
                                                        True)

        self.DNIC : cp.ComponentInput = self.add_input(self.ComponentName,
                                                    self.DirectNormalIrradiance,
                                                    lt.LoadTypes.Irradiance,
                                                    lt.Units.Wm2,
                                                    True)

        self.DNIextraC : cp.ComponentInput = self.add_input(self.ComponentName,
                                                         self.DirectNormalIrradianceExtra,
                                                         lt.LoadTypes.Irradiance,
                                                         lt.Units.Wm2,
                                                         True)

        self.DHIC: cp.ComponentInput = self.add_input(self.ComponentName,
                                                   self.DiffuseHorizontalIrradiance,
                                                   lt.LoadTypes.Irradiance,
                                                   lt.Units.Wm2,
                                                   True)

        self.GHIC: cp.ComponentInput = self.add_input(self.ComponentName,
                                                   self.GlobalHorizontalIrradiance,
                                                   lt.LoadTypes.Irradiance,
                                                   lt.Units.Wm2,
                                                   True)

        self.azimuthC : cp.ComponentInput = self.add_input(self.ComponentName,
                                                        self.Azimuth,
                                                        lt.LoadTypes.Any,
                                                        lt.Units.Degrees,
                                                        True)

        self.apparent_zenithC : cp.ComponentInput = self.add_input(self.ComponentName,
                                                                self.ApparentZenith,
                                                                lt.LoadTypes.Any,
                                                                lt.Units.Degrees,
                                                                True)

        self.wind_speedC: cp.ComponentInput = self.add_input(self.ComponentName,
                                                          self.WindSpeed,
                                                          lt.LoadTypes.Speed,
                                                          lt.Units.MeterPerSecond,
                                                          True)

        self.electricity_outputC : cp.ComponentOutput = self.add_output(self.ComponentName,
                                                             PVSystem.ElectricityOutput,
                                                             lt.LoadTypes.Electricity,
                                                             lt.Units.Watt,
                                                             False)

        self.add_default_connections(Weather, self.get_weather_default_connections())

    def get_weather_default_connections(self):
        log.information("setting weather default connections")
        connections = []
        weather_classname = Weather.get_classname()
        connections.append(cp.ComponentConnection(PVSystem.TemperatureOutside,weather_classname, Weather.TemperatureOutside))
        connections.append(cp.ComponentConnection(PVSystem.DirectNormalIrradiance,weather_classname, Weather.DirectNormalIrradiance))
        connections.append(cp.ComponentConnection(PVSystem.DirectNormalIrradianceExtra,weather_classname, Weather.DirectNormalIrradianceExtra))
        connections.append(cp.ComponentConnection(PVSystem.DiffuseHorizontalIrradiance,weather_classname, Weather.DiffuseHorizontalIrradiance))
        connections.append(cp.ComponentConnection(PVSystem.GlobalHorizontalIrradiance,weather_classname, Weather.GlobalHorizontalIrradiance))
        connections.append(cp.ComponentConnection(PVSystem.Azimuth,weather_classname, Weather.Azimuth))
        connections.append(cp.ComponentConnection(PVSystem.ApparentZenith,weather_classname, Weather.ApparentZenith))
        connections.append(cp.ComponentConnection(PVSystem.WindSpeed,weather_classname, Weather.WindSpeed))
        return connections

    def i_restore_state(self):
        pass

    def write_to_report(self):
        lines = []
        lines.append("Name: {}".format(self.ComponentName))
        lines.append("Power: {:3.0f} kWp".format(self.pvconfig.power*1E-3))
        lines.append("Module: {}".format(self.pvconfig.module_name))
        lines.append("Inverter: {}".format(self.pvconfig.inverter_name))
        return lines

    def i_simulate(self, timestep: int, stsv: cp.SingleTimeStepValues,  force_convergence: bool):

        if hasattr(self, "output"):
            #if(len(self.output) < timestep)
             #   raise Exception("Somehow the precalculated list of values for the PV system seems to be incorrect. Please delete the cache.")
            stsv.set_output_value(self.electricity_outputC, self.output[timestep] * self.pvconfig.power)
        else:
            DNI = stsv.get_input_value(self.DNIC)
            dni_extra = stsv.get_input_value(self.DNIextraC)
            DHI = stsv.get_input_value(self.DHIC)
            GHI = stsv.get_input_value(self.GHIC)
            azimuth = stsv.get_input_value(self.azimuthC)
            temperature = stsv.get_input_value(self.t_outC)
            wind_speed = stsv.get_input_value(self.wind_speedC)
            apparent_zenith = stsv.get_input_value(self.apparent_zenithC)

            #ac_power = self.simPhotovoltaic2(dni_extra=dni_extra,
            #                                 DNI=DNI,
            #                                 DHI=DHI,
            #                                 GHI=GHI,
            #                                 azimuth=azimuth,
            #                                 apparent_zenith=apparent_zenith,
            #                                 temperature=temperature,
            #                                 wind_speed=wind_speed)
            #ac_power = simPhotovoltaicSimple(
            #    dni_extra=dni_extra,
            #                                 DNI=DNI,
            #                                 DHI=DHI,
            #                                 GHI=GHI,
            #                                 azimuth=azimuth,
            #                                 apparent_zenith=apparent_zenith,
            #                                 temperature=temperature,
            #                                 wind_speed=wind_speed)

            ac_power = simPhotovoltaicFast(
                                            dni_extra=dni_extra,
                                            DNI=DNI,
                                            DHI=DHI,
                                            GHI=GHI,
                                            azimuth=azimuth,
                                            apparent_zenith=apparent_zenith,
                                            temperature=temperature,
                                            wind_speed=wind_speed)

            resultingvalue = ac_power * self.pvconfig.power
            # if you wanted to access the temperature forecast from the weather component:
            # val = self.simulation_repository.get_entry(Weather.Weather_Temperature_Forecast_24h)
            stsv.set_output_value(self.electricity_outputC, resultingvalue)
            self.data[timestep] = ac_power
            if timestep + 1 == self.data_length:
                database = pd.DataFrame(self.data, columns=["output"])

                database.to_csv(self.cache_filepath, sep=",", decimal=".", index=False)

        if self.my_simulation_parameters.system_config.predictive == True:
            last_forecast_timestep = int( timestep + 24 * 3600 / self.my_simulation_parameters.seconds_per_timestep )
            if ( last_forecast_timestep > len( self.output ) ):
                last_forecast_timestep = len( self.output )
            pvforecast = [ self.output[ t ] * self.pvconfig.power for t in range( timestep, last_forecast_timestep ) ]
            self.simulation_repository.set_entry( self.PV_Forecast_24h, pvforecast )

    def get_coordinates(self, location="Aachen", year=2019):
        """
        Reads a test reference year file and gets the GHI, DHI and DNI from it.

        Based on the tsib project @[tsib-kotzur] (Check header)

        Parameters
        -------
        try_num: int (default: 4)
            The region number of the test reference year.
        year: int (default: 2010)
            The year. Only data for 2010 and 2030 available
        """
        # get the correct file path
        filepath = os.path.join(utils.HISIMPATH["weather"][location])

        # get the geoposition
        with open(filepath + ".dat", encoding="utf-8") as fp:
            lines = fp.readlines()
            location_name = lines[0].split(maxsplit=2)[2].replace('\n', '')
            lat = float(lines[1][20:37])
            lon = float(lines[2][15:30])
        self.location = {"name": location_name, "latitude": lat, "longitude": lon}
        self.index = pd.date_range(
            "{}-01-01 00:00:00".format(year), periods=60*24*365, freq="T", tz="Europe/Berlin"
        )

    def i_save_state(self):
        pass

    def i_doublecheck(self, timestep: int, stsv: cp.SingleTimeStepValues):
        pass

    def build( self, load_module_data : bool, my_simulation_repository : Optional[ cp.SimRepository ] ):
        log.information(self.pvconfig.to_json())  # type: ignore
        file_exists, self.cache_filepath = utils.get_cache_file("PVSystem", self.pvconfig)

        if file_exists:
            self.output = pd.read_csv(self.cache_filepath, sep=',', decimal='.')['output'].tolist()
            if len(self.output) !=        self.my_simulation_parameters.timesteps:
                raise Exception("Reading the cached PV values seems to have failed. Expected "
                                + str(self.my_simulation_parameters.timesteps) + " values, but got " + str(len(self.output )))
        else:
            self.get_coordinates(location = self.pvconfig.location, year =  self.pvconfig.time)
            # Factor to guarantee peak power based on module with 250 Wh
            self.ac_power_factor = math.ceil( ( self.pvconfig.power * 1e3 ) / 250 )

            #when predictive control is activated, the PV simulation is run beforhand to make forecasting easier
            if self.my_simulation_parameters.system_config.predictive and my_simulation_repository is not None:
                #get yearly weather data from dictionary
                dni_extra = my_simulation_repository.get_entry( Weather.Weather_DirectNormalIrradianceExtra_yearly_forecast )
                DNI = my_simulation_repository.get_entry( Weather.Weather_DirectNormalIrradiance_yearly_forecast )
                DHI = my_simulation_repository.get_entry( Weather.Weather_DiffuseHorizontalIrradiance_yearly_forecast )
                GHI = my_simulation_repository.get_entry( Weather.Weather_GlobalHorizontalIrradiance_yearly_forecast )
                azimuth = my_simulation_repository.get_entry( Weather.Weather_Azimuth_yearly_forecast )
                apparent_zenith = my_simulation_repository.get_entry( Weather.Weather_ApparentZenith_yearly_forecast )
                temperature = my_simulation_repository.get_entry( Weather.Weather_TemperatureOutside_yearly_forecast )
                wind_speed = my_simulation_repository.get_entry( Weather.Weather_WindSpeed_yearly_forecast )

                x= [ ]
                for i in range( len( dni_extra ) ):
                    x.append( simPhotovoltaicFast( dni_extra[ i ], DNI[ i ], DHI[ i ], GHI[ i ], azimuth[ i ], apparent_zenith[ i ], temperature[ i ], wind_speed[ i ] ) )

                self.output = x

                database = pd.DataFrame( self.output, columns = [ "output" ] )

                database.to_csv( self.cache_filepath, sep=",", decimal=".", index=False )

            else:
                 self.data = [0] * self.my_simulation_parameters.timesteps
                 self.data_length = self.my_simulation_parameters.timesteps

        if self.my_simulation_parameters.system_config.predictive and my_simulation_repository is not None:
            my_simulation_repository.delete_entry( Weather.Weather_DirectNormalIrradianceExtra_yearly_forecast )
            my_simulation_repository.delete_entry( Weather.Weather_DirectNormalIrradiance_yearly_forecast )
            my_simulation_repository.delete_entry( Weather.Weather_DiffuseHorizontalIrradiance_yearly_forecast )
            my_simulation_repository.delete_entry( Weather.Weather_GlobalHorizontalIrradiance_yearly_forecast )
            my_simulation_repository.delete_entry( Weather.Weather_Azimuth_yearly_forecast )
            my_simulation_repository.delete_entry( Weather.Weather_ApparentZenith_yearly_forecast )
            my_simulation_repository.delete_entry( Weather.Weather_TemperatureOutside_yearly_forecast )
            my_simulation_repository.delete_entry( Weather.Weather_WindSpeed_yearly_forecast )

        self.modules = pd.read_csv(
            os.path.join(utils.HISIMPATH["photovoltaic"]["modules"]),
            index_col=0,
        )

        self.inverters = pd.read_csv(
            os.path.join(utils.HISIMPATH["photovoltaic"]["inverters"]),
            index_col=0,
        )

        self.temp_model = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS["sapm"]["open_rack_glass_glass"]

        # load the sandia data
        if load_module_data:
            # load module data online
            modules = pvlib.pvsystem.retrieve_sam(name="SandiaMod")
            self.module = modules[self.pvconfig.module_name]
            # get inverter data
            inverters = pvlib.pvsystem.retrieve_sam("cecinverter")
            self.inverter = inverters[self.pvconfig.inverter_name]
        else:
            # load module and inverter data from csv
            module = self.modules[self.pvconfig.module_name]
            self.module = pd.to_numeric(module, errors="coerce")

            inverter = self.inverters[self.pvconfig.inverter_name]
            self.inverter = pd.to_numeric(inverter, errors="coerce")
        #self.power = self.power
        #self.module_name =  module_name
        #self.inverter_name = inverter_name
        #self.integrateInverter = integrateInverter
        #self.simPhotovoltaicSimpleJit = simPhotovoltaicSimple

    def plot(self):
        # Plots ac_power. One day is represented by 1440 steps.
        #self.ac_power.iloc[0:7200].plot()
        plt.plot(self.data)
        plt.ylabel("Power [W]")
        plt.xlabel("Time")
        plt.show()

    def interpolate(self,pd_database,year):
        firstday = pd.Series([0.0], index=[
            pd.to_datetime(datetime.datetime(year-1, 12, 31, 23, 0), utc=True).tz_convert("Europe/Berlin")])
        lastday = pd.Series(pd_database[-1], index=[
            pd.to_datetime(datetime.datetime(year, 12, 31, 22, 59), utc=True).tz_convert("Europe/Berlin")])
        #pd_database = pd_database.append(firstday)
        pd_database = pd_database.append(lastday)
        pd_database = pd_database.sort_index()
        return pd_database.resample('1T').asfreq().interpolate(method='linear').tolist()

    def simPhotovoltaic2(
        self,
        dni_extra=None,
        DNI=None,
        DHI=None,
        GHI=None,
        azimuth=None,
        apparent_zenith=None,
        temperature=None,
        wind_speed=None,
        surface_tilt=30,
        surface_azimuth=180,
        albedo=0.2):
        """
        Simulates a defined PV array with the Sandia PV Array Performance Model.
        The implementation is done in accordance with following tutorial:
        https://github.com/pvlib/pvlib-python/blob/master/docs/tutorials/tmy_to_power.ipynb

        Based on the tsib project @[tsib-kotzur] (Check header)

        Parameters
        ----------
        tmy_data: pandas.DataFrame(), required
            Weatherfile in the format of a tmy file.
        surface_tilt: int or float, optional (default:30)
            Tilt angle of of the array in degree.
        surface_azimuth: int or float, optional (default:180)
            Azimuth angle of of the array in degree. 180 degree means south,
            90 degree east and 270 west.
        albedo: float, optional (default: 0.2)
            Reflection coefficient of the surrounding area.
        losses: float, optional (default: 0.1)
            Losses due to soiling, mismatch, diode connections, dc wiring etc.
        load_module_data: Boolean, optional (default: False)
            If True the module data base is loaded from the Sandia Website.
            Otherwise it is loaded from this relative path
                '\\profiles\\PV-Modules\\sandia_modules.csv'.
        module_name: str, optional (default:'Hanwha_HSL60P6_PA_4_250T__2013_')
            Module name. The string must be existens in Sandia Module database.
        integrateInverter: bool, optional (default: True)
            If an inverter shall be added to the simulation, providing the photovoltaic output after the inverter.
        inverter_name: str, optional (default: 'ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_')
            Type of inverter.

        Returns
        --------
        """
        # automatic pd time series in future pvlib version
        # calculate airmass
        airmass = pvlib.atmosphere.get_relative_airmass(apparent_zenith)
        # use perez model to calculate the plane of array diffuse sky radiation
        poa_sky_diffuse = pvlib.irradiance.perez(
            surface_tilt,
            surface_azimuth,
            DHI,
            np.float64(DNI),
            dni_extra,
            apparent_zenith,
            azimuth,
            airmass,
        )
        # calculate ground diffuse with specified albedo
        poa_ground_diffuse = pvlib.irradiance.get_ground_diffuse(
            surface_tilt, GHI, albedo=albedo
        )
        # calculate angle of incidence
        aoi = pvlib.irradiance.aoi(surface_tilt, surface_azimuth, apparent_zenith, azimuth)
        # calculate plane of array irradiance
        poa_irrad = pvlib.irradiance.poa_components(aoi, np.float64(DNI), poa_sky_diffuse, poa_ground_diffuse)
        # calculate pv cell and module temperature
        #temp_model = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS["sapm"]["open_rack_glass_glass"]
        pvtemps = pvlib.temperature.sapm_cell(poa_irrad["poa_global"], temperature, wind_speed, **self.temp_model)

        # calculate effective irradiance on pv module
        sapm_irr = pvlib.pvsystem.sapm_effective_irradiance(
            module=self.module,
            poa_direct=poa_irrad["poa_direct"],
            poa_diffuse=poa_irrad["poa_diffuse"],
            airmass_absolute=airmass,
            aoi=aoi,
        )
        # calculate pv performance
        sapm_out = pvlib.pvsystem.sapm(
            sapm_irr,
            module=self.module,
            temp_cell=pvtemps,
        )
        # calculate peak load of single module [W]
        peak_load = self.module.loc["Impo"] * self.module.loc["Vmpo"]
        ac_power = pd.DataFrame()
        if self.pvconfig.integrate_inverter:
            # calculate load after inverter
            iv_load = pvlib.inverter.sandia(inverter=self.inverter, v_dc=sapm_out["v_mp"], p_dc=sapm_out["p_mp"])
            ac_power = iv_load / peak_load
        else:
            # load in [kW/kWp]
            ac_power = sapm_out["p_mp"] / peak_load

        if math.isnan(ac_power):
            ac_power = 0.0

        #ac_power = ac_power * self.time_correction_factor
        #ac_power = ac_power

        #data = [DHI,
        #        DNI,
        #        GHI,
        #        dni_extra,
        #        aoi,
        #        apparent_zenith,
        #        azimuth,
        #        airmass,
        #        wind_speed]
        #if timestep % 60 == 0 and timestep < 1442:
        #    log.information(data)
        #    log.information("Timestep:{} , AcPower: {}".format(timestep, ac_power))

        return ac_power

def readTRY(location="Aachen", year=2010):
    """
    Reads a test reference year file and gets the GHI, DHI and DNI from it.

    Based on the tsib project @[tsib-kotzur] (Check header)

    Parameters
    -------
    try_num: int (default: 4)
        The region number of the test reference year.
    year: int (default: 2010)
        The year. Only data for 2010 and 2030 available
    """
    # get the correct file path
    filepath = os.path.join(utils.HISIMPATH["weather"][location])

    # get the geoposition
    with open(filepath + ".dat", encoding="utf-8") as fp:
        lines = fp.readlines()
        location_name = lines[0].split(maxsplit=2)[2].replace('\n', '')
        lat = float(lines[1][20:37])
        lon = float(lines[2][15:30])
    location = {"name": location_name, "latitude": lat, "longitude": lon}

    # check if time series data already exists as .csv with DNI
    if os.path.isfile(filepath + ".csv"):
        data = pd.read_csv(filepath + ".csv", index_col=0, parse_dates=True,sep=";",decimal=",")
        data.index = pd.to_datetime(data.index, utc=True).tz_convert("Europe/Berlin")
    # else read from .dat and calculate DNI etc.
    else:
        # get data
        data = pd.read_csv(
            filepath + ".dat", sep=r"\s+", skiprows=([i for i in range(0, 31)])
        )
        data.index = pd.date_range(
            "{}-01-01 00:00:00".format(year), periods=8760, freq="H", tz="Europe/Berlin"
        )
        data["GHI"] = data["D"] + data["B"]
        data = data.rename(columns={"D": "DHI", "t": "T", "WG": "WS"})

        # calculate direct normal
        data["DNI"] = calculateDNI(data["B"], lon, lat)
        # data["DNI"] = data["B"]

        # save as .csv
        #data.to_csv(filepath + ".csv",sep=";",decimal=",")
    return data, location

def calculateDNI(directHI, lon, lat, zenith_tol=87.0):
    """
    Calculates the direct NORMAL irradiance from the direct horizontal irradiance with the help of the PV lib.

    Based on the tsib project @[tsib-kotzur] (Check header)

    Parameters
    ----------
    directHI: pd.Series with time index
        Direct horizontal irradiance
    lon: float
        Longitude of the location
    lat: float
        Latitude of the location
    zenith_tol: float, optional
        Avoid cosines of values above a certain zenith angle of in order to avoid division by zero.

    Returns
    -------
    DNI: pd.Series
    """
    solarPos = pvlib.solarposition.get_solarposition(directHI.index, lat, lon)
    solarPos["apparent_zenith"][solarPos.apparent_zenith > zenith_tol] = zenith_tol
    DNI = directHI.div(solarPos["apparent_zenith"].apply(math.radians).apply(math.cos))
    DNI = DNI.fillna(0)
    if DNI.isnull().values.any():
        raise ValueError("Something went wrong...")
    return DNI


text:
#!/usr/bin/env python3
"""
Input-output routines

"""
import numpy as np
import sys

def save_xyzfile(fname, names, mat):
    """Save a coordinate matrix in a xyz format"""
    M = len(mat)
    s = "%i\nbla\n" % M
    for i in range(M):
        s += "%i\t%.10f\t%.10f\t%.10f\n" % \
                (names[i], mat[i, 0], mat[i, 1], mat[i, 2])
    open(fname, "w").write(s)

def read_xyzfile(fname):
    try:
        A = open(fname, "r").readlines()[2:]
    except FileNotFoundError:
        sys.exit("File %s not found." % fname)
    A = [line.split() for line in A]
    A = np.array(A).astype(float)
    if A.shape[1] != 4:
        sys.exit("Incorrect number of columns in %s." % fname)
    nm, xyz = A[:, 0], A[:, 1:]
    return nm, xyz

def parse_box(arg):
    L = [float(i) for i in arg.split()]
    if len(L) == 1:
        box = np.eye(3) * L[0]
    elif len(L) == 3:
        box = np.diag(L)
    else:
        sys.exit("L should have length 1 or 3.")
    return box

ext: py

max_stars_repo_name: petervanya/DPDsim

max_issues_repo_name: petervanya/DPDsim

max_forks_repo_name: petervanya/DPDsim

reason: import numpy


text:
"""
===============================================
vidgear library source-code is deployed under the Apache 2.0 License:

Copyright (c) 2019 Abhishek Thakur(@abhiTronix) <abhi.una12@gmail.com>

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
===============================================
"""

from vidgear.gears import NetGear
from vidgear.gears import VideoGear

import pytest
import cv2
import random
import tempfile
import os
import numpy as np
import traceback
from zmq.error import ZMQError
import logging as log

logger = log.getLogger('Test_netgear')

def return_testvideo_path():
	"""
	returns Test Video path
	"""
	path = '{}/Downloads/Test_videos/BigBuckBunny_4sec.mp4'.format(tempfile.gettempdir())
	return os.path.abspath(path)

def test_playback():
	"""
	Tests NetGear Bare-minimum network playback capabilities
	"""
	try:
		#open stream
		stream = VideoGear(source=return_testvideo_path()).start()
		#open server and client with default params
		client = NetGear(receive_mode = True)
		server = NetGear()
		#playback
		while True:
			frame_server = stream.read()
			if frame_server is None: break
			server.send(frame_server) #send
			frame_client = client.recv() #recv
		#clean resources
		stream.stop()
		server.close()
		client.close()
	except Exception as e:
		if isinstance(e, (ZMQError, ValueError)):
			logger.debug(traceback.print_tb(e.__traceback__))
		else:
			pytest.fail(str(e))

@pytest.mark.parametrize('pattern', [1, 2, 3]) #0:(zmq.PAIR,zmq.PAIR), 1:(zmq.REQ,zmq.REP), 2:(zmq.PUB,zmq.SUB) (#3 is incorrect value)
def test_patterns(pattern):
	"""
	Testing NetGear different messaging patterns
	"""
	#open stream
	try:
		stream = VideoGear(source=return_testvideo_path()).start()

		#define parameters
		options = {'flag' : 0, 'copy' : True, 'track' : False}
		client = NetGear(pattern = pattern, receive_mode = True, logging = True, **options)
		server = NetGear(pattern = pattern, logging = True, **options)
		#initialize
		frame_server = None
		#select random input frame from stream
		i = 0
		while (i < random.randint(10, 100)):
			frame_server = stream.read()
			i+=1
		#check if input frame is valid
		assert not(frame_server is None)
		#send frame over network
		server.send(frame_server)
		frame_client = client.recv()
		#clean resources
		stream.stop()
		server.close()
		client.close()
		#check if recieved frame exactly matches input frame
		assert np.array_equal(frame_server, frame_client)
	except Exception as e:
		if isinstance(e, (ZMQError, ValueError)):
			logger.error(traceback.print_tb(e.__traceback__))
		else:
			pytest.fail(str(e))

def test_compression():
	"""
	Testing NetGear's real-time frame compression capabilities
	"""
	try:
		#open streams
		stream = VideoGear(source=return_testvideo_path()).start()
		#define client parameters
		options = {'compression_param':cv2.IMREAD_COLOR} #read color image
		client = NetGear(pattern = 1, receive_mode = True, logging = True, **options)
		#define server parameters
		options = {'compression_format': '.jpg', 'compression_param':[cv2.IMWRITE_JPEG_OPTIMIZE, 20]} #JPEG compression
		server = NetGear(pattern = 1, logging = True, **options)
		#send over network
		while True:
			frame_server = stream.read()
			if frame_server is None: break
			server.send(frame_server)
			frame_client = client.recv()
		#clean resources
		stream.stop()
		server.close()
		client.close()
	except Exception as e:
		if isinstance(e, (ZMQError, ValueError)):
			logger.error(traceback.print_tb(e.__traceback__))
		else:
			pytest.fail(str(e))

test_data_class = [
	(0,1, tempfile.gettempdir(), True),
	(1,1, tempfile.gettempdir(), False)]
@pytest.mark.parametrize('pattern, security_mech, custom_cert_location, overwrite_cert', test_data_class)
def test_secure_mode(pattern, security_mech, custom_cert_location, overwrite_cert):
	"""
	Testing NetGear's Secure Mode
	"""
	try:
		#open stream
		stream = VideoGear(source=return_testvideo_path()).start()

		#define security mechanism
		options = {'secure_mode': security_mech, 'custom_cert_location': custom_cert_location, 'overwrite_cert': overwrite_cert}

		#define params
		server = NetGear(pattern = pattern, logging = True, **options)
		client = NetGear(pattern = pattern, receive_mode = True, logging = True, **options)
		#initialize
		frame_server = None
		#select random input frame from stream
		i = 0
		while (i < random.randint(10, 100)):
			frame_server = stream.read()
			i+=1
		#check input frame is valid
		assert not(frame_server is None)
		#send and recv input frame
		server.send(frame_server)
		frame_client = client.recv()
		#clean resources
		stream.stop()
		server.close()
		client.close()
		#check if recieved frame exactly matches input frame
		assert np.array_equal(frame_server, frame_client)
	except Exception as e:
		if isinstance(e, (ZMQError, ValueError)):
			logger.error(traceback.print_tb(e.__traceback__))
		else:
			pytest.fail(str(e))

@pytest.mark.parametrize('target_data', [[1,'string',['list']], {1:'apple', 2: 'cat'}])
def test_bidirectional_mode(target_data):
	"""
	Testing NetGear's Bidirectional Mode with different datatypes
	"""
	try:
		logger.debug('Given Input Data: {}'.format(target_data))

		#open strem
		stream = VideoGear(source=return_testvideo_path()).start()
		#activate bidirectional_mode
		options = {'bidirectional_mode': True}
		#define params
		client = NetGear(receive_mode = True, logging = True, **options)
		server = NetGear(logging = True, **options)
		#get frame from stream
		frame_server = stream.read()
		assert not(frame_server is None)

		#sent frame and data from server to client
		server.send(frame_server, message = target_data)
		#client recieves the data and frame and send its data
		server_data, frame_client = client.recv(return_data = target_data)
		#server recieves the data and cycle continues
		client_data = server.send(frame_server, message = target_data)

		#clean resources
		stream.stop()
		server.close()
		client.close()

		#logger.debug data recieved at client-end and server-end
		logger.debug('Data recieved at Server-end: {}'.format(server_data))
		logger.debug('Data recieved at Client-end: {}'.format(client_data))

		#check if recieved frame exactly matches input frame
		assert np.array_equal(frame_server, frame_client)
		#check if client-end data exactly matches server-end data
		assert client_data == server_data
	except Exception as e:
		if isinstance(e, (ZMQError, ValueError)):
			logger.error(traceback.print_tb(e.__traceback__))
		else:
			pytest.fail(str(e))

def test_multiserver_mode():
	"""
	Testing NetGear's Multi-Server Mode with three unique servers
	"""
	try:
		#open network stream
		stream = VideoGear(source=return_testvideo_path()).start()

		#define and activate Multi-Server Mode
		options = {'multiserver_mode': True}

		#define a single client
		client = NetGear(port = ['5556', '5557', '5558'], pattern = 1, receive_mode = True, logging = True, **options)
		#define client-end dict to save frames inaccordance with unique port
		client_frame_dict = {}

		#define three unique server
		server_1 = NetGear(pattern = 1, port = '5556', logging = True, **options) #at port `5556`
		server_2 = NetGear(pattern = 1, port = '5557', logging = True, **options) #at port `5557`
		server_3 = NetGear(pattern = 1, port = '5558', logging = True, **options) #at port `5558`

		#generate a random input frame
		frame_server = None
		i = 0
		while (i < random.randint(10, 100)):
			frame_server = stream.read()
			i+=1
		#check if input frame is valid
		assert not(frame_server is None)

		#send frame from Server-1 to client and save it in dict
		server_1.send(frame_server)
		unique_address, frame = client.recv()
		client_frame_dict[unique_address] = frame
		#send frame from Server-2 to client and save it in dict
		server_2.send(frame_server)
		unique_address, frame = client.recv()
		client_frame_dict[unique_address] = frame
		#send frame from Server-3 to client and save it in dict
		server_3.send(frame_server)
		unique_address, frame = client.recv()
		client_frame_dict[unique_address] = frame

		#clean all resources
		stream.stop()
		server_1.close()
		server_2.close()
		server_3.close()
		client.close()

		#check if recieved frames from each unique server exactly matches input frame
		for key in client_frame_dict.keys():
			assert np.array_equal(frame_server, client_frame_dict[key])

	except Exception as e:
		if isinstance(e, (ZMQError, ValueError)):
			logger.error(traceback.print_tb(e.__traceback__))
		else:
			pytest.fail(str(e))

max_stars_repo_name: cclauss/vidgear

max_issues_repo_name: cclauss/vidgear

max_forks_repo_name: cclauss/vidgear


text:
#!/usr/bin/env python

# Copyright (c) 2013, Carnegie Mellon University
# All rights reserved.
# Authors: Michael Koval <mkoval@cs.cmu.edu>
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
# - Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
# - Neither the name of Carnegie Mellon University nor the names of its
#   contributors may be used to endorse or promote products derived from this
#   software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

import numpy
import openravepy
from ..collision import (
    BakedRobotCollisionCheckerFactory,
    DefaultRobotCollisionCheckerFactory,
    SimpleRobotCollisionCheckerFactory,
)
from tsr import TSR, TSRChain
from ..util import SetTrajectoryTags
from .adapters import PlanToEndEffectorOffsetTSRAdapter
from .base import (
    Planner,
    LockedPlanningMethod,
    PlanningError,
    Tags,
    UnsupportedPlanningError,
    save_dof_limits
)
import contextlib

class CBiRRTPlanner(Planner):
    def __init__(self, robot_checker_factory=None, timelimit=5.):
        super(CBiRRTPlanner, self).__init__()

        self.timelimit = timelimit

        if robot_checker_factory is None:
            robot_checker_factory = DefaultRobotCollisionCheckerFactory

        self.robot_checker_factory = robot_checker_factory
        if isinstance(robot_checker_factory, SimpleRobotCollisionCheckerFactory):
            self._is_baked = False
        elif isinstance(robot_checker_factory, BakedRobotCollisionCheckerFactory):
            self._is_baked = True
        else:
            raise NotImplementedError(
                'CBiRRT only supports Simple and BakedRobotCollisionChecker.')

    def __str__(self):
        return 'CBiRRT'

    @LockedPlanningMethod
    def PlanToConfigurations(self, robot, goals, **kw_args):
        """
        Plan to multiple goal configurations with CBiRRT. This adds each goal
        in goals as a root node in the tree and returns a path that reaches one
        of the goals.
        @param robot
        @param goals list of goal configurations
        @return traj output path
        """
        kw_args.setdefault('smoothingitrs', 0)
        return self.Plan(robot, jointgoals=goals, **kw_args)

    @LockedPlanningMethod
    def PlanToConfiguration(self, robot, goal, **kw_args):
        """
        Plan to a single goal configuration with CBiRRT.
        @param robot
        @param goal goal configuration
        @return traj output path
        """
        kw_args.setdefault('smoothingitrs', 0)
        return self.Plan(robot, jointgoals=[goal], **kw_args)

    @LockedPlanningMethod
    def PlanToEndEffectorPose(self, robot, goal_pose, **kw_args):
        """
        Plan to a desired end-effector pose.
        @param robot
        @param goal desired end-effector pose
        @param psample probability of sampling a goal
        @return traj output path
        """
        manipulator_index = robot.GetActiveManipulatorIndex()
        goal_tsr = TSR(T0_w=goal_pose, manipindex=manipulator_index)
        tsr_chain = TSRChain(sample_goal=True, TSR=goal_tsr)

        kw_args.setdefault('psample', 0.1)
        kw_args.setdefault('smoothingitrs', 0)

        return self.Plan(robot, tsr_chains=[tsr_chain], **kw_args)

    @LockedPlanningMethod
    def PlanToEndEffectorOffset(self, robot, direction, distance, **kw_args):
        """
        Plan to a desired end-effector offset.

        @param robot
        @param direction unit vector in the direction of motion
        @param distance minimum distance in meters
        @param smoothingitrs number of smoothing iterations to run
        @return traj output path
        """
        chains = PlanToEndEffectorOffsetTSRAdapter.CreateTSRChains(
                    robot, direction, distance)
        return self.PlanToTSR(robot, tsr_chains=chains, **kw_args)

    @LockedPlanningMethod
    def PlanToTSR(self, robot, tsr_chains, smoothingitrs=100, **kw_args):
        """
        Plan to a goal specified as a list of TSR chains. CBiRRT supports an
        arbitrary list of start, goal, and/or constraint TSR chains. The path
        will be smoothed internally by CBiRRT if one or more constraint TSR
        chains are specified.
        @param robot
        @param tsr_chains list of TSR chains
        @param smoothingitrs number of smoothing iterations to run
        @return traj output path
        """
        is_constrained = False

        for chain in tsr_chains:
            if chain.sample_start or chain.sample_goal:
                kw_args.setdefault('psample', 0.1)

            if chain.constrain:
                is_constrained = True

        # Only smooth constrained trajectories.
        if not is_constrained:
            smoothingitrs = 0

        return self.Plan(
            robot,
            smoothingitrs=smoothingitrs,
            tsr_chains=tsr_chains,
            **kw_args
        )

    def Plan(self, robot, smoothingitrs=None, timelimit=None, allowlimadj=0,
             jointstarts=None, jointgoals=None, psample=None, tsr_chains=None,
             extra_args=None, **kw_args):
        """
        @param allowlimadj If True, adjust the joint limits to include
            the robot's start configuration
        """
        from openravepy import CollisionOptionsStateSaver, Robot, KinBody

        if timelimit is None:
            timelimit = self.timelimit

        if timelimit <= 0.:
            raise ValueError('Invalid value for "timelimit". Limit must be'
                             ' non-negative; got {:f}.'.format(timelimit))

        env = robot.GetEnv()
        problem = openravepy.RaveCreateProblem(env, 'CBiRRT')

        if problem is None:
            raise UnsupportedPlanningError('Unable to create CBiRRT module.')

        is_endpoint_deterministic = True
        is_constrained = False

        # TODO We may need this work-around because CBiRRT doesn't like it
        # when an IK solver other than GeneralIK is loaded (e.g. nlopt_ik).
        # self.ClearIkSolver(robot.GetActiveManipulator())

        env.LoadProblem(problem, robot.GetName())

        args =  ['RunCBiRRT']
        args += ['timelimit', str(timelimit)]

        # By default, CBiRRT interprets the DOF resolutions as an
        # L-infinity norm; this flag turns on the L-2 norm instead.
        args += ['bdofresl2norm', '1']
        args += ['steplength', '0.05999']

        if self._is_baked:
            args += ['bbakedcheckers', '1']

        if extra_args is not None:
            args += extra_args

        if smoothingitrs is not None:
            if smoothingitrs < 0:
                raise ValueError('Invalid number of smoothing iterations. '
                                 'Value must be non-negative; got  {:d}.'
                                 .format(smoothingitrs))

            args += ['smoothingitrs', str(smoothingitrs)]

        if allowlimadj is not None:
            args += ['allowlimadj', str(int(allowlimadj))]

        if psample is not None:
            if not (0 <= psample <= 1):
                raise ValueError('Invalid value for "psample". Value must be '
                                 'in the range [0, 1]; got {:f}.'
                                 .format(psample))

            args += ['psample', str(psample)]

        if jointstarts is not None:
            for start_config in jointstarts:
                if len(start_config) != robot.GetActiveDOF():
                    raise ValueError(
                        'Incorrect number of DOFs in start configuration;'
                        ' expected {:d}, got {:d}'.format(
                            robot.GetActiveDOF(), len(start_config)
                        )
                    )

                args += (['jointstarts'] +
                         self.serialize_dof_values(start_config))

        if jointgoals is not None:
            for goal_config in jointgoals:
                if len(goal_config) != robot.GetActiveDOF():
                    raise ValueError(
                        'Incorrect number of DOFs in goal configuration;'
                        ' expected {:d}, got {:d}'.format(
                            robot.GetActiveDOF(), len(goal_config)
                        )
                    )

                args += ['jointgoals'] + self.serialize_dof_values(goal_config)

            if len(jointgoals) > 1:
                is_endpoint_deterministic = False

        if tsr_chains is not None:
            for tsr_chain in tsr_chains:
                args += ['TSRChain', SerializeTSRChain(tsr_chain)]

                if tsr_chain.sample_goal:
                    is_endpoint_deterministic = False
                if tsr_chain.constrain:
                    is_constrained = True

        # FIXME: Why can't we write to anything other than cmovetraj.txt or
        # /tmp/cmovetraj.txt with CBiRRT?
        traj_path = 'cmovetraj.txt'
        args += ['filename', traj_path]
        args_str = ' '.join(args)

        # Bypass the context manager since CBiRRT does its own baking in C++.
        collision_checker = self.robot_checker_factory(robot)
        options = collision_checker.collision_options

        if tsr_chains is not None:
            mimicbodies = [env.GetKinBody(chain.mimicbodyname)
                           for chain in tsr_chains if chain.mimicbodyname is not 'NULL']
            mimicbody_savers = [
                mimicbody.CreateKinBodyStateSaver(KinBody.SaveParameters.LinkTransformation)
                for mimicbody in mimicbodies]
        else:
            mimicbody_savers = []

        with CollisionOptionsStateSaver(env.GetCollisionChecker(), options), \
            robot.CreateRobotStateSaver(Robot.SaveParameters.ActiveDOF |
                                        Robot.SaveParameters.LinkTransformation), \
            contextlib.nested(*mimicbody_savers), save_dof_limits(robot):
            response = problem.SendCommand(args_str, True)

        if not response.strip().startswith('1'):
            raise PlanningError('Unknown error: ' + response,
                deterministic=False)

        # Construct the output trajectory.
        with open(traj_path, 'rb') as traj_file:
            traj_xml = traj_file.read()
            traj = openravepy.RaveCreateTrajectory(env, '')
            traj.deserialize(traj_xml)

        # Tag the trajectory as non-determistic since CBiRRT is a randomized
        # planner. Additionally tag the goal as non-deterministic if CBiRRT
        # chose from a set of more than one goal configuration.
        SetTrajectoryTags(traj, {
            Tags.CONSTRAINED: is_constrained,
            Tags.DETERMINISTIC_TRAJECTORY: False,
            Tags.DETERMINISTIC_ENDPOINT: is_endpoint_deterministic,
        }, append=True)

        # Strip extraneous groups from the output trajectory.
        # TODO: Where are these groups coming from!?
        cspec = robot.GetActiveConfigurationSpecification('linear')
        openravepy.planningutils.ConvertTrajectorySpecification(traj, cspec)

        return traj

    def ClearIkSolver(self, manip):
        manip.SetIkSolver(None)
        manip.SetLocalToolTransform(manip.GetLocalToolTransform())
        manip.SetIkSolver(None)
        if manip.GetIkSolver() is not None:
            raise ValueError('Unable to clear IkSolver')

    @staticmethod
    def serialize_dof_values(dof_values):
        return [str(len(dof_values)),
                ' '.join([str(x) for x in dof_values])]

def SerializeTransform12Col(tm, format='%.5f'):
    return ' '.join([(format % (i,)) for i in tm[0:3, :].T.reshape(12)])

def SerializeArray(a, format='%.5f'):
    return ' '.join([(format % (i,)) for i in a.reshape(-1)])

def SerializeTSR(self):
    """
    Function for Serializing TSRs for CBIRRT.

    SerializeTSR(manipindex,bodyandlink,T0_w,Tw_e,Bw)

    Input:
    manipindex (int): the 0-indexed index of the robot's manipulator
    bodyandlink (str): body and link which is used as the 0 frame. Format
                       'body_name link_name'. For world frame, specify 'NULL'
    T0_w (double 4x4): transform matrix of the TSR's reference frame relative
                       to the 0 frame
    Tw_e (double 4x4): transform matrix of the TSR's offset frame relative to
                       the w frame
    Bw (double 1x12): bounds in x y z roll pitch yaw.
                      Format: [x_min, x_max, y_min, y_max ...]

    Output:
    outstring (str): string to use for SerializeTSRChain function
    """
    return '%d %s %s %s %s' % (self.manipindex, self.bodyandlink,
                               SerializeTransform12Col(self.T0_w),
                               SerializeTransform12Col(self.Tw_e),
                               SerializeArray(self.Bw))

def SerializeTSRChain(self):
    """
    Function for Serializing TSR Chains for CBIRRT.

    _SerializeTSRChain(bSampleFromChain, bConstrainToChain,
                       numTSRs, allTSRstring,
                       mimicbodyname, mimicbodyjoints)

    Input:
    bSampleStartFromChain (0/1): 1: Use this chain for sampling start configs
                                 0: Ignore for sampling starts
    bSampleGoalFromChain (0/1): 1: Use this chain for sampling goal configs
                                0: Ignore for sampling goals
    bConstrainToChain (0/1): 1: Use this chain for constraining configs
                             0: Ignore for constraining
    numTSRs (int): Number of TSRs in this chain (must be > 0)
    allTSRstring (str): string of concatenated TSRs from SerializeTSR.
                        Should be like [TSRstring 1 ' ' TSRstring2 ...]
    mimicbodyname (str): name of associated mimicbody for this chain
                         (NULL if none associated)
    mimicbodyjoints (int [1xn]): 0-indexed indices of mimicbody's joints that
                                 are mimiced (INCREASING AND CONSECUTIVE)

    Output:
    outstring (str): string to include in call to cbirrt planner
    """
    allTSRstring = ' '.join([SerializeTSR(tsr) for tsr in self.TSRs])
    numTSRs = len(self.TSRs)
    outstring = '%d %d %d' % (int(self.sample_start),
                              int(self.sample_goal),
                              int(self.constrain))
    outstring += ' %d %s' % (numTSRs, allTSRstring)
    outstring += ' ' + self.mimicbodyname
    if len(self.mimicbodyjoints) > 0:
        outstring += ' %d %s' % (len(self.mimicbodyjoints),
                                 SerializeArray(self.mimicbodyjoints))
    return outstring

max_stars_repo_name: MisoRobotics/prpy

max_issues_repo_name: MisoRobotics/prpy

max_forks_repo_name: MisoRobotics/prpy


text:
# -*- coding: utf-8 -*-
#
#  closures.py
#  pyobserver
#
#  Created by Alexander Rudy on 2014-07-19.
#  Copyright 2014 Alexander Rudy. All rights reserved.
#
"""
Handle LCH Closures
"""
from __future__ import (absolute_import, unicode_literals, division, print_function)

import re
import six
from six.moves import zip
import itertools
import warnings
import collections
import itertools
import numpy as np
import glob
import os, os.path
import datetime
import weakref

import astropy.units as u
from astropy.time import Time
from astropy.coordinates import FK5, AltAz, FK4
from astropy.utils.console import human_time
from astropy.utils.data import get_readable_fileobj

from .targets import Target
from .starlist import read_skip_comments

__all__ = ['Region', 'Closure', 'Opening', 'Regions']

class Region(Target):
    """Region for LCH Closures."""

    RADIUS = 2 * u.arcmin
    """Radius around this region covered by the closure."""

    def __init__(self, *args, **kwargs):
        super(Region, self).__init__(*args, **kwargs)
        self._openings = list()

    def contains(self, position):
        """Test whether a position is within this LCH region."""
        try:
            contained = (self.position.separation(location) > self.RADIUS)
        except (AttributeError, ValueError):
            return False
        return contained

    @property
    def openings(self):
        """A list of openings, in order."""
        return self._openings

    def add(self, opening, _no_sort=False):
        """Add an opening"""
        if opening not in self._openings:
            self._openings.append(opening)
        if not _no_sort:
            self._openings.sort(key=lambda o : o.start)

    def open(self, time):
        """Check if this window is open during a specific time."""
        for opening in self._openings:
            if (time > opening.start) and (time < opening.end):
                return True
        return False

    def closed(self, time):
        """Check if this window is closed during a specific time."""
        return (not self.open(time))

    @property
    def closures(self):
        """Get a closure list."""
        for open_start, open_end in self._iter_start_end():
            yield Closure(self, open_start.end, open_end.start)

    def _iter_start_end(self):
        """Iterate over pairs of openings, useful for iterating through closures, or through closures and openings."""
        if len(self.openings) == 0:
            return iter([])
        starts = itertools.islice(self.openings, len(self.openings)-1) # Take everything but the last element
        ends = itertools.islice(self.openings, 1, None) # Take everything but the first element
        return zip(starts, ends)

    @property
    def events(self):
        """Interleaved closures and openings."""
        for open_start, open_end in self._iter_start_end():
            yield open_start
            yield Closure(self, open_start.end, open_end.start)
        yield open_end

class WindowTimeProperty(object):
    """Descriptor for time properties."""
    def __init__(self, name):
        super(WindowTimeProperty, self).__init__()
        self.name = name

    @property
    def attr(self):
        """Attribute name."""
        return "_{}_{}".format(self.__class__.__name__, self.name)

    def __set__(self, obj, value):
        """Set the property."""
        setattr(obj, self.attr, Time(value))

    def __get__(self, obj, objtype=None):
        """Get the property"""
        if obj is None:
            return self
        if not hasattr(obj, self.attr):
            raise AttributeError("{:s}.{:s} has not been initialized".format(obj.__class__.__name__, self.name))
        return getattr(obj, self.attr)

class Window(object):
    """An LCH Window"""
    def __init__(self, region, start, end):
        super(Window, self).__init__()
        self._region = weakref.ref(region)
        self.start = start
        self.end = end

    start = WindowTimeProperty("start")
    end = WindowTimeProperty("end")

    def __repr__(self):
        """Representer for this LCH window object."""
        try:
            return "<{0:s} from {1.datetime:%H:%M:%S} to {2.datetime:%H:%M:%S} ({3:.0f})>".format(self.__class__.__name__, self.start, self.end, self.duration)
        except AttributeError:
            return super(Window, self).__repr__()

    @property
    def region(self):
        """Region property"""
        return self._region()

    @property
    def duration(self):
        """The duration of this closure."""
        return (self.end - self.start).to(u.s)

    def time_to(self, time = None, string = False):
        """Compute the time until this closure starts, from a given time.

        :param time: The reference time for this closure.
        :param string: Return the time to this closure as a six-character string.

        """
        time = Time.now() if time is None else time
        ttc = (self.start - time)
        if string:
            return human_time(ttc.to("s").value)
        return ttc

class Closure(Window):
    """An LCH closure window, when the laser cannot be propogated."""
    propogate = False

class Opening(Window):
    """An LCH Opening window, when the laser can be propogated."""
    propogate = True

    def __init__(self, region, start, end):
        super(Opening, self).__init__(region, start, end)
        self.region.add(self)

    @classmethod
    def from_openings_data(cls, region, data):
        """From already parsed openings data."""
        obj = cls(region, data["Start"], data["End"])
        if not np.abs(obj.duration - data["Opening"]) < 1.0 * u.second:
            raise ValueError("Duration of opening {} does not match parsed value: {}".format(obj, data["Opening"]))
        return obj

    @classmethod
    def from_openings(cls, region, text, date=None):
        """Parse a line of openings to create this object.

        :param region: The region to which this closure applied.
        :param text: The line of text for the closure.
        :param date: The date of the closure.

        """
        data = parse_openings_line(text, date)
        return cls.from_openings_data(region, data)

class _LCHParser(object):
    """Closure Parsing State"""

    def __init__(self, name, date=None, regioncls=None):
        super(_LCHParser, self).__init__()
        self._name = name
        self._date = Time(date)
        self._regioncls = regioncls if regioncls is not None else Regions
        self.reset()

    @property
    def date(self):
        """Return the active date."""
        return self._date

    def reset(self):
        """Reset the parser"""
        self._next_line = "header"
        self._this_line = None
        self._regions = self._regioncls()
        self._current_region = None
        self._line_number = 0

    def _warn(self, message):
        """Emit a formatted warning."""
        template = "[{name:s} line {ln:d}][{this:s}]{message:s}"
        warnings.warn(template.format(
            name = self._name,
            ln = self._line_number,
            this = self._this_line,
            message = message
        ))

    def header(self, text):
        """Handle the header line."""
        self._next_line = "blank"

    def blank(self, text):
        """Handle a blank line."""
        self._next_line = "starlist"

    def summary(self, text):
        """Handle a summary line."""
        self._next_line = "unknown"
        data = parse_summary_line(text)
        self._regions.verify(**data)

    def starlist(self, text):
        """Parse a starlist line."""
        self._next_line = "opening"
        self._current_region = Region.from_starlist(text)
        self._regions[self._current_region.name] = self._current_region

    def opening(self, text):
        """Parse an opening line."""
        self._next_line = "unknown"
        self._opening = Opening.from_openings(self._current_region, text, date=self._date)
        self._opening_text = text

    def handle_last_opening(self):
        """Handle the most recent opening."""
        if hasattr(self, '_opening_text'):
            final_closure = re.search(r"Closure\(sec\)[\s]+([\d]+)$", self._opening_text)
            if final_closure:
                next_start = self._opening.end + (int(final_closure.group(1)) * u.second)
                final_end = self._current_region.openings[0].start + 12 * u.hour
                self._opening = Opening(self._current_region, next_start, final_end)
            delattr(self, '_opening_text')

    def __call__(self, stream):
        """Parse a stream line."""
        self.reset()
        for line in stream:
            self.parse_line(line)
        regions = self._regions
        self.reset()
        return regions

    def parse_line(self, text):
        """Parse a string/file."""
        text = text.strip()
        self._line_number += 1
        if self._next_line == "unknown":
            if re.match(r"^[\s]*$", text):
                self.handle_last_opening()
                self._this_line = "blank"
                self.blank(text)
            elif re.match(r"^[\s]*First", text):
                self.handle_last_opening()
                self._this_line = "summary"
                self.summary(text)
            else:
                self._this_line = "opening"
                self.opening(text)
        else:
            self._this_line = self._next_line
            getattr(self, self._this_line)(text)

    @classmethod
    def parse_file(cls, filename, date=None):
        """Parse a whole closure file."""
        parser = cls(filename, date)
        with open(filename, 'r') as stream:
            return parser(stream)

class _LCHParserLick(_LCHParser):
    """Parser for Lick-style .lsm files."""

    def load_starlist(self, filename):
        """Load individual starlists."""
        for line in read_skip_comments(filename):
            region = Region.from_starlist(line)
            self._regions[region.name] = region

    def parse_lsmfile(self, filename):
        """Parse a single file."""
        region_name = os.path.splitext(os.path.basename(filename))[0]
        self._this_line = 'lsm file'
        if region_name not in self._regions:
            self._warn("Skipping region '{}', not found in starlist.".format(region_name))
            return
        self._current_region = self._regions[region_name]
        with open(filename, 'r') as stream:
            self._this_line = 'opening'
            for line in stream:
                opening = Opening.from_openings_data(self._current_region, parse_openings_compact_format(line, date=self.date))

    @staticmethod
    def glob_files(files):
        """Glob for filenames if nescessary."""
        if files is None:
            return glob.iglob(os.path.join(os.path.getcwd(),"*.lsm"))
        if isinstance(files, six.text_type):
            return glob.iglob(files)
        return iter(files)

    def __call__(self, starlist, files=None):
        """Parse a stream line."""
        self.reset()
        self.load_starlist(starlist)
        for fn, file in enumerate(self.glob_files(files)):
            self._line_number = fn
            self.parse_lsmfile(file)
        return self._regions

    @classmethod
    def parse_file(cls, starlist, filenames, date=None):
        """Parse a whole closure file."""
        parser = cls(starlist, date)
        parser(starlist, filenames)

_openings_re = re.compile(
    r"""
    ^[\s]*(?P<Start>[\d]{2}:[\d]{2}:[\d]{2}) # Start time
    [\s]*(?P<End>[\d]{2}:[\d]{2}:[\d]{2})    # End time
    [\s]*open\(min:sec\)\ (?P<Opening>[\d]{2,4}:[\d]{2}) # Opening duration
    [\s]*(?:Closure\(sec\)[\s]*(?P<Closure>[\d]+))? # Closure duration
    """,
    re.VERBOSE
)

def parse_openings_line(line, date=None):
    """Parse an openings line.

    :param line: The text line to parse as an opening.
    :param date: The date for this openings line.
    """

    # Setup the reference date
    date = Time.now() if date is None else date
    date.out_subfmt = 'date'
    date_str = date.iso

    # Match the regular expression.
    match = _openings_re.match(line)
    if not match:
        raise ValueError("Can't parse line as opening: '{}'".format(line))
    data = match.groupdict("")

    # Handle date types
    for time_key in "Start End".split():
        data[time_key] = Time("{date:s} {time:s}".format(date=date_str, time=data[time_key]), scale='utc')

    # Handle quantity types
    for delta_key in "Opening Closure".split():
        if delta_key in data and len(data[delta_key]):
            if ":" in data[delta_key]:
                minutes, seconds = map(float,data[delta_key].split(":"))
            else:
                minutes, seconds = 0, float(data[delta_key])
            data[delta_key] = minutes * u.minute + seconds * u.second

    return data

_openings_compact_re = re.compile(r"""
^[\s]*(?P<Start>[\d]{10})[\s]+(?P<End>[\d]{10})[\s]+(?P<Opening>[\d]+)[\s]+(?P<Closure>[\d]+)[\s]*$
""", re.VERBOSE)

def parse_openings_compact_format(line, date=None):
    """Parse openings which are in the short format."""
    # Setup the reference date
    date = Time.now() if date is None else date
    date.out_subfmt = 'date'
    date_str = date.iso

    # Match the regular expression.
    match = _openings_compact_re.match(line)
    if not match:
        raise ValueError("Can't parse line as opening: '{}'".format(line))
    data = match.groupdict("")
    # Handle date types
    for time_key in "Start End".split():
        data[time_key] = Time(datetime.datetime.utcfromtimestamp(int(data[time_key])), scale='utc')

    # Handle quantity types
    for delta_key in "Opening Closure".split():
        if delta_key in data and len(data[delta_key]):
            if ":" in data[delta_key]:
                minutes, seconds = map(float,data[delta_key].split(":"))
            else:
                minutes = 0
                seconds = float(data[delta_key])
            data[delta_key] = minutes * u.minute + seconds * u.second

    return data

_summary_line_re = re.compile(r"""
    ^[\s]*(?P<ordinal>[\w]+)[\s]+(?P<objects>[\d]+)[\s]+objects\:[\s]+ # Number of objects
    # Only one of the two following options will appear in a closure line.
    (?:(?:total\ time\ of\ closures:[\s]+(?P<time>[\d]+(?:\.[\d]+)?))| # Closure time
    (?:total\ number\ of\ closures:[\s]+(?P<number>[\d]+))) # Number of closures

    [\s]*$ # End the line.
    """, re.VERBOSE)

def parse_summary_line(line):
    """Parse a summary line."""
    match = _summary_line_re.match(line)
    if not match:
        raise ValueError("Can't parse line as summary: '{}'".format(line))
    data = match.groupdict()

    # Handle Time
    if data["time"] is not None:
        data["time"] = float(data["time"]) * u.second
    else:
        del data["time"]

    # Handle Number
    if data["number"] is not None:
        data["number"] = int(data["number"])
    else:
        del data["number"]

    data["objects"] = int(data["objects"])
    return data

def parse_closures_list(stream, name="<stream>", date=None, cls=_LCHParser):
    """Parse a closure list"""
    parser = cls(name, date)
    return parser(stream).values()

class Regions(collections.OrderedDict):
    """A dictionary interface for LCH regions.

    The dictionary keys are the region names, taken from the starlist formatted lines.
    The dictionary values are instances of :class:`Region`. Use the :meth:`Regions.parse`
    method to parse a closure file into region objects.

    """
    def __init__(self, *args, **kwargs):
        super(Regions, self).__init__(*args, **kwargs)

    def __repr__(self):
        """Represent this summary group."""
        repr_str = ["<{0:s}".format(self.__class__.__name__)]
        repr_str += ["N={0:d}".format(self.n_objects)]
        repr_str += ["({0:d} for {0:.0f})".format(self.n_closures, self.t_closures)]
        return " ".join(repr_str) + ">"

    def contains(self, location):
        """Check if this location is covered by any of the regions here."""
        return any(region.contains(location) for region in self.values())

    def open(self, location, time):
        """Check if this location and time is open for laser propagation.

        This method uses a conservative algorithm. Locations which are not
        contained in any LCH region are declared closed. Locations which are
        not declared open in all LCH regions where they fall are also declared
        closed. LCH region size is set by the :attr:`Region.RADIUS` attribute,
        and defaults to the 2 arcminute size which is standard at Keck.

        Parameters
        ----------
        location : :class:`~astropy.coordinates.SkyCoord`
            The location to check against all of the contained regions.
        time : :class:`~astropy.time.Time`
            The time to check against all contained regions.

        Returns
        -------
        open : bool
            Whether this location, at the specified time, is open for laser
            propogation.

        """
        contained = False
        open = True
        for region in self.values():
            this_contain = region.contains(location)
            contained |= this_contain
            if this_contain:
                this_open = region.open(time)
                if not this_open:
                    return False
                open &= this_open
        return (contained and open)

    def closed(self, location, time):
        """Check if this location and time is closed. Returns the opposite of :meth:`Regions.open`"""
        return not self.open(location, time)

    def to_starlist(self, filename, **kwargs):
        """Write these regions to a starlist file.

        Parameters
        ----------
        filename : str or fileobj
            Filename for writing

        """
        if isinstance(filename, io.IOBase) or hasattr(filename, 'write'):
            stream = filename
            self._to_starlist_stream(stream, **kwargs)
        else:
            with open(filename, mode) as stream:
                self._to_starlist_stream(stream, **kwargs)

    def _to_starlist_stream(self, file, **kwargs):
        """Create a full starlist."""
        for t in self.values():
            file.write(t.to_starlist(**kwargs))
            file.write("\n")

    @property
    def n_objects(self):
        """Number of regions."""
        return len(self)

    @property
    def n_closures(self):
        """Number of closures in all regions."""
        return sum([ len(region.openings) - 1 for region in self.values() ])

    @property
    def t_closures(self):
        """Total duration of closures in all regions."""
        return sum([ sum([ closure.duration.to('s').value for closure in region.closures ]) for region in self.values() ]) * u.second

    @classmethod
    def parse(cls, filename, date=None, _parser=_LCHParser, **kwargs):
        """Parse a closure file or stream into regions.

        Closure files contain starlists lines followed by listings
        of opening blocks when laser propogation is permitted. This
        object parses these into individual :class:`Region` objects
        which contain a list of :class:`Opening` periods.

        Parameters
        ----------
        filename : str or fileobj
            Filename or readable file object for the parser.
        date : :class:`~astropy.time.Time`
            String, or Time value which will be used to set the
            date for the closure file.

        Returns
        -------
        regions : :class:`Regions`
            The parsed regions.

        """
        date = Time.now() if date is None else date
        with get_readable_fileobj(filename) as stream:
            parser = _parser(stream.name, date, regioncls=cls)
            return parser(stream, **kwargs)

def upcoming_closures(regions, limit = 1 * u.hour, when = Time.now()):
    """Yield upcoming closures."""
    upcoming = []
    for region in regions:
        for closure in region.closures:
            till_start = closure.start - when
            till_finish = closure.end - when
            if (till_start > 0 * u.s and till_start < limit) or (till_finish > 0 * u.s and till_finish < limit):
                upcoming.append(closure)
    upcoming.sort(key = lambda o : o.start)
    return upcoming

max_stars_repo_name: alexrudy/KOPy

max_issues_repo_name: alexrudy/KOPy

max_forks_repo_name: alexrudy/KOPy

reason: import numpy,import astropy,from astropy


text:
import numpy as np
from skmultiflow.transform.base_transform import StreamTransform
from skmultiflow.utils.utils import get_dimensions

class OneHotToCategorical(StreamTransform):
    """ OneHotToCategorical

    Transform that receives a features matrix, with some features coded as
    one hot, and transform these features in basic int coded categorical
    features.

    Parameters
    ----------
    categorical_list: list of lists
        Each inner list contains all the attribute indexes that are associated with
        the same categorical feature.

    """
    def __init__(self, categorical_list):
        super().__init__()
        self.categorical_list = categorical_list

    def get_info(self):
        return 'OneHotToCategorical: categorical_list: ' + str(self.categorical_list)

    def transform(self, X):
        """ transform

        Transform one hot features in the X matrix into int coded
        categorical features.

        Parameters
        ----------
        X: numpy.ndarray of shape (n_samples, n_features)
            The sample or set of samples that should be transformed.

        Returns
        -------
        numpy.ndarray
            The transformed data.

        """
        r, c = get_dimensions(X)

        new_width = c
        for i in range(len(self.categorical_list)):
            new_width -= len(self.categorical_list[i]) - 1

        ret = np.zeros((0, new_width), dtype=X.dtype)
        for i in range(r):
            ret = np.concatenate((ret, self._transform(X[i, :], new_width)), axis=0)

        return ret

    def _transform(self, X, new_size):
        ret = np.zeros((1,new_size), dtype=X.dtype)
        list_index = 0
        i = 0
        new_i = 0
        while i < X.size:
            if i in self.categorical_list[list_index]:
                c = 0
                found = False
                for j in range(len(self.categorical_list[list_index])):
                    if X[i+j] == 1:
                        found = True
                    if not found:
                        c += 1
                i += len(self.categorical_list[list_index])
                list_index += 1
                ret[0][new_i] = c
                new_i += 1
            else:
                ret[0][new_i] = X[i]
                new_i += 1
                i += 1

        return ret

    def fit(self, X, y):
        return self

    def partial_fit_transform(self, X, y=None, classes=None):
        return self.transform(X)

    def partial_fit(self, X, y=None, classes=None):
        return self

max_stars_repo_name: AndreFCruz/scikit-multiflow

max_issues_repo_name: AndreFCruz/scikit-multiflow

max_forks_repo_name: AndreFCruz/scikit-multiflow


text:
import numpy as np
import pandas as pd
import sys

def duplicates(a):
	print(len(a), len(set(a)))
	return len(a) != len(set(a))

def intersection(a, b):
	return not set(a).isdisjoint(b)

trainfile = sys.argv[1]
testfile =  sys.argv[2]

train = pd.read_csv(trainfile)
trainsequences = list(train['sequence'])

test = pd.read_csv(testfile)
testsequences = list(test['sequence'])

print(duplicates(trainsequences))
print(duplicates(testsequences))
print(intersection(trainsequences, testsequences))

max_stars_repo_name: lelange/cu-ssp

max_issues_repo_name: lelange/cu-ssp

max_forks_repo_name: lelange/cu-ssp


text:
import re
import random
import time
import math
import numpy as np
import nltk
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
plt.style.use('default')

class AttentionModule(nn.Module):
    """Realizes a very simple attention module. Takes last state and word into account with two-layer feedforward"""

    def __init__(self, embedding_dim, state_dim, hidden_dim):
        super(AttentionModule, self).__init__()

        self.embedding_dim = embedding_dim
        self.hidden_dim = hidden_dim
        self.state_dim = state_dim
        self.output_layer = nn.Sequential(OrderedDict([
            ('fc1', nn.Linear(embedding_dim+state_dim, hidden_dim)),
            ('tanh', nn.Tanh()),
            ('dropout', nn.Dropout(0.25)),
            ('fc2', nn.Linear(hidden_dim, 1)),
            ('sigmoid', nn.Sigmoid())
        ]))

        self.reset_parameters()

    def reset_parameters(self):
        """This is PyTorch's default initialization method"""
        stdv = 1.0 / math.sqrt(self.hidden_dim)
        for weight in self.parameters():
            weight.data.uniform_(-stdv, stdv)

    def forward(self, word_embed, last_state):

        input_tensor = torch.cat([word_embed, last_state], dim=1)

        for layer in self.output_layer:
            input_tensor = layer(input_tensor)

        return input_tensor

class IterativeLSTMClassifier(nn.Module):
    """Encodes sentence with an LSTM and projects final hidden state"""

    def __init__(self, vocab_size, embedding_dim, hidden_dim, output_dim, vocab, iterations):
        super(IterativeLSTMClassifier, self).__init__()
        self.vocab = vocab
        self.hidden_dim = hidden_dim
        self.embedding_size = embedding_dim
        self.embed = nn.Embedding(vocab_size, embedding_dim, padding_idx=1)
        self.rnn = MyLSTMCell(embedding_dim, hidden_dim)
        self.attention_module = AttentionModule(embedding_dim, hidden_dim, embedding_dim)
        self.iterations = iterations

        self.output_layer = nn.Sequential(
            nn.Dropout(p=0.25),
            nn.Linear(hidden_dim, output_dim)
        )

    def forward(self, x):

        B = x.size(0)  # batch size (this is 1 for now, i.e. 1 single example)
        T = x.size(1)  # time (the number of words in the sentence)

        input_ = self.embed(x)

        # here we create initial hidden states containing zeros
        # we use a trick here so that, if input is on the GPU, then so are hx and cx
        hx = input_.new_zeros(B, self.rnn.hidden_size)
        cx = input_.new_zeros(B, self.rnn.hidden_size)

        # process input sentences one word/timestep at a time
        # input is batch-major, so the first word(s) is/are input_[:, 0]
        for iter_ind in range(self.iterations):
            if iter_ind > 0:
                # hx has dimension B,hidden_size -> need to be B,T,hidden_size
                exp_hx = hx.unsqueeze(dim=1).repeat([1, T, 1])
                attent_hx_input = exp_hx.view((B * T, self.rnn.hidden_size))
                attent_word_input = input_.view((B * T, self.embedding_size))
                print(attent_hx_input)
                print(attent_word_input)
                attent_scores = self.attention_module(attent_word_input, attent_hx_input)
                attent_scores = attent_scores.reshape(B, T, 1)

            outputs = []
            for i in range(T):
                hx_new, cx_new = self.rnn(input_[:, i], (hx, cx))
                if iter_ind == 0:
                    hx = hx_new
                    cx = cx_new
                else:
                    hx = attention_score[:, i] * hx_new + (1 - attention_score[:, i]) * hx
                    cx = attention_score[:, i] * cx_new + (1 - attention_score[:, i]) * cx
                outputs.append(hx)

            # if we have a single example, our final LSTM state is the last hx
            if B == 1:
                last_final_states = hx
            else:
                #
                # This part is explained in next section, ignore this else-block for now.
                #
                # we processed sentences with different lengths, so some of the sentences
                # had already finished and we have been adding padding inputs to hx
                # we select the final state based on the length of each sentence

                # two lines below not needed if using LSTM form pytorch
                outputs = torch.stack(outputs, dim=0)  # [T, B, D]
                outputs = outputs.transpose(0, 1).contiguous()  # [B, T, D]

                # to be super-sure we're not accidentally indexing the wrong state
                # we zero out positions that are invalid
                pad_positions = (x == 1).unsqueeze(-1)

                outputs = outputs.contiguous()
                outputs = outputs.masked_fill_(pad_positions, 0.)

                mask = (x != 1)  # true for valid positions [B, T]
                lengths = mask.sum(dim=1)  # [B, 1]

                indexes = (lengths - 1) + torch.arange(B, device=x.device, dtype=x.dtype) * T
                last_final_states = outputs.view(-1, self.hidden_dim)[indexes]  # [B, D]

            hx = last_final_states
            cx = input_.new_zeros(B, self.rnn.hidden_size)

        # we use the last hidden state to classify the sentence
        logits = self.output_layer(last_final_states)
        return logits

torch.cuda.empty_cache()
lstm_model = IterativeLSTMClassifier(
    len(v.w2i), 300, 168, len(t2i), v, 2)

# copy pre-trained vectors into embeddings table
with torch.no_grad():
    lstm_model.embed.weight.data.copy_(torch.from_numpy(vectors))
    lstm_model.embed.weight.requires_grad = False

print(lstm_model)
print_parameters(lstm_model)

lstm_model = lstm_model.to(device)

batch_size = 32
optimizer = optim.Adam(lstm_model.parameters(), lr=2e-4)

lstm_losses, lstm_accuracies = train_model(
    lstm_model, optimizer, num_iterations=30000,
    print_every=250, eval_every=250,
    batch_size=batch_size,
    batch_fn=get_minibatch,
    prep_fn=prepare_minibatch,
    eval_fn=evaluate)

max_stars_repo_name: sjoerdapp/colab

max_issues_repo_name: sjoerdapp/colab

max_forks_repo_name: sjoerdapp/colab


text:
# Similar to scratch7, but with the BAN Channel
from sys_simulator.general import db_to_power, power_to_db
from sys_simulator.q_learning.rewards import dis_reward_tensor_db
from scipy.spatial.distance import euclidean
from sys_simulator.channels import BANChannel
from sys_simulator.plots import plot_positions_actions_pie
from sys_simulator.q_learning.environments.completeEnvironment5dB \
    import CompleteEnvironment5dB
from sys_simulator.dqn.agents.dqnAgent import ExternalDQNAgent
from sys_simulator.dqn.externalDQNFramework import ExternalDQNFramework
from sys_simulator.parameters.parameters \
    import EnvironmentParameters, TrainingParameters, DQNAgentParameters
from matplotlib import pyplot as plt
import os
import torch
import numpy as np
import math

n_mues = 1  # number of mues
n_d2d = 2  # number of d2d pairs
n_rb = n_mues   # number of RBs
bs_radius = 500  # bs radius in m
rb_bandwidth = 180*1e3  # rb bandwidth in Hz
d2d_pair_distance = 50  # d2d pair distance in m
p_max = 23  # max tx power in dBm
noise_power = -116  # noise power per RB in dBm
bs_gain = 17    # macro bs antenna gain in dBi
user_gain = 4   # user antenna gain in dBi
sinr_threshold_mue = 6  # true mue sinr threshold in dB
mue_margin = 200  # mue sinr margin in dB
# conversions from dBm to dB
p_max = p_max - 30
noise_power = noise_power - 30
# q-learning parameters
STEPS_PER_EPISODE = 10
EPSILON_MIN = 0.01
EPSILON_DECAY = 100 * EPSILON_MIN / STEPS_PER_EPISODE
MAX_NUM_EPISODES = int(1.2/EPSILON_DECAY)
MAX_NUMBER_OF_AGENTS = 10
ALPHA = 0.05  # Learning rate
GAMMA = 0.98  # Discount factor
C = 80  # C constant for the improved reward function
TARGET_UPDATE = 10
REPLAY_MEMORY_SIZE = 10000
BATCH_SIZE = 256
HIDDEN_SIZE = 64
NUM_HIDDEN_LAYERS = 5
# channel parameters
CHANNEL_RND = False
# more parameters
cwd = os.getcwd()
# params objects
env_params = EnvironmentParameters(
    rb_bandwidth, d2d_pair_distance, p_max, noise_power,
    bs_gain, user_gain, sinr_threshold_mue, n_mues,
    n_d2d, n_rb, bs_radius, c_param=C, mue_margin=mue_margin
)
train_params = TrainingParameters(MAX_NUM_EPISODES, STEPS_PER_EPISODE)
agent_params = DQNAgentParameters(
    EPSILON_MIN, EPSILON_DECAY, 1,
    REPLAY_MEMORY_SIZE, BATCH_SIZE, GAMMA
)
# actions, rewards, environment, agent
reward_function = dis_reward_tensor_db
channel = BANChannel(rnd=CHANNEL_RND)
env = CompleteEnvironment5dB(env_params, reward_function, channel)
foo_agent = ExternalDQNAgent(agent_params, [1])
env.build_scenario([foo_agent])
env_state_size = env.get_state_size(foo_agent)
actions = power_to_db(np.linspace(
    db_to_power(p_max-20), db_to_power(p_max-10), 10
))
actions[0] = -1000
framework = ExternalDQNFramework(
    agent_params,
    env_state_size,
    len(actions),
    HIDDEN_SIZE,
    NUM_HIDDEN_LAYERS
)
framework.policy_net.load_state_dict(
    torch.load(f'{cwd}/models/dql/script32.pt')
)
pairs_positions = [
    (250, 0),
    (-250, 0),
    (0, 250),
    (0, -250)
]
mue_position = 250 / math.sqrt(2) * (1, 1)
n_agents = len(pairs_positions)
episode_steps = STEPS_PER_EPISODE

def calculate_interferences(env: CompleteEnvironment5dB):
    bs = env.bs
    mue = env.mue
    d2d_pairs = env.d2d_pairs
    txs = [mue]
    txs += [p[0] for p in d2d_pairs]
    rxs = [bs]
    rxs += [p[1] for p in d2d_pairs]
    interferences = np.zeros((len(txs), len(rxs)))
    for i, tx in enumerate(txs):
        for j, rx in enumerate(rxs):
            interf = \
                tx.tx_power + tx.gain \
                - env.channel.step(euclidean(tx.position, rx.position)) \
                + rx.gain
            interferences[i][j] = interf
    tx_labels = [d.id for d in txs]
    rx_labels = [d.id for d in rxs]
    return interferences, tx_labels, rx_labels

def run():
    # global actions
    agents = [ExternalDQNAgent(agent_params, actions)
              for _ in range(n_agents)]  # 1 agent per d2d tx
    env.set_scenario(pairs_positions, mue_position, agents)
    obs = [env.get_state(a) for a in agents]
    total_reward = 0.0
    actions_index = list()
    for _ in range(episode_steps):
        actions_index = list()
        for j, agent in enumerate(agents):
            aux = agent.act(framework, obs[j].float()).max(1)
            actions_index.append(aux[1].item())
            agent.set_action(aux[1].long(), agent.actions[aux[1]])
        next_obs, rewards = env.step(agents)
        obs = next_obs
        total_reward += sum(rewards)
    d2d_txs, d2d_rxs = zip(*env.d2d_pairs)
    # D2D interference on the MUE, in dB
    d2d_interferences = np.array([
        d.tx_power + env.params.user_gain + env.params.bs_gain -
        env.channel.step(d.distance_to_bs) for d in d2d_txs
    ])
    d2d_interferences_mag = db_to_power(d2d_interferences)
    d2d_total_interference = np.sum(d2d_interferences_mag)
    percentage_interferences = d2d_interferences_mag / d2d_total_interference
    interferences, tx_labels, rx_labels = calculate_interferences(env)
    if d2d_total_interference != 0:
        plot_positions_actions_pie(
            env.bs, env.mue, d2d_txs, d2d_rxs,
            actions_index, percentage_interferences,
            obs[0][0][4].item(), sinr_threshold_mue,
            env.reward, interferences, tx_labels, rx_labels
        )
    print_stuff(actions, env)
    plt.show()

def print_stuff(actions, env: CompleteEnvironment5dB):
    actions = [f'{i:.2f}' for i in actions]
    sinr_d2ds = [f'{d[0].sinr:.2f}' for d in env.d2d_pairs]
    print(f'MUE Tx Power [dBW]: {env.mue.tx_power:.2f}')
    print(f'D2D Power levels [dBW]: {actions}')
    print(f'D2D SINR [dB]: {sinr_d2ds}')
    print(f'D2D Spectral Efficiencies: {env.d2d_spectral_eff}')

if __name__ == '__main__':
    run()

reason: import numpy,from scipy


text:
#!/usr/bin/env python3

import sys
import numpy as np
from numpy import genfromtxt

dat = genfromtxt( sys.argv[ 1 ], delimiter=',')
fullmtx = dat.reshape(3,4)
mtx = fullmtx[:,0:3]
print( np.linalg.det( mtx ))

max_stars_repo_name: mzouink/template-building

max_issues_repo_name: mzouink/template-building

max_forks_repo_name: mzouink/template-building

reason: import numpy,from numpy