text:
#######################
# Header.R

from datetime import time
from operator import index
from os import path, times
import numpy as np
import pandas as pd
import os
import logging

from pathlib import Path

from pandas.core.reshape.merge import merge
from powergenome.util import regions_to_keep
from powergenome.us_state_abbrev import (state2abbr, abbr2state)
path_in = r"..\data\load_profiles_data\input" # fix

#read in state proportions
#how much state load should be distributed to GenXRegion
# pop = pd.read_parquet(path_in + "\GenX_State_Pop_Weight.parquet")
pop = pd.read_parquet(path_in + "\ipm_state_pop_weight_20210517.parquet")
states = pop.drop_duplicates(subset=["State"])["State"]
states_abb = list(map(state2abbr, states))
pop["State"] = list(map(state2abbr, pop["State"]))
states_eastern_abbr = ["ME","VT","NH","MA","RI","CT","NY","PA","NJ","DE","MD","DC","MI","IN","OH","KY","WV","VA","NC","SC","GA","FL"]
states_central_abbr = ["IL","MO","TN","AL","MS","WI","AR","LA","TX","OK","KS","NE","SD","ND","IA","MN"]
states_mountain_abbr = ["MT","WY","CO","NM","AZ","UT","ID"]
states_pacific_abbr = ["CA","NV","OR","WA"]
states_eastern = list(map(abbr2state, states_eastern_abbr))
states_central = list(map(abbr2state, states_central_abbr))
states_mountain = list(map(abbr2state, states_mountain_abbr))
states_pacific = list(map(abbr2state, states_pacific_abbr))

# some parameters
stated_states = ["New Jersey", "New York", "Virginia"]
# Date Jan 29, 2021
# (2) PA, NJ, VA, NY, MI all have EV and heat pump stocks from NZA DD case
# consistent with their economywide decarbonization goals.
# https://www.c2es.org/content/state-climate-policy/
# Date Feb 10, 2021
# Remove high electrification growth in PA and MI in stated policies;
# they dont have clean energy goals so kind of confusing/inconsistent to require high electrification in these states.
# So our new "Stated Policies" definition for electrification is states
# with BOTH economywide emissions goals + 100% carbon-free electricity standards
# = NY, NJ, VA.
stated_states_abbr = list(map(state2abbr, stated_states))
#years = ["2022", "2025", "2030", "2040", "2050"]
cases = ["current_policy", "stated_policy", "deep_decarbonization"]

running_sector = ['Residential','Residential', 'Commercial', 'Commercial','Transportation','Transportation','Transportation', 'Transportation']
running_subsector = ['space heating and cooling','water heating', 'space heating and cooling', 'water heating','light-duty vehicles','medium-duty trucks','heavy-duty trucks','transit buses']
Nsubsector = len(running_subsector)

logger = logging.getLogger(__name__)

#Define function for adjusting time-difference
def addhour(x):
    x += 1
    x = x.replace(8761, 1)
    return x

def SolveThreeUnknowns(a1, b1, c1, d1, a2, b2, c2, d2, a3, b3, c3, d3):
    D  = a1*b2*c3 + b1*c2*a3 + c1*a2*b3 - a1*c2*b3 - b1*a2*c3 - c1*b2*a3
    Dx = d1*b2*c3 + b1*c2*d3 + c1*d2*b3 - d1*c2*b3 - b1*d2*c3 - c1*b2*d3
    Dy = a1*d2*c3 + d1*c2*a3 + c1*a2*d3 - a1*c2*d3 - d1*a2*c3 - c1*d2*a3
    Dz = a1*b2*d3 + b1*d2*a3 + d1*a2*b3 - a1*d2*b3 - b1*a2*d3 - d1*b2*a3
    Sx = Dx/D
    Sy = Dy/D
    Sz = Dz/D
    d = {'Sx':Sx, 'Sy':Sy, 'Sz':Sz}
    return pd.DataFrame(d)

def SolveTwoUnknowns(a1, b1, c1, a2, b2, c2):
    D  = a1*b2 - a2*b1
    Dx = c1*b2 - c2*b1
    Dy = a1*c2 - a2*c1
    Sx = Dx/D
    Sy = Dy/D
    d = {'Sx':Sx, 'Sy':Sy}
    return pd.DataFrame(d)

def CreateOutputFolder(case_folder):
    path = case_folder / "extra_outputs"
    if not os.path.exists(path):
        os.makedirs(path)

######################################
# CreatingBaseLoad.R
def CreateBaseLoad(years, regions, output_folder, path_growthrate):
    path_processed = path_in
    path_result = output_folder.__str__()
    years = years
    regions = regions
    path_growthrate = path_growthrate

    ## Method 3: annually
    EFS_2020_LoadProf = pd.read_parquet(path_in + "\EFS_REF_load_2020.parquet")
    EFS_2020_LoadProf = pd.merge(EFS_2020_LoadProf, pop, on = ["State"])
    EFS_2020_LoadProf = EFS_2020_LoadProf.assign(weighted = EFS_2020_LoadProf["LoadMW"]*EFS_2020_LoadProf["State Prop"])
    EFS_2020_LoadProf = EFS_2020_LoadProf.groupby(["Year", "GenX.Region", "LocalHourID", "Sector", "Subsector"], as_index = False).agg({"weighted" : "sum"})

    # Read in 2019 Demand
    Original_Load_2019 = pd.read_parquet(path_in + "\ipm_load_curves_2019_EST.parquet")
    # Reorganize Demand
    Original_Load_2019 = Original_Load_2019.melt(id_vars="LocalHourID").rename(columns={"variable" : "GenX.Region", "value": "LoadMW_original"})
    Original_Load_2019 = Original_Load_2019.groupby(["LocalHourID"], as_index = False).agg({"LoadMW_original" : "sum"})

    ratio_A = Original_Load_2019["LoadMW_original"].sum() / EFS_2020_LoadProf["weighted"].sum()
    EFS_2020_LoadProf = EFS_2020_LoadProf.assign(weighted = EFS_2020_LoadProf["weighted"]*ratio_A)

    Base_Load_2019 = EFS_2020_LoadProf.rename(columns ={"weighted" : "LoadMW"})
    breakpoint()
    # Read in the Growth Rate
    GrowthRate = pd.read_parquet(path_in + "\ipm_growthrate_2019.parquet")
    try:
        GrowthRate = pd.read_parquet(path_growthrate)
    except:
        pass

    # Create Base loads
    Base_Load_2019 = Base_Load_2019[Base_Load_2019["GenX.Region"].isin(regions)]
    Base_Load_2019.loc[(Base_Load_2019["Sector"] == "Industrial") & (Base_Load_2019["Subsector"].isin(["process heat", "machine drives"])), "Subsector"] = "other"
    Base_Load_2019 = Base_Load_2019[Base_Load_2019["Subsector"] == "other"]
    Base_Load_2019 = Base_Load_2019.groupby(["Year", "LocalHourID", "GenX.Region", "Sector"], as_index= False).agg({'LoadMW' : 'sum'})
    Base_Load = Base_Load_2019
    for y in years:
        ScaleFactor = GrowthRate.assign(ScaleFactor = (1+GrowthRate["growth_rate"])**(int(y) - 2019)) \
            .drop(columns = "growth_rate")
        Base_Load_temp = pd.merge(Base_Load_2019, ScaleFactor, on = ["GenX.Region"])
        Base_Load_temp = Base_Load_temp.assign(Year = y, LoadMW = Base_Load_temp["LoadMW"]*Base_Load_temp["ScaleFactor"])\
                .drop(columns = "ScaleFactor")
        Base_Load = Base_Load.append(Base_Load_temp, ignore_index=True)
    Base_Load.to_parquet(path_result + "\Base_Load.parquet", index = False)
    del Base_Load, Base_Load_2019, Base_Load_temp, ScaleFactor,GrowthRate, Original_Load_2019

