text:
# LetsGetTogether.jl
# Uses SurfaceHopping module to combine all functions and start the simulation

using SurfaceHopping
using Base
using Distributions

# Initialise system
plotting = false;
# A
# global hbar_2ma2 = (1.05e-34)^2/(2*9.11e-31*(5.29e-11)^2)/100;
# global e2_4pie0 = (1.6e-19)^2*(8.987551787368e9)/100;
# A_0 = (-hbar_2ma2 +3*e2_4pie0/(1.0e-10))/1.6e-19; # [eV]
# A_1 = -(4*e2_4pie0/(R_0*1e-10))/1.6e-19; # [eV/A]
# A_2 = (2*e2_4pie0/(R_0*1e-10))/1.6e-19; # [eV/A^2]
alpha = 0.55; tau = 0.1;
K = 1.0; global M = 28.07; # 1amu/ps^2 = 0.00010375eV/A^2
#ExtField = 5.0
temp =9000
kb = 8.625000000000001e-5 # eV/Kelvin
R_0  = 4.0
global mean_last
F_min = -0.00; dF = 0.001; F_max = 0.00
#R_n = R_0 + sqrt(kb*temp/K); xr = abs(R_n); xl = 0; v_0 = sqrt(1.6e8/1.66*kb*temp/M);  # [A]
cycles = 1; T = 2e-2; dt = 1e-6; ddt = 1e-6 #
# dt = 1 = 0.1ns (use dt = 1e-10 to balance Angstrom)
n = Int(round(dt/ddt+1)); N = Int(round(T/(2*dt)))

# Set up empty arrays
KE = zeros(N); PE = zeros(N); TE = zeros(N); Ee = zeros(N);  En = zeros(N)
V_change = zeros(N); Empty = zeros(N); Empty2 = zeros(N); Empty3 = zeros(N); Amps = zeros(N);
Empty2 = complex(Empty2); runfreq = zeros(N); Amps2 = zeros(N); temporary = 0; FieldCount = 0
freqvsR = Dict{Float64,Int64}();RateField = Dict{Float64,Float64}();FieldFreqs = Dict{Float64,Int64}();

if Base.Filesystem.isdir("$temp, $F_min:$dF:$F_max")
    Base.Filesystem.rm("$temp, $F_min:$dF:$F_max", recursive = true)
end
Base.Filesystem.mkdir("$temp, $F_min:$dF:$F_max")
Base.Filesystem.cd("$temp, $F_min:$dF:$F_max")

# p1 = plot(apG, color = :blue, label = "Adiabatic PES for ground state", linewidth = 5)
# plot!(p1,apE, color = :blue, label = "Adiabatic PES for excited state", linewidth = 5)
# plot!(p1,dpL, color = :green, label = "Diabatic PES for H--H+", linewidth = 5)
# plot!(p1,dpR, color = :green, label = "Diabatic PES for H+--H", linewidth = 5)
# plot!(p1,title = "PES for avoided crossing", xlabel = "Bond Distance, Angstrom", ylabel = "Energy, eV", leg= ((0.37,0.98)))
#savefig(p1, "Potential Energy Surfaces")
phi_l = diabatic_state(0.);
phi_r(R) = diabatic_state(R);
dphi = diabatic_derivative()
# Calculate coupling term
J_lr(R) = overlap_phi(phi_l,phi_r(R))
# Generate diabatic potential energy surfaces
alpha_k(R) = alpha+0.0*R_0
dpL,dpR,dpX = diabatic_potential(alpha,alpha_k,tau,J_lr)
H_D(R) = H_diabatic(dpL(R),dpR(R),dpX(R))
#Generate adiabatic potential energy surfaces
#apg,ape = adiabatic_potential(dpl,dpr,J_lr)
H_aD(R) = H_adiabatic(H_D(R))
apG(R) = H_aD(R)[1]; apE(R) = H_aD(R)[4]

p2 = plot()
p3 = scatter()
Mean_pos = zeros(Int((F_max-F_min)/dF+1))
for ExtField in F_min:dF:F_max
    FieldCount+=1
    # Generate diabatic electonic states and form of their derivative
    #alpha_k(R) = alpha+ExtField*R

    println("External Field = ", ExtField*1000," mV/Angstrom")
    for l in 1:cycles
        #temp = temp*10
        State_Amplitudes = Array{Float64,}(5000)
        cl = Complex(1.0); cr = Complex(0.0)
        PESg = true
        R_0 = 4.0;
        #dR = rand(Normal(0,2*sqrt(kb*temp/K))); dv = rand(Normal(0,2*sqrt(1.6e8/1.66*kb*temp/M)))
        dR = 2*sqrt(kb*temp/K); dv = 2*sqrt(1.6e8/1.66*kb*temp/M)
        R_n = R_0 + dR; xr = R_n; xl = 0; v_0 = 0+dv;
        global count; count = 0
        for i in 1:N
            phi_r_new = phi_r(R_n)
            dphi = diabatic_derivative()
            J_lr_new = J_lr(R_n)
            dpl,dpr,dpx = diabatic_potential(alpha,alpha+ExtField*abs(R_n),tau,J_lr_new)
            H_d(R) = H_diabatic(dpl(R),dpr(R),dpx(R))
            H_ad(R) = H_adiabatic(H_d(R))
            apg(R) = H_ad(R)[1]; ape(R) = H_ad(R)[4]
            #apg,ape = adiabatic_potential(dpl,dpr,J_lr_new)
            # Determine which adiabatic potential energy surface is the active layer
            if PESg
                f = force(apg)
            else
                f = force(ape)
            end
            # Deteremine which site the electronic state is localised on
            # if PESg & R_n>0
            #     R_e = R_n
            # elseif PESg & R_n>0
            #     R_e = R_n
            # end
            # Perform classical propagation of the nuclei from t -> t+dt
            x,v = classical_propagation(2*dt, R_n, sign(R_n)*R_0, v_0, dt, f, M, K, 0.)
            dx = x[end] - x[1]; dx_l = -dx/2; dx_r = -dx_l

            # Generate necessary parameterd for quantum propagation
            #H_d = H_diabatic(dpl(R_n),dpr(R_n),J_lr_new)
            #H_ad = H_adiabatic(H_d)
            U_nk = diagonalise(H_d(R_n))
            #a_g, a_e = a_mn(U_nk, cl, cr);
            as = *(U_nk',[cl cr]')
            a_g = Complex(as[1]); a_e = Complex(as[2]);
            Empty2[1] = a_g; Empty2[2] = a_e
            #a_g = [a_gn; a_e = [a_en]
            d_ge, d_eg, d_gg, d_ee,ds = NACV(phi_l, phi_r_new, dphi, U_nk, R_n, dx_l, dx_r, 2*dt)
            # Perform quantum propagation of the electronic state from t -> t+dt
            for j in 2:2*n+1
                a_gn, a_en = new_a_mn(Empty2[2*j-3], Empty2[2*j-2], apg, ape, d_ge, d_eg, d_gg, d_ee, R_n, ddt)
                Empty2[j*2-1] = a_gn; Empty2[2*j] = a_en
            end
            a_g = Empty2[4*n+1]; a_e = Empty2[4*n+2]
            as[1] = a_g; as[2] = a_e
            #Empty2[i] = a_g

            #cl = (U_nk[1]*a_g + U_nk[3]*a_e); cr = (U_nk[2]*a_g + U_nk[4]*a_e)
            cs = *(as',U_nk'); cl = cs[1]; cr = cs[2]
            Amps[i] = abs(cl)^2; Amps2[i] = abs(a_g)^2
            #Ee[i] = abs(a_g)^2*apg(R_n) + abs(a_e)^2*ape(R_n)
            # Calculate the probability of hopping
            if PESg
                Empty3[i] = g_mn(a_g, a_e, d_eg, 2*dt);
                En[i] = apg(R_n)
                Ee[i] = apg(R_n)
                V_change[i] = ape(R_n) - Ee[i]
                #FROM g TO e
            else
                Empty3[i] = g_mn(a_e, a_g, d_ge, 2*dt);
                En[i] = ape(R_n)
                Ee[i] = ape(R_n)
                V_change[i] = apg(R_n) - Ee[i]
                 #FROM e TO g
            end
            #Empty3[i] = g
            if Empty3[i]>=0
                g = Empty3[i]
            else
                g = 0
                Empty3[i] = 0
            end

            KE[i] = 1.66/1.6e8*0.5*M*v[end]^2
            PE[i] = 0.5*K*(abs(x[end])-R_0)^2

            # Invoke surface hopping:
            chi = Distributions.Uniform()
            chi_rand = rand(chi)/1e3
            if chi_rand<g
                #println(chi_rand, "    ", g)
                if PESg
                    A = M*d_eg
                    B = M*(ds[3]^2+ds[4]^2)
                    d = (ds[3]+ds[4])
                else
                    A = M*d_ge
                    B = M*(ds[1]^2+ds[2]^2)
                    d = (ds[1]+ds[2])
                end
                #println(V_change[i] ,"      ", KE[i], "      ", 1.66/1.6e8*A^2/(2*B))
                if  V_change[i] < KE[i]
                    runfreq[i]+=1
                    if get(freqvsR, round(abs(R_n),1)+10*(ExtField-1), 0)==0
                        freqvsR[round(abs(R_n),1)+10*(ExtField-1)]=1
                    else
                        freqvsR[round(abs(R_n),1)+10*(ExtField-1)]+=1
                    end
                    v_new = sign(real(v[end]))*sqrt(complex(v[end]^2-2*(1.6e8/1.66*V_change[i])/M))#v[end]+d*A/B*(-1+sqrt(complex(1-2*(1.6e8/1.66*V_change[i]*B/A^2))))
                    #println(v_new-v[end])
                    v[end] = -real(v_new)

                    x[end] = -x[end]
                    if PESg
                        println(PESg," SWAP ", g)
                        PESg = false
                    else
                        println(PESg," SWAP ", g)
                        PESg = true
                    end
                end
            end
            # Update parameters for next iteration
            Empty[i] = R_n
            R_n = x[end]; v_0 = v[end]; xr = R_n; xl = 0;

            # if plotting
            #     psi_g, psi_e = adiabatic_states(phi_l, phi_r_new, U_nk)
            #     Psi_r(R) = real(a_g)*psi_g(R) + real(a_e)*psi_e(R)
            #     # plot(apg, color=:blue, -1.5,1.5)
            #     # plot!(ape, color=:blue)
            #     plot(phi_l,-2,7)
            #     plot!(phi_r_new)
            #     plot!(Psi_r)
            #     scatter!([0],[0], markershape=:circle)
            #     scatter!([R_n],[0], markershape=:circle)
            #     ylims!((-1.5,1.5))
            #     gui()
            # end

        end
        println(l)
        if l > 1
            mean_last = mean([mean(Amps),mean_last])
        else
            mean_last = mean(Amps)
        end

        if l == cycles
            APG = [apG(pos) for pos in Empty]
            APE = [apE(pos) for pos in Empty]
            writedlm("data, External Field = $ExtField", [APE,APG,En,Amps,Empty,PE,KE,Ee,V_change])
        end
        temporary+=sum(Empty3)
    end

    Mean_pos[FieldCount] = mean_last
    RateField[ExtField] = temporary/cyclesin
    num = 0

    # Max = sort(collect(freqvsR), by = tuple -> last(tuple))

    # if length(freqvsR)>0
    #     while Max[end][1]<0
    #         Max = Max[1:end-1]
    #     end
    #     Empty[Int(ExtField*10+51)] = Max[end][1]
    #     println(Max[end][1])
    # else
    #     Empty[Int(ExtField*10+51)] = 0
    # end

    for tuple in freqvsR
          num+=tuple[2]
    end
   FieldFreqs[ExtField] = num
end
Ee; V_change; Empty; Empty2; Empty3; runfreq;freqvsR;Amps;Amps2;RateField;
Mean_pos;

writedlm("Num_Hops_per_Field.txt", FieldFreqs)
writedlm("Rate_Hops_per_Field.txt", RateField)
writedlm("Freq_Hops_Histogram.txt", freqvsR)

#plot(KE); plot!(PE); plot!(TE); plot!(Ee)


text:
module Schrodinger
using LinearAlgebra, SparseArrays, Printf
using Base: tail, front, product, promote_eltype
using OrdinaryDiffEq, Optim

# Export submodules
export Gate

# Export types & functions
export Operator, Ket, Bra,
    data, dims, isnormalized, dimsmatch, dense, braket, isapproxhermitian,
    isunitary, isapproxunitary, isdensityop, checkdensityop,
    hermitian, symmetric,
    tensor, ptrace, expect, levelprobs, findstate,
    purity, fidelity, fidelity2, entanglement_fidelity, gate_fidelity,
    fock, basis, coherent, thermal, maxentangled, maxmixed, ket, @qb_str,
    qzero, qeye, create, destroy, numberop, projectorop,
    displacementop, squeezeop,
    rand_unitary,
    Liouvillian, SchrodingerEvo, LindbladEvo,
    Propagator, SchrodingerProp, LindbladProp,
    sesolve, mesolve, lsolve, psolve, psteady,
    grape, NormPSU, CoherentSubspaces,
    expim, gaussian, inner, scale, scale!

# re-export some LinearAlgebra methods
export eigen, eigs, eigvals, trace, normalize, normalize!

include("basepatch/v1.0.jl")
@static if VERSION < v"1.1"
    using Future: copy!
