text:
# Based on crying baby problem (Kochenderfer 2015)
# POMDP, Robust POMDP, and Robust ρPOMDP formulations

struct Baby2POMDP <: POMDP{Symbol, Symbol, Symbol}
    r_feed::Float64
    r_hungry::Float64
    discount::Float64
end

struct Baby2RPOMDP <: RPOMDP{Symbol, Symbol, Symbol}
    r_feed::Float64
    r_hungry::Float64
    discount::Float64
    usize::Float64
end

struct Baby2IPOMDP <: IPOMDP{Symbol, Symbol, Symbol}
    r_feed::Float64
    r_hungry::Float64
    discount::Float64
    inforeward::Vector{Vector{Float64}}
end

struct Baby2RIPOMDP <: RIPOMDP{Symbol, Symbol, Symbol}
    r_feed::Float64
    r_hungry::Float64
    discount::Float64
    usize::Float64
    inforeward::Vector{Vector{Float64}}
end

const states_baby2 = [:hungry, :full]
const actions_baby2 = [:feed, :nothing]
const observations_baby2 = [:quiet, :crying]
const inforeward_baby2 = [[1.0, -1.0], [-1.0, 1.0], [0.2, 0.2]]

Baby2POMDP(r_feed, r_hungry) = Baby2POMDP(r_feed, r_hungry, 0.9)
Baby2POMDP() = Baby2POMDP(-5.0, -10.0, 0.9)

Baby2RPOMDP(r_feed, r_hungry) = Baby2RPOMDP(r_feed, r_hungry, 0.9, 0.025)
Baby2RPOMDP(err) = Baby2RPOMDP(-5.0, -10.0, 0.9, err)
Baby2RPOMDP() = Baby2RPOMDP(-5.0, -10.0, 0.9, 0.025)

Baby2IPOMDP(r_feed, r_hungry) = Baby2IPOMDP(r_feed, r_hungry, 0.9, inforeward_baby2)
Baby2IPOMDP(alphas) = Baby2IPOMDP(-5.0, -10.0, 0.9, alphas)
Baby2IPOMDP() = Baby2IPOMDP(-5.0, -10.0, 0.9, inforeward_baby2)

Baby2RIPOMDP(r_feed, r_hungry) = Baby2RIPOMDP(r_feed, r_hungry, 0.9, 0.025, inforeward_baby2)
Baby2RIPOMDP(err::Float64, alphas::Vector{Vector{Float64}}) = Baby2RIPOMDP(-5.0, -10.0, 0.9, err, alphas)
Baby2RIPOMDP() = Baby2RIPOMDP(-5.0, -10.0, 0.9, 0.025, inforeward_baby2)

states(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}) = states_baby2
actions(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}) = actions_baby2
observations(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}) = observations_baby2

n_states(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}) = 2
n_actions(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}) = 2
n_observations(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}) = 2

state_index(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}, s::Symbol) = s == :hungry ? 1 : 2
action_index(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}, a::Symbol) = a == :feed ? 1 : 2
observation_index(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}, z::Symbol) = z == :quiet ? 1 : 2
obs_index(prob::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}, z::Symbol) = observation_index(prob, z)

initial_state_distribution(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}) = SparseCat([:hungry, :full], [0.0, 1.0])
initial_belief(::Union{Baby2POMDP, Baby2RPOMDP, Baby2IPOMDP, Baby2RIPOMDP}) = [0.0, 1.0]

# Transition functions
# Nominal transition function distributions
const pϵ_baby2 = 1e-6
const tdist_baby2 = [SparseCat([:hungry, :full], [0.0, 1.0]) SparseCat([:hungry, :full], [1.0, 0.0]);
                    SparseCat([:hungry, :full], [0.0, 1.0]) SparseCat([:hungry, :full], [0.1, 0.9])]
# Nominal transition function array tarray_baby2[s,a,s'] = Pr(s' | s ,a)
const tarray_baby2 = cat(3, [0.0 1.0; 0.0 0.1], [1.0 0.0; 1.0 0.9]) # P(hungry | s, a), P(full | s, a)

# Nominal transitions
transition(prob::Union{Baby2IPOMDP, Baby2POMDP}, s::Symbol, a::Symbol) = tdist_baby2[state_index(prob, s), action_index(prob, a)]
transition(prob::Union{Baby2IPOMDP, Baby2POMDP}, s::Symbol, a::Symbol, sp::Symbol) = tarray_baby2[state_index(prob, s), action_index(prob, a), state_index(prob, sp)]

# Robust transitions
# P(z,t|s,a) = P(z|s,a,t)P(t|s,a)
# PztLower = PzLower * PtLower
# PtLower = Plower(t|s,a) = t_lower[s,a,sp]
function transition(prob::Union{Baby2RPOMDP, Baby2RIPOMDP}, s::Symbol, a::Symbol)
    si, ai = state_index(prob, s), action_index(prob, a)
    plower = max.(tarray_baby2[si,ai,:] - prob.usize, 0.0 + pϵ_baby2)
    pupper = min.(tarray_baby2[si,ai,:] + prob.usize, 1.0 - pϵ_baby2)
    plower, pupper
end
function transition(prob::Union{Baby2RPOMDP, Baby2RIPOMDP}, s::Symbol, a::Symbol, sp::Symbol)
    si, ai, spi = state_index(prob, s), action_index(prob, a), state_index(prob, sp)
    plower = max.(tarray_baby2[si, ai, spi] - prob.usize, 0.0 + pϵ_baby2)
    pupper = min.(tarray_baby2[si, ai, spi] + prob.usize, 1.0 - pϵ_baby2)
    plower, pupper
end

# Nominal observation function distributions
const odist_baby2 = [SparseCat([:quiet, :crying], [0.2, 0.8])
                    SparseCat([:quiet, :crying], [0.9, 0.1])]

# Nominal observation function array oarray_baby2[a,sp,z] = Pr(z|a,sp)
const oarray_baby2 = cat(3, [0.2 0.9; 0.2 0.9], [0.8 0.1; 0.8 0.1])

# const oarray_baby2 = cat(3, [0.3 0.7; 0.3 0.7], [0.7 0.3; 0.7 0.3])

# Nominal observations
observation(prob::Union{Baby2IPOMDP, Baby2POMDP}, a::Symbol, sp::Symbol) = odist_baby2[state_index(prob, sp)]
observation(prob::Union{Baby2IPOMDP, Baby2POMDP}, a::Symbol, sp::Symbol, z::Symbol) = oarray_baby2[action_index(prob, a), state_index(prob, sp), observation_index(prob, z)]

