text:
################################################################################

@testset verbose = true "Write" begin

  ################################################################################

  # testing
  @testset "Write Hidden Markov model" begin

    ################################################################################

    # empty dict
    toWriteDc = Dict{String, Tuple{Array{Int64, 1}, Array{Array{Float64, 1}, 1}}}("Test1" => ([1, 2, 3], [[1., 2., 3.]]))

    # write
    writeHMM("out/hmm_", toWriteDc)

    # testing
    @test isfile("out/hmm_Test1_states.csv") == true
    @test isfile("out/hmm_Test1_traceb.csv") == true

    ################################################################################

    # read
    readStates = readdlm("out/hmm_Test1_states.csv")
    readTraceb = readdlm("out/hmm_Test1_traceb.csv")

    # paired values
    testStates = Array{Any, 2}(undef, 4, 1)
    testStates[:, 1] .= ["Test1"; 1; 2; 3]

    testTraceb = Array{Any, 2}(undef, 4, 1)
    testTraceb[:, 1] .= ["S1"; 1.; 2.; 3.]

    # testing
    @test readStates == testStates
    @test readTraceb == testTraceb

    ################################################################################

  end

  ################################################################################

  # testing
  @testset "Write Performance" begin

    ################################################################################

    # empty dict
    toWriteDc = Dict{String, Array{Float64, 2}}("Test1" => [1. 10.], "Test2" => [2. 20.], "Test3" => [3. 30.])

    # transform
    performanceTb = writePerformance(toWriteDc)

    # testing
    @test performanceTb == [["Electrode", "Test2", "Test1", "Test3"] ["Sensitivity", 2., 1., 3.] ["Specificity", 20., 10., 30.]]

    ################################################################################

    # write
    writePerformance("out/performace.csv", toWriteDc)

    # testing
    @test isfile("out/performace.csv") == true

    ################################################################################

    # read
    readPerformance = readdlm("out/performace.csv", ',')

    # testing
    @test readPerformance == [["Electrode", "Test2", "Test1", "Test3"] ["Sensitivity", 2., 1., 3.] ["Specificity", 20., 10., 30.]]

    ################################################################################

    # read
    readStates = readdlm("out/hmm_Test1_states.csv")
    readTraceb = readdlm("out/hmm_Test1_traceb.csv")

    # paired values
    testStates = Array{Any, 2}(undef, 4, 1)
    testStates[:, 1] .= ["Test1"; 1; 2; 3]

    testTraceb = Array{Any, 2}(undef, 4, 1)
    testTraceb[:, 1] .= ["S1"; 1.; 2.; 3.]

    # testing
    @test readStates == testStates
    @test readTraceb == testTraceb

    ################################################################################

  end

  ################################################################################

end

################################################################################

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
@get_units @with_kw struct SBM{T,N,M}
    # Model time step [s]
    Δt::T | "s"
    # Maximum number of soil layers
    maxlayers::Int | "-"
    # number of cells
    n::Int | "-"
    # Number of soil layers
    nlayers::Vector{Int} | "-"
    # Number of unsaturated soil layers
    n_unsatlayers::Vector{Int} | "-"
    # Fraction of river [-]
    riverfrac::Vector{T} | "-"
    # Saturated water content (porosity) [mm mm⁻¹]
    θₛ::Vector{T} | "mm mm-1"
    # Residual water content [mm mm⁻¹]
    θᵣ::Vector{T} | "mm mm-1"
    # Vertical hydraulic conductivity [mm Δt⁻¹] at soil surface
    kv₀::Vector{T} | "mm Δt-1"
    # Muliplication factor [-] applied to kv_z (vertical flow)
    kvfrac::Vector{SVector{N,T}} | "-"
    # Air entry pressure [cm] of soil (Brooks-Corey)
    hb::Vector{T} | "cm"
    # Soil thickness [mm]
    soilthickness::Vector{T}
    # Thickness of soil layers [mm]
    act_thickl::Vector{SVector{N,T}}
    # Cumulative sum of soil layers [mm], starting at soil surface (0)
    sumlayers::Vector{SVector{M,T}}
    # Infiltration capacity of the compacted areas [mm Δt⁻¹]
    infiltcappath::Vector{T} | "mm Δt-1"
    # Soil infiltration capacity [mm/Δt]
    infiltcapsoil::Vector{T} | "mm Δt-1"
    # Maximum leakage [mm/Δt] from saturated zone
    maxleakage::Vector{T} | "mm Δt-1"
    # Fraction of open water (excluding rivers) [-]
    waterfrac::Vector{T} | "-"
    # Fraction of compacted area  [-]
    pathfrac::Vector{T} | "-"
    # Rooting depth [mm]
    rootingdepth::Vector{T}
    # Controls how roots are linked to water table [-]
    rootdistpar::Vector{T} | "-"
    # Parameter [mm] controlling capillary rise
    capscale::Vector{T}
    # Multiplication factor [-] to correct
    et_reftopot::Vector{T} | "-"
    # Brooks-Corey power coefﬁcient [-] for each soil layer
    c::Vector{SVector{N,T}} | "-"
    # Stemflow [mm]
    stemflow::Vector{T}
    # Throughfall [mm]
    throughfall::Vector{T}
    # A scaling parameter [mm⁻¹] (controls exponential decline of kv₀)
    f::Vector{T} | "mm-1"
    # Amount of water in the unsaturated store, per layer [mm]
    ustorelayerdepth::Vector{SVector{N,T}}
    # Saturated store [mm]
    satwaterdepth::Vector{T}
    # Pseudo-water table depth [mm] (top of the saturated zone)
    zi::Vector{T}
    # Soilwater capacity [mm]
    soilwatercapacity::Vector{T}
    # Canopy storage [mm]
    canopystorage::Vector{T}
    # Maximum canopy storage [mm]
    cmax::Vector{T}
    # Canopy gap fraction [-]
    canopygapfraction::Vector{T} | "-"
    # Gash interception model parameter, ratio of the average evaporation from the
    # wet canopy [mm Δt⁻¹] and the average precipitation intensity [mm Δt⁻¹] on a saturated canopy
    e_r::Vector{T} | "-"
    # Precipitation [mm]
    precipitation::Vector{T}
    # Temperature [ᵒC]
    temperature::Vector{T} | "°C"
    # Potential evapotranspiration [mm]
    potential_evaporation::Vector{T}
    # Potential transpiration, open water, river and soil evaporation (after subtracting interception from potential_evaporation)
    pottrans_soil::Vector{T}
    # Transpiration [mm]
    transpiration::Vector{T}
    # Actual evaporation from unsaturated store [mm]
    ae_ustore::Vector{T}
    # Interception [mm]
    interception::Vector{T}
    # Soil evaporation from unsaturated store [mm]
    soilevap::Vector{T}
    # Soil evaporation from saturated store [mm]
    soilevapsat::Vector{T}
    # Actual capillary rise [mm]
    actcapflux::Vector{T}
    # Actual transpiration from saturated store [mm]
    actevapsat::Vector{T}
    # Total actual evapotranspiration [mm]
    actevap::Vector{T}
    # Runoff from river based on riverfrac [mm]
    runoff_river::Vector{T}
    # Runoff from land based on waterfrac [mm]
    runoff_land::Vector{T}
    # Actual evaporation from open water (land) [mm]
    ae_openw_l::Vector{T}
    # Actual evaporation from river [mm]
    ae_openw_r::Vector{T}
    # Net runoff from river [mm]
    net_runoff_river::Vector{T}
    # Water available for infiltration [mm]
    avail_forinfilt::Vector{T}
    # Actual infiltration into the unsaturated zone [mm]
    actinfilt::Vector{T}
    # Actual infiltration non-compacted fraction [mm]
    actinfiltsoil::Vector{T}
    # Actual infiltration compacted fraction [mm]
    actinfiltpath::Vector{T}
    # Infiltration capacity of the non-compacted areas [mm Δt⁻¹]
    infiltsoilpath::Vector{T} | "mm Δt-1"
    # Infiltration excess water [mm]
    infiltexcess::Vector{T}
    # Water that cannot infiltrate due to saturated soil (saturation excess) [mm]
    excesswater::Vector{T}
    # Water exfiltrating during saturation excess conditions [mm]
    exfiltsatwater::Vector{T}
    # Water exfiltrating from unsaturated store because of change in water table [mm]
    exfiltustore::Vector{T}
    # Excess water for non-compacted fraction [mm]
    excesswatersoil::Vector{T}
    # Excess water for compacted fraction [mm]
    excesswaterpath::Vector{T}
    # Total surface runoff from infiltration and saturation excess [mm]
    runoff::Vector{T}
    # Volumetric water content [mm mm⁻¹] per soil layer (including θᵣ and saturated zone)
    vwc::Vector{SVector{N,T}} | "mm mm-1"
    # Volumetric water content [%] per soil layer (including θᵣ and saturated zone)
    vwc_perc::Vector{SVector{N,T}} | "%"
    # Root water storage [mm] in unsaturated and saturated zone (excluding θᵣ)
    rootstore::Vector{T}
    # Volumetric water content [mm mm⁻¹] in root zone (including θᵣ and saturated zone)
    vwc_root::Vector{T} | "mm mm-1"
    # Volumetric water content [%] in root zone (including θᵣ and saturated zone)
    vwc_percroot::Vector{T} | "%"
    # Amount of available water in the unsaturated zone [mm]
    ustoredepth::Vector{T}
    # Downward flux from unsaturated to saturated zone [mm]
    transfer::Vector{T}
    # Net recharge to saturated store [mm]
    recharge::Vector{T}
    # Actual leakage from saturated store [mm]
    actleakage::Vector{T}
    ### Snow parameters ###
    # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹]
    cfmax::Vector{T} | "mm ᵒC-1 Δt-1"
    # Threshold temperature for snowfall [ᵒC]
    tt::Vector{T} | "ᵒC"
    # Threshold temperature interval length [ᵒC]
    tti::Vector{T} | "ᵒC"
    # Threshold temperature for snowmelt [ᵒC]
    ttm::Vector{T} | "ᵒC"
    # Water holding capacity as fraction of current snow pack [-]
    whc::Vector{T} | "-"
    # Soil temperature smooth factor [-]
    w_soil::Vector{T} | "-"
    # Controls soil infiltration reduction factor when soil is frozen [-]
    cf_soil::Vector{T} | "-"
    # Snow storage [mm]
    snow::Vector{T}
    # Liquid water content in the snow pack [mm]
    snowwater::Vector{T}
    # Snow melt + precipitation as rainfall [mm]
    rainfallplusmelt::Vector{T}
    # Threshold temperature for snowfall above glacier [ᵒC]
    g_tt::Vector{T} | "ᵒC"
    # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹] for glacier
    g_cfmax::Vector{T} | "mm ᵒC-1 Δt-1"
    # Fraction of the snowpack on top of the glacier converted into ice [-]
    g_sifrac::Vector{T} | "-"
    # Water within the glacier [mm]
    glacierstore::Vector{T}
    # Fraction covered by a glacier [-]
    glacierfrac::Vector{T} | "-"
    # Top soil temperature [ᵒC]
    tsoil::Vector{T} | "ᵒC"
    ## Interception related to leaf_area_index climatology ###
    # Specific leaf storage [mm]
    sl::Vector{T}
    # Storage woody part of vegetation [mm]
    swood::Vector{T}
    # Extinction coefficient [-] (to calculate canopy gap fraction)
    kext::Vector{T} | "-"
    # Leaf area index [m² m⁻²]
    leaf_area_index::Vector{T} | "m2 m-2"
    # Water level land [mm]
    waterlevel_land::Vector{T}
    # Water level river [mm]
    waterlevel_river::Vector{T}

    function SBM{T,N,M}(args...) where {T,N,M}
        equal_size_vectors(args)
        return new(args...)
    end
end

statevars(::SBM; snow = false) =
    snow ? (:satwaterdepth, :snow, :tsoil, :ustorelayerdepth, :snowwater, :canopystorage) :
    (:satwaterdepth, :ustorelayerdepth, :canopystorage)

function initialize_canopy(nc, config, inds)
    n = length(inds)
    # if leaf area index climatology provided use sl, swood and kext to calculate cmax, e_r and canopygapfraction
    if haskey(config.input.vertical, "leaf_area_index")
        # TODO confirm if leaf area index climatology is present in the NetCDF
        sl = ncread(
            nc,
            param(config, "input.vertical.specific_leaf");
            sel = inds,
            type = Float,
        )
        swood = ncread(
            nc,
            param(config, "input.vertical.storage_wood");
            sel = inds,
            type = Float,
        )
        kext = ncread(nc, param(config, "input.vertical.kext"); sel = inds, type = Float)
        cmax = fill(mv, n)
        e_r = fill(mv, n)
        canopygapfraction = fill(mv, n)
    else
        sl = fill(mv, n)
        swood = fill(mv, n)
        kext = fill(mv, n)
        # cmax, e_r, canopygapfraction only required when leaf area index climatology not provided
        cmax = ncread(
            nc,
            param(config, "input.vertical.cmax", nothing);
            sel = inds,
            defaults = 1.0,
            type = Float,
        )
        e_r = ncread(
            nc,
            param(config, "input.vertical.eoverr", nothing);
            sel = inds,
            defaults = 0.1,
            type = Float,
        )
        canopygapfraction = ncread(
            nc,
            param(config, "input.vertical.canopygapfraction", nothing);
            sel = inds,
            defaults = 0.1,
            type = Float,
        )
    end
    return cmax, e_r, canopygapfraction, sl, swood, kext
