# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using CSV
using DataFrames
using JSON
using JuMP
using OrderedCollections
using Random
using TimerOutputs
using Printf
using Gurobi
   
include("jumpext.jl")
include("dist.jl")
 
dict = OrderedDict
Time = Int
Random.seed!(42)

# Structs
# =========================================================================

Base.@kwdef mutable struct Component
    name::String
end

Base.@kwdef mutable struct Product
    name::String
    comp::Vector{Component}
end

Base.@kwdef mutable struct Center
    name::String
    latitude::Float64
    longitude::Float64
    prod_out::Vector{Product}
end

Base.@kwdef mutable struct Plant
    name::String
    #county::String
    latitude::Float64
    longitude::Float64
    prod_out::Vector{Product}
    prod_in::Product
end

Base.@kwdef mutable struct Emission
    name::String
end

Base.show(io::IO, p::Union{Component,Product,Center,Plant,Emission}) = print(io, p.name)

Base.@kwdef mutable struct Instance
    T::UnitRange{Int}
    centers::Vector{Center}
    plants::Vector{Plant}
    products::Vector{Product}
    emissions::Vector{Emission}
    alpha_mix::dict{Tuple{Plant,Product,Component,Component},Float64}
    alpha_plant_emission::dict{Tuple{Plant,Emission,Time},Float64}
    alpha_tr_emission::dict{Tuple{Product,Emission,Time},Float64}
    c_acq::dict{Tuple{Center,Product,Time},Float64}
    #c_interv::dict{Tuple{Center,Product,Time},Float64}
    c_center_disp::dict{Tuple{Center,Product,Time},Float64}
    c_emission::dict{Tuple{Emission,Time},Float64}
    c_fix::dict{Tuple{Plant,Time},Float64}
    c_open::dict{Tuple{Plant,Time},Float64}
    c_plant_disp::dict{Tuple{Plant,Product,Time},Float64}
    c_store::dict{Tuple{Center,Product,Time},Float64}
    c_tr::dict{Tuple{Product,Time},Float64}
    c_var::dict{Tuple{Plant,Time},Float64}
    m_cap::dict{Plant,Float64}
    m_center_disp::dict{Tuple{Product,Time},Float64}
    m_dist::dict{Tuple{Union{Center,Plant},Plant},Float64}
    m_emission::dict{Tuple{Emission,Time},Float64}
    m_init::dict{Tuple{Center,Product,Component,Time},Float64}
    m_plant_disp::dict{Tuple{Plant,Product,Time},Float64}
    m_store::dict{Tuple{Center,Product,Time},Float64}
    selected_edges::Union{Nothing,Set} = nothing
    max_cap::dict{Plant,Float64}
    max_OpC::dict{Plant,Float64}
    max_OpenC::dict{Plant,Float64}
    min_cap::dict{Plant,Float64}
    min_OpC::dict{Plant,Float64}
    min_OpenC::dict{Plant,Float64}
    unit_open::dict{Plant,Float64}
    unit_OpC::dict{Plant,Float64}
end


# Generate
# ==============================================================================

function generate_data()
    # Time window
    T = 1:2

    # Cities
    cities_a = dict(
        "Chicago" => [41.881832, -87.623177],
        "New York City" => [40.712776, -74.005974],
        "Los Angeles" => [34.052235, -118.243683],
        "Houston" => [29.760427, -95.369804],
        "Phoenix" => [33.448376, -112.074036],
        "Philadelphia" => [39.952583, -75.165222],
        "San Antonio" => [29.424122, -98.493629],
        "San Diego" => [32.715736, -117.161087],
        "Dallas" => [32.776664, -96.796988],
        "San Jose" => [37.338208, -121.886329],
    )

    cities_b = dict(
        "Chicago" => [41.881832, -87.623177],
        "Phoenix" => [33.448376, -112.074036],
        "Dallas" => [32.776664, -96.796988],
    )

    # Components
    film = Component("Film")
    paper = Component("Paper")
    cardboard = Component("Cardboard")

    # Products
    waste = Product(name="Waste", comp=[film, paper, cardboard])
    film_bale = Product(name="Film bale", comp=[film, paper, cardboard])
    cardboard_bale = Product(name="Cardboard bale", comp=[paper, cardboard])
    cardboard_sheets = Product(name="Cardboard sheets", comp=[cardboard])
    products = [waste, film_bale, cardboard_bale, cardboard_sheets]

    # Centers
    centers = [
        Center(
            name="Collection ($city_name)",
            latitude=city_lat,
            longitude=city_lon,
            prod_out=[waste],
        ) for (city_name, (city_lat, city_lon)) in cities_a
    ]

    # Plants
    plants_a = [
        Plant(
            name="MRF ($city_name)",
            latitude=city_lat,
            longitude=city_lon,
            prod_in=waste,
            prod_out=[film_bale, cardboard_bale],
        ) for (city_name, (city_lat, city_lon)) in cities_b
    ]
    plants_b = [
        Plant(
            name="Paper Mill ($city_name)",
            latitude=city_lat,
            longitude=city_lon,
            prod_in=cardboard_bale,
            prod_out=[cardboard_sheets],
        ) for (city_name, (city_lat, city_lon)) in cities_b
    ]
    plants = [plants_a; plants_b]

    # Emissions
    emissions = [Emission("CO2")]

    alpha_mix = dict(
        (p, r, cin, cout) => 0.0 for p in plants for cin in p.prod_in.comp for
        r in p.prod_out for cout in r.comp
    )
    for p in plants_a
        alpha_mix[p, film_bale, film, film] = 0.98
        alpha_mix[p, film_bale, paper, paper] = 0.02
        alpha_mix[p, film_bale, cardboard, cardboard] = 0.02
        alpha_mix[p, cardboard_bale, paper, paper] = 0.02
        alpha_mix[p, cardboard_bale, cardboard, cardboard] = 0.75
    end
    for p in plants_b
        alpha_mix[p, cardboard_sheets, cardboard, cardboard] = 0.95
    end
    alpha_plant_emission = dict((p, s, t) => 0.01 for p in plants, s in emissions, t in T)
    alpha_tr_emission = dict((r, s, t) => 0.01 for r in products, s in emissions, t in T)
    c_acq = dict((q, r, t) => 1.0 for q in centers for r in q.prod_out, t in T)
    c_center_disp = dict((q, r, t) => 0 for q in centers for r in q.prod_out for t in T)
    c_emission = dict((s, t) => 0.01 for s in emissions, t in T)
    c_fix = dict((p, t) => 1_000.0 for p in plants, t in T)
    c_open = dict((p, t) => 10_000.0 for p in plants, t in T)
    c_plant_disp = dict(
        (p, r, t) => (r == cardboard_sheets ? -100.0 : -10.0) for p in plants for
        r in p.prod_out, t in T
    )
    c_store = dict((q, r, t) => 1.0 for q in centers for r in q.prod_out, t in T)
    c_tr = dict((r, t) => 0.05 for r in products, t in T)
    c_var = dict((p, t) => 1.0 for p in plants, t in T)
    m_cap = dict(p => 50000.0 for p in plants)
    m_center_disp = dict((r, t) => 0.0 for r in products, t in T)
    metric = KnnDrivingDistance()
    m_dist = dict(
        (p, q) => _calculate_distance(
            p.latitude,
            p.longitude,
            q.latitude,
            q.longitude,
            metric,
        ) for p in [plants; centers], q in plants
    )
    m_emission = dict((s, t) => 1_000_000 for s in emissions, t in T)
    m_init = dict()
    for q in centers, r in q.prod_out
        ratio = dict(c => rand(1:10) for c in r.comp)
        total = dict(t => rand(1:1000) for t in T)
        for c in r.comp, t in T
            m_init[q, r, c, t] = ratio[c] * total[t]
        end
    end
    m_plant_disp =
        dict((p, r, t) => 1_000_000 for p in plants for r in p.prod_out for t in T)
    m_store = dict((q, r, t) => 1_000 for q in centers for r in q.prod_out, t in T)

    return Instance(;
        T,
        centers,
        plants,
        products,
        emissions,
        alpha_mix,
        alpha_plant_emission,
        alpha_tr_emission,
        c_acq,
        c_center_disp,
        c_emission,
        c_fix,
        c_open,
        c_plant_disp,
        c_store,
        c_tr,
        c_var,
        m_cap,
        m_center_disp,
        m_dist,
        m_emission,
        m_init,
        m_plant_disp,
        m_store,
    )