# Robust observations
# P(z,t|s,a) = P(z|s,a,t)P(t|s,a)
# PztLower = PzLower * PtLower
# PzLower = Plower(z|s,a,t) = o_lower[a,sp,z]
function observation(prob::Union{Baby2RPOMDP, Baby2RIPOMDP}, a::Symbol, sp::Symbol)
    ai, spi = action_index(prob, a), state_index(prob, sp)
    plower = max.(oarray_baby2[ai,spi,:] - prob.usize, 0.0 + pϵ_baby2)
    pupper = min.(oarray_baby2[ai,spi,:] + prob.usize, 1.0 - pϵ_baby2)
    plower, pupper
end
function observation(prob::Union{Baby2RPOMDP, Baby2RIPOMDP}, a::Symbol, sp::Symbol, z::Symbol)
    ai, spi, zi = action_index(prob, a), state_index(prob, sp), observation_index(prob, z)
    plower = max.(oarray_baby2[ai, spi, zi] - prob.usize, 0.0 + pϵ_baby2)
    pupper = min.(oarray_baby2[ai, spi, zi] + prob.usize, 1.0 - pϵ_baby2)
    plower, pupper
end

function dynamics(prob::Union{Baby2POMDP,Baby2IPOMDP})
    ns = n_states(prob)
    nz = n_observations(prob)
    na = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na)
    plower = Array{Float64}(ns,nz,ns,na)
    pupper = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        pt = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pz = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        p[tind, zind, sind, aind] = pz * pt
    end
    p
end

function dynamics(prob::Union{Baby2RPOMDP,Baby2RIPOMDP})
    ns = n_states(prob)
    nz = n_observations(prob)
    na = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na)
    plower = Array{Float64}(ns,nz,ns,na)
    pupper = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        ptl, ptu = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pzl, pzu = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        plower[tind, zind, sind, aind] = pzl * ptl
        pupper[tind, zind, sind, aind] = pzu * ptu
    end
    plower, pupper
end

function reward(prob::Union{Baby2POMDP, Baby2RPOMDP}, s::Symbol, a::Symbol)
    r = 0.0
    if s == :hungry
        r += prob.r_hungry
    end
    if a == :feed
        r += prob.r_feed
    end
    return r
end
reward(prob::Union{Baby2POMDP,Baby2RPOMDP}, b::Vector{Float64}, s::Symbol, a::Symbol, sp::Symbol) = reward(prob,s,a)

function reward(prob::Union{Baby2IPOMDP, Baby2RIPOMDP}, b::Vector{Float64}, a::Symbol)
    rmax = -Inf
    for α in prob.inforeward
        rmax = max(rmax, dot(α, b))
    end
    rmax
end
reward(prob::Union{Baby2IPOMDP,Baby2RIPOMDP}, b::Vector{Float64}, s::Symbol, a::Symbol, sp::Symbol) = reward(prob,b,a)

function rewardalpha(prob::Union{Baby2IPOMDP, Baby2RIPOMDP}, b::Vector{Float64}, a::Symbol)
    rmax = -Inf
    ralpha = nothing
    for α in prob.inforeward
        rnext = dot(α, b)
        if rnext > rmax
            rmax = rnext
            ralpha = copy(α)
        end
    end
    ralpha
end

discount(p::Union{Baby2POMDP,Baby2RPOMDP,Baby2IPOMDP,Baby2RIPOMDP}) = p.discount

# Simulation Functions
function initial_state(prob::Union{Baby2POMDP,Baby2IPOMDP,Baby2RPOMDP,Baby2RIPOMDP}, rng::AbstractRNG)
    d = initial_state_distribution(prob)
    return rand(rng, d)
end

function generate_sor(prob::Union{Baby2POMDP,Baby2IPOMDP}, b, s, a, rng::AbstractRNG)
    sp = rand(rng, transition(prob, s, a))
    o = rand(rng, observation(prob, a, sp))
    r = reward(prob, b, s, a, sp)
    sp, o, r
end

function generate_sor(prob::Union{Baby2RPOMDP,Baby2RIPOMDP}, b, s, a, rng::AbstractRNG)
    tdist = SparseCat(states(prob), psample(transition(prob, s, a)...))
    sp = rand(rng, tdist)
    odist = SparseCat(observations(prob), psample(observation(prob, a, sp)...))
    o = rand(rng, odist)
    r = reward(prob, b, s, a, sp)
    sp, o, r
end


text:
function encode_char(input::Char)
    current_string = ""
    value = Int(input)
    if iszero(value)
        current_string = "❤"
    else
        if value >= 200
            copies, value = divrem(value, 200)
            current_string *= "🫂" ^ copies
        end
        if value >= 50
            copies, value = divrem(value, 50)
            current_string *= "💖" ^ copies
        end
        if value >= 10
            copies, value = divrem(value, 10)
            current_string *= "✨" ^ copies
        end
        if value >= 5
            copies, value = divrem(value, 5)
            current_string *= "🥺" ^ copies
        end
        if !iszero(value)
            current_string *= "," ^ value
        end
    end

    return current_string * "👉👈"
end

encode_bottom(input::String) = join(map(encode_char, collect(input)))

function decode_byte(encoded_byte::AbstractString, emoji_to_value::Dict{String,Int})
    # Decodes a string of bottom characters into the single character they represent

    bottom_chars = emoji_to_value |> keys |> collect
    return map(char -> count(char, encoded_byte) * emoji_to_value[char], bottom_chars) |> sum |> Char
end

function decode_bottom(input::String)
    emoji_to_value = Dict{String,Int}("🫂" => 200, "💖" => 50, "✨" => 10, "🥺" => 5, "," => 1, "❤" => 0)
    byte_separator = "👉👈"
    all_bottom_chars = [emoji_to_value |> keys |> collect ; byte_separator]

    if sum(map(char -> count(char, input) * length(char), all_bottom_chars)) ≠ length(input)
        error("Input string contains a character not in the Bottom specification.")
    end

    if !endswith(input, byte_separator)
        error("The input string should end with the byte separator.")
    end

    input_to_decode = split(input, byte_separator)[1:end-1]
    decoded_chars = decode_byte.(input_to_decode, Ref(emoji_to_value))

    return join(decoded_chars)
end

function write_output(data::String, filename::String)
    open(filename, "w") do output_file
        write(output_file, data)
    end
end

function main()
    if length(ARGS) != 3
        println("Usage: $(PROGRAM_FILE) [--bottomify/--regress] input_file output_file")
        return
    end

    input_file = ARGS[2]
    output_file = ARGS[3]

    if ARGS[1] == "--bottomify"
        input_text = read(input_file, String)
        result = encode_bottom(input_text)
        write_output(result, output_file)
    elseif ARGS[1] == "--regress"
        input_text = string(readchomp(input_file))
        result = decode_bottom(input_text)
        write_output(result, output_file)
    else
        println("Invalid option specified (expected --bottomify or --regress).")
    end