end

function initialize_sbm(nc, config, riverfrac, inds)

    Δt = Second(config.timestepsecs)
    config_thicknesslayers = get(config.model, "thicknesslayers", Float[])
    if length(config_thicknesslayers) > 0
        thicknesslayers = SVector(Tuple(push!(Float.(config_thicknesslayers), mv)))
        sumlayers = pushfirst(cumsum(thicknesslayers), 0.0)
        maxlayers = length(thicknesslayers) # max number of soil layers
    else
        maxlayers = 1
    end

    n = length(inds)

    cfmax =
        ncread(
            nc,
            param(config, "input.vertical.cfmax", nothing);
            sel = inds,
            defaults = 3.75653,
            type = Float,
        ) .* (Δt / basetimestep)
    tt = ncread(
        nc,
        param(config, "input.vertical.tt", nothing);
        sel = inds,
        defaults = 0.0,
        type = Float,
    )
    tti = ncread(
        nc,
        param(config, "input.vertical.tti", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float,
    )
    ttm = ncread(
        nc,
        param(config, "input.vertical.ttm", nothing);
        sel = inds,
        defaults = 0.0,
        type = Float,
    )
    whc = ncread(
        nc,
        param(config, "input.vertical.whc", nothing);
        sel = inds,
        defaults = 0.1,
        type = Float,
    )
    w_soil =
        ncread(
            nc,
            param(config, "input.vertical.w_soil", nothing);
            sel = inds,
            defaults = 0.1125,
            type = Float,
        ) .* (Δt / basetimestep)
    cf_soil = ncread(
        nc,
        param(config, "input.vertical.cf_soil", nothing);
        sel = inds,
        defaults = 0.038,
        type = Float,
    )
    # glacier parameters
    g_tt = ncread(
        nc,
        param(config, "input.vertical.g_tt", nothing);
        sel = inds,
        defaults = 0.0,
        type = Float,
        fill = 0.0,
    )
    g_cfmax =
        ncread(
            nc,
            param(config, "input.vertical.g_cfmax", nothing);
            sel = inds,
            defaults = 3.0,
            type = Float,
            fill = 0.0,
        ) .* (Δt / basetimestep)
    g_sifrac = ncread(
        nc,
        param(config, "input.vertical.g_sifrac", nothing);
        sel = inds,
        defaults = 0.001,
        type = Float,
        fill = 0.0,
    )
    glacierfrac = ncread(
        nc,
        param(config, "input.vertical.glacierfrac", nothing);
        sel = inds,
        defaults = 0.0,
        type = Float,
        fill = 0.0,
    )
    glacierstore = ncread(
        nc,
        param(config, "input.vertical.glacierstore", nothing);
        sel = inds,
        defaults = 5500.0,
        type = Float,
        fill = 0.0,
    )
    # soil parameters
    θₛ = ncread(
        nc,
        param(config, "input.vertical.θₛ", nothing);
        sel = inds,
        defaults = 0.6,
        type = Float,
    )
    θᵣ = ncread(
        nc,
        param(config, "input.vertical.θᵣ", nothing);
        sel = inds,
        defaults = 0.01,
        type = Float,
    )
    kv₀ =
        ncread(
            nc,
            param(config, "input.vertical.kv₀", nothing);
            sel = inds,
            defaults = 3000.0,
            type = Float,
        ) .* (Δt / basetimestep)
    f = ncread(
        nc,
        param(config, "input.vertical.f", nothing);
        sel = inds,
        defaults = 0.001,
        type = Float,
    )
    hb = ncread(
        nc,
        param(config, "input.vertical.hb", nothing);
        sel = inds,
        defaults = 10.0,
        type = Float,
    )
    soilthickness = ncread(
        nc,
        param(config, "input.vertical.soilthickness", nothing);
        sel = inds,
        defaults = 2000.0,
        type = Float,
    )
    infiltcappath =
        ncread(
            nc,
            param(config, "input.vertical.infiltcappath", nothing);
            sel = inds,
            defaults = 10.0,
            type = Float,
        ) .* (Δt / basetimestep)
    infiltcapsoil =
        ncread(
            nc,
            param(config, "input.vertical.infiltcapsoil", nothing);
            sel = inds,
            defaults = 100.0,
            type = Float,
        ) .* (Δt / basetimestep)
    maxleakage =
        ncread(
            nc,
            param(config, "input.vertical.maxleakage", nothing);
            sel = inds,
            defaults = 0.0,
            type = Float,
        ) .* (Δt / basetimestep)

    c = ncread(
        nc,
        param(config, "input.vertical.c", nothing);
        sel = inds,
        defaults = 10.0,
        type = Float,
        dimname = :layer,
    )
    if size(c, 1) != maxlayers
        parname = param(config, "input.vertical.c")
        size1 = size(c, 1)
        error("$parname needs a layer dimension of size $maxlayers, but is $size1")
    end
    kvfrac = ncread(
        nc,
        param(config, "input.vertical.khfrac", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float,
        dimname = :layer,
    )
    if size(kvfrac, 1) != maxlayers
        parname = param(config, "input.vertical.khfrac")
        size1 = size(kvfrac, 1)
        error("$parname needs a layer dimension of size $maxlayers, but is $size1")
    end

    # fraction open water and compacted area (land cover)
    waterfrac = ncread(
        nc,
        param(config, "input.vertical.waterfrac", nothing);
        sel = inds,
        defaults = 0.0,
        type = Float,
    )
    pathfrac = ncread(
        nc,
        param(config, "input.vertical.pathfrac", nothing);
        sel = inds,
        defaults = 0.01,
        type = Float,
    )

    # vegetation parameters
    rootingdepth = ncread(
        nc,
        param(config, "input.vertical.rootingdepth", nothing);
        sel = inds,
        defaults = 750.0,
        type = Float,
    )
    rootdistpar = ncread(
        nc,
        param(config, "input.vertical.rootdistpar", nothing);
        sel = inds,
        defaults = -500.0,
        type = Float,
    )
    capscale = ncread(
        nc,
        param(config, "input.vertical.capscale", nothing);
        sel = inds,
        defaults = 100.0,
        type = Float,
    )
    et_reftopot = ncread(
        nc,
        param(config, "input.vertical.et_reftopot", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float,
    )

    cmax, e_r, canopygapfraction, sl, swood, kext = initialize_canopy(nc, config, inds)

    # these are filled in the loop below
    # TODO see if we can replace this approach
    nlayers = zeros(Int, n)
    act_thickl = zeros(Float, maxlayers, n)
    s_layers = zeros(Float, maxlayers + 1, n)

    for i = 1:n
        if length(config_thicknesslayers) > 0
            act_thickl_, nlayers_ =
                set_layerthickness(soilthickness[i], sumlayers, thicknesslayers)
            s_layers_ = pushfirst(cumsum(act_thickl_), 0.0)

            nlayers[i] = nlayers_
            act_thickl[:, i] = act_thickl_
            s_layers[:, i] = s_layers_
        else
            nlayers[i] = 1
            act_thickl[:, i] = SVector(soilthickness[i])
            s_layers[:, i] = pushfirst(cumsum(SVector(soilthickness[i])), 0.0)
        end
    end

    # needed for derived parameters below
    act_thickl = svectorscopy(act_thickl, Val{maxlayers}())
    θₑ = θₛ .- θᵣ
    soilwatercapacity = soilthickness .* θₑ
    satwaterdepth = 0.85 .* soilwatercapacity # cold state value for satwaterdepth
    zi = max.(0.0, soilthickness .- satwaterdepth ./ θₑ) # cold state value for zi

    # copied to array of sarray below
    vwc = fill(mv, maxlayers, n)
    vwc_perc = fill(mv, maxlayers, n)

    sbm = SBM{Float,maxlayers,maxlayers + 1}(
        Δt = tosecond(Δt),
        maxlayers = maxlayers,
        n = n,
        nlayers = nlayers,
        n_unsatlayers = fill(0, n),
        riverfrac = riverfrac,
        θₛ = θₛ,
        θᵣ = θᵣ,
        kv₀ = kv₀,
        kvfrac = svectorscopy(kvfrac, Val{maxlayers}()),
        hb = hb,
        soilthickness = soilthickness,
        act_thickl = act_thickl,
        sumlayers = svectorscopy(s_layers, Val{maxlayers + 1}()),
        infiltcappath = infiltcappath,
        infiltcapsoil = infiltcapsoil,
        maxleakage = maxleakage,
        waterfrac = max.(waterfrac .- riverfrac, Float(0.0)),
        pathfrac = pathfrac,
        rootingdepth = rootingdepth,
        rootdistpar = rootdistpar,
        capscale = capscale,
        et_reftopot = et_reftopot,
        c = svectorscopy(c, Val{maxlayers}()),
        stemflow = fill(mv, n),
        throughfall = fill(mv, n),
        f = f,
        ustorelayerdepth = zero(act_thickl),
        satwaterdepth = satwaterdepth,
        zi = zi,
        soilwatercapacity = soilwatercapacity,
        canopystorage = zeros(Float, n),
        cmax = cmax,
        canopygapfraction = canopygapfraction,
        e_r = e_r,
        precipitation = fill(mv, n),
        temperature = fill(mv, n),
        potential_evaporation = fill(mv, n),
        pottrans_soil = fill(mv, n),
        transpiration = fill(mv, n),
        ae_ustore = fill(mv, n),
        interception = fill(mv, n),
        soilevap = fill(mv, n),
        soilevapsat = fill(mv, n),
        actcapflux = fill(mv, n),
        actevapsat = fill(mv, n),
        actevap = fill(mv, n),
        runoff_river = fill(mv, n),
        runoff_land = fill(mv, n),
        ae_openw_l = fill(mv, n),
        ae_openw_r = fill(mv, n),
        avail_forinfilt = fill(mv, n),
        actinfilt = fill(mv, n),
        actinfiltsoil = fill(mv, n),
        actinfiltpath = fill(mv, n),
        infiltsoilpath = fill(mv, n),
        infiltexcess = fill(mv, n),
        excesswater = fill(mv, n),
        exfiltsatwater = fill(mv, n),
        exfiltustore = fill(mv, n),
        excesswatersoil = fill(mv, n),
        excesswaterpath = fill(mv, n),
        runoff = fill(mv, n),
        net_runoff_river = fill(mv, n),
        vwc = svectorscopy(vwc, Val{maxlayers}()),
        vwc_perc = svectorscopy(vwc_perc, Val{maxlayers}()),
        rootstore = fill(mv, n),
        vwc_root = fill(mv, n),
        vwc_percroot = fill(mv, n),
        ustoredepth = fill(mv, n),
        transfer = fill(mv, n),
        recharge = fill(mv, n),
        actleakage = fill(mv, n),
        # snow parameters
        cfmax = cfmax,
        tt = tt,
        tti = tti,
        ttm = ttm,
        whc = whc,
        w_soil = w_soil,
        cf_soil = cf_soil,
        snow = zeros(Float, n),
        snowwater = zeros(Float, n),
        rainfallplusmelt = fill(mv, n),
        tsoil = fill(Float(10.0), n),
        # glacier parameters
        g_tt = g_tt,
        g_sifrac = g_sifrac,
        g_cfmax = g_cfmax,
        glacierstore = glacierstore,
        glacierfrac = glacierfrac,
        # Interception related to climatology (leaf_area_index)
        sl = sl,
        swood = swood,
        kext = kext,
        leaf_area_index = fill(mv, n),
        waterlevel_land = fill(mv, n),
        waterlevel_river = zeros(Float, n), #set to zero to account for cells outside river domain
    )

    return sbm

end

function update_until_snow(sbm::SBM, config)

    do_lai = haskey(config.input.vertical, "leaf_area_index")
    modelglacier = get(config.model, "glacier", false)::Bool
    modelsnow = get(config.model, "snow", false)::Bool

    @threads for i = 1:sbm.n
        if do_lai
            cmax = sbm.sl[i] * sbm.leaf_area_index[i] + sbm.swood[i]
            canopygapfraction = exp(-sbm.kext[i] * sbm.leaf_area_index[i])
            ewet =
                (1.0 - exp(-sbm.kext[i] * sbm.leaf_area_index[i])) *
                sbm.potential_evaporation[i]
            e_r =
                sbm.precipitation[i] > 0.0 ?
                min(0.25, ewet / max(0.0001, sbm.precipitation[i])) : 0.0
        else
            cmax = sbm.cmax[i]
            canopygapfraction = sbm.canopygapfraction[i]
            e_r = sbm.e_r[i]
        end

        potential_evaporation = sbm.potential_evaporation[i] * sbm.et_reftopot[i]
        # should we include tempcor in SBM?
        # potential_evaporation = PotenEvap #??

        if Second(sbm.Δt) >= Hour(23)
            throughfall, interception, stemflow, canopystorage = rainfall_interception_gash(
                cmax,
                e_r,
                canopygapfraction,
                sbm.precipitation[i],
                sbm.canopystorage[i],
                potential_evaporation,
            )
            pottrans_soil = max(0.0, potential_evaporation - interception) # now in mm
        else
            netinterception, throughfall, stemflow, leftover, interception, canopystorage =
                rainfall_interception_modrut(
                    sbm.precipitation[i],
                    potential_evaporation,
                    sbm.canopystorage[i],
                    canopygapfraction,
                    cmax,
                )
            pottrans_soil = max(0.0, leftover)  # now in mm
            interception = netinterception
        end

        if modelsnow
            tsoil = sbm.tsoil[i] + sbm.w_soil[i] * (sbm.temperature[i] - sbm.tsoil[i])
            snow, snowwater, snowmelt, rainfallplusmelt, snowfall = snowpack_hbv(
                sbm.snow[i],
                sbm.snowwater[i],
                throughfall + stemflow,
                sbm.temperature[i],
                sbm.tti[i],
                sbm.tt[i],
                sbm.ttm[i],
                sbm.cfmax[i],
                sbm.whc[i],
            )
        end
        # update the outputs and states
        sbm.e_r[i] = e_r
        sbm.cmax[i] = cmax
        sbm.canopygapfraction[i] = canopygapfraction
        sbm.canopystorage[i] = canopystorage
        sbm.interception[i] = interception
        sbm.stemflow[i] = stemflow
        sbm.throughfall[i] = throughfall
        sbm.pottrans_soil[i] = pottrans_soil
        if modelsnow
            sbm.snow[i] = snow
            sbm.snowwater[i] = snowwater
            sbm.tsoil[i] = tsoil
            sbm.rainfallplusmelt[i] = rainfallplusmelt
        end
    end
end

function update_until_recharge(sbm::SBM, config)

    # start dummy variables (should be generated from model reader and from Config.jl TOML)
    soilinfreduction = get(config.model, "soilinfreduction", false)::Bool
    modelglacier = get(config.model, "glacier", false)::Bool
    modelsnow = get(config.model, "snow", false)::Bool
    transfermethod = get(config.model, "transfermethod", false)::Bool
    ust = get(config.model, "whole_ust_available", false)::Bool # should be removed from optional setting and code?

    @threads for i = 1:sbm.n
        if modelsnow
            rainfallplusmelt = sbm.rainfallplusmelt[i]
            if modelglacier
                # Run Glacier module and add the snowpack on-top of it.
                # Estimate the fraction of snow turned into ice (HBV-light).
                # Estimate glacier melt.

                snow, snow2glacier, glacierstore, glaciermelt = glacier_hbv(
                    sbm.glacierfrac[i],
                    sbm.glacierstore[i],
                    sbm.snow[i],
                    sbm.temperature[i],
                    sbm.g_tt[i],
                    sbm.g_cfmax[i],
                    sbm.g_sifrac[i],
                    Second(sbm.Δt),
                )
                # Convert to mm per grid cell and add to snowmelt
                glaciermelt = glaciermelt * sbm.glacierfrac[i]
                rainfallplusmelt = rainfallplusmelt + glaciermelt

            end
        else
            rainfallplusmelt = sbm.stemflow[i] + sbm.throughfall[i]
        end

        avail_forinfilt = rainfallplusmelt
        ustoredepth = sum(@view sbm.ustorelayerdepth[i][1:sbm.nlayers[i]])

        runoff_river = min(1.0, sbm.riverfrac[i]) * avail_forinfilt
        runoff_land = min(1.0, sbm.waterfrac[i]) * avail_forinfilt
        avail_forinfilt = max(avail_forinfilt - runoff_river - runoff_land, 0.0)

        rootingdepth = min(sbm.soilthickness[i] * 0.99, sbm.rootingdepth[i])

        ae_openw_r = min(
            sbm.waterlevel_river[i] * sbm.riverfrac[i],
            sbm.riverfrac[i] * sbm.pottrans_soil[i],
        )
        ae_openw_l = min(
            sbm.waterlevel_land[i] * sbm.waterfrac[i],
            sbm.waterfrac[i] * sbm.pottrans_soil[i],
        )

        restevap = sbm.pottrans_soil[i] - ae_openw_r - ae_openw_l

        # evap available for soil evaporation and transpiration
        potsoilevap = restevap * sbm.canopygapfraction[i]
        pottrans = restevap * (1.0 - sbm.canopygapfraction[i])

        # Calculate the initial capacity of the unsaturated store
        ustorecapacity = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i] - ustoredepth

        # Calculate the infiltration flux into the soil column
        infiltsoilpath, infiltsoil, infiltpath, soilinf, pathinf, infiltexcess =
            infiltration(
                avail_forinfilt,
                sbm.pathfrac[i],
                sbm.cf_soil[i],
                sbm.tsoil[i],
                sbm.infiltcapsoil[i],
                sbm.infiltcappath[i],
                ustorecapacity,
                modelsnow,
                soilinfreduction,
            )

        usl, n_usl = set_layerthickness(sbm.zi[i], sbm.sumlayers[i], sbm.act_thickl[i])
        z = cumsum(usl)
        usld = sbm.ustorelayerdepth[i]

        ast = 0.0
        soilevapunsat = 0.0
        if n_usl > 0
            # Using the surface infiltration rate, calculate the flow rate between the
            # different soil layers that contain unsaturated storage assuming gravity
            # based flow only, estimate the gravity based flux rate to the saturated zone
            # (ast) and the updated unsaturated storage for each soil layer.
            if transfermethod && sbm.maxlayers == 1
                ustorelayerdepth = sbm.ustorelayerdepth[i][1] + infiltsoilpath
                kv_z = sbm.kvfrac[i][1] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.zi[i])
                ustorelayerdepth, ast = unsatzone_flow_sbm(
                    ustorelayerdepth,
                    sbm.soilwatercapacity[i],
                    sbm.satwaterdepth[i],
                    kv_z,
                    usl[1],
                    sbm.θₛ[i],
                    sbm.θᵣ[i],
                )
                usld = setindex(usld, ustorelayerdepth, m)
            else
                for m = 1:n_usl
                    l_sat = usl[m] * (sbm.θₛ[i] - sbm.θᵣ[i])
                    kv_z = sbm.kvfrac[i][m] * sbm.kv₀[i] * exp(-sbm.f[i] * z[m])
                    ustorelayerdepth =
                        m == 1 ? sbm.ustorelayerdepth[i][m] + infiltsoilpath :
                        sbm.ustorelayerdepth[i][m] + ast
                    ustorelayerdepth, ast =
                        unsatzone_flow_layer(ustorelayerdepth, kv_z, l_sat, sbm.c[i][m])
                    usld = setindex(usld, ustorelayerdepth, m)
                end
            end

            # then evapotranspiration from layers
            # Calculate saturation deficit
            saturationdeficit = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i]

            # First calculate the evaporation of unsaturated storage into the
            # atmosphere from the upper layer.
            if sbm.maxlayers == 1
                soilevapunsat =
                    potsoilevap * min(1.0, saturationdeficit / sbm.soilwatercapacity[i])
            else
                # In case only the most upper soil layer contains unsaturated storage
                if n_usl == 1
                    # Check if groundwater level lies below the surface
                    soilevapunsat =
                        potsoilevap *
                        min(1.0, usld[1] / (sbm.zi[i] * (sbm.θₛ[i] - sbm.θᵣ[i])))
                else
                    # In case first layer contains no saturated storage
                    soilevapunsat =
                        potsoilevap *
                        min(1.0, usld[1] / (usl[1] * ((sbm.θₛ[i] - sbm.θᵣ[i]))))
                end
            end
            # Ensure that the unsaturated evaporation rate does not exceed the
            # available unsaturated moisture
            soilevapunsat = min(soilevapunsat, usld[1])
            # Update the additional atmospheric demand
            potsoilevap = potsoilevap - soilevapunsat
            usld = setindex(usld, usld[1] - soilevapunsat, 1)
        end
        transfer = ast

        if sbm.maxlayers == 1
            soilevapsat = 0.0
        else
            if n_usl == 0 || n_usl == 1
                soilevapsat =
                    potsoilevap *
                    min(1.0, (sbm.act_thickl[i][1] - sbm.zi[i]) / sbm.act_thickl[i][1])
                soilevapsat = min(
                    soilevapsat,
                    (sbm.act_thickl[i][1] - sbm.zi[i]) * (sbm.θₛ[i] - sbm.θᵣ[i]),
                )
            else
                soilevapsat = 0.0
            end
        end
        soilevap = soilevapunsat + soilevapsat
        satwaterdepth = sbm.satwaterdepth[i] - soilevapsat

        # transpiration from saturated store
        wetroots = scurve(sbm.zi[i], rootingdepth, Float(1.0), sbm.rootdistpar[i])
        actevapsat = min(pottrans * wetroots, satwaterdepth)
        satwaterdepth = satwaterdepth - actevapsat
        restpottrans = pottrans - actevapsat

        # actual transpiration from ustore
        actevapustore = 0.0
        for k = 1:n_usl
            ustorelayerdepth, actevapustore, restpottrans = acttransp_unsat_sbm(
                rootingdepth,
                usld[k],
                sbm.sumlayers[i][k],
                restpottrans,
                actevapustore,
                sbm.c[i][k],
                usl[k],
                sbm.θₛ[i],
                sbm.θᵣ[i],
                sbm.hb[i],
                ust,
            )
            usld = setindex(usld, ustorelayerdepth, k)
        end

        # check soil moisture balance per layer
        du = 0.0
        for k = n_usl:-1:1
            du = max(0.0, usld[k] - usl[k] * (sbm.θₛ[i] - sbm.θᵣ[i]))
            usld = setindex(usld, usld[k] - du, k)
            if k > 1
                usld = setindex(usld, usld[k-1] + du, k - 1)
            end
        end

        actinfilt = infiltsoilpath - du
        excesswater = avail_forinfilt - infiltsoilpath - infiltexcess + du

        # Separation between compacted and non compacted areas (correction with the satflow du)
        # This is required for D-Emission/Delwaq
        if infiltsoil + infiltpath > 0.0
            actinfiltsoil = infiltsoil - du * infiltsoil / (infiltpath + infiltsoil)
            actinfiltpath = infiltpath - du * infiltpath / (infiltpath + infiltsoil)
        else
            actinfiltsoil = 0.0
            actinfiltpath = 0.0
        end
        excesswatersoil = max(soilinf - actinfiltsoil, 0.0)
        excesswaterpath = max(pathinf - actinfiltpath, 0.0)

        actcapflux = 0.0
        if n_usl > 0
            ksat = sbm.kvfrac[i][n_usl] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.zi[i])
            ustorecapacity =
                sbm.soilwatercapacity[i] - satwaterdepth - sum(@view usld[1:sbm.nlayers[i]])
            maxcapflux = max(0.0, min(ksat, actevapustore, ustorecapacity, satwaterdepth))

            if sbm.zi[i] > rootingdepth
                capfluxscale =
                    sbm.capscale[i] / (sbm.capscale[i] + sbm.zi[i] - rootingdepth) *
                    (sbm.Δt / tosecond(basetimestep))
            else
                capfluxscale = 0.0
            end
            capflux = maxcapflux * capfluxscale

            netcapflux = capflux
            for k = n_usl:-1:1
                toadd =
                    min(netcapflux, max(usl[k] * (sbm.θₛ[i] - sbm.θᵣ[i]) - usld[k], 0.0))
                usld = setindex(usld, usld[k] + toadd, k)
                netcapflux = netcapflux - toadd
                actcapflux = actcapflux + toadd
            end
        end
        deepksat = sbm.kv₀[i] * exp(-sbm.f[i] * sbm.soilthickness[i])
        deeptransfer = min(satwaterdepth, deepksat)
        actleakage = max(0.0, min(sbm.maxleakage[i], deeptransfer))

        # recharge (mm) for saturated zone
        recharge = (transfer - actcapflux - actleakage - actevapsat - soilevapsat)
        transpiration = actevapsat + actevapustore
        actevap = soilevap + transpiration + ae_openw_r + ae_openw_l

        # update the outputs and states
        sbm.n_unsatlayers[i] = n_usl
        sbm.net_runoff_river[i] = runoff_river - ae_openw_r
        sbm.avail_forinfilt[i] = avail_forinfilt
        sbm.actinfilt[i] = actinfilt
        sbm.infiltexcess[i] = infiltexcess
        sbm.recharge[i] = recharge
        sbm.transpiration[i] = transpiration
        sbm.soilevap[i] = soilevap
        sbm.soilevapsat[i] = soilevapsat
        sbm.ae_openw_r[i] = ae_openw_r
        sbm.ae_openw_l[i] = ae_openw_l
        sbm.runoff_land[i] = runoff_land
        sbm.runoff_river[i] = runoff_river
        sbm.actevapsat[i] = actevapsat
        sbm.actevap[i] = actevap
        sbm.ae_ustore[i] = actevapustore
        sbm.ustorelayerdepth[i] = usld
        sbm.transfer[i] = transfer
        sbm.actcapflux[i] = actcapflux
        sbm.actleakage[i] = actleakage
        sbm.actinfiltsoil[i] = actinfiltsoil
        sbm.actinfiltpath[i] = actinfiltpath
        sbm.excesswater[i] = excesswater
        sbm.excesswatersoil[i] = excesswatersoil
        sbm.excesswaterpath[i] = excesswaterpath
        sbm.rainfallplusmelt[i] = rainfallplusmelt
        sbm.infiltsoilpath[i] = infiltsoilpath
        sbm.satwaterdepth[i] = satwaterdepth
        if modelsnow
            if modelglacier
                sbm.snow[i] = snow
                sbm.glacierstore[i] = glacierstore
            end
        end
    end