#####################################
# Add_Electrification.R
def AddElectrification(years, regions, electrification, output_folder, path_stock):
    path_processed = path_in
    path_result = output_folder.__str__()
    path_stock = path_stock
    years = years
    electrification = electrification
    regions = regions
    #Creating Time-series

    SCENARIO_STOCK = pd.read_parquet(path_processed + "\SCENARIO_STOCK.parquet")
    SCENARIO_STOCK = SCENARIO_STOCK[(SCENARIO_STOCK["YEAR"].isin(years)) & (SCENARIO_STOCK["SCENARIO"].isin(electrification))]
    SCENARIO_STOCK_temp = pd.DataFrame()
    for year, case in zip(years, electrification):
        SCENARIO_STOCK_temp = SCENARIO_STOCK_temp.append(SCENARIO_STOCK[(SCENARIO_STOCK["YEAR"] == year) & (SCENARIO_STOCK["SCENARIO"] == case)])
    SCENARIO_STOCK = SCENARIO_STOCK_temp
    del SCENARIO_STOCK_temp

    try:
        CUSTOM_STOCK = pd.read_parquet(path_stock)
        CUSTOM_STOCK = CUSTOM_STOCK[(CUSTOM_STOCK["YEAR"].isin(years)) & (CUSTOM_STOCK["SCENARIO"].isin(electrification))]
        SCENARIO_STOCK = SCENARIO_STOCK.append(CUSTOM_STOCK)
    except:
        pass

    #Method 1 Calculate from Type1 and Type 2
    for i in range(0, Nsubsector):
        timeseries = pd.read_parquet(path_processed + "\\" + running_sector[i] + "_" + running_subsector[i] + "_Incremental_Factor.parquet")
        timeseries = timeseries[["State", "Year", "LocalHourID", "Unit", "Factor_Type1", "Factor_Type2" ]]
        stock_temp = SCENARIO_STOCK[(SCENARIO_STOCK["SECTOR"] == running_sector[i]) & (SCENARIO_STOCK["SUBSECTOR"] == running_subsector[i])]
        stock_temp = stock_temp[["SCENARIO", "STATE", "YEAR", "AGG_STOCK_TYPE1", "AGG_STOCK_TYPE2"]].rename(columns={"STATE" : "State", "YEAR" : "Year"})
        years_pd = pd.Series(years)
        IF_years = pd.Series(timeseries["Year"].unique())
        for year in years_pd:
            exists = year in IF_years.values
            if not exists:
                diff = np.array(IF_years - year)
                index = diff[np.where(diff <= 0)].argmax()
                year_approx = IF_years[index]
                timeseries_temp = timeseries[timeseries["Year"] == year_approx]
                timeseries_temp["Year"] = year
                logger.warning("No incremental factor available for year " + str(year) + ": using factors from year " + str(year_approx) + ".")
                timeseries = timeseries.append(timeseries_temp)

        timeseries = pd.merge(timeseries, stock_temp, on = ["State", "Year"])
        timeseries = timeseries.assign(LoadMW = timeseries["AGG_STOCK_TYPE1"]*timeseries["Factor_Type1"] + timeseries["AGG_STOCK_TYPE2"]*timeseries["Factor_Type2"])
        timeseries = timeseries[["SCENARIO", "State", "Year", "LocalHourID", "LoadMW"]].dropna()
        timeseries.to_parquet(path_result + "\\" + running_sector[i] + "_" + running_subsector[i] + "_Scenario_Timeseries_Method1.parquet", index = False)
    del timeseries, stock_temp

    ##########################
    # Read in time series and combine them
    Method = "Method1"
    Res_SPH = pd.read_parquet(path_result + "\Residential_space heating and cooling_Scenario_Timeseries_" + Method + ".parquet")
    Res_SPH = Res_SPH.rename(columns={"LoadMW" : "Res_SPH_LoadMW"})
    Res_SPH_sum = Res_SPH
    Res_SPH_sum = Res_SPH.groupby(["SCENARIO", "State", "Year"], as_index = False)["Res_SPH_LoadMW"].agg({"Total_Res_SPH_TWh" : "sum"})
    Res_SPH_sum["Total_Res_SPH_TWh"] = 10**-6*Res_SPH_sum["Total_Res_SPH_TWh"]

    Res_WH = pd.read_parquet(path_result + "\Residential_water heating_Scenario_Timeseries_" + Method +".parquet")
    Res_WH = Res_WH.rename(columns ={"LoadMW" : "Res_WH_LoadMW"})
    Res_WH_sum = Res_WH
    Res_WH_sum = Res_WH.groupby(["SCENARIO", "State", "Year"], as_index = False)["Res_WH_LoadMW"].agg({"Total_Res_WH_TWh" : "sum"})
    Res_WH_sum["Total_Res_WH_TWh"] = 10**-6*Res_WH_sum["Total_Res_WH_TWh"]

    Com_SPH = pd.read_parquet(path_result + "\Commercial_space heating and cooling_Scenario_Timeseries_" + Method +".parquet")
    Com_SPH = Com_SPH.rename(columns={"LoadMW" : "Com_SPH_LoadMW"})
    Com_SPH_sum = Com_SPH
    Com_SPH_sum = Com_SPH.groupby(["SCENARIO", "State", "Year"], as_index = False)["Com_SPH_LoadMW"].agg({"Total_Com_SPH_TWh" : "sum"})
    Com_SPH_sum["Total_Com_SPH_TWh"] = 10**-6*Com_SPH_sum["Total_Com_SPH_TWh"]

    Com_WH = pd.read_parquet(path_result + "\Commercial_water heating_Scenario_Timeseries_" + Method +".parquet")
    Com_WH = Com_WH.rename(columns ={"LoadMW" : "Com_WH_LoadMW"})
    Com_WH_sum = Com_WH
    Com_WH_sum = Com_WH.groupby(["SCENARIO", "State", "Year"], as_index = False)["Com_WH_LoadMW"].agg({"Total_Com_WH_TWh" : "sum"})
    Com_WH_sum["Total_Com_WH_TWh"] = 10**-6*Com_WH_sum["Total_Com_WH_TWh"]

    Trans_LDV = pd.read_parquet(path_result + "\Transportation_light-duty vehicles_Scenario_Timeseries_" + Method +".parquet")
    Trans_LDV = Trans_LDV.rename(columns ={"LoadMW" : "LDV_LoadMW"})
    Trans_LDV_sum = Trans_LDV
    Trans_LDV_sum = Trans_LDV.groupby(["SCENARIO", "State", "Year"], as_index = False)["LDV_LoadMW"].agg({"Total_Trans_LDV_TWh" : "sum"})
    Trans_LDV_sum["Total_Trans_LDV_TWh"] = 10**-6*Trans_LDV_sum["Total_Trans_LDV_TWh"]

    Trans_MDV = pd.read_parquet(path_result + "\Transportation_medium-duty trucks_Scenario_Timeseries_" + Method +".parquet")
    Trans_MDV = Trans_MDV.rename(columns ={"LoadMW" : "MDV_LoadMW"})
    Trans_MDV_sum = Trans_MDV
    Trans_MDV_sum = Trans_MDV.groupby(["SCENARIO", "State", "Year"], as_index = False)["MDV_LoadMW"].agg({"Total_Trans_MDV_TWh" : "sum"})
    Trans_MDV_sum["Total_Trans_MDV_TWh"] = 10**-6*Trans_MDV_sum["Total_Trans_MDV_TWh"]

    Trans_HDV = pd.read_parquet(path_result + "\Transportation_heavy-duty trucks_Scenario_Timeseries_" + Method +".parquet")
    Trans_HDV = Trans_HDV.rename(columns ={"LoadMW" : "HDV_LoadMW"})
    Trans_HDV_sum = Trans_HDV
    Trans_HDV_sum = Trans_HDV.groupby(["SCENARIO", "State", "Year"], as_index = False)["HDV_LoadMW"].agg({"Total_Trans_HDV_TWh" : "sum"})
    Trans_HDV_sum["Total_Trans_HDV_TWh"] = 10**-6*Trans_HDV_sum["Total_Trans_HDV_TWh"]

    Trans_BUS = pd.read_parquet(path_result + "\Transportation_transit buses_Scenario_Timeseries_" + Method +".parquet")
    Trans_BUS = Trans_BUS.rename(columns ={"LoadMW" : "BUS_LoadMW"})
    Trans_BUS_sum = Trans_BUS
    Trans_BUS_sum = Trans_BUS.groupby(["SCENARIO", "State", "Year"], as_index = False)["BUS_LoadMW"].agg({"Total_Trans_BUS_TWh" : "sum"})
    Trans_BUS_sum["Total_Trans_BUS_TWh"] = 10**-6*Trans_BUS_sum["Total_Trans_BUS_TWh"]

    # Total_SPH = pd.merge(Res_SPH, Com_SPH, on = ["SCENARIO","State","Year","LocalHourID"])
    # Total_SPH = Total_SPH.assign(LoadMW = Total_SPH['Res_SPH_LoadMW'] + Total_SPH["Com_SPH_LoadMW"]).drop(columns = ["Res_SPH_LoadMW", "Com_SPH_LoadMW"])
    # Total_WH = pd.merge(Res_WH, Com_WH, on = ["SCENARIO","State","Year","LocalHourID"])
    # Total_WH = Total_WH.assign(LoadMW = Total_WH['Res_WH_LoadMW'] + Total_WH["Com_WH_LoadMW"]).drop(columns = ["Res_WH_LoadMW", "Com_WH_LoadMW"])
    # Total_MHDV = pd.merge(Trans_MDV, Trans_HDV, on = ["SCENARIO", "State", "Year", "LocalHourID"])
    # Total_MHDV = pd.merge(Total_MHDV, Trans_BUS, on = ["SCENARIO", "State", "Year", "LocalHourID"])
    # Total_MHDV = Total_MHDV.assign(LoadMW = Total_MHDV["MDV_LoadMW"] + Total_MHDV["HDV_LoadMW"] + Total_MHDV["BUS_LoadMW"]).drop(columns=["MDV_LoadMW", "HDV_LoadMW", "BUS_LoadMW"])

    #Total_SPH.to_parquet(path_or_buf = path_result + "\Total_SPH_" + Method + ".parquet", index=False)
    #Total_WH.to_parquet(path_or_buf = path_result + "\Total_WH_" + Method + ".parquet", index = False)
    #Trans_LDV.to_parquet(path_or_buf = path_result + "\Trans_LDV_" + Method + ".parquet", index = False)
    #Total_MHDV.to_parquet(path_or_buf = path_result + "\Total_MHDV_" + Method + ".parquet", index = False)

    # Total_SPH_Sum = Total_SPH.groupby(["SCENARIO", "State", "Year"], as_index=False)["LoadMW"].sum()
    # Total_SPH_Sum = Total_SPH_Sum.rename(columns = {"LoadMW" : "TotalMW"})
    # Total_WH_Sum = Total_WH.groupby(["SCENARIO", "State", "Year"], as_index=False)["LoadMW"].sum()
    # Total_WH_Sum = Total_WH_Sum.rename(columns = {"LoadMW" : "TotalMW"})
    # Trans_LDV_Sum = Trans_LDV.groupby(["SCENARIO", "State", "Year"], as_index=False)["LDV_LoadMW"].sum()
    # Trans_LDV_Sum = Trans_LDV_Sum.rename(columns = {"LDV_LoadMW" : "TotalMW"})
    # Total_MHDV_Sum = Total_MHDV.groupby(["SCENARIO", "State", "Year"], as_index=False)["LoadMW"].sum()
    # Total_MHDV_Sum = Total_MHDV_Sum.rename(columns = {"LoadMW" : "TotalMW"})

    #Total_SPH_Sum.to_parquet(path_or_buf= path_result + "\Total_SPH_sum_" + Method + ".parquet", index = False)
    #Total_WH_Sum.to_parquet(path_or_buf= path_result + "\Total_WH_sum_" + Method + ".parquet", index = False)
    #Trans_LDV_Sum.to_parquet(path_or_buf= path_result + "\Trans_LDV_sum_" + Method + ".parquet", index = False)
    #Total_MHDV_Sum.to_parquet(path_or_buf= path_result + "\Total_MHDV_sum_" + Method + ".parquet", index = False)

    del  \
        Res_SPH_sum, Res_WH_sum, Com_SPH_sum, Trans_LDV_sum, \
             Trans_MDV_sum, Trans_HDV_sum, Trans_BUS_sum

    ################
    #Distribute Load to GenX.Region
    Method = 'Method1'
    subsectors = [Res_SPH, Res_WH, Com_SPH, Com_WH, Trans_LDV, Trans_MDV, Trans_HDV, Trans_BUS]
    subsector_names = ['Res_SPH', 'Res_WH', 'Com_SPH', 'Com_WH', 'Trans_LDV', 'Trans_MDV', 'Trans_HDV', 'Trans_BUS']
    column_names = ['Res_SPH', 'Res_WH', 'Com_SPH', 'Com_WH', 'LDV', 'MDV', 'HDV', 'BUS']
    j = 0

    for i in subsectors:
        temp = pd.merge(i, pop, on = ["State"], how = 'left')
        temp = temp.assign(weighted = temp[column_names[j] + "_" + "LoadMW"] * temp["State Prop"]).groupby(["SCENARIO", "Year", "LocalHourID", "GenX.Region"], as_index=False)["weighted"].sum()\
            .rename(columns={"weighted" : column_names[j] + "_MW"})
        temp.to_parquet(path_result + "\\" + subsector_names[j] + "_By_region.parquet", index = False)
        j = j + 1
    del temp, subsectors, Res_SPH, Res_WH, Com_SPH, Com_WH, Trans_LDV, Trans_MDV, Trans_HDV, Trans_BUS

    #Total_WH = pd.read_parquet(path_result + "\Total_WH_" + Method + ".parquet")
    # Total_WH = pd.merge(Total_WH, pop, on = ["State"], how = 'left')
    # Total_WH = Total_WH.assign(weighted = Total_WH["LoadMW"] * Total_WH["State Prop"]).groupby(["SCENARIO", "Year", "LocalHourID", "GenX.Region"], as_index=False)["weighted"].sum()\
    #     .rename(columns={"weighted" : "water_heat_MW"})
    # Total_WH = Total_WH[Total_WH["GenX.Region"].isin(regions)]
    # Total_WH.to_parquet(path_result + "\Total_WH_By_region.parquet", index = False)
    # del Total_WH

    # #Total_SPH = pd.read_parquet(path_result + "\Total_SPH_" + Method + ".parquet")
    # Total_SPH = pd.merge(Total_SPH, pop, on = ["State"], how = 'left')
    # Total_SPH = Total_SPH.assign(weighted = Total_SPH["LoadMW"] * Total_SPH["State Prop"]).groupby(["SCENARIO", "Year", "LocalHourID", "GenX.Region"], as_index=False)["weighted"].sum()\
    #     .rename(columns={"weighted" : "space_heat_MW"})

    # Total_SPH.to_parquet(path_result + "\Total_SPH_By_region.parquet", index = False)
    # del Total_SPH

    # #Total_LDEV = pd.read_parquet(path_result + "\Trans_LDV_" + Method + ".parquet")
    # Total_LDEV = Trans_LDV
    # Total_LDEV = pd.merge(Total_LDEV, pop, on = ["State"], how = 'left')
    # Total_LDEV = Total_LDEV.assign(weighted = Total_LDEV["LDV_LoadMW"] * Total_LDEV["State Prop"]).groupby(["SCENARIO", "Year", "LocalHourID", "GenX.Region"], as_index=False)["weighted"].sum()\
    #     .rename(columns={"weighted" : "LDEV_MW"})
    # Total_LDEV.to_parquet(path_result + "\Total_LDEV_By_region.parquet", index = False)
    # del Total_LDEV

    # #Total_MHDV = pd.read_parquet(path_result + "\Total_MHDV_" + Method + ".parquet")
    # Total_MHDV = pd.merge(Total_MHDV, pop, on = ["State"], how = 'left')
    # Total_MHDV = Total_MHDV.assign(weighted = Total_MHDV["LoadMW"] * Total_MHDV["State Prop"]).groupby(["SCENARIO", "Year", "LocalHourID", "GenX.Region"], as_index=False)["weighted"].sum()\
    #     .rename(columns={"weighted" : "MHBEV_MW"})
    # Total_MHDV.to_parquet(path_result + "\Total_MHDV_By_region.parquet", index = False)
    # del Total_MHDV

    ######
    #Construct Total Load

    Base_Load = pd.read_parquet(path_result + "\Base_Load.parquet")
    Base_Load = Base_Load.rename(columns= {"LoadMW" : "Base_MW"})

    Base_Load.loc[(Base_Load["Sector"] == "Commercial"), "Subsector"] = "Base_Com_other"
    Base_Load.loc[(Base_Load["Sector"] == "Residential"), "Subsector"] = "Base_Res_other"
    Base_Load.loc[(Base_Load["Sector"] == "Transportation"), "Subsector"] = "Base_Trans_other"
    Base_Load.loc[(Base_Load["Sector"] == "Industrial"), "Subsector"] = "Base_Ind"

    Base_Load = Base_Load.drop(columns = ["Sector"])
    Base_Load = Base_Load.pivot_table(index = [Base_Load.index.values, 'LocalHourID', "GenX.Region", "Year"], columns = 'Subsector', values = 'Base_MW')\
        .reset_index(['LocalHourID', "GenX.Region", "Year"]).fillna(0)
    Base_Load = Base_Load.groupby(['LocalHourID', 'GenX.Region', 'Year'], as_index=False).agg({'Base_Com_other' : 'sum', 'Base_Res_other' : 'sum', \
        'Base_Trans_other' : 'sum', 'Base_Ind' : 'sum'})

    Res_WH = pd.read_parquet(path_result + "\Res_WH_By_region.parquet")
    Com_WH = pd.read_parquet(path_result + "\Com_WH_By_region.parquet")
    Res_SPH = pd.read_parquet(path_result + "\Res_SPH_By_region.parquet")
    Com_SPH = pd.read_parquet(path_result + "\Com_SPH_By_region.parquet")
    Trans_LDV = pd.read_parquet(path_result + "\Trans_LDV_By_region.parquet")
    Trans_MDV = pd.read_parquet(path_result + "\Trans_MDV_By_region.parquet")
    Trans_HDV = pd.read_parquet(path_result + "\Trans_HDV_By_region.parquet")
    #Trans_BUS = pd.read_parquet(path_result + "\Trans_BUS_By_region.parquet")

    Total_Load = pd.merge(Res_WH, Com_WH, on = ["SCENARIO","Year","LocalHourID","GenX.Region"])
    Total_Load = pd.merge(Total_Load, Res_SPH, on = ["SCENARIO","Year","LocalHourID","GenX.Region"])
    Total_Load = pd.merge(Total_Load, Com_SPH, on = ["SCENARIO","Year","LocalHourID","GenX.Region"])
    Total_Load = pd.merge(Total_Load, Trans_LDV, on = ["SCENARIO","Year","LocalHourID","GenX.Region"])
    Total_Load = pd.merge(Total_Load, Trans_MDV, on = ["SCENARIO","Year","LocalHourID","GenX.Region"])
    Total_Load = pd.merge(Total_Load, Trans_HDV, on = ["SCENARIO","Year","LocalHourID","GenX.Region"])
    #Total_Load = pd.merge(Total_Load, Trans_BUS, on = ["SCENARIO","Year","LocalHourID","GenX.Region"])
    Total_Load = pd.merge(Total_Load, Base_Load, on = ["Year","LocalHourID","GenX.Region"])
    Total_Load = Total_Load[Total_Load["Year"].isin(years)]
    del Base_Load, Res_WH, Com_WH, Res_SPH, Com_SPH, Trans_LDV, Trans_MDV, Trans_HDV

    Total_Load = Total_Load[(Total_Load["Year"].isin(years)) & (Total_Load["GenX.Region"].isin(regions))]
    Total_Load.to_parquet(path_result + "\Total_load_by_region_full.parquet", index = False)

    #Total_Load_plot = Total_Load.assign(Total = Total_Load["Base_Com_SPH"] + Total_Load["Base_Res_SPH"] + Total_Load["Base_Res_WH"] + Total_Load["Base_Com_WH"] \
    #    + Total_Load["Base_Trans_LDV"] + Total_Load["Base_Trans_MDV"] + Total_Load["Base_Trans_HDV"] + Total_Load["Res_SPH_MW"] + Total_Load['Res_WH_MW']\
    #    + Total_Load['LDV_MW'] + Total_Load["MDV_MW"] + Total_Load["HDV_MW"] + Total_Load["Com_WH_MW"] + Total_Load["Com_SPH_MW"])\
    #        .melt(id_vars=["SCENARIO", "Year","LocalHourID","GenX.Region"])

    #Total_Load_sum = Total_Load_plot.groupby(["SCENARIO", "Year", "GenX.Region", "variable"], as_index = False)["value"].sum()
    #Total_Load_sum = Total_Load_sum.rename(columns={"value" : "AnnualTWh"})
    #Total_Load_sum["AnnualTWh"] = round(10**-6*Total_Load_sum["AnnualTWh"], 2)
    #del Total_Load_plot
    #Total_Load_sum_by_state = pd.merge(Total_Load_sum, pop, on = ["GenX.Region"], how = 'left')
    #Total_Load_sum_by_state = Total_Load_sum_by_state.assign(weighted = Total_Load_sum_by_state["AnnualTWh"] * Total_Load_sum_by_state["Zone Prop"])\
    #    .groupby(["SCENARIO", "State", "Year", "variable"])["weighted"].sum().rename(columns={"weighted" : "AnnualTWh"})
    #Total_Load_sum_by_state.to_parquet(path_result + "\Total_MWh_by_state.parquet", index = False)
    #del Total_Load_sum

    #Total_Load = Total_Load.assign(Total_MW = Total_Load["water_heat_MW"] + Total_Load["space_heat_MW"] + Total_Load["LDEV_MW"] + Total_Load["MHBEV_MW"] + Total_Load["Base_MW"])\
    #        .drop(columns = ["water_heat_MW", "space_heat_MW", "LDEV_MW", "MHBEV_MW", "Base_MW", "SCENARIO", "Year"]) \
    #        .pivot_table(index = 'time_index', columns = 'region', values = 'Total_MW')

    return Total_Load

