text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Remade on Sun 30 May 2021 22:32

@author: Laurens Stoop - l.p.stoop@uu.nl

"""

#%%
# =============================================================================
# Dependencies
# =============================================================================

## Importing modules
import xarray as xr
import numpy as np
# import regionmask
# import geopandas as gpd
import datetime
import pandas as pd
# import matplotlib.pyplot as plt
import os.path

# Select the years to run
years = np.array([
            # '1950', '1951', '1952',
            # '1953', '1954', '1955',
            # '1956', '1957', '1958',
            # '1959', '1960', '1961',
            # '1962', '1963', '1964',
            # '1965', '1966', '1967',
            # '1968', '1969', '1970',
            # '1971', '1972', '1973',
            # '1974', '1975', '1976',
            # '1977', '1978',
            # '1979', '1980', '1981',
            # '1982', '1983', '1984',
            # '1985', '1986', '1987',
            # '1988', '1989', '1990',
            # '1991', '1992', '1993',
            # '1994', '1995', '1996',
            # '1997', '1998', '1999',

            # '2000', '2001', '2002',
            # '2003', '2004', '2005',
            # '2006', '2007', '2008',
            # '2009', '2010', '2011',
            # '2012', '2013',
            '2014',
            '2015', '2016', '2017',
            '2018', '2019',
            #'2020'
        ])

# Versions of the TYNDP
CD_TYNDP_input = np.array([
            'DE_2030',
            'DE_2040',
            'GA_2030',
            'GA_2040',
            'NT_2025',
            'NT_2030',
            'NT_2040'
        ])

# Set the path for the data
PATH_TO_TYNDP = '/media/DataStager1/Other/CapacityDistribution/TYNDP/Origin/'

# Read NetCDF
FOLDER_ERA5_CF_NUTS0 = '/media/DataStager2/ERA5-EU_CF/NUTS0/'
FOLDER_EV_NUTS0 = '/media/DataStager2/ERA5-EU_EV/NUTS0/'

print('NOTIFY: Initialization is complete, Skynet active')
#%%
# =============================================================================
# Load in the datafiles with the capacity distributions
# =============================================================================

# Select the distribution to run over
# capacity_distribution = CD_TYNDP_input[0]

for capacity_distribution in CD_TYNDP_input:

    print('NOTIFY: Working on Distribution '+capacity_distribution+'!')

    # Read in the Capacity Distribution from the TYNDP
    df_cd_tyndp = pd.read_csv(PATH_TO_TYNDP+'TYNDP-'+capacity_distribution+'.csv', sep=';', names=['Code', 'Type', 'Capacity'])

    #%%
    # =============================================================================
    # Load in the NUTS0 data
    # =============================================================================

    # Set the year to run over
    for year in years:

        print('NOTIFY: Working on year '+year+'!')

        # Load in the NetCDF
        ds_cf_solar = xr.open_dataset(FOLDER_ERA5_CF_NUTS0+'ERA5-EU_CF-NUTS0_solar_'+str(year)+'.nc') #, chunks = {'time': 8760})
        ds_cf_windoff = xr.open_dataset(FOLDER_ERA5_CF_NUTS0+'ERA5-EU_CF-NUTS0_windoff_'+str(year)+'.nc') #, chunks = {'time': 8760})
        ds_cf_windon = xr.open_dataset(FOLDER_ERA5_CF_NUTS0+'ERA5-EU_CF-NUTS0_windon_'+str(year)+'.nc') #, chunks = {'time': 8760})

        #%%
        # =============================================================================
        # Multiply the capacity distribution with the capacity factor
        # =============================================================================

        # Set a new dataset for energy variables
        ds_ev_solar = xr.Dataset()
        ds_ev_windoff = xr.Dataset()
        ds_ev_windon = xr.Dataset()

        # The easy countries
        countries = [
            'BA', 'BE', 'DE',
            'EE', 'ES', 'FI',
            'HR', 'IE',
            'LT', 'LU',
            'LV', 'NL', 'SI',
            'AT', 'CH', 'CZ',
            'HU', 'PL',
            'PT', 'SK', 'TR',
            'AL', 'BG', 'CY',
            'ME', 'MK',
            'RO', 'RS',
            'MT', # to small but used offshore data for onshore figures (small capacity)
            'FR', # without the island of corsica
            'UK', # without the region of northern ireland

            ]

        # The easy countries
        for country in countries:

            #%% working on solar
            # Define the capacity installed
            country_cap_distr_solar = df_cd_tyndp.where(df_cd_tyndp['Code'] == country).where(df_cd_tyndp['Type'] == 'Solar PV').dropna()['Capacity'].values

            # If this capacity is not defined, do not calculate
            if country_cap_distr_solar.size == 0:
                print('There is no solar capacity for '+country)

            # Apply the wind offshore capacity distribution
            else:
                # apply the cap distribution
                ds_ev_solar[country] = country_cap_distr_solar * ds_cf_solar[country+'00']

            #%% working on onshore wind
            # Define the capacity installed
            country_cap_distr_windon = df_cd_tyndp.where(df_cd_tyndp['Code'] == country).where(df_cd_tyndp['Type'] == 'Onshore Wind').dropna()['Capacity'].values

            # If this capacity is not defined, do not calculate
            if country_cap_distr_windon.size == 0:
                print('There is no onshore wind capacity for '+country)

            # Apply the wind offshore capacity distribution
            else:
                # apply the cap distribution
                ds_ev_windon[country] = country_cap_distr_windon * ds_cf_windon[country+'00']

            #%% working on offshore wind
            # Define the capacity installed
            country_cap_distr_windoff = df_cd_tyndp.where(df_cd_tyndp['Code'] == country).where(df_cd_tyndp['Type'] == 'Offshore Wind').dropna()['Capacity'].values

            # If this capacity is not defined, do not calculate
            if country_cap_distr_windoff.size == 0:
                print('There is no offshore capacity for '+country)

            # Apply the wind offshore capacity distribution
            else:
                ds_ev_windoff[country] = country_cap_distr_windoff * ds_cf_windoff[country+'00_OFF']

        #%%
        # =============================================================================
        # Working on Greece (just a naming issue)
        # =============================================================================
        country_cap_distr_windoff = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'EL').where(df_cd_tyndp['Type'] == 'Offshore Wind').dropna()['Capacity'].values

        # If this capacity is not defined, do not calculate
        if country_cap_distr_windoff.size == 0:
            print('There is no offshore capacity for Greece')
        else:
            ds_ev_windoff['EL'] = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'EL').where(df_cd_tyndp['Type'] == 'Offshore Wind').dropna()['Capacity'].values * ds_cf_windoff['GR00_OFF']

        ds_ev_windon['EL'] = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'EL').where(df_cd_tyndp['Type'] == 'Onshore Wind').dropna()['Capacity'].values * ds_cf_windon['GR00']
        ds_ev_solar['EL'] = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'EL').where(df_cd_tyndp['Type'] == 'Solar PV').dropna()['Capacity'].values * ds_cf_solar['GR00']

        # =============================================================================
        # Working on Ukraine, without the oblast Moekatsjevo
        # =============================================================================
        ds_ev_windon['UA'] = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'UA').where(df_cd_tyndp['Type'] == 'Onshore Wind').dropna()['Capacity'].values * ds_cf_windon['UA02']

        # =============================================================================
        # Working on Sweden, capacity devided based on CF's of regions
        # =============================================================================

        # Determin the devisors
        SE_devisor_windoff = ds_cf_windoff['SE01_OFF'].mean(...).values + ds_cf_windoff['SE02_OFF'].mean(...).values + ds_cf_windoff['SE03_OFF'].mean(...).values + ds_cf_windoff['SE04_OFF'].mean(...).values
        SE_devisor_windon = ds_cf_windon['SE01'].mean(...).values + ds_cf_windon['SE02'].mean(...).values + ds_cf_windon['SE03'].mean(...).values + ds_cf_windon['SE04'].mean(...).values
        SE_devisor_solar = ds_cf_solar['SE01'].mean(...).values + ds_cf_solar['SE02'].mean(...).values + ds_cf_solar['SE03'].mean(...).values + ds_cf_solar['SE04'].mean(...).values

        # Find the capacities
        SE_cap_off = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'SE').where(df_cd_tyndp['Type'] == 'Offshore Wind').dropna()['Capacity'].values[0]
        SE_cap_on = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'SE').where(df_cd_tyndp['Type'] == 'Offshore Wind').dropna()['Capacity'].values[0]
        SE_cap_sol = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'SE').where(df_cd_tyndp['Type'] == 'Solar PV').dropna()['Capacity'].values[0]

        # Calculate the regional energy generation for offshore wind
        SE01_off = SE_cap_off * ds_cf_windoff['SE01_OFF'] * ds_cf_windoff['SE01_OFF'].mean(...).values * SE_devisor_windoff**-1
        SE02_off = SE_cap_off * ds_cf_windoff['SE02_OFF'] * ds_cf_windoff['SE02_OFF'].mean(...).values * SE_devisor_windoff**-1
        SE03_off = SE_cap_off * ds_cf_windoff['SE03_OFF'] * ds_cf_windoff['SE03_OFF'].mean(...).values * SE_devisor_windoff**-1
        SE04_off = SE_cap_off * ds_cf_windoff['SE04_OFF'] * ds_cf_windoff['SE04_OFF'].mean(...).values * SE_devisor_windoff**-1
        # Sum over all regions to obtain national figures
        ds_ev_windoff['SE'] = SE01_off + SE02_off + SE03_off + SE04_off

        # Calculate the regional energy generation for onshore wind
        SE01_on = SE_cap_on * ds_cf_windon['SE01'] * ds_cf_windon['SE01'].mean(...).values * SE_devisor_windon**-1
        SE02_on = SE_cap_on * ds_cf_windon['SE02'] * ds_cf_windon['SE02'].mean(...).values * SE_devisor_windon**-1
        SE03_on = SE_cap_on * ds_cf_windon['SE03'] * ds_cf_windon['SE03'].mean(...).values * SE_devisor_windon**-1
        SE04_on = SE_cap_on * ds_cf_windon['SE04'] * ds_cf_windon['SE04'].mean(...).values * SE_devisor_windon**-1
        # Sum over all regions to obtain national figures
        ds_ev_windon['SE'] = SE01_on + SE02_on + SE03_on + SE04_on

        # Calculate the regional energy generation for solar PV
        SE01_sol = SE_cap_sol * ds_cf_solar['SE01'] * ds_cf_solar['SE01'].mean(...).values * SE_devisor_solar**-1
        SE02_sol = SE_cap_sol * ds_cf_solar['SE02'] * ds_cf_solar['SE02'].mean(...).values * SE_devisor_solar**-1
        SE03_sol = SE_cap_sol * ds_cf_solar['SE03'] * ds_cf_solar['SE03'].mean(...).values * SE_devisor_solar**-1
        SE04_sol = SE_cap_sol * ds_cf_solar['SE04'] * ds_cf_solar['SE04'].mean(...).values * SE_devisor_solar**-1
        # Sum over all regions to obtain national figures
        ds_ev_solar['SE'] = SE01_sol + SE02_sol + SE03_sol + SE04_sol

        # =============================================================================
        # Working on Norway
        # =============================================================================

        # Determin the devisors
        NO_dev_on = ds_cf_windon['NOM1'].mean(...).values + ds_cf_windon['NON1'].mean(...).values + ds_cf_windon['NOS0'].mean(...).values
        NO_dev_sol = ds_cf_solar['NOM1'].mean(...).values + ds_cf_solar['NON1'].mean(...).values + ds_cf_solar['NOS0'].mean(...).values

        # Find the capacities
        NO_cap_on = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'NO').where(df_cd_tyndp['Type'] == 'Onshore Wind').dropna()['Capacity'].values[0]
        NO_cap_sol = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'NO').where(df_cd_tyndp['Type'] == 'Solar PV').dropna()['Capacity'].values[0]

        # Calculate the regional energy generation for offshore wind
        # If this capacity is not defined, do not calculate
        if df_cd_tyndp.where(df_cd_tyndp['Code'] == 'NO').where(df_cd_tyndp['Type'] == 'Offshore Wind').dropna()['Capacity'].values.size == 0:
            print('There is no offshore capacity for Norway')
        else:
            NO_cap_off = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'NO').where(df_cd_tyndp['Type'] == 'Offshore Wind').dropna()['Capacity'].values[0]
            NO_dev_off = ds_cf_windoff['NOM1_OFF'].mean(...).values + ds_cf_windoff['NON1_OFF'].mean(...).values
            NOM1_off = NO_cap_off * ds_cf_windoff['NOM1_OFF'] * ds_cf_windoff['NOM1_OFF'].mean(...).values * NO_dev_off**-1
            NON1_off = NO_cap_off * ds_cf_windoff['NON1_OFF'] * ds_cf_windoff['NON1_OFF'].mean(...).values * NO_dev_off**-1
            # Sum over all regions to obtain national figures
            ds_ev_windoff['NO'] = NOM1_off + NON1_off

        # Calculate the regional energy generation for onshore wind
        NOM1_on = NO_cap_on * ds_cf_windon['NOM1'] * ds_cf_windon['NOM1'].mean(...).values * NO_dev_on**-1
        NON1_on = NO_cap_on * ds_cf_windon['NON1'] * ds_cf_windon['NON1'].mean(...).values * NO_dev_on**-1
        NOS0_on = NO_cap_on * ds_cf_windon['NOS0'] * ds_cf_windon['NOS0'].mean(...).values * NO_dev_on**-1
        # Sum over all regions to obtain national figures
        ds_ev_windon['NO'] = NOM1_on + NON1_on + NOS0_on

        # Calculate the regional energy generation for solar PV
        NOM1_sol = NO_cap_sol * ds_cf_solar['NOM1'] * ds_cf_solar['NOM1'].mean(...).values * NO_dev_sol**-1
        NON1_sol = NO_cap_sol * ds_cf_solar['NON1'] * ds_cf_solar['NON1'].mean(...).values * NO_dev_sol**-1
        NOS0_sol = NO_cap_sol * ds_cf_solar['NOS0'] * ds_cf_solar['NOS0'].mean(...).values * NO_dev_sol**-1
        # Sum over all regions to obtain national figures
        ds_ev_solar['NO'] = NOM1_sol + NON1_sol + NOS0_sol

        # =============================================================================
        # Working on Denmark
        # =============================================================================

        # Determin the devisors
        DK_dev_off = ds_cf_windoff['DKE1_OFF'].mean(...).values + ds_cf_windoff['DKKF_OFF'].mean(...).values + ds_cf_windoff['DKW1_OFF'].mean(...).values
        DK_dev_on = ds_cf_windon['DKE1'].mean(...).values + ds_cf_windon['DKW1'].mean(...).values
        DK_dev_sol = ds_cf_solar['DKE1'].mean(...).values + ds_cf_solar['DKW1'].mean(...).values

        # Find the capacities
        DK_cap_off = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'DK').where(df_cd_tyndp['Type'] == 'Offshore Wind').dropna()['Capacity'].values[0]
        DK_cap_on = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'DK').where(df_cd_tyndp['Type'] == 'Onshore Wind').dropna()['Capacity'].values[0]
        DK_cap_sol = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'DK').where(df_cd_tyndp['Type'] == 'Solar PV').dropna()['Capacity'].values[0]

        # Calculate the regional energy generation for offshore wind
        DKE1_off = DK_cap_off * ds_cf_windoff['DKE1_OFF'] * ds_cf_windoff['DKE1_OFF'].mean(...).values * DK_dev_off**-1
        DKKF_off = DK_cap_off * ds_cf_windoff['DKKF_OFF'] * ds_cf_windoff['DKKF_OFF'].mean(...).values * DK_dev_off**-1
        DKW1_off = DK_cap_off * ds_cf_windoff['DKW1_OFF'] * ds_cf_windoff['DKW1_OFF'].mean(...).values * DK_dev_off**-1
        # Sum over all regions to obtain national figures
        ds_ev_windoff['DK'] = DKE1_off + DKKF_off + DKW1_off

        # Calculate the regional energy generation for onshore wind
        DKE1_on = DK_cap_on * ds_cf_windon['DKE1'] * ds_cf_windon['DKE1'].mean(...).values * DK_dev_on**-1
        DKW1_on = DK_cap_on * ds_cf_windon['DKW1'] * ds_cf_windon['DKW1'].mean(...).values * DK_dev_on**-1
        # Sum over all regions to obtain national figures
        ds_ev_windon['DK'] = DKE1_on + DKW1_on

        # Calculate the regional energy generation for solar PV
        DKE1_sol = DK_cap_sol * ds_cf_solar['DKE1'] * ds_cf_solar['DKE1'].mean(...).values * DK_dev_sol**-1
        DKW1_sol = DK_cap_sol * ds_cf_solar['DKW1'] * ds_cf_solar['DKW1'].mean(...).values * DK_dev_sol**-1
        # Sum over all regions to obtain national figures
        ds_ev_solar['DK'] = DKE1_sol + DKW1_sol

        # =============================================================================
        # Working on Italy
        # =============================================================================

        # Determin the devisors
        IT_dev_on = ds_cf_windon['ITCA'].mean(...).values + ds_cf_windon['ITCN'].mean(...).values + ds_cf_windon['ITCS'].mean(...).values + ds_cf_windon['ITN1'].mean(...).values + ds_cf_windon['ITS1'].mean(...).values + ds_cf_windon['ITSA'].mean(...).values + ds_cf_windon['ITSI'].mean(...).values
        IT_dev_sol = ds_cf_solar['ITCA'].mean(...).values + ds_cf_solar['ITCN'].mean(...).values + ds_cf_solar['ITCS'].mean(...).values + ds_cf_solar['ITN1'].mean(...).values + ds_cf_solar['ITS1'].mean(...).values + ds_cf_solar['ITSA'].mean(...).values + ds_cf_solar['ITSI'].mean(...).values

        # Find the capacities
        IT_cap_on = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'IT').where(df_cd_tyndp['Type'] == 'Onshore Wind').dropna()['Capacity'].values[0]
        IT_cap_sol = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'IT').where(df_cd_tyndp['Type'] == 'Solar PV').dropna()['Capacity'].values[0]

        # Calculate the regional energy generation for offshore wind
        # If this capacity is not defined, do not calculate
        if df_cd_tyndp.where(df_cd_tyndp['Code'] == 'IT').where(df_cd_tyndp['Type'] == 'Offshore Wind').dropna()['Capacity'].values.size == 0:
            print('There is no offshore capacity for Italy')
        else:

            IT_dev_off = ds_cf_windoff['ITCA_OFF'].mean(...).values + ds_cf_windoff['ITCN_OFF'].mean(...).values + ds_cf_windoff['ITCS_OFF'].mean(...).values + ds_cf_windoff['ITN1_OFF'].mean(...).values + ds_cf_windoff['ITS1_OFF'].mean(...).values + ds_cf_windoff['ITSA_OFF'].mean(...).values + ds_cf_windoff['ITSI_OFF'].mean(...).values
            IT_cap_off = df_cd_tyndp.where(df_cd_tyndp['Code'] == 'IT').where(df_cd_tyndp['Type'] == 'Offshore Wind').dropna()['Capacity'].values[0]
            ITCA_off = IT_cap_off * ds_cf_windoff['ITCA_OFF'] * ds_cf_windoff['ITCA_OFF'].mean(...).values * IT_dev_off**-1
            ITCN_off = IT_cap_off * ds_cf_windoff['ITCN_OFF'] * ds_cf_windoff['ITCN_OFF'].mean(...).values * IT_dev_off**-1
            ITCS_off = IT_cap_off * ds_cf_windoff['ITCS_OFF'] * ds_cf_windoff['ITCS_OFF'].mean(...).values * IT_dev_off**-1
            ITN1_off = IT_cap_off * ds_cf_windoff['ITN1_OFF'] * ds_cf_windoff['ITN1_OFF'].mean(...).values * IT_dev_off**-1
            ITS1_off = IT_cap_off * ds_cf_windoff['ITS1_OFF'] * ds_cf_windoff['ITS1_OFF'].mean(...).values * IT_dev_off**-1
            ITSA_off = IT_cap_off * ds_cf_windoff['ITSA_OFF'] * ds_cf_windoff['ITSA_OFF'].mean(...).values * IT_dev_off**-1
            ITSI_off = IT_cap_off * ds_cf_windoff['ITSI_OFF'] * ds_cf_windoff['ITSI_OFF'].mean(...).values * IT_dev_off**-1
            # Sum over all regions to obtain national figures
            ds_ev_windoff['IT'] = ITCA_off + ITCN_off + ITCS_off + ITN1_off + ITS1_off + ITSA_off + ITSI_off

        # Calculate the regional energy generation for onshore wind
        ITCA_on = IT_cap_on * ds_cf_windon['ITCA'] * ds_cf_windon['ITCA'].mean(...).values * IT_dev_on**-1
        ITCN_on = IT_cap_on * ds_cf_windon['ITCN'] * ds_cf_windon['ITCN'].mean(...).values * IT_dev_on**-1
        ITCS_on = IT_cap_on * ds_cf_windon['ITCS'] * ds_cf_windon['ITCS'].mean(...).values * IT_dev_on**-1
        ITN1_on = IT_cap_on * ds_cf_windon['ITN1'] * ds_cf_windon['ITN1'].mean(...).values * IT_dev_on**-1
        ITS1_on = IT_cap_on * ds_cf_windon['ITS1'] * ds_cf_windon['ITS1'].mean(...).values * IT_dev_on**-1
        ITSA_on = IT_cap_on * ds_cf_windon['ITSA'] * ds_cf_windon['ITSA'].mean(...).values * IT_dev_on**-1
        ITSI_on = IT_cap_on * ds_cf_windon['ITSI'] * ds_cf_windon['ITSI'].mean(...).values * IT_dev_on**-1
        # Sum over all regions to obtain national figures
        ds_ev_windon['IT'] = ITCA_on + ITCN_on + ITCS_on + ITN1_on + ITS1_on + ITSA_on + ITSI_on

        # Calculate the regional energy generation for solar PV
        ITCA_sol = IT_cap_sol * ds_cf_solar['ITCA'] * ds_cf_solar['ITCA'].mean(...).values * IT_dev_sol**-1
        ITCN_sol = IT_cap_sol * ds_cf_solar['ITCN'] * ds_cf_solar['ITCN'].mean(...).values * IT_dev_sol**-1
        ITCS_sol = IT_cap_sol * ds_cf_solar['ITCS'] * ds_cf_solar['ITCS'].mean(...).values * IT_dev_sol**-1
        ITN1_sol = IT_cap_sol * ds_cf_solar['ITN1'] * ds_cf_solar['ITN1'].mean(...).values * IT_dev_sol**-1
        ITS1_sol = IT_cap_sol * ds_cf_solar['ITS1'] * ds_cf_solar['ITS1'].mean(...).values * IT_dev_sol**-1
        ITSA_sol = IT_cap_sol * ds_cf_solar['ITSA'] * ds_cf_solar['ITSA'].mean(...).values * IT_dev_sol**-1
        ITSI_sol = IT_cap_sol * ds_cf_solar['ITSI'] * ds_cf_solar['ITSI'].mean(...).values * IT_dev_sol**-1
        # Sum over all regions to obtain national figures
        ds_ev_solar['IT'] = ITCA_sol + ITCN_sol + ITCS_sol + ITN1_sol + ITS1_sol + ITSA_sol + ITSI_sol

        #%%
        # =============================================================================
        # Time to save the data
        # =============================================================================

        # Setting the general dataset attributes
        ds_ev_windoff.attrs.update(
            author = 'Laurens Stoop UU/KNMI/TenneT',
            variables = 'Wind offshore electricity generation for a certain region',
            units = 'MWh',
            created = datetime.datetime.today().strftime('%d-%m-%Y'),
            region_definition = 'ENTSO-E StudyZones at national level aggregated',
            CapacityDistribution = 'TYNDP-'+capacity_distribution,
            data_source = 'Capacity factors based on ERA5 reanalysis data, contains modified Copernicus Climate Change Service information [31-05-2021]'
            )

        # copy most and update partially
        ds_ev_windon.attrs = ds_ev_windoff.attrs
        ds_ev_windon.attrs.update(
             variables = 'Wind onshore electricity generation for a certain region',
            )

        ds_ev_solar.attrs = ds_ev_windoff.attrs
        ds_ev_solar.attrs.update(
             variables = 'Solar PV electricity generation for a certain region',
            )

        # Saving the files as NetCDF
        ds_ev_windoff.to_netcdf(FOLDER_EV_NUTS0+capacity_distribution+'/ERA5-EU_EV_TYNDP-'+capacity_distribution+'_WOF_'+str(year)+'.nc', encoding={'time':{'units':'days since 1900-01-01'}})
        ds_ev_windon.to_netcdf(FOLDER_EV_NUTS0+capacity_distribution+'/ERA5-EU_EV_TYNDP-'+capacity_distribution+'_WON_'+str(year)+'.nc', encoding={'time':{'units':'days since 1900-01-01'}})
        ds_ev_solar.to_netcdf(FOLDER_EV_NUTS0+capacity_distribution+'/ERA5-EU_EV_TYNDP-'+capacity_distribution+'_SPV_'+str(year)+'.nc', encoding={'time':{'units':'days since 1900-01-01'}})

        # Converting ot Pandas
        df_windoff = ds_ev_windoff.to_pandas()
        df_windon = ds_ev_windon.to_pandas()
        df_solar = ds_ev_solar.to_pandas()

        # Saving as CSV
        df_windoff.to_csv(FOLDER_EV_NUTS0+'csv/ERA5-EU_EV_TYNDP-'+capacity_distribution+'_WOF_'+str(year)+'.csv')
        df_windon.to_csv(FOLDER_EV_NUTS0+'csv/ERA5-EU_EV_TYNDP-'+capacity_distribution+'_WON_'+str(year)+'.csv')
        df_solar.to_csv(FOLDER_EV_NUTS0+'csv/ERA5-EU_EV_TYNDP-'+capacity_distribution+'_SPV_'+str(year)+'.csv')

max_stars_repo_name: laurensstoop/EnergyVariables

max_issues_repo_name: laurensstoop/EnergyVariables

max_forks_repo_name: laurensstoop/EnergyVariables


text:
# Copyright 2018 The dm_control Authors.
#
# 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.
# ============================================================================

"""Tests for `dm_control.mjcf.physics`."""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import copy
import os

# Internal dependencies.
from absl.testing import absltest
from absl.testing import parameterized
from dm_control import mjcf
from dm_control.mjcf import physics as mjcf_physics
from dm_control.mujoco.wrapper import mjbindings
import mock
import numpy as np
import six
from six.moves import cPickle
from six.moves import range
from six.moves import zip

mjlib = mjbindings.mjlib

ARM_MODEL = os.path.join(os.path.dirname(__file__), 'test_assets/robot_arm.xml')

class PhysicsTest(parameterized.TestCase):
  """Tests for `mjcf.Physics`."""

  def setUp(self):
    super(PhysicsTest, self).setUp()
    self.model = mjcf.from_path(ARM_MODEL)
    self.physics = mjcf.Physics.from_xml_string(
        self.model.to_xml_string(), assets=self.model.get_assets())
    self.random = np.random.RandomState(0)

  def sample_elements(self, namespace, single_element):
    all_elements = self.model.find_all(namespace)
    if single_element:
      # A single randomly chosen element from this namespace.
      elements = self.random.choice(all_elements)
      full_identifiers = elements.full_identifier
    else:
      # A random permutation of all elements in this namespace.
      elements = self.random.permutation(all_elements)
      full_identifiers = [element.full_identifier for element in elements]
    return elements, full_identifiers

  def test_construct_and_reload_from_mjcf_model(self):
    physics = mjcf.Physics.from_mjcf_model(self.model)
    physics.data.time = 1.
    physics.reload_from_mjcf_model(self.model)
    self.assertEqual(physics.data.time, 0.)

  @parameterized.parameters(
      # namespace, single_element
      ('geom', True),
      ('geom', False),
      ('joint', True),
      ('joint', False))
  def test_id(self, namespace, single_element):
    elements, full_identifiers = self.sample_elements(namespace, single_element)
    actual = self.physics.bind(elements).element_id
    if single_element:
      expected = self.physics.model.name2id(full_identifiers, namespace)
    else:
      expected = [self.physics.model.name2id(name, namespace)
                  for name in full_identifiers]
    np.testing.assert_array_equal(expected, actual)

  def assertCanGetAndSetBindingArray(
      self, binding, attribute_name, named_indexer, full_identifiers):

    # Read the values using the binding attribute.
    actual = getattr(binding, attribute_name)

    # Read them using the normal named indexing machinery.
    expected = named_indexer[full_identifiers]
    np.testing.assert_array_equal(expected, actual)

    # Assign an array of unique values to the attribute.
    expected = np.arange(actual.size, dtype=actual.dtype).reshape(actual.shape)
    setattr(binding, attribute_name, expected)

    # Read the values back using the normal named indexing machinery.
    actual = named_indexer[full_identifiers]
    np.testing.assert_array_equal(expected, actual)

  @parameterized.parameters(
      # namespace, attribute_name, model_or_data, field_name, single_element
      ('geom', 'xpos', 'data', 'geom_xpos', True),
      ('geom', 'xpos', 'data', 'geom_xpos', False),
      ('joint', 'qpos', 'data', 'qpos', True),
      ('joint', 'qpos', 'data', 'qpos', False),
      ('site', 'rgba', 'model', 'site_rgba', True),
      ('site', 'rgba', 'model', 'site_rgba', False),
      ('sensor', 'sensordata', 'data', 'sensordata', True),
      ('sensor', 'sensordata', 'data', 'sensordata', False))
  def test_attribute_access(self, namespace, attribute_name, model_or_data,
                            field_name, single_element):

    elements, full_identifiers = self.sample_elements(namespace, single_element)
    named_indexer = getattr(getattr(self.physics.named, model_or_data),
                            field_name)
    binding = self.physics.bind(elements)
    self.assertCanGetAndSetBindingArray(
        binding, attribute_name, named_indexer, full_identifiers)

  @parameterized.parameters(
      # namespace, attribute_name, model_or_data, field_name, single_element,
      # column_index
      ('geom', 'pos', 'model', 'geom_pos', True, None),
      ('geom', 'pos', 'model', 'geom_pos', False, None),
      ('geom', 'pos', 'model', 'geom_pos', True, 1),
      ('geom', 'pos', 'model', 'geom_pos', False, 1),
      ('geom', 'pos', 'model', 'geom_pos', True, 'y'),
      ('geom', 'pos', 'model', 'geom_pos', False, 'y'),
      ('geom', 'pos', 'model', 'geom_pos', True, slice(0, None, 2)),
      ('geom', 'pos', 'model', 'geom_pos', False, slice(0, None, 2)),
      ('geom', 'pos', 'model', 'geom_pos', True, [0, 2]),
      ('geom', 'pos', 'model', 'geom_pos', False, [0, 2]),
      ('geom', 'pos', 'model', 'geom_pos', True, ['x', 'z']),
      ('geom', 'pos', 'model', 'geom_pos', False, ['x', 'z']),
      ('joint', 'qpos', 'data', 'qpos', True, None),
      ('joint', 'qpos', 'data', 'qpos', False, None))
  def test_indexing(self, namespace, attribute_name, model_or_data,
                    field_name, single_element, column_index):

    elements, full_identifiers = self.sample_elements(namespace, single_element)
    named_indexer = getattr(getattr(self.physics.named, model_or_data),
                            field_name)
    binding = self.physics.bind(elements)

    if column_index is not None:
      binding_index = (attribute_name, column_index)
      try:
        named_index = np.ix_(full_identifiers, column_index)
      except ValueError:
        named_index = (full_identifiers, column_index)
    else:
      binding_index = attribute_name
      named_index = full_identifiers

    # Read the values by indexing the binding.
    actual = binding[binding_index]

    # Read them using the normal named indexing machinery.
    expected = named_indexer[named_index]
    np.testing.assert_array_equal(expected, actual)

    # Write an array of unique values into the binding.
    expected = np.arange(actual.size, dtype=actual.dtype).reshape(actual.shape)
    binding[binding_index] = expected

    # Read the values back using the normal named indexing machinery.
    actual = named_indexer[named_index]
    np.testing.assert_array_equal(expected, actual)

  def test_bind_mocap_body(self):
    pos = [1, 2, 3]
    quat = [1, 0, 0, 0]
    model = mjcf.RootElement()
    # Bodies are non-mocap by default.
    normal_body = model.worldbody.add('body', pos=pos, quat=quat)
    mocap_body = model.worldbody.add('body', pos=pos, quat=quat, mocap='true')
    physics = mjcf.Physics.from_xml_string(model.to_xml_string())

    binding = physics.bind(mocap_body)
    np.testing.assert_array_equal(pos, binding.mocap_pos)
    np.testing.assert_array_equal(quat, binding.mocap_quat)

    new_pos = [4, 5, 6]
    new_quat = [0, 1, 0, 0]
    binding.mocap_pos = new_pos
    binding.mocap_quat = new_quat
    np.testing.assert_array_equal(
        new_pos, physics.named.data.mocap_pos[mocap_body.full_identifier])
    np.testing.assert_array_equal(
        new_quat, physics.named.data.mocap_quat[mocap_body.full_identifier])

    with self.assertRaises(AttributeError):
      _ = physics.bind(normal_body).mocap_pos

    with six.assertRaisesRegex(
        self,
        ValueError,
        'Cannot bind to a collection containing multiple element types'):
      physics.bind([mocap_body, normal_body])

  def test_bind_worldbody(self):
    expected_mass = 10
    model = mjcf.RootElement()
    child = model.worldbody.add('body')
    child.add('geom', type='sphere', size=[0.1], mass=expected_mass)
    physics = mjcf.Physics.from_mjcf_model(model)
    mass = physics.bind(model.worldbody).subtreemass
    self.assertEqual(mass, expected_mass)

  def test_caching(self):
    all_joints = self.model.find_all('joint')

    original = self.physics.bind(all_joints)
    cached = self.physics.bind(all_joints)
    self.assertIs(cached, original)

    different_order = self.physics.bind(all_joints[::-1])
    self.assertIsNot(different_order, original)

    # Reloading the `Physics` instance should clear the cache.
    self.physics.reload_from_xml_string(
        self.model.to_xml_string(), assets=self.model.get_assets())
    after_reload = self.physics.bind(all_joints)
    self.assertIsNot(after_reload, original)

  def test_exceptions(self):
    joint = self.model.find_all('joint')[0]
    geom = self.model.find_all('geom')[0]
    with six.assertRaisesRegex(
        self,
        ValueError,
        'Cannot bind to a collection containing multiple element types'):
      self.physics.bind([joint, geom])

    with six.assertRaisesRegex(self, ValueError, 'cannot be bound to physics'):
      mjcf.physics.Binding(self.physics, 'invalid_namespace', 'whatever')

    binding = self.physics.bind(joint)
    with six.assertRaisesRegex(self, AttributeError, 'does not have attribute'):
      getattr(binding, 'invalid_attribute')

  def test_dirty(self):
    self.physics.forward()
    self.assertFalse(self.physics.is_dirty)

    joints, _ = self.sample_elements('joint', single_element=False)
    sites, _ = self.sample_elements('site', single_element=False)

    # Accessing qpos shouldn't trigger a recalculation.
    _ = self.physics.bind(joints).qpos
    self.assertFalse(self.physics.is_dirty)

    # Reassignments to qpos should cause the physics to become dirty.
    site_xpos_before = copy.deepcopy(self.physics.bind(sites).xpos)
    self.physics.bind(joints).qpos += 0.5
    self.assertTrue(self.physics.is_dirty)

    # Accessing stuff in mjModel shouldn't trigger a recalculation.
    _ = self.physics.bind(sites).pos
    self.assertTrue(self.physics.is_dirty)

    # Accessing stuff in mjData should trigger a recalculation.
    actual_sites_xpos_after = copy.deepcopy(self.physics.bind(sites).xpos)
    self.assertFalse(self.physics.is_dirty)
    self.assertFalse((actual_sites_xpos_after == site_xpos_before).all())

    # Automatic recalculation should render `forward` a no-op here.
    self.physics.forward()
    expected_sites_xpos_after = self.physics.bind(sites).xpos
    np.testing.assert_array_equal(actual_sites_xpos_after,
                                  expected_sites_xpos_after)

    # `forward` should not be called on subsequent queries to xpos.
    with mock.patch.object(
        self.physics, 'forward',
        side_effect=self.physics.forward) as mock_forward:
      _ = self.physics.bind(sites).xpos
      mock_forward.assert_not_called()

  @parameterized.parameters(True, False)
  def test_assign_while_dirty(self, assign_via_slice):
    actuators = self.model.find_all('actuator')
    if assign_via_slice:
      self.physics.bind(actuators).ctrl[:] = 0.75
    else:
      self.physics.bind(actuators).ctrl = 0.75
    self.assertTrue(self.physics.is_dirty)
    self.physics.step()
    self.assertTrue(self.physics.is_dirty)
    sensors = self.model.find_all('sensor')
    if assign_via_slice:
      self.physics.bind(sensors).sensordata[:] = 12345
    else:
      self.physics.bind(sensors).sensordata = 12345
    self.assertFalse(self.physics.is_dirty)
    np.testing.assert_array_equal(
        self.physics.bind(sensors).sensordata,
        [12345] * len(self.physics.bind(sensors).sensordata))

  def test_setitem_on_binding_attr(self):
    bodies, _ = self.sample_elements('body', single_element=False)
    xfrc_binding = self.physics.bind(bodies).xfrc_applied

    xfrc_binding[:, 1] = list(range(len(bodies)))
    for i, body in enumerate(bodies):
      self.assertEqual(xfrc_binding[i, 1], i)
      self.assertEqual(
          self.physics.named.data.xfrc_applied[body.full_identifier][1], i)

    xfrc_binding[:, 1] *= 2
    for i, body in enumerate(bodies):
      self.assertEqual(xfrc_binding[i, 1], 2 * i)
      self.assertEqual(
          self.physics.named.data.xfrc_applied[body.full_identifier][1], 2 * i)

    xfrc_binding[[1, 3, 5], 2] = 42
    for i, body in enumerate(bodies):
      actual_value = (
          self.physics.named.data.xfrc_applied[body.full_identifier][2])
      if i in [1, 3, 5]:
        self.assertEqual(actual_value, 42)
      else:
        self.assertNotEqual(actual_value, 42)

    # Bind to a single element.
    single_binding = self.physics.bind(bodies[0]).xfrc_applied
    single_binding[:2] = 55
    np.testing.assert_array_equal(single_binding[:2], [55, 55])
    np.testing.assert_array_equal(
        self.physics.named.data.xfrc_applied[bodies[0].full_identifier][:2],
        [55, 55])

  def test_empty_binding(self):
    binding = self.physics.bind([])
    self.assertEqual(binding.xpos.shape, (0,))
    with self.assertRaisesWithLiteralMatch(
        ValueError, 'Cannot assign a value to an empty binding.'):
      binding.xpos = 5

  @parameterized.parameters([('data', 'act'), ('data', 'act_dot')])
  def test_actuator_state_binding(self, model_or_data, attribute_name):

    def make_model_with_mixed_actuators():
      actuators = []
      is_stateful = []
      root = mjcf.RootElement()
      body = root.worldbody.add('body')
      body.add('geom', type='sphere', size=[0.1])
      slider = body.add('joint', type='slide', name='slide_joint')
      # Third-order `general` actuator.
      actuators.append(
          root.actuator.add(
              'general', dyntype='integrator', biastype='affine',
              dynprm=[1, 0, 0], joint=slider, name='general_act'))
      is_stateful.append(True)
      # Cylinder actuators are also third-order.
      actuators.append(
          root.actuator.add('cylinder', joint=slider, name='cylinder_act'))
      is_stateful.append(True)
      # A second-order actuator, added after the third-order actuators. The
      # actuators will be automatically reordered in the generated XML so that
      # the second-order actuator comes first.
      actuators.append(
          root.actuator.add('velocity', joint=slider, name='velocity_act'))
      is_stateful.append(False)
      return root, actuators, is_stateful

    model, actuators, is_stateful = make_model_with_mixed_actuators()
    physics = mjcf.Physics.from_mjcf_model(model)
    binding = physics.bind(actuators)
    named_indexer = getattr(
        getattr(physics.named, model_or_data), attribute_name)
    stateful_actuator_names = [
        actuator.full_identifier
        for actuator, stateful in zip(actuators, is_stateful) if stateful]
    self.assertCanGetAndSetBindingArray(
        binding, attribute_name, named_indexer, stateful_actuator_names)

  def test_bind_stateless_actuators_only(self):
    actuators = []
    root = mjcf.RootElement()
    body = root.worldbody.add('body')
    body.add('geom', type='sphere', size=[0.1])
    slider = body.add('joint', type='slide', name='slide_joint')
    actuators.append(
        root.actuator.add('velocity', joint=slider, name='velocity_act'))
    actuators.append(
        root.actuator.add('motor', joint=slider, name='motor_act'))
    # `act` should be an empty array if there are no stateful actuators.
    physics = mjcf.Physics.from_mjcf_model(root)
    self.assertEqual(physics.bind(actuators).act.shape, (0,))

  def make_simple_model(self):
    def add_submodel(root):
      body = root.worldbody.add('body')
      geom = body.add('geom', type='ellipsoid', size=[0.1, 0.2, 0.3])
      site = body.add('site', type='sphere', size=[0.1])
      return body, geom, site
    root = mjcf.RootElement()
    add_submodel(root)
    add_submodel(root)
    return root

  def quat2mat(self, quat):
    result = np.empty(9, dtype=np.double)
    mjlib.mju_quat2Mat(result, np.asarray(quat))
    return result

  @parameterized.parameters(['body', 'geom', 'site'])
  def test_write_to_pos(self, entity_type):
    root = self.make_simple_model()
    entity1, entity2 = root.find_all(entity_type)
    physics = mjcf.Physics.from_mjcf_model(root)
    first = physics.bind(entity1)
    second = physics.bind(entity2)

    # Initially both entities should be 'sameframe'
    self.assertEqual(first.sameframe, 1)
    self.assertEqual(second.sameframe, 1)

    # Assigning to `pos` should disable 'sameframe' only for that entity.
    new_pos = (0., 0., 0.1)
    first.pos = new_pos
    self.assertEqual(first.sameframe, 0)
    self.assertEqual(second.sameframe, 1)
    # `xpos` should reflect the new position.
    np.testing.assert_array_equal(first.xpos, new_pos)

    # Writing into the `pos` array should also disable 'sameframe'.
    new_x = -0.1
    pos_array = second.pos
    pos_array[0] = new_x
    self.assertEqual(second.sameframe, 0)
    # `xpos` should reflect the new position.
    self.assertEqual(second.xpos[0], new_x)

  @parameterized.parameters(['body', 'geom', 'site'])
  def test_write_to_quat(self, entity_type):
    root = self.make_simple_model()
    entity1, entity2 = root.find_all(entity_type)
    physics = mjcf.Physics.from_mjcf_model(root)
    first = physics.bind(entity1)
    second = physics.bind(entity2)

    # Initially both entities should be 'sameframe'
    self.assertEqual(first.sameframe, 1)
    self.assertEqual(second.sameframe, 1)

    # Assigning to `quat` should disable 'sameframe' only for that entity.
    new_quat = (0., 0., 0., 1.)
    first.quat = new_quat
    self.assertEqual(first.sameframe, 0)
    self.assertEqual(second.sameframe, 1)
    # `xmat` should reflect the new quaternion.
    np.testing.assert_allclose(first.xmat, self.quat2mat(new_quat))

    # Writing into the `quat` array should also disable 'sameframe'.
    new_w = -1.
    quat_array = second.quat
    quat_array[0] = new_w
    self.assertEqual(second.sameframe, 0)
    # `xmat` should reflect the new quaternion.
    np.testing.assert_allclose(second.xmat, self.quat2mat(quat_array))

  def test_write_to_ipos(self):
    root = self.make_simple_model()
    entity1, entity2 = root.find_all('body')
    physics = mjcf.Physics.from_mjcf_model(root)
    first = physics.bind(entity1)
    second = physics.bind(entity2)

    # Initially both bodies should be 'simple' and 'sameframe'
    self.assertEqual(first.simple, 1)
    self.assertEqual(first.sameframe, 1)
    self.assertEqual(second.simple, 1)
    self.assertEqual(second.sameframe, 1)

    # Assigning to `ipos` should disable 'simple' and 'sameframe' only for that
    # body.
    new_ipos = (0., 0., 0.1)
    first.ipos = new_ipos
    self.assertEqual(first.simple, 0)
    self.assertEqual(first.sameframe, 0)
    self.assertEqual(second.simple, 1)
    self.assertEqual(second.sameframe, 1)
    # `xipos` should reflect the new position.
    np.testing.assert_array_equal(first.xipos, new_ipos)

    # Writing into the `ipos` array should also disable 'simple' and
    # 'sameframe'.
    new_x = -0.1
    ipos_array = second.ipos
    ipos_array[0] = new_x
    self.assertEqual(second.simple, 0)
    self.assertEqual(second.sameframe, 0)
    # `xipos` should reflect the new position.
    self.assertEqual(second.xipos[0], new_x)

  def test_write_to_iquat(self):
    root = self.make_simple_model()
    entity1, entity2 = root.find_all('body')
    physics = mjcf.Physics.from_mjcf_model(root)
    first = physics.bind(entity1)
    second = physics.bind(entity2)

    # Initially both bodies should be 'simple' and 'sameframe'
    self.assertEqual(first.simple, 1)
    self.assertEqual(first.sameframe, 1)
    self.assertEqual(second.simple, 1)
    self.assertEqual(second.sameframe, 1)

    # Assigning to `iquat` should disable 'simple' and 'sameframe' only for that
    # body.
    new_iquat = (0., 0., 0., 1.)
    first.iquat = new_iquat
    self.assertEqual(first.simple, 0)
    self.assertEqual(first.sameframe, 0)
    self.assertEqual(second.simple, 1)
    self.assertEqual(second.sameframe, 1)
    # `ximat` should reflect the new quaternion.
    np.testing.assert_allclose(first.ximat, self.quat2mat(new_iquat))

    # Writing into the `iquat` array should also disable 'simple' and
    # 'sameframe'.
    new_w = -0.1
    iquat_array = second.iquat
    iquat_array[0] = new_w
    self.assertEqual(second.simple, 0)
    self.assertEqual(second.sameframe, 0)
    # `ximat` should reflect the new quaternion.
    np.testing.assert_allclose(second.ximat, self.quat2mat(iquat_array))

  @parameterized.parameters([dict(order='C'), dict(order='F')])
  def test_copy_synchronizing_array_wrapper(self, order):
    root = self.make_simple_model()
    physics = mjcf.Physics.from_mjcf_model(root)
    xpos_view = physics.bind(root.find_all('body')).xpos
    xpos_view_copy = xpos_view.copy(order=order)

    np.testing.assert_array_equal(xpos_view, xpos_view_copy)
    self.assertFalse(np.may_share_memory(xpos_view, xpos_view_copy),
                     msg='Original and copy should not share memory.')
    self.assertIs(type(xpos_view_copy), np.ndarray)

    # Check that `order=` is respected.
    if order == 'C':
      self.assertTrue(xpos_view_copy.flags.c_contiguous)
      self.assertFalse(xpos_view_copy.flags.f_contiguous)
    elif order == 'F':
      self.assertFalse(xpos_view_copy.flags.c_contiguous)
      self.assertTrue(xpos_view_copy.flags.f_contiguous)

    # The copy should be writeable.
    self.assertTrue(xpos_view_copy.flags.writeable)
    new_value = 99.
    xpos_view_copy[0, -1] = new_value
    self.assertEqual(xpos_view_copy[0, -1], new_value)

  def test_error_when_pickling_synchronizing_array_wrapper(self):
    root = self.make_simple_model()
    physics = mjcf.Physics.from_mjcf_model(root)
    xpos_view = physics.bind(root.find_all('body')).xpos
    with self.assertRaisesWithLiteralMatch(
        NotImplementedError,
        mjcf_physics._PICKLING_NOT_SUPPORTED.format(type=type(xpos_view))):
      cPickle.dumps(xpos_view)

if __name__ == '__main__':
  absltest.main()

ext: py

reason: import numpy


text:
import os

import numpy as np
import h5py

import thimbles as tmb

resources_dir = os.path.join(os.path.dirname(tmb.__file__), "resources")
standards_file = os.path.join(resources_dir, "standards.h5")

class StandardsDatabase(object):

    def __init__(self, standards_file=standards_file):
        pass

    def extract_standard(self):
        pass

    def inject_standard(self):
        pass

    def match_standard(self, spectra):
        wv_bounds = [(s.wv[0], s.wv[-1]) for s in spectra]
        min_wv = np.min(wv_bounds)
        max_wv = np.max(wv_bounds)

max_stars_repo_name: quidditymaster/thimbles

max_issues_repo_name: quidditymaster/thimbles

max_forks_repo_name: quidditymaster/thimbles


text:
import tensorflow as tf
import zutils.tf_math_funcs as tmf
import numpy as np
import zutils.pt_utils as ptu
import prettytensor as pt
from net_modules.common import *

from net_modules.hourglass import hourglass

from net_modules.auto_struct.keypoint_decoder import Factory as BaseFactory

class Factory(BaseFactory):

    @staticmethod
    def pt_defaults_scope_value():
        return {
            # non-linearity
            'activation_fn': default_activation.current_value,
            # for batch normalization
            'batch_normalize': True,
            'learned_moments_update_rate': 0.0003,
            # 'learned_moments_update_rate': 1.,
            'variance_epsilon': 0.001,
            'scale_after_normalization': True
        }

    default_patch_feature_dim = 8

    def __init__(self, recon_dist_param_num=1, options=None):
        super().__init__(recon_dist_param_num, options)
        if "image_channels" in options:
            self.image_channels = options["image_channels"]
        else:
            self.image_channels = 3

    def image_size(self):
        return 64, 64

    def feature2image(self, feature_map):
        output_channels = 3*self.recon_dist_param_num
        hgd = [
            {"type": "conv2d", "depth": 64},
            {"type": "skip", "layer_num": 2},
            {"type": "pool", "pool": "max", "kernel": 2, "stride": 2},  # 32x32
            {"type": "conv2d", "depth": 128},
            {"type": "skip", "layer_num": 2},
            {"type": "pool", "pool": "max", "kernel": 2, "stride": 2},  # 16x16
            {"type": "conv2d", "depth": 256},
            {"type": "skip", "layer_num": 2},
            {"type": "pool", "pool": "max", "kernel": 2, "stride": 2},  # 8x8
            {"type": "conv2d", "depth": 512},
            {"type": "skip", "layer_num": 2},
            {"type": "pool", "pool": "max", "kernel": 2, "stride": 2},  # 4x4
            {"type": "conv2d", "depth": 512},
        ]
        with pt.defaults_scope(**self.pt_defaults_scope_value()):
            output_tensor = hourglass(
                feature_map, hgd,
                net_type=self.options["hourglass_type"] if "hourglass_type" in self.options else None,
                extra_highlevel_feature=None
            )
            output_tensor = (
                pt.wrap(output_tensor).
                deconv2d(3, 32, stride=2).
                conv2d(3, output_channels, activation_fn=None)
            ).tensor
        return output_tensor

    def input_feature_dim(self):
        return 64

    def bg_feature2image(self, bg_feature):
        batch_size = tmf.get_shape(bg_feature)[0]
        with pt.defaults_scope(**self.pt_defaults_scope_value()):
            return (
                pt.wrap(bg_feature).
                conv2d(3, 512).                 # 2
                deconv2d(3, 512, stride=2).     # 4
                deconv2d(3, 256, stride=2).     # 8
                deconv2d(3, 256, stride=2).     # 16
                deconv2d(3, 128, stride=2).     # 32
                deconv2d(3, 64, stride=2).      # 64
                deconv2d(3, 32, stride=2).      # 128
                conv2d(3, 3*self.recon_dist_param_num, activation_fn=None).tensor
            )


text:
"""
Convert each frame in a anime video to 256 x 256 images
"""
import cv2
import os
import numpy as np
import argparse
from moviepy.video.io.VideoFileClip import VideoFileClip
from tqdm import tqdm

def parse_args():
    parser = argparse.ArgumentParser()
    parser.add_argument('--video-path', type=str, default='/home/ubuntu/Downloads/kimetsu_yaiba.mp4')
    parser.add_argument('--save-path', type=str, default='./script/test_crop')
    parser.add_argument('--max-image', type=int, default=1800)
    parser.add_argument('--start', type=int, default=0)
    parser.add_argument('--end', type=int, default=0)
    parser.add_argument('--image-size', type=int, default=256)

    return parser.parse_args()

class VideoConverter:
    def __init__(self, video_path, save_dir, max_image=1600, start=0, end=0, image_size=256):
        os.makedirs(save_dir, exist_ok=True)

        assert os.path.isfile(video_path), f'{video_path} must be a video file'

        self.save_dir = save_dir
        self.video_path = video_path
        self.max_image = max_image
        self.start = start
        self.end = end or None
        self.counter = 0
        self.image_size = image_size

    def crop_and_save(self, image):
        '''
        Crop a large image into smaller images of given size

        @Returns:
            - True if counter reach max_image
        '''

        height, width, _ = image.shape

        # Image can be divided into rows * cols sub-images
        rows = height // self.image_size
        cols = width // self.image_size

        for r in range(rows):
            for c in range(cols):
                start_x = r * self.image_size
                end_x = start_x + self.image_size
                start_y = c * self.image_size
                end_y = start_y + self.image_size
                # print(f'x: [{start_x}:{end_x}], y: [{start_y}:{end_y}]')
                sub_im = image[start_x: end_x, start_y: end_y, :]

                if np.std(sub_im) > 25.0:
                    save_path = os.path.join(self.save_dir, f'{self.counter}.jpg')
                    cv2.imwrite(save_path, sub_im)
                    self.counter += 1

                if self.counter == self.max_image:
                    return True

        return False

    def process(self):
        '''
        Process video
        '''
        video_clip = VideoFileClip(self.video_path)
        if self.start or self.end:
            video_clip = video_clip.subclip(self.start, self.end)

        video_clip = video_clip.set_fps(video_clip.fps // 10)

        total_frames = round(video_clip.fps * video_clip.duration)
        print(f'Processing video {self.video_path}, {total_frames} frames, size: {video_clip.size}')

        for frame in tqdm(video_clip.iter_frames()):
            # It's better if we resizing before crop to keep the image looks more 'scene'
            h, w, _ = frame.shape
            aspect_ratio = w / h
            ratio = h / (self.image_size * 2)
            w /= ratio
            h = w / aspect_ratio

            frame = cv2.resize(frame[...,::-1], (int(h) , int(w)))
            if self.crop_and_save(frame):
                break

        print(f'Saved {self.counter} images to {self.save_dir}')

if __name__ == '__main__':

    args = parse_args()
    converter = VideoConverter(
        args.video_path,
        args.save_path,
        max_image=args.max_image,
        start=args.start,
        end=args.end,
        image_size=args.image_size
        )

    converter.process()


text:
# -*- coding: utf-8 -*-

# MIT License
#
# Copyright (c) 2017 Moritz Lode
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

import os.path
import sys
import codecs
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

from youtubeAPICrawler.database import *

class ThinDiGraph(nx.DiGraph):
    all_edge_dict = {'weight': 1}
    def single_edge_dict(self):
        return self.all_edge_dict
    edge_attr_dict_factory = single_edge_dict

if os.path.isfile('yt.adjlist'):
    G=nx.read_adjlist('yt.adjlist', create_using=ThinDiGraph())
else:
    db = YTDatabase() # db should be configured to the 16k ytdatabase
    G=ThinDiGraph()

    with db._session_scope(True) as session:

        print '#', session.query(Channel).count() # Number of channels in db

        channels = session.query(Channel).all()

        for channel in channels:
            G.add_node(channel.id)

            for link in channel.featured:
                G.add_edge(link.channelID, link.featuredChannelID)

    nx.write_adjlist(G,"yt.adjlist")

print G.number_of_nodes()
print G.size()

pr = nx.pagerank_scipy(G)

import operator

sorted_pr = sorted(pr.items(), key=operator.itemgetter(1), reverse=True)

sorted_deg = sorted( [(n, G.degree(n)) for n in G.nodes()], key=operator.itemgetter(1), reverse=True)

sorted_degpr = sorted( [(n, G.degree(n)*pr[n]) for n in G.nodes()], key=operator.itemgetter(1), reverse=True)

for ch in sorted_pr[:5]:
    print ch[0], 'pr:', ch[1], 'deg:', G.degree(ch[0]), 'out:', len(G.out_edges(ch[0])), 'in:', len(G.in_edges(ch[0]))

print '...'
for ch in sorted_pr[len(sorted_pr)-5:len(sorted_pr)]:
    print ch[0], 'pr:', pr[ch[0]], 'deg:', G.degree(ch[0]), 'out:', len(G.out_edges(ch[0])), 'in:', len(G.in_edges(ch[0]))

print '\n'

for ch in sorted_deg[:5]:
    print ch[0], 'deg:', ch[1], 'out:', len(G.out_edges(ch[0])), 'in:', len(G.in_edges(ch[0]))

print '...'
for ch in sorted_deg[len(sorted_deg)-5:len(sorted_deg)]:
    print ch[0], 'deg:', ch[1], 'out:', len(G.out_edges(ch[0])), 'in:', len(G.in_edges(ch[0]))

print '\n'

for ch in sorted_degpr[:5]:
    print ch[0], 'deg*pr:', ch[1], 'out:', len(G.out_edges(ch[0])), 'in:', len(G.in_edges(ch[0]))

print '...'
for ch in sorted_degpr[len(sorted_degpr)-5:len(sorted_degpr)]:
    print ch[0], 'deg*pr:', ch[1], 'out:', len(G.out_edges(ch[0])), 'in:', len(G.in_edges(ch[0]))

print '\n'

#with open('page_rank_sampled_channel_id.json', 'wb') as channelid:
#    channelid.write(json.dumps([id[0] for id in sorted_pr[:800]]))

largest = max(nx.strongly_connected_components(G), key=len)
print 'largest strong connected subgraph:', len(largest), '\n'

# trying to get all bidirectional connected nodes only
print '############ BI-DIR ###########\n'

bidirectnodes = []
BG = ThinDiGraph()

for nodeu in G.nodes_iter():
    for nodep in G.neighbors(nodeu):
        if G.has_edge(nodeu, nodep) and G.has_edge(nodep, nodeu):
            BG.add_node(nodeu)
            BG.add_node(nodep)
            BG.add_edge(nodeu, nodep)
            BG.add_edge(nodep, nodeu)

print BG.number_of_nodes()
print BG.size()

UG = BG.to_undirected(True) # is the same as BG only undirected
print UG.number_of_nodes()
print UG.size()

with open('bidi_channel_graph_ids.json', 'wb') as sfile:
    sfile.write(json.dumps(UG.nodes()))

nx.write_adjlist(BG,"yt_bidi.adjlist")

prbi = nx.pagerank_scipy(BG)

import operator
sorted_prbi = sorted(prbi.items(), key=operator.itemgetter(1), reverse=True)

sorted_deg = sorted( [(n, BG.degree(n)) for n in BG.nodes()],key=operator.itemgetter(1), reverse=True)

print sorted_prbi[:5]

print '\n'

print sorted_prbi[len(sorted_prbi)-5:len(sorted_prbi)]

for ch in sorted_prbi[:5]:
    print ch[0], 'pr:', ch[1], 'deg:', BG.degree(ch[0]), 'out:', len(BG.out_edges(ch[0])), 'in:', len(BG.in_edges(ch[0]))

for ch in sorted_prbi[len(sorted_prbi)-5:len(sorted_prbi)]:
    print ch[0], 'pr:', prbi[ch[0]], 'deg:', BG.degree(ch[0]), 'out:', len(BG.out_edges(ch[0])), 'in:', len(BG.in_edges(ch[0]))

strongsubgraphs = [Gc for Gc in sorted(nx.strongly_connected_component_subgraphs(BG),key=len, reverse=True)]

print 'biggest strongly connected subgraphs:'
for cs in strongsubgraphs[:5]:
    print 'nodes:', cs.number_of_nodes(),'edges:', cs.size()

#app = Viewer(SBGnn)
#app.mainloop()

#pos=nx.spring_layout(SBG)
#nx.draw(SBG,pos, edge_cmap=plt.cm.Reds)
#plt.show()

max_stars_repo_name: alxgrh/CATANA

max_issues_repo_name: christiannkoch/CATANA

max_forks_repo_name: lodemo/CATANA

reason: import numpy,import networkx


text:
# pylint: disable=too-many-statements

"""Joints Analysis: Supplementary material of MonStereo"""

import json
import os
from collections import defaultdict

import numpy as np
import matplotlib.pyplot as plt

from ..utils import find_cluster, average
from ..visuals.figures import get_distances
from ..prep.transforms import COCO_KEYPOINTS

def joints_variance(joints, clusters, dic_ms):
    # CLUSTERS = ('3', '5', '7', '9', '11', '13', '15', '17', '19', '21', '23', '25', '27', '29', '31', '49')
    BF = 0.54 * 721
    phase = 'train'
    methods = ('pifpaf', 'mask')
    dic_fin = {}

    for method in methods:
        dic_var = defaultdict(lambda: defaultdict(list))
        dic_joints = defaultdict(list)
        dic_avg = defaultdict(lambda: defaultdict(float))
        path_joints = joints + '_' + method + '.json'

        with open(path_joints, 'r') as f:
            dic_jo = json.load(f)

        for idx, keypoint in enumerate(dic_jo[phase]['kps']):
            # if dic_jo[phase]['names'][idx] == '005856.txt' and dic_jo[phase]['Y'][idx][2] > 14:
            #     aa = 4
            assert len(keypoint) < 2
            kps = np.array(keypoint[0])[:, :17]
            kps_r = np.array(keypoint[0])[:, 17:]
            disps = kps[0] - kps_r[0]
            zz = dic_jo[phase]['Y'][idx][2]
            disps_3 = get_variance(kps, kps_r, zz)
            disps_8 = get_variance_conf(kps, kps_r, num=8)
            disps_4 = get_variance_conf(kps, kps_r, num=4)
            disp_gt = BF / zz
            clst = find_cluster(zz, clusters)   # 4 = '3'    35 = '31'     42 = 2 = 'excl'
            dic_var['std_d'][clst].append(disps.std())
            errors = np.minimum(30, np.abs(zz - BF / disps))
            dic_var['mean_dev'][clst].append(min(30, abs(zz - BF / np.median(disps))))
            dic_var['mean_3'][clst].append(min(30, abs(zz - BF / disps_3.mean())))
            dic_var['mean_8'][clst].append(min(30, abs(zz - BF / np.median(disps_8))))
            dic_var['mean_4'][clst].append(min(30, abs(zz - BF / np.median(disps_4))))
            arg_best = np.argmin(errors)
            conf = np.mean((kps[2][arg_best], kps_r[2][arg_best]))
            dic_var['mean_best'][clst].append(np.min(errors))
            dic_var['conf_best'][clst].append(conf)
            dic_var['conf'][clst].append(np.mean((np.mean(kps[2]), np.mean(kps_r[2]))))
            # dic_var['std_z'][clst].append(zzs.std())
            for ii, el in enumerate(disps):
                if abs(el-disp_gt) < 1:
                    dic_var['rep'][clst].append(1)
                    dic_joints[str(ii)].append(1)
                else:
                    dic_var['rep'][clst].append(0)
                    dic_joints[str(ii)].append(0)

        for key in dic_var:
            for clst in clusters[:-1]:  # 41 needs to be excluded (36 = '31')
                dic_avg[key][clst] = average(dic_var[key][clst])
        dic_fin[method] = dic_avg
        for key in dic_joints:
            dic_fin[method]['joints'][key] = average(dic_joints[key])
        dic_fin['monstereo'] = {clst: dic_ms[clst]['mean'] for clst in clusters[:-1]}
    variance_figures(dic_fin, clusters)

def get_variance(kps, kps_r, zz):

    thresh = 0.5 - zz / 100
    disps_2 = []
    disps = kps[0] - kps_r[0]
    arg_disp = np.argsort(disps)[::-1]

    for idx in arg_disp[1:]:
        if kps[2][idx] > thresh and kps_r[2][idx] > thresh:
            disps_2.append(disps[idx])
        if len(disps_2) >= 3:
            return np.array(disps_2)
    return disps

def get_variance_conf(kps, kps_r, num=8):

    disps_conf = []
    confs = (kps[2, :] + kps_r[2, :]) / 2
    disps = kps[0] - kps_r[0]
    arg_disp = np.argsort(confs)[::-1]

    for idx in arg_disp[:num]:
        disps_conf.append(disps[idx])
    return np.array(disps_conf)

def variance_figures(dic_fin, clusters):
    """Predicted confidence intervals and task error as a function of ground-truth distance"""

    dir_out = 'docs'
    x_min = 3
    x_max = 43
    y_min = 0
    y_max = 1

    plt.figure(0)
    plt.xlabel("Ground-truth distance [m]")
    plt.title("Repeatability by distance")
    plt.xlim(x_min, x_max)
    plt.ylim(y_min, y_max)
    plt.grid(linewidth=0.2)

    xxs = get_distances(clusters)
    yys_p = [el for _, el in dic_fin['pifpaf']['rep'].items()]
    yys_m = [el for _, el in dic_fin['mask']['rep'].items()]
    plt.plot(xxs, yys_p, marker='s', label="PifPaf")
    plt.plot(xxs, yys_m, marker='o', label="Mask R-CNN")
    plt.tight_layout()
    plt.legend()
    path_fig = os.path.join(dir_out, 'repeatability.png')
    plt.savefig(path_fig)
    print("Figure of repeatability saved in {}".format(path_fig))

    plt.figure(1)
    plt.xlabel("Ground-truth distance [m]")
    plt.ylabel("[m]")
    plt.title("Depth error")
    plt.grid(linewidth=0.2)
    y_min = 0
    y_max = 2.7
    plt.ylim(y_min, y_max)
    yys_p = [el for _, el in dic_fin['pifpaf']['mean_dev'].items()]
    # yys_m = [el for _, el in dic_fin['mask']['mean_dev'].items()]
    yys_p_3 = [el for _, el in dic_fin['pifpaf']['mean_3'].items()]
    yys_p_8 = [el for _, el in dic_fin['pifpaf']['mean_8'].items()]
    yys_p_4 = [el for _, el in dic_fin['pifpaf']['mean_4'].items()]
    # yys_m_3 = [el for _, el in dic_fin['mask']['mean_3'].items()]
    yys_ms = [el for _, el in dic_fin['monstereo'].items()]
    yys_p_best = [el for _, el in dic_fin['pifpaf']['mean_best'].items()]
    plt.plot(xxs, yys_p_4, marker='o', linestyle=':', label="PifPaf (highest 4)")
    plt.plot(xxs, yys_p, marker='+', label="PifPaf (median)")
    # plt.plot(xxs, yys_m, marker='o', label="Mask R-CNN (median")
    plt.plot(xxs, yys_p_3, marker='s', linestyle='--', label="PifPaf (closest 3)")
    plt.plot(xxs, yys_p_8, marker='*', linestyle=':', label="PifPaf (highest 8)")
    plt.plot(xxs, yys_ms, marker='^', label="MonStereo")
    plt.plot(xxs, yys_p_best, marker='o', label="PifPaf (best)")
    # plt.plot(xxs, yys_m_3, marker='o', color='r', label="Mask R-CNN (closest 3)")
    # plt.plot(xxs, yys_mon, marker='o', color='b', label="Our MonStereo")

    plt.legend()
    plt.tight_layout()
    path_fig = os.path.join(dir_out, 'mean_deviation.png')
    plt.savefig(path_fig)
    print("Figure of mean deviation saved in {}".format(path_fig))

    plt.figure(2)
    plt.xlabel("Ground-truth distance [m]")
    plt.ylabel("Pixels")
    plt.title("Standard deviation of joints disparity")
    yys_p = [el for _, el in dic_fin['pifpaf']['std_d'].items()]
    yys_m = [el for _, el in dic_fin['mask']['std_d'].items()]
    # yys_p_z = [el for _, el in dic_fin['pifpaf']['std_z'].items()]
    # yys_m_z = [el for _, el in dic_fin['mask']['std_z'].items()]
    plt.plot(xxs, yys_p, marker='s', label="PifPaf")
    plt.plot(xxs, yys_m, marker='o', label="Mask R-CNN")
    # plt.plot(xxs, yys_p_z, marker='s', color='b', label="PifPaf (meters)")
    # plt.plot(xxs, yys_m_z, marker='o', color='r', label="Mask R-CNN (meters)")

    plt.grid(linewidth=0.2)
    plt.legend()
    path_fig = os.path.join(dir_out, 'std_joints.png')
    plt.savefig(path_fig)
    print("Figure of standard deviation of joints by distance in {}".format(path_fig))

    plt.figure(3)
    # plt.style.use('ggplot')
    width = 0.35
    xxs = np.arange(len(COCO_KEYPOINTS))
    yys_p = [el for _, el in dic_fin['pifpaf']['joints'].items()]
    yys_m = [el for _, el in dic_fin['mask']['joints'].items()]
    plt.bar(xxs, yys_p, width, color='C0', label='Pifpaf')
    plt.bar(xxs + width, yys_m, width, color='C1', label='Mask R-CNN')
    plt.ylim(0, 1)

    plt.xlabel("Keypoints")
    plt.title("Repeatability by keypoint type")

    plt.xticks(xxs + width / 2, xxs)
    plt.legend(loc='best')
    path_fig = os.path.join(dir_out, 'repeatability_2.png')
    plt.savefig(path_fig)
    plt.close('all')
    print("Figure of standard deviation of joints by keypointd in {}".format(path_fig))

    plt.figure(4)
    plt.xlabel("Ground-truth distance [m]")
    plt.ylabel("Confidence")
    plt.grid(linewidth=0.2)
    xxs = get_distances(clusters)
    yys_p_conf = [el for _, el in dic_fin['pifpaf']['conf'].items()]
    yys_p_conf_best = [el for _, el in dic_fin['pifpaf']['conf_best'].items()]
    yys_m_conf = [el for _, el in dic_fin['mask']['conf'].items()]
    yys_m_conf_best = [el for _, el in dic_fin['mask']['conf_best'].items()]
    plt.plot(xxs, yys_p_conf_best, marker='s', color='lightblue', label="PifPaf (best)")
    plt.plot(xxs, yys_p_conf, marker='s', color='b', label="PifPaf (mean)")
    plt.plot(xxs, yys_m_conf_best, marker='^', color='darkorange', label="Mask (best)")
    plt.plot(xxs, yys_m_conf, marker='o', color='r', label="Mask R-CNN (mean)")
    plt.legend()
    plt.tight_layout()
    path_fig = os.path.join(dir_out, 'confidence.png')
    plt.savefig(path_fig)
    print("Figure of confidence saved in {}".format(path_fig))

max_stars_repo_name: galvinw/fairmotdocker

max_issues_repo_name: galvinw/fairmotdocker

max_forks_repo_name: galvinw/fairmotdocker


text:
from octis.models.model import AbstractModel
import numpy as np
from gensim.models import nmf
import gensim.corpora as corpora
import octis.configuration.citations as citations
import octis.configuration.defaults as defaults

class NMF(AbstractModel):

    def __init__(self, num_topics=100, chunksize=2000, passes=1, kappa=1.0, minimum_probability=0.01, w_max_iter=200,
                 w_stop_condition=0.0001, h_max_iter=50, h_stop_condition=0.001, eval_every=10, normalize=True,
                 random_state=None, use_partitions=True):
        """
        Initialize NMF model

        Parameters
        ----------
        num_topics (int, optional) – Number of topics to extract.

        chunksize (int, optional) – Number of documents to be used in each
        training chunk.

        passes (int, optional) – Number of full passes over the
        training corpus. Leave at default passes=1 if your input
        is an iterator.

        kappa (float, optional) – Gradient descent step size.
        Larger value makes the model train faster, but could
        lead to non-convergence if set too large.

        minimum_probability – If normalize is True, topics with
        smaller probabilities are filtered out. If normalize is False,
        topics with smaller factors are filtered out. If set to None,
        a value of 1e-8 is used to prevent 0s.

        w_max_iter (int, optional) – Maximum number of iterations to
        train W per each batch.

        w_stop_condition (float, optional) – If error difference gets less
        than that, training of W stops for the current batch.

        h_max_iter (int, optional) – Maximum number of iterations to train
        h per each batch.

        h_stop_condition (float) – If error difference gets less than that,
        training of h stops for the current batch.

        eval_every (int, optional) – Number of batches after which l2 norm
        of (v - Wh) is computed. Decreases performance if set too low.

        normalize (bool or None, optional) – Whether to normalize the result.

        random_state ({np.random.RandomState, int}, optional) – Seed for
        random generator. Needed for reproducibility.
        """
        super().__init__()
        self.hyperparameters["num_topics"] = num_topics
        self.hyperparameters["chunksize"] = chunksize
        self.hyperparameters["passes"] = passes
        self.hyperparameters["kappa"] = kappa
        self.hyperparameters["minimum_probability"] = minimum_probability
        self.hyperparameters["w_max_iter"] = w_max_iter
        self.hyperparameters["w_stop_condition"] = w_stop_condition
        self.hyperparameters["h_max_iter"] = h_max_iter
        self.hyperparameters["h_stop_condition"] = h_stop_condition
        self.hyperparameters["eval_every"] = eval_every
        self.hyperparameters["normalize"] = normalize
        self.hyperparameters["random_state"] = random_state
        self.use_partitions = use_partitions

        self.id2word = None
        self.id_corpus = None
        self.update_with_test = False

    def info(self):
        """
        Returns model informations
        """
        return {
            "citation": citations.models_NMF,
            "name": "NMF, Non-negative Matrix Factorization"
        }

    def hyperparameters_info(self):
        """
        Returns hyperparameters informations
        """
        return defaults.NMF_gensim_hyperparameters_info

    def partitioning(self, use_partitions, update_with_test=False):
        """
        Handle the partitioning system to use and reset the model to perform
        new evaluations

        Parameters
        ----------
        use_partitions: True if train/set partitioning is needed, False
                        otherwise
        update_with_test: True if the model should be updated with the test set,
                          False otherwise
        """
        self.use_partitions = use_partitions
        self.update_with_test = update_with_test
        self.id2word = None
        self.id_corpus = None

    def train_model(self, dataset, hyperparameters=None, top_words=10):
        """
        Train the model and return output

        Parameters
        ----------
        dataset : dataset to use to build the model
        hyperparameters : hyperparameters to build the model
        top_words : if greather than 0 returns the most significant words
                 for each topic in the output
                 Default True

        Returns
        -------
        result : dictionary with up to 3 entries,
                 'topics', 'topic-word-matrix' and
                 'topic-document-matrix'
        """
        if hyperparameters is None:
            hyperparameters = {}
        if self.use_partitions:
            partition = dataset.get_partitioned_corpus(use_validation=False)
        else:
            partition = [dataset.get_corpus(), []]

        if self.id2word is None:
            self.id2word = corpora.Dictionary(dataset.get_corpus())
        if self.id_corpus is None:
            self.id_corpus = [self.id2word.doc2bow(
                document) for document in partition[0]]

        hyperparameters["corpus"] = self.id_corpus
        hyperparameters["id2word"] = self.id2word
        self.hyperparameters.update(hyperparameters)

        self.trained_model = nmf.Nmf(**self.hyperparameters)

        result = {}

        result["topic-word-matrix"] = self.trained_model.get_topics()

        if top_words > 0:
            topics_output = []
            for topic in result["topic-word-matrix"]:
                top_k = np.argsort(topic)[-top_words:]
                top_k_words = list(reversed([self.id2word[i] for i in top_k]))
                topics_output.append(top_k_words)
            result["topics"] = topics_output

        result["topic-document-matrix"] = self._get_topic_document_matrix()

        if self.use_partitions:
            new_corpus = [self.id2word.doc2bow(
                document) for document in partition[1]]
            if self.update_with_test:
                self.trained_model.update(new_corpus)
                self.id_corpus.extend(new_corpus)

                result["test-topic-word-matrix"] = self.trained_model.get_topics()

                if top_words > 0:
                    topics_output = []
                    for topic in result["test-topic-word-matrix"]:
                        top_k = np.argsort(topic)[-top_words:]
                        top_k_words = list(
                            reversed([self.id2word[i] for i in top_k]))
                        topics_output.append(top_k_words)
                    result["test-topics"] = topics_output

                result["test-topic-document-matrix"] = self._get_topic_document_matrix()
            else:
                result["test-topic-document-matrix"] = self._get_topic_document_matrix(new_corpus)
        return result

    def _get_topics_words(self, topk):
        """
        Return the most significative words for each topic.
        """
        topic_terms = []
        for i in range(self.hyperparameters["num_topics"]):
            topic_words_list = []
            for word_tuple in self.trained_model.get_topic_terms(i, topk):
                topic_words_list.append(self.id2word[word_tuple[0]])
            topic_terms.append(topic_words_list)
        return topic_terms

    def _get_topic_document_matrix(self, test_corpus=None):
        """
        Return the topic representation of the
        corpus
        """
        doc_topic_tuples = []

        if test_corpus is None:
            for document in self.id_corpus:
                doc_topic_tuples.append(
                    self.trained_model.get_document_topics(document, minimum_probability=0))
        else:
            for document in test_corpus:
                doc_topic_tuples.append(
                    self.trained_model.get_document_topics(document, minimum_probability=0))
        topic_document = np.zeros((
            self.hyperparameters["num_topics"],
            len(doc_topic_tuples)))

        for ndoc in range(len(doc_topic_tuples)):
            document = doc_topic_tuples[ndoc]
            for topic_tuple in document:
                topic_document[topic_tuple[0]][ndoc] = topic_tuple[1]
        return topic_document

max_issues_repo_name: mpalenciaolivar/OCTIS

max_forks_repo_name: mpalenciaolivar/OCTIS


text:
import gym
import pylab
import random
import numpy as np
from collections import deque
from keras.layers import Dense
from keras.optimizers import Adam
from keras.models import Sequential

EPISODES = 4000

class DQNAgent:
    def __init__(self, state_size, action_size):
        # Cartpole이 학습하는 것을 보려면 "True"로 바꿀 것
        self.render = True

        # state와 action의 크기를 가져와서 모델을 생성하는데 사용함
        self.state_size = state_size
        self.action_size = action_size

        # Cartpole DQN 학습의 Hyper parameter 들
        # deque를 통해서 replay memory 생성
        self.discount_factor = 0.99
        self.learning_rate = 0.001
        self.epsilon = 1.0
        self.epsilon_min = 0.005
        self.epsilon_decay = (self.epsilon - self.epsilon_min) / 50000
        self.batch_size = 64
        self.train_start = 1000
        self.memory = deque(maxlen=10000)

        # 학습할 모델과 타겟 모델을 생성
        self.model = self.build_model()
        self.target_model = self.build_model()
        # 학습할 모델을 타겟 모델로 복사 --> 타겟 모델의 초기화(weight를 같게 해주고 시작해야 함)
        self.update_target_model()

    # Deep Neural Network를 통해서 Q Function을 근사
    # state가 입력, 각 행동에 대한 Q Value가 출력인 모델을 생성
    def build_model(self):
        model = Sequential()
        model.add(Dense(32, input_dim=self.state_size, activation='relu', kernel_initializer='he_uniform'))
        model.add(Dense(16, activation='relu', kernel_initializer='he_uniform'))
        model.add(Dense(self.action_size, activation='linear', kernel_initializer='he_uniform'))
        model.summary()
        model.compile(loss='mse', optimizer=Adam(lr=self.learning_rate))
        return model

    # 일정한 시간 간격마다 타겟 모델을 현재 학습하고 있는 모델로 업데이트
    def update_target_model(self):
        self.target_model.set_weights(self.model.get_weights())

    # 행동의 선택은 현재 네트워크에 대해서 epsilon-greedy 정책을 사용
    def get_action(self, state):
        if np.random.rand() <= self.epsilon:
            return random.randrange(self.action_size)
        else:
            q_value = self.model.predict(state)
            return np.argmax(q_value[0])

    # <s,a,r,s'>을 replay_memory에 저장함
    def replay_memory(self, state, action, reward, next_state, done):
        if action == 2:
            action = 1
        self.memory.append((state, action, reward, next_state, done))
        if self.epsilon > self.epsilon_min:
            self.epsilon -= self.epsilon_decay
        # print(len(self.memory))

    # replay memory에서 batch_size 만큼의 샘플들을 무작위로 뽑아서 학습
    def train_replay(self):
        if len(self.memory) < self.train_start:
            return
        batch_size = min(self.batch_size, len(self.memory))
        mini_batch = random.sample(self.memory, batch_size)

        update_input = np.zeros((batch_size, self.state_size))
        update_target = np.zeros((batch_size, self.action_size))

        for i in range(batch_size):
            state, action, reward, next_state, done = mini_batch[i]
            target = self.model.predict(state)[0]

            # 큐러닝에서와 같이 s'에서의 최대 Q Value를 가져옴. 단, 타겟 모델에서 가져옴
            if done:
                target[action] = reward
            else:
                target[action] = reward + self.discount_factor * \
                                          np.amax(self.target_model.predict(next_state)[0])
            update_input[i] = state
            update_target[i] = target

        # 학습할 정답인 타겟과 현재 자신의 값의 minibatch를 만들고 그것으로 한 번에 모델 업데이트
        self.model.fit(update_input, update_target, batch_size=batch_size, epochs=1, verbose=0)

    # 저장한 모델을 불러옴
    def load_model(self, name):
        self.model.load_weights(name)

    # 학습된 모델을 저장함
    def save_model(self, name):
        self.model.save_weights(name)

if __name__ == "__main__":
    # CartPole-v1의 경우 500 타임스텝까지 플레이가능
    env = gym.make('MountainCar-v0')
    # 환경으로부터 상태와 행동의 크기를 가져옴
    state_size = env.observation_space.shape[0]
    #action_size = env.action_space.n
    action_size = 2
    # DQN 에이전트의 생성
    agent = DQNAgent(state_size, action_size)
    agent.load_model("./save_model/MountainCar_DQN.h5")
    scores, episodes = [], []

    for e in range(EPISODES):
        done = False
        score = 0
        state = env.reset()
        state = np.reshape(state, [1, state_size])
        print(state)

        # 액션 0(좌), 1(아무것도 안함), 3(아무것도 하지 않는 액션을 하지 않기 위한 fake_action 선언
        fake_action = 0

        # 같은 액션을 4번하기 위한 카운터
        action_count = 0

        while not done:
            if agent.render:
                env.render()

            # 현재 상태에서 행동을 선택하고 한 스텝을 진행
            action_count = action_count + 1

            if action_count == 4:
                action = agent.get_action(state)
                action_count = 0

                if action == 0:
                    fake_action = 0
                elif action == 1:
                    fake_action = 2

            # 선택한 액션으로 1 step을 시행한다
            next_state, reward, done, info = env.step(fake_action)
            next_state = np.reshape(next_state, [1, state_size])
            # 에피소드를 끝나게 한 행동에 대해서 -100의 패널티를 줌
            #reward = reward if not done else -100

            # <s, a, r, s'>을 replay memory에 저장
            agent.replay_memory(state, fake_action, reward, next_state, done)
            # 매 타임스텝마다 학습을 진행
            agent.train_replay()
            score += reward
            state = next_state

            if done:
                env.reset()
                # 매 에피소드마다 학습하는 모델을 타겟 모델로 복사
                agent.update_target_model()

                # 각 에피소드마다 cartpole이 서있었던 타임스텝을 plot
                scores.append(score)
                episodes.append(e)
                #pylab.plot(episodes, scores, 'b')
                #pylab.savefig("./save_graph/MountainCar_DQN.png")
                print("episode:", e, "  score:", score, "  memory length:", len(agent.memory),
                      "  epsilon:", agent.epsilon)

        # 50 에피소드마다 학습 모델을 저장
        if e % 50 == 0:
             agent.save_model("./save_model/MountainCar_DQN.h5")


text:
import os
from setuptools import setup, Extension, Command
from Cython.Distutils import build_ext
from Cython.Build import cythonize
import numpy as np
from os.path import join as pjoin

mpi_compile_args = os.popen("mpic++ --showme:compile").read().strip().split(' ')
mpi_link_args = os.popen("mpic++ --showme:link").read().strip().split(' ')

class cuda_setup():

    def find_in_path(name, path):
        """Find a file in a search path"""

        # Adapted fom http://code.activestate.com/recipes/52224
        for dir in path.split(os.pathsep):
            binpath = pjoin(dir, name)
            if os.path.exists(binpath):
                return os.path.abspath(binpath)
        return None

    def locate_cuda(self):
        """Locate the CUDA environment on the system
        Returns a dict with keys 'home', 'nvcc', 'include', and 'lib64'
        and values giving the absolute path to each directory.
        Starts by looking for the CUDAHOME env variable. If not found,
        everything is based on finding 'nvcc' in the PATH.
        """
        CUDA = self.locate_cuda()

        # First check if the CUDAHOME env variable is in use
        if 'CUDAHOME' in os.environ:
            home = os.environ['CUDAHOME']
            nvcc = pjoin(home, 'bin', 'nvcc')
        else:
            # Otherwise, search the PATH for NVCC
            nvcc = self.find_in_path('nvcc', os.environ['PATH'])
            if nvcc is None:
                raise EnvironmentError('The nvcc binary could not be '
                                       'located in your $PATH. Either add it to your path, '
                                       'or set $CUDAHOME')
            home = os.path.dirname(os.path.dirname(nvcc))

        cudaconfig = {'home': home, 'nvcc': nvcc,
                      'include': pjoin(home, 'include'),
                      'lib64': pjoin(home, 'lib64')}
        for k, v in iter(cudaconfig.items()):
            if not os.path.exists(v):
                raise EnvironmentError('The CUDA %s path could not be '
                                       'located in %s' % (k, v))
        return cudaconfig

    try:
        numpy_include = np.get_include()
    except AttributeError:
        numpy_include = np.get_numpy_include()

    def get_module(self):
        CUDA = self.locate_cuda()
        return Extension('da', sources=['da.cu', 'src/da.py'],
                         include_dirs=[np.get_include(), CUDA['include']],
                         library_dirs=['/usr/include/boost/', CUDA['lib64']],
                         language='c++',
                         libraries=['/usr/include/boost/chrono', 'cudart', '/usr/lib/x86_64-linux-gnu/openmpi/lib'],
                         extra_compile_args={'gcc': [],
                                             'nvcc': [
                                                 '-arch=sm_61', '--ptxas-options=-v', '-c',
                                                 '--compiler-options', "'-fPIC'", '-lmpi_cxx', '-lmpi', '-Xcompiler',
                                                 '-fopenmp'
                                             ], 'openmp': "-fopenmp", 'mpi': ' '.join(mpi_compile_args)
                                             },  # (["-fopenmp"] + mpi_compile_args),
                         runtime_library_dirs=[CUDA['lib64']],
                         extra_link_args=(["-fopenmp"] + mpi_link_args)
                         )

# Run the customize_compiler
class custom_build_ext(build_ext):
    def build_extensions(self):
        self.customize_compiler_for_nvcc(self.compiler)
        build_ext.build_extensions(self)

def customize_compiler_for_nvcc(self):
    """Inject deep into distutils to customize how the dispatch
    to gcc/nvcc works.
    If you subclass UnixCCompiler, it's not trivial to get your subclass
    injected in, and still have the right customizations (i.e.
    distutils.sysconfig.customize_compiler) run on it. So instead of going
    the OO route, I have this. Note, it's kindof like a wierd functional
    subclassing going on.
    """

    # Tell the compiler it can processes .cu
    self.src_extensions.append('.cu')

    # Save references to the default compiler_so and _comple methods
    default_compiler_so = self.compiler_so
    super = self._compile
    CUDA = self.locate_cuda()

    # Now redefine the _compile method. This gets executed for each
    # object but distutils doesn't have the ability to change compilers
    # based on source extension: we add it.
    def _compile(obj, src, ext, cc_args, extra_postargs, pp_opts):
        if os.path.splitext(src)[1] == '.cu':
            # use the cuda for .cu files
            self.set_executable('compiler_so', CUDA['nvcc'])
            # use only a subset of the extra_postargs, which are 1-1
            # translated from the extra_compile_args in the Extension class
            # /usr/lib/x86_64-linux-gnu/openmpi/include/
            # -I/usr/lib/x86_64-linux-gnu/openmpi/include/openmpi -I/usr/lib/x86_64-linux-gnu/openmpi/include -pthread
            print(extra_postargs['mpi'])
            postargs = extra_postargs['nvcc'] + mpi_compile_args
        else:
            postargs = extra_postargs['gcc'] + extra_postargs['openmp'] + mpi_compile_args

        super(obj, src, ext, cc_args, postargs, pp_opts)
        # Reset the default compiler_so, which we might have changed for cuda
        self.compiler_so = default_compiler_so

    # Inject our redefined _compile method into the class
    self._compile = _compile


text:
from matplotlib import pyplot as plt
import pandas as pd
import numpy as np

df = pd.read_excel("comparison.xlsx")

E_eed = df["E_eed"]
E_gpot = df["E_gpot"]
E_comsol = df["E_comsol"]

bins = np.arange(-2.6, 1.2, 0.2)

print(E_comsol-E_eed)
print(E_comsol-E_gpot)

plt.hist(E_comsol-E_eed, bins, fc="None", ec="b")
#plt.hist(E_comsol-E_gpot, bins, fc="None", ec="r")
#plt.gca().set_xticks(bins)


text:
import requests
import sys
import os
import logging
import coloredlogs
import zipfile
import pandas as pd
import shutil
from scipy import spatial
from typing import List

class LatLongResolver():
    def __init__(self, cities_model_filename: str):
        if not os.path.exists(cities_model_filename):
            self.build_cities_dataset(cities_model_filename)
        self.cities_filename = cities_model_filename
        self.pd = pd.read_csv(self.cities_filename)
        self.init()

    def init(self):
        # build several trees based on population
        self.psmall = self.pd[(self.pd.population >= 5000) & (self.pd.population < 20000)]
        self.pmedium = self.pd[(self.pd.population >= 20000) & (self.pd.population < 1000000)]
        self.plarge = self.pd[(self.pd.population >= 1000000) & (self.pd.population < 10000000)]

        self.t5000 = spatial.KDTree(list(zip(self.psmall.latitude,
                                             self.psmall.longitude)))

        self.t6000 = spatial.KDTree(list(zip(self.pmedium.latitude,
                                             self.pmedium.longitude)))

        self.t7000 = spatial.KDTree(list(zip(self.plarge.latitude,
                                             self.plarge.longitude)))

    def nearest_helper(self, dataset: pd.DataFrame, tree: spatial.KDTree, latitude: float, longitude: float) -> List[str]:
        distance, index = tree.query([(latitude, longitude)])
        result = dataset.iloc[index]
        country = result.Country.item()
        city = result.asciiname.item()
        latitude = result.latitude.item()
        longitude = result.longitude.item()
        return [city, country]

    def build_cities_dataset(self, model_filename: str) -> pd.DataFrame:
        if os.path.exists(model_filename):
            return pd.read_csv(model_filename)

        download_url = 'http://download.geonames.org/export/dump/'

        if not os.path.exists('tmp'):
            os.mkdir('tmp')

        logging.info('downloading {}'.format(download_url + 'cities500.zip'))
        r = requests.get(download_url + 'cities500.zip')
        with open('tmp/cities500.zip', 'wb') as f:
            f.write(r.content)

        logging.info('downloading {}'.format(download_url + 'countryInfo.txt'))
        r2 = requests.get(download_url + 'countryInfo.txt')
        with open('tmp/countryInfo.txt', 'wb') as f2:
            f2.write(r2.content)

        with zipfile.ZipFile('tmp/cities500.zip', 'r') as zip_ref:
            zip_ref.extractall('tmp')

        cities = pd.read_csv('tmp/cities500.txt', sep='\t',
                             names=['geonameid', 'name', 'asciiname', 'alternatenames',
                                    'latitude', 'longitude', 'feature class', 'feature code', 'country code', 'cc2',
                                    'admin1 code', 'admin2 code', 'admin3 code', 'admin4 code',
                                    'population', 'elevation', 'dem', 'timezone', 'modification date'])

        # remove all the comment junk at the beginning of the file
        text = ''
        with open('tmp/countryInfo.txt', 'r') as f3:
            text = f3.read()
            text = text[text.index('#ISO') + 1:]
        with open('tmp/countryInfo_rewritten.txt', 'w') as f4:
            f4.write(text)

        countries = pd.read_csv('tmp/countryInfo_rewritten.txt', sep='\t')

        shaped_cities = cities[['asciiname', 'latitude', 'longitude', 'country code', 'population']]
        shaped_countries = countries[['ISO', 'Country']]
        result = shaped_cities.merge(shaped_countries, left_on='country code', right_on='ISO')

        result.to_csv(model_filename, header=True, index=False)
        shutil.rmtree('tmp')

        return result

    def nearest(self, latitude: float, longitude: float) -> List[str]:
        result = [self.nearest_helper(self.psmall, self.t5000, latitude, longitude),
                  self.nearest_helper(self.pmedium, self.t6000, latitude, longitude),
                  self.nearest_helper(self.plarge, self.t7000, latitude, longitude)]

        flat_list = [item for sublist in result for item in sublist]
        return list(set(flat_list))

coloredlogs.install(level='INFO')

max_stars_repo_name: joelpob/photos-gallery

max_issues_repo_name: joelpob/photos-gallery

max_forks_repo_name: joelpob/photos-gallery

reason: from scipy


text:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as spst

import proppa
from utilities import gillespie,normalise_trace,extract_states

def sample_dist(dist):
    try:
        return dist.rvs()
    except AttributeError:
        ind = np.random.random_integers(low=0,high=len(dist)-1)
        return dist[ind]
#   alternatively:
#    if 'rvs' in dir(dist):
#        return dist.rvs()
#    else:
#        ind = np.random.random_integers(len(dist)-1)
#        return dist[ind]

def debug_dist(dists,n_samples=5000):
    n_params = len(dists)
    samples = [dists[i].rvs(size=n_samples) for i in range(n_params)]
    for i in range(n_params):
        plt.hist(samples[i])
        plt.show()

def sample_paths(model,t_final=None,dists=None,n_paths=100):
    # retrieve rate functions, initial state, stop time and updates from model
    abstract_rate_funcs = model.reaction_functions()
    updates = model.updates
    init_state = model.init_state
    if dists is None:
        print("Using priors")
        dists = [p.rhs.to_distribution() for p in model.uncertain]
    if t_final is None: # use the final observation time
        t_final = model.obs[-1][0]

    paths = []
    i = 0
    while i < n_paths:
        params = [sample_dist(d) for d in dists]
        rate_funcs = [r(params) for r in abstract_rate_funcs]
        path = gillespie(rate_funcs,t_final,init_state,updates)
        paths.append(path)
        i = i + 1
    return paths

if __name__ == "__main__":
    # toy model and distribution
    model = proppa.load_model('predPreyTest.proppa')
    model.numerize()
    t_f = 1
    paths = sample_paths(model,t_f,n_paths=1000)
#    fixed_values = [5,5,10,10]
#    dists = [spst.rv_discrete(values=(v,1)) for v in fixed_values]
#    paths = sample_paths(model,t_f,n_paths=1000,dists=dists)

    for p in paths:
        p_t, p_y = zip(*[(pp[0],pp[1][0]) for pp in p])
        plt.step(p_t,p_y)
    plt.title('Paths of prey (%d SSA runs)' % len(paths))
    plt.show()

    plot_times = np.linspace(0,t_f,101)
    norm_paths = [normalise_trace(p,plot_times) for p in paths]
    norm_prey = [ [e[0] for e in extract_states(p)] for p in norm_paths]
    avg_path = np.average(norm_prey,axis=0)
    std_path = np.std(norm_prey,axis=0)

    upper_line = avg_path + std_path
#    lower_line = avg_path - std_path
    lower_line = np.maximum(avg_path - std_path,0)
    plt.plot(plot_times,avg_path,lw=2)
    plt.plot(plot_times,upper_line,'k--',lw=2)
    plt.plot(plot_times,lower_line,'k--',lw=2)
    plt.fill_between(plot_times,upper_line,lower_line,color='grey',alpha='0.5')
    plt.title('Average prey (+/- 1 std)')
    plt.show()

    min_prey = np.min(norm_prey,axis=0)
    max_prey = np.max(norm_prey,axis=0)
    plt.plot(plot_times,avg_path)
    plt.plot(plot_times,min_prey,'k--')
    plt.plot(plot_times,max_prey,'k--')
    plt.title('Mean, min and max prey (%d SSA runs)' % len(paths))
    plt.show()

reason: import numpy,import scipy


text:
import csv
import json
import numpy
import tensorflow as tf

class Song:
	def __init__(self,number, artist,name, dance ,energy, loud, duration):
		self.number=number
		self.artist=artist
		self.name=name
		self.dance=dance
		self.energy=energy
		self.loud=loud
		self.duration=duration

def Wishman_score(song,user):
	da=song.dance/user.dance
	en=song.energy/user.energy
	lo=song.loud/user.loud
	du=song.duration/user.duration

	wishman=da+(en+(lo+du))
	return wishman

#gathering data and sorting it
songs=[]
with open("songs.csv","r") as f:
	data=csv.reader(f)
	header=next(data)
	size=1
	for line in data:
		artist=line[2]
		name=line[1]
		dance=float(line[3])
		energy=float(line[4])
		loud=float(line[6])
		duration=float(line[14])
		song=Song(size,artist,name,dance,energy, loud,duration)
		songs.append(song)
		size=size+1

#crete initial play list
bluff=1

#start running loop
while(bluff>0):
	i=bluff-1
	user=Song(i,songs[i].artist, songs[i].name, songs[i].dance, songs[i].energy, songs[i].loud, songs[i].duration)
	print("currently playing:   " + user.artist + " - " + user.name)

	#create recomendation play list
	print("recomended play list:")
	wishmans=[]
	for song in songs:
		wishmans.append(Wishman_score(song,user))

	#picking the songs with the best wishmans
	#ideal song has wishmn score of 5
	#recomendation should pick songs most close to 5

	closeness=[]
	for wishman in wishmans:
		if(wishman<4):
			closeness.append(4-wishman)
		else:
			closeness.append(wishman-4)

	#picking up the songs
	song_numbers=[]
	x=0
	#those greater close, 6 songs
	sub1=closeness
	while(x<7):
		j=1
		number=0
		last=100
		for n in sub1:
			if((n<last)&(n!=-1)):
				last=n
				number=j
			else:
				pass
			j=j+1
		song_numbers.append(number)
		x=x+1
		#remove the found number
		j=1
		sub2=[]
		for k in sub1:
			if(j==number):
				sub2.append(-1)
			else:
				sub2.append(k)
			j=j+1
		sub1=sub2

	for n in song_numbers:
		print(str(songs[n-1].number) + "  " + songs[n-1].artist + " - " + songs[n-1].name)

	print("Please select song number or 0 to quit")
	bluff=int(input())

max_stars_repo_name: mershcap/project

max_issues_repo_name: mershcap/project

max_forks_repo_name: mershcap/project