end

function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)

    @threads for i = 1:sbm.n
        usl, n_usl = set_layerthickness(zi[i], sbm.sumlayers[i], sbm.act_thickl[i])
        # exfiltration from ustore
        usld = sbm.ustorelayerdepth[i]
        exfiltustore = 0.0
        for k = sbm.n_unsatlayers[i]:-1:1
            if k <= n_usl
                exfiltustore = max(0, usld[k] - usl[k] * (sbm.θₛ[i] - sbm.θᵣ[i]))
            else
                exfiltustore = usld[k]
            end
            usld = setindex(usld, usld[k] - exfiltustore, k)
            if k > 1
                usld = setindex(usld, usld[k-1] + exfiltustore, k - 1)
            end
        end

        ustoredepth = sum(@view usld[1:n_usl])

        runoff = max(
            exfiltustore +
            exfiltsatwater[i] +
            sbm.excesswater[i] +
            sbm.runoff_land[i] +
            sbm.infiltexcess[i] - sbm.ae_openw_l[i],
            0.0,
        )

        # volumetric water content per soil layer and root zone
        vwc = sbm.vwc[i]
        vwc_perc = sbm.vwc_perc[i]
        for k = 1:sbm.nlayers[i]
            if k <= n_usl
                vwc = setindex(
                    vwc,
                    (usld[k] + (sbm.act_thickl[i][k] - usl[k]) * (sbm.θₛ[i] - sbm.θᵣ[i])) / sbm.act_thickl[i][k] + sbm.θᵣ[i],
                    k,
                )
            else
                vwc = setindex(vwc, sbm.θₛ[i], k)
            end
            vwc_perc = setindex(vwc_perc, (vwc[k] / sbm.θₛ[i]) * 100.0, k)
        end

        rootstore_unsat = 0
        for k = 1:n_usl
            rootstore_unsat =
                rootstore_unsat +
                min(1.0, (max(0.0, sbm.rootingdepth[i] - sbm.sumlayers[i][k]) / usl[k])) *
                usld[k]
        end

        rootstore_sat = max(0.0, sbm.rootingdepth[i] - zi[i]) * (sbm.θₛ[i] - sbm.θᵣ[i])
        rootstore = rootstore_sat + rootstore_unsat
        vwc_root = rootstore / sbm.rootingdepth[i] + sbm.θᵣ[i]
        vwc_percroot = (vwc_root / sbm.θₛ[i]) * 100.0

        satwaterdepth = (sbm.soilthickness[i] - zi[i]) * (sbm.θₛ[i] - sbm.θᵣ[i])

        # update the outputs and states
        sbm.n_unsatlayers[i] = n_usl
        sbm.ustorelayerdepth[i] = usld
        sbm.ustoredepth[i] = ustoredepth
        sbm.satwaterdepth[i] = satwaterdepth
        sbm.exfiltsatwater[i] = exfiltsatwater[i]
        sbm.exfiltustore[i] = exfiltustore
        sbm.runoff[i] = runoff
        sbm.vwc[i] = vwc
        sbm.vwc_perc[i] = vwc_perc
        sbm.rootstore[i] = rootstore
        sbm.vwc_root[i] = vwc_root
        sbm.vwc_percroot[i] = vwc_percroot
        sbm.zi[i] = zi[i]
    end