def MakeLoadProfiles(settings, case_folder):
    #path_processed = r"C:\Users\ritib\Dropbox\Project_LoadConstruction\data\processed"
    #path_result = r"C:\Users\ritib\Dropbox\Project_LoadConstruction\data\result"
    CreateOutputFolder(case_folder)

    output_folder = case_folder / "extra_outputs"

    years = []
    regions = []
    electrification = []
    # scenarios = pd.DataFrame()
    for year in settings:
        for cases, _settings in settings[year].items():
            years.append(_settings["model_year"])
            regions = regions +  regions_to_keep(_settings)[0]
            electrification.append(_settings["NZA_electrification"])
            try:
                path_stock = str(_settings["input_folder"]) + "\\" + _settings["custom_stock"]
                path_growthrate = str(_settings["input_folder"]) + "\\" + _settings["custom_growthrate"]
            except:
                path_stock = ""
                path_growthrate = ""
            # scenarios
    #years = list(set(years))
    regions = list(set(regions))
    #electrification = list(set(electrification))
    CreateBaseLoad(years, regions, output_folder, path_growthrate)
    return AddElectrification(years, regions, electrification, output_folder, path_stock)

def FilterTotalProfile(settings, total_load):
    TotalLoad = total_load
    settings = settings
    TotalLoad = TotalLoad[TotalLoad["Year"] == settings["model_year"]]
    TotalLoad = TotalLoad.assign(TotalMW = TotalLoad['Res_WH_MW'] + TotalLoad['Com_WH_MW'] + TotalLoad['Res_SPH_MW'] + TotalLoad['Com_SPH_MW']\
        + TotalLoad['LDV_MW'] + TotalLoad['MDV_MW'] + TotalLoad['HDV_MW'] + TotalLoad['Base_Com_other']\
             + TotalLoad['Base_Res_other'] + TotalLoad['Base_Trans_other'] + TotalLoad['Base_Ind']).drop(columns = ['Res_WH_MW', 'Com_WH_MW', 'Res_SPH_MW', 'Com_SPH_MW', 'LDV_MW', 'MDV_MW', 'HDV_MW',\
                    'Base_Com_other', 'Base_Res_other', 'Base_Trans_other', 'Base_Ind'])
    TotalLoad = TotalLoad[TotalLoad["SCENARIO"] == settings["NZA_electrification"]]
    TotalLoad = TotalLoad.drop(columns = ["SCENARIO", "Year"]).rename(columns = {"LocalHourID" : "time_index", "GenX.Region" : "region"})

    keep_regions, region_agg_map = regions_to_keep(settings)

    TotalLoad.loc[:, "region_id_epaipm"] = TotalLoad.loc[:, "region"]
    TotalLoad.loc[
        TotalLoad.region_id_epaipm.isin(region_agg_map.keys()), "region"
    ] = TotalLoad.loc[
        TotalLoad.region_id_epaipm.isin(region_agg_map.keys()), "region_id_epaipm"
    ].map(
        region_agg_map
    )

    logger.info("Aggregating load curves in grouped regions")
    TotalLoad = TotalLoad.groupby(["region", "time_index"]).sum()

    TotalLoad = TotalLoad.unstack(level=0)
    TotalLoad.columns = TotalLoad.columns.droplevel()
    #TotalLoad = TotalLoad.pivot_table(index = 'time_index', columns = 'region', values = 'TotalMW')
    return TotalLoad