end

main()

ext: jl

max_stars_repo_name: bottom-software-foundation/bottom-jl

max_issues_repo_name: bottom-software-foundation/bottom-jl

max_forks_repo_name: bottom-software-foundation/bottom-jl


text:
"""
  run_mcmc( x_::Array{Float64,1} , f , threshold::Float64 , G::Array{Float64,2} , h::Array{Float64,1} , conf::DecorrelationConfig)

runs mcmc steps which sample from the probability distribution f_density, restricted to
the space {x|f_cost(x)<threshold}

# Note
if all steps fail, then the returned value of x (xf) will be NaN
"""
function run_mcmc( x_::Array{Float64,1} , x_f_density::Float64 , f_density , f_cost , threshold::Float64 , G::Array{Float64,2} , h::Array{Float64,1} , conf::DecorrelationConfig)
    d = length(x_)

    dd_xx       = NaN * ones(d,conf.steps+1)
    dd_xx_fdens = NaN * ones(conf.steps+1)
    dd_xx_fcost = NaN * ones(conf.steps+1)

    x  = copy(x_)

    use_density = !isa(f_density,Void)

    xfdens::Float64 = NaN
    if(use_density)
        if(isnan(x_f_density))
            xfdens = f_density(x)
        else
            xfdens = x_f_density
        end
    end

    hrconf = HRStepConfig()
    tdata_hr  = HRSamplingTrackingData()
    tdata_mh  = MHSamplingTrackingData()
    steps_data = Vector{Int64}(conf.steps)

    # config:
    mh_max_tries = 8

    cnt_successes = 0
    cnt_fails = 0
    rv_all = mvnrnd(x,conf.cov_estimate,conf.steps)
    step_types_all = sample_proportional([conf.p_hr;conf.p_mh],conf.steps)

    xfcost::Float64 = NaN

    for zi in 1:conf.steps
        if(step_types_all[zi]==1)
            # HR STEP
            x_candidate_a = x - 4*rv_all[:,zi]  # times 4, because we want to be outside of viable space, if possible
            x_candidate_b = x + 4*rv_all[:,zi]
            #print("HRStep:\n")
            #print([x , x_candidate_a , x_candidate_b ])
            hr_result = hr_step( f_density , f_cost , threshold , G , h , x , x_candidate_a, x_candidate_b , hrconf ,  tdata_hr  )
            if( hr_result[1]>0 )
                cnt_successes += 1
                steps_data[zi] = +1
                x  = hr_result[2]
                # NOTE: following part will work only after upgrades to hr_step function..
                xfdens = hr_result[3]
                xfcost = hr_result[4]
                dd_xx[:,cnt_successes] = x
                dd_xx_fdens[cnt_successes] = xfdens
                dd_xx_fcost[cnt_successes] = xfcost
            else
                cnt_fails += 1
                steps_data[zi] = -1
            end
        end
        if(step_types_all[zi]==2)
            #MH STEP
            mh_result = mh_step( f_density , f_cost , threshold , G , h , x , xfdens , conf.cov_estimate , mh_max_tries , tdata_mh  )
            if(mh_result[1])
                cnt_successes += 1
                steps_data[zi] = +2
                x  = mh_result[2]
                xfdens = mh_result[3]
                xfcost = mh_result[4]
                dd_xx[:,cnt_successes] = x
                dd_xx_fdens[cnt_successes] = xfdens
                dd_xx_fcost[cnt_successes] = xfcost
            else
                cnt_fails += 1
                steps_data[zi] = -2
            end
        end
    end

    if(cnt_successes >= 1)
        # extract the successful steps data:
        dd_xx   = dd_xx[:,1:cnt_successes]
        dd_xx_fdens = dd_xx_fdens[1:cnt_successes]
        dd_xx_fcost = dd_xx_fcost[1:cnt_successes]
        # extract xv data:
        xv = dd_xx - [ x_ dd_xx[:,1:(size(dd_xx,2)-1)] ]
        # create dc:
        dd = DecorrelationData(xv , cnt_successes , cnt_fails , steps_data , dd_xx,dd_xx_fdens,dd_xx_fcost,tdata_hr,tdata_mh)
        return DecorrelationResult(x,xfdens,xfcost,dd)
    else
        print("mcmc failed..\n")
        print("step_types: \n")
        print(step_types_all)
        return DecorrelationResult(x,NaN,NaN,DecorrelationData( zeros(d,0) , cnt_successes , cnt_fails , steps_data ,zeros(d,0),[NaN],[NaN],tdata_hr,tdata_mh ))
    end
end

struct AdaptiveMCMCResult
    x::Vector{Float64}
    xfdens::Float64
    xfcost::Float64
    dd::DecorrelationData
    ecov_raw::Dict{Int,Array{Float64,2}}
    ecov_regularized::Dict{Int,Array{Float64,2}}
end