end

ext: jl


text:
@static if Base.VERSION > v"1.3-" && !(Sys.iswindows() && Sys.WORD_SIZE == 32)
    using OpenSSL_jll: libcrypto
else
    isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")) ||
        error("PDFIO not properly installed. Please run Pkg.build(\"PDFIO\")")
    include("../deps/deps.jl")
end

using Base: SecretBuffer, SecretBuffer!
import Base: copy

export
    BIO,
    read_cddate,
    openssl_error,
    DigestContext,
    CipherContext,
    set_padding,
    reset,
    init,
    update!,
    close,
    CertStore,
    set_params!,
    Cert,
    is_self_signed,
    find_cert,
    get_info,
    PKey,
    CMSContentInfo,
    CMSSignedInfo,
    get_signer,
    get_signer_info,
    get_signer_info_signing_time,
    get_signer_info_timestamp,
    harvest_certs,
    CMS_BINARY,
    CMS_NO_SIGNER_CERT_VERIFY,
    CMS_CADES,
    verify,
    PKCS1SignedInfo,
    get_rsa_digest,
    BigNum,
    crypto_random,
    read_pkcs12,
    decrypt

const V_ASN1_UNDEF                    = Cint(-1)
const V_ASN1_EOC                      = Cint(0)
const V_ASN1_BOOLEAN                  = Cint(1)
const V_ASN1_INTEGER                  = Cint(2)
const V_ASN1_BIT_STRING               = Cint(3)
const V_ASN1_OCTET_STRING             = Cint(4)
const V_ASN1_NULL                     = Cint(5)
const V_ASN1_OBJECT                   = Cint(6)
const V_ASN1_OBJECT_DESCRIPTOR        = Cint(7)
const V_ASN1_EXTERNAL                 = Cint(8)
const V_ASN1_REAL                     = Cint(9)
const V_ASN1_ENUMERATED               = Cint(10)
const V_ASN1_UTF8STRING               = Cint(12)
const V_ASN1_SEQUENCE                 = Cint(16)
const V_ASN1_SET                      = Cint(17)
const V_ASN1_NUMERICSTRING            = Cint(18)
const V_ASN1_PRINTABLESTRING          = Cint(19)
const V_ASN1_T61STRING                = Cint(20)
const V_ASN1_TELETEXSTRING            = Cint(20)
const V_ASN1_VIDEOTEXSTRING           = Cint(21)
const V_ASN1_IA5STRING                = Cint(22)
const V_ASN1_UTCTIME                  = Cint(23)
const V_ASN1_GENERALIZEDTIME          = Cint(24)
const V_ASN1_GRAPHICSTRING            = Cint(25)
const V_ASN1_ISO64STRING              = Cint(26)
const V_ASN1_VISIBLESTRING            = Cint(26)
const V_ASN1_GENERALSTRING            = Cint(27)
const V_ASN1_UNIVERSALSTRING          = Cint(28)
const V_ASN1_BMPSTRING                = Cint(30)

function openssl_error(ret)
    ret <= 0 || return nothing
    errmsg = "Error in Crypto Library: "
    bio = BIO()
    ccall((:ERR_print_errors, libcrypto), Cvoid, (Ptr{Cvoid}, ), bio.data)
    data = read(bio)
    msg = transcode(String, data)
    error(errmsg*"\n\n"*msg)
end

function crypto_random(n::Int)
    b = Vector{UInt8}(undef, n)
    ret = ccall((:RAND_bytes, libcrypto), Cint,
                (Ptr{Cuchar}, Cint), pointer(b), Cint(n))
    openssl_error(ret)
    return b
end

mutable struct BigNum
    n::Ptr{Cvoid}
    function BigNum(b::Vector{UInt8})
        bn = ccall((:BN_bin2bn, libcrypto), Ptr{Cvoid},
                   (Ptr{Cuchar}, Cint, Ptr{Cvoid}),
                   pointer(b), Cint(length(b)), C_NULL)
        bn == C_NULL && error("Cannot instantiate BigNum.")
        this = new(bn)
        finalizer(x->ccall((:BN_free, libcrypto),
                           Cvoid, (Ptr{Cvoid}, ), x.n), this)
        return this
    end
end

Base.mod(n1::BigNum, n2::UInt) =
    ccall((:BN_mod_word, libcrypto), Culonglong,
          (Ptr{Cvoid}, Culonglong), n1.n, Culonglong(n2))

#====

    BIO

=====#

# Memory BIO from OpenSSL BIO interface as a Julia IO
mutable struct BIO <: IO
    data::Ptr{Cvoid}
    function BIO(bio::Ptr{Cvoid})
        this = new(bio)
        finalizer(x->ccall((:BIO_vfree, libcrypto), Cvoid,
                           (Ptr{Cvoid}, ), x.data), this)
    end
end

# R or W File BIO
function BIO(fname::String, mode::String="r")
    bio = ccall((:BIO_new_file, libcrypto), Ptr{Cvoid},
                (Ptr{Cstring}, Ptr{Cstring}), pointer(fname), pointer(mode))
    bio == C_NULL && openssl_error(0)
    return BIO(bio)
end

# Readonly memory BIO
function BIO(membuf::Vector{UInt8})
    bio = ccall((:BIO_new_mem_buf, libcrypto), Ptr{Cvoid},
                (Ptr{Cvoid}, Cint), pointer(membuf), length(membuf))
    bio == C_NULL && openssl_error(0)
    return BIO(bio)
end

# R/W Memory BIO
function BIO()
    buf_bio = ccall((:BIO_s_mem, libcrypto), Ptr{Cvoid}, ())
    buf_bio == C_NULL && openssl_error(0)
    bio = ccall((:BIO_new, libcrypto), Ptr{Cvoid}, (Ptr{Cvoid}, ), buf_bio)
    bio == C_NULL && openssl_error(0)
    return BIO(bio)
end

function Base.read(bio::BIO)
    n, size = Csize_t(4096), Csize_t(4096)
    out, buf = Vector{UInt8}(undef, 0), Vector{UInt8}(undef, size)
    while n >= size
        n = ccall((:BIO_read, libcrypto), Cint,
                  (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
                  bio.data, pointer(buf), size)
        n > 0 && append!(out, @view buf[1:n])
    end
    return out
end

function get_date_from_utctime(d::String)
    # Remove fraction seconds if there
    d[end] == 'Z' || error(E_INVALID_DATE)
    a = split(d, '.')
    length(a) > 1 && (d = a[1]*"Z")
    pfx = length(d) == 15 ? "D:" :
          parse(Int, d[1:2]) < 50 ? "D:20" : "D:19"
    return CDDate(pfx*d)
end

read_cddate(bio::BIO) =
    bio |> read |> pointer |> unsafe_string |> get_date_from_utctime

#====

    DigestContext

====#

mutable struct DigestContext
    data::Ptr{Cvoid}
    md::Ptr{Cvoid}
    function DigestContext(algo::String)
        md = ccall((:EVP_get_digestbyname, libcrypto), Ptr{Cvoid},
                   (Ptr{Cstring}, ), pointer(algo))
        md == C_NULL && error("Unable to find message digest algorithm "*algo)
        data = ccall((:EVP_MD_CTX_new, libcrypto), Ptr{Cvoid}, ())
        data == C_NULL && error("Unable to create digest context")
        this = new(data, md)
        finalizer(x->ccall((:EVP_MD_CTX_free, libcrypto), Cvoid,
                           (Ptr{Cvoid}, ), x.data), this)
        init(this)
        return this
    end
end

function init(ctx::DigestContext)
    ret = ccall((:EVP_DigestInit_ex, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), ctx.data, ctx.md, C_NULL)
    openssl_error(ret)
    return nothing
end

function Base.reset(ctx::DigestContext)
    ret = ccall((:EVP_MD_CTX_reset, libcrypto), Cint, (Ptr{Cvoid}, ), ctx.data)
    openssl_error(ret)
    return nothing
end

function update!(ctx::DigestContext, v::AbstractVector{UInt8})
    ret = ccall((:EVP_DigestUpdate, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cuchar}, Cint),
                ctx.data, pointer(v), length(v))
    openssl_error(ret)
    return nothing
end

update!(ctx::DigestContext, s::SecretBuffer) =
    update!(ctx, (@view s.data[1:s.size]))

const EVP_MAX_MD_SIZE = Cint(64)