end

# Write
# ==============================================================================

function write_json(data, filename)
    json = dict()
    json["parameters"] = dict("time horizon (years)" => data.T.stop)
    json["products"] = dict(
        r.name => dict(
            "components" => [c.name for c in r.comp],
            "disposal limit (tonne)" => [data.m_center_disp[r, t] for t in data.T],
            "transportation cost (\$/km/tonne)" => [data.c_tr[r, t] for t in data.T],
            "transportation emissions (tonne/km/tonne)" => dict(
                s => [data.alpha_tr_emission[r, s, t] for t in data.T] for
                s in data.emissions
            ),
        ) for r in data.products
    )
    json["centers"] = dict(
        q.name => dict(
            "latitude" => q.latitude,
            "longitude" => q.longitude,
            "output" => dict(
                r.name => dict(
                    "initial amount (tonne)" => [
                        data.m_init[q, r, c, t] for c in r.comp, t in data.T
                    ],
                    "disposal cost (\$/tonne)" =>
                        [data.c_center_disp[q, r, t] for t in data.T],
                    "storage cost (\$/tonne)" =>
                        [data.c_store[q, r, t] for t in data.T],
                    "storage limit (tonne)" => [data.m_store[q, r, t] for t in data.T],
                    "acquisition cost (\$/tonne)" =>
                        [data.c_acq[q, r, t] for t in data.T],
                ) for r in q.prod_out
            ),
        ) for q in data.centers
    )
    json["plants"] = dict(
        p.name => dict(
            "latitude" => p.latitude,
            "longitude" => p.longitude,
            "input" => p.prod_in.name,
            "output" => dict(
                r.name => dict(
                    "output matrix" => [
                        data.alpha_mix[p, r, c_in, c_out] for
                        c_in in p.prod_in.comp, c_out in r.comp
                    ],
                    "disposal limit (tonne)" =>
                        [data.m_plant_disp[p, r, t] for t in data.T],
                    "disposal cost (\$/tonne)" =>
                        [data.c_plant_disp[p, r, t] for t in data.T],
                ) for r in p.prod_out
            ),
            
            "Min_capacity" => data.min_cap[p],
            "Min_FixOperatingCost" => data.min_OpC[p],
            "Min_OpeningCost" => data.min_OpenC[p],
            "Max_capacity" => data.max_cap[p],
            "Max_FixOperatingCost" => data.max_OpC[p],
            "Max_OpeningCost" => data.max_OpenC[p],
            "fixed operating cost (\$)" => [data.c_fix[p, t] for t in data.T],
            "variable operating cost (\$/tonne)" => [data.c_var[p, t] for t in data.T],
            "opening cost (\$)" => [data.c_open[p, t] for t in data.T],
            "capacity (tonne)" => data.m_cap[p],
            "emissions (tonne/tonne)" => dict(
                s => [data.alpha_plant_emission[p, s, t] for t in data.T] for
                s in data.emissions
            ),
        ) for p in data.plants
    )
    json["emissions"] = dict(
        s.name => dict(
            "penalty (\$/tonne)" => [data.c_emission[s, t] for t in data.T],
            "limit (tonne)" => [data.m_emission[s, t] for t in data.T],
        ) for s in data.emissions
    )
    open(filename, "w") do io
        JSON.print(io, json, 2)
    end
end