max_stars_repo_name: ritibhandarkar/PowerGenome

max_issues_repo_name: ritibhandarkar/PowerGenome

max_forks_repo_name: ritibhandarkar/PowerGenome


text:
"""
Target Problem:
---------------
* A classifier for the diagnosis of Autism Spectrum Disorder (ASD)

Proposed Solution (Machine Learning Pipeline):
----------------------------------------------
* Standard Scaling -> PCA -> Decision Tree Classifier

Input to Proposed Solution:
---------------------------
* Directories of training and testing data in csv file format
* These two types of data should be stored in n x m pattern in csv file format.

  Typical Example:
  ----------------
  n x m samples in training csv file (n number of samples, m - 1 number of features, ground truth labels at last column)
  k x s samples in testing csv file (k number of samples, s number of features)

* These data set files are ready by load_data() function.
* For comprehensive information about input format, please check the section
  "Data Sets and Usage Format of Source Codes" in README.md file on github.

Output of Proposed Solution:
----------------------------
* Predictions generated by learning model for testing set
* They are stored in "submission.csv" file.

Code Owner:
-----------
* Copyright © Team 4. All rights reserved.
* Copyright © Istanbul Technical University, Learning From Data Spring 2019. All rights reserved. """

import csv
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.tree import DecisionTreeClassifier

np.random.seed(1)  # Anchoring randomization during training step

def load_data():

    """
    The method reads train and test file to obtain train and test data.
    Then, it splits train data into two parts which are features and labels.
    :return: features and labels of train data and test data are returned.

    """

    train_data = pd.read_csv('train.csv')
    test_data = pd.read_csv('test.csv')

    all_train = train_data.iloc[:, :].values
    y_train = train_data.iloc[:, 595].values  # labels
    x_train = train_data.iloc[:, 0:595].values  # features

    x_test = test_data.iloc[:, 0:595].values

    return x_train, y_train, x_test

def standardization(x_train, x_test):

    """
    The method performs standard scaling on training and testing data.
    When doing this, only train data is included in training phase (fitting operation).

    Parameters
    ----------
    x_train: features of train data
    x_test: features of test data

    """

    sc = StandardScaler()
    x_train_std = sc.fit_transform(x_train)
    x_test_std = sc.transform(x_test)

    return x_train_std, x_test_std

def dim_red(x_train, x_test):

    """
    The method reduces the dimension of training and testing data by using PCA.
    When doing this, only training data is included in training phase (fitting operation).

    Parameters
    ----------
    x_train: features of scaled training data
    x_test: features of scaled testing data

    """

    pca = PCA(n_components=15)
    x_train_red = pca.fit_transform(x_train)
    x_test_red = pca.transform(x_test)

    return x_train_red, x_test_red

def decision_tree(criterion_name, x_train, y_train, x_test):

    """
    The method creates a decision tree learning model, trains it by using train data and generate predictions for
        testing data.
    Then, predicted labels are returned.

    Parameters
    ----------
    criterion_name: it specifies which function is used to measure quality of a split
    x_train: features of scaled and reduced training set
    y_train: labels of scaled and reduced training set
    x_test: features of scaled and reduced testing set

    """

    dtc = DecisionTreeClassifier(criterion=criterion_name, max_depth=4,
                                 max_features=3, random_state=3, max_leaf_nodes=2)
    dtc.fit(x_train, y_train)
    y_pred = dtc.predict(x_test)

    return y_pred

def write_output(y_pred):

    fields = ['ID', 'Predicted']
    filename = "submission.csv"
    rows = list()
    with open(filename, 'w',newline="") as csvfile:
        csvwriter = csv.writer(csvfile)
        csvwriter.writerow(fields)
        for i in range(len(y_pred)):
            rows.append([i+1, y_pred[i]])
        csvwriter.writerows(rows)

# ********** MAIN PROGRAM ********** #

x_train, y_train, x_test = load_data()
x_train_std, x_test_std = standardization(x_train, x_test)
x_train_red, x_test_red = dim_red(x_train_std, x_test_std)

y_pred = decision_tree('entropy', x_train_red, y_train, x_test_red)
write_output(y_pred)

ext: py

max_stars_repo_path: Team 4/main.py

max_stars_repo_name: basiralab/BrainNet-ML-ToolBox

max_issues_repo_path: Team 4/main.py

max_issues_repo_name: aycagokdag/BrainNet-ML-ToolBox

max_forks_repo_path: Team 4/main.py

max_forks_repo_name: aycagokdag/BrainNet-ML-ToolBox

reason: import numpy


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

from torch.autograd import Variable
from onnx import numpy_helper

import io
import onnx
import os
import shutil
import torch
import traceback

import test_pytorch_common
import test_onnx_common
from common_nn import module_tests
from test_nn import new_module_tests

# Take a test case (a dict) as input, return the test name.
def get_test_name(testcase):
    if "fullname" in testcase:
        return "test_" + testcase["fullname"]

    test_name = "test_" + testcase["constructor"].__name__
    if "desc" in testcase:
        test_name += "_" + testcase["desc"]
    return test_name

# Take a test case (a dict) as input, return the input for the module.
def gen_input(testcase):
    if "input_size" in testcase:
        return Variable(torch.randn(*testcase["input_size"]))
    elif "input_fn" in testcase:
        input = testcase["input_fn"]()
        if isinstance(input, Variable):
            return input
        return Variable(testcase["input_fn"]())

def gen_module(testcase):
    if "constructor_args" in testcase:
        args = testcase["constructor_args"]
        module = testcase["constructor"](*args)
        module.train(False)
        return module
    module = testcase["constructor"]()
    module.train(False)
    return module

def convert_tests(testcases, sets=1):
    print("Collect {} test cases from PyTorch.".format(len(testcases)))
    failed = 0
    ops = set()
    for t in testcases:
        test_name = get_test_name(t)
        module = gen_module(t)
        try:
            input = gen_input(t)
            f = io.BytesIO()
            torch.onnx._export(module, input, f)
            onnx_model = onnx.load_from_string(f.getvalue())
            onnx.checker.check_model(onnx_model)
            onnx.helper.strip_doc_string(onnx_model)
            output_dir = os.path.join(test_onnx_common.pytorch_converted_dir, test_name)

            if os.path.exists(output_dir):
                shutil.rmtree(output_dir)
            os.makedirs(output_dir)
            with open(os.path.join(output_dir, "model.onnx"), "wb") as file:
                file.write(onnx_model.SerializeToString())

            for i in range(sets):
                output = module(input)
                data_dir = os.path.join(output_dir, "test_data_set_{}".format(i))
                os.makedirs(data_dir)

                for index, var in enumerate([input]):
                    tensor = numpy_helper.from_array(var.data.numpy())
                    with open(os.path.join(data_dir, "input_{}.pb".format(index)), "wb") as file:
                        file.write(tensor.SerializeToString())
                for index, var in enumerate([output]):
                    tensor = numpy_helper.from_array(var.data.numpy())
                    with open(os.path.join(data_dir, "output_{}.pb".format(index)), "wb") as file:
                        file.write(tensor.SerializeToString())
                input = gen_input(t)
        except:
            traceback.print_exc()
            failed += 1

    print("Collect {} test cases from PyTorch repo, failed to export {} cases.".format(
        len(testcases), failed))
    print("PyTorch converted cases are stored in {}.".format(test_onnx_common.pytorch_converted_dir))

if __name__ == '__main__':
    testcases = module_tests + new_module_tests
    convert_tests(testcases)

max_stars_repo_name: fumihwh/onnx-fb-universe

max_issues_repo_name: fumihwh/onnx-fb-universe

max_forks_repo_name: fumihwh/onnx-fb-universe


text:
"""
This framework is based on SSD_tensorlow(https://github.com/balancap/SSD-Tensorflow)
"""
import math
import numpy as np
import tensorflow as tf
import tf_extended as tfe
from nets import custom_layers
from nets import textbox_common
from collections import namedtuple
import tensorflow.contrib.slim as slim

# =========================================================================== #
# Text class definition.
# =========================================================================== #
# Define a namedtuple data structure to keep the params of TextboxNet class.
TextboxParams = namedtuple('TextboxParameters',['img_shape',
												'num_classes',
												'feat_layers',
												'feat_shapes',
												'anchor_ratios',
												'anchor_sizes',
												'anchor_steps', #辅助放缩搜索网格中心点的位置
												'normalizations',
												'prior_scaling'
])