function Base.close(ctx::DigestContext)
    md_value = Vector{UInt8}(undef, EVP_MAX_MD_SIZE)
    md_len = Ref(EVP_MAX_MD_SIZE)
    ret = ccall((:EVP_DigestFinal_ex, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cuchar}, Ptr{Cint}),
                ctx.data, md_value, md_len)
    openssl_error(ret)
    resize!(md_value, md_len[])
    return md_value
end

#====

    CipherContext

====#

mutable struct CipherContext
    data::Ptr{Cvoid}
    ca::Ptr{Cvoid}
    function CipherContext(algo::String, key::SecretBuffer,
                           iv::SecretBuffer, isencrypt::Bool)
        ca = ccall((:EVP_get_cipherbyname, libcrypto), Ptr{Cvoid},
                   (Ptr{Cstring}, ), pointer(algo))
        if ca == C_NULL
            ca = (algo == "aes_128_cbc") ?
                ccall((:EVP_aes_128_cbc, libcrypto), Ptr{Cvoid}, ()) :
                (algo == "aes_256_cbc") ?
                ccall((:EVP_aes_256_cbc, libcrypto), Ptr{Cvoid}, ()) :
                (algo == "aes_256_ecb") ?
                ccall((:EVP_aes_256_ecb, libcrypto), Ptr{Cvoid}, ()) :
                error("Unable to find message digest algorithm "*algo)
        end
        data = ccall((:EVP_CIPHER_CTX_new, libcrypto), Ptr{Cvoid}, ())
        data == C_NULL && error("Unable to create cipher context")
        this = new(data, ca)
        finalizer(x->ccall((:EVP_CIPHER_CTX_free, libcrypto), Cvoid,
                           (Ptr{Cvoid}, ), x.data), this)
        init(this, key, iv, isencrypt)
        return this
    end
end

function init(ctx::CipherContext, key::SecretBuffer,
              iv::SecretBuffer, isencrypt::Bool)
    enc = Cint(isencrypt ? 1 : 0)
    piv = iv.size > 0 ? pointer(iv.data) : C_NULL
    ret = ccall((:EVP_CipherInit_ex, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
                 Ptr{Cuchar}, Ptr{Cuchar}, Cint),
                ctx.data, ctx.ca, C_NULL, pointer(key.data), piv, enc)
    openssl_error(ret)
    ret = ccall((:EVP_CIPHER_CTX_set_key_length, libcrypto), Cint,
                (Ptr{Cvoid}, Cint), ctx.data, key.size)
    openssl_error(ret)
    return nothing
end

function set_padding(ctx::CipherContext, pad_size::Int)
    ret = ccall((:EVP_CIPHER_CTX_set_padding, libcrypto), Cint,
                (Ptr{Cvoid}, Cint), ctx.data, Cint(pad_size))
    openssl_error(ret)
end

function Base.reset(ctx::CipherContext)
    ret = ccall((:EVP_CIPHER_CTX_reset, libcrypto), Cint,
                (Ptr{Cvoid}, ), ctx.data)
    openssl_error(ret)
end

const EVP_MAX_BLOCK_LENGTH = Cint(32)

function update!(ctx::CipherContext, indata::AbstractVector{UInt8})
    inlen = Cint(length(indata))
    out, outlen = Vector{UInt8}(undef, inlen + EVP_MAX_BLOCK_LENGTH),
                  Ref(inlen + EVP_MAX_BLOCK_LENGTH)
    ret = ccall((:EVP_CipherUpdate, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cuchar}, Ptr{Cint}, Ptr{Cuchar}, Cint),
                ctx.data, out, outlen, pointer(indata), inlen)
    openssl_error(ret)
    return resize!(out, outlen[])
end

update!(ctx::CipherContext, s::SecretBuffer) =
    update!(ctx, (@view s.data[1:s.size]))

function Base.close(ctx::CipherContext)
    c_value = Vector{UInt8}(undef, EVP_MAX_BLOCK_LENGTH)
    c_len = Ref(EVP_MAX_BLOCK_LENGTH)
    ret = ccall((:EVP_CipherFinal_ex, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cuchar}, Ptr{Cint}),
                ctx.data, c_value, c_len)
    openssl_error(ret)
    resize!(c_value, c_len[])
    return c_value
end

#====

    CertStore

=====#
mutable struct CertStore
    data::Ptr{Cvoid}
    function CertStore()
        store   = ccall((:X509_STORE_new,   libcrypto), Ptr{Cvoid}, ())
        this = new(store)
        finalizer(x->ccall((:X509_STORE_free, libcrypto),
                           Cvoid, (Ptr{Cvoid}, ), x.data), this)

        flookup = ccall((:X509_LOOKUP_file, libcrypto), Ptr{Cvoid}, ())
        lookup  = ccall((:X509_STORE_add_lookup, libcrypto), Ptr{Cvoid},
                        (Ptr{Cvoid}, Ptr{Cvoid}), store, flookup)

        cacerts   = joinpath(@__DIR__, "..", "data", "certs", "cacerts.pem")
        adoberoot = joinpath(@__DIR__, "..", "data", "certs", "adoberoot.pem")
        # Ensure a cacerts.pem file is always there.
        isfile(cacerts) || cp(adoberoot, cacerts)

        ret = ccall((:X509_STORE_load_locations, libcrypto), Cint,
                    (Ptr{Cvoid}, Ptr{Cstring}, Ptr{Cstring}),
                    store, transcode(UInt8, cacerts), C_NULL)
        openssl_error(ret)

        ccall((:X509_STORE_set_default_paths, libcrypto),
              Cvoid, (Ptr{Cvoid}, ), store)
        return this
    end
end

const X509_PURPOSE_SSL_CLIENT         = Cint(1)
const X509_PURPOSE_SSL_SERVER         = Cint(2)
const X509_PURPOSE_NS_SSL_SERVER      = Cint(3)
const X509_PURPOSE_SMIME_SIGN         = Cint(4)
const X509_PURPOSE_SMIME_ENCRYPT      = Cint(5)
const X509_PURPOSE_CRL_SIGN           = Cint(6)
const X509_PURPOSE_ANY                = Cint(7)
const X509_PURPOSE_OCSP_HELPER        = Cint(8)
const X509_PURPOSE_TIMESTAMP_SIGN     = Cint(9)

function set_params!(store::CertStore, d::Dict)
    param = ccall((:X509_VERIFY_PARAM_new, libcrypto), Ptr{Cvoid}, ())
    depth   = get(d, :depth, 16)
    ccall((:X509_VERIFY_PARAM_set_depth, libcrypto), Cvoid,
          (Ptr{Cvoid}, Cint), param, depth)
    purpose = get(d, :purpose, X509_PURPOSE_ANY)
        ccall((:X509_VERIFY_PARAM_set_purpose, libcrypto), Cvoid,
              (Ptr{Cvoid}, Cint), param, purpose)
    atepoch = get(d, :atepoch, nothing)
    atepoch !== nothing &&
        ccall((:X509_VERIFY_PARAM_set_time, libcrypto), Cvoid,
              (Ptr{Cvoid}, Cint), param, atepoch)
    ccall((:X509_STORE_set1_param, libcrypto), Cvoid,
          (Ptr{Cvoid}, Ptr{Cvoid}), store.data, param)
    ccall((:X509_VERIFY_PARAM_free, libcrypto), Cvoid, (Ptr{Cvoid}, ), param)
    return nothing
end

# abstract type SignedInfo end
mutable struct PKey
    data::Ptr{Cvoid}
    function PKey(p::Ptr{Cvoid}, clean::Bool=true)
        this = new(p)
        finalizer(x->ccall((:EVP_PKEY_free, libcrypto),
                           Cvoid, (Ptr{Cvoid}, ), x.data), this)
        return this
    end
end

function copy(pkey::PKey)
    ret = ccall((:EVP_PKEY_up_ref, libcrypto), Cint, (Ptr{Cvoid}, ), pkey.data)
    openssl_error(ret)
    return pkey
end

struct PKCS1SignedInfo
    sdata::Vector{UInt8}
    function PKCS1SignedInfo(sd::Vector{UInt8})
        str = ccall((:d2i_ASN1_OCTET_STRING, libcrypto), Ptr{Cvoid},
                    (Ptr{Cvoid}, Ptr{Ptr{UInt8}}, Clong),
                    C_NULL, Ref(pointer(sd)), length(sd))
        str == C_NULL && error("Unable to read PKCS#1 SignedInfo data")
        len = ccall((:ASN1_STRING_length, libcrypto), Cint, (Ptr{Cvoid}, ), str)
        p = ccall((:ASN1_STRING_get0_data, libcrypto), Ptr{Cuchar},
                  (Ptr{Cvoid}, ), str)
        sdata = Vector{UInt8}(undef, len)
        unsafe_copyto!(pointer(sdata), p, len)
        new(sdata)
    end
end

mutable struct CMSContentInfo
    cms::Ptr{Cvoid}
    detached::Bool
    function CMSContentInfo(contents::Vector{UInt8},
                            detached::Bool = false)
        cms = ccall((:d2i_CMS_ContentInfo, libcrypto), Ptr{Cvoid},
                    (Ptr{Cvoid}, Ptr{Ptr{UInt8}}, Clong),
                    C_NULL, Ref(pointer(contents)), length(contents))
        cms == C_NULL && error("Unable to read CMS SignedInfo data")
        this = new(cms, detached)
        finalizer(x->ccall((:CMS_ContentInfo_free, libcrypto),
                           Cvoid, (Ptr{Cvoid}, ), x.cms), this)
        return this
    end
end

const CMSSignedInfo = CMSContentInfo

function Base.show(io::IO, si::CMSSignedInfo)
    bio = BIO()
    ret = ccall((:CMS_ContentInfo_print_ctx, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Ptr{Cvoid}),
                bio.data, si.cms, 2, C_NULL)
    openssl_error(ret)
    return print(io, unsafe_string(pointer(read(bio))))
end

const CMS_TEXT                   = Cint(0x1)
const CMS_NOCERTS                = Cint(0x2)
const CMS_NO_CONTENT_VERIFY      = Cint(0x4)
const CMS_NO_ATTR_VERIFY         = Cint(0x8)
const CMS_NOSIGS                 = CMS_NO_CONTENT_VERIFY +
                                   CMS_NO_ATTR_VERIFY
const CMS_NOINTERN               = Cint(0x10)
const CMS_NO_SIGNER_CERT_VERIFY  = Cint(0x20)
const CMS_NOVERIFY               = Cint(0x20)
const CMS_DETACHED               = Cint(0x40)
const CMS_BINARY                 = Cint(0x80)
const CMS_NOATTR                 = Cint(0x100)
const CMS_NOSMIMECAP             = Cint(0x200)
const CMS_NOOLDMIMETYPE          = Cint(0x400)
const CMS_CRLFEOL                = Cint(0x800)
const CMS_STREAM                 = Cint(0x1000)
const CMS_NOCRL                  = Cint(0x2000)
const CMS_PARTIAL                = Cint(0x4000)
const CMS_REUSE_DIGEST           = Cint(0x8000)
const CMS_USE_KEYID              = Cint(0x10000)
const CMS_DEBUG_DECRYPT          = Cint(0x20000)
const CMS_KEY_PARAM              = Cint(0x40000)
const CMS_ASCIICRLF              = Cint(0x80000)
const CMS_CADES                  = Cint(0x100000)

function verify(si::CMSSignedInfo, store::CertStore,
                hdata::Vector{UInt8}, flags::Cint)
    ret = ccall((:CMS_verify, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
                 Ptr{Cvoid}, Ptr{Cvoid}, Cint),
                si.cms, C_NULL, store.data, BIO(hdata).data, C_NULL, flags)
    openssl_error(ret)
    return ret == 1
end

# Can be used only after validation
function get_signer(si::CMSSignedInfo)
    certs = ccall((:CMS_get0_signers, libcrypto), Ptr{Cvoid},
                  (Ptr{Cvoid}, ), si.cms)
    certs == C_NULL && openssl_error(0)
    num = ccall((:OPENSSL_sk_num, libcrypto), Cint, (Ptr{Cvoid}, ), certs)
    num > 1 && error("More than one signer in CMSSignedInfo")
    cert = ccall((:OPENSSL_sk_value, libcrypto), Ptr{Cvoid},
                 (Ptr{Cvoid}, Cint), certs, 0)
    return Cert(cert, si)
end

function get_signer_info(cmsinfo::CMSSignedInfo)
    cms = cmsinfo.cms
    certs = ccall((:CMS_get1_certs, libcrypto), Ptr{Cvoid},
                   (Ptr{Cvoid}, ), cms)
    num = ccall((:OPENSSL_sk_num, libcrypto), Cint, (Ptr{Cvoid}, ), certs)
    openssl_error(num)

    sinfos = ccall((:CMS_get0_SignerInfos, libcrypto), Ptr{Cvoid},
                   (Ptr{Cvoid}, ), cms)
    nsi = ccall((:OPENSSL_sk_num, libcrypto), Cint, (Ptr{Cvoid}, ), sinfos)
    @assert nsi == 1 "PDF specification requires only one SignerInfo"
    si = ccall((:OPENSSL_sk_value, libcrypto), Ptr{Cvoid},
               (Ptr{Cvoid}, Cint), sinfos, 0)
    return si