function run_adaptive_mcmc( steps::Int , x_::Array{Float64,1} , x_f_density::Float64 , f_density , f_cost , threshold::Float64 , G::Array{Float64,2} , h::Array{Float64,1} , conf_mcmc::MCMCConfig ; steps_between_cov_estim::Int = 5)

    output_level = 1

    d    = length(x_)
    xi   = copy(x_)
    xi_f = x_f_density

    # create decorrelation data vector
    dd::Vector{DecorrelationData} = Vector{DecorrelationData}()

    # create structures to save ecov info:
    data_ecov_raw = Dict{Int,Array{Float64,2}}()
    data_ecov_regularized = Dict{Int,Array{Float64,2}}()

    # Performed steps
    Xv = zeros(length(xi),0)

    # Init decorr_result struct
    sh_dcr::DecorrelationResult = DecorrelationResult()

    #
    total_steps = 0
    while(total_steps < steps)
        # 3. decorrelate
        # 3.1 estimate covariance..
        ecov_raw = eye(d)
        ecov     = NaN * ones(d,d)
        if( size(Xv,2)>4*d )
            if(size(Xv,2)>conf_mcmc.cov_estim.n_considered)
                ecov_raw = cov(Xv[ : , (end-conf_mcmc.cov_estim.n_considered+1):end ]' )
            else
                ecov_raw = cov(Xv')
            end
            sv_ecov = svd(ecov_raw)[2]
            ecov = ecov_raw + diagm(ones(d))*conf_mcmc.cov_estim.blowup_abs + diagm(ones(d)) * maximum(sv_ecov) * conf_mcmc.cov_estim.blowup_rel_maxsv + diagm(ones(d)) * minimum(sv_ecov) * conf_mcmc.cov_estim.blowup_rel_minsv

            if(output_level>=2)
                print("\necov_raw:\n$(ecov_raw)\n")
                print("ecov_regularized:\n$(ecov)\n")
                print("ecov_raw_sv: $(sv_ecov')\n")
                #print("ecov : blowup_abs: $(conf.mcmc.cov_estim.blowup_abs) , blowup_rel_max_sv: $(conf.mcmc.cov_estim.blowup_rel_maxsv*maximum(sv_ecov))) , blowup_rel_min_sv: $(conf.mcmc.cov_estim.blowup_rel_minsv*minimum(sv_ecov)) \n")
                print( @sprintf("ecov : blowup_abs: %6.3f  , blowup_rel_maxsv: %6.3f   , blowup_rel_minsv: %6.3f  \n", (conf_mcmc.cov_estim.blowup_abs) , (conf_mcmc.cov_estim.blowup_rel_maxsv*maximum(sv_ecov)) , (conf_mcmc.cov_estim.blowup_rel_minsv*minimum(sv_ecov)) ) )
            end
        else
            ecov = ecov_raw[:,:]
        end
        data_ecov_raw[total_steps]          = ecov_raw
        data_ecov_regularized[total_steps]  = ecov

        #threshold::Float64 = threshold
        steps_i            = min( steps-total_steps , steps_between_cov_estim )
        dc = DecorrelationConfig( steps_i , ecov , conf_mcmc.p_hr , conf_mcmc.p_mh )
        #sh_dcr = run_mcmc( xi , NaN, f_density , f_cost, threshold , G, h , dc )
        sh_dcr = run_mcmc( xi , xi_f , f_density , f_cost, threshold , G, h , dc )
        total_steps += sh_dcr.dd.steps_success

        # add decorrelation data
        push!(dd,sh_dcr.dd)
        # add performed steps
        if(sh_dcr.dd.steps_success>1)
            Xv  = [Xv    sh_dcr.dd.xv ]
        end

        xi   = sh_dcr.x
        xi_f = sh_dcr.xfdens
        #xi_fcost = sh_dcr.xfcost
    end

    # collapse the dds into one..
    dd_collapsed = reduce( Base.hcat , dd )

    return AdaptiveMCMCResult( sh_dcr.x , sh_dcr.xfdens , sh_dcr.xfcost , dd_collapsed , data_ecov_raw , data_ecov_regularized )
end

        #sh_dcr::Array{DecorrelationResult,1} = Array{DecorrelationResult}(size(X0,2))
        # parallel / not parallel
        #if(false)
        #    # 3.2 decorrelate samples..
        #    sh_dcr_af = Array{Future}(size(X0,2))
        #    for zi in 1:size(X0,2)
        #        sh_dcr_af[zi] = @spawn decorrelate( X0[:,zi] , f, threshold , G, h , dc )
        #    end
        #    for zi in 1:size(X0,2)
        #        sh_dcr[zi] = fetch(sh_dcr_af[zi])
        #    end
        #else
        #    for zi in 1:size(X0,2)
        #    sh_dcr[zi] = run_mcmc( X0[:,zi] , NaN, Void(), f, threshold , G, h , dc )#decorrelate( X0[:,zi] , f, threshold , G, h , dc )
        #        #sh_dcr[zi] = decorrelate( X0[:,zi] , f, threshold , G, h , dc )
        #        #if(sh_dcr[zi].)
        #    end
        #end

        #X_new   = Array{Float64}(d,0)
        #X_new_f = Array{Float64}(0)
        #Xv_new  = Array{Float64}(d,0)
        ##for zi in 1:length(sh_dcr)
        #    push!(dd,sh_dcr[zi].dd)
        #    if(sh_dcr[zi].dd.steps_success>1)
        #        X_new   = [ X_new sh_dcr[zi].x ]
        #        Xv_new  = [ Xv_new  sh_dcr[zi].dd.xv ]
        #        X_new_f = [ X_new_f ; sh_dcr[zi].xfcost ]
        #        total_steps += sh_dcr[zi].dd.steps_success
        #    end
        ##end
        #X   = [X     X_new   ]
        #Xv  = [Xv    Xv_new  ]
        #X_f = [X_f ; X_new_f ]
    #end
#end


text:
export KnightTourDrawing

"""
`KnightTourDrawing(r,c)` illustrates a Knight's tour of an r-by-c
chessboard. Throws an error if no such tour exists.
"""
function KnightTourDrawing(r::Int = 6, c::Int = 6)
    G = Knight(r, c)
    @info "Searching for a Hamiltonian cycle in a $r-by-$c Knight's move graph"
    @time h = hamiltonian_cycle(G)
    @info "Finished"

    if length(h) == 0
        throw(error("There is no such Knight's move tour"))
    end

    T = eltype(G)
    H = SimpleGraph{T}()
    for v in G.V
        add!(H, v)
    end

    for k = 1:NV(H)-1
        add!(H, h[k], h[k+1])
    end

    add!(H, h[1], h[end])

    xy = getxy(H)

    for v in H.V
        xy[v] = collect(v)
    end

    embed(H, xy)

    draw(H)

    for a = 0:c
        plot!([0.5, r + 0.5], [a + 0.5, a + 0.5], color = "black", linestyle = :dash)
    end

    for b = 0:r
        plot!([b + 0.5, b + 0.5], [0.5, c + 0.5], color = "black", linestyle = :dash)
    end

    finish()
end


text:
"""
    Uniform( xmin, xmax, nx )

Uniform nodes

"""
struct Uniform <: NodesType

    nx    :: Int64
    xmin  :: Float64
    xmax  :: Float64
    xk    :: Vector{Float64}

    function Uniform( xmin, xmax, nx )

        xk  = zeros(Float64, nx)
        xk .= range(-1, stop=1, length=nx);
        xk .= ((xk .+ 1) * (xmax - xmin)) /2 .+ xmin
        new( nx, xmin, xmax, xk )

    end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BlockWeighting(sides)
    BlockWeighting(side₁, side₂, ...)

A weighting method that assigns weights to points in spatial object
based on blocks of given `sides`. The number `n` of points inside a
block determines the weights `1/n` of these points.
"""
struct BlockWeighting{Dim,T} <: WeightingMethod
  sides::SVector{Dim,T}
end

BlockWeighting(sides::NTuple) = BlockWeighting(SVector(sides))
BlockWeighting(sides::Vararg) = BlockWeighting(SVector(sides))

weight(object, method::BlockWeighting) =
  weight(domain(object), method)

function weight(domain::Domain, method::BlockWeighting)
  p = partition(domain, BlockPartition(method.sides))

  weights = Vector{Float64}(undef, nelements(domain))
  for s in subsets(p)
    weights[s] .= 1 / length(s)
  end

  GeoWeights(domain, weights)
end


text:
using Test
using HDF5
using DSGE

path = dirname(@__FILE__)

m = AnSchorfheide()
Γ0, Γ1, C, Ψ, Π = eqcond(m)
stake = 1 + 1e-6
#G1, C, impact, fmat, fwt, ywt, gev, eu, loose = gensys(Γ0, Γ1, C, Ψ, Π, stake)
G1, C, impact, eu = gensys(Γ0, Γ1, C, Ψ, Π, stake)

file = "$path/../reference/gensys.h5"
G1_exp = h5read(file, "G1_gensys")
C_exp = h5read(file,"C_gensys")
impact_exp = h5read(file, "impact")
eu_exp = h5read(file, "eu")

@testset "Check gensys outputs match reference" begin
    @test @test_matrix_approx_eq G1_exp G1
    @test @test_matrix_approx_eq C_exp C
    @test @test_matrix_approx_eq impact_exp impact

    @test isequal(eu_exp, eu)
end

G1, C, impact, eu = gensys(Γ0, Γ1, C, Ψ, Π)

# trigger indeterminacy error
m = Model1002()
Γ0, Γ1, C, Ψ, Π = eqcond(m)
Γ1[m.equilibrium_conditions[:eq_z], m.endogenous_states[:y_t]] = 1.
Γ1[m.equilibrium_conditions[:eq_z], m.endogenous_states[:y_f_t]] = 1.

@testset "Check indeterminacy zeros warning" begin
    @test_logs (:warn, "Indeterminacy: 1 loose endogenous error(s)") gensys(Γ0, Γ1, C, Ψ, Π)
end

# trigger coincident zeros
m = Model1002()
Γ0, Γ1, C, Ψ, Π = eqcond(m)
Γ1[1,:] .= 0.
Γ0[1,:] .= 0.
Ψ[1,:] .= 0.
Π[1,:] .= 0.
@testset "Check output coincident zeros warning" begin
    @test_logs (:warn, "Coincident zeros. Indeterminacy and/or nonexistence.") gensys(Γ0, Γ1, C, Ψ, Π)
    @info "After this message, a warning will be intentionally printed."
    G1, tmpC, impact, eu = gensys(Γ0, Γ1, C, Ψ, Π)
    @test eu == vec([-2, -2])
    @test isempty(G1)
    @test isempty(tmpC)
    @test isempty(impact)
end

Γ1[1,1] = 1e3
@testset "Check nonexistence warning" begin
    @test_logs (:warn, "Nonexistence: number of unstable roots exceeds number of jump variables") (:warn, "Indeterminacy: 1 loose endogenous error(s)") gensys(Γ0, Γ1, C, Ψ, Π)
end

Random.seed!(1793)
F = schur!(complex(rand(3,3)), complex(rand(3,3)))
F.T[1,1] = abs(F.S[1,1]) * (1 + 1e-3)
Γ0, Γ1, ~, ~, ~ = eqcond(m)
Fm = schur!(complex(Γ0), complex(Γ1))
@testset "Check new_div function" begin
    @test DSGE.new_div(F) == (1+5e-4)
    @test DSGE.new_div(Fm) == 1.01
end


text:
using Pkg
Pkg.activate(".")

# ] instantiate

using Revise, StableDQMC, GenericSVD, JacobiSVD, LinearAlgebra, Plots, Statistics, LaTeXStrings, DataStructures
using BenchmarkTools, StatsPlots, DataFrames

gr();
#PyPlot.plt.style.use("publication_tex")
default(
    framestyle = :grid,
    grid = true,
    guidefont=font(15),
#     framestyle = :box,
#     grid = false,
#     guidefont=font(18),
    legend = true,
    size = (1.6*360, 360),
    xtickfont=font(13),
    ytickfont=font(13),
    legendfont=font(12),
    linewidth = 1.5
);

# default(
#     framestyle = :box,
#     legend = true,
#     grid = false,
#     size = (1.5*250, 250),
#     xtickfont=font(8),
#     ytickfont=font(8),
#     guidefont=font(11),
#     legendfont=font(8),
#     linewidth = 1
# );

# Some basic tests
B = rand(5,5) # regular matrix
N = 10
@assert isapprox(calc_Bchain(B, N)[1], B^N)
@assert isapprox(Matrix(calc_Bchain_svd(B, N; svdalg = gesdd)[1]), B^N)
@assert isapprox(Matrix(calc_Bchain_svd(B, N; svdalg =  gesvd)[1]), B^N)
@assert isapprox(Matrix(calc_Bchain_svd(B, N; svdalg =  genericsvd)[1]), B^N)
@assert isapprox(Matrix(calc_Bchain_qr(B, N)[1]), B^N)
# @assert isapprox(calc_Bchain_stabilized(B, N, gesvj)[1], B^N)

n = 8
t = -1
μ = -0.1
Δτ = 0.1
T = diagm(-1 => fill(-t, n-1), 1 => fill(-t, n-1), 0 => fill(μ, n))

# PBC
T[1,end] = -t
T[end,1] = -t

# B from spin-fermion model
#B = [0.975173 0.0971936 0.192475 0.0191836 0.0 0.0 0.0 0.0; 0.0971936 0.975173 0.0191836 0.192475 0.0 0.0 0.0 0.0; 0.192475 0.0191836 0.975173 0.0971936 0.0 0.0 0.0 0.0; 0.0191836 0.192475 0.0971936 0.975173 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.975173 -0.192475 -0.0971936 0.0191836; 0.0 0.0 0.0 0.0 -0.192475 0.975173 0.0191836 -0.0971936; 0.0 0.0 0.0 0.0 -0.0971936 0.0191836 0.975173 -0.192475; 0.0 0.0 0.0 0.0 0.0191836 -0.0971936 -0.192475 0.975173]

B = exp(-Δτ .* T)

# Compare to exact discretized energies
ks = LinRange(-pi, pi, n+1)[1:end-1]
ϵ(k) = -2*t*cos(k) + μ
sort(ϵ.(ks))

plot(ks, sort(ϵ.(ks)), marker=true, label="exact")
plot!(ks, eigvals(T), marker=true, label="eigenvalues")

cond(B)

kappa(N) = exp(4abs(t)*N*Δτ)

cond(B)

kappa(1)

estimate_breakdown(sigdigits) = sigdigits * log(10) / (4 * abs(t) * Δτ)

estimate_breakdown(17)

Bbig = BigFloat.(B)
N = 500

setprecision(100) do
    global svs_genericsvd = calc_Bchain_svd(Bbig, N; svdalg = genericsvd)[2];
end
svs_qr_udt = calc_Bchain_qr(B, N)[2];
svs_gesvd  = calc_Bchain_svd(B, N; svdalg = gesvd)[2];
svs_gesdd  = calc_Bchain_svd(B, N; svdalg = gesdd)[2];
svs_gesvj  = calc_Bchain_svd(B, N; svdalg = gesvj)[2];

cond(exp.(svs_qr_udt[end,:]))

kappa(N)

betas = range(1, N * 0.1, length=N)
plot(betas, svs_genericsvd, color="#CA0020", labels = reshape(["exact", "", "","", "","", "", ""], (1,8)))
plot!(betas, svs_qr_udt, color="#E66101", labels = reshape(["QR", "", "","", "","", "", ""], (1,8)))
plot!(betas, svs_gesvj, color="#FDB863", labels = reshape(["SVD (Jacobi)", "", "","", "","", "", ""], (1,8)))
plot!(betas, svs_gesvd, color="#008837", labels = reshape(["SVD", "", "","", "","", "", ""], (1,8)))
plot!(betas, svs_gesdd, color="#5E3C99", labels = reshape(["SVD (D\\&C)", "", "","", "","", "", ""], (1,8)))

ylabel!("log singular values")
xlabel!(L"inverse temperature $\beta$")
ylims!(-160, 100)
xlims!(0., 50.)
savefig("../paper/figures/decomp_comparison_simple_mkl.pdf")

svs_naive = calc_Bchain(B, N)[2];

betas = range(1, N * 0.1, length=N)
# plot(svs_genericsvd, color="#CA0020", labels = reshape(["exact", "", "","", "","", "", ""], (1,8)))
plot(betas, svs_qr_udt, color="#E66101", labels = reshape(["stable", "", "","", "","", "", ""], (1,8)))
plot!(betas, svs_naive, color=:green, labels = reshape(["naive", "", "","", "","", "", ""], (1,8)))
plot!(1:40, fill(-17, 40), color=:darkgray, linestyle=:dash, label = "Float64")

ylabel!("log singular values")
xlabel!(L"inverse temperature $\beta$")
ylims!(-80, 80)
xlims!((0., 40.))
savefig("../paper/figures/naive_vs_stable_mkl.pdf")

using DoubleFloats
Bdb64 = Double64.(B)
svs_naive = calc_Bchain(Bdb64, N)[2];

betas = range(1, N * 0.1, length=N)
# plot(svs_genericsvd, color="#CA0020", labels = reshape(["exact", "", "","", "","", "", ""], (1,8)))
plot(betas, svs_qr_udt, color="#E66101", labels = reshape(["stable", "", "","", "","", "", ""], (1,8)))
plot!(betas, svs_naive, color=:green, labels = reshape(["naive", "", "","", "","", "", ""], (1,8)))
plot!(1:40, fill(-34, 40), color=:darkgray, linestyle=:dash, label = "Double64")

ylabel!("log singular values")
xlabel!(L"inverse temperature $\beta$")
ylims!(-80, 80)
xlims!((0., 40.))
savefig("../paper/figures/naive_vs_stable_db64_mkl.pdf")

function accuracy_greens(B; svdinversion = inv_one_plus)
    Bbig = BigFloat.(B)

    Δs_qr    = SortedDict{Int, Float64}()
    Δs_gesvd = SortedDict{Int, Float64}()
    Δs_gesdd = SortedDict{Int, Float64}()
    Δs_gesvj = SortedDict{Int, Float64}()

    for beta in range(5, 40, step=5)
        N = Int(beta * 10)

        local G_genericsvd
        setprecision(1000) do
            F_genericsvd = calc_Bchain_svd(Bbig, N; svdalg = genericsvd)[1];
            G_genericsvd = inv(I + Matrix(F_genericsvd))
        end

        # Bchain
        F_qr_udt = calc_Bchain_qr(B, N)[1];
        F_gesvd  = calc_Bchain_svd(B, N; svdalg = gesvd)[1];
        F_gesdd  = calc_Bchain_svd(B, N; svdalg = gesdd)[1];
        F_gesvj  = calc_Bchain_svd(B, N; svdalg = gesvj)[1];

        # inversion
        G_qr_udt = inv_one_plus(F_qr_udt)
        G_gesvd  = svdinversion(F_gesvd)
        G_gesdd  = svdinversion(F_gesdd)
        G_gesvj  = svdinversion(F_gesvj)

        Δs_qr[beta]    = maximum(abs.(G_qr_udt - G_genericsvd))
        Δs_gesvd[beta] = maximum(abs.(G_gesvd - G_genericsvd))
        Δs_gesdd[beta] = maximum(abs.(G_gesdd - G_genericsvd))
        Δs_gesvj[beta] = maximum(abs.(G_gesvj - G_genericsvd))
    end

    return Δs_qr, Δs_gesvd, Δs_gesdd, Δs_gesvj
end

Δs_qr, Δs_gesvd, Δs_gesdd, Δs_gesvj = accuracy_greens(B);

betas = collect(keys(Δs_qr))
p = plot(betas, log.(10, values(Δs_qr)), color="#E66101", marker=true, label="QR")
plot!(betas, log.(10, values(Δs_gesvd)), color="#008837", marker=true, label="SVD")
plot!(betas, log.(10, values(Δs_gesdd)), color="#5E3C99", marker=true, label="SVD (D\\&C)")
plot!(betas, log.(10, values(Δs_gesvj)), color="#FDB863", marker=true, label="SVD (Jacobi)")

xlabel!(L"inverse temperature $\beta$")
ylabel!(L"\log(\textrm{max}(\textrm{abs}(G - G_{\textrm{exact}})))")
# ylims!(0., 1e-3)
savefig("../paper/figures/accuracy_greens_svd_regularinv_mkl.pdf")

Δs_qr, Δs_gesvd, Δs_gesdd, Δs_gesvj = accuracy_greens(B, svdinversion = inv_one_plus_loh);

betas = collect(keys(Δs_qr))
p = plot(betas, log.(10, values(Δs_qr)), color="#E66101", marker=true, label="QR")
plot!(betas, log.(10, values(Δs_gesvd)), color="#008837", marker=true, label="SVD")
plot!(betas, log.(10, values(Δs_gesdd)), color="#5E3C99", marker=true, label="SVD (D\\&C)")
plot!(betas, log.(10, values(Δs_gesvj)), color="#FDB863", marker=true, label="SVD (Jacobi)")

xlabel!(L"inverse temperature $\beta$")
ylabel!(L"\log(\textrm{max}(\textrm{abs}(G - G_{\textrm{exact}})))")
savefig("../paper/figures/accuracy_greens_svd_loh_mkl.pdf")

absdiff(x,y) = abs(x-y)
reldiff(x, y) = abs((x - y) / (x + y))

function accuracy_determinant(B; svdinversion = svd_inv_one_plus)
    Bbig = BigFloat.(B)

    Δs_abs_qr    = SortedDict{Int, Float64}()
    Δs_abs_gesvd = SortedDict{Int, Float64}()
    Δs_abs_gesdd = SortedDict{Int, Float64}()
    Δs_abs_gesvj = SortedDict{Int, Float64}()

    Δs_rel_qr    = SortedDict{Int, Float64}()
    Δs_rel_gesvd = SortedDict{Int, Float64}()
    Δs_rel_gesdd = SortedDict{Int, Float64}()
    Δs_rel_gesvj = SortedDict{Int, Float64}()

    for beta in range(5, 40, step=5)
        N = Int(beta * 10)

        local det_genericsvd
        setprecision(1000) do
            F_genericsvd = calc_Bchain_svd(Bbig, N; svdalg = genericsvd)[1];
            det_genericsvd = prod(svd_inv_one_plus(F_genericsvd).S)
        end

        # Bchain
        F_qr_udt   = calc_Bchain_qr(B, N)[1];
        F_gesvd    = calc_Bchain_svd(B, N; svdalg = gesvd)[1];
        F_gesdd    = calc_Bchain_svd(B, N; svdalg = gesdd)[1];
        F_gesvj    = calc_Bchain_svd(B, N; svdalg = gesvj)[1];

        # inversion + det extraction
        det_qr_udt = prod(udt_inv_one_plus(F_qr_udt).D)
        det_gesvd  = prod(svdinversion(F_gesvd; svdalg = gesvd!).S)
        det_gesdd  = prod(svdinversion(F_gesdd; svdalg = gesdd!).S)
        det_gesvj  = prod(svdinversion(F_gesvj; svdalg = StableDQMC.gesvj!).S)

        Δs_abs_qr[beta]    = absdiff(det_qr_udt, det_genericsvd)
        Δs_abs_gesvd[beta] = absdiff(det_gesvd, det_genericsvd)
        Δs_abs_gesdd[beta] = absdiff(det_gesdd, det_genericsvd)
        Δs_abs_gesvj[beta] = absdiff(det_gesvj, det_genericsvd)

        Δs_rel_qr[beta]    = reldiff(det_qr_udt, det_genericsvd)
        Δs_rel_gesvd[beta] = reldiff(det_gesvd, det_genericsvd)
        Δs_rel_gesdd[beta] = reldiff(det_gesdd, det_genericsvd)
        Δs_rel_gesvj[beta] = reldiff(det_gesvj, det_genericsvd)
    end

    return Δs_abs_qr, Δs_abs_gesvd, Δs_abs_gesdd, Δs_abs_gesvj,
            Δs_rel_qr, Δs_rel_gesvd, Δs_rel_gesdd, Δs_rel_gesvj
end

Δs_abs_qr, Δs_abs_gesvd, Δs_abs_gesdd, Δs_abs_gesvj,
Δs_rel_qr, Δs_rel_gesvd, Δs_rel_gesdd, Δs_rel_gesvj = accuracy_determinant(B, svdinversion = svd_inv_one_plus)

betas = collect(keys(Δs_rel_qr))

p = plot(betas, log.(10, values(Δs_rel_qr)), color="#E66101", marker=true, label="QR")
plot!(betas, log.(10, values(Δs_rel_gesvd)), color="#008837", marker=true, label="SVD")
plot!(betas, log.(10, values(Δs_rel_gesdd)), color="#5E3C99", marker=true, label="SVD (D\\&C)")
plot!(betas, log.(10, values(Δs_rel_gesvj)), color="#FDB863", marker=true, label="SVD (Jacobi)")

xlabel!(L"inverse temperature $\beta$")
ylabel!(L"\log(\textrm{reldiff}(\textrm{det}G - \textrm{det}G_{\textrm{exact}}))")
savefig("../paper/figures/accuracy_det_svd_regularinv_mkl.pdf")

Δs_abs_qr, Δs_abs_gesvd, Δs_abs_gesdd, Δs_abs_gesvj,
Δs_rel_qr, Δs_rel_gesvd, Δs_rel_gesdd, Δs_rel_gesvj = accuracy_determinant(B, svdinversion = svd_inv_one_plus_loh)

betas = collect(keys(Δs_rel_qr))
p = plot(betas, log.(10, values(Δs_rel_qr)), color="#E66101", marker=true, label="QR")
plot!(betas, log.(10, values(Δs_rel_gesvd)), color="#008837", marker=true, label="SVD")
plot!(betas, log.(10, values(Δs_rel_gesdd)), color="#5E3C99", marker=true, label="SVD (D\\&C)")
plot!(betas, log.(10, values(Δs_rel_gesvj)), color="#FDB863", marker=true, label="SVD (Jacobi)")

xlabel!(L"inverse temperature $\beta$")
ylabel!(L"\log(\textrm{reldiff}(\textrm{det}G - \textrm{det}G_{\textrm{exact}}))")
savefig("../paper/figures/accuracy_det_svd_loh_mkl.pdf")

function benchmark_decomps(; evals=3)
    Ns = Int.(range(100, 500, step=100))
    df = DataFrame(N=Int[], qr=Float64[], gesdd=Float64[], gesvd=Float64[], gesvj=Float64[])

    for N in Ns
        x = rand(ComplexF64, N,N)
        t_qr = @belapsed udt($x) evals=evals
        t_gesdd = @belapsed gesdd($x) evals=evals
        t_gesvd = @belapsed gesvd($x) evals=evals
        t_gesvj = @belapsed gesvj($x) evals=1
        push!(df, [N, t_qr, t_gesdd, t_gesvd, t_gesvj])
    end

    return df
end

#df_decomps = benchmark_decomps()

#@df df_decomps plot(:N, [:gesvd ./ :qr, :gesdd ./ :qr, :gesvj ./ :qr],
#            marker=true,
#            color=permutedims(["#008837", "#5E3C99", "#FDB863"]),
#            label=permutedims(["SVD", "SVD (D\\&C)", "SVD (Jacobi)"]),
##             yscale=:log10
#            ylims=(0,29)
#           )
#@df df_decomps plot!(:N, fill(1, length(:N)), color=:grey, linestyle=:dash, label="QR")
#xlabel!(L"N of $N\times N$ square matrix")
#ylabel!("slowdown factor rel. to QR")
#savefig("../paper/figures/benchmark_decomps_mkl.pdf")

function benchmark_inv_one_plus(; evals=3)
    Ns = Int.(range(100, 500, step=100))
    df = DataFrame(N=Int[], qr=Float64[], gesdd=Float64[], gesvd=Float64[], gesvj=Float64[])

    for N in Ns
        x = rand(ComplexF64, N,N)
        u = udt(x)
        s = svd(x)

        t_qr = @belapsed udt_inv_one_plus($u) evals=evals
        t_gesdd = @belapsed svd_inv_one_plus($s, svdalg = $gesdd!) evals=evals
        t_gesvd = @belapsed svd_inv_one_plus($s, svdalg = $gesvd!) evals=evals
        t_gesvj = @belapsed svd_inv_one_plus($s, svdalg = $StableDQMC.gesvj!) evals=evals
        push!(df, [N, t_qr, t_gesdd, t_gesvd, t_gesvj])
    end

    return df
end

#df_greens = benchmark_inv_one_plus(evals=1)

#@df df_greens plot(:N, [:gesvd ./ :qr, :gesdd ./ :qr, :gesvj ./ :qr],
#            marker=true,
#            color=permutedims(["#008837", "#5E3C99", "#FDB863"]), # #E66101 for QR
#            label=permutedims(["SVD", "SVD (D\\&C)", "SVD (Jacobi)"]),
##             yscale=:log10
##             ylims=(0,50)
#           )
#@df df_greens plot!(:N, fill(1, length(:N)), color=:grey, linestyle=:dash, label="QR")
#xlabel!(L"number of sites $N$")
#ylabel!("slowdown factor rel. to QR")
#savefig("../paper/figures/benchmark_fact_inv_one_plus_mkl.pdf")

function benchmark_inv_one_plus_loh(; evals=3)
    Ns = Int.(range(100, 500, step=100))
    df = DataFrame(N=Int[], qr=Float64[], gesdd=Float64[], gesvd=Float64[], gesvj=Float64[])

    for N in Ns
        x = rand(ComplexF64, N,N)
        u = udt(x)
        s = svd(x)

        t_qr = @belapsed udt_inv_one_plus_loh($u) evals=evals
        t_gesdd = @belapsed svd_inv_one_plus_loh($s, svdalg = $gesdd!) evals=evals
        t_gesvd = @belapsed svd_inv_one_plus_loh($s, svdalg = $gesvd!) evals=evals
        t_gesvj = @belapsed svd_inv_one_plus_loh($s, svdalg = $StableDQMC.gesvj!) evals=evals
        push!(df, [N, t_qr, t_gesdd, t_gesvd, t_gesvj])
    end

    return df
end

#df_loh = benchmark_inv_one_plus_loh()

#@df df_loh plot(:N, [:gesvd ./ :qr, :gesdd ./ :qr, :gesvj ./ :qr],
#            marker=true,
#            color=permutedims(["#008837", "#5E3C99", "#FDB863"]), # #E66101 for QR
#            label=permutedims(["SVD", "SVD (D\\&C)", "SVD (Jacobi)"]),
##             yscale=:log10
##             ylims=(0,50)
#           )
#@df df_loh plot!(:N, fill(1, length(:N)), color=:grey, linestyle=:dash, label="QR")
#xlabel!(L"number of sites $N$")
#ylabel!("slowdown factor rel. to QR")
#savefig("../paper/figures/benchmark_fact_inv_one_plus_loh_mkl.pdf")

# G(slice, 1) = [B(slice, 1)^-1 + B(beta, slice)]^-1
function accuracy_tdgf(B; loh = false)
    Bbig = BigFloat.(B)

#     betas = range(5, 40, step=5)
    beta = 40
    M = Int(beta * 10)
    slices = 1:10:M

    df = DataFrame(slice=Int[], qr=Float64[], gesdd=Float64[], gesvd=Float64[], gesvj=Float64[])

    inv_sum_method = loh ? inv_sum_loh : inv_sum

    for s in slices
        local TDGF_genericsvd
        setprecision(1000) do
            TDGF_genericsvd = calc_tdgf_svd(Bbig, s, M-s, inv_sum_method = inv_sum_method, svdalg = genericsvd!)
        end

        TDGF_qr = calc_tdgf_qr(B, s, M-s, inv_sum_method = inv_sum_method)
        TDGF_gesdd = calc_tdgf_svd(B, s, M-s, inv_sum_method = inv_sum_method, svdalg = gesdd!)
        TDGF_gesvd = calc_tdgf_svd(B, s, M-s, inv_sum_method = inv_sum_method, svdalg = gesvd!)
        TDGF_gesvj = calc_tdgf_svd(B, s, M-s, inv_sum_method = inv_sum_method, svdalg = StableDQMC.gesvj!)

        Δ_qr    = mean(abs, TDGF_qr - TDGF_genericsvd)
        Δ_gesdd = mean(abs, TDGF_gesdd - TDGF_genericsvd)
        Δ_gesvd = mean(abs, TDGF_gesvd - TDGF_genericsvd)
        Δ_gesvj = mean(abs, TDGF_gesvj - TDGF_genericsvd)

        push!(df, [s, Δ_qr, Δ_gesdd, Δ_gesvd, Δ_gesvj])
    end

    return df
end

df_tdgf = accuracy_tdgf(B; loh = false)

@df df_tdgf plot(:slice .* 0.1, (x -> log.(10, x)).([:qr, :gesvd, :gesdd, :gesvj]),
            marker=true,
            color=permutedims(["#E66101", "#008837", "#5E3C99", "#FDB863"]), # #E66101 for QR
            label=permutedims(["QR", "SVD", "SVD (D\\&C)", "SVD (Jacobi)"]),
#             xlims=(0,10)
           )
xlabel!(L"imaginary time $\tau$")
ylabel!(L"\log(\textrm{max}(\textrm{abs}(G(\tau, 0) - G_{\textrm{exact}}(\tau, 0))))")
savefig("../paper/figures/accuracy_tdgf_mkl.pdf")

df_tdgf_loh = accuracy_tdgf(B; loh = true)

@df df_tdgf_loh plot(:slice .* 0.1, (x -> log.(10, x)).([:qr, :gesvd, :gesdd, :gesvj]),
            marker=true,
            color=permutedims(["#E66101", "#008837", "#5E3C99", "#FDB863"]), # #E66101 for QR
            label=permutedims(["QR", "SVD", "SVD (D\\&C)", "SVD (Jacobi)"]),
#             xlims=(0,10)
           )
xlabel!(L"imaginary time $\tau$")
ylabel!(L"\log(\textrm{max}(\textrm{abs}(G(\tau, 0) - G_{\textrm{exact}}(\tau, 0))))")
savefig("../paper/figures/accuracy_tdgf_loh_mkl.pdf")


text:
# Abstract Code
# =============

# The existence of .codeword::Vector{DNAKmer{k}} is assumed.
abstract AbstractCode{k} <: AbstractVector{DNAKmer{k}}

function Base.getindex(code::AbstractCode, i::Integer)
    return code.codewords[i]
end

function Base.size(code::AbstractCode)
    return size(code.codewords)
end

function Base.filter{k}(f::Function, code::AbstractCode{k})
    return filter!(f, deepcopy(code))
end

function Base.filter!{k}(f::Function, code::AbstractCode{k})
    filter!(f, code.codewords)
    return code
end