class TextboxNet(object):
	"""
	Implementation of the Textbox 384 network.
	The default features layers with 384x384 image input are:
	  conv4  ==> 48 x 48
	  conv7  ==> 24 x 24
	  conv8  ==> 12 x 12
	  conv9  ==> 6 x 6
	  conv10 ==> 4 x 4
	  conv11 ==> 2 x 2
	The default image size used to train this network is 384x384.
	Stage 1: 384x384， Stage 2: 768x768.
	"""
	default_params = TextboxParams(
		img_shape=(384, 384),
		num_classes=2,
		feat_layers=['conv4', 'conv7', 'conv8', 'conv9', 'conv10', 'conv11'], # 要进行检测的特征图name
		feat_shapes=[(48, 48), (24, 24), (12, 12), (6, 6), (4, 4), (2, 2)], # 特征图大小
		# Anchors 48 = 8 ratios * 6 sizes, each textbox layer 8 = 8 ratios * 1 size
		#TODO: Different anchor ratios setting from the original paper (1,2,3,5,1/2,1/3,1/5)->(2,1/2,3,1/3,4,1/4,5,1/5).
		anchor_ratios=[
			[2.0, 1. / 2, 3.0, 1. / 3, 4.0, 1. / 4, 5., 1. / 5],
			[2.0, 1. / 2, 3.0, 1. / 3, 4.0, 1. / 4, 5., 1. / 5],
			[2.0, 1. / 2, 3.0, 1. / 3, 4.0, 1. / 4, 5., 1. / 5],
			[2.0, 1. / 2, 3.0, 1. / 3, 4.0, 1. / 4, 5., 1. / 5],
			[2.0, 1. / 2, 3.0, 1. / 3, 4.0, 1. / 4, 5., 1. / 5],
			[2.0, 1. / 2, 3.0, 1. / 3, 4.0, 1. / 4, 5., 1. / 5],
		], # 特征图先验框所采用的长宽比
		anchor_sizes=[
			(30.,60.),
			(30.,90.),
			(90.,150.),
			(150., 210.),
			(210., 270.),
			(270., 330.)
		], # 不同特征图的先验框尺度（第一个值是s_k，第2个值是s_k+1）
		# 辅助放缩搜索网格中心点的位置
		anchor_steps=[8, 16, 32, 64, 100, 300], # 特征图单元大小 （image_shape / feature_map_shape），用来计算单元的center,
												# 即该层前所有卷积步长stride的乘积， 48*8 = 24*16 = 32*12 = 64*6 = 384 image_size
		normalizations=[20, -1, -1, -1, -1, -1], # L2 Normalization
		prior_scaling=[0.1, 0.1, 0.2, 0.2])  # offset calculate variances

	def __init__(self, params=None):
		"""
		Init the Textbox net with some parameters. Use the default ones
		if none provided.
		"""
		if isinstance(params, TextboxParams):
			self.params = params
		else:
			self.params = self.default_params

	# ======================================================================= #
	def net(self,
			inputs,
			is_training=True,
			dropout_keep_prob=0.5,
			reuse=None,
			scope='text_box_384',
			update_feat_shapes=False):
		"""
		Text network definition.
		"""
		# r = [predictions, localisations, logits, end_points, shape_list(if update_feat_shapes = True)]
		r = text_net(
			inputs,
			feat_layers=self.params.feat_layers,
			anchor_sizes=self.params.anchor_sizes,
			anchor_ratios=self.params.anchor_ratios,
			normalizations=self.params.normalizations,
			is_training=is_training,
			dropout_keep_prob=dropout_keep_prob,
			reuse=reuse,
			scope=scope,
			update_feat_shapes=update_feat_shapes)
		# Update feature shapes (try at least!)
		# feat_shape: [(48, 48), (24, 24), (12, 12), (6, 6), (4, 4), (2, 2)]
		if update_feat_shapes:
			#TODO: r[0]：各选中层预测结果，predictions of r[-1]: 各选中层的feature shape?
			# 获取各个中间层shape（不含0维），如果含有None则返回默认的feat_shapes
			shapes = textboxes_feat_shapes_from_net(r[-1],
													self.params.feat_shapes)
			self.params = self.params._replace(feat_shapes=shapes)
		print(self.default_params.feat_shapes)
		return r

	def arg_scope(self, weight_decay=0.0005, data_format='NHWC'):
		"""Get the defined Network arg_scope Encapsulation.
		"""
		return ssd_arg_scope(weight_decay, data_format=data_format)

	def arg_scope_caffe(self, caffe_scope):
		"""Get the defined Caffe arg_scope Encapsulation used for weights importing.
		"""
		return ssd_arg_scope_caffe(caffe_scope)

	# ======================================================================= #

	def anchors(self, img_shape, dtype=np.float32):
		"""Compute the default anchor boxes, given an image shape.
		"""
		return textbox_anchor_all_layers(
			img_shape,
			self.params.feat_shapes,
			self.params.anchor_ratios,
			self.params.anchor_sizes,
			self.params.anchor_steps,
			0.5, dtype)

	def bboxes_encode(self,
					  glabels,
					  bboxes,
					  anchors,
					  gxs,
					  gys,
					  scope='text_bboxes_encode'):
		"""Encode labels and bounding boxes.
		"""
		return textbox_common.tf_text_bboxes_encode(
			glabels,
			bboxes,
			anchors,
			gxs,
			gys,
			matching_threshold=0.5,
			prior_scaling=self.params.prior_scaling,
			dtype=tf.float32,
			scope=scope)

	def losses(self,
			   logits,
			   localisations,
			   glabels,
			   glocalisations,
			   gscores,
			   match_threshold=0.5,
			   negative_ratio=3.,
			   alpha=0.2,
			   label_smoothing=0.,
			   batch_size=16,
			   scope='txt_losses'):
		"""Define the SSD network losses.
		"""
		return ssd_losses(
			logits,
			localisations,
			glabels,
			glocalisations,
			gscores,
			match_threshold=match_threshold,
			negative_ratio=negative_ratio,
			alpha=alpha,
			label_smoothing=label_smoothing,
			batch_size=batch_size,
			scope=scope)

def textboxes_feat_shapes_from_net(predictions, default_shapes=None):
	"""Try to obtain the feature shapes from the prediction layers.
	Args:
	  predictions: a list of feature shapes return from text_net()
	Return:
	  list of feature shapes. Default values if predictions shape not fully determined.
	"""
	feat_shapes = []
	for l in predictions:
		print(l)
		shape = l.get_shape().as_list()[1:3]
		# shape = tuple(l[1:3])

		if None in shape:
			return default_shapes
		else:
			feat_shapes.append(shape)
	return feat_shapes

def text_multibox_layer(layer,
						inputs,
						anchor_size,
						anchor_ratio,
						normalization=-1):
	"""
	Construct a multibox layer, return a class and localization predictions.
	The most different between textbox and ssd is the prediction shape
	where textbox has prediction score shape (48,48,2,6) and location has shape (48,48,2,6,12)
	besise,the kernel for every layer same
	"""
	net = inputs
	if normalization > 0:
		net = custom_layers.l2_normalization(net, scaling=True)

	num_classes = 2
	# Number of anchors, 8+2.
	num_anchors = len(anchor_ratio) + len(anchor_size)
	# Vertical offsets.
	num_prior_per_location = num_anchors * 2

	# Location prediction.
	# Location 4+8 coordinates, [minimum horizontal bounding box offsets-4] + [quadrilateral bounding box offsets-8] = 12.
	loc_pred = slim.conv2d(net, num_prior_per_location*12, [3, 5], activation_fn=None,    # 3x5 conv used in textbox layer.
						   padding='SAME', scope='conv_loc')
	# TODO： different from ssd300
	# transform to NHWC
	#loc_pred = custom_layers.channel_to_last(loc_pred)
	loc_pred = slim.flatten(loc_pred)
	l_shape = loc_pred.shape
	batch_size = l_shape[0]
	loc_pred = tf.reshape(loc_pred, [batch_size, -1, 12])
	# NHW(num_anchors+4) --> [NHW, num_anchors, 4]
	#loc_pred = tf.reshape(loc_pred, tensor_shape(loc_pred, 4)[:-1]+[num_prior_per_location, 12])

	# Class prediction.
	cls_pred = slim.conv2d(net, num_prior_per_location*num_classes, [3, 5], activation_fn=None,
						   padding='SAME',scope='conv_cls')
	# transform to NHWC
	#cls_pred = custom_layers.channel_to_last(cls_pred)
	l_shape = tf.shape(cls_pred)
	cls_pred = slim.flatten(cls_pred)
	cls_pred = tf.reshape(cls_pred, [batch_size, -1 ,2])
	# NHW(num_anchors+classes) --> [NHW, num_anchors, classes]
	#cls_pred = tf.reshape(cls_pred, tensor_shape(cls_pred, 4)[:-1]+[num_prior_per_location, num_classes])

	return cls_pred, loc_pred, l_shape

def tensor_shape(x, rank=3):
	"""Returns the dimensions of a tensor.
	Args:
	  x: A N-D Tensor of shape.
	  rank:
	Returns:
	  A list of dimensions. Dimensions that are statically known are python
		integers,otherwise they are integer scalar tensors.
	"""
	if x.get_shape().is_fully_defined():
		return x.get_shape().as_list()
	else:
		static_shape = x.get_shape().with_rank(rank).as_list()
		dynamic_shape = tf.unstack(tf.shape(x), rank)
		return [s if s is not None else d
				for s, d in zip(static_shape, dynamic_shape)]

def text_net(inputs,
			 feat_layers=TextboxNet.default_params.feat_layers,
			 anchor_sizes=TextboxNet.default_params.anchor_sizes,
			 anchor_ratios = TextboxNet.default_params.anchor_ratios,
			 normalizations=TextboxNet.default_params.normalizations,
			 is_training=True,
			 dropout_keep_prob=0.5,
			 reuse=None,
			 scope='text_box_384',
			 update_feat_shapes=False):
	"""
	Define the backbone (13 original vgg layers + 10 extra conv layers to extract multi-scale feature maps form as SSD)
	of the textboxes and the neck net -- 6 textbox layers.
	:param inputs: input image size
	:param feat_layers: feature map which connect to the textbox layer
	:param anchor_sizes: multi-scale anchor sizes
	:param anchor_ratios: multi-scale anchor aspect ratios [2.0, 1. / 2, 3.0, 1. / 3, 4.0, 1. / 4, 5., 1. / 5]
	:param normalizations:
	:param is_training: train or not
	:param dropout_keep_prob:
	:param reuse:
	:param scope:
	:param update_feat_shapes:
	:return: [predictions, localisations, logits, end_points, shape_list(if update_feat_shapes=True)]
	"""
	# End_points collect relevant activations for external use.
	end_points = {}
	with tf.variable_scope(scope, 'text_box_300', [inputs], reuse=reuse):  # 300*300 384*384
		######################################
		# 前五个 Blocks，首先照搬 VGG16 架构   #
		# 注意这里使用 end_points 标注中间结果  #
		######################################
		# ——————————————————Original VGG-16 blocks (total 13 conv layers)———————————————————————
		# Block 1.
		net = slim.repeat(inputs, 2, slim.conv2d, 64, [3, 3], scope='conv1')  # 300 384
		end_points['conv1'] = net
		net = slim.max_pool2d(net, [2, 2], scope='pool1')  # 150
		# Block 2.
		net = slim.repeat(net, 2, slim.conv2d, 128, [3, 3], scope='conv2')  # 150 192
		end_points['conv2'] = net
		net = slim.max_pool2d(net, [2, 2], scope='pool2')  # 75
		# Block 3.
		net = slim.repeat(net, 3, slim.conv2d, 256, [3, 3], scope='conv3')  # 75 81
		end_points['conv3'] = net
		net = slim.max_pool2d(net, [2, 2], scope='pool3')  # 38
		# Block 4.
		end_point = 'conv4'
		net = slim.repeat(net, 3, slim.conv2d, 512, [3, 3], scope='conv4')  # 38 40
		end_points[end_point] = net
		net = slim.max_pool2d(net, [2, 2], scope='pool4')  # 19
		# Block 5.
		net = slim.repeat(net, 3, slim.conv2d, 512, [3, 3], scope='conv5')  # 19
		end_points['conv5'] = net
		net = slim.max_pool2d(net, [3, 3], stride=1, scope='pool5')  # 19  Pooling size 2 -> 3

		####################################
		# 后六个 Blocks，使用额外卷积层      #
		####################################
		# ————————————Additional SSD blocks.——————————————————————
		# Block 6: let's dilate the hell out of it!  dilation -> 6
		net = slim.conv2d(net, 1024, [3, 3], rate=6, scope='conv6')  # 19
		end_points['conv6'] = net
		# Block 7: 1x1 conv. Because the fuck.
		end_point = 'conv7'
		net = slim.conv2d(net, 1024, [1, 1], scope='conv7')  # 19
		end_points[end_point] = net

		# Block 8/9/10/11: 1x1 and 3x3 convolutions stride 2 (except lasts).
		end_point = 'conv8'
		with tf.variable_scope(end_point):
			net = slim.conv2d(net, 256, [1, 1], scope='conv1x1')
			net = custom_layers.pad2d(net, pad=(1, 1))
			net = slim.conv2d(net, 512, [3, 3], stride=2, scope='conv3x3', padding='VALID')
		end_points[end_point] = net  # 10

		end_point = 'conv9'
		with tf.variable_scope(end_point):
			net = slim.conv2d(net, 128, [1, 1], scope='conv1x1')
			net = custom_layers.pad2d(net, pad=(1, 1))
			net = slim.conv2d(net, 256, [3, 3], stride=2, scope='conv3x3', padding='VALID')
		end_points[end_point] = net # 5

		end_point = 'conv10'
		with tf.variable_scope(end_point):
			net = slim.conv2d(net, 128, [1, 1], scope='conv1x1')
			net = slim.conv2d(net, 256, [3, 3], scope='conv3x3', padding='VALID')
		end_points[end_point] = net  # 3

		end_point = 'conv11'
		with tf.variable_scope(end_point):
			net = slim.conv2d(net, 128, [1, 1], scope='conv1x1')
			net = slim.conv2d(net, 256, [3, 3], scope='conv3x3', padding='VALID')
		end_points[end_point] = net  #

		######################################
		# 每个中间层 end_points 返回中间结果   #
		# 将各层预测结果存入列表，返回给优化函数 #
		######################################
		# Prediction and localisations layers.
		predictions = []
		logits = []
		localisations = []
		shape_list = []
		# feat_layers=['conv4', 'conv7', 'conv8', 'conv9', 'conv10', 'conv11']
		for i, layer in enumerate(feat_layers):
			with tf.variable_scope(layer + '_box'):
				cls, loc, shape = text_multibox_layer(layer,
											 end_points[layer],
											 anchor_sizes[i],
											 anchor_ratios[i],
											 normalizations[i])
			prediction_fn = slim.softmax
			# Prediction of conference and bbox location of each textbox layer.
			predictions.append(prediction_fn(cls))
			logits.append(cls)
			localisations.append(loc)
			shape_list.append(shape)

		if update_feat_shapes is True:
			return predictions, localisations, logits, end_points, shape_list
		else:
			return predictions, localisations, logits, end_points