end

include("quobj/types.jl")
include("quobj/basicmethods.jl")
include("math/basemath.jl")
include("math/eigen.jl")
include("math/special.jl")
include("math/superops.jl")
include("math/norms.jl")
include("math/ptrace.jl")
include("dynamics/liouvillian.jl")
include("dynamics/propagator.jl")
include("dynamics/constructors.jl")
include("dynamics/interface.jl")
include("control/grape.jl")
include("control/normpsu.jl")
include("control/coherentsubspaces.jl")
include("control/amplitudelimiter.jl")
include("tomography/statetomo.jl")
include("tomography/processtomo.jl")
include("misc/utils.jl")
include("misc/indexing.jl")
include("misc/tupletools.jl")
include("misc/inner.jl")
include("misc/kron.jl")
include("misc/sparsevec.jl")
include("misc/checks.jl")
include("misc/hermfact.jl")
include("library/operators.jl")
include("library/states.jl")
include("library/random.jl")
include("library/drivefuns.jl")
include("library/constants.jl")
include("library/gates.jl")

include("quobj/gpu.jl")

end

ext: jl


text:
module AMR

export sample_function

include("utils.jl")
include("sampling.jl")

end # module


text:
struct TravisToken end

base_url(t::TravisToken) = "https://api.travis-ci.org"
headers(t::TravisToken) = Dict(
    "Accept" => "application/vnd.travis-ci.2+json",
    "User-Agent" => "Travis/1.6.8")

function get_travis_token(github_token)
    @info "getting travis token"
    json_parse(talk_to(HTTP.post, TravisToken(),
        "/auth/github",
        Dict("github_token" => github_token)
    ))["access_token"]
end


text:
module FractionalFlow

using Roots, Dierckx, PyPlot, JFVM, ProgressMeter, JSON
import Optim
# import GR
include("CoreyFunctions.jl")
include("water_flood.jl")
include("low_sal_flood.jl")
include("water_solvent_flood.jl")
include("tertiary_flood.jl")
include("water_flood_fvm.jl")
include("water_flood_fvm_upwind.jl")
include("polymer_flood_fvm_upwind.jl")
include("water_flood_upwind_BL.jl")
include("fractional_flow_cases.jl")
include("HTMLElements.jl")
include("fractional_flow_input.jl")

CF = CoreyFunctions # a shorter name

# types:
"""
Core relative permeability parameters
    krw0::Float64
    kro0::Float64
    swc::Float64
    sor::Float64
    nw::Float64
    no::Float64
"""
struct CoreyRelativePermeability
    krw0::Float64
    kro0::Float64
    swc::Float64
    sor::Float64
    nw::Float64
    no::Float64
end

"""
The properties of the core are stored in this structure
    length::Float64
    diameter::Float64
    porosity::Float64
    permeability::Float64
"""
struct CoreProperties
    length::Float64
    diameter::Float64
    porosity::Float64
    permeability::Float64
    surface_area::Float64
    matrix_density::Float64
end

"""
species concentration
    species::AbstractString
    molality::Float64
"""
struct SpeciesConcentration
    species::AbstractString
    concentration::Float64
end

"""
Defines a brine solution, phreeqc style
    name::AbstractString
    number::Integer
    temperature::Float64
    pressure::Float64
    unit::AbstractString
    salinity::Array{SpeciesConcentration, 1}
    pH::Float64
"""
struct Brine
    name::AbstractString
    number::Integer
    temperature::Float64
    pressure::Float64
    unit::AbstractString
    salinity::Array{SpeciesConcentration, 1}
    pH::Float64
    density::Float64
    viscosity::Float64
end

struct Oil
    density::Float64
    viscosity::Float64
    acid_number::Float64
    base_number::Float64
    asphaltene::Float64
end

"""
A structure for storing the core flooding experimental conditions
    injection_velocity::Float64        # m/s
    injected_pore_volume::Float64
    back_pressure::Float64             # Pa
    initial_water_saturation::Float64
    injected_water_saturation::Float64
"""
struct CoreFlooding
    injection_velocity::Float64 # m/s
    injected_pore_volume::Float64
    back_pressure::Float64
    initial_water_saturation::Float64
    injected_water_saturation::Float64
end

"""
a structure that stores the viscosity of the fluids
"""
struct Fluids
    oil_viscosity::Float64
    water_viscosity::Float64
end

"""
A line to be used for storing the shock line results
"""
struct Line
    start_point::Array{Float64,1}
    end_point::Array{Float64,1}
end

"""
A structure that stores the fractional flow curves and the final solution
    fractional_flow_functions::Array{Function, 1}
    shock_line::Array{Float64, 2}
    recovery_pv::Array{Float64, 2}
    recovery_time::Array{Float64, 2}
    saturation_profile_xt::Array{Float64, 2}
"""
struct FracFlowResults
    fractional_flow_functions::Array{Function, 1}
    shock_lines::Array{Line, 1}
    recovery_pv::Array{Float64, 2}
    recovery_time::Array{Float64, 2}
    saturation_profile_xt::Array{Float64, 2}
    tracer_profile_xt::Array{Float64, 2}
    dp_pv::Array{Float64, 2}
    dp_time::Array{Float64, 2}
    water_cut_pv::Array{Float64, 2}
    water_cut_time::Array{Float64, 2}
end

"""
Exp_data stores the core flooding experimental data, recovery factor and pressure drop
    t_exp_dp: time [s] at which the pressure drop is measured
    dp_exp: pressure drop [Pa] measurements
    t_exp_R: time at which recovery factor is measured
    R_exp: recovery factor
"""
struct Exp_data
    t_exp_dp
    dp_exp
    t_exp_R
    R_exp
end

# functions

"""
from old julia
"""
function linspace(a,b,c)
    return range(a, stop=b, length=c)
end

function indmin(a)
    return argmin(a)
end

function indmax(a)
    return argmax(a)
end

"""
visualize(rel_perm::CoreyRelativePermeability)
This function visualizes the relative permeability curve that is defined
by a CoreyRelativePermeability structure
"""
function visualize(rel_perms::CoreyRelativePermeability; label="")
    # convert the rel_perm structure to the rel_perm functions
    krw, kro, dkrwdsw, dkrodsw = rel_perm_functions(rel_perms)
    sw = linspace(0,1, 100)
    plot(sw, krw.(sw), label="krw-$label")
    plot(sw, kro.(sw), label="kro-$label")
    xlabel("Sw")
    ylabel("Relative permeability")
    legend()
end

"""
tabulate the fluids viscosity
"""
function print_fluids(fluids::Fluids; title = nothing)
    HTMLElements.table([fluids.oil_viscosity fluids.water_viscosity],
        header_row=[:mu_water, :mu_oil], title=title)
end

"""
tabulate the core properties
"""
function print_core_properties(core_props::CoreProperties; title=nothing)
    core_param = [core_props.length core_props.diameter core_props.porosity core_props.permeability]
    HTMLElements.table(core_param, header_row=[:length, :diameter, :porosity, :permeability], title=title)
end

"""
tabulate the ralative permeability table in a readable format
"""
function print_relperm(rel_perm::CoreyRelativePermeability; title=nothing)
    rel_perm_params = [rel_perm.krw0 rel_perm.kro0 rel_perm.nw rel_perm.no rel_perm.swc rel_perm.sor]
    HTMLElements.table(rel_perm_params, header_row = [:krw0, :kro0, :nw, :no, :Swc, :Sor], title=title)
end

"""
tabulates the core flooding conditions
"""
function print_core_flood(core_flood::CoreFlooding; title=nothing)
    core_flood_param = [core_flood.injection_velocity core_flood.injected_pore_volume core_flood.initial_water_saturation core_flood.injected_water_saturation]
    header_row = [:u_inj_m_s, :pv_injected, :Sw_init, :Sw_inj]
    HTMLElements.table(core_flood_param, header_row = header_row, title=title)
end

"""
This function visualizes the fractional flow curves
"""
function visualize(rel_perms::CoreyRelativePermeability, fluids::Fluids; label="")
    # convert the rel_perm structure to the rel_perm functions
    fw, dfw = fractional_flow_function(rel_perms, fluids)
    sw = linspace(0,1, 100)
    subplot(2,1,1)
    plot(sw, fw.(sw), label="fw-$label")
    xlabel("Sw [-]")
    ylabel("Fractional flow [-]")
    legend()
    subplot(2,1,2)
    plot(sw, dfw.(sw), label="dfw/dSw-$label")
    xlabel("Sw [-]")
    ylabel("Fractional flow derivative [-]")
    legend()
end

"""
Only visualizes the solution procedure
"""
function visualize_solution(wf_res::FracFlowResults)
    # plot the fractional flow functions and the solution procedure
    fw = wf_res.fractional_flow_functions
    sw = linspace(0, 1, 100)
    for f in fw
        plot(sw, f.(sw))
    end
    shock_lines = wf_res.shock_lines
    for sl in shock_lines
        plot(sl.start_point[1], sl.start_point[2], "ro")
        plot(sl.end_point[1], sl.end_point[2], "ro")
        plot([sl.start_point[1], sl.end_point[1]], [sl.start_point[2], sl.end_point[2]])
    end
    xlabel("Water saturation [-]")
    ylabel("water fractional flow [-]")
    legend()
end

"""
visualize the saturation and concentration profiles
"""
function visualize_profiles(wf_res::FracFlowResults)
    plot(wf_res.saturation_profile_xt[:,1], wf_res.saturation_profile_xt[:,2], label = "Sw")
    plot(wf_res.tracer_profile_xt[:,1], wf_res.tracer_profile_xt[:,2], label = "tracer")
    xlabel("xD/tD [-]")
    ylabel("Water saturation")
    legend()
end

"""
visualize the fractional flow results
"""
function visualize(wf_res::FracFlowResults)
    # plot the fractional flow functions and the solution procedure
    fw = wf_res.fractional_flow_functions
    sw = linspace(0, 1, 100)
    figure()
    for f in fw
        plot(sw, f.(sw))
    end
    shock_lines = wf_res.shock_lines
    for sl in shock_lines
        plot(sl.start_point[1], sl.start_point[2], "ro")
        plot(sl.end_point[1], sl.end_point[2], "ro")
        plot([sl.start_point[1], sl.end_point[1]], [sl.start_point[2], sl.end_point[2]])
    end
    xlabel("Water saturation [-]")
    ylabel("water fractional flow [-]")
    # plot the recovery factors (pv)
    figure()
    plot(wf_res.recovery_pv[:,1], wf_res.recovery_pv[:,2])
    xlabel("Pore volume [-]")
    ylabel("Recovery factor")
    figure()
    plot(wf_res.dp_pv[:,1], wf_res.dp_pv[:,2])
    xlabel("Pore volume [-]")
    ylabel("Pressure drop [Pa]")
    figure()
    plot(wf_res.water_cut_pv[:,1], wf_res.water_cut_pv[:,2])
    xlabel("Pore volume [-]")
    ylabel("Water cut [-]")
    # plot the recovery factors (time)
    # figure()
    # plot(wf_res.recovery_time[:,1], wf_res.recovery_time[:,2])
    # xlabel("time [s]")
    # ylabel("Recovery factor")
    # plot the saturation profile
    figure()
    plot(wf_res.saturation_profile_xt[:,1], wf_res.saturation_profile_xt[:,2], label = "Sw")
    plot(wf_res.tracer_profile_xt[:,1], wf_res.tracer_profile_xt[:,2], label = "tracer")
    xlabel("xD/tD [-]")
    ylabel("Water saturation")
    legend()
end

"""
corey_rel_perm = oil_water_rel_perms(krw0=0.4, kro0=0.9,
    swc=0.15, sor=0.2, nw=2.0, no = 2.0)
This functions defines a CoreyRelativePermeability structure
with a set of predefined values
"""
function oil_water_rel_perms(;krw0=0.4, kro0=0.9,
    swc=0.15, sor=0.2, nw=2.0, no = 2.0)
    return CoreyRelativePermeability(krw0, kro0, swc, sor, nw, no)
end

function oil_water_fluids(;mu_water=1e-3, mu_oil=2e-3)
    return Fluids(mu_oil, mu_water)
end

function core_flooding(;u_inj=1.15e-5, pv_inject=5.0, p_back=1e5, sw_init=0.2, sw_inj=1.0, rel_perms=[])
    swc = sw_init
    if typeof(rel_perms)==CoreyRelativePermeability
        if rel_perms.swc>sw_init
            swc = rel_perms.swc
        end
    end
    return CoreFlooding(u_inj, pv_inject, p_back, swc, sw_inj)
end

function core_properties(;L=0.15, D=0.03, φ=0.3, k=1e-12, a=2000, ρ=2700)
    return CoreProperties(L, D, φ, k, a, ρ)
end

"""
rel_perm_functions(rel_perm::CoreyRelativePermeability)
returns:
     krw, kro, dkrwdsw, dkrodsw
"""
function rel_perm_functions(rel_perm::CoreyRelativePermeability)
    kro0 = rel_perm.kro0
    krw0 = rel_perm.krw0
    sor  = rel_perm.sor
    swc = rel_perm.swc
    no = rel_perm.no
    nw = rel_perm.nw

    kro = sw -> CF.kro(sw, kro0, sor, swc, no)
    krw = sw -> CF.krw(sw, krw0, sor, swc, nw)
    dkrwdsw = sw -> CF.dkrwdsw(sw, krw0, sor, swc, nw)
    dkrodsw = sw -> CF.dkrodsw(sw, kro0, sor, swc, no)

    return krw, kro, dkrwdsw, dkrodsw