# Read
# ==============================================================================
function read_json(filename; max_centers=Inf, max_plants=Inf)::Instance
    json = JSON.parsefile(filename)
    T = 1:json["parameters"]["time horizon (years)"]
    centers = []
    components_by_name = dict()
    emissions = []
    emissions_by_name = dict()
    plants = []
    products = []
    products_by_name = dict()
    alpha_mix = dict()
    alpha_plant_emission = dict()
    alpha_tr_emission = dict()
    c_acq = dict()
    #c_interv = dict()
    c_center_disp = dict()
    c_emission = dict()
    c_fix = dict()
    c_open = dict()
    c_plant_disp = dict()
    c_store = dict()
    c_tr = dict()
    c_var = dict()
    m_cap = dict()
    m_center_disp = dict()
    m_emission = dict()
    m_init = dict()
    m_plant_disp = dict()
    m_store = dict()
    max_cap = dict()
    max_OpC = dict()
    max_OpenC = dict()
    min_cap = dict()
    min_OpC = dict()
    min_OpenC = dict()
    unit_open = dict()
    unit_OpC = dict()

    @timeit "Read: Emissions" begin
        for (emission_name, emission_data) in json["emissions"]
            s = Emission(emission_name)
            emissions_by_name[emission_name] = s
            push!(emissions, s)
            for t in T
                c_emission[s, t] = emission_data["penalty (\$/tonne)"][t]
                m_emission[s, t] = emission_data["limit (tonne)"][t]
            end
        end
    end

    @timeit "Read: Products" begin
        for (name, prod_data) in json["products"]
            print(name)
            comp = []
            for (comp_name) in prod_data["components"]
                if comp_name ∉ keys(components_by_name)
                    print(comp_name)
                    components_by_name[comp_name] = Component(comp_name)
                end
                push!(comp, components_by_name[comp_name])
            end
            r = Product(name, comp)
            products_by_name[name] = r
            push!(products, r)
            for t in T
                m_center_disp[r, t] = prod_data["disposal limit (tonne)"][t]
                c_tr[r, t] = prod_data["transportation cost (\$/km/tonne)"][t]
            end
            for (s_name, s_data) in prod_data["transportation emissions (tonne/km/tonne)"]
                s = emissions_by_name[s_name]
                for t in T
                    alpha_tr_emission[r, s, t] = s_data[t]
                end
            end
        end
    end

    @timeit "Read: Centers" begin
        for (name, center_data) in json["centers"]
            if length(centers) >= max_centers
                @warn "Maximum number of centers reached. Skipping remaining ones."
                break
            end
            latitude = center_data["latitude"]
            longitude = center_data["longitude"]
            prod_out = [products_by_name[r] for r in keys(center_data["output"])]
            q = Center(; name, latitude, longitude, prod_out)
            push!(centers, q)
            for r in prod_out, t in T
                c_acq[q, r, t] =
                    center_data["output"][r.name]["acquisition cost (\$/tonne)"][t]
                #c_interv[q, r, t] =
                 #   center_data["output"][r.name]["intervention costs ($)"][t]
                c_center_disp[q, r, t] =
                    center_data["output"][r.name]["disposal cost (\$/tonne)"][t]
                c_store[q, r, t] =
                    center_data["output"][r.name]["storage cost (\$/tonne)"][t]
                m_store[q, r, t] = center_data["output"][r.name]["storage limit (tonne)"][t]
                for (c_idx, c) in enumerate(r.comp)
                    m_init[q, r, c, t] =
                        center_data["output"][r.name]["initial amount (tonne)"][t][c_idx]
                end
            end
        end
    end

    @timeit "Read: Plants" begin
        for (plant_name, plant_data) in json["plants"]
            if length(plants) >= max_plants
                @warn "Maximum number of plants reached. Skipping remaining ones."
                break
            end
            #county = plant_data["county"]
            latitude = plant_data["latitude"]
            longitude = plant_data["longitude"]
            prod_in = products_by_name[plant_data["input"]]
            prod_out = [products_by_name[r] for r in keys(plant_data["output"])]
            p = Plant(plant_name, latitude, longitude, prod_out, prod_in)
            push!(plants, p)
            m_cap[p] = plant_data["capacity (tonne)"]
            min_cap[p] = plant_data["Min_capacity"]
            min_OpC[p] = plant_data["Min_FixOperatingCost"]
            min_OpenC[p] = plant_data["Min_OpeningCost"]
            max_cap[p] = plant_data["Max_capacity"]
            max_OpC[p] = plant_data["Max_FixOperatingCost"]
            max_OpenC[p] = plant_data["Max_OpeningCost"]

            #Compute unit cost for expansion calculations----------------------
            #unit opening cost
            if max_cap[p] > min_cap[p]
                unit_open[p] = (max_OpenC[p] - min_OpenC[p])/(max_cap[p] - min_cap[p])
            else
                unit_open[p] = (max_OpenC[p]/max_cap[p])
            end
            #unit operating cost
            if max_cap[p] > min_cap[p]
                unit_OpC[p] = (max_OpC[p] - min_OpC[p])/(max_cap[p] - min_cap[p])
            else
                unit_OpC[p] = (max_OpC[p]/max_cap[p])
            end
            #------------------------------------------------------------------
            
            for t in T
                c_fix[p, t] = plant_data["fixed operating cost (\$)"][t]
                c_var[p, t] = plant_data["variable operating cost (\$/tonne)"][t]
                c_open[p, t] = plant_data["opening cost (\$)"][t]
            end
            for r in prod_out,
                (cin_idx, c_in) in enumerate(prod_in.comp),
                (cout_idx, c_out) in enumerate(r.comp)

                alpha_mix[p, r, c_in, c_out] =
                    plant_data["output"][r.name]["output matrix"][cout_idx][cin_idx]
            end
            for r in prod_out, t in T
                c_plant_disp[p, r, t] =
                    plant_data["output"][r.name]["disposal cost (\$/tonne)"][t]
                m_plant_disp[p, r, t] =
                    plant_data["output"][r.name]["disposal limit (tonne)"][t]
            end
            for (s_name, s_data) in plant_data["emissions (tonne/tonne)"]
                s = emissions_by_name[s_name]
                for t in T
                    alpha_plant_emission[p, s, t] = s_data[t]
                end
            end
        end
    end

    @timeit "Calculate distances" begin
        metric = KnnDrivingDistance()
        m_dist = dict(
            (p, q) => _calculate_distance(
                p.latitude,
                p.longitude,
                q.latitude,
                q.longitude,
                metric,
            ) for p in [plants; centers], q in plants
        )
    end

    return Instance(;
        T,
        centers,
        plants,
        products,
        emissions,
        alpha_mix,
        alpha_plant_emission,
        alpha_tr_emission,
        c_acq,
        #c_interv,
        c_center_disp,
        c_emission,
        c_fix,
        c_open,
        c_plant_disp,
        c_store,
        c_tr,
        c_var,
        m_cap,
        m_center_disp,
        m_dist,
        m_emission,
        m_init,
        m_plant_disp,
        m_store,
        max_cap,
        max_OpC,
        max_OpenC,
        min_cap,
        min_OpC,
        min_OpenC,
        unit_open,
        unit_OpC
    )
end

# Multi-period Heuristic
# ==============================================================================