def textbox_anchor_one_layer(img_shape, feat_shape, ratios, size, step, offset=0.5, dtype=np.float32):
	"""
		Computer TextBoxes++ default anchor boxes for one selected feature layer.
		Each feature point has 12 default textboxes (6 default boxes + 6 offsets vertical boxes).
		Determine the relative position grid of the centers and the relative width and height.

	:param: img_shape: Image shape, used for computing height, width relatively to the former;
	:param: feat_shape: Feature shape of each layer which connect to textbox layer, used for computing relative position grids;
	:param: ratios: Anchor Aspect ratios to use on these features;
	:param: size: Anchor size;
	:param: step:
	:param:offset: Anchor vertical offset.

	:return: y, x, h, w: Relative x and y grids, and height and width of a anchor associated to a textbox layer.
	"""
	# Follow the papers scheme
	# 12 ahchors boxes with out sk' = sqrt(sk * sk+1), 生成feat_shape中HW对应的网格坐标
	y, x = np.mgrid[0:feat_shape[0], 0:feat_shape[1]] + offset

	# vertical offset, step*feat_shape 约等于img_shape，这使得网格点坐标介于0~1，放缩一下即可到图像大小
	y = y.astype(dtype) * step / img_shape[0]
	x = x.astype(dtype) * step / img_shape[1]
	y_offset = (y.astype(dtype) + offset) * step / img_shape[0]
	x_offset = x

	# (38, 38, 2) origin anchor + offset anchor
	x_out = np.stack((x, x_offset), -1)
	y_out = np.stack((y, y_offset), -1)

	# Expand dims to support easy broadcasting
	y_out = np.expand_dims(y_out, axis=-1)
	x_out = np.expand_dims(x_out, axis=-1)

	# Compute relative height and width, 8 ratios + 2 sizes.
	num_anchors = len(ratios) + len(size)
	# shape: (10,)
	h = np.zeros((num_anchors,), dtype=dtype)
	w = np.zeros((num_anchors,), dtype=dtype)

	# Add first anchor boxes with ratio=1 (smaller square).
	# e.g., h[0] = 30 / 300, w[0] = 30 / 300.
	h[0] = size[0] / img_shape[0]
	w[0] = size[0] / img_shape[1]
	di = 1

	if len(size) > 1:
		# Add last anchor boxes with ratio=1' (bigger square).
		# e.g., h[1] = sqrt(30 * 60) / 300, w[1] = sqrt(30 * 60) / 300.
		h[1] = math.sqrt(size[0] * size[1]) / img_shape[0]
		w[1] = math.sqrt(size[0] * size[1]) / img_shape[1]
		di += 1

	for i, r in enumerate(ratios):
		# Add the other anchors.
		# h[k] = Sk / sqrt(ar), w[k] = sk * sqrt(ar). e.g., h[2] = 30 / 300 / sqrt(2.0), w[2] = 30 / 300 * sqrt(2.0).
		h[i+di] = size[0] / img_shape[0] / math.sqrt(r)
		w[i+di] = size[0] / img_shape[1] * math.sqrt(r)

	xmin = x_out - w/2
	ymin = y_out - h/2
	xmax = x_out + w/2
	ymax = y_out + h/2

	xmin = xmin.reshape([xmin.shape[0], xmin.shape[1], -1], order='F').reshape(-1)
	ymin = ymin.reshape([ymin.shape[0], ymin.shape[1], -1], order='F').reshape(-1)
	xmax = xmax.reshape([xmax.shape[0], xmax.shape[1], -1], order='F').reshape(-1)
	ymax = ymax.reshape([ymax.shape[0], ymax.shape[1], -1], order='F').reshape(-1)

	return xmin, ymin, xmax, ymax

## produce anchor for all layers
def ssd_arg_scope(weight_decay=0.0005, data_format='NHWC'):
	"""Defines the VGG arg scope, 约束网络中的超参数设定.

	Args:
	  weight_decay: The l2 regularization coefficient.

	Returns:
	  An arg_scope.
	"""
	with slim.arg_scope(
		[slim.conv2d, slim.fully_connected],
			activation_fn=tf.nn.relu,
			weights_regularizer=slim.l2_regularizer(weight_decay),
			weights_initializer=tf.contrib.layers.xavier_initializer(),
			biases_initializer=tf.zeros_initializer()):
		with slim.arg_scope(
			[slim.conv2d, slim.max_pool2d],
				padding='SAME',
				data_format=data_format):
			with slim.arg_scope(
				[custom_layers.pad2d,
				 custom_layers.l2_normalization,
				 custom_layers.channel_to_last],
				data_format=data_format) as sc:
				return sc

def textbox_anchor_all_layers(img_shape, layers_shape, anchor_ratios, anchor_sizes, anchor_steps,
							  offset=0.5, dtype=np.float32):
	"""
		Compute anchor boxes for all feature layers.
	:param img_shape:
	:param layers_shape:
	:param anchor_ratios:
	:param anchor_sizes:
	:param anchor_steps:
	:param offset: vertical offset which used to cover all text and
			makes the anchors dense in the vertical orientation.
	:param dtype:
	:return: the anchors list of each layer
	"""
	layers_anchors = []
	# layers_shape: [(48, 48), (24, 24), (12, 12), (6, 6), (4, 4), (2, 2)]
	for i, s in enumerate(layers_shape):
		anchor_bboxes = textbox_anchor_one_layer(img_shape, s,
												 anchor_ratios[i],
												 anchor_sizes[i],
												 anchor_steps[i],
												 offset=offset, dtype=dtype)
		layers_anchors.append(anchor_bboxes)
	return layers_anchors

# =========================================================================== #
# Caffe scope: importing weights at initialization.
# =========================================================================== #
def ssd_arg_scope_caffe(caffe_scope):
	"""Caffe scope definition.

	Args:
	  caffe_scope: Caffe scope object with loaded weights.

	Returns:
	  An arg_scope.
	"""
	# Default network arg scope.
	with slim.arg_scope(
		[slim.conv2d],
			activation_fn=tf.nn.relu,
			weights_initializer=caffe_scope.conv_weights_init(),
			biases_initializer=caffe_scope.conv_biases_init()):
		with slim.arg_scope([slim.fully_connected], activation_fn=tf.nn.relu):
			with slim.arg_scope(
				[custom_layers.l2_normalization],
					scale_initializer=caffe_scope.l2_norm_scale_init()):
				with slim.arg_scope(
					[slim.conv2d, slim.max_pool2d], padding='SAME') as sc:
					return sc

# =========================================================================== #
# Text loss function.
# =========================================================================== #
def ssd_losses(logits, localisations,
			   glabels, glocalisations, gscores,
			   match_threshold=0.5,
			   negative_ratio=3.,
			   alpha=0.2,
			   label_smoothing=0.,
			   batch_size=16,
			   scope=None):
	"""
	Loss functions for training the textboxes++ network.
	This function defines the different loss components of the SSD, and adds them to the TF loss collection.

	Arguments:
	  logits: (list of) predictions logits Tensors;                x
	  localisations: (list of) predictions localisations Tensors;  l
	  glocalisations: (list of) groundtruth localisations Tensors; g
	  gscores: (list of) groundtruth score / conference Tensors;   c
	"""
	# from ssd loss
	with tf.name_scope(scope, 'txt_losses'):
		lshape = tfe.get_shape(logits[0], 5)
		num_classes = lshape[-1]
		batch_size = batch_size

		l_cross_pos = []
		l_cross_neg = []
		l_loc = []

		# Flatten out all vectors!
		flogits = []
		fgscores = []
		flocalisations = []
		fglocalisations = []
		fglabels = []

		# Transform all matrix to vector.
		for i in range(len(logits)):
			flogits.append(tf.reshape(logits[i], [-1, num_classes]))
			fgscores.append(tf.reshape(gscores[i], [-1]))
			fglabels.append(tf.reshape(glabels[i], [-1]))
			flocalisations.append(tf.reshape(localisations[i], [-1, 12]))
			fglocalisations.append(tf.reshape(glocalisations[i], [-1, 12]))

		# And concat the crap!
		glabels = tf.concat(fglabels, axis=0)
		logits = tf.concat(flogits, axis=0)                      # x
		gscores = tf.concat(fgscores, axis=0)                   # c
		localisations = tf.concat(flocalisations, axis=0)      # l
		glocalisations = tf.concat(fglocalisations, axis=0)   # g
		dtype = logits.dtype

		# Compute positive matching mask...
		pmask = gscores > match_threshold  # positive mask
		# pmask = tf.concat(axis=0, values=[pmask[:tf.argmax(gscores, axis=0)], [True], pmask[tf.argmax(gscores, axis=0) + 1:]])
		ipmask = tf.cast(pmask, tf.int32)  # int positive mask
		fpmask = tf.cast(pmask, dtype)  # float positive mask
		n_positives = tf.reduce_sum(fpmask)  # calculate all number

		# Hard negative mining...
		no_classes = tf.cast(pmask, tf.int32)
		predictions = slim.softmax(logits)  #
		nmask = tf.logical_and(tf.logical_not(pmask), gscores > -0.5)  #
		fnmask = tf.cast(nmask, dtype)
		nvalues = tf.where(nmask, predictions[:, 0], 1. - fnmask)
		nvalues_flat = tf.reshape(nvalues, [-1])
		# Number of negative entries to select.
		max_neg_entries = tf.cast(tf.reduce_sum(fnmask), tf.int32)
		n_neg = tf.cast(negative_ratio * n_positives, tf.int32) + batch_size
		n_neg = tf.minimum(n_neg, max_neg_entries)

		val, idxes = tf.nn.top_k(-nvalues_flat, k=n_neg)
		max_hard_pred = -val[-1]
		# Final negative mask.
		nmask = tf.logical_and(nmask, nvalues < max_hard_pred)
		fnmask = tf.cast(nmask, dtype)

		# Add cross-entropy loss.
		# logits [batch_size, num_classes] labels [batch_size] ~ 0,num_class
		with tf.name_scope('cross_entropy_pos'):
			loss = tf.nn.sparse_softmax_cross_entropy_with_logits(
				logits=logits, labels=glabels)
			loss = tf.div(
				tf.reduce_sum(loss * fpmask), batch_size, name='value')
			tf.losses.add_loss(loss)
			l_cross_pos.append(loss)

		with tf.name_scope('cross_entropy_neg'):
			loss = tf.nn.sparse_softmax_cross_entropy_with_logits(
				logits=logits, labels=no_classes)
			loss = tf.div(
				tf.reduce_sum(loss * fnmask), batch_size, name='value')
			tf.losses.add_loss(loss)
			l_cross_neg.append(loss)

		# Add localization loss: smooth L1, L2, ...
		with tf.name_scope('localization'):
			# Weights Tensor: positive mask + random negative.
			weights = tf.expand_dims(alpha * fpmask, axis=-1)
			# localisations = tf.Print(localisations, [localisations, tf.shape(localisations)], "pre is:         ", summarize=20)
			# glocalisations = tf.Print(glocalisations, [glocalisations,  tf.shape(glocalisations)], "gt is :         ",summarize=20)
			loss = custom_layers.abs_smooth(localisations - glocalisations)
			loss = tf.div(
				tf.reduce_sum(loss * weights), batch_size, name='value')
			tf.losses.add_loss(loss)
			l_loc.append(loss)


