using CSV
using Dates
using DataFrames
using Statistics
using LinearAlgebra
using Graphs

# path to data folder 
network = "elec_s_47_ec_lcopt_Co2L-2H"

# construct network
loads = CSV.read("$(@__DIR__)/$network/loads.csv", DataFrame)
lines = CSV.read("$(@__DIR__)/$network/lines.csv", DataFrame)
g = Graph() 
add_vertices!(g, length(nodes))
srcs = indexin(lines.bus0, loads.bus)
dsts = indexin(lines.bus1, loads.bus)
for (u,v) in zip(srcs, dsts)
    add_edge!(g, u,v)
end

## Time series 

# time stamps
snapshots = CSV.read("$(@__DIR__)/$network/snapshots.csv", DataFrame)
snapshots = DateTime.(snapshots.snapshot, "yyyy-mm-dd HH:MM:SS")

# consumption
loads_ts = CSV.read("$(@__DIR__)/$network/loads-p_set.csv", DataFrame)
loads_ts[!, 1] = snapshots
rename!(loads_ts, :Column1 => :snapshot)

# production 
generators = CSV.read("$(@__DIR__)/$network/generators.csv", DataFrame)
generators_ts = CSV.read("$(@__DIR__)/$network/generators-p_max_pu.csv", DataFrame)
generators_ts[!, 1] = snapshots
rename!(generators_ts, :Column1 => :snapshot)


## Renewables
wind = generators[occursin.("wind", generators.carrier), :]
wind_ts = generators_ts[!, ["snapshot", wind.name...]]

solar = generators[occursin.("solar", generators.carrier), :]
solar_ts = generators_ts[!, ["snapshot", solar.name...]]

# aggregate over nodes  
renewables = vcat(wind, solar)
renewables_ts = generators_ts[!, ["snapshot", wind.name..., solar.name...]]
for b in unique(renewables.bus)
    idx = renewables[ renewables.bus .== b, : ].name
    agg = sum.(eachrow(renewables_ts[!, idx]))
    for c in idx 
        select!(renewables_ts, Not(c))
    end
    renewables_ts[!, b] = agg
end


## Calculate averages 
avg_load  = mean( sum.(eachrow(loads_ts[!, 2:end])) )
avg_prod  = mean( sum.(eachrow(generators_ts[!, 2:end])) )
avg_wind  = mean( sum.(eachrow(wind_ts[!, 2:end])) )
avg_solar = mean( sum.(eachrow(solar_ts[!, 2:end])) )
avg_renewables = mean( sum.(eachrow(renewables_ts[!, 2:end])) )

println("Penetration")
println("\tWind\t", avg_wind / avg_load * 100,  "%")
println("\tSolar\t", avg_solar / avg_load * 100, "%")
println("\tRenewables\t", avg_renewables / avg_load * 100, "%")
println("\tAll sources\t", avg_prod / avg_load * 100, "%")
# [Q]: why is this number so low? how much is going into storage?


## Calculate power injection and transmission distribution factor
# mismatch 
mismatch = renewables_ts[!,2:end] .- loads_ts[!,2:end]
Δ⁺ = max.(0,         mismatch)
Δ⁻ = max.(0, (-1) .* mismatch)

# backup energy
b = Δ⁺ .- (mean.(eachcol(loads_ts[!,2:end])) / avg_load)' .* Δ⁻
# [Q]: what about first term? here sum(Δ⁺)/sum(Δ⁺) was left out since this would be 0/0 

# power injection
p = hcat(DataFrame(:snapshot => snapshots), mismatch .- b)

# power transmission distribution factor
Ω = diagm(1 ./ lines.x)
K = incidence_matrix(g, oriented=true)
B = diagm( ones(nv(g)) ) #???
H = Ω * K' * pinv(B)
# [Q]: what is the nodal succeptance matrix B?

# power flow
f = Matrix( hcat([ H * [p[i, 2:end]...]  for i=1:size(p,1)]... )')

## Time averages injection and flow
avg_p = mean.(eachcol(p[!, 2:end]))
avg_f = mean(f, dims=1)[1,:]