end

function get_signer_info_signing_time(si::Ptr{Cvoid})
    loc = ccall((:CMS_signed_get_attr_by_NID, libcrypto), Cint,
                (Ptr{Cvoid}, Cint, Cint), si, NID_pkcs9_signingTime, Cint(-1))
    loc < 0 && return nothing
    attr = ccall((:CMS_signed_get_attr, libcrypto), Ptr{Cvoid},
                 (Ptr{Cvoid}, Cint), si, loc)
    st_asn = ccall((:X509_ATTRIBUTE_get0_data, libcrypto), Ptr{Cvoid},
                   (Ptr{Cvoid}, Cint, Cint, Ptr{Cvoid}),
                   si, 0, V_ASN1_IA5STRING, C_NULL)
    st_asn == C_NULL && return nothing
    bio = BIO()
    ret = ccall((:ASN1_STRING_print, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}), bio.data, st_asn)
    openssl_error(ret)
    return read_cddate(bio)
end

function get_signer_info_timestamp(si::Ptr{Cvoid})
    loc = ccall((:CMS_unsigned_get_attr_by_NID, libcrypto), Cint,
                (Ptr{Cvoid}, Cint, Cint),
                si, NID_id_smime_aa_timeStampToken, Cint(-1))
    loc < 0 && return nothing
    attr = ccall((:CMS_unsigned_get_attr, libcrypto), Ptr{Cvoid},
                 (Ptr{Cvoid}, Cint), si, loc)
    attr == C_NULL && return nothing
    dptr = ccall((:X509_ATTRIBUTE_get0_data, libcrypto), Ptr{Cvoid},
                 (Ptr{Cvoid}, Cint, Cint, Ptr{Cvoid}),
                 attr, 0, V_ASN1_SEQUENCE, C_NULL)
    xder = Ref(Ptr{Cuchar}(C_NULL))
    len = ccall((:ASN1_STRING_length, libcrypto), Cint, (Ptr{Cvoid}, ), dptr)
    ptr = ccall((:ASN1_STRING_get0_data, libcrypto), Ptr{Cuchar},
                (Ptr{Cvoid}, ), dptr)
    xptr = Ref(ptr)
    ts = ccall((:d2i_CMS_ContentInfo, libcrypto), Ptr{Cvoid},
               (Ptr{Cvoid}, Ptr{Ptr{Cuchar}}, Clong), C_NULL, xptr, len)
    ts == C_NULL && return nothing
    eobj = ccall((:CMS_get0_eContentType, libcrypto), Ptr{Cvoid},
                   (Ptr{Cvoid}, ), ts)
    nid = ccall((:OBJ_obj2nid, libcrypto), Cint, (Ptr{Cvoid}, ), eobj)
    nid != NID_id_smime_ct_TSTInfo && return nothing
    ptinfo = ccall((:CMS_get0_content, libcrypto), Ptr{Ptr{Cvoid}},
                   (Ptr{Cvoid}, ), ts)
    tinfo = unsafe_load(ptinfo)
    tinfo == C_NULL && return nothing
    len  = ccall((:ASN1_STRING_length, libcrypto), Cint, (Ptr{Cvoid}, ), tinfo)
    cbuf = ccall((:ASN1_STRING_get0_data, libcrypto), Ptr{Cuchar},
                 (Ptr{Cuchar}, ), tinfo)
    tst_info = ccall((:d2i_TS_TST_INFO, libcrypto), Ptr{Cvoid},
                     (Ptr{Cvoid}, Ptr{Ptr{Cuchar}}, Clong),
                     C_NULL, Ref(cbuf), len)
    gentime = ccall((:TS_TST_INFO_get_time, libcrypto), Ptr{Cvoid},
                    (Ptr{Cvoid}, ), tst_info)
    bio = BIO()
    ret = ccall((:ASN1_STRING_print, libcrypto), Cint,
                    (Ptr{Cvoid}, Ptr{Cvoid}), bio.data, gentime)
    openssl_error(ret)
    return read_cddate(bio)
end

function harvest_certs(cmsinfo::CMSSignedInfo)
    cms = cmsinfo.cms
    certs = ccall((:CMS_get1_certs, libcrypto), Ptr{Cvoid},
                   (Ptr{Cvoid}, ), cms)
    num = ccall((:OPENSSL_sk_num, libcrypto), Cint, (Ptr{Cvoid}, ), certs)
    openssl_error(num)

    si = get_signer_info(cmsinfo)
    xcert = Ref(C_NULL)
    ccall((:CMS_SignerInfo_get0_algs, libcrypto), Cvoid,
          (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Ptr{Cvoid}),
          si, C_NULL, xcert, C_NULL, C_NULL)
    signed_cert = Cert(xcert[], nothing)

    certs_info = Vector{Dict{Symbol, Any}}()
    for i = 1:num
        cert = ccall((:OPENSSL_sk_value, libcrypto), Ptr{Cvoid},
                     (Ptr{Cvoid}, Cint), certs, i-1)
        # Take ownership of cert clean up
        pdc = Cert(cert)
        push!(certs_info, get_info(pdc))
        if xcert[] == C_NULL
            ret = ccall((:CMS_SignerInfo_cert_cmp, libcrypto), Cint,
                        (Ptr{Cvoid}, Ptr{Cvoid}), si, cert)
            ret == 0 && (signed_cert = pdc)
        end
    end
    ccall((:OPENSSL_sk_free, libcrypto), Cvoid, (Ptr{Cvoid}, ), certs)
    return certs_info, signed_cert
end

mutable struct Cert
    data::Ptr{Cvoid}
    lockref::Union{CMSSignedInfo, Nothing}
end

# Manage life cycle also. Will be freed along with the Julia object
function Cert(cert::Ptr{Cvoid})
    this = Cert(cert, nothing)
    finalizer(x->ccall((:X509_free, libcrypto), Cvoid,
                       (Ptr{Cvoid}, ), x.data), this)
    return this
end

function Cert(cbytes::Vector{UInt8})
    cdata = ccall((:d2i_X509, libcrypto), Ptr{Cvoid},
                  (Ptr{Cvoid}, Ptr{Ptr{UInt8}}, Clong),
                  C_NULL, Ref(pointer(cbytes)), length(cbytes))
    cdata == C_NULL && error("Invalid certificate data")
    return Cert(cdata)
end

Base.:(==)(c1::Cert, c2::Cert) = !(c1 < c2 || c2 < c1)
Base.isless(c1::Cert, c2::Cert) =
    ccall((:X509_cmp, libcrypto), Cint,
          (Ptr{Cvoid}, Ptr{Cvoid}), c1.data, c2.data) < 0

function is_self_signed(c::Cert)
    cert = c.data
    ns  = ccall((:X509_get_subject_name, libcrypto), Ptr{Cvoid},
                (Ptr{Cvoid}, ), cert)
    ni  = ccall((:X509_get_issuer_name,  libcrypto), Ptr{Cvoid},
                (Ptr{Cvoid}, ), cert)
    ret = ccall((:X509_NAME_cmp, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}), ns, ni)
    return ret == 0
end

const ASN1_STRFLGS_ESC_2253           = Culong(1)
const ASN1_STRFLGS_ESC_CTRL           = Culong(2)
const ASN1_STRFLGS_ESC_MSB            = Culong(4)
const ASN1_STRFLGS_ESC_QUOTE          = Culong(8)
const ASN1_STRFLGS_UTF8_CONVERT       = Culong(0x10)
const ASN1_STRFLGS_DUMP_UNKNOWN       = Culong(0x100)
const ASN1_STRFLGS_DUMP_DER           = Culong(0x200)

const ASN1_STRFLGS_RFC2253 =
    ASN1_STRFLGS_ESC_2253 +
    ASN1_STRFLGS_ESC_CTRL +
    ASN1_STRFLGS_ESC_MSB  +
    ASN1_STRFLGS_UTF8_CONVERT +
    ASN1_STRFLGS_DUMP_UNKNOWN +
    ASN1_STRFLGS_DUMP_DER

const XN_FLAG_FN_SN                   = Culong(0)
const XN_FLAG_SEP_CPLUS_SPC           = Culong(2 << 16)
const XN_FLAG_SPC_EQ                  = Culong(1 << 23)

const XN_FLAG_ONELINE =
    ASN1_STRFLGS_ESC_2253 +
    ASN1_STRFLGS_ESC_CTRL +
    ASN1_STRFLGS_ESC_MSB +
    ASN1_STRFLGS_ESC_QUOTE +
    ASN1_STRFLGS_UTF8_CONVERT +
    ASN1_STRFLGS_DUMP_UNKNOWN +
    ASN1_STRFLGS_DUMP_DER +
    XN_FLAG_SEP_CPLUS_SPC +
    XN_FLAG_SPC_EQ +
    XN_FLAG_FN_SN

const NID_commonName  = Cint(13)
const NID_pkcs9_signingTime = Cint(52)
const NID_id_smime_aa_timeStampToken = Cint(225)
const NID_id_smime_ct_TSTInfo = Cint(207)

const X509_V_OK = Cint(0)

function get_info(c::Cert)
    info = Dict{Symbol, Any}()
    ns  = ccall((:X509_get_subject_name, libcrypto), Ptr{Cvoid},
                (Ptr{Cvoid}, ), c.data)
    ni  = ccall((:X509_get_issuer_name,  libcrypto), Ptr{Cvoid},
                (Ptr{Cvoid}, ), c.data)
    bio = BIO()
    ret = ccall((:X509_NAME_print_ex, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Culong),
                bio.data, ns, 0, XN_FLAG_ONELINE)
    openssl_error(ret)
    info[:subject] = unsafe_string(pointer(read(bio)))
    ret = ccall((:X509_NAME_print_ex, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Culong),
                bio.data, ni, 0, XN_FLAG_ONELINE)
    openssl_error(ret)
    info[:issuer]  = unsafe_string(pointer(read(bio)))

    lastpos = Cint(-1)

    while true
        lastpos = ccall((:X509_NAME_get_index_by_NID, libcrypto), Cint,
                        (Ptr{Cvoid}, Cint, Cint), ns, NID_commonName, lastpos)
        lastpos == -1 && break
        e = ccall((:X509_NAME_get_entry, libcrypto), Ptr{Cvoid},
                  (Ptr{Cvoid}, Cint), ns, lastpos)
        ed = ccall((:X509_NAME_ENTRY_get_data, libcrypto), Ptr{Cvoid},
                   (Ptr{Cvoid}, ), e)
        ret = ccall((:ASN1_STRING_print_ex, libcrypto), Cint,
                    (Ptr{Cvoid}, Ptr{Cvoid}, Cint),
                    bio.data, ed, ASN1_STRFLGS_RFC2253)
        openssl_error(ret)
        cn = unsafe_string(pointer(read(bio)))
        cns = get!(info, :cn, CDTextString[])
        push!(cns, cn)
    end

    x509_notBefore = ccall((:X509_get0_notBefore, libcrypto), Ptr{Cvoid},
                           (Ptr{Cvoid}, ), c.data)
    ret = ccall((:ASN1_STRING_print, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}), bio.data, x509_notBefore)
    openssl_error(ret)
    info[:notBefore] = read_cddate(bio)

    x509_notAfter  = ccall((:X509_get0_notAfter,  libcrypto), Ptr{Cvoid},
                           (Ptr{Cvoid}, ), c.data)
    ret = ccall((:ASN1_STRING_print, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}), bio.data, x509_notAfter)
    openssl_error(ret)
    info[:notAfter]  = read_cddate(bio)
    info[:text] = string(c)
    return info
end

# Flags for X509_print_ex()

const X509_FLAG_COMPAT                = Cint(0)
const X509_FLAG_NO_HEADER             = Cint(1)
const X509_FLAG_NO_VERSION            = Cint(1 << 1)
const X509_FLAG_NO_SERIAL             = Cint(1 << 2)
const X509_FLAG_NO_SIGNAME            = Cint(1 << 3)
const X509_FLAG_NO_ISSUER             = Cint(1 << 4)
const X509_FLAG_NO_VALIDITY           = Cint(1 << 5)
const X509_FLAG_NO_SUBJECT            = Cint(1 << 6)
const X509_FLAG_NO_PUBKEY             = Cint(1 << 7)
const X509_FLAG_NO_EXTENSIONS         = Cint(1 << 8)
const X509_FLAG_NO_SIGDUMP            = Cint(1 << 9)
const X509_FLAG_NO_AUX                = Cint(1 << 10)
const X509_FLAG_NO_ATTRIBUTES         = Cint(1 << 11)
const X509_FLAG_NO_IDS                = Cint(1 << 12)