end

function fractional_flow_function(rel_perms, fluids)
    krw, kro, dkrwdsw, dkrodsw = rel_perm_functions(rel_perms)
    muw = fluids.water_viscosity
    muo = fluids.oil_viscosity

    fw = sw -> ((krw(sw)/muw)/(krw(sw)/muw+kro(sw)/muo))
    dfw = sw -> ((dkrwdsw(sw)/muw*(krw(sw)/muw+kro(sw)/muo)-
    (dkrwdsw(sw)/muw+dkrodsw(sw)/muo)*krw(sw)/muw)/
    (kro(sw)/muo+krw(sw)/muw)^2)

    return fw, dfw
end

function tangent_line_saturation(rel_perms, fluids, sw_fw_point)
    eps1 = 1e-5 # distance from the end point saturations
    swc = rel_perms.swc
    sor = rel_perms.sor
    sw0 = sw_fw_point[1]
    fw0 = sw_fw_point[2]
    fw, dfw = fractional_flow_function(rel_perms, fluids)

    f_shock = sw -> (1-((fw(sw)-fw0)/(sw-sw0))/dfw(sw))
    sw_tmp = collect(linspace(swc+eps1, 1-sor-eps1, 1000))
    # ind_min = indmin(abs.(f_shock.(sw_tmp)))
    ind_max = indmax(dfw.(sw_tmp))
    sw_max = sw_tmp[ind_max]
    res = Optim.optimize(x->-dfw(x), swc, 1-sor, Optim.GoldenSection())
    sw_max = res.minimizer
    sw_tmp = linspace(sw_max, 1-sor-eps1, 100) # sw_tmp[ind_max:end]
    f_tmp = f_shock.(sw_tmp)
    # println(f_tmp)
    sw_right = 1-sor-eps1
    try
        sw_right = find(f_tmp.<0.0)[1]
    catch
        sw_right = 1-sor-eps1
        # @info "difficulty finding the shock front saturation range!"
    end
    ind_min = indmin(abs.(f_shock.(sw_tmp)))
    # ind_shock = indmax(abs.((fw.(sw_tmp[ind_max:end])-fw0)./(sw_tmp[ind_max:end]-sw0)))
    sw_shock_est = sw_tmp[ind_min]
    # println(sw_shock_est)
    sw_shock = 0.0
    try
        if f_shock(sw_max)*f_shock(sw_right)>0
            sw_shock = fzero(f_shock, sw_shock_est) #  ftol = 1e-10, xtol = 1e-7)
            if (abs(f_shock(sw_shock))>abs(f_shock(sw_shock_est))) || (sw_shock<sw_max) || (sw_shock>1-sor-eps1)
                sw_shock = sw_shock_est
            end
        else
            sw_shock = fzero(f_shock, [sw_max,sw_right])
            if abs(f_shock(sw_shock))>abs(f_shock(sw_shock_est))
                sw_shock = sw_shock_est
            end
        end
    catch
        sw_shock = sw_shock_est
        # @info "shock front saturation is estimated: $sw_shock, error is $(f_shock(sw_shock))"
    end
    # println(abs(f_shock(sw_shock_est)))
    # println(abs(f_shock(sw_shock)))
    return sw_shock
end

function cross_point_saturation(fw, rel_perms, point1, point2; sw_left=0, sw_right=0)
    # create the equation for line
    sor = rel_perms.sor
    swc = rel_perms.swc
    sw_tmp = collect(linspace(swc+eps(), 1-sor-eps(), 1000))
    m = (point2[2]-point1[2])/(point2[1]-point1[1])
    b = point1[2]-m*point1[1]
    f = sw -> m*sw+b-fw(sw)
    f_tmp = f.(sw_tmp)
    sw_neg = sw_tmp[indmin(f_tmp)]
    sw_pos = sw_tmp[indmax(f_tmp)]
    if sw_right==0
        sw_left, sw_right = minmax(sw_neg, sw_pos)
    end
    sw_cross = sw_tmp[indmin(abs.(f_tmp))]
    try
        if f(sw_left)*f(sw_right)>0
            sw_cross = fzero(f, sw_cross)
        else
            sw_cross = fzero(f, [sw_left, sw_right])
        end
    catch
        # @info "Cross point is estimated: $sw_cross, error is $(f(sw_cross))"
    end
    return sw_cross
end

"""
find the saturation at the outlet of the core at the end of injection
"""
function outlet_saturation(pv_inj, sw_shock, dfw, rel_perms)
    eps1 = 1e-2
    sor = rel_perms.sor
    sw_tmp = collect(linspace(sw_shock+eps(), 1-sor-eps1, 1000))
    f = sw -> (pv_inj-1/dfw(sw)) # find the outlet saturation at the end of injection
    f_tmp = f.(sw_tmp)
    sw_left = sw_shock+eps()
    sw_right = 1-sor-eps1
    # println(f_tmp)
    sw_max = 0.0
    try
        sw_max = sw_tmp[find(f_tmp.<0.0)[1]]
    catch
        sw_max = sw_tmp[indmin(f_tmp)]
    end
    try
        if f(sw_left)*f(sw_right)>0
            sw_max = fzero(f_sw, [sw_left, sw_right])
        else
            sw_max = fzero(f, [sw_left, sw_right])
        end
    catch
        sw_max = sw_tmp[indmin(abs.(f_tmp))]
        # @info "Outlet saturation is estimated: $sw_max, error is $(f(sw_max))"
    end
    return sw_max
end

function trapz(x,y)
    # copied from https://github.com/hwborchers/NumericalMath.jl/blob/master/src/integrate.jl
        n=length(x)
        if (length(y) != n)
            error("Vectors 'x', 'y' must be of same length")
        end
        r=0
        for i in 2:n
            r+=(x[i] - x[i-1])*(y[i] + y[i-1])
        end
        return r/2.0
end

"""
function interp1d_linear(x,y, xi) does a linear interpolation using the
arrays x and y, based on the specified value xi
"""
function interp1d_linear(x, y, xi)
    if xi<=x[1]
        return y[1]
    elseif xi>=x[end]
        return y[end]
    else
        for i in 1:length(x)-1
            if x[i]<=xi<=x[i+1]
                m = (y[i+1]-y[i])/(x[i+1]-x[i])
                return m*(xi-x[i])+y[i]
            end
        end
    end
end

"""
myeps returns 1e-15 as the value of eps; for some reason the internal
eps() does npt work properly for me:
example:
sw+eps()==sw returns true which is not what I want.
"""
function myeps()
    return 1e-15
end

end # module

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
using JLD, BenchmarkTools

text_X_orig = load("test_X_orig.jld", "test_X_orig")

@benchmark load("test_X_orig.jld", "test_X_orig")

d = load("dataset.jld")
d


text:
abstract type Port end

mutable struct Input
    sys::NLStateSpace
    id
    subscription
    Input(sys; id=1, subscription=nothing) = new(sys, id, subscription)
end

mutable struct Output
    sys::NLStateSpace
    id
    subscribers::Array{Input}
    Output(sys; id=1, subscribers=[]) = new(sys, id, subscribers)
end

function connect(in::Input, out::Output)
    out.sys.f = combine(out.sys.f, in.sys.f, in.sys.h)
    out.sys.h = combine(out.sys.h, in.sys.h)
end

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
using Documenter, Shuffle

DocMeta.setdocmeta!(
    Shuffle,
    :DocTestSetup,
    quote
        using Shuffle
        import Random
        using Random: MersenneTwister
    end;
    recursive=true
)