function compress(original::Instance)::Instance
    T = original.T
    alpha_plant_emission = Dict(
        (p, s, 1) => mean([original.alpha_plant_emission[p, s, t] for t in T]) for
        p in original.plants, s in original.emissions
    )
    alpha_tr_emission = Dict(
        (r, s, 1) => mean([original.alpha_tr_emission[r, s, t] for t in T]) for
        r in original.products, s in original.emissions
    )
    c_acq = Dict(
        (q, r, 1) => mean([original.c_acq[q, r, t] for t in T]) for
        q in original.centers for r in q.prod_out
    )
    #c_interv = Dict(
    #    (q, r, 1) => mean([original.c_interv[q, r, t] for t in T]) for
    #    q in original.centers for r in q.prod_out
    #)
    c_center_disp = Dict(
        (q, r, 1) => mean([original.c_center_disp[q, r, t] for t in T]) for
        q in original.centers for r in q.prod_out
    )
    c_emission = Dict(
        (s, 1) => mean([original.c_emission[s, t] for t in T]) for s in original.emissions
    )
    c_fix = Dict((p, 1) => sum([original.c_fix[p, t] for t in T]) for p in original.plants)
    c_open =
        Dict((p, 1) => mean([original.c_open[p, t] for t in T]) for p in original.plants)
    c_plant_disp = Dict(
        (p, r, 1) => mean([original.c_plant_disp[p, r, t] for t in T]) for
        p in original.plants for r in p.prod_out
    )
    c_store = Dict(
        (q, r, 1) => mean([original.c_store[q, r, t] for t in T]) for
        q in original.centers for r in q.prod_out
    )
    c_tr = Dict((r, 1) => mean([original.c_tr[r, t] for t in T]) for r in original.products)
    c_var = Dict((p, 1) => mean([original.c_var[p, t] for t in T]) for p in original.plants)
    m_cap = Dict(p => original.m_cap[p] * length(T) for p in original.plants)
    max_cap = Dict(p => original.max_cap[p] * length(T) for p in original.plants)
    min_cap = Dict(p => original.min_cap[p] * length(T) for p in original.plants)
   
    max_OpC = Dict(p => original.max_OpC[p] * length(T) for p in original.plants)
    min_OpC = Dict(p => original.min_OpC[p] * length(T) for p in original.plants)
    
    max_OpenC = Dict(p => original.max_OpenC[p] * length(T) for p in original.plants)
    min_OpenC = Dict(p => original.min_OpenC[p] * length(T) for p in original.plants)
    
    unit_open = Dict(p => original.unit_open[p] * length(T) for p in original.plants)
    #unit_open = Dict(p => ((max_OpenC[p] - min_OpenC[p])/(max_cap[p] - min_cap[p])) for p in original.plants)
    #unit_open = Dict(p => original.unit_open[p] for p in original.plants)
    
    unit_OpC = Dict(p => original.unit_OpC[p] * length(T) for p in original.plants)
    #unit_OpC =  Dict(p => ((max_OpC[p] - min_OpC[p])/(max_cap[p] - min_cap[p])) for p in original.plants)
    #unit_OpC = Dict(p => original.unit_OpC[p] for p in original.plants)
    
    m_center_disp = Dict(
        (r, 1) => sum([original.m_center_disp[r, t] for t in T]) for r in original.products
    )
    m_emission = Dict(
        (s, 1) => sum([original.m_emission[s, t] for t in T]) for s in original.emissions
    )
    m_init = Dict(
        (q, r, c, 1) => length(T) * maximum([original.m_init[q, r, c, t] for t in T]) for
        q in original.centers for r in q.prod_out for c in r.comp
    ) 
    m_plant_disp = Dict(
        (p, r, 1) => length(T) * minimum([original.m_plant_disp[p, r, t] for t in T]) for
        p in original.plants for r in p.prod_out
    )
    m_store = Dict(
        (q, r, 1) => length(T) * minimum([original.m_store[q, r, t] for t in T]) for
        q in original.centers for r in q.prod_out
    )
    #=
    m_init = Dict(
        (q, r, c, 1) => sum([original.m_init[q, r, c, t] for t in T]) for
        q in original.centers for r in q.prod_out for c in r.comp
    )
    m_plant_disp = Dict(
        (p, r, 1) => sum([original.m_plant_disp[p, r, t] for t in T]) for
        p in original.plants for r in p.prod_out
    )
    m_store = Dict(
        (q, r, 1) => sum([original.m_store[q, r, t] for t in T]) for
        q in original.centers for r in q.prod_out
    )
    =#
    return Instance(;
        T=1:1,
        centers=original.centers,
        plants=original.plants,
        products=original.products,
        emissions=original.emissions,
        alpha_mix=original.alpha_mix,
        alpha_plant_emission,
        alpha_tr_emission,
        c_acq,
        #c_interv,
        c_center_disp,
        c_emission,
        c_fix,
        c_open,
        c_plant_disp,
        c_store,
        c_tr,
        c_var,
        m_cap,
        m_center_disp,
        m_dist=original.m_dist,
        m_emission,
        m_init,
        m_plant_disp,
        m_store,
        max_cap,
        max_OpC,
        max_OpenC,
        min_cap,
        min_OpC,
        min_OpenC,
        unit_open,
        unit_OpC
    )
end

function benchmark_compress(filename, optimizer; max_centers=[Inf], max_plants=[Inf])
    stats = []
    for mc in max_centers, mp in max_plants
        # Solve original
        orig = read_json(filename; max_centers=mc, max_plants=mp)
        reset_timer!()
        _, stats_orig = solve(orig; optimizer)
        stats_orig["Filename"] = filename
        stats_orig["Method"] = "Original"
        push!(stats, stats_orig)

        # Solve compressed
        compressed = compress(orig)
        reset_timer!()
        model_comp, stats_comp = solve(compressed; optimizer)
        stats_comp["Filename"] = filename
        stats_comp["Method"] = "Compressed"
        push!(stats, stats_comp)

    end

    return DataFrame(stats)
end


# Run
# ==============================================================================

function generate_json()
    @info "Generating data"
    data = generate_data()

    @info "Writing JSON file"
    write_json(data, "output-3/case.json")
end

function solve(filename, optimizer; max_centers=Inf, max_plants=Inf, heuristic=false)
    reset_timer!()
    @timeit "Read JSON" begin
        data = read_json(filename; max_centers, max_plants)
        write_json(data, "output-3/BaselinePlastics_capExTestCase10C7P.json")
    end
    return solve(data; optimizer=optimizer, output_dir=dirname(filename), heuristic)
    print_timer()
end