function show(io::IO, cert::Cert)
    bio = BIO()
    ccall((:X509_print_ex, libcrypto), Cint,
          (Ptr{Cvoid}, Ptr{Cvoid}, Culong, Culong),
          bio.data, cert.data, XN_FLAG_ONELINE, X509_FLAG_COMPAT)
    ccall((:PEM_write_bio_X509, libcrypto), Cint,
          (Ptr{Cvoid}, Ptr{Cvoid}), bio.data, cert.data)
    return print(io, unsafe_string(pointer(read(bio))))
end

const X509_LU_X509 = Cint(1)

function find_cert(s::CertStore, c::Cert)
    os = ccall((:X509_STORE_get0_objects, libcrypto), Ptr{Cvoid},
               (Ptr{Cvoid}, ), s.data)
    ns = ccall((:X509_get_subject_name, libcrypto), Ptr{Cvoid},
               (Ptr{Cvoid}, ), c.data)
    o  = ccall((:X509_OBJECT_retrieve_by_subject, libcrypto), Ptr{Cvoid},
               (Ptr{Cvoid}, Cint, Ptr{Cvoid}), os, X509_LU_X509, ns)
    o == C_NULL && return false
    c1 = ccall((:X509_OBJECT_get0_X509, libcrypto),
               Ptr{Cvoid}, (Ptr{Cvoid}, ), o)
    c1 == C_NULL && return false
    return Cert(c1, nothing) == c
end

function verify(store::CertStore, cert::Cert)
    ctx = ccall((:X509_STORE_CTX_new, libcrypto), Ptr{Cvoid}, ())
    try
        ret = ccall((:X509_STORE_CTX_init, libcrypto), Cint,
                    (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
                    ctx, store.data, cert.data, C_NULL)
        openssl_error(ret)

        ret = ccall((:X509_verify_cert, libcrypto), Cint, (Ptr{Cvoid}, ), ctx)
        if ret <= 0
            depth = ccall((:X509_STORE_CTX_get_error_depth, libcrypto), Cint,
                          (Ptr{Cvoid}, ), ctx)
            err   = ccall((:X509_STORE_CTX_get_error, libcrypto), Cint,
                          (Ptr{Cvoid}, ), ctx)
            cestr = ccall((:X509_verify_cert_error_string, libcrypto),
                          Cstring, (Cint, ), err)
            errstr = unsafe_string(cestr)
            return false, errstr
        else
            chain = ccall((:X509_STORE_CTX_get1_chain, libcrypto), Ptr{Cvoid},
                          (Ptr{Cvoid}, ), ctx)
            num = ccall((:OPENSSL_sk_num, libcrypto), Cint,
                        (Ptr{Cvoid}, ), chain)
            vrf_chain = Dict{Symbol, Any}[]
            for i = 0:(num-1)
                cert = ccall((:OPENSSL_sk_value, libcrypto), Ptr{Cvoid},
                             (Ptr{Cvoid}, Cint), chain, i)
                pdc = Cert(cert, nothing)
                push!(vrf_chain, get_info(pdc))
                ccall((:X509_free, libcrypto), Cvoid, (Ptr{Cvoid},), cert)
            end
            ccall((:OPENSSL_sk_free, libcrypto), Cvoid, (Ptr{Cvoid}, ), chain)
            return true, vrf_chain
        end
    finally
        ccall((:X509_STORE_CTX_free, libcrypto), Cvoid, (Ptr{Cvoid}, ), ctx)
    end
end

function PKey(cert::Cert)
    kdata = ccall((:X509_get_pubkey, libcrypto), Ptr{Cvoid},
                  (Ptr{Cvoid}, ), cert.data)
    kdata == C_NULL && error("Cannot extract public key from certificate")
    return PKey(kdata)
end

Base.:(==)(k1::PKey, k2::PKey) =
    ccall((:EVP_PKEY_cmp, libcrypto), Cint,
          (Ptr{Cvoid}, Ptr{Cvoid}), k1.data, k2.data) == 1

const RSA_PKCS1_PADDING = Cint(1)

function get_rsa_digest(si::PKCS1SignedInfo, cert::Cert)
    pkey = PKey(cert)
    sdata = si.sdata

    rsa = ccall((:EVP_PKEY_get0_RSA, libcrypto), Ptr{Cvoid},
                (Ptr{Cvoid}, ), pkey.data)
    rsa == C_NULL && error("Invalid RSA certificate")

    len = ccall((:RSA_size, libcrypto), Cint, (Ptr{Cvoid}, ), rsa) - 11
    out = Vector{UInt8}(undef, len)
    len = ccall((:RSA_public_decrypt, libcrypto), Cint,
                (Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cvoid}, Cint),
                length(sdata), pointer(sdata), pointer(out), rsa,
                RSA_PKCS1_PADDING)
    openssl_error(len)

    # X509_SIG is recommendeded to be used for DigestInfo
    sig = ccall((:d2i_X509_SIG, libcrypto), Ptr{Cvoid},
                (Ptr{Cvoid}, Ptr{Ptr{Cuchar}}, Clong),
                C_NULL, Ref(pointer(out)), len)
    try
        xalgor, xdigest = Ref(C_NULL), Ref(C_NULL)
        ccall((:X509_SIG_get0, libcrypto), Cvoid,
              (Ptr{Cvoid}, Ptr{Ptr{Cvoid}}, Ptr{Ptr{Cvoid}}),
              sig, xalgor, xdigest)

        # Digest data extraction from X509_SIG
        c_digest = ccall((:ASN1_STRING_get0_data, libcrypto), Ptr{Cuchar},
                         (Ptr{Cvoid}, ), xdigest[])
        c_len = ccall((:ASN1_STRING_length, libcrypto), Cint,
                      (Ptr{Cvoid}, ), xdigest[])
        xtract_digest = Vector{UInt8}(undef, c_len)
        unsafe_copyto!(pointer(xtract_digest), c_digest, c_len)

        # Compute digest from PDF file data applying the hash algorithm
        # available in the PKCS#1 signature
        xoid, xvtyp, xpval = Ref(C_NULL), Ref(Cint(0)), Ref(C_NULL)
        ccall((:X509_ALGOR_get0, libcrypto), Cvoid,
              (Ptr{Ptr{Cvoid}}, Ptr{Cint}, Ptr{Ptr{Cvoid}}, Ptr{Cvoid}),
              xoid, xvtyp, xpval, xalgor[])

        nid = ccall((:OBJ_obj2nid, libcrypto), Cint, (Ptr{Cvoid}, ), xoid[])
        cstr_sn  = ccall((:OBJ_nid2sn, libcrypto), Cstring, (Cint, ), nid)
        sn = unsafe_string(cstr_sn)
        return xtract_digest, sn
    finally
        ccall((:X509_SIG_free, libcrypto), Cvoid, (Ptr{Cvoid}, ), sig)
    end
end

function read_pkcs12(fn::AbstractString, pw::SecretBuffer)
    data = read(fn)
    p12 = ccall((:d2i_PKCS12, libcrypto), Ptr{Cvoid},
                (Ptr{Cvoid}, Ptr{Ptr{Cuchar}}, Clong),
                C_NULL, Ref(pointer(data)), length(data))
    p12 == C_NULL && error("Unable to read $fn")

    xkey, xcert, xca = Ref(C_NULL), Ref(C_NULL), Ref(C_NULL)

    ret = ccall((:PKCS12_parse, libcrypto), Cint,
                (Ptr{Cvoid}, Ptr{Cstring},
                 Ptr{Ptr{Cvoid}}, Ptr{Ptr{Cvoid}}, Ptr{Ptr{Cvoid}}),
                p12, pointer(pw.data), xkey, xcert, xca)
    openssl_error(ret)
    return Cert(xcert[]), PKey(xkey[])
end

const CMS_RECIPINFO_TRANS = Cint(0)

function decrypt(ci::CMSContentInfo, key::PKey, cert::Cert,
                 detached::Vector{UInt8}, flags::Cint)
    ris = ccall((:CMS_get0_RecipientInfos, libcrypto), Ptr{Cvoid},
                (Ptr{Cvoid}, ), ci.cms)

    nri = ccall((:OPENSSL_sk_num, libcrypto), Cint, (Ptr{Cvoid}, ), ris)

    found = false
    for i = 1:nri
        ri = ccall((:OPENSSL_sk_value, libcrypto), Ptr{Cvoid},
                   (Ptr{Cvoid}, Cint), ris, Cint(i-1))
        if CMS_RECIPINFO_TRANS == ccall((:CMS_RecipientInfo_type, libcrypto),
                                        Cint, (Ptr{Cvoid}, ), ri)
            if ccall((:CMS_RecipientInfo_ktri_cert_cmp, libcrypto), Cint,
                     (Ptr{Cvoid}, Ptr{Cvoid}), ri, cert.data) == 0
                found = true
                kcpy = copy(key)
                ret = ccall((:CMS_RecipientInfo_set0_pkey, libcrypto), Cint,
                            (Ptr{Cvoid}, Ptr{Cvoid}), ri, kcpy.data)
                openssl_error(ret)
                ret = ccall((:CMS_RecipientInfo_decrypt, libcrypto), Cint,
                            (Ptr{Cvoid}, Ptr{Cvoid}), ci.cms, ri)
                openssl_error(ret)
                break
            end
        end
    end

    if found
        out  = BIO()
        dbio = ci.detached ? BIO(detached).data : C_NULL
        ret = ccall((:CMS_decrypt, libcrypto), Cint,
                    (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
                     Ptr{Cvoid}, Ptr{Cvoid}, Cint),
                    ci.cms, C_NULL, C_NULL, dbio, out.data, flags)
        openssl_error(ret)
        return read(out)
    else
        return nothing
    end
end


text:
using LinearAlgebra
using Dates
using Plots
using SparseArrays
using Distributions

include("truncated_svd.jl")
include("resetqr.jl")
include("cresetqr.jl")

"""
Return the rate of convergence of different trials (Figure 3.6 and 3.7)

# Arguments
- `min_rows::Integer`: minimum number of rows the generated matrix can have in each trial
- `min_cols::Integer`:minimum number of columns the generated matrix can have in each trial
- `min_density::Integer`: minimum density the generated matrix can have in each trial. It must be > 0 and < 1. We suggest to not use too small densities.
- `k::Integer`: rank of the UV decomposition
- `num_tests::Integer`: number of trials to perform
- `max_iterations::Integer`: max_iteration of CRESETQR

# Returns
- `x::Array{Float64,1}`: vector for x axis in the plot
- `y::Array{Float64,2}`: matrix of size <num_tests> x <max_iterations>. Row i contains the rate of convergence values for the i-th trial.

# Plots
plot(x, y',  lw=3, xlabel = "iteration", ylabel = "R",legend=:false)

"""
function multiple_convergence_test(min_rows, min_cols, min_density, k, num_tests=25, max_iterations=100)

    #NOTE: all the trials must end with 1000 iterations. Therefore choose the parameters so that it will happen.

    plotly()

    y = zeros(num_tests, max_iterations)
    j = 1

    for i = 1:num_tests
        m = rand(min_rows:750)
        n = rand(min_cols:1000)
        scaleA = rand(1:5)
        densityA = rand(Uniform(min_density,1))

        #generate matrix A
        A = Array(sprand(Float64,m,n,densityA))*scaleA

        _, y1, _, _ = test_convergence_same_starting_point(A,k,max_iterations)

        y[j,:] .= y1[:,1]

        j = j + 1

    end

    x = [1:1:max_iterations;]

    #To generate charts Figure 3.6 and 3.7
    #plot(x, y',  lw=3, xlabel = "iteration", ylabel = "R",legend=:false)

    return x, y
end