makedocs(
    sitename = "Shuffle.jl",
    modules = [Shuffle],
    pages = [
        "Home" => "index.md"
        "API Reference" => "reference.md"
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

deploydocs(repo = "github.com/Luapulu/Shuffle.jl.git")


text:
mutable struct DensityMatrix{T, d} <: QuantumVariable{T}
    head::Symbol
    id_hash::UInt64
    size::Tuple{Int, Int}
    value::Convex.ValueOrNothing
    vexity::Convex.Vexity
    function DensityMatrix{T}(d::Int) where {T <: Number}
        this = new{T, d}(:DensityMatrix, 0, (d,d), nothing, Convex.AffineVexity())
        this.id_hash = objectid(this)
        Convex.id_to_variables[this.id_hash] = this
        this
    end
    DensityMatrix(d::Int) = DensityMatrix{ComplexF64}(d)
end

Convex.constraints(ρ::DensityMatrix) = [ ρ ⪰ 0, tr(ρ) == 1 ]
Convex.sign(::DensityMatrix) = Convex.ComplexSign()
Convex.vartype(::DensityMatrix) = Convex.ContVar


text:
using PastaQ
using LinearAlgebra
using Test

@testset "linear array" begin
  N = 10
  couplings = PastaQ.lineararray(N)
  @test length(couplings) == 2
  @test length(couplings[1]) == N ÷ 2
  @test length(couplings[2]) == (N ÷ 2 - 1)

  for cycle in couplings
    q = collect(1:N)
    for bond in cycle
      @test bond[1] in q
      deleteat!(q, findfirst(x -> x == bond[1], q))
      @test bond[2] in q
      deleteat!(q, findfirst(x -> x == bond[2], q))
    end
  end

  N = 11
  couplings = PastaQ.lineararray(N)
  @test length(couplings) == 2
  @test length(couplings[1]) == N ÷ 2
  @test length(couplings[2]) == N ÷ 2
  for cycle in couplings
    q = collect(1:N)
    for bond in cycle
      @test bond[1] in q
      deleteat!(q, findfirst(x -> x == bond[1], q))
      @test bond[2] in q
      deleteat!(q, findfirst(x -> x == bond[2], q))
    end
  end
end

@testset "square array" begin
  Lx = 4
  Ly = 4

  couplings = PastaQ.squarearray(Lx, Ly)
  @test length(couplings) == 4
  @test length(couplings[1]) == Lx * (Ly ÷ 2)
  @test length(couplings[2]) == (Lx ÷ 2) * (Ly ÷ 2)
  @test length(couplings[3]) == Lx * (Ly ÷ 2)
  @test length(couplings[4]) == (Lx ÷ 2) * (Ly ÷ 2)

  N = Lx * Ly
  for cycle in couplings
    q = collect(1:N)
    for bond in cycle
      @test bond[1] in q
      deleteat!(q, findfirst(x -> x == bond[1], q))
      @test bond[2] in q
      deleteat!(q, findfirst(x -> x == bond[2], q))
    end
  end

  couplings = PastaQ.squarearray(Lx, Ly; rotated=true)
  for cycle in couplings
    q = collect(1:N)
    for bond in cycle
      @test bond[1] in q
      deleteat!(q, findfirst(x -> x == bond[1], q))
      @test bond[2] in q
      deleteat!(q, findfirst(x -> x == bond[2], q))
    end
  end

  Lx = 5
  Ly = 5

  couplings = PastaQ.squarearray(Lx, Ly)
  @test length(couplings) == 4
  @test length(couplings[1]) == Lx * (Ly ÷ 2)
  @test length(couplings[2]) == (Lx ÷ 2) * Ly
  @test length(couplings[3]) == Lx * (Ly ÷ 2)
  @test length(couplings[4]) == Lx * (Ly ÷ 2)

  N = Lx * Ly
  for cycle in couplings
    q = collect(1:N)
    for bond in cycle
      @test bond[1] in q
      deleteat!(q, findfirst(x -> x == bond[1], q))
      @test bond[2] in q
      deleteat!(q, findfirst(x -> x == bond[2], q))
    end
  end

  couplings = PastaQ.squarearray(Lx, Ly; rotated=true)
  for cycle in couplings
    q = collect(1:N)
    for bond in cycle
      @test bond[1] in q
      deleteat!(q, findfirst(x -> x == bond[1], q))
      @test bond[2] in q
      deleteat!(q, findfirst(x -> x == bond[2], q))
    end
  end
end

@testset "random couplints" begin
  N = 20
  R = 10

  for d in 1:10
    q = collect(1:N)
    bonds = PastaQ.randomcouplings(N, R)
    for bond in bonds
      @test bond[1] in q
      deleteat!(q, findfirst(x -> x == bond[1], q))
      @test bond[2] in q
      deleteat!(q, findfirst(x -> x == bond[2], q))
    end
  end
end


text:
module AutoOffload

using LinearAlgebra, AbstractFFTs, FFTW

@static if Base.find_package("CuArrays") !== nothing
    using CuArrays
    if Float64(CuArrays.CUDAdrv.totalmem(first(CuArrays.CUDAdrv.devices()))) > 1e9
        @info("CUDA support found, automatic GPU acceleration will be enabled.")
        const GPU_SUPPORTED = true
        const AUTO_GPU_SIZE = 100
        cuify(x) = CuArrays.CuArray(x)

        # Piracy, should get upstreamed
        function Base.ldiv!(x::CuArrays.CuArray,_qr::CuArrays.CUSOLVER.CuQR,b::CuArrays.CuArray)
          _x = UpperTriangular(_qr.R) \ (_qr.Q' * reshape(b,length(b),1))
          x .= vec(_x)
        end
    else
        @info("CUDA support not found, GPU acceleration will not be available.")
        const GPU_SUPPORTED = false
        const AUTO_GPU_SIZE = 100
        cuify(x) = x
    end
else
    @info("CUDA support not found, GPU acceleration will not be available.")
    const GPU_SUPPORTED = false
    const AUTO_GPU_SIZE = 100
    cuify(x) = x
end

function accelerated_mul!(X,A,B)
    if GPU_SUPPORTED && size(A,1) > AUTO_GPU_SIZE && b isa Array
        A,B = cuify(A),cuify(B)
        X .= Array(A*B)
    else
        mul!(X,A,B)
    end
end

function accelerated_ldiv!(X,A,B)
    if GPU_SUPPORTED && size(A,1) > AUTO_GPU_SIZE && b isa Array
        _X,A,B = cuify(X),cuify(A),cuify(B)
        X .= Array(ldiv!(_X,A,B))
    else
        ldiv!(X,A,B)
    end
end

function accelerated_factorize!(A)
    if GPU_SUPPORTED && size(A,1) > AUTO_GPU_SIZE && b isa Array
        qr(cuify(A))
    else
        factorize(A)
    end
end

function accelerated_fft!(A)
    if GPU_SUPPORTED && size(A,1) > AUTO_GPU_SIZE && b isa Array
        A .= Array(fft!(cuify(A)))
    else
        FFTW.fft!(A)
    end
end

export accelerated_mul!, accelerated_ldiv!, accelerated_factorize!,
       accelerated_fft!

module Pirate

    # Somehow take over mul!, *, etc to autooffload to GPU?

end

end # module


text:
using ImageInTerminal, ImageCore
using TestImages, ImageTransformations, CoordinateTransformations
using Rotations, OffsetArrays
using SparseArrays
using Test, ReferenceTests

function ensurecolor(f, args...)
    old_color = Base.have_color
    try
        Core.eval(Base, :(have_color = true))
        return @inferred f(args...)
    finally
        Core.eval(Base, :(have_color = $old_color))
    end
end

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

# define some test images
gray_square = colorview(Gray, N0f8[0. 0.3; 0.7 1])
gray_square_alpha = colorview(GrayA, N0f8[0. 0.3; 0.7 1], N0f8[1 0.7; 0.3 0])
gray_line = colorview(Gray, N0f8[0., 0.3, 0.7, 1])
gray_line_alpha = colorview(GrayA, N0f8[0., 0.3, 0.7, 1], N0f8[1, 0.7, 0.3, 0])
rgb_line = colorview(RGB, range(0, stop=1, length=20), zeroarray, range(1, stop=0, length=20))
rgb_line_4d = repeat(repeat(rgb_line', 1, 1, 1, 1), 1, 1, 2, 2)

camera_man = testimage("camera")
lighthouse = testimage("lighthouse")
toucan = testimage("toucan")
lena = testimage("lena_color_256")

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

tests = [
    "tst_colorant2ansi.jl",
    "tst_encodeimg.jl",
    "tst_imshow.jl",
    "tst_baseshow.jl",
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end


text:
using Test
using DaprProto

@testset "DaprProto" begin
end


text:
""" This file executes the CM with the previously calculated, new load profiles"""

using Pkg
if isfile("Project.toml") && isfile("Manifest.toml")
    Pkg.activate(".")
end

###################
# Terminal ########
###################

    usedayperiod = true
    startday = 1
    lastday = 365
    starthour = 1
    lasthour = 8760
	trm = 0	# reduces - or at negative value, increases - available capacity of transmission lines
	use221 = true	# whether the fictional generators with mc 221 are used

    datapath = "data/"
	newpath = "NewData/"
	inputpath = "output_2030/Output_LoadShift/Sens_quart/"
	outputpath = "output_2030/Output_DSOM/DSOM_Sens_quart_trm_0/"

   if usedayperiod == true
        timeperiod = string((startday-1)*24+1)*":"*string((lastday*24))
   else
        timeperiod = string(starthour)*":"*string(lasthour)
   end
    timeperiod_str = replace(string(timeperiod), ":" => "-")

    println()
    # println("Case >" * case * "< selected.")
    println("Calculating for time period >" * timeperiod * "<.")
    println()

###################
# Run CM ##########
###################

    println("Loading required packages: ProgressMeter, CSV, DataFrames, Dates, JuMP, Gurobi.")
    using ProgressMeter
    using CSV
    using DataFrames
    using Dates
    using JuMP
    using Gurobi
    println("Loading of packages complete.")

    include("model/src/Sesam_ptdf_Loadshift.jl") #WHEN USING PTDF
    # include("model/src/Sesam.jl") WHEN USING NO PTDF
    using .Sesam

###################
# Data ############
###################

	println("Reading and processing Data from CSV...")

    # nodes_df_30 =   # Dataframe mit ID, Lat, longitude
    nodes_df_30 = CSV.read(joinpath(newpath, "nodeCoordinates.csv"))
    nodes_df_30 = rename!(nodes_df_30, :node => :id)[:, 1:3]
#	nodes_df_30[:id] = nodes_df_30[:id] .+1

    # load_df_30 =    # Dataframe mit systemload und knotenweisen shares über alle stunden
    load_df_30 = CSV.read(joinpath(newpath, "demandHourly.csv"))
    rename!(load_df_30, 1 => :systemload)
    b = CSV.read(joinpath(newpath, "loadShare.csv"))
#   	b[:Column1] = b[:Column1] .+1
    for i in 1:nrow(b)
        load_df_30[Symbol(b[i, 1])] = repeat([b[i, 2]], 8760)
    end

    exchange_df_30 = CSV.read(joinpath(datapath, "exchange.csv"))

    # lines_df_30 =    # Dataframe mit id, from, to, voltage, reactance, resistance, circuits, pmax
    lines_df_30 = CSV.read(joinpath(newpath, "lines.csv"))
#	lines_df_30[:id] = lines_df_30[:id] .+1
#	lines_df_30[:from] = lines_df_30[:from] .+1
#	lines_df_30[:to] = lines_df_30[:to] .+1

    # powerplants_df_30 = # Dataframe mit unit, node, technology, fuel, capacity, efficiency, emission, varcost, rup, rdn, latitude, longitude
    temp = CSV.read(joinpath(newpath, "convGenData.csv"))
#	temp[:node] = temp[:node] .+ 1
    powerplants_df_30 = DataFrame(unit = String[], node = Int64[], technology = String[], fuel = String[], capacity = Float64[], efficiency = Float64[],
                                    emission = Float64[], varcost = Float64[], rup = Float64[], rdn = Float64[])
	use221 == true ? techs = names(temp)[2:end] : techs = names(temp)[2:end-1]
    for i in 1:nrow(temp)
        for class in names(temp)[2:end]
            if !ismissing(temp[class][i]) && temp[class][i] > 0
                pp_id = "NODE"*string(temp[i,1])*"C"*string(class)
                cap = temp[class][i]
                emission = 0
                efficiency = 1
                varcost = 0
                rup = 0
                rdn = 0
                push!(powerplants_df_30, [pp_id temp[i, 1] string(class) string(class) cap efficiency emission varcost rup rdn])
            end
        end
    end

    # avail_powerplants_df_30
    avail_powerplants_df_30 = DataFrame()
    for i in powerplants_df_30.unit
        avail_powerplants_df_30[Symbol(i)] = ones(8760)
    end

    # mustrun_powerplants_df_30 = # Dataframe mit überall nullen pro generator
    mustrun_powerplants_df_30 = DataFrame()
    for i in powerplants_df_30.unit
        mustrun_powerplants_df_30[Symbol(i)] = zeros(8760)
    end

    # fuelcost_df_30
    fuelcost_df_30 = DataFrame()
    for i in names(temp)[2:end]
        fuelcost_df_30[i] = [parse(Int, string(i))]
    end
    fuelcost_df_30[:CO2] = [0]

    # renewables_df_30 = # Dataframe mit Node, WindOnshore, WindOffshore, SolarPV, RoR, Geothermal
    solar_df_30 = CSV.read(joinpath(newpath, "solarHourly.csv"))
    solarshare_df_30 = CSV.read(joinpath(newpath, "solarShareNodes.csv"))
#	solarshare_df_30[:node] = solarshare_df_30[:node] .+1
    wind_df_30 = CSV.read(joinpath(newpath, "windHourly.csv"))
    windshare_df_30 = CSV.read(joinpath(newpath, "windShareNodes.csv"))
#	windshare_df_30[:node] = windshare_df_30[:node] .+ 1

    renewables_df_30 = DataFrame(node=nodes_df_30[:id])
    renewables_df_30[:SolarPV] = [maximum(solar_df_30[:solarhourly]*solarshare_df_30[i,2]) for i in 1:nrow(solarshare_df_30)]
    renewables_df_30[:Wind] = [maximum(wind_df_30[:windhourly]*windshare_df_30[i,2]) for i in 1:nrow(windshare_df_30)]
    renewables_df_30

    # avail_solar_pv_30 = # Dataframe mit nodalen, stündlichen relativen Verfügbarkeiten

    # avail_solar_pv_df
    avail_solar_pv_df_30 = DataFrame()
    for n in 1:nrow(solarshare_df_30)
        temparray = []
        for h in 1:8760
            renewables_df_30[:SolarPV][n] > 0 ? x = (solarshare_df_30[n,2]*solar_df_30[h,1])/renewables_df_30[n,2] : x = 0
            push!(temparray, x)
        end
        avail_solar_pv_df_30[Symbol(n)] = temparray
    end
	avail_solar_pv_df_30

    avail_wind_df_30 = DataFrame()
    for n in 1:nrow(windshare_df_30)
        temparray = []
        for h in 1:8760
            renewables_df_30[:Wind][n] > 0 ? x = (windshare_df_30[n,2]*wind_df_30[h,1])/renewables_df_30[n,3] : x = 0
            push!(temparray, x)
        end
        avail_wind_df_30[Symbol(n)] = temparray
    end
    avail_wind_df_30

    renewables_avail_dict_30 = Dict(:SolarPV => avail_solar_pv_df_30,
        :Wind => avail_wind_df_30)

		println("Reading of data complete.")
		println()

###################
# Model setup #####
###################

		println("Preparing model technologies, nodes, and grid...")
		prog = Progress(6)

	# Create sets
		# Create set of nodes and lines
		nodes = Nodes(nodes_df_30, load_df_30, exchange_df_30, 1000, 0)
			next!(prog)
		lines = Lines(lines_df_30, 1000, trm)
			next!(prog)
		# Create set of power plants	ED_L_S_opt, ED_p
		pp = PowerPlants(powerplants_df_30, avail_powerplants_df_30, mustrun_powerplants_df_30, fuelcost_df_30)
			next!(prog)
		# Create set of renewables
		res = Renewables(renewables_df_30, renewables_avail_dict_30)
			next!(prog)
		# Create a set of storages
		storages_df = CSV.read(joinpath(datapath, "storages.csv"))
			next!(prog)
		storages_df[:storage] = zeros(32)
		storages_df[:power] = zeros(32)
		storages = Storages(storages_df)
			next!(prog)
		sum(sum((res.infeed[i]) for i in keys(res.infeed)))
	# Create Dict Zone to Gen
		N = nodes.id
		P = pp.unit
		R = res.unit
		S = storages.unit
		Node2PP = Dict((i, []) for i in nodes.id)
		for i in keys(pp.node)
			push!(Node2PP[pp.node[i]], i)
		end
		Node2Res = Dict((i, []) for i in nodes.id)
		for i in keys(res.node)
			push!(Node2Res[res.node[i]], i)
		end
		Node2Stor = Dict((i, []) for i in nodes.id)
		for i in keys(storages.node)
			push!(Node2Stor[storages.node[i]], i)
		end

	println("Preparation complete.")
	println()

	println("Loading new load profiles.")
	Newloads_df = CSV.read(joinpath(inputpath, "Newload.csv"))

	newloads_dict = Dict((Newloads_df[r,1] => Newloads_df[r, 2:end]) for r in 1:nrow(Newloads_df))
	nodes.load = newloads_dict

		session_time = Dates.now()
		session = Dates.format(session_time, "yyyy-mm-dd__HH-MM-SS")*"_trm_" * string(trm)

	# Create session folder for output
		mkdir(joinpath(outputpath, session*"__"*timeperiod_str))

	prog = Progress(length(dayslicer(timeperiod)))

	### Parameter outputs
		println("Exporting parameters.")
		println("Exporting merit order...")
		# Merit order
		output_ED_merit_order_pp = DataFrame(unit=Symbol[],fuel=String[],capacity=Float64[],mc=Float64[])

		try
			for g in pp.unit
					mc = pp.mc[g]
					fuel = pp.fuel[g]
					capacity = pp.capacity[g]
					push!(output_ED_merit_order_pp, [g fuel capacity mc])
			end

			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_merit_order_pp.csv"), output_ED_merit_order_pp)
		catch e
			println("The marginal costs for power plants are not available.")
		end

		output_ED_merit_order_res = DataFrame(unit=Symbol[],fuel=Symbol[],capacity=Float64[],mc=Float64[])

		try
			for r in res.unit
					mc = 0
					fuel = res.fuel[r]
					capacity = res.capacity[r]
					push!(output_ED_merit_order_res, [r fuel capacity mc])
			end

			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_merit_order_res.csv"), output_ED_merit_order_res)
		catch e
			println("The marginal costs for renewables are not available.")
		end

		output_ED_merit_order_storages = DataFrame(unit=Symbol[],fuel=String[],capacity=Float64[],mc=Float64[])

		try
			for s in storages.unit
					mc = 0
					fuel = storages.fuel[s]
					capacity = storages.power[s]
					push!(output_ED_merit_order_storages, [s fuel capacity mc])
			end

			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_merit_order_storages.csv"), output_ED_merit_order_storages)
		catch e
			println("The marginal costs for renewables are not available.")
		end

	println("Export complete.")
	println()

### Initial run of first dayslice
	println("Initial run. Calculating first dayslice: " * string(dayslicer(timeperiod)[1]))
	println()
	println()

	ED_mod, ED_P_opt, ED_P_R_opt, ED_P_R_opt_dist, ED_P_S_opt, ED_D_S_opt,
		ED_L_S_opt, ED_price = EconomicDispatch80(dayslicer(timeperiod)[1], nodes, pp, res, storages)

	CM_mod, CM_ΔP_up_opt, CM_ΔP_dn_opt, CM_ΔP_R_up_opt, CM_ΔP_R_dn_opt,
		CM_Θ_opt, CM_P_flow_opt, CM_price, CM_P_gen_lost_opt, CM_P_load_lost_opt, CM_P_S_up_opt = CongestionManagementPHS80(dayslicer(timeperiod)[1],
		nodes, lines, pp, res, storages, ED_price, ED_P_opt, ED_P_R_opt_dist,
		ED_P_S_opt, ED_D_S_opt, ED_L_S_opt)

	next!(prog)

		C_V_RD_N = DataFrame(time=Int64[], node=Int64[], CostPP=Float64[], CostRes=Float64[], CostStor=Float64[], CostLostLoad=Float64[], CostLostGen=Float64[], CostTotal=Float64[],
			VolumePP=Float64[], VolumeRes=Float64[], VolumeStor=Float64[], VolumeLostLoad=Float64[], VolumeLostGen=Float64[], VolumeTotal=Float64[])
			for t in dayslicer(timeperiod)[1]
				for n in N
					if isempty(Node2PP[n])
						costpp = 0
					else
						costpp = sum((pp.mc[g][1]*CM_ΔP_up_opt[t, g] + (ED_price[t] - pp.mc[g][]) * CM_ΔP_dn_opt[t, g]) for g in Node2PP[n])
					end
					if isempty(Node2Res[n])
						costres = 0
					else
						costres = sum(((ED_price[t]) * CM_ΔP_R_dn_opt[t, res]) for res in Node2Res[n])
					end
					if isempty(Node2Stor[n])
						coststor = 0
					else
						coststor = sum(((ED_price[t]/storages.efficiency[stor]) * CM_P_S_up_opt[t, stor]) for stor in Node2Stor[n])
					end

					costlostload = 1000*CM_P_load_lost_opt[t, n]
					costlostgen = 1000*CM_P_gen_lost_opt[t, n] # ÜBERALL NULL!
					totalcost = costpp + costres + coststor + costlostload + costlostgen

					if isempty(Node2PP[n])
						volumepp = 0
					else
						volumepp = sum((CM_ΔP_up_opt[t, g] - CM_ΔP_dn_opt[t, g]) for g in Node2PP[n])
					end
					if isempty(Node2Res[n])
						volumeres = 0
					else
						volumeres = sum(-(CM_ΔP_R_dn_opt[t, res]) for res in Node2Res[n])
					end
					if isempty(Node2Stor[n])
						volumestor = 0
					else
						volumestor = sum((CM_P_S_up_opt[t, stor]) for stor in Node2Stor[n])
					end

					volumelostload = CM_P_load_lost_opt[t, n]
					volumelostgen = CM_P_gen_lost_opt[t, n]
					totalvolume = volumepp + volumeres + volumestor + volumelostload + volumelostgen

					push!(C_V_RD_N, [t n costpp costres coststor costlostload costlostgen totalcost volumepp volumeres volumestor volumelostload volumelostgen totalvolume])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "C_V_RD_N__"*timeperiod_str*".csv"), C_V_RD_N)

	# Save ED Data
		# Dispatchable power plants
		output_ED_P_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])

		try
			for t in dayslicer(timeperiod)[1]
				for g in pp.unit
					val = ED_P_opt[t, g]
					fuel =  pp.fuel[g]
					push!(output_ED_P_opt, [t g fuel val pp.node[g]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P.csv"), output_ED_P_opt)
		catch e
			println("ED_P_opt does not exist.")
		end

		# RES generation
		output_ED_P_R_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for r in res.unit
					val = ED_P_R_opt[t, r]
					fuel =  res.fuel[r]
					push!(output_ED_P_R_opt, [t r fuel val res.node[r]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P_R.csv"), output_ED_P_R_opt)
		catch e
			println("ED_P_R_opt does not exist.")
		end

		# Maximum available RES infeed from input data
		output_MaxRES = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for r in res.unit
					output = res.infeed[r][t]
					fuel =  res.fuel[r]
					push!(output_MaxRES, [t r fuel output res.node[r]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "MaxRES.csv"), output_MaxRES)
		catch e
			println("res.infeed does not exist.")
		end

		# Generation from storages
		output_ED_P_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for s in storages.unit
					val = ED_P_S_opt[t, s]
					fuel =  storages.fuel[s]
					push!(output_ED_P_S_opt, [t s fuel val storages.node[s]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P_S.csv"), output_ED_P_S_opt)
		catch e
			println("ED_P_S_opt does not exist.")
		end

		# Electricity demand from storages
		output_ED_D_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for s in storages.unit
					val = -1*ED_D_S_opt[t, s]
					fuel =  storages.fuel[s]
					push!(output_ED_D_S_opt, [t s fuel val storages.node[s]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_D_S.csv"), output_ED_D_S_opt)
		catch e
			println("ED_D_S_opt does not exist.")
		end

		# Storage level
		output_ED_L_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],storage=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for s in storages.unit
					val = ED_L_S_opt[t, s]
					fuel =  storages.fuel[s]
					push!(output_ED_L_S_opt, [t s fuel val storages.node[s]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_L_S.csv"), output_ED_L_S_opt)
		catch e
			println("ED_L_S_opt does not exist.")
		end

		# Load
		output_load = DataFrame(time=Int64[],load=Float64[])
		try
			for t in dayslicer(timeperiod)[1]
				load = nodes.systemload[t]
				push!(output_load, [t load])
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "load.csv"), output_load)
		catch e
			println("nodes.systemload does not exist.")
		end

		# Prices
		output_ED_price = DataFrame(time=Int64[],price=Float64[])
		try
			for t in dayslicer(timeperiod)[1]
				price =  ED_price[t]
				push!(output_ED_price, [t price])
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_price.csv"), output_ED_price)
		catch e
			println("ED_price does not exist.")
		end

		# Must-run
		output_pmin = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for g in pp.unit
					val = pp.pmin[g][t]
					fuel =  pp.fuel[g]
					push!(output_pmin, [t g fuel val pp.node[g]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "pmin.csv"), output_pmin)
		catch e
			println("pp.pmin does not exist.")
		end

		### Model parameters
		output_ED_info = DataFrame(timeperiod=String[],
										termination_status=String[],
										primal_status=String[],
										dual_status=String[],
										objective_value=Float64[],
										solve_time=Float64[])

		push!(output_ED_info, [string(dayslicer(timeperiod)[1]) string(termination_status(ED_mod)) string(primal_status(ED_mod)) string(dual_status(ED_mod)) objective_value(ED_mod) solve_time(ED_mod)])
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_info.csv"), output_ED_info)

		# Nodal exchange
		output_exchange = DataFrame(time=Int64[],node=Int64[],output=Float64[])
		try
			for t in dayslicer(timeperiod)[1]
				for j in nodes.id
					push!(output_exchange, [t j nodes.exchange[j][t]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "exchange.csv"), output_exchange)
		catch e
			println("output_exchange does not exist.")
		end

###################
# CM ##############
###################
# Save CM Data
	# Dispatchable power plants
	output_CM_ΔP_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for g in pp.unit
				val = CM_ΔP_up_opt[t, g]
				fuel =  pp.fuel[g]
				push!(output_CM_ΔP_up_opt, [t g fuel val pp.node[g]])
			end
		end
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_up.csv"), output_CM_ΔP_up_opt)
	catch e
		println("CM_ΔP_up_opt does not exist.")
	end

	output_CM_ΔP_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for g in pp.unit
				val = -1*CM_ΔP_dn_opt[t, g]
				fuel =  pp.fuel[g]
				push!(output_CM_ΔP_dn_opt, [t g fuel val pp.node[g]])
			end
		end
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_dn.csv"), output_CM_ΔP_dn_opt)
	catch e
		println("CM_ΔP_dn_opt does not exist.")
	end

	# RES generation
	output_CM_ΔP_R_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for r in res.unit
				val = CM_ΔP_R_up_opt[t, r]
				fuel =  res.fuel[r]
				push!(output_CM_ΔP_R_up_opt, [t r fuel val res.node[r]])
			end
		end
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_R_up.csv"), output_CM_ΔP_R_up_opt)
	catch e
		println("CM_ΔP_R_up_opt does not exist.")
	end

	output_CM_ΔP_R_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for r in res.unit
				val = -1*CM_ΔP_R_dn_opt[t, r]
				fuel =  res.fuel[r]
				push!(output_CM_ΔP_R_dn_opt, [t r fuel val res.node[r]])
			end
		end
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_R_dn.csv"), output_CM_ΔP_R_dn_opt)
	catch e
		println("CM_ΔP_R_dn_opt does not exist.")
	end

	### Model parameters
	output_CM_info = DataFrame(timeperiod=String[],
									   termination_status=String[],
									   primal_status=String[],
									   dual_status=String[],
									   objective_value=Float64[],
									   solve_time=Float64[])

	push!(output_CM_info, [string(dayslicer(timeperiod)[1]) string(termination_status(CM_mod)) string(primal_status(CM_mod)) string(dual_status(CM_mod)) objective_value(CM_mod) solve_time(CM_mod)])
	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_info.csv"), output_CM_info)

	### Power flow results
	output_CM_P_flow_opt = DataFrame(time=Int64[],
									 line=Int64[],
									 from=Int64[],
									 to=Int64[],
									 pmax=Float64[],
									 pflow=Float64[])

	 try
 		for t in dayslicer(timeperiod)[1]
 			for l in lines.id
				from = parse(Int, string(lines.from[l]))
				to = parse(Int, string(lines.to[l]))
				pmax = lines.pmax[l]
				pflow = CM_P_flow_opt[t, l]
 				push!(output_CM_P_flow_opt, [t l from to pmax pflow])
 			end
 		end
 		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_flow.csv"), output_CM_P_flow_opt)
 	catch e
 		println("CM_P_flow_opt does not exist.")
 	end

	# Lost generation
	output_CM_P_gen_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
	try
		for t in dayslicer(timeperiod)[1]
			for j in nodes.id
				output = CM_P_gen_lost_opt[t, j]
				push!(output_CM_P_gen_lost, [t, j, output])
			end
		end
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_gen_lost.csv"), output_CM_P_gen_lost)
	catch e
		println("CM_P_gen_lost does not exist.")
	end

	# Lost load
	output_CM_P_load_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
	try
		for t in dayslicer(timeperiod)[1]
			for j in nodes.id
				output = CM_P_load_lost_opt[t, j]
				push!(output_CM_P_load_lost, [t, j, output])
			end
		end
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_load_lost.csv"), output_CM_P_load_lost)
	catch e
		println("CM_P_load_lost does not exist.")
	end

	output_CM_price = DataFrame(time=Int64[],node=Int64[],price=Float64[])
	 try
		 for t in dayslicer(timeperiod)[1]
			 for j in nodes.id
				 price =  CM_price[j, t]
				 push!(output_CM_price, [t j price])
			 end
		 end
		 CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_price.csv"), output_CM_price)
	 catch e
		 println("CM_price does not exist.")
	 end

	# PHS output increase
	output_CM_P_S_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for s in storages.unit
				val = CM_P_S_up_opt[t, s]
				fuel =  storages.fuel[s]
				push!(output_CM_P_S_up_opt, [t s fuel val storages.node[s]])
			end
		end
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_S_up.csv"), output_CM_P_S_up_opt)
	catch e
		println("CM_P_S_up_opt does not exist.")
	end

### Remaining periods: Initialising first hour with last hour of previous day
# 	and append data

	for hours in dayslicer(timeperiod)[2:end]
		hours_extended = [hours[1]-1; hours]

	println()
	println()
	println("Calculating economic dispatch: " * string(hours))

	global ED_mod, ED_P_opt, ED_P_R_opt, ED_P_R_opt_dist, ED_P_S_opt, ED_D_S_opt,
		ED_L_S_opt, ED_price = EconomicDispatch80(hours, nodes, pp, res, storages)

	    ###################
    # ED ##############
	    ###################

	    # Dispatchable power plants
	    output_ED_P_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])

	    try
	    	for t in hours
	    		for g in pp.unit
	    			val = ED_P_opt[t, g]
	    			fuel =  pp.fuel[g]
	    			push!(output_ED_P_opt, [t g fuel val pp.node[g]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P.csv"), output_ED_P_opt, append = true)
	    catch e
	    	println("ED_P_opt does not exist.")
	    end

	    # RES generation
	    output_ED_P_R_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for r in res.unit
	    			val = ED_P_R_opt[t, r]
	    			fuel =  res.fuel[r]
	    			push!(output_ED_P_R_opt, [t r fuel val res.node[r]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P_R.csv"), output_ED_P_R_opt, append = true)
	    catch e
	    	println("ED_P_R_opt does not exist.")
	    end

	    # Maximum available RES infeed from input data
	    output_MaxRES = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for r in res.unit
	    			output = res.infeed[r][t]
	    			fuel =  res.fuel[r]
	    			push!(output_MaxRES, [t r fuel output res.node[r]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "MaxRES.csv"), output_MaxRES, append = true)
	    catch e
	    	println("res.infeed does not exist.")
	    end

	    # Generation from storages
	    output_ED_P_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for s in storages.unit
	    			val = ED_P_S_opt[t, s]
	    			fuel =  storages.fuel[s]
	    			push!(output_ED_P_S_opt, [t s fuel val storages.node[s]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P_S.csv"), output_ED_P_S_opt, append = true)
	    catch e
	    	println("ED_P_S_opt does not exist.")
	    end

	    # Electricity demand from storages
	    output_ED_D_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for s in storages.unit
	    			val = -1*ED_D_S_opt[t, s]
	    			fuel =  storages.fuel[s]
	    			push!(output_ED_D_S_opt, [t s fuel val storages.node[s]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_D_S.csv"), output_ED_D_S_opt, append = true)
	    catch e
	    	println("ED_D_S_opt does not exist.")
	    end

	    # Storage level
	    output_ED_L_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],storage=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for s in storages.unit
	    			val = ED_L_S_opt[t, s]
	    			fuel =  storages.fuel[s]
	    			push!(output_ED_L_S_opt, [t s fuel val storages.node[s]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_L_S.csv"), output_ED_L_S_opt, append = true)
	    catch e
	    	println("ED_L_S_opt does not exist.")
	    end

	    # Load
	    output_load = DataFrame(time=Int64[],load=Float64[])
	    try
	    	for t in hours
	    		load = nodes.systemload[t]
	    		push!(output_load, [t load])
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "load.csv"), output_load, append = true)
	    catch e
	    	println("nodes.systemload does not exist.")
	    end

	    # Prices
	    output_ED_price = DataFrame(time=Int64[],price=Float64[])
	    try
	    	for t in hours
	    		price =  ED_price[t]
	    		push!(output_ED_price, [t price])
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_price.csv"), output_ED_price, append = true)
	    catch e
	    	println("ED_price does not exist.")
	    end

	    # Must-run
	    output_pmin = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for g in pp.unit
	    			val = pp.pmin[g][t]
	    			fuel =  pp.fuel[g]
	    			push!(output_pmin, [t g fuel val pp.node[g]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "pmin.csv"), output_pmin, append = true)
	    catch e
	    	println("pp.pmin does not exist.")
	    end

		### Model parameters
		output_ED_info = DataFrame(timeperiod=String[],
										   termination_status=String[],
										   primal_status=String[],
										   dual_status=String[],
										   objective_value=Float64[],
										   solve_time=Float64[])

		push!(output_ED_info, [string(hours) string(termination_status(ED_mod)) string(primal_status(ED_mod)) string(dual_status(ED_mod)) objective_value(ED_mod) solve_time(ED_mod)])
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_info.csv"), output_ED_info, append = true)

		###################
# CM ##############
		###################

		println()
		println()
		println("Calculating congestion management: " * string(hours))

			global CM_mod, CM_ΔP_up_opt, CM_ΔP_dn_opt, CM_ΔP_R_up_opt, CM_ΔP_R_dn_opt,
				CM_Θ_opt, CM_P_flow_opt, CM_price, CM_P_gen_lost_opt, CM_P_load_lost_opt, CM_P_S_up_opt = CongestionManagementPHS80(hours, nodes, lines, pp, res, storages, ED_price,
				ED_P_opt,  ED_P_R_opt_dist, ED_P_S_opt, ED_D_S_opt, ED_L_S_opt)

		C_V_RD_N = DataFrame(time=Int64[], node=Int64[], CostPP=Float64[], CostRes=Float64[], CostStor=Float64[], CostLostLoad=Float64[], CostLostGen=Float64[], CostTotal=Float64[],
			VolumePP=Float64[], VolumeRes=Float64[], VolumeStor=Float64[], VolumeLostLoad=Float64[], VolumeLostGen=Float64[], VolumeTotal=Float64[])
	    for t in hours
			for n in N
				if isempty(Node2PP[n])
					costpp = 0
				else
					costpp = sum((pp.mc[g][1]*CM_ΔP_up_opt[t, g] + (ED_price[t] - pp.mc[g][]) * CM_ΔP_dn_opt[t, g]) for g in Node2PP[n])
				end
				if isempty(Node2Res[n])
					costres = 0
				else
					costres = sum(((ED_price[t]) * CM_ΔP_R_dn_opt[t, res]) for res in Node2Res[n])
				end
				if isempty(Node2Stor[n])
					coststor = 0
				else
					coststor = sum(((ED_price[t]/storages.efficiency[stor]) * CM_P_S_up_opt[t, stor]) for stor in Node2Stor[n])
				end

				costlostload = 1000*CM_P_load_lost_opt[t, n]
				costlostgen = 1000*CM_P_gen_lost_opt[t, n] # ÜBERALL NULL!
				totalcost = costpp + costres + coststor + costlostload + costlostgen

				if isempty(Node2PP[n])
					volumepp = 0
				else
					volumepp = sum((CM_ΔP_up_opt[t, g] - CM_ΔP_dn_opt[t, g]) for g in Node2PP[n])
	    #            VolumePP = sum((pp.mc[g][1]*CM_ΔP_up_opt[t, g] + (ED_price[t] - pp.mc[g][]) * CM_ΔP_dn_opt[t, g]) for g in Node2PP[n])
				end
				if isempty(Node2Res[n])
					volumeres = 0
				else
					volumeres = sum(-(CM_ΔP_R_dn_opt[t, res]) for res in Node2Res[n])
				end
				if isempty(Node2Stor[n])
					volumestor = 0
				else
					volumestor = sum((CM_P_S_up_opt[t, stor]) for stor in Node2Stor[n])
				end

				volumelostload = CM_P_load_lost_opt[t, n]
				volumelostgen = CM_P_gen_lost_opt[t, n]
				totalvolume = volumepp + volumeres + volumestor + volumelostload + volumelostgen

				push!(C_V_RD_N, [t n costpp costres coststor costlostload costlostgen totalcost volumepp volumeres volumestor volumelostload volumelostgen totalvolume])
			end
		end
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "C_V_RD_N__"*timeperiod_str*".csv"), C_V_RD_N, append = true)

			# Dispatchable power plants
			output_CM_ΔP_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
			try
				for t in hours
					for g in pp.unit
						val = CM_ΔP_up_opt[t, g]
						fuel =  pp.fuel[g]
						push!(output_CM_ΔP_up_opt, [t g fuel val pp.node[g]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_up.csv"), output_CM_ΔP_up_opt, append = true)
			catch e
				println("CM_ΔP_up_opt does not exist.")
			end

			output_CM_ΔP_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
			try
				for t in hours
					for g in pp.unit
						val = -1*CM_ΔP_dn_opt[t, g]
						fuel =  pp.fuel[g]
						push!(output_CM_ΔP_dn_opt, [t g fuel val pp.node[g]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_dn.csv"), output_CM_ΔP_dn_opt, append = true)
			catch e
				println("CM_ΔP_dn_opt does not exist.")
			end

			# RES generation
			output_CM_ΔP_R_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
			try
				for t in hours
					for r in res.unit
						val = CM_ΔP_R_up_opt[t, r]
						fuel =  res.fuel[r]
						push!(output_CM_ΔP_R_up_opt, [t r fuel val res.node[r]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_R_up.csv"), output_CM_ΔP_R_up_opt, append = true)
			catch e
				println("CM_ΔP_R_up_opt does not exist.")
			end

			output_CM_ΔP_R_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
			try
				for t in hours
					for r in res.unit
						val = -1*CM_ΔP_R_dn_opt[t, r]
						fuel =  res.fuel[r]
						push!(output_CM_ΔP_R_dn_opt, [t r fuel val res.node[r]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_R_dn.csv"), output_CM_ΔP_R_dn_opt, append = true)
			catch e
				println("CM_ΔP_R_dn_opt does not exist.")
			end

			### Model parameters
			output_CM_info = DataFrame(timeperiod=String[],
											   termination_status=String[],
											   primal_status=String[],
											   dual_status=String[],
											   objective_value=Float64[],
											   solve_time=Float64[])

			push!(output_CM_info, [string(hours) string(termination_status(CM_mod)) string(primal_status(CM_mod)) string(dual_status(CM_mod)) objective_value(CM_mod) solve_time(CM_mod)])
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_info.csv"), output_CM_info, append = true)

			### Power flow results
			output_CM_P_flow_opt = DataFrame(time=Int64[],
											 line=Int64[],
											 from=Int64[],
											 to=Int64[],
											 pmax=Float64[],
											 pflow=Float64[])

			try
			   for t in hours
				   for l in lines.id
					   from = parse(Int, string(lines.from[l]))
					   to = parse(Int, string(lines.to[l]))
					   pmax = lines.pmax[l]
					   pflow = CM_P_flow_opt[t, l]
					   push!(output_CM_P_flow_opt, [t l from to pmax pflow])
				   end
			   end
			   CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_flow.csv"), output_CM_P_flow_opt, append = true)
		   catch e
			   println("CM_P_flow_opt does not exist.")
		   end

		   # Lost generation
		   output_CM_P_gen_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
		   try
			   for t in hours
				   for j in nodes.id
					   output = CM_P_gen_lost_opt[t, j]
					   push!(output_CM_P_gen_lost, [t, j, output])
				   end
			   end
			   CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_gen_lost.csv"), output_CM_P_gen_lost, append = true)
		   catch e
			   println("CM_P_gen_lost does not exist.")
		   end

		   # Lost load
		   output_CM_P_load_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
		   try
			   for t in hours
				   for j in nodes.id
					   output = CM_P_load_lost_opt[t, j]
					   push!(output_CM_P_load_lost, [t, j, output])
				   end
			   end
			   CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_load_lost.csv"), output_CM_P_load_lost, append = true)
		   catch e
			   println("CM_P_load_lost does not exist.")
		   end

		   # Noral Price
		   output_CM_price = DataFrame(time=Int64[],node=Int64[],price=Float64[])
	   	    try
	   	    	for t in hours
					for j in nodes.id
		   	    		price =  CM_price[j, t]
		   	    		push!(output_CM_price, [t j price])
					end
	   	    	end
	   	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_price.csv"), output_CM_price, append = true)
	   	    catch e
	   	    	println("CM_price does not exist.")
	   	    end

			output_CM_P_S_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
			try
				for t in hours
					for s in storages.unit
						val = CM_P_S_up_opt[t, s]
						fuel =  storages.fuel[s]
						push!(output_CM_P_S_up_opt, [t s fuel val storages.node[s]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_S_up.csv"), output_CM_P_S_up_opt, append = true)
			catch e
				println("CM_P_S_up_opt does not exist.")
			end

			# Nodal exchange
			output_exchange = DataFrame(time=Int64[],node=Int64[],output=Float64[])
			try
				for t in hours
					for j in nodes.id
						push!(output_exchange, [t j nodes.exchange[j][t]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "exchange.csv"), output_exchange, append = true)
			catch e
				println("output_exchange does not exist.")
			end

	next!(prog)

end

println()
println("Calculation complete.")


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

precompile(Tuple{typeof(Base.empty!), Base.Dict{Int64, Union{Distributed.Worker, Distributed.LocalProcess}}})
precompile(Tuple{typeof(Core.Compiler.isbits), Distributed.DefaultClusterManager})
precompile(Tuple{typeof(Distributed.init_worker), String, Distributed.DefaultClusterManager})
precompile(Tuple{typeof(Distributed.local_remotecall_thunk), typeof(Distributed.set_valid_processes), Tuple{Array{Int64, 1}}, Array{Any, 1}})
precompile(Tuple{typeof(Distributed.remote_do), typeof(Distributed.set_valid_processes), Distributed.Worker, Array{Int64, 1}})
precompile(Tuple{typeof(Distributed.remote_do), typeof(Distributed.set_valid_processes), Distributed.LocalProcess, Array{Int64, 1}})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remote_do")), Array{Any, 1}, typeof(Distributed.remote_do), typeof(Distributed.set_valid_processes), Distributed.Worker, Array{Int64, 1}})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remote_do")), Array{Any, 1}, typeof(Distributed.remote_do), typeof(Distributed.set_valid_processes), Distributed.LocalProcess, Array{Int64, 1}})
precompile(Tuple{typeof(Distributed.default_addprocs_params)})
precompile(Tuple{typeof(Distributed.topology), Symbol})
precompile(Tuple{typeof(Base.popfirst!), Array{Distributed.WorkerConfig, 1}})
precompile(Tuple{typeof(Distributed.workers)})
precompile(Tuple{getfield(Distributed, Symbol("#kw#addprocs_locked")), Array{Any, 1}, typeof(Distributed.addprocs_locked), Distributed.SSHManager})
precompile(Tuple{typeof(Distributed.check_addprocs_args), Array{Any, 1}})
precompile(Tuple{Type{Distributed.SSHManager}, Array{Any, 1}})
precompile(Tuple{getfield(Distributed, Symbol("#kw#addprocs")), Array{Any, 1}, typeof(Distributed.addprocs), Distributed.SSHManager})
precompile(Tuple{getfield(Distributed, Symbol("#kw#addprocs_locked")), Array{Any, 1}, typeof(Distributed.addprocs_locked), Distributed.LocalManager})
precompile(Tuple{getfield(Distributed, Symbol("#kw#addprocs")), Array{Any, 1}, typeof(Distributed.addprocs), Distributed.LocalManager})
precompile(Tuple{typeof(Distributed.check_master_connect)})
precompile(Tuple{typeof(Distributed.terminate_all_workers)})
precompile(Tuple{typeof(Distributed.local_remotecall_thunk), typeof(Base.exit), Tuple{}, Array{Any, 1}})
precompile(Tuple{typeof(Distributed.remote_do), typeof(Base.exit), Distributed.Worker})
precompile(Tuple{typeof(Distributed.remote_do), typeof(Base.exit), Distributed.LocalProcess})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remote_do")), Array{Any, 1}, typeof(Distributed.remote_do), typeof(Base.exit), Distributed.Worker})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remote_do")), Array{Any, 1}, typeof(Distributed.remote_do), typeof(Base.exit), Distributed.LocalProcess})
precompile(Tuple{typeof(Distributed.set_worker_state), Distributed.Worker, Distributed.WorkerState})
precompile(Tuple{typeof(Distributed.set_worker_state), Distributed.LocalProcess, Distributed.WorkerState})
precompile(Tuple{getfield(Distributed, Symbol("#kw#rmprocs")), Array{Any, 1}, typeof(Distributed.rmprocs), Array{Int64, 1}})
precompile(Tuple{typeof(Distributed.interrupt), Array{Int64, 1}})
precompile(Tuple{typeof(Distributed.flush_gc_msgs)})
precompile(Tuple{typeof(Distributed.addprocs), Int64})
precompile(Tuple{Type{Distributed.WorkerConfig}})
precompile(Tuple{typeof(Distributed.launch), Distributed.LocalManager, Base.Dict{Any, Any}, Array{Distributed.WorkerConfig, 1}, Base.Condition})
precompile(Tuple{typeof(Distributed.start_worker), Base.PipeEndpoint, String})
precompile(Tuple{typeof(Distributed.socket_reuse_port)})
precompile(Tuple{typeof(Distributed.flush_gc_msgs)})
precompile(Tuple{typeof(Distributed.disable_nagle), Sockets.TCPServer})
precompile(Tuple{typeof(Distributed.next_tunnel_port)})
precompile(Tuple{typeof(Base._delete!), Base.Dict{Int64, Union{Distributed.Worker, Distributed.LocalProcess}}, Int64})
precompile(Tuple{typeof(Distributed.send_msg_), Distributed.Worker, Distributed.MsgHeader, Distributed.JoinPGRPMsg, Bool})
precompile(Tuple{typeof(Distributed.send_msg_now), Distributed.Worker, Distributed.MsgHeader, Distributed.JoinPGRPMsg})
precompile(Tuple{typeof(Distributed.connect_w2w), Int64, Distributed.WorkerConfig})
precompile(Tuple{typeof(Distributed.create_worker), Distributed.LocalManager, Distributed.WorkerConfig})
precompile(Tuple{typeof(Distributed.setup_launched_worker), Distributed.LocalManager, Distributed.WorkerConfig, Array{Int64, 1}})
precompile(Tuple{typeof(Sockets.connect), Distributed.LocalManager, Int64, Distributed.WorkerConfig})
precompile(Tuple{typeof(Distributed.read_worker_host_port), Base.Pipe})
precompile(Tuple{typeof(Distributed.parse_connection_info), String})
precompile(Tuple{typeof(Distributed.connect_to_worker), Base.SubString{String}, Int16})
precompile(Tuple{typeof(Distributed.process_messages), Sockets.TCPSocket, Sockets.TCPSocket, Bool})
precompile(Tuple{getfield(Core, Symbol("#kw#Type")), Array{Any, 1}, Type{Distributed.Worker}, Int64, Sockets.TCPSocket, Sockets.TCPSocket, Distributed.LocalManager})
precompile(Tuple{typeof(Distributed.worker_id_from_socket), Sockets.TCPSocket})
precompile(Tuple{Type{Distributed.ClusterSerializer{Sockets.TCPSocket}}, Sockets.TCPSocket})
precompile(Tuple{typeof(Distributed.send_msg_), Distributed.Worker, Distributed.MsgHeader, Distributed.ResultMsg, Bool})
precompile(Tuple{typeof(Distributed.send_msg_now), Distributed.Worker, Distributed.MsgHeader, Distributed.ResultMsg})
precompile(Tuple{typeof(Distributed.register_worker_streams), Distributed.Worker})
precompile(Tuple{typeof(Distributed.register_worker_streams), Distributed.LocalProcess})
precompile(Tuple{Type{Distributed.ClusterSerializer{Sockets.TCPSocket}}, Sockets.TCPSocket})
precompile(Tuple{typeof(Distributed.worker_id_from_socket), Sockets.TCPSocket})
precompile(Tuple{typeof(Base.convert), Type{Distributed.ClusterSerializer{I} where I<:IO}, Distributed.ClusterSerializer{Sockets.TCPSocket}})
precompile(Tuple{typeof(Base.convert), Type{Distributed.ClusterManager}, Distributed.LocalManager})
precompile(Tuple{typeof(Base.convert), Type{Distributed.WorkerConfig}, Distributed.WorkerConfig})
precompile(Tuple{typeof(Base.get), Base.Dict{Any, Any}, Distributed.RRID, Bool})
precompile(Tuple{typeof(Base.ht_keyindex), Base.Dict{Any, Any}, Distributed.RRID})
precompile(Tuple{typeof(Distributed.local_remotecall_thunk), typeof(Distributed.rmprocs), Tuple{Int64}, Array{Any, 1}})
precompile(Tuple{typeof(Distributed.remote_do), typeof(Distributed.rmprocs), Distributed.Worker, Int64})
precompile(Tuple{typeof(Distributed.remote_do), typeof(Distributed.rmprocs), Distributed.LocalProcess, Int64})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remote_do")), Array{Any, 1}, typeof(Distributed.remote_do), typeof(Distributed.rmprocs), Distributed.Worker, Int64})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remote_do")), Array{Any, 1}, typeof(Distributed.remote_do), typeof(Distributed.rmprocs), Distributed.LocalProcess, Int64})
precompile(Tuple{Type{Distributed.ResultMsg}, Distributed.RemoteException})
precompile(Tuple{Type{Distributed.ResultMsg}, Symbol})
precompile(Tuple{typeof(Distributed.send_msg_now), Sockets.TCPSocket, Distributed.MsgHeader, Distributed.ResultMsg})
precompile(Tuple{typeof(Base._delete!), Base.Dict{Int64, Union{Distributed.Worker, Distributed.LocalProcess}}, Int64})
precompile(Tuple{typeof(Distributed.def_rv_channel)})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Any, Any}, Distributed.RemoteValue, Distributed.RRID})
precompile(Tuple{typeof(Base.ht_keyindex2!), Base.Dict{Any, Any}, Distributed.RRID})
precompile(Tuple{typeof(Base._setindex!), Base.Dict{Any, Any}, Distributed.RemoteValue, Distributed.RRID, Int64})
precompile(Tuple{typeof(Base.notify), Base.Condition, Distributed.ProcessExitedException, Bool, Bool})
precompile(Tuple{typeof(Distributed.process_messages), Sockets.TCPSocket, Sockets.TCPSocket, Bool})
precompile(Tuple{typeof(Base.pop!), Base.Dict{Int64, Union{Distributed.Worker, Distributed.LocalProcess}}, Int64, Nothing})
precompile(Tuple{typeof(Distributed.send_connection_hdr), Distributed.Worker, Bool})
precompile(Tuple{typeof(Distributed.deregister_worker), Distributed.ProcessGroup, Int64})
precompile(Tuple{typeof(Distributed.process_hdr), Sockets.TCPSocket, Bool})
precompile(Tuple{typeof(Distributed.deserialize_msg), Distributed.ClusterSerializer{Sockets.TCPSocket}})
precompile(Tuple{typeof(Distributed.null_id), Distributed.RRID})
precompile(Tuple{typeof(Distributed.deliver_result), Sockets.TCPSocket, Symbol, Distributed.RRID, Distributed.RemoteException})
precompile(Tuple{typeof(Distributed.disable_nagle), Sockets.TCPSocket})
precompile(Tuple{typeof(Distributed.message_handler_loop), Sockets.TCPSocket, Sockets.TCPSocket, Bool})
precompile(Tuple{typeof(Distributed.process_tcp_streams), Sockets.TCPSocket, Sockets.TCPSocket, Bool})
precompile(Tuple{Type{Distributed.JoinPGRPMsg}, Int64, Array{Union{Tuple{Any, Int64}, Tuple{Tuple{}, Any, Bool}}, 1}, Symbol, Bool})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Distributed.JoinPGRPMsg})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Int64})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Array{Union{Tuple{Any, Int64}, Tuple{Tuple{}, Any, Bool}}, 1}})
precompile(Tuple{typeof(Serialization.should_send_whole_type), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{Any}})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Core.SimpleVector})
precompile(Tuple{typeof(Serialization.serialize_type_data), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{Any}, Bool})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{Any}})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Tuple{Int64}})
precompile(Tuple{typeof(Serialization.serialize_cycle), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{Union{Tuple{Any, Int64}, Tuple{Tuple{}, Any, Bool}}}})
precompile(Tuple{typeof(Serialization.serialize_type), Distributed.ClusterSerializer{Sockets.TCPSocket}, DataType})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Bool})
precompile(Tuple{typeof(Distributed.serialize_global_from_main), Distributed.ClusterSerializer{Sockets.TCPSocket}, Symbol})
precompile(Tuple{typeof(Serialization.serialize_mod_names), Distributed.ClusterSerializer{Sockets.TCPSocket}, Module})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Symbol})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Module})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Array{Any, 1}})
precompile(Tuple{typeof(Serialization.serialize_cycle), Distributed.ClusterSerializer{Sockets.TCPSocket}, Core.TypeName})
precompile(Tuple{typeof(Serialization.serialize_typename), Distributed.ClusterSerializer{Sockets.TCPSocket}, Core.TypeName})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Array{Symbol, 1}})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Core.TypeName})
precompile(Tuple{typeof(Serialization.serialize_cycle_header), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{Union{Tuple{Any, Int64}, Tuple{Tuple{}, Any, Bool}}}})
precompile(Tuple{typeof(Serialization.serialize_any), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{Union{Tuple{Any, Int64}, Tuple{Tuple{}, Any, Bool}}}})
precompile(Tuple{typeof(Distributed.send_msg_), Distributed.Worker, Distributed.MsgHeader, Distributed.ResultMsg, Bool})
precompile(Tuple{typeof(Distributed.send_msg_now), Distributed.Worker, Distributed.MsgHeader, Distributed.ResultMsg})
precompile(Tuple{typeof(Distributed.local_remotecall_thunk), typeof(Distributed.rmprocs), Tuple{Int64}, Array{Any, 1}})
precompile(Tuple{typeof(Distributed.remote_do), typeof(Distributed.rmprocs), Distributed.Worker, Int64})
precompile(Tuple{typeof(Distributed.remote_do), typeof(Distributed.rmprocs), Distributed.LocalProcess, Int64})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remote_do")), Array{Any, 1}, typeof(Distributed.remote_do), typeof(Distributed.rmprocs), Distributed.Worker, Int64})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remote_do")), Array{Any, 1}, typeof(Distributed.remote_do), typeof(Distributed.rmprocs), Distributed.LocalProcess, Int64})
precompile(Tuple{Type{Distributed.ResultMsg}, Distributed.RemoteException})
precompile(Tuple{Type{Distributed.ResultMsg}, Symbol})
precompile(Tuple{typeof(Distributed.send_msg_now), Sockets.TCPSocket, Distributed.MsgHeader, Distributed.ResultMsg})
precompile(Tuple{typeof(Base.pop!), Base.Dict{Int64, Union{Distributed.Worker, Distributed.LocalProcess}}, Int64, Nothing})
precompile(Tuple{typeof(Distributed.deregister_worker), Distributed.ProcessGroup, Int64})
precompile(Tuple{typeof(Distributed.process_hdr), Sockets.TCPSocket, Bool})
precompile(Tuple{typeof(Distributed.null_id), Distributed.RRID})
precompile(Tuple{typeof(Distributed.deliver_result), Sockets.TCPSocket, Symbol, Distributed.RRID, Distributed.RemoteException})
precompile(Tuple{typeof(Distributed.disable_nagle), Sockets.TCPSocket})
precompile(Tuple{typeof(Distributed.message_handler_loop), Sockets.TCPSocket, Sockets.TCPSocket, Bool})
precompile(Tuple{typeof(Distributed.process_tcp_streams), Sockets.TCPSocket, Sockets.TCPSocket, Bool})
precompile(Tuple{typeof(Serialization.deserialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{Union}})
precompile(Tuple{typeof(Serialization.deserialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{Module}})
precompile(Tuple{typeof(Serialization.deserialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{Core.SimpleVector}})
precompile(Tuple{Type{Distributed.JoinPGRPMsg}, Int64, Array{Union{Tuple{Any, Int64}, Tuple{Tuple{}, Any, Bool}}, 1}, Symbol, Bool})
precompile(Tuple{typeof(Distributed.handle_msg), Distributed.JoinPGRPMsg, Distributed.MsgHeader, Sockets.TCPSocket, Sockets.TCPSocket, Base.VersionNumber})
precompile(Tuple{Type{Distributed.WorkerConfig}})
precompile(Tuple{typeof(Distributed.send_msg_), Distributed.Worker, Distributed.MsgHeader, Distributed.JoinCompleteMsg, Bool})
precompile(Tuple{typeof(Distributed.send_msg_now), Distributed.Worker, Distributed.MsgHeader, Distributed.JoinCompleteMsg})
precompile(Tuple{getfield(Core, Symbol("#kw#Type")), Array{Any, 1}, Type{Distributed.Worker}, Int64, Sockets.TCPSocket, Sockets.TCPSocket, Distributed.DefaultClusterManager})
precompile(Tuple{typeof(Distributed.register_worker_streams), Distributed.Worker})
precompile(Tuple{typeof(Distributed.register_worker_streams), Distributed.LocalProcess})
precompile(Tuple{typeof(Base.convert), Type{Distributed.ClusterSerializer{I} where I<:IO}, Distributed.ClusterSerializer{Sockets.TCPSocket}})
precompile(Tuple{typeof(Base.convert), Type{Distributed.ClusterManager}, Distributed.DefaultClusterManager})
precompile(Tuple{typeof(Base.convert), Type{Distributed.WorkerConfig}, Distributed.WorkerConfig})
precompile(Tuple{typeof(Distributed.send_connection_hdr), Distributed.Worker, Bool})
precompile(Tuple{typeof(Distributed.manage), Distributed.LocalManager, Int64, Distributed.WorkerConfig, Symbol})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Distributed.JoinCompleteMsg})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Int64})
precompile(Tuple{typeof(Serialization.deserialize), Distributed.ClusterSerializer{Sockets.TCPSocket}})
precompile(Tuple{typeof(Serialization.deserialize_cycle), Distributed.ClusterSerializer{Sockets.TCPSocket}, Expr})
precompile(Tuple{typeof(Serialization.handle_deserialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Int32})
precompile(Tuple{typeof(Serialization.deserialize_array), Distributed.ClusterSerializer{Sockets.TCPSocket}})
precompile(Tuple{typeof(Serialization.deserialize_datatype), Distributed.ClusterSerializer{Sockets.TCPSocket}})
precompile(Tuple{typeof(Serialization.deserialize_expr), Distributed.ClusterSerializer{Sockets.TCPSocket}, Int64})
precompile(Tuple{typeof(Serialization.deserialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{Int64}})
precompile(Tuple{Type{Distributed.JoinCompleteMsg}, Int64, Int64})
precompile(Tuple{typeof(Distributed.handle_msg), Distributed.JoinCompleteMsg, Distributed.MsgHeader, Sockets.TCPSocket, Sockets.TCPSocket, Base.VersionNumber})
precompile(Tuple{typeof(Base.hash), Distributed.RemoteChannel{Base.Channel{Any}}, UInt64})
precompile(Tuple{typeof(Base.ht_keyindex), Base.Dict{WeakRef, Nothing}, Distributed.RemoteChannel{Base.Channel{Any}}})
precompile(Tuple{typeof(Distributed.remotecall_fetch), typeof(Distributed.put_ref), Distributed.Worker, Distributed.RRID, Distributed.WorkerPool})
precompile(Tuple{typeof(Distributed.remotecall_fetch), typeof(Distributed.put_ref), Distributed.LocalProcess, Distributed.RRID, Distributed.WorkerPool})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remotecall_fetch")), Array{Any, 1}, typeof(Distributed.remotecall_fetch), typeof(Distributed.put_ref), Distributed.Worker, Distributed.RRID, Distributed.WorkerPool})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remotecall_fetch")), Array{Any, 1}, typeof(Distributed.remotecall_fetch), typeof(Distributed.put_ref), Distributed.LocalProcess, Distributed.RRID, Distributed.WorkerPool})
precompile(Tuple{typeof(Base.finalizer), Distributed.RemoteChannel{Base.Channel{Any}}, typeof(Distributed.finalize_ref)})
precompile(Tuple{typeof(Distributed.test_existing_ref), Distributed.RemoteChannel{Base.Channel{Any}}})
precompile(Tuple{Type{Distributed.RemoteChannel{T} where T<:Base.AbstractChannel}, Int64})
precompile(Tuple{Type{Distributed.WorkerPool}})
precompile(Tuple{typeof(Distributed.default_worker_pool)})
precompile(Tuple{typeof(Distributed.call_on_owner), typeof(Distributed.put_ref), Distributed.RemoteChannel{Base.Channel{Any}}, Distributed.WorkerPool})
precompile(Tuple{typeof(Distributed.put_ref), Distributed.RRID, Distributed.WorkerPool})
precompile(Tuple{typeof(Base.put!), Distributed.RemoteValue, Distributed.WorkerPool})
precompile(Tuple{typeof(Base.put!), Base.Channel{Any}, Distributed.WorkerPool})
precompile(Tuple{typeof(Base.put_buffered), Base.Channel{Any}, Distributed.WorkerPool})
precompile(Tuple{typeof(Base.put_unbuffered), Base.Channel{Any}, Distributed.WorkerPool})
precompile(Tuple{typeof(Base.push!), Distributed.WorkerPool, Int64})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Distributed.RemoteDoMsg})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, typeof(Distributed.set_valid_processes)})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Tuple{Array{Int64, 1}}})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Array{Int64, 1}})
precompile(Tuple{typeof(Serialization.deserialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{typeof(Distributed.set_valid_processes)}})
precompile(Tuple{typeof(Distributed.handle_msg), Distributed.RemoteDoMsg, Distributed.MsgHeader, Sockets.TCPSocket, Sockets.TCPSocket, Base.VersionNumber})
precompile(Tuple{typeof(Distributed.set_valid_processes), Array{Int64, 1}})
precompile(Tuple{typeof(Distributed._rmprocs), Array{Int64, 1}, Float64})
precompile(Tuple{typeof(Base.kill), Distributed.LocalManager, Int64, Distributed.WorkerConfig})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, typeof(Base.exit)})
precompile(Tuple{typeof(Serialization.serialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Tuple{}})
precompile(Tuple{typeof(Serialization.deserialize), Distributed.ClusterSerializer{Sockets.TCPSocket}, Type{typeof(Base.exit)}})
precompile(Tuple{typeof(Distributed.finalize_ref), Distributed.RemoteChannel{Base.Channel{Any}}})
precompile(Tuple{typeof(Distributed.send_del_client), Distributed.RemoteChannel{Base.Channel{Any}}})
precompile(Tuple{typeof(Base.:(==)), Distributed.RemoteChannel{Base.Channel{Any}}, Distributed.RemoteChannel{Base.Channel{Any}}})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remotecall_fetch")), Array{Any, 1}, typeof(Distributed.remotecall_fetch), typeof(Base.open), Distributed.LocalProcess, typeof(Base.read), String})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remotecall_fetch")), Array{Any, 1}, typeof(Distributed.remotecall_fetch), typeof(Base.open), Distributed.Worker, typeof(Base.read), String})
precompile(Tuple{getfield(Distributed, Symbol("#kw#remote_do")), Array{Any, 1}, typeof(Distributed.remote_do), typeof(Base.exit), Distributed.Worker})
precompile(Tuple{getfield(Distributed, Symbol("#kw#rmprocs")), Array{Any, 1}, typeof(Distributed.rmprocs), Array{Int64, 1}})
precompile(Tuple{Type{Distributed.Future}, Int64})
precompile(Tuple{typeof(Distributed.flush_gc_msgs), Distributed.Worker})
precompile(Tuple{typeof(Distributed.remote_do), typeof(Base.exit), Distributed.Worker})
precompile(Tuple{typeof(Distributed.send_del_client), Distributed.Future})
precompile(Tuple{typeof(Distributed.send_msg), Distributed.Worker, Distributed.MsgHeader, Distributed.CallMsg{:call}})
precompile(Tuple{typeof(Distributed.send_msg_), Distributed.Worker, Distributed.MsgHeader, Distributed.CallMsg{:call}, Bool})
precompile(Tuple{typeof(Distributed.send_msg), Distributed.Worker, Distributed.MsgHeader, Distributed.CallMsg{:call_fetch}})
precompile(Tuple{typeof(Distributed.send_msg_), Distributed.Worker, Distributed.MsgHeader, Distributed.CallMsg{:call_fetch}, Bool})
precompile(Tuple{typeof(Distributed.send_msg), Distributed.Worker, Distributed.MsgHeader, Distributed.RemoteDoMsg})
precompile(Tuple{typeof(Distributed.send_msg_), Distributed.Worker, Distributed.MsgHeader, Distributed.RemoteDoMsg, Bool})
precompile(Tuple{typeof(Distributed.terminate_all_workers)})
precompile(Tuple{typeof(Distributed.test_existing_ref), Distributed.Future})
precompile(Tuple{typeof(Base.finalizer), Distributed.Future, typeof(Distributed.finalize_ref)})
precompile(Tuple{typeof(Base.hash), Distributed.Future, UInt64})
precompile(Tuple{typeof(Base.ht_keyindex), Base.Dict{WeakRef, Nothing}, Distributed.Future})
precompile(Tuple{Type{Union{}}, Distributed.RRID})