function solve(data::Instance; optimizer, output_dir=nothing, heuristic=false)

    if heuristic
        @info "Solving compressed instance..."
        comp_data = compress(data)
        #print("data compressed")
        comp_model, comp_stats = solve(comp_data; optimizer, heuristic=false)
        #print("compressed model solved")
        
        # Filter plants
        selected_plants = Plant[]
        for p in comp_data.plants
            if value(comp_model[:x_open][p, 1]) > 0.5
                push!(selected_plants, p)
            end
        end
        @info "Selected $(length(selected_plants)) out of $(length(comp_data.plants)) plants"

          
        # Filter edges
        selected_edges = Set()
        for (src, dst, r) in comp_model[:E]
            if value(comp_model[:y_total][src, dst, r, 1]) > 0.5
                push!(selected_edges, (src, dst, r))
            end
        end
        @info "Selected $(length(selected_edges)) out of $(length(comp_model[:E])) transportation edges"
        

        data = Instance(;
            data.T,
            data.centers,
            plants=selected_plants,
            data.products,
            data.emissions,
            data.alpha_mix,
            data.alpha_plant_emission,
            data.alpha_tr_emission,
            data.c_acq,
            #data.c_interv,
            data.c_center_disp,
            data.c_emission,
            data.c_fix,
            data.c_open,
            data.c_plant_disp,
            data.c_store,
            data.c_tr,
            data.c_var,
            data.m_cap,
            data.m_center_disp,
            data.m_dist,
            data.m_emission,
            data.m_init,
            data.m_plant_disp,
            data.m_store,
            selected_edges,
            data.max_cap,
            data.max_OpC,
            data.max_OpenC,
            data.min_cap,
            data.min_OpC,
            data.min_OpenC,
            data.unit_open,
            data.unit_OpC
            
        )
    end


    T = data.T
    centers = data.centers
    plants = data.plants
    products = data.products
    emissions = data.emissions

    #model = Model(optimizer)
    model = Model(Gurobi.Optimizer)
    set_optimizer_attribute(model, "MIPGap", 0.005)
    init_time = time()

    # Graph
    # -------------------------------------------------------------------------
    @timeit "Build graph" begin
        E = []
        E_in = dict(src => [] for src in plants)
        E_out = dict(src => [] for src in plants ∪ centers)
        function push_edge!(src, dst, r)
            e = (src, dst, r)
            if data.selected_edges !== nothing && e ∉ data.selected_edges
                return
            end
            push!(E, e)
            push!(E_out[src], (dst, r))
            push!(E_in[dst], (src, r))
        end
        for r in products
            # Plant to plant
            for p1 in plants
                r ∈ p1.prod_out || continue
                for p2 in plants
                    p1 != p2 || continue
                    r == p2.prod_in || continue
                    push_edge!(p1, p2, r)
                end
            end
            # Center to plant
            for q in centers
                r ∈ q.prod_out || continue
                for p in plants
                    r == p.prod_in || continue
                    push_edge!(q, p, r)
                end
            end
        end
    end
    model[:E] = E
    model[:E_in] = E_in
    model[:E_out] = E_out

    @printf("Building optimization problem with:\n")
    @printf("    %8d plants\n", length(plants))
    @printf("    %8d centers\n", length(centers))
    @printf("    %8d products\n", length(products))
    @printf("    %8d time periods\n", length(T))
    @printf("    %8d transportation edges\n", length(E))

    available = Dict((r, t) => 0.0 for r in data.products, t in T)
    for q in centers, t in T
        for r in q.prod_out, s in r.comp
            available[r, t] += data.m_init[q, r, s, t]
        end
    end

    capacity = Dict(r => 0.0 for r in data.products)
    for p in plants
        capacity[p.prod_in] += data.max_cap[p]
    end

    for r in products, t in T
        if available[r, t] > capacity[r]
            @warn "Not enough capacity to process $(r.name) at time $t: $(available[r,t]) > $(capacity[r])"
        end
    end


    # Decision variables
    # -------------------------------------------------------------------------
    @timeit "Model: Add variables" begin
        @timeit "y" begin
            y = _init(model, :y)
            for (q, p, r) in E, c in r.comp, t in T
                y[q, p, r, c ,t] = @variable(model, lower_bound = 0)
            end
        end
        @timeit "y_total" begin
            y_total = _init(model, :y_total)
            for (q, p, r) in E, t in T
                y_total[q, p, r, t] = @variable(model, lower_bound = 0)
            end
        end
        @timeit "z_center_disp" begin
            z_center_disp = _init(model, :z_center_disp)
            for q in centers, r in q.prod_out, c in r.comp, t in T
                z_center_disp[q, r, c, t] = @variable(model, lower_bound = 0)
            end
        end
        @timeit "z_center_disp_total" begin
            z_center_disp_total = _init(model, :z_center_disp_total)
            for q in centers, r in q.prod_out, t in T
                z_center_disp_total[q, r, t] = @variable(model, lower_bound = 0)
            end
        end
        @timeit "z_store" begin
            z_store = _init(model, :z_store)
            for q in centers, r in q.prod_out, c in r.comp, t in T
                if t == T.stop
                    z_store[q, r, c, t] = 0.0
                else
                    z_store[q, r, c, t] = @variable(model, lower_bound = 0)
                end
            end
        end
        @timeit "z_store_total" begin
            z_store_total = _init(model, :z_store_total)
            for q in centers, r in q.prod_out
                z_store_total[q, r, 0] = 0.0
                for t in T
                    if t == T.stop
                        z_store_total[q, r, t] = 0.0
                    else
                        z_store_total[q, r, t] = @variable(model, lower_bound = 0)
                    end
                end
            end
        end
        #=
        @timeit "x_open" begin
            x_open = _init(model, :x_open)
            for p in plants
                x_open[p, 0] = 0
                for t in T
                    x_open[p, t] = @variable(model, binary = true)
                end
            end
        end
        =#
        
        @timeit "x_open" begin
            x_open = _init(model, :x_open)
            for p in plants
                if data.m_cap[p]>0
                    x_open[p,0] = 1
                else
                    x_open[p,0] = 0
                end 
                for t in T
                    x_open[p, t] = @variable(model, binary = true)
                end
            end
        end
        
        @timeit "x_send" begin
            x_send = _init(model, :x_send)
            for p in plants, (q, r) in E_out[p], t in T
                x_send[p, q, r, t] = @variable(model, binary = true)
            end
        end
        @timeit "x_disp" begin
            x_disp = _init(model, :x_disp)
            for p in plants, r in p.prod_out, t in T
                x_disp[p, r, t] = @variable(model, binary = true)
            end
        end
        @timeit "z_prod" begin
            z_prod = _init(model, :z_prod)
            for p in plants, r in p.prod_out, c in r.comp, t in T
                z_prod[p, r, c, t] = @variable(model, lower_bound = 0)
            end
        end
        @timeit "z_plant_disp" begin
            z_plant_disp = _init(model, :z_plant_disp)
            for p in plants, r in p.prod_out, c in r.comp, t in T
                z_plant_disp[p, r, c, t] = @variable(model, lower_bound = 0)
            end
        end
        @timeit "z_tr_emissions" begin
            z_tr_emissions = _init(model, :z_tr_emissions)
            for (q, p, r) in E, s in emissions, t in T
                z_tr_emissions[q, p, r, s, t] = @variable(model, lower_bound = 0)
            end
        end
        @timeit "z_plant_emissions" begin
            z_plant_emissions = _init(model, :z_plant_emissions)
            for p in plants, s in emissions, t in T
                z_plant_emissions[p, s, t] = @variable(model, lower_bound = 0)
            end
        end

        ########################## FOR CAPACITY EXPANSION ##############################################################
        #Binary variable - if a new plant is built
        @timeit "x_build" begin
            x_build = _init(model, :x_build)
            for p in plants, t in T
                x_build[p, t] = @variable(model, binary = true)
            end
        end
            
        #Plant capacity expansion Variable
        @timeit "z_capExp" begin
            z_capExp = _init(model, :z_capExp)
            for p in plants, t in T
                z_capExp[p, t] = @variable(model, lower_bound = 0, upper_bound = (data.max_cap[p] - data.m_cap[p]))      
            end
        end
        
        #Total plant capacity Variable
        @timeit "z_capTot" begin
            z_capTot = _init(model, :z_capTot)
            for p in plants
                z_capTot[p, 0] = data.m_cap[p]
                for t in T
                    z_capTot[p, t] = @variable(model, lower_bound = 0, upper_bound = data.max_cap[p])  
                end
            end
        end
        ##############################################################################################################
    end

    # Objective function
    # -------------------------------------------------------------------------
    @timeit "Model: Objective function" begin
        obj = AffExpr()

        # Center disposal
        @timeit "c_center_disp" begin
            for q in centers, r in q.prod_out, t in T
                add_to_expression!(
                    obj,
                    data.c_center_disp[q, r, t],
                    z_center_disp_total[q, r, t],
                )
            end
        end

        # Center acquisition
        @timeit "c_acq" begin
            for q in centers, r in q.prod_out, c in r.comp, t in T
                add_to_expression!(obj, data.m_init[q, r, c, t] * data.c_acq[q, r, t]) #+data.c_interv[q,r,t])
            end
        end

        # Center storage
        @timeit "c_store" begin
            for q in centers, r in q.prod_out, t in T
                add_to_expression!(obj, data.c_store[q, r, t], z_store_total[q, r, t])
            end
        end

        # Transportation, variable operating cost
        @timeit "c_tr, c_var" begin
            for (q, p, r) in E
                dist = data.m_dist[q, p]
                #y_qpr = y[q, p, r]
                for c in r.comp, t in T
                    add_to_expression!(obj, dist * data.c_tr[r, t], y[q, p, r,c, t])
                    add_to_expression!(obj, data.c_var[p, t], y[q, p, r, c, t])
                end
            end
        end

        # Transportation emissions
        @timeit "c_emission" begin
            for (q, p, r) in E, s in emissions, t in T
                add_to_expression!(
                    obj,
                    data.c_emission[s, t],
                    z_tr_emissions[q, p, r, s, t],
                )
            end
        end

        # Fixed cost
        @timeit "c_fix" begin
            for p in plants, t in T
                add_to_expression!(obj, data.c_fix[p, t], x_open[p, t])
            end
        end

        #=
        # Opening cost
        @timeit "c_open" begin
            for p in plants, t in T
                add_to_expression!(obj, data.c_open[p, t], x_open[p, t] - x_open[p, t-1])
            end
        end
        =#

        # Plant disposal
        @timeit "c_plant_disp" begin
            for p in plants, r in p.prod_out, c in r.comp, t in T
                add_to_expression!(
                    obj,
                    data.c_plant_disp[p, r, t],
                    z_plant_disp[p, r, c, t],
                )
            end
        end

        # Plant emissions
        @timeit "c_emission" begin
            for p in plants, s in emissions, t in T
                add_to_expression!(obj, data.c_emission[s, t], z_plant_emissions[p, s, t])
            end
        end

        ######################## FOR CAPACITY EXPANSION ####################################################
        
        #Plant opening cost
        @timeit "c_open" begin
            for p in plants, t in T
                add_to_expression!(
                    obj,
                    data.min_OpenC[p],
                    x_build[p, t]
                )
            end
        end
    
        #Plant Expansion cost
        @timeit "c_exp" begin
            for p in plants, t in T
                add_to_expression!(
                    obj,
                    data.unit_open[p],
                    (z_capExp[p, t])
                )
            end
        end
    
        #Plant additional fixed operating cost due to expansion
        @timeit "c_fix_exp" begin
            for p in plants, t in T
                add_to_expression!(
                    obj,
                    data.unit_OpC[p],
                    sum(z_capExp[p, i] for i = 1:t)
                )
            end
        end
        
        ####################################################################################################

        @objective(model, Min, obj)
    end

    # Constraints
    # -------------------------------------------------------------------------
    @timeit "Model: Constraints" begin
        @timeit "eq_balance" begin
            eq_balance = _init(model, :eq_balance)
            for q in centers, r in q.prod_out, t in T
                eq_balance[q, r, t] = @constraint(
                    model,
                    sum(y_total[q, p, r2, t] for (p, r2) in E_out[q] if r == r2) +
                    z_store_total[q, r, t] + sum(z_center_disp[q, r, c, t] for c in r.comp) ==
                    sum(data.m_init[q, r, c, t] for c in r.comp) + z_store_total[q, r, t-1]
                )
            end
        end
        ratio(q, r, c, t) =
            (data.m_init[q, r, c, t] / sum(data.m_init[q, r, d, t] for d in r.comp))

        @timeit "eq_y_total" begin
            eq_y_total = _init(model, :eq_y_total)
            for (q, p, r) in E, t in T
                eq_y_total[q, p, r, t] = @constraint(
                    model,
                    y_total[q, p, r, t] == sum(y[q, p, r, c, t] for c in r.comp)
                )
            end
        end
        @timeit "eq_split_y" begin
            eq_split_y = _init(model, :eq_split_y)
            for (q, p, r) in E, c in r.comp, t in T
                if q ∉ centers
                    continue
                end
                eq_split_y[q, p, r, c, t] = @constraint(
                    model,
                    y[q, p, r, c, t] == ratio(q, r, c, t) * y_total[q, p, r, t]
                )
            end
        end
        @timeit "eq_split_center_disp" begin
            eq_split_center_disp = _init(model, :eq_split_center_disp)
            for q in centers, r in q.prod_out, c in r.comp, t in T
                eq_split_center_disp[q, r, c, t] = @constraint(
                    model,
                    z_center_disp[q, r, c, t] ==
                    ratio(q, r, c, t) * z_center_disp_total[q, r, t]
                )
            end
        end
        @timeit "eq_split_store" begin
            eq_split_store = _init(model, :eq_split_store)
            for q in centers, r in q.prod_out, c in r.comp, t in T
                eq_split_store[q, r, c, t] = @constraint(
                    model,
                    z_store[q, r, c, t] == ratio(q, r, c, t) * z_store_total[q, r, t]
                )
            end
        end
        
        @timeit "eq_center_disposal" begin
            eq_center_disposal = _init(model, :eq_center_disposal)
            for r in products, t in T
                centers_r = [q for q in centers if r ∈ q.prod_out]
                if isempty(centers_r)
                    continue
                end
                eq_center_disposal[r, t] = @constraint(
                    model,
                    sum(z_center_disp_total[q, r, t] for q in centers_r) <=
                    data.m_center_disp[r, t]
                )
            end
        end
       
        @timeit "eq_center_storage" begin
            eq_center_storage = _init(model, :eq_center_storage)
            for q in centers, r in q.prod_out, t in T
                eq_center_storage[q, r, t] =
                    @constraint(model, z_store_total[q, r, t] <= data.m_store[q, r, t])
            end
        end
        @timeit "eq_tr_emissions" begin
            eq_tr_emissions = _init(model, :eq_tr_emissions)
            for (q, p, r) in E, s in emissions, t in T
                eq_tr_emissions[q, p, r, s, t] = @constraint(
                    model,
                    z_tr_emissions[q, p, r, s, t] ==
                    data.m_dist[q, p] *
                    data.alpha_tr_emission[r, s, t] *
                    y_total[q, p, r, t]
                )
            end
        end
            
        #=
        @timeit "eq_plant_capacity" begin
            eq_plant_capacity = _init(model, :eq_plant_capacity)
            for p in plants, t in T
                eq_plant_capacity[p, t] = @constraint(
                    model,
                    sum(y_total[q, p, r][t] for (q, r) in E_in[p]) <=
                    data.m_cap[p] * x_open[p, t]
                )
            end
        end
        =#
            
        @timeit "eq_plant_prod" begin
            eq_plant_prod = _init(model, :eq_plant_prod)
            for p in plants, r_out in p.prod_out, c_out in r_out.comp, t in T
                eq_plant_prod[p, r_out, c_out, t] = @constraint(
                    model,
                    z_prod[p, r_out, c_out, t] == sum(
                        data.alpha_mix[p, r_out, c_in, c_out] * y[q, p, r_in, c_in, t]
                        for (q, r_in) in E_in[p] for
                        c_in in r_in.comp if data.alpha_mix[p, r_out, c_in, c_out] > 0
                    )
                )
            end
        end
        @timeit "eq_plant_disp" begin
            eq_plant_disp = _init(model, :eq_plant_disp)
            for p in plants, r in p.prod_out, t in T
                eq_plant_disp[p, r, t] = @constraint(
                    model,
                    sum(z_plant_disp[p, r, c, t] for c in r.comp) <=
                    data.m_plant_disp[p, r, t] * x_disp[p, r, t]
                )
            end
        end
        @timeit "eq_plant_emissions" begin
            eq_plant_emissions = _init(model, :eq_plant_emissions)
            for p in plants, s in emissions, t in T
                eq_plant_emissions[p, s, t] = @constraint(
                    model,
                    z_plant_emissions[p, s, t] ==
                    sum(y_total[q, p, r, t] for (q, r) in E_in[p]) *
                    data.alpha_plant_emission[p, s, t]
                )
            end
        end
        @timeit "eq_emissions_limit" begin
            eq_emissions_limit = _init(model, :eq_emissions_limit)
            for s in emissions, t in T
                eq_emissions_limit[s, t] = @constraint(
                    model,
                    sum(z_plant_emissions[p, s, t] for p in plants) +
                    sum(z_tr_emissions[q, p, r, s, t] for (q, p, r) in E) <=
                    data.m_emission[s, t]
                )
            end
        end
        
        @timeit "eq_plant_remains_open" begin
            eq_plant_remains_open = _init(model, :eq_plant_remains_open)
            for p in plants, t in T
                eq_plant_remains_open[p, t] =
                    @constraint(model, x_open[p, t] >= x_open[p, t-1])
            end
        end
        
        #=
        @timeit "eq_plant_single_dest" begin
            eq_plant_single_dest = _init(model, :eq_plant_single_dest)
            for p in plants, r in p.prod_out, t in T
                eq_plant_single_dest[p, r, t] = @constraint(
                    model,
                    sum(x_send[p, q, r, t] for (q, r2) in E_out[p] if r == r2) +
                    x_disp[p, r, t] <= 1
                )
            end
        end
        
        @timeit "eq_plant_send_limit" begin
            eq_plant_send_limit = _init(model, :eq_plant_send_limit)
            for p in plants, (q, r) in E_out[p], t in T
                eq_plant_send_limit[p, q, r, t] = @constraint(
                    model,
                    y_total[p, q, r][t] <= data.m_cap[q] * x_send[p, q, r, t]
                )
            end
        end
        =#
            
        @timeit "eq_plant_balance" begin
            eq_plant_balance = _init(model, :eq_plant_balance)
            for p in plants, r in p.prod_out, c in r.comp, t in T
                eq_plant_balance[p, r, c, t] = @constraint(
                    model,
                    z_prod[p, r, c, t] ==
                    z_plant_disp[p, r, c, t] +
                    sum(y[p, q, r, c, t] for (q, r2) in E_out[p] if r == r2)
                )
            end
        end

        ############################################### FOR CAPACITY EXPANSION ###########################################
        #Total amount send to a plant should satisfy the capacity constraints of the plant
        
        @timeit "eq_plant_capacity" begin
            eq_plant_capacity = _init(model, :eq_plant_capacity)
            for p in plants, t in T
                eq_plant_capacity[p, t] = @constraint(
                    model,
                    sum(y_total[q, p, r, t] for (q, r) in E_in[p]) <=
                    z_capTot[p, t] 
                )
            end
        end
        
        # total plant capacity
        @timeit "eq_plant_totalcapacity" begin
            eq_plant_totalcapacity = _init(model, :eq_plant_totalcapacity)
            for p in plants, t in T
                eq_plant_totalcapacity[p, t] = @constraint(
                    model,
                    z_capTot[p, t] == z_capTot[p, t-1] + z_capExp[p, t] + (data.min_cap[p] * x_build[p, t])
                )
            end
        end
         
        #Expansion capacity limitation constraints for plants
        @timeit "eq_plant_Expcapacity" begin
            eq_plant_ExpCapacity = _init(model, :eq_plant_Expcapacity)
            for p in plants, t in T
                eq_plant_ExpCapacity[p, t] = @constraint(
                    model,
                    z_capExp[p, t] <= (data.max_cap[p] - data.m_cap[p]) * x_open[p, t]
                )
            end
        end
        
        #Plant Feasibility constraint - A plant is open iff it was open at t-1or new plant is built at t
        @timeit "eq_plant_feasibility" begin
            eq_plant_feasibility = _init(model, :eq_plant_feasibility)
            for p in plants, t in T
                eq_plant_feasibility[p, t] = @constraint(
                    model,
                    x_open[p, t] == x_open[p, t-1] + x_build[p, t]
                )
            end
        end
        
        #Plant Initial opening condition Constraints
        @timeit "eq_plant_opencondition" begin
            eq_plant_opencondition = _init(model, :eq_plant_opencondition)
            for p in plants, t in T
                eq_plant_opencondition[p, t] = @constraint(
                    model,
                    x_build[p, t] <= 1 - x_open[p, t-1] 
                )
            end
        end
        
        #Total capacity can only increase or stay the same
        @timeit "eq_plant_totalcap_continuity" begin
            eq_plant_totalcap_continuity = _init(model, :eq_plant_totalcap_continuity)
            for p in plants, t in T
                #if t>1
                    eq_plant_totalcap_continuity[p, t] = @constraint(
                        model,
                        z_capTot[p, t] >= z_capTot[p, t-1] 
                    )
                #end
            end
        end
        ################################################################################################################
            
    end

    # Optimize
    # -------------------------------------------------------------------------
     _set_names!(model)
    
    model_build_time = time() - init_time
    #write_to_file(model, "model.lp")
    #return
    #error("STOP")
    
    optimize!(model)

    stats = OrderedDict{String,Any}(
        "Plants" => length(plants),
        "Centers" => length(centers),
        "Products" => length(products),
        "Time periods" => length(T),
        "Model build time (s)" => model_build_time,
        "Variables" => num_variables(model),
        "Constraints" =>
            num_constraints(model, count_variable_in_set_constraints=false),
        "Objective Value" => objective_value(model),
        "Solve time (s)" => solve_time(model),
    )
    
    write_to_file(model, "model.lp")
   
    #open("model.lp", "w") do f
    #print(f, model)
    #end
    #write_to_file(model, "model.ilp")
   
    #=
    #Add check for composition consistency
    #-----------------------------------------------------------------------------
    pdt_ratio = dict((p, r, c, t) => 0.0 for p in plants, r in p.prod_out, c in r.comp, t in data.T)
    transp_ratio = dict((p, q, r, c, t) => 0.0 for p in plants, (q,r) in E_out[p], c in r.comp, t in T)
    
    #compute the proportion of components produced in the plants
    for p in plants, r in p.prod_out, c in r.comp, t in T
        pdt_ratio[p,r,c,t] = value(z_prod[p,r,c,t])/sum(value(z_prod[p,r,d,t]) for d in r.comp)
    end
    #compute the proportion of components transported from a plant to another plant
    for p in plants, (q,r) in E_out[p], c in r.comp, t in T
        transp_ratio[p,q,r,c,t] = y[q,p,r][c,t]/sum(value(y[q,p,r][d,t]) for d in r.comp)
    end
    for p in plants, r in p.prod_out, c in r.comp, t in T
        if (pdt_ratio[p,r,c,t] - transp_ratio[q,p,r,c,t] for (q, r2) in E_out[p] if r == r2) >0.1
            @warn "composition inconsistency detected for $r.name transported from $p.name to $q.name at time $t"
        end
    end
    =#
    if output_dir === nothing
        return model, stats
    end

    # Report: Transportation
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."source" = String[]
    df."source longitude (deg)" = Float64[]
    df."source latitude (deg)" = Float64[]
    df."destination" = String[]
    df."destination longitude (deg)" = Float64[]
    df."destination latitude (deg)" = Float64[]
    df."product" = String[]
    df."component" = String[]
    df."time" = Int[]
    df."distance (km)" = Float64[]
    df."amount sent (tonne)" = Float64[]
    df."transportation cost (\$)" = Float64[]
    df."variable operating cost (\$)" = Float64[]
    for (q, p, r) in E, c in r.comp, t in T
        if value(y[q, p, r, c, t]) ≈ 0
            continue
        end
        push!(
            df,
            [
                q.name,
                q.longitude,
                q.latitude,
                p.name,
                p.longitude,
                p.latitude,
                r.name,
                c.name,
                t,
                round(data.m_dist[q, p], digits=2),
                round(value(y[q, p, r, c, t]), digits=2),
                round(
                    data.m_dist[q, p] * data.c_tr[r, t] * value(y[q, p, r, c, t]),
                    digits=2,
                ),
                round(data.c_var[p, t] * value(y[q, p, r, c, t]), digits=2),
            ],
        )
    end
    CSV.write("$output_dir/transp.csv", df)

    # Report: Centers
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."center" = String[]
    df."product" = String[]
    df."component" = String[]
    df."time" = Int[]
    df."amount available (tonne)" = Float64[]
    df."amount sent (tonne)" = Float64[]
    df."amount stored (tonne)" = Float64[]
    df."amount disposed (tonne)" = Float64[]
    df."acquisition cost (\$)" = Float64[]
    df."storage cost (\$)" = Float64[]
    df."disposal cost (\$)" = Float64[]
    for q in centers, r in q.prod_out, c in r.comp, t in T
        push!(
            df,
            [
                q.name,
                r.name,
                c.name,
                t,
                round(data.m_init[q, r, c, t], digits=2),
                round(
                    sum(value(y[q, p, r, c, t]) for (p, r2) in E_out[q] if r == r2; init=0),
                    digits=2,
                ),
                round(value(z_store[q, r, c, t]), digits=2),
                round(value(z_center_disp[q, r, c, t]), digits=2),
                round(data.m_init[q, r, c, t] * data.c_acq[q, r, t], digits=2),
                round(data.c_store[q, r, t] * value(z_store[q, r, c, t]), digits=2),
                round(
                    data.c_center_disp[q, r, t] * value(z_center_disp[q, r, c, t]),
                    digits=2,
                ),
            ],
        )
    end
    CSV.write("$output_dir/centers.csv", df)

    # Report: Plants
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."plant" = String[]
    df."Plant latitude (deg)" = Float64[]
    df."Plant longitude (deg)" = Float64[]
    df."Maximum Capacity" = Float64[]
    df."time(year)" = Int[]
    df."is open?" = Float64[]
    df."is newly built" = Float64[]
    df."Initial Capacity" = Float64[]
    df."Expanded Capacity" = Float64[]
    df."Total Capacity" = Float64[]
    df."opening cost (\$)" = Float64[]
    df."Expansion Cost" = Float64[]
    df."fixed operating cost (\$)" = Float64[]
    df."Exp fixed op cost" = Float64[]
    df."Total fixed operating cost" = Float64[]
    for p in plants, t in T
        push!(
            df,
            [
                p.name,
                p.latitude,
                p.longitude,
                data.max_cap[p],
                t,
                value(x_open[p, t]),
                value(x_build[p, t]),
                data.m_cap[p],
                value(z_capExp[p, t]),
                value(z_capTot[p, t]),
                #data.c_open[p, t] * (value(x_open[p, t]) - value(x_open[p, t-1])),
                data.c_open[p, t] * value(x_build[p, t]),
                data.unit_open[p] * value(z_capExp[p, t]),
                data.c_fix[p, t] * value(x_open[p, t]),
                data.unit_OpC[p] * sum(value(z_capExp[p, i]) for i in 1:t),
                ((data.c_fix[p, t] * value(x_open[p, t])) + (data.unit_OpC[p] * sum(value(z_capExp[p, i]) for i in 1:t))),
            ]
        )
    end
    CSV.write("$output_dir/plants.csv", df)

    # Report: Plant Outputs
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."plant" = String[]
    #df."county" = String[]
    df."Plant Longitude (deg)" = Float64[]
    df."Plant Latitude (deg)" = Float64[]
    df."product" = String[]
    df."component" = String[]
    df."time" = Int[]
    df."amount produced (tonne)" = Float64[]
    df."amount disposed (tonne)" = Float64[]
    df."amount sent (tonne)" = Float64[]
    df."disposal cost (\$)" = Float64[]
    for p in plants, r in p.prod_out, c in r.comp, t in T
        if value(z_prod[p, r, c, t]) ≈ 0
            continue
        end
        push!(
            df,
            [
                p.name,
                #p.county,
                p.longitude,
                p.latitude,
                r.name,
                c.name,
                t,
                round(value(z_prod[p, r, c, t]), digits=2),
                round(value(z_plant_disp[p, r, c, t]), digits=2),
                round(
                    sum(
                        value(y[p, q, r, c, t]) for (q, r2) in E_out[p] if r == r2;
                        init=0.0,
                    ),
                    digits=2,
                ),
                round(
                    data.c_plant_disp[p, r, t] * value(z_plant_disp[p, r, c, t]),
                    digits=2,
                ),
            ],
        )
    end
    CSV.write("$output_dir/plant-outputs.csv", df)

    # Report: Plant Emissions
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."plant" = String[]
    df."emission" = String[]
    df."time" = Int[]
    df."amount emitted (tonne)" = Float64[]
    df."emission cost (\$)" = Float64[]
    for p in plants, s in emissions, t in T
        push!(
            df,
            [
                p.name,
                s.name,
                t,
                round(value(z_plant_emissions[p, s, t]), digits=2),
                round(
                    data.c_emission[s, t] * value(z_plant_emissions[p, s, t]),
                    digits=2,
                ),
            ],
        )
    end
    CSV.write("$output_dir/plant-emissions.csv", df)

    # Report: Transportation Emissions
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."source" = String[]
    df."destination" = String[]
    df."product" = String[]
    df."emission" = String[]
    df."time" = Int[]
    df."distance (km)" = Float64[]
    df."amount sent (tonne)" = Float64[]
    df."amount emitted (tonne)" = Float64[]
    df."emission cost (\$)" = Float64[]
    for (q, p, r) in E, s in emissions, t in T
        if value(y_total[q, p, r, t]) ≈ 0
            continue
        end
        push!(
            df,
            [
                q.name,
                p.name,
                r.name,
                s.name,
                t,
                round(data.m_dist[q, p], digits=2),
                round(value(y_total[q, p, r, t]), digits=2),
                round(value(z_tr_emissions[q, p, r, s, t]), digits=2),
                round(
                    data.c_emission[s, t] * value(z_tr_emissions[q, p, r, s, t]),
                    digits=2,
                ),
            ],
        )
    end
    CSV.write("$output_dir/transp-emissions.csv", df)

    return model, stats
end