"""
To generate rate of convergence, relative_error_x and relative_error_f plots

# Arguments
- `A::Array{Float64,2}`: matrix of size m x n
- `k::Integer`: rank of the UV decomposition
- `max_iterations::Integer=100`: max_iteration of CRESETQR
- `p::Integer=1`: p factor in rate of convergence

# Returns
- `x::Array{Float64,1}`: vector for x axis in the plot
- `y1::Array{Float64,1}`: rate of convergence values during the execution
- `y2::Array{Float64,1}`: relative_error_x values during the execution
- `y3::Array{Float64,1}`: relative_error_f values during the execution

# Plots
- Plot rate of convergence (Figure 3.5):
    - plot(x, y1, lw=3, xlabel = "iteration", ylabel = "R", legend=:false)
- Plot relative_error_x (Figure 3.4a):
    - plot(x, y2, lw=3, yaxis=:log10, ylabel = "log(relative_error_x)", xlabel = "iteration", legend=:false)
- Plot relative_error_f (Figure 3.4b):
    - plot(x, y3, lw=3, yaxis=:log10, ylabel = "log(relative_error_f)", xlabel = "iteration", legend=:false)
"""
function test_convergence1(A, k, max_iterations=100, p=1)
    plotly()
    #generate starting point
    V = rand(k, size(A)[2])
    solution = approx_matrix_svd(A, k)
    #exec cresetqr
    _, _, errors, residuals = cresetqr(A, k; initV=V, opt=solution, max_iterations=max_iterations, min_error = 1e-16, print_error=true, save_errors=true)

    #Generate the chart
    y1 = zeros((length(errors), 1))
    y2 = zeros((length(residuals), 1))
    y3 = zeros((length(residuals), 1))

    #Compute the optimal value using SVD
    optimal_f = norm(A - solution)
    x = [1:1:max(length(errors),1);]

    #Fill y1 axis (rate of convergence)
    for i = 2:length(errors)
        y1[i,1] = (errors[i] - optimal_f)/(errors[i-1] - optimal_f)^p
    end

    #fill y2 axis (relative_error_x)
    for i = 1:length(errors)
        y2[i,1] = (residuals[i])/norm(solution)
    end

    #fill y3 axis (relative_error_f)
    for i = 1:length(errors)
        y3[i,1] = (errors[i] - optimal_f)/optimal_f
    end

    #To generate charts

    #Plot rate of convergence (Figure 3.5)
    #plot(x, y1, lw=3, xlabel = "iteration", ylabel = "R", legend=:false)

    #plot relative_error_x (Figure 3.4a)
    #plot(x, y2, lw=3, yaxis=:log10, ylabel = "log(relative_error_x)", xlabel = "iteration", legend=:false)

    #plot relative_error_f (Figure 3.4b)
    #plot(x, y3, lw=3, yaxis=:log10, ylabel = "log(relative_error_f)", xlabel = "iteration", legend=:false)

    return x, y1, y2, y3
end

"""
To generate plot 3.8 (in the report)

# Arguments
- `A::Array{Float64,2}`: matrix of size m x n
- `k::Integer`: rank of the UV decomposition
- `max_iterations::Integer=100`: max_iteration of CRESETQR
- `p::Integer=1`: p factor in rate of convergence

# Returns
- `x::Array{Float64,1}`: denominator of rate of convergence during the execution
- `y::Array{Float64,1}`: numerator of rate of convergence during the execution

# Plots
plot(x, y, lw=3, scale=:log10, xlabel = "log(f(i-1) - f(*))", ylabel = "log(f(i) - f(*))",legend=:false)

"""
function test_convergence2(A, k, max_iterations=100, p=1)
    plotly()
    #generate starting point
    V = rand(k, size(A)[2])
    solution = approx_matrix_svd(A, k)

    #exec cresetqr
    _, _, errors, _ = cresetqr(A, k; initV=V, opt=solution, max_iterations=max_iterations, min_error = 1e-12, print_error=true, save_errors=true)

    #Generate the chart

    x = zeros(length(errors)-1)
    y = zeros(length(errors)-1)

    #Compute the optimal value using SVD
    optimal_f = norm(A - solution)

    #Fill y1 axis (rate of convergence)
    for i = 1:length(errors)-1
        y[i] = (errors[i+1] - optimal_f)
    end

    for i = 1:length(errors)-1
        x[i] = (errors[i] - optimal_f)^p
    end

    #To generate charts
    #plot(x, y, lw=3, scale=:log10, xlabel = "log(f(i-1) - f(*))", ylabel = "log(f(i) - f(*))",legend=:false)

    return x, y
end


text:
struct AdvancedData{T,V}
    G::SparseMatrixCSC{T,V}
    cc::Vector{Vector{V}}
    nodemap::Matrix{V}
    polymap::Matrix{V}
    hbmeta::RasterMeta
    sources::Vector{T}
    grounds::Vector{T}
    source_map::Matrix{T} # Needed for one to all mode
    finite_grounds::Vector{T}
    check_node::V
    src::V
    cellmap::Matrix{T}
end

function raster_advanced(T, V, cfg)::Matrix{T}

    # Load raster data
    rasterdata = load_raster_data(T, V, cfg)

    # Get flags
    flags = get_raster_flags(cfg)

    # Generate advanced
    advanced_data = compute_advanced_data(rasterdata, flags)

    # Send to main kernel
    v, _ = advanced_kernel(advanced_data, flags, cfg)

    v
end

function compute_advanced_data(data::RasData{T,V},
                        flags)::AdvancedData{T,V} where {T,V}

    # Data
    cellmap = data.cellmap
    polymap = data.polymap
    points_rc = data.points_rc
    included_pairs = data.included_pairs
    hbmeta = data.hbmeta
    source_map = data.source_map

    # Flags
    avg_res = flags.avg_res
    four_neighbors = flags.four_neighbors
    policy = flags.policy

    # Nodemap and graph construction
    nodemap = construct_node_map(cellmap, polymap)
    A = construct_graph(cellmap, nodemap, avg_res, four_neighbors)
    G = laplacian(A)

    # Connected Components
    cc = connected_components(SimpleWeightedGraph(G))

    # Advanced mode specific stuff
    sources, grounds, finite_grounds =
            get_sources_and_grounds(data, flags, G, nodemap)

    AdvancedData(G, cc, nodemap, polymap, hbmeta,
                sources, grounds, source_map,
                finite_grounds, V(-1), V(0), cellmap)
end

function get_sources_and_grounds(data, flags, G, nodemap)

    # Data
    source_map = data.source_map
    ground_map = data.ground_map

    _get_sources_and_grounds(source_map, ground_map, flags, G, nodemap)
end

function _get_sources_and_grounds(source_map, ground_map,
                                  flags, G, nodemap::Matrix{V}, override_policy = :none) where V
    # Flags
    is_raster = flags.is_raster
    grnd_file_is_res = flags.grnd_file_is_res
    policy = override_policy == :none ? flags.policy : override_policy

    # Initialize sources and grounds
    sources = zeros(eltype(G), size(G, 1))
    grounds = zeros(eltype(G), size(G, 1))

    if is_raster
        (i1, j1, v1) = begin _I = findall(!iszero, source_map); getindex.(_I, 1), getindex.(_I, 2), source_map[_I] end
        (i2, j2, v2) = begin _I = findall(!iszero, ground_map); getindex.(_I, 1), getindex.(_I, 2), ground_map[_I] end
        for i = 1:size(i1, 1)
            v = V(nodemap[i1[i], j1[i]])
            if v != 0
                sources[v] += v1[i]
            end
        end
        for i = 1:size(i2, 1)
            v = V(nodemap[i2[i], j2[i]])
            if v != 0
                grounds[v] += v2[i]
            end
        end
    else
        if grnd_file_is_res
            ground_map[:,2] = 1 ./ ground_map[:,2]
        end
        sources[V.(source_map[:,1])] = source_map[:,2]
        grounds[V.(ground_map[:,1])] = ground_map[:,2]
    end
    sources, grounds, finitegrounds =
        resolve_conflicts(sources, grounds, policy)
end

function resolve_conflicts(sources::Vector{T},
                            grounds::Vector{T}, policy) where T

    finitegrounds = similar(sources)
    l = size(sources, 1)

    finitegrounds = map(x -> x < T(Inf) ? x : T(0.), grounds)
    if count(x -> x != 0, finitegrounds) == 0
        finitegrounds = T.([-9999.])
    end

    conflicts = falses(l)
    for i = 1:l
        conflicts[i] = sources[i] != 0 && grounds[i] != 0
    end

    if any(conflicts)
        if policy == :rmvsrc
            sources[findall(x->x!=0,conflicts)] .= 0
        elseif policy == :rmvgnd
            grounds[findall(x->x!=0,conflicts)] .= 0
        elseif policy == :rmvall
            sources[findall(x->x!=0,conflicts)] .= 0
        end
    end

    infgrounds = map(x -> x == Inf, grounds)
    infconflicts = map((x,y) -> x > 0 && y > 0, infgrounds, sources)
    grounds[infconflicts] .= 0

    sources, grounds, finitegrounds
end

function advanced_kernel(data::AdvancedData{T,V}, flags, cfg)::Tuple{Matrix{T},Matrix{T}} where {T,V}

    # Data
    G = data.G
    nodemap = data.nodemap
    polymap = data.polymap
    hbmeta = data.hbmeta
    sources = data.sources
    grounds = data.grounds
    finitegrounds = data.finite_grounds
    cc = data.cc
    src = data.src
    check_node = data.check_node
    source_map = data.source_map # Need it for one to all mode
    cellmap = data.cellmap

    # Flags
    is_raster = flags.is_raster
    is_alltoone = flags.is_alltoone
    is_onetoall = flags.is_onetoall
    write_v_maps = flags.outputflags.write_volt_maps
    write_c_maps = flags.outputflags.write_cur_maps
    write_cum_cur_map_only = flags.outputflags.write_cum_cur_map_only

    volt = zeros(eltype(G), size(nodemap))
    ind = findall(x->x!=0,nodemap)
    f_local = Vector{eltype(G)}()
    solver_called = false
    voltages = Vector{eltype(G)}()
    outvolt = alloc_map(hbmeta)
    outcurr = alloc_map(hbmeta)

    for c in cc

        if check_node != -1 && !(check_node in c)
            continue
        end

        # a_local = laplacian(G[c, c])
        a_local = G[c,c]
        s_local = sources[c]
        g_local = grounds[c]

        if sum(s_local) == 0 || sum(g_local) == 0
            continue
        end

        if finitegrounds != [-9999.]
            f_local = finitegrounds[c]
        else
            f_local = finitegrounds
        end

        voltages = multiple_solver(cfg, a_local, s_local, g_local, f_local)
        local_nodemap = construct_local_node_map(nodemap, c, polymap)
        solver_called = true

        if write_v_maps && is_raster
            accum_voltages!(outvolt, voltages, local_nodemap, hbmeta)
        end
        if write_c_maps && is_raster
            accum_currents!(outcurr, voltages, cfg, a_local, voltages, f_local, local_nodemap, hbmeta)
        end

        for i in eachindex(volt)
            _i = Int(local_nodemap[i])
            _i != 0 && (volt[i] = voltages[_i])
        end
    end

    name = src == 0 ? "" : "_$(V(src))"
    if write_v_maps
        if !is_raster
            write_volt_maps(name, voltages, FullGraph(G, cellmap), flags, cfg)
        else
            write_aagrid(outvolt, name, cfg, hbmeta, cellmap, voltage = true)
        end
    end

    if write_c_maps && !write_cum_cur_map_only
        if !is_raster
            write_cur_maps(name, voltages, FullGraph(G, cellmap), finitegrounds, flags, cfg)
        else
            write_aagrid(outcurr, name, cfg, hbmeta, cellmap)
        end
    end

    if !is_raster
        v = [collect(1:size(G, 1))  voltages]
        return v, outcurr
    end

    scenario = cfg["scenario"]
    if !solver_called
        ret = Matrix{T}(undef,1,1)
        ret[1] = -1
        return ret, outcurr
    end

    if is_onetoall
        idx = findall(x->x!=0,source_map)
        val = volt[idx] ./ source_map[idx]
        if val[1] ≈ 0
            ret = Matrix{T}(undef,1,1)
            ret[1] = -1
            return ret, outcurr
        else
            ret = Matrix{T}(undef,length(val),1)
            ret[:,1] = val
            return ret, outcurr
        end
    elseif is_alltoone
        ret = Matrix{T}(undef,1,1)
        ret[1] = 0
        return ret, outcurr
    end

    return volt, outcurr
end

function multiple_solver(cfg, a::SparseMatrixCSC{T,V}, sources, grounds, finitegrounds) where {T,V}

    asolve = deepcopy(a)
    if finitegrounds[1] != -9999
        # asolve = a + spdiagm(finitegrounds, 0, size(a, 1), size(a, 1))
        asolve = a + spdiagm(0 => finitegrounds)
    end

    infgrounds = findall(x -> x == Inf, grounds)
    deleteat!(sources, infgrounds)
    dst_del = V[]
    append!(dst_del, infgrounds)
    r = collect(1:size(a, 1))
    deleteat!(r, dst_del)
    asolve = asolve[r, r]

    t1 = @elapsed M = aspreconditioner(smoothed_aggregation(asolve))
    csinfo("Time taken to construct preconditioner = $t1 seconds")
    t1 = @elapsed volt = solve_linear_system(cfg, asolve, sources, M)
    csinfo("Time taken to solve linear system = $t1 seconds")

    # Replace the inf with 0
    voltages = zeros(eltype(a), length(volt) + length(infgrounds))
    k = 1
    for i = 1:size(voltages, 1)
        if i in infgrounds
            voltages[i] = 0
        else
            #voltages[i] = volt[1][k]
            voltages[i] = volt[k]
            k += 1
        end
    end
    voltages
end

struct FullGraph{T,V}
    matrix::SparseMatrixCSC{T,V}
    cc::Vector{V}
    local_nodemap::Matrix{V}
    hbmeta::RasterMeta
    cellmap::Matrix{T}
end
FullGraph(G::SparseMatrixCSC{T,V}, cellmap) where {T,V} = FullGraph(G, collect(V, 1:size(G,1)),
                            Matrix{V}(undef,0,0), RasterMeta(), cellmap)