max_stars_repo_name: musm/julia

max_issues_repo_name: musm/julia

max_forks_repo_name: musm/julia


text:
function inline_consts!(ir::IR)
  env = Dict()
  for (v, st) in ir
    if st.type isa Union{Const,Node}
      delete!(ir, v)
      env[v] = st.type.value
    else
      env[v] = v
    end
  end
  return IRTools.prewalk!(x -> get(env, x, x), ir)
end

function partials!(ir::IR)
  slots = IdDict()
  slot(k, T) = Base.@get!(slots, k, Slot(gensym(:partial), T))
  for (v, st) in ir
    ex = st.expr
    if iscall(ex, setindex!)
      exprtype(ir, ex.args[2]) isa Partial{<:Dict} || continue
      key = exprtype(ir, ex.args[2]), exprtype(ir, ex.args[4])::Const
      T = key[1].value[key[2].value]
      insert!(ir, v, Expr(:(=), slot(key, T), ex.args[3]))
    elseif iscall(ex, getindex)
      exprtype(ir, ex.args[2]) isa Partial{<:Dict} || continue
      key = exprtype(ir, ex.args[2]), exprtype(ir, ex.args[3])::Const
      T = key[1].value[key[2].value]
      ir[v] = slot(key, T)
    end
  end
  for (v, st) in ir
    st.type isa Partial{<:Dict} && delete!(ir, v)
  end
  return ir