text:
import numpy as np
import matplotlib.pyplot as plt
from .PlotData import PlotData

def PlotChain(Stations,Date,ut=None,fig=None,comp=['Bx','By','Bz','Bm'],high=None,low=None,useytitle=False,nolegend=False,figsize=(6,4),ylim=None,coords='hdz',RemoveMean=True):
	'''
	'''
	#create the plot window if it doesn't already exist
	if fig is None:
		fig = plt
		fig.figure(figsize=figsize)

	#now work out the total number of stations to plot
	ns = np.size(Stations)

	#loop through each plot
	for i in range(0,ns):
		nox = i < (ns-1)
		ax = PlotData(Stations[i],Date,ut=ut,fig=fig,maps=[1,ns,0,i],comp=comp,high=high,low=low,nox=nox,useytitle=useytitle,nolegend=nolegend,coords=coords,RemoveMean=RemoveMean)
		if not ylim is None:
			ax.set_ylim(ylim)

	#remove space between subplots
	fig.subplots_adjust(hspace=0.0)


text:
import numpy as np
import cv2
import glob

# termination criteria
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)

# prepare object points, like (0,0,0), (1,0,0), (2,0,0) ....,(6,5,0)
objp = np.zeros((6*7,3), np.float32)
objp[:,:2] = np.mgrid[0:7,0:6].T.reshape(-1,2)

# Arrays to store object points and image points from all the images.
objpoints = [] # 3d point in real world space
imgpoints = [] # 2d points in image plane.

images = glob.glob('calib_images/*.jpg')

for fname in images:
    img = cv2.imread(fname)
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

    # Find the chess board corners
    ret, corners = cv2.findChessboardCorners(gray, (7,6),None)

    # If found, add object points, image points (after refining them)
    if ret == True:
        objpoints.append(objp)

        corners2 = cv2.cornerSubPix(gray,corners,(11,11),(-1,-1),criteria)
        imgpoints.append(corners2)

        # Draw and display the corners
        img = cv2.drawChessboardCorners(img, (7,6), corners2,ret)

ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1],None,None)

# Function to draw the axis
def draw(img, corners, imgpts):
    corner = tuple(corners[0].ravel())
    img = cv2.line(img, corner, tuple(imgpts[0].ravel()), (255, 0, 0), 5)
    img = cv2.line(img, corner, tuple(imgpts[1].ravel()), (0, 255, 0), 5)
    img = cv2.line(img, corner, tuple(imgpts[2].ravel()), (0, 0, 255), 5)
    return img

criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
objp = np.zeros((6*7,3), np.float32)
objp[:,:2] = np.mgrid[0:7,0:6].T.reshape(-1,2)

axis = np.float32([[3,0,0], [0,3,0], [0,0,-3]]).reshape(-1,3)

for fname in glob.glob('calib_images/left03.jpg'):
    img = cv2.imread(fname)
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    ret, corners = cv2.findChessboardCorners(gray, (7,6),None)

    if ret == True:
        corners2 = cv2.cornerSubPix(gray,corners,(11,11),(-1,-1),criteria)

        # Find the rotation and translation vectors.
        _,rvecs, tvecs, inliers = cv2.solvePnPRansac(objp, corners2, mtx, dist)
        # project 3D points to image plane
        imgpts, jac = cv2.projectPoints(axis, rvecs, tvecs, mtx, dist)

        img = draw(img,corners2,imgpts)
        cv2.imshow('img',img)
        k = cv2.waitKey(0) & 0xff

cv2.destroyAllWindows()

max_stars_repo_name: RishabhBrajabasi/laserSmart


text:
# Runs after normalization and per_person_ratio_and_factor and pre_plot_aggregation.
import shutil
from pathlib import Path
import itertools
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import collections
from scipy.optimize import minimize_scalar

cols1 = ['F1_' + str(i) for i in range(3, 20, 2)]
cols2 = ['F2_' + str(i) for i in range(3, 20, 2)]
kCols = cols1 + cols2

def PlotWithSlices(df, data_name, output_dir):
    for group_name in ['Gender', 'AgeGroup', 'Family1', 'Family2', 'Family3', 'Family4', 'Education1', 'Career1', 'Career2', 'Language1', 'Word']:
        grouped_df = df.groupby([group_name])[kCols].mean()
        # grouped_df.to_csv(output_dir / (data_name + '_' + group_name + '_raw.csv'), index=True)
        full_group_name = '@'.join([data_name, group_name])
        fig = plt.figure()
        ax = fig.add_subplot(111, projection='3d')
        print(group_name)
        z_label = grouped_df.index.to_numpy().tolist()
        print(z_label)
        cmap = plt.get_cmap('viridis')
        colors = cmap(np.linspace(0, 1, len(z_label)))
        for key in z_label:
            x = np.arange(0, 9)
            color = colors[z_label.index(key)]
            z = z_label.index(key)
            mdf = grouped_df.loc[key]
            y1 = mdf[cols1].to_numpy(dtype='float')
            y2 = mdf[cols2].to_numpy(dtype='float')
            coeff1 = np.polyfit(x, y1, 4)
            coeff2 = np.polyfit(x, y2, 4)
            line1 = np.poly1d(coeff1)
            line2 = np.poly1d(coeff2)
            line1dd = np.polyder(line1, 2)
            line2dd = np.polyder(line2, 2)
            line1dd_max = minimize_scalar(-line1dd,
                                          bounds=(0, 8), method='bounded')
            line2dd_max = minimize_scalar(-line2dd,
                                          bounds=(0, 8), method='bounded')
            inflection1 = line1dd_max.x
            inflection2 = line2dd_max.x
            inflection1y = line1(inflection1)
            inflection2y = line2(inflection2)
            ax.plot(x, y1, zs=z, zdir='x', c=color, label='F1', linewidth=3.0)
            ax.plot(x, y2, zs=z, zdir='x', c=color, label='F2')
            ax.plot([inflection1, inflection1], [inflection1y-100, inflection1y+100], zs=z, zdir='x', c='black')
            ax.plot([inflection2, inflection2], [inflection2y-100, inflection2y+100], zs=z, zdir='x', c='black')

        ax.set(xticks=range(len(z_label)), xticklabels=z_label)
        plt.title(full_group_name)
        plt.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
        plt.savefig(output_dir / (full_group_name + '.png'),
                    bbox_inches="tight")
        plt.clf()
        plt.cla()

def PlotNoSlice(df, full_group_name, output_dir):
    x = np.arange(0, 9)
    y1 = df[cols1].to_numpy(dtype='float')
    y2 = df[cols2].to_numpy(dtype='float')
    coeff1 = np.polyfit(x, y1, 4)
    coeff2 = np.polyfit(x, y2, 4)
    line1 = np.poly1d(coeff1)
    line2 = np.poly1d(coeff2)
    line1dd = np.polyder(line1, 2)
    line2dd = np.polyder(line2, 2)
    line1dd_max = minimize_scalar(-line1dd,
                                  bounds=(0, 8), method='bounded')
    line2dd_max = minimize_scalar(-line2dd,
                                  bounds=(0, 8), method='bounded')
    inflection1 = line1dd_max.x
    inflection2 = line2dd_max.x

    # Plot f1/f2
    plt.plot(x, y1, 'o')
    plt.plot(x, y2, 'x')
    plt.plot(x, line1(x), label='F1 fitted')
    plt.plot(x, line2(x), label='F2 fitted')
    plt.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
    plt.title(full_group_name)
    plt.savefig(output_dir / (full_group_name + '.fitted.png'),
                bbox_inches="tight")
    plt.clf()
    plt.cla()
    # Plot deriv and inflection
    plt.plot(x, line1dd(x), label='F1 2nd deriv')
    plt.plot(x, line2dd(x), label='F2 2nd deriv')
    plt.axvline(x=inflection1, linestyle=':', label='F1 break')
    plt.axvline(x=inflection2, linestyle='-.', label='F2 break')
    plt.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
    plt.title(full_group_name)
    plt.savefig(output_dir / (full_group_name + '.break.png'),
                bbox_inches="tight")
    plt.clf()
    plt.cla()

def SlicePlotDataS(df, output_dir):
    cols1 = ['F1_' + str(i) for i in range(3, 20, 2)]
    cols2 = ['F2_' + str(i) for i in range(3, 20, 2)]
    kCols = cols1 + cols2
    matched_rows = []
    sa_a1_sb_a1 = df[df['IsSba2'] == 'No']
    #sa_a1_sb_a1.to_csv(output_dir / 'sa_a1_sb_a1_raw.csv', index=False)
    sa_a1_sb_a1_mean = sa_a1_sb_a1.groupby(['Pos'])[kCols].mean()
    #sa_a1_sb_a1_mean.to_csv(output_dir / 'sa_a1_sb_a1_mean.csv', index=True)
    PlotNoSlice(sa_a1_sb_a1_mean.iloc[0], 'sa_a1_sb_a1_a', output_dir)
    PlotNoSlice(sa_a1_sb_a1_mean.iloc[1], 'sa_a1_sb_a1_b', output_dir)
    sa_a1_sb_a2 = df[df['IsSba2'] == 'Yes']
    #sa_a1_sb_a2.to_csv(output_dir / 'sa_a1_sb_a2_raw.csv', index=False)
    sa_a1_sb_a2_mean = sa_a1_sb_a2.groupby(['Pos'])[kCols].mean()
    #sa_a1_sb_a2_mean.to_csv(output_dir / 'sa_a1_sb_a2_mean.csv', index=True)
    PlotNoSlice(sa_a1_sb_a2_mean.iloc[0], 'sa_a1_sb_a2_a', output_dir)
    PlotNoSlice(sa_a1_sb_a2_mean.iloc[1], 'sa_a1_sb_a2_b', output_dir)

    matched_rows = []
    for _, row in df.iterrows():
        comps = row['Filename'].split('_')
        lang = comps[0]
        pos = comps[4]
        if lang == 'S' and pos == 'b' and row['Annotation'] == 'a2':
            matched_rows.append(row)
    input_df = pd.DataFrame(matched_rows)
    PlotWithSlices(input_df, 'all_s_sb_a2', output_dir)

def SlicePlotDataM(df, output_dir):
    m_sb_a1 = df[df['IsSba2'] == 'No']
    PlotWithSlices(m_sb_a1, 'm_sb_a1', output_dir)
    m_sb_a1_mean = m_sb_a1.groupby(['IsSba2'])[kCols].mean()
    PlotNoSlice(m_sb_a1_mean.iloc[0], 'm_sb_a1', output_dir)

    m_sb_a2 = df[df['IsSba2'] == 'Yes']
    PlotWithSlices(m_sb_a2, 'm_sb_a2', output_dir)
    m_sb_a2_mean = m_sb_a2.groupby(['IsSba2'])[kCols].mean()
    PlotNoSlice(m_sb_a2_mean.iloc[0], 'm_sb_a2', output_dir)

input_base_dir = Path('./analysis/output/')
output_base_dir = Path('./analysis/output/break/')
shutil.rmtree(output_base_dir, ignore_errors=True)
output_base_dir.mkdir(parents=True, exist_ok=True)

df = pd.read_csv(input_base_dir / 'S_all_plot_raw_data.csv')
SlicePlotDataS(df, output_base_dir)

df = pd.read_csv(input_base_dir / 'M_all_plot_raw_data.csv')
SlicePlotDataM(df, output_base_dir)

max_stars_repo_name: weishi/praat

max_issues_repo_name: weishi/praat

max_forks_repo_name: weishi/praat

reason: import numpy,from scipy


text:
import matplotlib.pyplot as plt
import numpy as np
import nengo
import nengo_gui
nengo.log(level='info')

from run_nengo import build_target_layer, load_network

import base64
import PIL
import cStringIO

backend = 'nengo_ocl'

loadfile = 'checkpoints/ilsvrc2012-lif-48'
layers, data, dp = load_network(loadfile)

# image = data[0][0]
# image = np.transpose(image, (1, 2, 0))
# print(image.min(), image.max())

# plt.imshow((image / 255. + 0.5).clip(0, 1))
# plt.show()

# --- build network in Nengo
model = nengo.Network()

presentation_time = 0.2

model.config[nengo.Connection].synapse = nengo.synapses.Lowpass(0.0)
# model.config[nengo.Connection].synapse = nengo.synapses.Alpha(0.003)

outputs = build_target_layer(
    'logprob', layers, data, model, pt=presentation_time)

label_names = dp.batch_meta['label_names']
vocab_names = [
    (name.split(',')[0] if ',' in name else name).upper().replace(' ', '_')
    for i, name in enumerate(label_names)]
vocab_vectors = np.eye(len(vocab_names))

# number duplicates
unique = set()
duplicates = []
for name in vocab_names:
    if name in unique:
        duplicates.append(name)
    else:
        unique.add(name)

duplicates = {name: 0 for name in duplicates}
for i, name in enumerate(vocab_names):
    if name in duplicates:
        vocab_names[i] = '%s%d' % (name, duplicates[name])
        duplicates[name] += 1

with model:
    vocab = nengo.spa.Vocabulary(len(vocab_names))
    for name, vector in zip(vocab_names, vocab_vectors):
        vocab.add(name, vector)

    config = nengo.Config(nengo.Ensemble)
    config[nengo.Ensemble].neuron_type = nengo.Direct()
    with config:
        output = nengo.spa.State(len(vocab_names), subdimensions=10, vocab=vocab)
    nengo.Connection(outputs['probs'], output.input, transform=1.0, synapse=nengo.synapses.Alpha(0.01))

    # --- image display
    input_image = outputs['data']
    input_shape = data[0].shape[1:]

    def display_func(t, x):
        values = x.reshape(input_shape)
        values = values.transpose((1,2,0))
        # values = values * 255.
        values = (values / 255. + 0.5).clip(0, 1) * 255.

        values = values.astype('uint8')

        png = PIL.Image.fromarray(values)
        buffer = cStringIO.StringIO()
        png.save(buffer, format="PNG")
        img_str = base64.b64encode(buffer.getvalue())

        display_func._nengo_html_ = '''
            <svg width="100%%" height="100%%" viewbox="0 0 100 100">
            <image width="100%%" height="100%%"
                   xlink:href="data:image/png;base64,%s"
                   style="image-rendering: pixelated;">
            </svg>''' % (''.join(img_str))

    display_node = nengo.Node(display_func, size_in=input_image.size_out)
    nengo.Connection(input_image, display_node, synapse=None)

for obj in model.all_nodes + model.all_ensembles:
    obj.label = ''

# from nengo.utils.builder import remove_passthrough_nodes
# model = remove_passthrough_nodes(model)

# test whole network
# with model:
#     # xp = nengo.Probe(outputs['data'])
#     yp = nengo.Probe(outputs['probs'])
#     zp = nengo.Probe(outputs['logprob'])

#     spikes_p = {name: nengo.Probe(outputs[name])
#                 for name in layers if layers[name]['type'] == 'neuron'}

# nengo_gui.GUI(__file__, backend=backend).start()


text:
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 02 19:49:34 2015

@author: jeffrey.f.bryant
"""

#!/usr/bin/env python
import cv2
import numpy as np
import math

# local modules
import video

def nothing(x):
    pass

def sortCirclesBySize(circles):
    """
    Sort a list of circles by size;
    The imput is a list of [x,y,r] items where r is the size
    The return is the list sorted
    """
    sortedCircles = sorted(circles,key=lambda circles: circles[2])
    return sortedCircles

def medianColor(image,circle):
    """
    Compute te average color by creating a rectangle inside the circle,
    converting it to hsv colorspace and then computing the average
    hue, saturation and brightness
    """

    # Get a box inside the circle
    x = circle[0]
    y = circle[1]
    r = circle[2]
    span = int(r * .707)
    x1 = x-span
    x2 = x+span+1
    y1 = y-span
    y2 = y+span+1
    #img1 = image[x1:x2, y1:y2]
    img1 = image[y1:y2, x1:x2]
    if np.size(img1) > 4:
        cv2.imshow('inside',img1)

    #Compute the average color in the HSV colorspace
    hsv = cv2.cvtColor(img1, cv2.COLOR_BGR2HSV)
    hsv = np.ravel(hsv)
    try:
        hmed = np.median(hsv[0::3])
        smed = np.median(hsv[1::3])
        vmed = np.median(hsv[2::3])
    except:
        print hsv.size
        hmed = 255
        smed = 255
        vmed = 255

    # Display the average color in a window
    colorImage = np.zeros(64*64*3,np.uint8)
    colorImage.resize(64,64,3)
    colorImage[:,:,0] = hmed
    colorImage[:,:,1] = smed
    colorImage[:,:,2] = vmed
    dispImage = cv2.cvtColor(colorImage,cv2.COLOR_HSV2BGR)
    cv2.imshow('Avg Color',dispImage)

    return hmed,smed,vmed

def removeNestedCircles(circles):
    return circles

def persistantCircles(circles):
    return

def processLines(frame,cimg,edges):
#            line1 = cv2.getTrackbarPos('Line1','Detected')
#            line2 = cv2.getTrackbarPos('Line2','Detected')
#            line3 = cv2.getTrackbarPos('Line3','Detected')

    line1 = 50
    line2 = 5
    line3 = 1
            # Find and draw the first 100 lines
    lines = cv2.HoughLinesP(edges,line3, math.pi/180.0, 3, np.array([]), line1, line2)
    if lines is not None:
        a,b,c = lines.shape
        lim = b
        if lim > 100: lim = 100
        for i in range(lim):
            cv2.line(cimg, (lines[i][0][0], lines[i][0][1]), (lines[i][0][2], lines[i][0][3]), (0, 0, 255), 3, cv2.LINE_AA)

def processCircles(frame,cimg,img,cir1,cir2):
         # Find and draw the circles

        circles = cv2.HoughCircles(img,cv2.HOUGH_GRADIENT,1,20,
                                param1=cir1,param2=cir2,minRadius=70,maxRadius=600)

        if circles is not None:

           circles = np.round(circles[0, :]).astype("int")
           scircles = sortCirclesBySize(circles)

           h,s,v = medianColor(cimg,scircles[0])
           print h,s,v

           # loop over the (x, y) coordinates and radius of the circles
           for (x, y, r) in scircles:
               cv2.circle(cimg, (x, y), r, (0, 255, 0), 4)
               cv2.rectangle(cimg, (x - 5, y - 5), (x + 5, y + 5), (0, 128, 255), -1)

def main():

    showLines = False
    showCircles = True
    cam = video.create_capture(0, fallback='synth:bg=../data/baboon.jpg:class=chess:noise=0.05')

    cv2.namedWindow('Detected')

    # create trackbars for thresholds
    cv2.createTrackbar('Cir1','Detected',250,1000,nothing)
    cv2.createTrackbar('Cir2','Detected',30,100,nothing)

#    cv2.createTrackbar('Line1','Detected',50,255,nothing)
#    cv2.createTrackbar('Line2','Detected',5,100,nothing)
#    cv2.createTrackbar('Line3','Detected',1,10,nothing)

    k = 0;
    while True:
        flag, frame = cam.read()
        # cv2.imshow('camera', frame)

        gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

        img = cv2.medianBlur(gray, 5)
        cimg = frame.copy() # numpy
        #cv2.imshow('Blur',img)

        cir1 = cv2.getTrackbarPos('Cir1','Detected')
        cir2 = cv2.getTrackbarPos('Cir2','Detected')

        edges = cv2.Canny(img,cir1,cir2)

        if showLines:
            processLines(frame,cimg,edges)

        if showCircles:
            processCircles(frame,cimg,img,cir1,cir2)

        cv2.imshow('Detected',cimg)
        cv2.imshow('Edges',edges)

        # print k
        k = k + 1

        ch = 0xFF & cv2.waitKey(1)
        if ch == 27:
            break

    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()