text:
module Absorption_heat_pump_district_heating

    using JuMP
    using Cbc

    using CSV
    using DelimitedFiles
    csv_sep = ','

    function district_solver(input_folder, output_folder, tech_infos, network_id, print_function = print)
#         example_file = readdlm(string(input_folder, "\\example_file.csv"), csv_sep)
#         write_results_example(output_folder, network_id, example_file)
# 	end
        m = Model(with_optimizer(Cbc.Optimizer, logLevel=1, ratioGap = 1e+1))

        ####### Load distributions #########
        DEM=readdlm(string(input_folder, "\\DEM_time.csv"), csv_sep)

        ST_specific=readdlm(string(input_folder, "\\ST_specific_time.csv"), csv_sep)

        ST_specific_2=readdlm(string(input_folder, "\\ST_specific_time_2.csv"), csv_sep)

        ST_specific_seasonal=readdlm(string(input_folder, "\\ST_specific_time_seasonal.csv"), csv_sep)

        #28/01/2020:  from tech_infos
        #Electricity_price=readdlm(string(input_folder, "\\Electricity_price_time.csv"), csv_sep)

        Heat_exchanger_specific=readdlm(string(input_folder, "\\Heat_exchanger_specific_time.csv"), csv_sep)

        eta_HP = readdlm(string(input_folder, "\\eta_HP_1.csv"), csv_sep)
        eta_HP_2 = readdlm(string(input_folder, "\\eta_HP_2.csv"), csv_sep)
        eta_HP_3 = readdlm(string(input_folder, "\\eta_HP_3.csv"), csv_sep)

        eta_HP_waste_heat_I_1 = readdlm(string(input_folder, "\\eta_HP_I_1.csv"), csv_sep)
        eta_HP_waste_heat_I_2 = readdlm(string(input_folder, "\\eta_HP_I_1.csv"), csv_sep)
        eta_HP_waste_heat_I_3 = readdlm(string(input_folder, "\\eta_HP_I_1.csv"), csv_sep)

        eta_HP_waste_heat_II_1=readdlm(string(input_folder, "\\eta_HP_II_1.csv"), csv_sep)
        eta_HP_waste_heat_II_2=readdlm(string(input_folder, "\\eta_HP_II_2.csv"), csv_sep)
        eta_HP_waste_heat_II_3=readdlm(string(input_folder, "\\eta_HP_II_3.csv"), csv_sep)

        eta_HP_waste_heat_III_1=readdlm(string(input_folder, "\\eta_HP_III_1.csv"), csv_sep)
        eta_HP_waste_heat_III_2=readdlm(string(input_folder, "\\eta_HP_III_2.csv"), csv_sep)
        eta_HP_waste_heat_III_3=readdlm(string(input_folder, "\\eta_HP_III_3.csv"), csv_sep)

        eta_HP_waste_heat_seasonal=readdlm(string(input_folder, "\\COP_heat_pump_temperature_group_seasonal.csv"), csv_sep)
        #eta_HP_waste_heat_seasonal=eta_HP_waste_heat_seasonal[1:4380]

        Available_waste_heat_heat_pump_I_1=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_1.csv"), csv_sep)
        Available_waste_heat_heat_pump_I_2=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_2.csv"), csv_sep)
        Available_waste_heat_heat_pump_I_3=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_3.csv"), csv_sep)

        Available_waste_heat_heat_pump_II_1=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_II_1.csv"), csv_sep)
        Available_waste_heat_heat_pump_II_2=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_II_2.csv"), csv_sep)
        Available_waste_heat_heat_pump_II_3=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_II_3.csv"), csv_sep)

        Available_waste_heat_heat_pump_III_1=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_III_1.csv"), csv_sep)
        Available_waste_heat_heat_pump_III_2=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_III_2.csv"), csv_sep)
        Available_waste_heat_heat_pump_III_3=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_III_3.csv"), csv_sep)

        Available_waste_heat_heat_pump_seasonal=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_seasonal.csv"), csv_sep)
        #Available_waste_heat_heat_pump_seasonal=Available_waste_heat_heat_pump_seasonal[1:4380]

        Available_waste_heat_heat_pump_absorption=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_absorption.csv"), csv_sep)
        #Available_waste_heat_heat_pump_absorption=Available_waste_heat_heat_pump_absorption[1:4380]

        ORC_waste_heat_availability=readdlm(string(input_folder, "\\ORC_waste_heat_availability.csv"), csv_sep)

        ####### Technology input data ########

        ## Heat only boiler ##

        # Installed capacity [MW]
        P_max_HOB = tech_infos["P_max_HOB"] # 90
        # Technical minimum [MW]
        Tecnical_minimum_HOB= tech_infos["Tecnical_minimum_HOB"] # 0.0
        P_min_HOB = Tecnical_minimum_HOB*P_max_HOB
        # Efficiency [-]
        eta_HOB = tech_infos["eta_HOB"] # 0.8
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HOB= tech_infos["Percentage_ramp_up_down_HOB"] # 0.99
        ramp_up_down_HOB = P_max_HOB*Percentage_ramp_up_down_HOB
        # Variable cost [EUR/MWh]
        cost_var_HOB = tech_infos["cost_var_HOB"] # 5.4
        # Fuel cost [EUR/MWh]
        cost_fuel_HOB = tech_infos["cost_fuel_HOB"] # 25

        ## Heat only boiler 2 ##

        # Installed capacity [MW]
        P_max_HOB_2 = tech_infos["P_max_HOB_2"] # 0
        # Technical minimum [MW]
        Tecnical_minimum_HOB_2= tech_infos["Tecnical_minimum_HOB_2"] # 0.0
        P_min_HOB_2 = Tecnical_minimum_HOB_2*P_max_HOB_2
        # Efficiency [-]
        eta_HOB_2 = tech_infos["eta_HOB_2"] # 0.8
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HOB_2= tech_infos["Percentage_ramp_up_down_HOB_2"] # 0.8
        ramp_up_down_HOB_2 = P_max_HOB_2*Percentage_ramp_up_down_HOB_2
        # Variable cost [EUR/MWh]
        cost_var_HOB_2 = tech_infos["cost_var_HOB_2"] # 5.4
        # Fuel cost [EUR/MWh]
        cost_fuel_HOB_2 = tech_infos["cost_fuel_HOB_2"] # 12.5

        ## Heat only boiler 3 ##

        # Installed capacity [MW]
        P_max_HOB_3 = tech_infos["P_max_HOB_3"] # 0
        # Technical minimum [MW]
        Tecnical_minimum_HOB_3= tech_infos["Tecnical_minimum_HOB_3"] # 0.8
        P_min_HOB_3 = Tecnical_minimum_HOB_3*P_max_HOB_3
        # Efficiency [-]
        eta_HOB_3 = tech_infos["eta_HOB_3"] # 0.8
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HOB_3= tech_infos["Percentage_ramp_up_down_HOB_3"] # 0.8
        ramp_up_down_HOB_3 = P_max_HOB_3*Percentage_ramp_up_down_HOB_3
        # Variable cost [EUR/MWh]
        cost_var_HOB_3 = tech_infos["cost_var_HOB_3"] # 3.4
        # Fuel cost [EUR/MWh]
        cost_fuel_HOB_3 = tech_infos["cost_fuel_HOB_3"] # 10

        ## Electrical heater ##

        # Installed capacity [MW]
        P_max_EH = tech_infos["P_max_EH"] # 0
        # Technical minimum [MW]
        Technical_minimum_EH= tech_infos["Technical_minimum_EH"] # 0.1
        P_min_EH = P_max_EH*Technical_minimum_EH
        # Efficiency [-]
        eta_EH = tech_infos["eta_EH"] # 0.98
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_EH= tech_infos["Percentage_ramp_up_down_EH"] # 0.95
        ramp_up_down_EH = Percentage_ramp_up_down_EH*P_max_EH
        # Variable cost [EUR/MWh]
        cost_var_EH = tech_infos["cost_var_EH"] # 0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_EH = tech_infos["cost_fuel_EH"] # Electricity_price

        ## Electrical heater 2 ##

        # Installed capacity [MW]
        P_max_EH_2 = tech_infos["P_max_EH_2"] # 0
        # Technical minimum [MW]
        Technical_minimum_EH_2= tech_infos["Technical_minimum_EH_2"] # 0.15
        P_min_EH_2 = P_max_EH_2*Technical_minimum_EH_2
        # Efficiency [-]
        eta_EH_2 = tech_infos["eta_EH_2"] # 0.98
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_EH_2=tech_infos["Percentage_ramp_up_down_EH_2"] # 0.9
        ramp_up_down_EH_2 = Percentage_ramp_up_down_EH_2*P_max_EH_2
        # Variable cost [EUR/MWh]
        cost_var_EH_2 = tech_infos["cost_var_EH_2"] # 0.7
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_EH_2 = tech_infos["cost_fuel_EH_2"] # Electricity_price

        ## Electrical heater 3 ##

        # Installed capacity [MW]
        P_max_EH_3 = tech_infos["P_max_EH_3"] # 0
        # Technical minimum [MW]
        Technical_minimum_EH_3= tech_infos["Technical_minimum_EH_3"] # 0.1
        P_min_EH_3 = P_max_EH_3*Technical_minimum_EH_3
        # Efficiency [-]
        eta_EH_3 = tech_infos["eta_EH_3"] # 0.98
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_EH_3= tech_infos["Percentage_ramp_up_down_EH_3"] # 0.99
        ramp_up_down_EH_3 = Percentage_ramp_up_down_EH_3*P_max_EH_3
        # Variable cost [EUR/MWh]
        cost_var_EH_3 = tech_infos["cost_var_EH_3"] # 0.8
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_EH_3 = tech_infos["cost_fuel_EH_3"] # Electricity_price

        ## Heat pump ##
        # Installed capacity [MW]
        P_max_HP = tech_infos["P_max_HP"] # 0
        # Technical minimum [MW]
        Technical_minimum_HP= tech_infos["Technical_minimum_HP"] # 0.1
        P_min_HP = Technical_minimum_HP*P_max_HP
        # Efficiency [-]
        eta_HP = eta_HP
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP= tech_infos["Percentage_ramp_up_down_HP"] # 0.95
        ramp_up_down_HP = Percentage_ramp_up_down_HP*P_max_HP
        # Variable cost [EUR/MWh]
        cost_var_HP = tech_infos["cost_var_HP"] # 0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP = tech_infos["cost_fuel_HP"] # Electricity_price

        ## Heat pump 2 ##

        # Installed capacity [MW]
        P_max_HP_2 = tech_infos["P_max_HP_2"] # 0
        # Technical minimum [MW]
        Technical_minimum_HP_2= tech_infos["Technical_minimum_HP_2"] # 0.1
        P_min_HP_2 = Technical_minimum_HP_2*P_max_HP_2
        # Efficiency [-]
        eta_HP_2 = eta_HP_2
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_2= tech_infos["Percentage_ramp_up_down_HP_2"] # 0.95
        ramp_up_down_HP_2 = Percentage_ramp_up_down_HP_2*P_max_HP_2
        # Variable cost [EUR/MWh]
        cost_var_HP_2 = tech_infos["cost_var_HP_2"] # 0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_2 = tech_infos["cost_fuel_HP_2"] # Electricity_price

        ## Heat pump 3 ##
        # Installed capacity [MW]
        P_max_HP_3 = tech_infos["P_max_HP_3"] # 0
        # Technical minimum [MW]
        Technical_minimum_HP_3= tech_infos["Technical_minimum_HP_3"] # 0.1
        P_min_HP_3 = Technical_minimum_HP_3*P_max_HP_3
        # Efficiency [-]
        eta_HP_3 = eta_HP_3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_3= tech_infos["Percentage_ramp_up_down_HP_3"] # 0.95
        ramp_up_down_HP_3 = Percentage_ramp_up_down_HP_3*P_max_HP_3
        # Variable cost [EUR/MWh]
        cost_var_HP_3 = tech_infos["cost_var_HP_3"] # 0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_3 = tech_infos["cost_fuel_HP_3"] # Electricity_price

        ## Absorption heat pump ##

        # Installed capacity [MW]
        P_max_HP_absorption = tech_infos["P_max_HP_absorption"] # 30
        # Technical minimum [MW]
        Technical_minimum_HP_absorption= tech_infos["Technical_minimum_HP_absorption"] # 0.0
        P_min_HP_absorption = Technical_minimum_HP_absorption*P_max_HP_absorption
        # Efficiency [-]
        eta_HP_absorption = tech_infos["eta_HP_absorption"] # 1.4
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_absorption=tech_infos["Percentage_ramp_up_down_HP_absorption"] # 0.95
        ramp_up_down_HP_absorption = Percentage_ramp_up_down_HP_absorption*P_max_HP_absorption
        # Variable cost [EUR/MWh]
        cost_var_HP_absorption = tech_infos["cost_var_HP_absorption"] # 0.5
        # Fuel cost [EUR/MWh]
        cost_fuel_HP_absorption = tech_infos["cost_fuel_HP_absorption"] # 25

        ## Cogeneration ##

        # Installed capacity [MW]
        P_max_CHP = tech_infos["P_max_CHP"] # 0
        # Technical minimum [MW]
        Technical_minimum_CHP= tech_infos["Technical_minimum_CHP"] # 0.4
        P_min_CHP = Technical_minimum_CHP*P_max_CHP
        # Efficiency [-]
        eta_CHP_th = tech_infos["eta_CHP_th"] # 0.65
        # Power-to-heat ratio [-]
        cb_CHP = tech_infos["cb_CHP"] # 0.3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_CHP= tech_infos["Percentage_ramp_up_down_CHP"] # 0.2
        ramp_up_down_CHP = Percentage_ramp_up_down_CHP*P_max_CHP
        # Variable cost [EUR/MWh]
        cost_var_CHP = tech_infos["cost_var_CHP"] # 3.9
        # Fuel cost [EUR/MWh]
        cost_fuel_CHP = tech_infos["cost_fuel_CHP"] # 20
        # Price paid to CHP for electrcity production
        CHP_electricity_price = tech_infos["CHP_electricity_price"] # Electricity_price

        ## Cogeneration 2 ##

        # Installed capacity [MW]
        P_max_CHP_2 = tech_infos["P_max_CHP_2"] # 0
        # Technical minimum [MW]
        Technical_minimum_CHP_2= tech_infos["Technical_minimum_CHP_2"] # 0.4
        P_min_CHP_2 = Technical_minimum_CHP_2*P_max_CHP_2
        # Efficiency [-]
        eta_CHP_th_2 = tech_infos["eta_CHP_th_2"] # 0.65
        # Power-to-heat ratio [-]
        cb_CHP_2 = tech_infos["cb_CHP_2"] # 0.3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_CHP_2= tech_infos["Percentage_ramp_up_down_CHP_2"] # 0.2
        ramp_up_down_CHP_2 = Percentage_ramp_up_down_CHP_2*P_max_CHP_2
        # Variable cost [EUR/MWh]
        cost_var_CHP_2 = tech_infos["cost_var_CHP_2"] # 3.9
        # Fuel cost [EUR/MWh]
        cost_fuel_CHP_2 = tech_infos["cost_fuel_CHP_2"] # 15
        # Price paid to CHP for electrcity production
        CHP_electricity_price_2 = tech_infos["CHP_electricity_price_2"] # Electricity_price

        ## ORC Cogeneration ##

        # THERMAL installed capacity [MW]
        P_max_CHP_ORC =tech_infos["P_max_CHP_ORC"] #0
        # Technical minimum [MW]
        Technical_minimum_CHP_ORC=tech_infos["Technical_minimum_CHP_ORC"] #0
        P_min_CHP_ORC = Technical_minimum_CHP_ORC*P_max_CHP_ORC
        # Electrical efficiency [-]
        eta_CHP_el = tech_infos["eta_CHP_el"]
        # Power-to-heat ratio [-]
        cb_CHP_ORC =tech_infos["cb_CHP_ORC"]
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_CHP_ORC=tech_infos["Percentage_ramp_up_down_CHP_ORC"]
        ramp_up_down_CHP_ORC = Percentage_ramp_up_down_CHP_ORC*P_max_CHP_ORC
        # Variable cost [EUR/MWh]
        cost_var_CHP_ORC = tech_infos["cost_var_CHP_ORC"]
        # Price paid to CHP for electrcity production
        CHP_ORC_electricity_price = tech_infos["CHP_ORC_electricity_price"] # Electricity_price
        # ORC Thermal availability
        ORC_thermal_availability=(ORC_waste_heat_availability*eta_CHP_el)/cb_CHP_ORC

        ## Solar thermal ##

        # Installed solar thermal collector area [m2]
        A_ST= tech_infos["A_ST"] # 0
        # Variable cost [EUR/MWh]
        cost_var_ST = tech_infos["cost_var_ST"] # 0.5
        # Solar thermal collector prodution [MW]
        P_max_ST=A_ST*ST_specific/1000000

        ## Solar thermal 2 ##

        # Installed solar thermal collector area [m2]
        A_ST_2= tech_infos["A_ST_2"] # 0
        # Variable cost [EUR/MWh]
        cost_var_ST_2 = tech_infos["cost_var_ST_2"] # 0.5
        # Solar thermal collector prodution [MW]
        P_max_ST_2=A_ST_2*ST_specific_2/1000000

        ## Seasonal solar thermal ##
        # Installed solar thermal collector area [m2]
        A_ST_seasonal= tech_infos["A_ST_seasonal"] # 0
        # Variable cost [EUR/MWh]
        cost_var_ST_seasonal = tech_infos["cost_var_ST_seasonal"] # 0.5
        # Solar thermal collector prodution [MW]
        P_max_ST_seasonal=A_ST_seasonal*ST_specific_seasonal/1000000

        ## Waste heat - heat exchangers ##

        # Installed capacity [MW]
        P_HEX_capacity = tech_infos["P_HEX_capacity"] # 0
        P_max_HEX=P_HEX_capacity*Heat_exchanger_specific
        # Technical minimum [MW]
        P_min_HEX= tech_infos["P_min_HEX"] # 0
        # Variable cost [EUR/MWh]
        cost_var_HEX = tech_infos["cost_var_HEX"] # 0.25

        ## Waste heat - heat pumps 1 ##

        # Installed capacity [MW]
        P_max_HP_waste_heat_I_1 = tech_infos["P_max_HP_waste_heat_I_1"] # 0
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_I_1 = tech_infos["Technical_minimum_HP_waste_heat_I_1"] # 0.0
        P_min_HP_waste_heat_I_1 = Technical_minimum_HP_waste_heat_I_1*P_max_HP_waste_heat_I_1
        # Efficiency [-]
        eta_HP_waste_heat_I_1 = eta_HP_waste_heat_I_1
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_I_1= tech_infos["Percentage_ramp_up_down_HP_waste_heat_I_1"] # 0.95
        ramp_up_down_HP_waste_heat_I_1 = Percentage_ramp_up_down_HP_waste_heat_I_1*P_max_HP_waste_heat_I_1
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_I_1 = tech_infos["cost_var_HP_waste_heat_I_1"] # 0.7
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_I_1 = tech_infos["cost_fuel_HP_waste_heat_I_1"] # Electricity_price

        ## Waste heat - heat pumps 2 ##

        # Installed capacity [MW]
        P_max_HP_waste_heat_I_2 = tech_infos["P_max_HP_waste_heat_I_2"] # 0
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_I_2 = tech_infos["Technical_minimum_HP_waste_heat_I_2"] # 0.0
        P_min_HP_waste_heat_I_2 = Technical_minimum_HP_waste_heat_I_2*P_max_HP_waste_heat_I_2
        # Efficiency [-]
        eta_HP_waste_heat_2 = eta_HP_waste_heat_I_2
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_I_2= tech_infos["Percentage_ramp_up_down_HP_waste_heat_I_2"] # 0.95
        ramp_up_down_HP_waste_heat_I_2 = Percentage_ramp_up_down_HP_waste_heat_I_2*P_max_HP_waste_heat_I_2
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_I_2 = tech_infos["cost_var_HP_waste_heat_I_2"] # 1.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_I_2 = tech_infos["cost_fuel_HP_waste_heat_I_2"] # Electricity_price

        ## Waste heat - heat pumps 3##

        # Installed capacity [MW]
        P_max_HP_waste_heat_I_3 = tech_infos["P_max_HP_waste_heat_I_3"] # 0
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_I_3 = tech_infos["Technical_minimum_HP_waste_heat_I_3"] # 0.0
        P_min_HP_waste_heat_I_3 = Technical_minimum_HP_waste_heat_I_3*P_max_HP_waste_heat_I_3
        # Efficiency [-]
        eta_HP_waste_heat_I_3 = eta_HP_waste_heat_I_3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_I_3= tech_infos["Percentage_ramp_up_down_HP_waste_heat_I_3"] # 0.95
        ramp_up_down_HP_waste_heat_I_3 = Percentage_ramp_up_down_HP_waste_heat_I_3*P_max_HP_waste_heat_I_3
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_I_3 = tech_infos["cost_var_HP_waste_heat_I_3"] # 1.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_I_3 = tech_infos["cost_fuel_HP_waste_heat_I_3"] # Electricity_price

        ## Waste heat - heat pumps - temperature group II ##

        # Heat pump II_1 #

        # Installed capacity [MW]
        P_max_HP_waste_heat_II_1 = tech_infos["P_max_HP_waste_heat_II_1"]
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_II_1 =tech_infos["Technical_minimum_HP_waste_heat_II_1"]
        P_min_HP_waste_heat_II_1 = Technical_minimum_HP_waste_heat_II_1*P_max_HP_waste_heat_II_1
        # Efficiency [-]
        eta_HP_waste_heat_II_1 = eta_HP_waste_heat_II_1
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_II_1=tech_infos["Percentage_ramp_up_down_HP_waste_heat_II_1"]
        ramp_up_down_HP_waste_heat_II_1 = Percentage_ramp_up_down_HP_waste_heat_II_1*P_max_HP_waste_heat_II_1
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_II_1 =tech_infos["cost_var_HP_waste_heat_II_1"]
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_II_1 = tech_infos["cost_fuel_HP_waste_heat_II_1"] # Electricity_price

        # Heat pump II_2 #

        # Installed capacity [MW]
        P_max_HP_waste_heat_II_2 =tech_infos["P_max_HP_waste_heat_II_2"]
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_II_2 =tech_infos["Technical_minimum_HP_waste_heat_II_2"]
        P_min_HP_waste_heat_II_2 = Technical_minimum_HP_waste_heat_II_2*P_max_HP_waste_heat_II_2
        # Efficiency [-]
        eta_HP_waste_heat_II_2 = eta_HP_waste_heat_II_2
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_II_2=tech_infos["Percentage_ramp_up_down_HP_waste_heat_II_2"]
        ramp_up_down_HP_waste_heat_II_2 = Percentage_ramp_up_down_HP_waste_heat_II_2*P_max_HP_waste_heat_II_2
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_II_2 =tech_infos["cost_var_HP_waste_heat_II_2"]
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_II_2 = tech_infos["cost_fuel_HP_waste_heat_II_2"] # Electricity_price

        # Heat pump II_3 #

        # Installed capacity [MW]
        P_max_HP_waste_heat_II_3 = tech_infos["P_max_HP_waste_heat_II_3"]
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_II_3 =tech_infos["Technical_minimum_HP_waste_heat_II_3"]
        P_min_HP_waste_heat_II_3 = Technical_minimum_HP_waste_heat_II_3*P_max_HP_waste_heat_II_3
        # Efficiency [-]
        eta_HP_waste_heat_II_3 = eta_HP_waste_heat_II_3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_II_3=tech_infos["Percentage_ramp_up_down_HP_waste_heat_II_3"]
        ramp_up_down_HP_waste_heat_II_3 = Percentage_ramp_up_down_HP_waste_heat_II_3*P_max_HP_waste_heat_II_3
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_II_3 = tech_infos["cost_var_HP_waste_heat_II_3"]
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_II_3 = tech_infos["cost_fuel_HP_waste_heat_II_3"] # Electricity_price

        ## Waste heat - heat pumps - temperature group III ##

        # Heat pump III_1 #

        # Installed capacity [MW]
        P_max_HP_waste_heat_III_1 = tech_infos["P_max_HP_waste_heat_III_1"]
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_III_1 =tech_infos["Technical_minimum_HP_waste_heat_III_1"]
        P_min_HP_waste_heat_III_1 = Technical_minimum_HP_waste_heat_III_1*P_max_HP_waste_heat_III_1
        # Efficiency [-]
        eta_HP_waste_heat_III_1 = eta_HP_waste_heat_III_1
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_III_1=tech_infos["Percentage_ramp_up_down_HP_waste_heat_III_1"]
        ramp_up_down_HP_waste_heat_III_1 = Percentage_ramp_up_down_HP_waste_heat_III_1*P_max_HP_waste_heat_III_1
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_III_1 = tech_infos["cost_var_HP_waste_heat_III_1"]
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_III_1 = tech_infos["cost_fuel_HP_waste_heat_III_1"] # Electricity_price

        # Heat pump III_2 #

        # Installed capacity [MW]
        P_max_HP_waste_heat_III_2 = tech_infos["P_max_HP_waste_heat_III_2"]
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_III_2 =tech_infos["Technical_minimum_HP_waste_heat_III_2"]
        P_min_HP_waste_heat_III_2 = Technical_minimum_HP_waste_heat_III_2*P_max_HP_waste_heat_III_2
        # Efficiency [-]
        eta_HP_waste_heat_III_2 = eta_HP_waste_heat_III_2
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_III_2=tech_infos["Percentage_ramp_up_down_HP_waste_heat_III_2"]
        ramp_up_down_HP_waste_heat_III_2 = Percentage_ramp_up_down_HP_waste_heat_III_2*P_max_HP_waste_heat_III_2
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_III_2 = tech_infos["cost_var_HP_waste_heat_III_2"]
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_III_2 = tech_infos["cost_fuel_HP_waste_heat_III_2"] # Electricity_price

        # Heat pump III_3 #

        # Installed capacity [MW]
        P_max_HP_waste_heat_III_3 = tech_infos["P_max_HP_waste_heat_III_3"]
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_III_3 =tech_infos["Technical_minimum_HP_waste_heat_III_3"]
        P_min_HP_waste_heat_III_3 = Technical_minimum_HP_waste_heat_III_3*P_max_HP_waste_heat_III_3
        # Efficiency [-]
        eta_HP_waste_heat_III_3 = eta_HP_waste_heat_III_3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_III_3=tech_infos["Percentage_ramp_up_down_HP_waste_heat_III_3"]
        ramp_up_down_HP_waste_heat_III_3 = Percentage_ramp_up_down_HP_waste_heat_III_3*P_max_HP_waste_heat_III_3
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_III_3 = tech_infos["cost_var_HP_waste_heat_III_3"]
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_III_3 = tech_infos["cost_fuel_HP_waste_heat_III_3"] # Electricity_price

        ## Seasonal waste heat heat pumps ##

        # Installed capacity [MW]
        P_max_HP_waste_heat_seasonal = tech_infos["P_max_HP_waste_heat_seasonal"] # 0
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_seasonal = tech_infos["Technical_minimum_HP_waste_heat_seasonal"] # 0.0
        P_min_HP_waste_heat_seasonal = Technical_minimum_HP_waste_heat_seasonal*P_max_HP_waste_heat_seasonal
        # Efficiency [-]
        eta_HP_waste_heat_seasonal = eta_HP_waste_heat_seasonal
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_seasonal= tech_infos["Percentage_ramp_up_down_HP_waste_heat_seasonal"] # 0.95
        ramp_up_down_HP_waste_heat_seasonal = Percentage_ramp_up_down_HP_waste_heat_seasonal*P_max_HP_waste_heat_seasonal
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_seasonal = tech_infos["cost_var_HP_waste_heat_seasonal"] # 0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_seasonal = tech_infos["cost_fuel_HP_waste_heat_seasonal"] # Electricity_price

        ## Waste heat absorption heat pump ##

        # Installed capacity [MW]
        P_max_HP_waste_heat_absorption = tech_infos["P_max_HP_waste_heat_absorption"] # 30
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_absorption = tech_infos["Technical_minimum_HP_waste_heat_absorption"] # 0.0
        P_min_HP_waste_heat_absorption = Technical_minimum_HP_waste_heat_absorption*P_max_HP_waste_heat_absorption
        # Efficiency [-]
        eta_HP_waste_heat_absorption = tech_infos["eta_HP_waste_heat_absorption"] # 1.5
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_absorption= tech_infos["Percentage_ramp_up_down_HP_waste_heat_absorption"] # 0.99
        ramp_up_down_HP_waste_heat_absorption = Percentage_ramp_up_down_HP_waste_heat_absorption*P_max_HP_waste_heat_absorption
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_absorption = tech_infos["cost_var_HP_waste_heat_absorption"] # 0.5
        # Fuel cost [EUR/MWh]
        cost_fuel_HP_waste_heat_absorption = tech_infos["cost_fuel_HP_waste_heat_absorption"] # 25

        ## Thermal energy storage ##
        # Thermal energy storage size [MWh]
        TES_size = tech_infos["TES_size"] # 50
        # Minimum percentage of state of charge of thermal energy storage [-]
        SOC_min = tech_infos["SOC_min"] # 0.15
        # Thermal energy storage state of charge at the start/end of the year [-]
        TES_start_end = tech_infos["TES_start_end"] # 0.5
        # Thermal storage charge/discharge capacity [MW]
        TES_charge_discharge_time = tech_infos["TES_charge_discharge_time"] # 5
        TES_in_out_max = TES_size/TES_charge_discharge_time
        TES_loss= tech_infos["TES_loss"] # 0.1/(24*7)

        ## Seasonal thermal energy storage ##

        # Thermal energy storage size [MWh]
        TES_size_seasonal = tech_infos["TES_size_seasonal"] # 0
        # Minimum percentage of state of charge of thermal energy storage [-]
        SOC_min_seasonal = tech_infos["SOC_min_seasonal"] # 0.0
        # Thermal energy storage state of charge at the start/end of the year [-]
        TES_start_end_seasonal = tech_infos["TES_start_end_seasonal"] # 0.0
        # Thermal storage charge/discharge capacity [MW]
        TES_charge_discharge_time_seasonal = tech_infos["TES_charge_discharge_time_seasonal"] # 2000
        TES_in_out_max_seasonal = TES_size_seasonal/TES_charge_discharge_time_seasonal
        TES_loss_seasonal= tech_infos["TES_loss_seasonal"] # 0.01/(24*7)

        ######## Constraints and variables ########

        ### Variables ###

        # Technologies' operation
        @variable(m, P_max_HOB >= HOB[1:length(DEM)] >= 0)
        @variable(m, P_max_HOB_2 >= HOB_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HOB_3 >= HOB_3[1:length(DEM)] >= 0)

        @variable(m, P_max_EH >= EH[1:length(DEM)] >= 0)
        @variable(m, P_max_EH_2 >= EH_2[1:length(DEM)] >= 0)
        @variable(m, P_max_EH_3 >= EH_3[1:length(DEM)] >= 0)

        @variable(m, P_max_HP >= HP[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_2 >= HP_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_3 >= HP_3[1:length(DEM)] >= 0)

        @variable(m, P_max_HP_absorption >= HP_absorption[1:length(DEM)] >= 0)

        @variable(m, P_max_CHP >= CHP[1:length(DEM)] >= 0)
        @variable(m, P_max_CHP_2 >= CHP_2[1:length(DEM)] >= 0)

        @variable(m, P_max_CHP_ORC >= CHP_ORC[1:length(DEM)] >= 0)

        @variable(m, ST[1:length(DEM)] >= 0)
        @variable(m, ST_2[1:length(DEM)] >= 0)
        @variable(m, ST_seasonal[1:length(DEM)] >= 0)

        @variable(m, HEX[1:length(DEM)] >= 0)

        @variable(m, P_max_HP_waste_heat_I_1 >= HP_waste_heat_I_1[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_waste_heat_I_2 >= HP_waste_heat_I_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_waste_heat_I_3 >= HP_waste_heat_I_3[1:length(DEM)] >= 0)

        @variable(m, P_max_HP_waste_heat_II_1 >= HP_waste_heat_II_1[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_waste_heat_II_2 >= HP_waste_heat_II_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_waste_heat_II_3 >= HP_waste_heat_II_3[1:length(DEM)] >= 0)

        @variable(m, P_max_HP_waste_heat_III_1 >= HP_waste_heat_III_1[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_waste_heat_III_2 >= HP_waste_heat_III_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_waste_heat_III_3 >= HP_waste_heat_III_3[1:length(DEM)] >= 0)

        @variable(m, P_max_HP_waste_heat_seasonal >= HP_waste_heat_seasonal[1:length(DEM)] >= 0)

        @variable(m, P_max_HP_waste_heat_absorption >= HP_waste_heat_absorption[1:length(DEM)] >= 0)
        @variable(m, TES_in_out_max >=TES[1:length(DEM)]>= -TES_in_out_max)
        @variable(m, TES_size >= SOC[1:length(DEM)] >= TES_size*SOC_min)

        @variable(m, 0 >=TES_seasonal[1:length(DEM)]>= -TES_in_out_max_seasonal)
        @variable(m, TES_size_seasonal >= SOC_seasonal[1:length(DEM)] >= TES_size*SOC_min_seasonal)

        #@variable(m, maximum(DEM) >= Error[1:length(DEM)] >= 0)
        @variable(m, 0 >= Error[1:length(DEM)] >= 0)

        # Binary variables, on/off

        @variable(m, OnOffHOB[1:length(DEM)], Bin)
        @variable(m, OnOffHOB_2[1:length(DEM)], Bin)
        @variable(m, OnOffHOB_3[1:length(DEM)], Bin)

        @variable(m, OnOffEH[1:length(DEM)], Bin)
        @variable(m, OnOffEH_2[1:length(DEM)], Bin)
        @variable(m, OnOffEH_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP[1:length(DEM)], Bin)
        @variable(m, OnOffHP_2[1:length(DEM)], Bin)
        @variable(m, OnOffHP_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP_absorption[1:length(DEM)], Bin)

        @variable(m, OnOffCHP[1:length(DEM)], Bin)
        @variable(m, OnOffCHP_2[1:length(DEM)], Bin)

        @variable(m, OnOffCHP_ORC[1:length(DEM)], Bin)

        @variable(m, OnOffHP_waste_heat_I_1[1:length(DEM)], Bin)
        @variable(m, OnOffHP_waste_heat_I_2[1:length(DEM)], Bin)
        @variable(m, OnOffHP_waste_heat_I_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP_waste_heat_II_1[1:length(DEM)], Bin)
        @variable(m, OnOffHP_waste_heat_II_2[1:length(DEM)], Bin)
        @variable(m, OnOffHP_waste_heat_II_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP_waste_heat_III_1[1:length(DEM)], Bin)
        @variable(m, OnOffHP_waste_heat_III_2[1:length(DEM)], Bin)
        @variable(m, OnOffHP_waste_heat_III_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP_waste_heat_seasonal[1:length(DEM)], Bin)

        @variable(m, OnOffHP_waste_heat_absorption[1:length(DEM)], Bin)

        ### Constraints ###
        # Supply + TES = DEM
        @constraint(m, HOB+HOB_2+HOB_3+EH+EH_2+EH_3+HP+HP_2+HP_3+HP_absorption+CHP+CHP_2+CHP_ORC+ST+ST_2+HEX+HP_waste_heat_I_1+HP_waste_heat_I_2+HP_waste_heat_I_3+HP_waste_heat_II_1+HP_waste_heat_II_2+HP_waste_heat_II_3+HP_waste_heat_III_1+HP_waste_heat_III_2+HP_waste_heat_III_3+HP_waste_heat_absorption-TES-TES_seasonal + Error .== DEM)

        # Thermal storage constraints
        t=2:(length(DEM)-1)
        @constraint(m, SOC[t]-SOC[t .- 1]-TES[t]+SOC[t]*TES_loss .== 0)
        @constraint(m, SOC[length(DEM)-1]+TES[length(DEM)] == TES_size*TES_start_end)
        @constraint(m, SOC[1]-TES[1] == TES_size*TES_start_end)
        @constraint(m, SOC[length(DEM)] == TES_size*TES_start_end)

        # Seasonal thermal storage constraints
        @constraint(m, SOC_seasonal[t]-SOC_seasonal[t .- 1]-TES_seasonal[t]+SOC_seasonal[t]*TES_loss_seasonal - HP_waste_heat_seasonal[t]-ST_seasonal[t].== 0)
        @constraint(m, SOC_seasonal[length(DEM)-1]+TES_seasonal[length(DEM)]+HP_waste_heat_seasonal[length(DEM)]+ST_seasonal[length(DEM)] == TES_size_seasonal*TES_start_end_seasonal)
        @constraint(m, SOC_seasonal[1]-TES_seasonal[1]+HP_waste_heat_seasonal[1]+ST_seasonal[1] == TES_size_seasonal*TES_start_end_seasonal)
        @constraint(m, SOC_seasonal[length(DEM)] == TES_size_seasonal*TES_start_end_seasonal)

        # Solar thermal collector constraints
        #@constraint(m, ST .== P_max_ST)
        @constraint(m, ST .<= P_max_ST)
        #@constraint(m, ST_2 .== P_max_ST_2)
        @constraint(m, ST_2 .<= P_max_ST_2)

        #@constraint(m, ST_seasonal .== P_max_ST_seasonal)
        @constraint(m, ST_seasonal .<= P_max_ST_seasonal)

        # Waste heat constraints - can't utilize more than available
        @constraint(m, HEX .<=P_max_HEX)

        @constraint(m, HP_waste_heat_I_1 .<= Available_waste_heat_heat_pump_I_1)
        @constraint(m, HP_waste_heat_I_2 .<= Available_waste_heat_heat_pump_I_2)
        @constraint(m, HP_waste_heat_I_3 .<= Available_waste_heat_heat_pump_I_3)

        @constraint(m, HP_waste_heat_II_1 .<= Available_waste_heat_heat_pump_II_1)
        @constraint(m, HP_waste_heat_II_2 .<= Available_waste_heat_heat_pump_II_2)
        @constraint(m, HP_waste_heat_II_3 .<= Available_waste_heat_heat_pump_II_3)

        @constraint(m, HP_waste_heat_III_1 .<= Available_waste_heat_heat_pump_III_1)
        @constraint(m, HP_waste_heat_III_2 .<= Available_waste_heat_heat_pump_III_2)
        @constraint(m, HP_waste_heat_III_3 .<= Available_waste_heat_heat_pump_III_3)

        @constraint(m, HP_waste_heat_seasonal .<= Available_waste_heat_heat_pump_seasonal)

        @constraint(m, HP_waste_heat_absorption .<= Available_waste_heat_heat_pump_absorption)

        @constraint(m, CHP_ORC.<= ORC_thermal_availability)

        # OnOff constraints - technical minimum
        @constraint(m, OnOffHOB*P_min_HOB-HOB.<=0)
        @constraint(m, -OnOffHOB*P_max_HOB+HOB.<=0)
        @constraint(m, OnOffHOB_2*P_min_HOB_2-HOB_2.<=0)
        @constraint(m, -OnOffHOB_2*P_max_HOB_2+HOB_2.<=0)
        @constraint(m, OnOffHOB_3*P_min_HOB_3-HOB_3.<=0)
        @constraint(m, -OnOffHOB_3*P_max_HOB_3+HOB_3.<=0)

        @constraint(m, OnOffEH*P_min_EH-EH.<=0)
        @constraint(m, -OnOffEH*P_max_EH+EH.<=0)
        @constraint(m, OnOffEH_2*P_min_EH_2-EH_2.<=0)
        @constraint(m, -OnOffEH_2*P_max_EH_2+EH_2.<=0)
        @constraint(m, OnOffEH_3*P_min_EH_3-EH_3.<=0)
        @constraint(m, -OnOffEH_3*P_max_EH_3+EH_3.<=0)

        @constraint(m, OnOffHP*P_min_HP-HP.<=0)
        @constraint(m, -OnOffHP*P_max_HP+HP.<=0)
        @constraint(m, OnOffHP_2*P_min_HP_2-HP_2.<=0)
        @constraint(m, -OnOffHP_2*P_max_HP_2+HP_2.<=0)
        @constraint(m, OnOffHP_3*P_min_HP_3-HP_3.<=0)
        @constraint(m, -OnOffHP_3*P_max_HP_3+HP_3.<=0)

        @constraint(m, OnOffHP_absorption*P_min_HP_absorption-HP_absorption.<=0)
        @constraint(m, -OnOffHP_absorption*P_max_HP_absorption+HP_absorption.<=0)

        @constraint(m, OnOffCHP*P_min_CHP-CHP.<=0)
        @constraint(m, -OnOffCHP*P_max_CHP+CHP.<=0)
        @constraint(m, OnOffCHP_2*P_min_CHP_2-CHP_2.<=0)
        @constraint(m, -OnOffCHP_2*P_max_CHP_2+CHP_2.<=0)

        @constraint(m, OnOffCHP_ORC*P_min_CHP_ORC-CHP_ORC.<=0)
        @constraint(m, -OnOffCHP_ORC*P_max_CHP_ORC+CHP_ORC.<=0)

        @constraint(m, OnOffHP_waste_heat_I_1*P_min_HP_waste_heat_I_1-HP_waste_heat_I_1.<=0)
        @constraint(m, -OnOffHP_waste_heat_I_1*P_max_HP_waste_heat_I_1+HP_waste_heat_I_1.<=0)
        @constraint(m, OnOffHP_waste_heat_I_2*P_min_HP_waste_heat_I_2-HP_waste_heat_I_2.<=0)
        @constraint(m, -OnOffHP_waste_heat_I_2*P_max_HP_waste_heat_I_2+HP_waste_heat_I_2.<=0)
        @constraint(m, OnOffHP_waste_heat_I_3*P_min_HP_waste_heat_I_3-HP_waste_heat_I_3.<=0)
        @constraint(m, -OnOffHP_waste_heat_I_3*P_max_HP_waste_heat_I_3+HP_waste_heat_I_3.<=0)

        @constraint(m, OnOffHP_waste_heat_II_1*P_min_HP_waste_heat_II_1-HP_waste_heat_II_1.<=0)
        @constraint(m, -OnOffHP_waste_heat_II_1*P_max_HP_waste_heat_II_1+HP_waste_heat_II_1.<=0)
        @constraint(m, OnOffHP_waste_heat_II_2*P_min_HP_waste_heat_II_2-HP_waste_heat_II_2.<=0)
        @constraint(m, -OnOffHP_waste_heat_II_2*P_max_HP_waste_heat_II_2+HP_waste_heat_II_2.<=0)
        @constraint(m, OnOffHP_waste_heat_II_3*P_min_HP_waste_heat_II_3-HP_waste_heat_II_3.<=0)
        @constraint(m, -OnOffHP_waste_heat_II_3*P_max_HP_waste_heat_II_3+HP_waste_heat_II_3.<=0)

        @constraint(m, OnOffHP_waste_heat_III_1*P_min_HP_waste_heat_III_1-HP_waste_heat_III_1.<=0)
        @constraint(m, -OnOffHP_waste_heat_III_1*P_max_HP_waste_heat_III_1+HP_waste_heat_III_1.<=0)
        @constraint(m, OnOffHP_waste_heat_III_2*P_min_HP_waste_heat_III_2-HP_waste_heat_III_2.<=0)
        @constraint(m, -OnOffHP_waste_heat_III_2*P_max_HP_waste_heat_III_2+HP_waste_heat_III_2.<=0)
        @constraint(m, OnOffHP_waste_heat_III_3*P_min_HP_waste_heat_III_3-HP_waste_heat_III_3.<=0)
        @constraint(m, -OnOffHP_waste_heat_III_3*P_max_HP_waste_heat_III_3+HP_waste_heat_III_3.<=0)

        @constraint(m, OnOffHP_waste_heat_seasonal*P_min_HP_waste_heat_seasonal-HP_waste_heat_seasonal.<=0)
        @constraint(m, -OnOffHP_waste_heat_seasonal*P_max_HP_waste_heat_seasonal+HP_waste_heat_seasonal.<=0)

        @constraint(m, OnOffHP_waste_heat_absorption*P_min_HP_waste_heat_absorption-HP_waste_heat_absorption.<=0)
        @constraint(m, -OnOffHP_waste_heat_absorption*P_max_HP_waste_heat_absorption+HP_waste_heat_absorption.<=0)

        # Ramp-up-down constraint for each technology
        i=2:(length(DEM)-1)

        @constraint(m, HOB[i .- 1]-HOB[i].<=ramp_up_down_HOB*OnOffHOB[i]+(1 .- OnOffHOB[i])*P_min_HOB)
        @constraint(m, HOB[i]-HOB[i .- 1].<=ramp_up_down_HOB*OnOffHOB[i .- 1]+(1 .- OnOffHOB[i .- 1])*P_min_HOB)
        @constraint(m, HOB_2[i .- 1]-HOB_2[i].<=ramp_up_down_HOB_2*OnOffHOB_2[i]+(1 .- OnOffHOB_2[i])*P_min_HOB_2)
        @constraint(m, HOB_2[i]-HOB_2[i .- 1].<=ramp_up_down_HOB_2*OnOffHOB_2[i .- 1]+(1 .- OnOffHOB_2[i .- 1])*P_min_HOB_2)
        @constraint(m, HOB_3[i .- 1]-HOB_3[i].<=ramp_up_down_HOB_3*OnOffHOB_3[i]+(1 .- OnOffHOB_3[i])*P_min_HOB_3)
        @constraint(m, HOB_3[i]-HOB_3[i .- 1].<=ramp_up_down_HOB_3*OnOffHOB_3[i .- 1]+(1 .- OnOffHOB_3[i .- 1])*P_min_HOB_3)

        @constraint(m, EH[i .- 1]-EH[i].<=ramp_up_down_EH*OnOffEH[i]+(1 .- OnOffEH[i])*P_min_EH)
        @constraint(m, EH[i]-EH[i .- 1].<=ramp_up_down_EH*OnOffEH[i .- 1]+(1 .- OnOffEH[i .- 1])*P_min_EH)
        @constraint(m, EH_2[i .- 1]-EH_2[i].<=ramp_up_down_EH_2*OnOffEH_2[i]+(1 .- OnOffEH_2[i])*P_min_EH_2)
        @constraint(m, EH_2[i]-EH_2[i .- 1].<=ramp_up_down_EH_2*OnOffEH_2[i .- 1]+(1 .- OnOffEH_2[i .- 1])*P_min_EH_2)
        @constraint(m, EH_3[i .- 1]-EH_3[i].<=ramp_up_down_EH_3*OnOffEH_3[i]+(1 .- OnOffEH_3[i])*P_min_EH_3)
        @constraint(m, EH_3[i]-EH_3[i .- 1].<=ramp_up_down_EH_3*OnOffEH_3[i .- 1]+(1 .- OnOffEH_3[i .- 1])*P_min_EH_3)
        @constraint(m, HP[i .- 1]-HP[i].<=ramp_up_down_HP*OnOffHP[i]+(1 .- OnOffHP[i])*P_min_HP)
        @constraint(m, HP[i]-HP[i .- 1].<=ramp_up_down_HP*OnOffHP[i .- 1]+(1 .- OnOffHP[i .- 1])*P_min_HP)
        @constraint(m, HP_2[i .- 1]-HP_2[i].<=ramp_up_down_HP_2*OnOffHP_2[i]+(1 .- OnOffHP_2[i])*P_min_HP_2)
        @constraint(m, HP_2[i]-HP_2[i .- 1].<=ramp_up_down_HP_2*OnOffHP_2[i .- 1]+(1 .- OnOffHP_2[i .- 1])*P_min_HP_2)
        @constraint(m, HP_3[i .- 1]-HP_3[i].<=ramp_up_down_HP_3*OnOffHP_3[i]+(1 .- OnOffHP_3[i])*P_min_HP_3)
        @constraint(m, HP_3[i]-HP_3[i .- 1].<=ramp_up_down_HP_3*OnOffHP_3[i .- 1]+(1 .- OnOffHP_3[i .- 1])*P_min_HP_3)

        @constraint(m, HP_absorption[i .- 1]-HP_absorption[i].<=ramp_up_down_HP_absorption*OnOffHP_absorption[i]+(1 .- OnOffHP_absorption[i])*P_min_HP_absorption)
        @constraint(m, HP_absorption[i]-HP_absorption[i .- 1].<=ramp_up_down_HP_absorption*OnOffHP_absorption[i .- 1]+(1 .- OnOffHP_absorption[i .- 1])*P_min_HP_absorption)

        @constraint(m, CHP_ORC[i .- 1]-CHP_ORC[i].<=ramp_up_down_CHP_ORC*OnOffCHP_ORC[i]+(1 .- OnOffCHP_ORC[i])*P_min_CHP_ORC)
        @constraint(m, CHP_ORC[i]-CHP_ORC[i .- 1].<=ramp_up_down_CHP_ORC*OnOffCHP_ORC[i .- 1]+(1 .- OnOffCHP_ORC[i .- 1])*P_min_CHP_ORC)

        @constraint(m, CHP[i .- 1]-CHP[i].<=ramp_up_down_CHP*OnOffCHP[i]+(1 .- OnOffCHP[i])*P_min_CHP)
        @constraint(m, CHP[i]-CHP[i .- 1].<=ramp_up_down_CHP*OnOffCHP[i .- 1]+(1 .- OnOffCHP[i .- 1])*P_min_CHP)
        @constraint(m, CHP_2[i .- 1]-CHP_2[i].<=ramp_up_down_CHP_2*OnOffCHP_2[i]+(1 .- OnOffCHP_2[i])*P_min_CHP_2)
        @constraint(m, CHP_2[i]-CHP_2[i .- 1].<=ramp_up_down_CHP_2*OnOffCHP_2[i .- 1]+(1 .- OnOffCHP_2[i .- 1])*P_min_CHP_2)

        @constraint(m, HP_waste_heat_I_1[i .- 1]-HP_waste_heat_I_1[i].<=ramp_up_down_HP_waste_heat_I_1*OnOffHP_waste_heat_I_1[i]+(1 .- OnOffHP_waste_heat_I_1[i])*P_min_HP_waste_heat_I_1)
        @constraint(m, HP_waste_heat_I_1[i]-HP_waste_heat_I_1[i .- 1].<=ramp_up_down_HP_waste_heat_I_1*OnOffHP_waste_heat_I_1[i .- 1]+(1 .- OnOffHP_waste_heat_I_1[i .- 1])*P_min_HP_waste_heat_I_1)
        @constraint(m, HP_waste_heat_I_2[i .- 1]-HP_waste_heat_I_2[i].<=ramp_up_down_HP_waste_heat_I_2*OnOffHP_waste_heat_I_2[i]+(1 .- OnOffHP_waste_heat_I_2[i])*P_min_HP_waste_heat_I_2)
        @constraint(m, HP_waste_heat_I_2[i]-HP_waste_heat_I_2[i .- 1].<=ramp_up_down_HP_waste_heat_I_2*OnOffHP_waste_heat_I_2[i .- 1]+(1 .- OnOffHP_waste_heat_I_2[i .- 1])*P_min_HP_waste_heat_I_2)
        @constraint(m, HP_waste_heat_I_3[i .- 1]-HP_waste_heat_I_3[i].<=ramp_up_down_HP_waste_heat_I_3*OnOffHP_waste_heat_I_3[i]+(1 .- OnOffHP_waste_heat_I_3[i])*P_min_HP_waste_heat_I_3)
        @constraint(m, HP_waste_heat_I_3[i]-HP_waste_heat_I_3[i .- 1].<=ramp_up_down_HP_waste_heat_I_3*OnOffHP_waste_heat_I_3[i .- 1]+(1 .- OnOffHP_waste_heat_I_3[i .- 1])*P_min_HP_waste_heat_I_3)

        @constraint(m, HP_waste_heat_II_1[i .- 1]-HP_waste_heat_II_1[i].<=ramp_up_down_HP_waste_heat_II_1*OnOffHP_waste_heat_II_1[i]+(1 .- OnOffHP_waste_heat_II_1[i])*P_min_HP_waste_heat_II_1)
        @constraint(m, HP_waste_heat_II_1[i]-HP_waste_heat_II_1[i .- 1].<=ramp_up_down_HP_waste_heat_II_1*OnOffHP_waste_heat_II_1[i .- 1]+(1 .- OnOffHP_waste_heat_II_1[i .- 1])*P_min_HP_waste_heat_II_1)
        @constraint(m, HP_waste_heat_II_2[i .- 1]-HP_waste_heat_II_2[i].<=ramp_up_down_HP_waste_heat_II_2*OnOffHP_waste_heat_II_2[i]+(1 .- OnOffHP_waste_heat_II_2[i])*P_min_HP_waste_heat_II_2)
        @constraint(m, HP_waste_heat_II_2[i]-HP_waste_heat_II_2[i .- 1].<=ramp_up_down_HP_waste_heat_II_2*OnOffHP_waste_heat_II_2[i .- 1]+(1 .- OnOffHP_waste_heat_II_2[i .- 1])*P_min_HP_waste_heat_II_2)
        @constraint(m, HP_waste_heat_II_3[i .- 1]-HP_waste_heat_II_3[i].<=ramp_up_down_HP_waste_heat_II_3*OnOffHP_waste_heat_II_3[i]+(1 .- OnOffHP_waste_heat_II_3[i])*P_min_HP_waste_heat_II_3)
        @constraint(m, HP_waste_heat_II_3[i]-HP_waste_heat_II_3[i .- 1].<=ramp_up_down_HP_waste_heat_II_3*OnOffHP_waste_heat_II_3[i .- 1]+(1 .- OnOffHP_waste_heat_II_3[i .- 1])*P_min_HP_waste_heat_II_3)

        @constraint(m, HP_waste_heat_III_1[i .- 1]-HP_waste_heat_III_1[i].<=ramp_up_down_HP_waste_heat_III_1*OnOffHP_waste_heat_III_1[i]+(1 .- OnOffHP_waste_heat_III_1[i])*P_min_HP_waste_heat_III_1)
        @constraint(m, HP_waste_heat_III_1[i]-HP_waste_heat_III_1[i .- 1].<=ramp_up_down_HP_waste_heat_III_1*OnOffHP_waste_heat_III_1[i .- 1]+(1 .- OnOffHP_waste_heat_III_1[i .- 1])*P_min_HP_waste_heat_III_1)
        @constraint(m, HP_waste_heat_III_2[i .- 1]-HP_waste_heat_III_2[i].<=ramp_up_down_HP_waste_heat_III_2*OnOffHP_waste_heat_III_2[i]+(1 .- OnOffHP_waste_heat_III_2[i])*P_min_HP_waste_heat_III_2)
        @constraint(m, HP_waste_heat_III_2[i]-HP_waste_heat_III_2[i .- 1].<=ramp_up_down_HP_waste_heat_III_2*OnOffHP_waste_heat_III_2[i .- 1]+(1 .- OnOffHP_waste_heat_III_2[i .- 1])*P_min_HP_waste_heat_III_2)
        @constraint(m, HP_waste_heat_III_3[i .- 1]-HP_waste_heat_III_3[i].<=ramp_up_down_HP_waste_heat_III_3*OnOffHP_waste_heat_III_3[i]+(1 .- OnOffHP_waste_heat_III_3[i])*P_min_HP_waste_heat_III_3)
        @constraint(m, HP_waste_heat_III_3[i]-HP_waste_heat_III_3[i .- 1].<=ramp_up_down_HP_waste_heat_III_3*OnOffHP_waste_heat_III_3[i .- 1]+(1 .- OnOffHP_waste_heat_III_3[i .- 1])*P_min_HP_waste_heat_III_3)

        @constraint(m, HP_waste_heat_seasonal[i .- 1]-HP_waste_heat_seasonal[i].<=ramp_up_down_HP_waste_heat_seasonal*OnOffHP_waste_heat_seasonal[i]+(1 .- OnOffHP_waste_heat_seasonal[i])*P_min_HP_waste_heat_seasonal)
        @constraint(m, HP_waste_heat_seasonal[i]-HP_waste_heat_seasonal[i .- 1].<=ramp_up_down_HP_waste_heat_seasonal*OnOffHP_waste_heat_seasonal[i .- 1]+(1 .- OnOffHP_waste_heat_seasonal[i .- 1])*P_min_HP_waste_heat_seasonal)

        @constraint(m, HP_waste_heat_absorption[i .- 1]-HP_waste_heat_absorption[i].<=ramp_up_down_HP_waste_heat_absorption*OnOffHP_waste_heat_absorption[i]+(1 .- OnOffHP_waste_heat_absorption[i])*P_min_HP_waste_heat_absorption)
        @constraint(m, HP_waste_heat_absorption[i]-HP_waste_heat_absorption[i .- 1].<=ramp_up_down_HP_waste_heat_absorption*OnOffHP_waste_heat_absorption[i .- 1]+(1 .- OnOffHP_waste_heat_absorption[i .- 1])*P_min_HP_waste_heat_absorption)

        ###### Objective function ######

        # Minimize total running costs

        @objective(m,Min, sum(Error)*1e2
                          + sum(HOB)*(cost_var_HOB .+ cost_fuel_HOB/eta_HOB)
                          + sum(HOB_2)*(cost_var_HOB_2 .+ cost_fuel_HOB_2/eta_HOB_2)
                          + sum(HOB_3)*(cost_var_HOB_3 .+ cost_fuel_HOB_3/eta_HOB_3)
                          + sum(EH.*(cost_var_EH .+ cost_fuel_EH/eta_EH))
                          + sum(EH_2.*(cost_var_EH_2 .+ cost_fuel_EH_2/eta_EH_2))
                          + sum(EH_3.*(cost_var_EH_3 .+ cost_fuel_EH_3/eta_EH_3))
                          + sum(HP.*(cost_var_HP .+ cost_fuel_HP./eta_HP))
                          + sum(HP_2.*(cost_var_HP_2 .+ cost_fuel_HP_2./eta_HP_2))
                          + sum(HP_3.*(cost_var_HP_3 .+ cost_fuel_HP_3./eta_HP_3))
                          + sum(HP_absorption.*(cost_var_HP_absorption .+ cost_fuel_HP_absorption./eta_HP_absorption))
                          + sum(CHP)*(cost_var_CHP .+ cost_fuel_CHP/eta_CHP_th)
                          + sum(CHP_2)*(cost_var_CHP_2 .+ cost_fuel_CHP_2/eta_CHP_th_2)
                          + sum(CHP_ORC)*(cost_var_CHP_ORC)
                          + sum(ST)*cost_var_ST
                          + sum(ST_2)*cost_var_ST_2
                          + sum(ST_seasonal)*cost_var_ST_seasonal
                          + sum(HEX)*cost_var_HEX
                          + sum(HP_waste_heat_I_1.*(cost_var_HP_waste_heat_I_1 .+ cost_fuel_HP_waste_heat_I_1./eta_HP_waste_heat_I_1))
                          + sum(HP_waste_heat_I_2.*(cost_var_HP_waste_heat_I_2 .+ cost_fuel_HP_waste_heat_I_2./eta_HP_waste_heat_I_2))
                          + sum(HP_waste_heat_I_3.*(cost_var_HP_waste_heat_I_3 .+ cost_fuel_HP_waste_heat_I_3./eta_HP_waste_heat_I_3))
                          + sum(HP_waste_heat_II_1.*(cost_var_HP_waste_heat_II_1 .+ cost_fuel_HP_waste_heat_II_1./eta_HP_waste_heat_II_1))
                          + sum(HP_waste_heat_II_2.*(cost_var_HP_waste_heat_II_2 .+ cost_fuel_HP_waste_heat_II_2./eta_HP_waste_heat_II_2))
                          + sum(HP_waste_heat_II_3.*(cost_var_HP_waste_heat_II_3 .+ cost_fuel_HP_waste_heat_II_3./eta_HP_waste_heat_II_3))
                          + sum(HP_waste_heat_III_1.*(cost_var_HP_waste_heat_III_1 .+ cost_fuel_HP_waste_heat_III_1./eta_HP_waste_heat_III_1))
                          + sum(HP_waste_heat_III_2.*(cost_var_HP_waste_heat_III_2 .+ cost_fuel_HP_waste_heat_III_2./eta_HP_waste_heat_III_2))
                          + sum(HP_waste_heat_III_3.*(cost_var_HP_waste_heat_III_3 .+ cost_fuel_HP_waste_heat_III_3./eta_HP_waste_heat_III_3))
                          + sum(HP_waste_heat_seasonal.*(cost_var_HP_waste_heat_seasonal .+ cost_fuel_HP_waste_heat_seasonal./eta_HP_waste_heat_seasonal))
                          + sum(HP_waste_heat_absorption.*(cost_var_HP_waste_heat_absorption .+ cost_fuel_HP_waste_heat_absorption./eta_HP_waste_heat_absorption))
                          - sum(CHP*cb_CHP.*CHP_electricity_price)
                          - sum(CHP_2*cb_CHP_2.*CHP_electricity_price_2)
                          - sum(CHP_ORC*cb_CHP_ORC.*CHP_ORC_electricity_price))

        print_function("Julia code for district heating: solve(m)...")

        optimize!(m)
        status = termination_status(m)

        print_function(string("Exit with status:", status))
        print_function("Julia code for district heating: END OF SCRIPT!")
        if primal_status(m) == MOI.NO_SOLUTION
            print_function("No solution found.")
            return
        end

        ##### Result visualization ####

        Result_HOB=value.(HOB)
        Result_HOB_2=value.(HOB_2)
        Result_HOB_3=value.(HOB_3)

        Result_EH=value.(EH)
        Result_EH_2=value.(EH_2)
        Result_EH_3=value.(EH_3)

        Result_HP=value.(HP)
        Result_HP_2=value.(HP_2)
        Result_HP_3=value.(HP_3)
        Result_HP_absorption=value.(HP_absorption)

        Result_CHP=value.(CHP)
        Result_CHP_2=value.(CHP_2)

        Result_CHP_ORC=value.(CHP_ORC)

        Result_ST=value.(ST)
        Result_ST_2=value.(ST_2)
        Result_ST_seasonal=value.(ST_seasonal)

        Result_HEX=value.(HEX)

        Result_HP_waste_heat_I_1=value.(HP_waste_heat_I_1)
        Result_HP_waste_heat_I_2=value.(HP_waste_heat_I_2)
        Result_HP_waste_heat_I_3=value.(HP_waste_heat_I_3)

        Result_HP_waste_heat_II_1=value.(HP_waste_heat_II_1)
        Result_HP_waste_heat_II_2=value.(HP_waste_heat_II_2)
        Result_HP_waste_heat_II_3=value.(HP_waste_heat_II_3)

        Result_HP_waste_heat_III_1=value.(HP_waste_heat_III_1)
        Result_HP_waste_heat_III_2=value.(HP_waste_heat_III_2)
        Result_HP_waste_heat_III_3=value.(HP_waste_heat_III_3)

        Result_HP_waste_heat_seasonal=value.(HP_waste_heat_seasonal)
        Result_HP_waste_heat_absorption=value.(HP_waste_heat_absorption)

        Result_SOC=value.(SOC)
        Result_SOC_seasonal=value.(SOC_seasonal)
        Result_Error=value.(Error)

        #load duration for each technology #

        sort_Result_HOB=sort(Result_HOB, rev=true)
        sort_Result_HOB_2=sort(Result_HOB_2, rev=true)
        sort_Result_HOB_3=sort(Result_HOB_3, rev=true)
        sort_Result_EH=sort(Result_EH, rev=true)
        sort_Result_EH_2=sort(Result_EH_2, rev=true)
        sort_Result_EH_3=sort(Result_EH_3, rev=true)
        sort_Result_HP=sort(Result_HP, rev=true)
        sort_Result_HP_2=sort(Result_HP_2, rev=true)
        sort_Result_HP_3=sort(Result_HP_3, rev=true)
        sort_Result_HP_absorption=sort(Result_HP_absorption, rev=true)
        sort_Result_CHP=sort(Result_CHP, rev=true)
        sort_Result_CHP_2=sort(Result_CHP_2, rev=true)
        sort_Result_CHP_ORC=sort(Result_CHP_ORC, rev=true)
        sort_Result_ST=sort(Result_ST, rev=true)
        sort_Result_ST_2=sort(Result_ST_2, rev=true)
        sort_Result_ST_seasonal=sort(Result_ST_seasonal, rev=true)
        sort_Result_HEX=sort(Result_HEX, rev=true)
        sort_Result_HP_waste_heat_I_1=sort(Result_HP_waste_heat_I_1, rev=true)
        sort_Result_HP_waste_heat_I_2=sort(Result_HP_waste_heat_I_2, rev=true)
        sort_Result_HP_waste_heat_I_3=sort(Result_HP_waste_heat_I_3, rev=true)
        sort_Result_HP_waste_heat_seasonal=sort(Result_HP_waste_heat_seasonal, rev=true)
        sort_Result_HP_waste_heat_absorption=sort(Result_HP_waste_heat_absorption, rev=true)
        sort_Result_Error=sort(Result_Error, rev=true)

        sort_Result_HP_waste_heat_II_1=sort(Result_HP_waste_heat_II_1, rev=true)
        sort_Result_HP_waste_heat_II_2=sort(Result_HP_waste_heat_II_2, rev=true)
        sort_Result_HP_waste_heat_II_3=sort(Result_HP_waste_heat_II_3, rev=true)

        sort_Result_HP_waste_heat_III_1=sort(Result_HP_waste_heat_III_1, rev=true)
        sort_Result_HP_waste_heat_III_2=sort(Result_HP_waste_heat_III_2, rev=true)
        sort_Result_HP_waste_heat_III_3=sort(Result_HP_waste_heat_III_3, rev=true)

        ## Export results ##
        writedlm(string(output_folder, "\\Result_HOB_",network_id ,".csv"), Result_HOB, csv_sep)
        writedlm(string(output_folder, "\\Result_HOB_2_",network_id ,".csv"), Result_HOB_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HOB_3_",network_id ,".csv"), Result_HOB_3, csv_sep)

        writedlm(string(output_folder, "\\Result_EH_",network_id ,".csv"), Result_EH, csv_sep)
        writedlm(string(output_folder, "\\Result_EH_2_",network_id ,".csv"), Result_EH_2, csv_sep)
        writedlm(string(output_folder, "\\Result_EH_3_",network_id ,".csv"), Result_EH_3, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_1_",network_id ,".csv"), Result_HP, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_2_",network_id ,".csv"), Result_HP_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_3_",network_id ,".csv"), Result_HP_3, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_absorption_",network_id ,".csv"), Result_HP_absorption, csv_sep)

        writedlm(string(output_folder, "\\Result_CHP_",network_id ,".csv"), Result_CHP, csv_sep)
        writedlm(string(output_folder, "\\Result_CHP_2_",network_id ,".csv"), Result_CHP_2, csv_sep)

        writedlm(string(output_folder, "\\Result_CHP_ORC_",network_id ,".csv"), Result_CHP_ORC, csv_sep)

        writedlm(string(output_folder, "\\Result_ST_",network_id ,".csv"), Result_ST, csv_sep)
        writedlm(string(output_folder, "\\Result_ST_2_",network_id ,".csv"), Result_ST_2, csv_sep)
        writedlm(string(output_folder, "\\Result_ST_seasonal_",network_id ,".csv"), Result_ST_seasonal, csv_sep)

        writedlm(string(output_folder, "\\Result_HEX_",network_id ,".csv"), Result_HEX, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_I_1_",network_id ,".csv"), Result_HP_waste_heat_I_1, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_I_2_",network_id ,".csv"), Result_HP_waste_heat_I_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_I_3_",network_id ,".csv"), Result_HP_waste_heat_I_3, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_II_1_",network_id ,".csv"), Result_HP_waste_heat_II_1, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_II_2_",network_id ,".csv"), Result_HP_waste_heat_II_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_II_3_",network_id ,".csv"), Result_HP_waste_heat_II_3, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_III_1_",network_id ,".csv"), Result_HP_waste_heat_III_1, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_III_2_",network_id ,".csv"), Result_HP_waste_heat_III_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_III_3_",network_id ,".csv"), Result_HP_waste_heat_III_3, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_waste_heat_seasonal_",network_id ,".csv"), Result_HP_waste_heat_seasonal, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_waste_heat_absorption_",network_id ,".csv"), Result_HP_waste_heat_absorption, csv_sep)
        writedlm(string(output_folder, "\\Result_DEM_",network_id ,".csv"), DEM, csv_sep)
        writedlm(string(output_folder, "\\Result_SOC_",network_id ,".csv"), Result_SOC, csv_sep)
        writedlm(string(output_folder, "\\Result_SOC_seasonal_",network_id ,".csv"), Result_SOC_seasonal, csv_sep)
        writedlm(string(output_folder, "\\Result_Error_",network_id ,".csv"), Result_Error, csv_sep)
    end

    export district_solver

    function write_results_example(output_folder, network_id, example_file)
        writedlm(string(output_folder, "\\Result_HOB_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HOB_2_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HOB_3_",network_id ,".csv"), example_file, csv_sep)

        writedlm(string(output_folder, "\\Result_EH_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_EH_2_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_EH_3_",network_id ,".csv"), example_file, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_2_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_3_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_absorption_",network_id ,".csv"), example_file, csv_sep)

        writedlm(string(output_folder, "\\Result_CHP_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_CHP_2_",network_id ,".csv"), example_file, csv_sep)

        writedlm(string(output_folder, "\\Result_CHP_ORC_",network_id ,".csv"), example_file, csv_sep)

        writedlm(string(output_folder, "\\Result_ST_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_ST_2_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_ST_seasonal_",network_id ,".csv"), example_file, csv_sep)

        writedlm(string(output_folder, "\\Result_HEX_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_I_1_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_I_2_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_I_3_",network_id ,".csv"), example_file, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_II_1_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_II_2_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_II_3_",network_id ,".csv"), example_file, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_III_1_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_III_2_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_III_3_",network_id ,".csv"), example_file, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_waste_heat_seasonal_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_waste_heat_absorption_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_DEM_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_SOC_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_SOC_seasonal_",network_id ,".csv"), example_file, csv_sep)
        writedlm(string(output_folder, "\\Result_Error_",network_id ,".csv"), example_file, csv_sep)
    end
end

max_stars_repo_name: Planheat/Planheat-Tool

max_issues_repo_name: Planheat/Planheat-Tool

max_forks_repo_name: Planheat/Planheat-Tool


text:
module Decoder

using DSP, FFTW
using Match

export decode

# Returns two vectors containing
# the indexes where samples with
# information (aka: different from silences)
# start and end
#
function find_data_samples(x, window, get_threshold)

	if isempty(x)
		return [], []
	end

	# Normalize samples
	x = x ./ maximum(abs.(x))

	# Calculate the envolvent
	y = conv(window, x)

	# Keep the samples beginning in N/2
	# because the convolution shift it in that amount
	y1 = y[div(length(window), 2):end]

	# Threshold
	thres = ones(length(y1)) .* get_threshold(y1)

	# Detecting tones and silences
	pulses = map(x -> Int(x), y1 .> thres)

	# To detect where the
	# start and end of each pulse are
	differ = diff(pulses)

	# Indexes where the pulses start
	pulses_start = map(x -> x + 1, findall(x -> x == 1, differ))

	# Indexes where the pulses end
	pulses_end = findall(x -> x == -1, differ)

	return pulses_start, pulses_end
end

function get_time_window(fs)

	# The samples per period
	# The longest period is for the smallest freq
	# In this case: 697 Hz ~ 500 Hz
	samples_period = Int(0.005 * fs)

	# Window size is 3 periods
	N = samples_period * 3

	h = ones(N) ./ N

	return h
end

function get_pulse_window()

	N = 120
	h = tukey(N, 0.5) ./ (N/4)

	return h
end

function get_freq_window()

	N = 10
	h = tukey(N, 0.5) ./ (N/2)

	return h
end

function time_threshold(y)
	return 0.0125
end

function pulse_threshold(y)
	return 0.85 * maximum(y)
end

function freq_threshold(y)
	return 0.15 * maximum(y)
end

# Given a frequency it matches to a
# (start, end) interval it returns
# the frequency that belongs to it
get_freq(f) = @match f begin
	600:733 => 697
	744:811 => 770
	812:896 => 852
	897:1020 => 941
	1100:1272 => 1209
	1273:1406 => 1336
	1407:1555 => 1477
	1556:1700 => 1633
	_ => 0
end

get_symbol(f1, f2) = @match (f1, f2) begin

	(941, 1336) || (1336, 941) => "0"
	(697, 1209) || (1209, 697) => "1"
	(697, 1336) || (1336, 697) => "2"
	(697, 1477) || (1477, 697) => "3"
	(770, 1209) || (1209, 770) => "4"
	(770, 1336) || (1336, 770) => "5"
	(770, 1477) || (1477, 770) => "6"
	(852, 1209) || (1209, 852) => "7"
	(852, 1336) || (1336, 852) => "8"
	(852, 1477) || (1477, 852) => "9"
	(941, 1209) || (1209, 941) => "*"
	(941, 1477) || (1477, 941) => "#"
	(697, 1633) || (1633, 697) => "A"
	(770, 1633) || (1633, 770) => "B"
	(852, 1633) || (1633, 852) => "C"
	(941, 1633) || (1633, 941) => "D"
	_ => ""
end

# Receives a list of tuples containing
# the start and end of the interval
# that corresponds to a given frequency
#
# Returns the symbol it belongs to
function decode_symbol(freqs)

	if isempty(freqs)
		return ""
	end

	freqs = map(f -> get_freq(div(f[1] + f[2], 2)), freqs)

	if length(freqs) < 2
		return ""
	end

	return get_symbol(freqs[1], freqs[2])
end

# Find frequency peaks and filter
# only the ones that are not in the range
# of dialing tones
function analyze_freqs(X, fs)

	freqs = []

	if isempty(X)
		return freqs
	end

	starts, ends = find_data_samples(X,
					 get_freq_window(),
					 freq_threshold)

	samples_per_freq = Int(floor(fs/(2*length(X))))

	starts = starts .* samples_per_freq
	ends = ends .* samples_per_freq

	for i in range(1, length(starts))

		# Filter frequencies only in the
		# range of (500 Hz, 1700 Hz)
		if (500 <= starts[i] <= 1700 &&
			ends[i] <= 1800 &&
			1 <= ends[i] - starts[i] < 250)

			push!(freqs, (starts[i], ends[i]))
		end
	end

	return freqs
end

function analyze_pulse(pulse, fs)

	e = pulse .^ 2

	starts, ends = find_data_samples(e,
					 get_pulse_window(),
					 pulse_threshold)

	symbols = []

	for i in range(1, length(starts))

		# Because of the convolution with
		# the window, in some cases the 'ends'
		# indexes may be outside the range of
		# the pulse.
		if ends[i] > length(pulse)
			p = pulse[starts[i]:end]
		else
			p = pulse[starts[i]:ends[i]]
		end

		X = abs.(fft(p))
		X = X[1:div(end,2)]

		freqs = analyze_freqs(X, fs)
		sym = decode_symbol(freqs)

		append!(symbols, sym)
	end

	return symbols
end

function decode(s, fs)

	s = vec(s)

	e = s .^ 2

	symbols = []

	starts, ends = find_data_samples(e,
					 get_time_window(fs),
					 time_threshold)

	# Analyze all data chunks looking for
	# the "dialing" frequencies
	for i in range(1, length(starts))

		pulse = s[starts[i]:ends[i]]

		syms = analyze_pulse(pulse, fs)

		append!(symbols, syms)
	end

	return symbols
end

end # module

ext: jl

max_stars_repo_name: PatricioIribarneCatella/TP-Telefono

max_issues_repo_name: PatricioIribarneCatella/TP-Telefono

max_forks_repo_name: PatricioIribarneCatella/TP-Telefono


text:
mutable struct CachedCovariate{T<:AbstractFloat} <: AbstractCovariate{T}
    cache::Union{Vector{T}, Nothing}
    basecovariate::AbstractCovariate{T}
    lockobj::Threads.TatasLock
end

Base.length(var::CachedCovariate) = length(var.basecovariate)

getname(var::CachedCovariate) = getname(var.basecovariate)

function CachedCovariate(basecovariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    CachedCovariate{T}(nothing, basecovariate, Threads.TatasLock())
end

function cache(basecovariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    CachedCovariate{T}(nothing, basecovariate, Threads.TatasLock())
end

function slice(covariate::CachedCovariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:AbstractFloat}
    basecov = covariate.basecovariate
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    if isa(basecov, CachedCovariate) || isa(basecov, Covariate)
        slice(basecov, fromobs, toobs, slicelength)
    else
        lockobj = covariate.lockobj
        lock(lockobj)
        try
            if covariate.cache === nothing
                slices = slice(basecov, 1, length(basecov), SLICELENGTH)
                cachedata = Vector{T}(undef, length(basecov))
                fold(0, slices) do offset, slice
                    n = length(slice)
                    view(cachedata, (1 + offset):(n + offset)) .= slice
                    offset += n
                end
                covariate.cache = cachedata
            end
        finally
            unlock(lockobj)
        end
        slice(covariate.cache, fromobs, toobs, slicelength)
    end
end

function Base.map(covariate::CachedCovariate, dataframe::AbstractDataFrame)
    CachedCovariate(map(covariate.basecovariate, dataframe))
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
using Documenter, GenTransportMCMC

makedocs(
  sitename="GenTransportMCMC",
  modules=[GenTransportMCMC],
  pages=[
    "Home" => "index.md",
    "Getting Started" => "getting_started.md",
    "Mathematical Details" => "mathematical_details.md",
    "API" => [
      "Inference" => "ref/inference.md",
      "Transport" => "ref/transport.md",
    ],
  ],
)

deploydocs(
  repo = "github.com/mvarble/GenTransportMCMC.jl.git",
  target = "build",
  devbranch = "main",
)

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
const _SELF_ROOT = dirname(@__DIR__)

_load_gitignore() = read(joinpath(_SELF_ROOT, ".gitignore"), String)
_default_user() = LibGit2.getconfig("user.name", "")

_default_plugins() = [
    # allowed
    PkgTemplates.Git(ignore = [_load_gitignore()], manifest = true, branch = "main"),
    PkgTemplates.ProjectFile(),
    PkgTemplates.SrcDir(),
    PkgTemplates.Tests(),
    PkgTemplates.Readme(),
    PkgTemplates.License(),
    PkgTemplates.Codecov(),
    PkgTemplates.GitHubActions(;coverage = true),

    # disallowed
    !PkgTemplates.CompatHelper,
    !PkgTemplates.TagBot
]

function _cp(src_, dst_; kwargs...)
    mkpath(dirname(src_))
    mkpath(dirname(dst_))
    cp(src_, dst_; kwargs...)
end

function _generate(pkgname;
        user = _default_user(), julia = v"1.6.0"
    )

    dir = Pkg.devdir()
    pkgdir = joinpath(dir, replace(pkgname, ".jl" => ""))

    plugins = _default_plugins()
    t = PkgTemplates.Template(;user, julia, plugins, dir)
    PkgTemplates.generate(t, pkgname)

    # copy tagged-release.yml
    _cp(
        joinpath(_SELF_ROOT, ".github/workflows/tagged-release.yml"),
        joinpath(pkgdir, ".github/workflows/tagged-release.yml");
        force = true
    )

    return pkgdir
end

function run_generate_pkg(argv=ARGS)

    argset = ArgParse.ArgParseSettings()
    ArgParse.@add_arg_table! argset begin
        "pkgname"
            help="The package name"
            required = true
        "--code", "-c"
            help = "Will open the project after creation"
            action = :store_true
    end

    parsed_args = ArgParse.parse_args(argv, argset)
    pkgname = parsed_args["pkgname"]
    code = parsed_args["code"]

    ## ---------------------------------------------------------
    pkgdir = _generate(pkgname)

    ## ---------------------------------------------------------
    if code
        run(Cmd(["code", pkgdir]); wait=true)
    end

end


text:
using Atum
using Atum.Euler

using Test
using StaticArrays: SVector
using LinearAlgebra: norm

function vortex(law, x⃗, t)
  FT = eltype(law)
  c⃗ = SVector(FT(5), FT(0))
  u⃗₀ = SVector(FT(1), FT(0))
  β = FT(5)
  r⃗ = x⃗ - c⃗ - u⃗₀ * t
  f = β * exp(1 - r⃗' * r⃗)

  ρ = (1 - (γ(law) - 1) * f ^ 2 / (16 * γ(law) * π ^ 2)) ^ (1 / (γ(law) - 1))
  u⃗ = u⃗₀ + f / 2π * SVector(-r⃗[2], r⃗[1])
  ρu⃗ = ρ * u⃗
  p = ρ ^ γ(law)
  ρe = Euler.energy(law, ρ, ρu⃗, p)

  SVector(ρ, ρu⃗..., ρe)
end

function run(A, FT, N, K)
  Nq = N + 1

  law = EulerLaw{FT, 2}()

  cell = LobattoCell{FT, A}(Nq, Nq)
  vx = range(FT(0), stop=FT(20), length=K+1)
  vy = range(FT(-5), stop=FT(5), length=K+1)
  grid = brickgrid(cell, (vx, vy); periodic=(true, true))

  dg = DGSEM(; law, cell, grid, numericalflux = RusanovFlux())

  cfl = FT(1 // 4)
  dt = cfl * step(vy) / N
  timeend = FT(5)

  q = vortex.(Ref(law), points(grid), FT(0))

  odesolver = LSRK54(dg, q, dt)
  solve!(q, timeend, odesolver)

  qexact = vortex.(Ref(law), points(grid), timeend)
  errf = map(components(q), components(qexact)) do f, fexact
    sqrt(sum(dg.MJ .* (f .- fexact) .^ 2))
  end
  norm(errf)
end

let
  A = Array
  FT = Float64
  N = 4

  nlevels = 4
  errors = zeros(FT, nlevels)

  for l in 1:nlevels
    K = 5 * 2 ^ (l - 1)
    errf = run(A, FT, N, K)
    errors[l] = errf
  end
  rates = log2.(errors[1:(nlevels-1)] ./ errors[2:nlevels])
  @test rates[end] ≈ N + 1 atol = 0.25
end


text:
import InstrumentConfig: initialize, terminate

"""
An instrument is a generic device with which you can take and read measurements

All instruments can be initialized and terminated. For more
information on how to connect to a supported instrument:
```
help> initialize
```

You can use the help feature on any supported instrument group.

For more information on the actual instruments you can initialize use julia's help feature on one of the supported instrument groups.

For example lets say you want to learn about oscilloscopes.
```
julia> ?
help> Oscilloscope

  Supported Instruments:
  ≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡

    •  AgilentDSOX4034A
```

Pick one of the supported instruments for more information on
how to use it and for all its available functions:
```
help> AgilentDSOX4034A
```

# Supported Instrument Groups
- `Oscilloscope`
- `MultiMeter`
- `PowerSupply`
- `WaveformGenerator`
- `ImpedanceAnalyzer`
- `XYZStage`

"""
abstract type Instrument end
"""
# Supported Instruments
- [`AgilentDSOX4024A`](@ref)
- `AgilentDSOX4034A`
"""
abstract type Oscilloscope <: Instrument end

"""
# Supported Instruments
- KeysightDMM34465A
"""
abstract type MultiMeter <: Instrument end

"""
# Supported Instruments
- `AgilentE36312A`
- `VersatilePower`
- `PS310`
"""
abstract type PowerSupply <: Instrument end

# Maybe store ips in a config file and it dynamically shows you address?
"""
# Supported Instruments:
- `Keysight33612A`: Default ip ~ "10.1.30.36"
"""
abstract type WaveformGenerator <: Instrument end

"""
# Supported Instruments:
- `Agilent4294A`
- `Agilent4395A`

"""
abstract type ImpedanceAnalyzer <: Instrument end

mutable struct Instr{ T <: Instrument } <: Instrument
    model::Union{Type, T}
    address::String
    sock::TCPSocket
    connected::Bool
end

function CreateTcpInstr(model, address)
    Instr{model}(model, address, TCPSocket(), false)
end

function Base.show(io::IO, ::MIME"text/plain", i::TcpInstruments.Instr)
    model = i.model isa DataType ? i.model : typeof(i.model)
    println("TcpInstruments.Instr{$(i.model)}")
    println("    Group: $(supertype(model))")
    println("    Model: $(model)")
    println("    Address: $(i.address)")
    println("    Connected: $(i.connected)")
end

"""
    initialize(model, address)

Initializes a connection to the instrument at the given (input) IP address.

# Arguments
- `model`: They device type you are connecting to. Use `help> Instrument` to see available options.
- `address::String`: The ip address of the device. Ex. "10.3.30.23"

# Keywords
- `GPIB_ID::Int`: The GPIB interface ID of your device. This is optional and doesn't need to be set unless you are using a prologix controller to control it remotely.
"""
function initialize(model::Type{T}, address; GPIB_ID=-1) where T <: Instrument
    instr_h = CreateTcpInstr(model, address)
    connect!(instr_h)
    remote_mode(instr_h)
    if GPIB_ID >= 0
        set_prologix_chan(instr_h, GPIB_ID)
    end
    return instr_h
end

function initialize(model::Type{T}) where T <: Instrument
    data = nothing
    try
        data = get_config()[string(model)]
    catch e
        error("""
        $(string(model)) was not found in your .tcp_instruments.yml file.
        To update to the latest version:
        `create_config()`

        Otherwise please add it to your config file or
        specify an ip address:
        `initialize($(string(model)), "10.1.30.XX")`

        """)
    end
    if data isa String
        return initialize(model, data)
    end

    address = get(data, "address", "")
    gpib = get(data, "gpib", "")

    if isempty(gpib)
        return initialize(model, address)
    end
    return initialize(model, address, GPIB_ID=gpib)
end

"""
    terminate(instr)

Closes the TCP connection.
"""
function terminate(instr::Instrument)
    close!(instr)
    local_mode(instr)
end

reset(obj) = write(obj, "*RST")

remote_mode(obj)   = nothing
local_mode(obj) = nothing
set_prologix_chan(obj, chan) = write(obj, "++addr $chan")
get_prologix_chan(obj) = query(obj, "++addr")
"""
    info(instr_h)
Asks an instrument to print model number and other device info.
"""
info(obj) = query(obj, "*IDN?")

import Base.write, Base.read

function connect!(instr::Instrument)
    instr.connected && error("Cannot connect. Instrument is already connected!")
	SCPI_port = 5025
	host,port = split_str_into_host_and_port(instr.address)
	port == 0 && (port = SCPI_port)
	instr.sock = connect(host,port)
	instr.connected = true
end

function close!(instr::Instrument)::Bool
    !instr.connected && error("Cannot disconnect. Instrument is not connected!")
	close(instr.sock)
	instr.connected = false
end

function write(instr::Instrument, message::AbstractString)
    !instr.connected && error("Instrument is not connected, cannot write to it!")
	println(instr.sock, message)
end

function read(instr::Instrument)
    !instr.connected && error("Instrument is not connected, cannot read from it!")
	return rstrip(readline(instr.sock), ['\r', '\n'])
end

"""
Writes a message to a device then listens for and returns any output from
the device.

This is a blocking procedure and will block until a response is received from the device or
till it has been blocking for longer than the designated `timeout` time after which an
error will be thrown.

# Arguments
- `instr::Instrument`: Any instrument that supports being written to and read from
- `message::AbstractString`: The message to be sent to the device before listening for a response
- `timeout`: _Optional flag_ ~ How long to try and listen for a response before giving up and throwing an error. The default time is 2.8 seconds. _Note_: if timeout is set to 0 then this will turn off the timeout functionality and `query` may listen/block indefinitely for a response
"""
function query(instr::Instrument, message::AbstractString; timeout=2.8)
    write(instr, message)
    if timeout == 0
        return read(instr)
    end
    proc = @spawn read(instr)
    start_clock = time()
    while (time() - start_clock) < timeout
        if proc.state != :runnable
            break
        end
        sleep(0.05)
    end
    if proc.state == :runnable
        schedule(proc, ErrorException("Query timed out"), error=true)
        error("Query timed out")
    end
    return fetch(proc)
end

"""
Differs from [query](@ref) in that it will return a Float64 and not a String
"""
f_query(obj, ins; timeout=0.5) = parse(Float64, query(obj, ins; timeout=timeout))

"""
Differs from [query](@ref) in that it will return an Int64 and not a String
"""
i_query(obj, ins; timeout=0.5) = parse(Int64, query(obj, ins; timeout=timeout))


text:
using GLVisualize, GLAbstraction

if !isdefined(:runtests)
	window = glscreen()
    timesignal = loop(linspace(0f0,1f0,360))
end

# create a glsl fragment shader
parametric_func = frag"""
    uniform float arg1; // you can add arbitrary uniforms and supply them via the keyword args
    float function(float x) {
     return arg1*sin(1/tan(x));
   }
"""
# _view the function on a 1700x800 pixel plane
paremetric = visualize(parametric_func, arg1=timesignal, dimensions=(1700, 800))
_view(paremetric, window, camera=:orthographic_pixel)

if !isdefined(:runtests)
	renderloop(window)
end


text:
# High level MCMC and AIS inference functions.  These functions
# construct the star and galaxy log likelihood and prior functions, and
# run either Slice-Sampling or Slice-Sampling-within-AIS

"""
Runs Annealed Importance Sampling on both the Star and Galaxy posteriors
to estimate the marginal likelihood of each model and posterior samples
simultaneously
"""
function run_ais(entry::CatalogEntry,
                 imgs::Vector,
                 patches::Array{ImagePatch, 2},
                 background_images::Array{Array{Float64, 2}, 1},
                 pos_delta::Array{Float64, 1}=[2., 2.];
                 num_samples::Int=2,
                 num_temperatures::Int=50,
                 print_skip::Int=20,
                 num_samples_per_chain::Int=25)
    Log.info("\nRunning AIS on with patch size $(imgs[1].H) x $(imgs[1].W)")
    Log.info("  catalog type: $(entry.is_star ? "star" : "galaxy")")
    Log.info("  num images  : $(length(imgs))")
    n_active = sum([sum(p.active_pixel_bitmap) for p in patches[1, :]]) / length(patches[1,:])
    Log.info("  num active pixels per patch: $(n_active)")
    #imgs, pos_delta, num_samples, print_skip, num_temperatures, num_samples_per_chain = patch_images, [1., 1.], 3, 1, 10, 25

    ###################
    # run star MCMC   #
    ###################
    star_loglike, star_logprior, star_logpost, sample_star_prior, ra_lim, dec_lim, uniform_to_deg, deg_to_uniform =
          MCMC.make_star_inference_functions(imgs, entry;
                                        patches=patches[1, :],
                                        background_images=background_images,
                                        pos_delta=pos_delta)

    th_cat = [log.(entry.star_fluxes)..., deg_to_uniform(entry.pos)...]
    th_rand = sample_star_prior()
    Log.info("star loglike at CATALOG vs PRIOR : ", star_loglike(th_cat), ", ", star_loglike(th_rand))
    Log.info("star logprior at CATALOG vs PRIOR : ", star_logprior(th_cat), ", ", star_logprior(th_rand))
    star_schedule = MCMC.sigmoid_schedule(num_temperatures; rad=4)
    res_star = MCMC.ais_slicesample(star_logpost, star_logprior,
                                    sample_star_prior;
                                    schedule=star_schedule,
                                    num_samps=num_samples,
                                    num_samples_per_step=1)
    star_chains, star_chain_lls = [], []
    for i in 1:size(res_star[:zsamps], 2)
        Log.info("    star chain ", i, " of ", size(res_star[:zsamps], 2))
        star_chain, star_lls = MCMC.slicesample_chain(star_logpost,
            res_star[:zsamps][:,1], num_samples_per_chain;
            print_skip=Int64(num_samples_per_chain/5), verbose=false)
        push!(star_chains, star_chain)
        push!(star_chain_lls, star_lls)
    end
    res_star[:zsamps]    = transpose(vcat(star_chains...))
    res_star[:zsamp_lls] = vcat(star_chain_lls...)
    lnZ  = res_star[:lnZ]
    lnZs = res_star[:lnZ_bootstrap]
    lo, hi = percentile(lnZs, 2.5), percentile(lnZs, 97.5)
    Log.info(@sprintf "STAR AIS estimate : %6.3f [%6.3f, %6.3f]\n" lnZ lo hi)
    Log.info(@sprintf "  CI width : %6.5f \n" (hi-lo))

    ####################
    # run galaxy AIS   #
    ####################
    gal_loglike, gal_logprior, gal_logpost, sample_gal_prior, ra_lim, dec_lim, uniform_to_deg, deg_to_uniform =
          MCMC.make_gal_inference_functions(imgs, entry;
                                        patches=patches[1, :],
                                        background_images=background_images,
                                        pos_delta=pos_delta)
    th_cat = MCMC.parameters_from_catalog(entry; is_star=false)
    th_rand = sample_gal_prior()
    Log.info("gal loglike  at CATALOG vs PRIOR : ", gal_loglike(th_cat), ", ", gal_loglike(th_rand))
    Log.info("gal logprior at CATALOG vs PRIOR : ", gal_logprior(th_cat), ", ", gal_logprior(th_rand))
    gal_schedule = MCMC.sigmoid_schedule(num_temperatures; rad=4)
    res_gal = MCMC.ais_slicesample(gal_logpost, gal_logprior,
                                   sample_gal_prior;
                                   schedule=gal_schedule,
                                   num_samps=num_samples,
                                   num_samples_per_step=1)
    gal_chains, gal_chain_lls = [], []
    for i in 1:size(res_gal[:zsamps], 2)
        Log.info("    gal chain ", i, " of ", size(res_gal[:zsamps], 2))
        gal_chain, gal_lls = MCMC.slicesample_chain(gal_logpost,
            res_gal[:zsamps][:,1], num_samples_per_chain;
            print_skip=Int64(num_samples_per_chain/5), verbose=false)
        push!(gal_chains, gal_chain)
        push!(gal_chain_lls, gal_lls)
    end
    res_gal[:zsamps]    = transpose(vcat(gal_chains...))
    res_gal[:zsamp_lls] = vcat(gal_chain_lls...)
    lnZ = res_gal[:lnZ]
    lnZs = res_gal[:lnZ_bootstrap]
    lo, hi = percentile(lnZs, 2.5), percentile(lnZs, 97.5)
    Log.info(@sprintf "GAL AIS estimate : %6.3f [%6.3f, %6.3f]\n" lnZ lo hi)
    Log.info(@sprintf "  CI width : %6.5f \n" (hi-lo))

    #########################################################
    # Compute prob star vs gal based on marginal likelihood #
    #########################################################
    type_chain  = zeros(length(res_gal[:lnZ_bootstrap]))
    # is_star = [.28, .72] vs [.999, .001] vs [.5, .5]
    lnprob_a    = log(.28)
    lnprob_nota = log(.72)
    for n in 1:length(res_gal[:lnZ_bootstrap])
        # normalizing constant is ln p(data | star)
        # so p(star | data) \propto p(data | star) p(star)
        lnprob_star = res_star[:lnZ_bootstrap][n] + lnprob_a
        lnprob_gal  = res_gal[:lnZ_bootstrap][n] + lnprob_nota
        lnsum = Model.logsumexp([lnprob_star, lnprob_gal]) # normalize
        type_chain[n] = lnprob_star - lnsum
    end
    ave_pstar = Model.logsumexp(type_chain) - log(length(type_chain))
    Log.info("  source p-star = ", exp(ave_pstar))

    ####################################################################
    # convert positions to RA/Dec and organize chains into dataframes  #
    ####################################################################
    for n in 1:size(res_gal[:zsamps], 2)
        res_gal[:zsamps][6:7,n]  = uniform_to_deg(res_gal[:zsamps][6:7,n])
        res_star[:zsamps][6:7,n] = uniform_to_deg(res_star[:zsamps][6:7,n])
    end
    star_chain = MCMC.samples_to_dataframe(transpose(res_star[:zsamps]), is_star=true)
    gal_chain  = MCMC.samples_to_dataframe(transpose(res_gal[:zsamps]), is_star=false)

    # store objid (for concatenation)
    mcmc_results = Dict("star_samples" => star_chain,
                        "star_lls"     => res_star[:zsamp_lls], #[res_star[:lnZsamps],
                        "star_bootstrap"=> res_star[:lnZ_bootstrap],
                        "gal_samples"  => gal_chain,
                        "gal_lls"      => res_gal[:zsamp_lls], #res_gal[:lnZsamps],
                        "gal_bootstrap" => res_gal[:lnZ_bootstrap],
                        "type_samples" => type_chain,
                        "ave_pstar"    => ave_pstar)
    return mcmc_results
end

"""
Run single source MCMC chain
"""
function run_mcmc(entry::CatalogEntry,
                  imgs::Vector,
                  patches::Array{ImagePatch, 2},
                  background_images::Array{Array{Float64, 2}, 1},
                  pos_delta::Array{Float64, 1}=[1., 1.];
                  num_samples::Int=500,
                  print_skip::Int=20)
    Log.info("\nRunning mcmc on entry with patch size ",
            imgs[1].H, "x", imgs[1].W)
    Log.info("  catalog type: ", entry.is_star ? "star" : "galaxy")

    # position log prior --- same for both star and galaxy (constrains to a
    # small window around the existing catalog location
    #imgs, pos_delta, num_samples, print_skip = patch_images, [1., 1.], 200, 20
    pos_logprior, ra_lim, dec_lim = MCMC.make_location_prior(
        imgs[1], entry.pos; pos_pixel_delta=[2., 2.])

    ###################
    # run star MCMC   #
    ###################
    star_loglike, constrain_pos, unconstrain_pos =
        MCMC.make_star_loglike(imgs, entry.pos;
                               patches=patches[1,:],
                               background_images=background_images,
                               pos_delta=pos_delta)

    function star_logprior(th)
        lnfluxes, pos = th[1:5], th[6:end]
        return MCMC.logflux_logprior(lnfluxes; is_star=true) + pos_logprior(pos)
    end

    function star_logpost(th)
        return star_loglike(th) + star_logprior(th)
    end

    # log likelihood at data generating parameters, and random prior
    th_cat = [log.(entry.star_fluxes)..., entry.pos...]
    #th_cat = MCMC.parameters_from_catalog(entry, unconstrain_pos; is_star=true)
    #th_rand = [th_cat[1:5]..., [.001, .001]...]
    th_rand = th_cat + .0001*randn(length(th_cat))
    Log.info("loglike at true initial position: ", star_loglike(th_cat))
    Log.info("loglike at random prior position: ", star_loglike(th_rand))
    Log.info("logprior at true position:        ", star_logprior(th_cat))

    # draw MCMC samples
    star_chain, star_lls = MCMC.slicesample_chain(star_logpost, th_cat,
        num_samples; print_skip=print_skip, verbose=false)
    num_warmup = Int(round(.25 * num_samples))
    star_chain = star_chain[num_warmup:end, :]
    star_lls   = star_lls[num_warmup:end]

    ####################
    # run galaxy MCMC  #
    ####################
    gal_loglike, constrain_pos, unconstrain_pos =
        MCMC.make_gal_loglike(imgs, entry.pos;
                              patches=patches[1,:],
                              background_images=background_images,
                              pos_delta=pos_delta)
    gal_logprior = MCMC.make_gal_logprior()

    # test at catalog initialized position vs shifted --- eye test
    th_cat = MCMC.parameters_from_catalog(entry; is_star=false)
    th_rand = th_cat + .00001*randn(length(th_cat))
    Log.info("gal ll at Catalog vs Shifted : ",
        gal_loglike(th_cat), " vs. ", gal_loglike(th_rand))
    Log.info("gal lnprior at true initial pos:  ", gal_logprior(th_cat))
    th_cat[1:5] = th_rand[1:5]

    function gal_logpost(th)
        # catch illegal values in the prior
        llprior = gal_logprior(th)
        if llprior < -1e100
            return llprior
        end
        return gal_loglike(th; print_params=false) + llprior
    end

    # draw MCMC samples
    gal_chain, gal_lls = MCMC.slicesample_chain(
        gal_logpost, th_cat, num_samples; print_skip=print_skip, verbose=false)
    gal_chain = gal_chain[num_warmup:end, :]
    gal_lls   = gal_lls[num_warmup:end]

    #############################################
    # Compute prob star vs gal based on samples #
    #############################################
    type_chain  = zeros(length(gal_lls))
    lnprob_a    = log(.5)
    lnprob_nota = log(1.-.5)
    for n in 1:length(gal_lls)
        # a == 1 full joint prob ln p(pixels | theta^star, a=star) p(theta^star) p(theta^gal) p(a)
        lnjoint_star = star_lls[n] + gal_logprior(gal_chain[n, :]) + lnprob_a

        # a == 0 (galaxy) full joint
        lnjoint_gal = gal_lls[n] + star_logprior(star_chain[n, :]) + lnprob_nota

        # compute log prob of type = star (1)
        lnsum = Model.logsumexp([lnjoint_star, lnjoint_gal])
        type_chain[n] = lnjoint_star - lnsum
    end
    ave_pstar = Model.logsumexp(type_chain) - log(length(gal_lls))
    Log.info("  source p-star = ", exp(ave_pstar))

    ####################################################################
    # convert positions to RA/Dec and organize chains into dataframes  #
    ####################################################################
    star_chain = MCMC.samples_to_dataframe(star_chain; is_star=true)
    gal_chain = MCMC.samples_to_dataframe(gal_chain; is_star=false)

    # store objid (for concatenation)
    mcmc_results = Dict("star_samples" => star_chain,
                        "star_lls"     => star_lls,
                        "gal_samples"  => gal_chain,
                        "gal_lls"      => gal_lls,
                        "type_samples" => type_chain)
    return mcmc_results
end

"""
Turn MCMC results into a single dataframe row that summarizes the
posterior distribution
"""
function summarize_samples(results_dict)
    stardf   = results_dict["star_samples"]
    galdf    = results_dict["gal_samples"]
    star_row = samples_to_dataframe_row(stardf; is_star=true)
    gal_row  = samples_to_dataframe_row(galdf; is_star=false)
    return star_row, gal_row
end

"""
Consolidate samples into results dataframes (single entry per source)
"""
function consolidate_samples(reslist; objids=nothing)
    # give each source a unique label
    if objids==nothing
        objids = ["samp_$(i)" for i in 1:length(reslist)]
    end

    # loop through each result sample list and summarize
    joint_summary, star_summary, gal_summary = [], [], []
    for i in 1:length(reslist)
        r = reslist[i]

        # compute star and gal summaries
        star_row, gal_row = summarize_samples(r)
        star_row[:objid], gal_row[:objid] = objids[i], objids[i]
        push!(star_summary, star_row)
        push!(gal_summary, gal_row)

        # compute star average
        pstar = exp(r["ave_pstar"])
        srow  = (pstar > .5) ? star_row : gal_row
        srow[:, :is_star] = [pstar]
        push!(joint_summary, srow)
    end
    joint_summary = vcat(joint_summary...)
    star_summary  = vcat(star_summary...)
    gal_summary   = vcat(gal_summary...)
    return star_summary, gal_summary, joint_summary
end


text:
struct TIVAlignmentResult{D,S,T,N,M,F<:AffineMap,G<:AbstractAffineMap,H<:AbstractAffineMap,V<:AbstractGMM{D,S},W<:AbstractGMM{D,T},X<:AbstractGMM{D,S},Y<:AbstractGMM{D,T}} <: AlignmentResults
    x::X
    y::Y
    upperbound::T
    lowerbound::T
    tform::F
    tform_params::NTuple{N,T}
    obj_calls::Int
    num_splits::Int
    num_blocks::Int
    stagnant_splits::Int
    rotation_result::GMAlignmentResult{D,S,T,M,G,V,W}
    translation_result::GMAlignmentResult{D,S,T,M,H,X,Y}
end

"""
    tgmm = tivgmm(gmm::IsotropicGMM, c=Inf)
    tgmm = tivgmm(mgmm::MultiGMM, c=Inf)

Returns a new `IsotropicGMM` or `MultiGMM` containing up to `c*length(gmm)` translation invariant vectors (TIVs) connecting Gaussian means in `gmm`.
TIVs are chosen to maximize length multiplied by the weights of the connected distributions.

See [Li et. al. (2019)](https://arxiv.org/abs/1812.11307) for a description of TIV construction.
"""
function tivgmm(gmm::AbstractIsotropicGMM, c=Inf)
    t = numbertype(gmm)
    npts, ndims = size(gmm)
    n = ceil(c*npts)
    if npts^2 < n
        n = npts^2
    end
    scores = fill(zero(t), npts, npts)
    for i=1:npts
        for j = i+1:npts
            scores[i,j] = scores[j,i] = norm(gmm.gaussians[i].μ-gmm.gaussians[j].μ) * √(gmm.gaussians[i].ϕ * gmm.gaussians[j].ϕ)
        end
    end

    tivgaussians = IsotropicGaussian{ndims,t}[]
    order = sortperm(vec(scores), rev=true)
    for idx in order[1:Int(n)]
        i = Int(floor((idx-1)/npts)+1)
        j = mod(idx-1, npts)+1
        x, y = gmm.gaussians[i], gmm.gaussians[j]
        push!(tivgaussians, IsotropicGaussian(x.μ-y.μ, √(x.σ*y.σ), √(x.ϕ*y.ϕ), x.dirs))
    end
    return IsotropicGMM(tivgaussians)
end

function tivgmm(mgmm::AbstractIsotropicMultiGMM, c=Inf)
    gmms = Dict{Symbol, IsotropicGMM{dims(mgmm),numbertype(mgmm)}}()
    for key in keys(mgmm.gmms)
        push!(gmms, Pair(key, tivgmm(mgmm.gmms[key],c)))
    end
    return IsotropicMultiGMM(gmms)
end

# fit a plane to a set of points, returning the normal vector
function planefit(pts)
    decomp = GenericLinearAlgebra.svd(pts .- sum(pts, dims=2))
    dist, nvecidx = findmin(decomp.S)
    return decomp.U[:, nvecidx], dist
end

function planefit(gmm::AbstractIsotropicGMM, R)
    ptsmat = fill(zero(numbertype(gmm)), 3, length(gmm))
    for (i,gauss) in enumerate(gmm.gaussians)
        ptsmat[:,i] = gauss.μ
    end
    return planefit(R * ptsmat)
end

function planefit(mgmm::AbstractIsotropicMultiGMM, R)
    len = sum([length(gmm) for gmm in values(mgmm.gmms)])
    ptsmat = fill(zero(numbertype(mgmm)), 3, len)
    idx = 1
    for gmm in values(mgmm.gmms)
        for gauss in gmm.gaussians
            ptsmat[:,idx] = gauss.μ
            idx += 1
        end
    end
    return planefit(R * ptsmat)
end

"""
    result = tiv_gogma_align(gmmx, gmmy, cx, cy; kwargs...)

Finds the globally optimal translation for alignment between two isotropic Gaussian mixtures, `gmmx`
and `gmmy`, using a modified GOGMA algorithm that performs rotational and translational optimization separately
by making use of translation invariant vectors (TIVs).

`cx` and `cy` are the ratios between number of TIVs used to represent a GMM and the number of Gaussians it contains,
for `gmmx` and `gmmy` respectively, during rotational alignment.

For details about keyword arguments, see `gogma_align()`.
"""
function tiv_gogma_align(gmmx::AbstractGMM, gmmy::AbstractGMM, cx=Inf, cy=Inf; kwargs...)
    t = promote_type(numbertype(gmmx),numbertype(gmmy))
    p = t(π)
    z = zero(t)

    # align TIVs
    tivgmmx, tivgmmy = tivgmm(gmmx, cx), tivgmm(gmmy, cy)
    rot_res = rot_gogma_align(tivgmmx, tivgmmy; kwargs...)
    rotblock = Block(((-p,p), (-p,p), (-p,p)), rot_res.tform_params, z, z)
    rotscore, rotpos = local_align(tivgmmx, tivgmmy, rotblock, objfun=rot_alignment_objective)

    # spin the moving tivgmm around to check for a better rotation (helps when the Gaussians are largely coplanar)
    R = AngleAxis(rot_res.tform_params...)
    spinvec, dist = planefit(tivgmmx, R)
    spinblock = Block(((-p,p), (-p,p), (-p,p)), rot_to_axis(AngleAxis(π*spinvec...) * R), z, z)
    spinscore, spinrotpos = local_align(tivgmmx, tivgmmy, spinblock, objfun=rot_alignment_objective)
    if spinscore < rotscore
        rotpos = spinrotpos
    end

    # perform translation alignment
    trl_res = trl_gogma_align(gmmx, gmmy; rot=rotpos, kwargs...)

    # perform local alignment in the full transformation space
    pos = NTuple{6, t}([rotpos..., trl_res.tform_params...])
    trlim = translation_limit(gmmx, gmmy)
    localblock = Block(((-p,p), (-p,p), (-p,p), (-trlim,trlim), (-trlim,trlim), (-trlim,trlim)), pos, trl_res.upperbound, typemin(t))
    min, bestpos = local_align(gmmx, gmmy, localblock)

    return TIVAlignmentResult(gmmx, gmmy, min, trl_res.lowerbound, AffineMap(bestpos...), bestpos,
                             rot_res.obj_calls+trl_res.obj_calls, rot_res.num_splits+trl_res.num_splits,
                             rot_res.num_blocks+trl_res.num_blocks, rot_res.stagnant_splits+trl_res.stagnant_splits,
                             rot_res, trl_res)
end


text:
#
# LocalFilters.jl --
#
# Local filters for Julia arrays.
#
#------------------------------------------------------------------------------
#
# This file is part of the `LocalFilters.jl` package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2018, Éric Thiébaut.
#

isdefined(Base, :__precompile__) && __precompile__(true)

module LocalFilters

export
    Neighborhood,
    RectangularBox,
    bilateralfilter!,
    bilateralfilter,
    bottom_hat,
    closing!,
    closing,
    convolve!,
    convolve,
    dilate!,
    dilate,
    erode!,
    erode,
    localextrema!,
    localextrema,
    localfilter!,
    localfilter,
    localmean!,
    localmean,
    opening!,
    opening,
    strel,
    top_hat

import Base: eltype, ndims, size, length, first, last, tail,
    getindex, setindex!, convert, reverse

# Deal with compatibility issues.
using Compat
const USE_CARTESIAN_RANGE = !isdefined(Base, :CartesianIndices)
@static if USE_CARTESIAN_RANGE
    import Base: CartesianRange
    import Compat: CartesianIndices
else
    import Base: CartesianIndices
end
@static if isdefined(Base, :axes)
    import Base: axes
else
    import Base: indices
    const axes = indices
end
@static if VERSION < v"0.7.0-alpha"
    Base.Tuple(index::CartesianIndex) = index.I
end

function localfilter end
function localfilter! end

include("types.jl")
include("basics.jl")
include("filters.jl")
include("morphology.jl")
include("separable.jl")
include("bilateral.jl")
import .BilateralFilter: bilateralfilter!, bilateralfilter

end