end

function trimblocks!(ir)
  i = 1
  while i <= length(blocks(ir))
    bl = block(ir, i)
    # assume no implicit branches
    # don't bother with blocks that have args (for now)
    if isempty(bl) && length(branches(bl)) == 1 && isempty(arguments(bl))
      for a in predecessors(bl), i in 1:length(branches(a))
        br = branches(bl)[1]
        branches(a)[i].block == bl.id || continue
        branches(a)[i] = Branch(branches(a)[i], block = br.block, args = br.args)
      end
      deleteblock!(ir, i)
    else
      i += 1
    end
  end
  return ir
end

cleanup!(ir) =
  ir |> inline_consts! |> partials! |> ssa! |> prune! |> renumber


text:
using ConvexTests, GLPK

@info "Starting GLPK tests"

do_tests("GLPK", GLPK.Optimizer; exclude = [r"exp", r"sdp", r"socp"])

max_stars_repo_name: ericphanson/ConvexTests

max_issues_repo_name: ericphanson/ConvexTests

max_forks_repo_name: ericphanson/ConvexTests


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Collection of sources required to build GMP
name = "GMP"
version = v"6.2.0"

sources = [
    ArchiveSource("https://gmplib.org/download/gmp/gmp-$(version).tar.bz2",
                  "f51c99cb114deb21a60075ffb494c1a210eb9d7cb729ed042ddb7de9534451ea"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gmp-*

# Include Julia-carried patches
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/gmp_alloc_overflow_func.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/gmp-exception.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/gmp-apple-arm64.patch

flags=(--enable-cxx --enable-shared --disable-static)

# On x86_64 architectures, build fat binary
if [[ ${proc_family} == intel ]]; then
    flags+=(--enable-fat)
fi

autoreconf
./configure --prefix=$prefix --build=${MACHTYPE} --host=${target} ${flags[@]}

make -j${nproc}
make install

# On Windows, we need to make sure that the non-versioned dll names exist too
if [[ ${target} == *mingw* ]]; then
    cp -v ${libdir}/libgmp-*.dll "${libdir}/libgmp.dll"
    cp -v ${libdir}/libgmpxx-*.dll "${libdir}/libgmpxx.dll"
fi

# GMP is dual-licensed, install all license files
install_license COPYING*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis([supported_platforms(); Platform("aarch64", "macos")])

# The products that we will ensure are always built
products = [
    LibraryProduct("libgmp", :libgmp),
    LibraryProduct("libgmpxx", :libgmpxx),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6", julia_compat="1.6")

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil