text:
using DifferentialEquations
using DataFrames
using CSV
using LaTeXStrings
using Statistics
using GLM
using StatsPlots
using Optim
using ForwardDiff
using ColorSchemes
using Random, Distributions
using Plots, ColorSchemes
using GLM, Measures
using LightGraphs, SimpleWeightedGraphs#, GraphPlot
using Compose, Cairo, Fontconfig
using LsqFit, Dates
import Gadfly
g = Gadfly
#
pyplot(grid=false)
#
fnt = Plots.font("Arial")
default(titlefont=fnt, guidefont=fnt, tickfont=fnt, legendfont=fnt)

crec = parse(Colorant, "#F5EBCE")
crel = parse(Colorant, "#DEE6F0")

################################################################################
# Model:
# A <-> R <-> I
# ka: R -> A
# kr: A -> R
# ki: R -> I
# kp: I -> R
################################################################################

function three_state_model(dc, c, k, t)
    dc[1] = k['a'] * c[2] - k['s'] * c[1]
    dc[2] = k['s'] * c[1] - k['a'] * c[2] - k['i'] * c[2] + k['p'] * c[3]
    dc[3] = k['i'] * c[2] - k['p'] * c[3]
end

function isoutofdomain(c, k, t)
    if (c[1]<=1) && (c[2]<=1) && (c[3]<=1) && (c[1]>=0) && (c[2]>=0) && (c[3]>=0)
        return false
    else
        return true
    end
end

function run_recruitment(k_time_pairs, get_irrev=true, c0=[1.0; 0.0; 0.0])
    t = fill(Float64, 0)
    u = fill(Float64[], 0)

    for pair in k_time_pairs
        k = pair[1]
        times = pair[2]
        prob = ODEProblem(three_state_model, c0, times, k)
        sol = solve(prob, isoutofdomain=isoutofdomain, saveat=0.001)
        t = [t; sol.t]
        u = [u; sol.u]
        c0 = sol.u[end]
    end
    as = [x[1] for x in u]
    rs = [x[2] for x in u]
    is = [x[3] for x in u]
    if get_irrev
        data = [as rs is]
    else
        data = [as rs]
    end

    return [t, data]
end

function plot_recruitment(k_time_pairs, grey_pairs, blue_pairs, filename)
    sol = run_recruitment(k_time_pairs)

    p1 = plot()
    for pair in grey_pairs
        p1 = vspan!(p1, pair, label="", color=crec)
    end
    for pair in blue_pairs
        p1 = vspan!(p1, pair, label="", color=crel)
    end
    p1 = plot!(p1,
               sol[1], sol[2],
               label=["Active (A)" "Reversibly Silent (R)" "Irreversibly Silent (I)"],
               color = [:blue :orange :red],
               xlabel="Days",
               ylabel="Fraction of Cells",
               yticks = 0:0.2:1.0,
               linewidth=4,
               legendfontsize=13,
               size=(600,300))
    savefig(p1, filename)
end

################################################################################
### FIGURE 3A
################################################################################

k_recruit = Dict('a' => 0.1,
                 's' => 1,
                 'i' => 0.2,
                 'p' => 0)

k_release = Dict('a' => 0.1,
                 's' => 0,
                 'i' => 0,
                 'p' => 0)
kt_pairs = [[k_recruit, (0.0, 5.0)],
            [k_release, (5.0, 15.0)]]
gp = [[0.0, 5.0]]
bp = [[5.0, 15.0]]
fname = "fig_3a.pdf"

plot_recruitment(kt_pairs, gp, bp, fname)

################################################################################
### FIGURE 3B-C
################################################################################

function get_pulsed_silencing_curve(s=200)
    pulse_pairs = [[k_recruit, (0.0, 1.0)],
                   [k_release, (1.0, 2.0)],
                   [k_recruit, (2.0, 3.0)],
                   [k_release, (3.0, 4.0)],
                   [k_recruit, (4.0, 5.0)],
                   [k_release, (5.0, 6.0)],
                   [k_recruit, (6.0, 7.0)],
                   [k_release, (7.0, 8.0)],
                   [k_recruit, (8.0, 9.0)],
                   [k_release, (9.0, 10.0)],
                   [k_release, (10.0, 15.0)]
                   ]
    gpp = [[0.0, 1.0], [2.0, 3.0], [4.0, 5.0], [6.0, 7.0], [8.0, 9.0]]
    bpp = [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0], [7.0, 8.0], [9.0, 15.0]]

    psol = run_recruitment(pulse_pairs)

    pp = plot()
    for pair in gpp
        ymin = 0.005
        ymax = 1.05
        pp = plot!(pp,
                   [pair[1], pair[1], pair[2], pair[2]],
                   [ymin, ymax, ymax, ymin],
                   color = crec,
                   linealpha = 0,
                   linewidth = 0,
                   label = "",
                   seriestype = :shape)
        #  pp = vspan!(pp, pair, label="", linecolor=crec, fillcolor=crec, linealpha=0, fillalpha=0)
    end
    for pair in bpp
        ymin = 0.005
        ymax = 1.05
        pp = plot!(pp,
                   [pair[1], pair[1], pair[2], pair[2]],
                   [ymin, ymax, ymax, ymin],
                   color = crel,
                   linealpha = 0,
                   linewidth = 0,
                   label = "",
                   seriestype = :shape)
        #  pp = vspan!(pp, pair, label="", linecolor=crel, fillcolor=crel, linealpha=0)
    end
    #  for pair in gpp
    #      pp = vspan!(pp, pair, label="", color=crec)
    #  end
    #  for pair in bpp
    #      pp = vspan!(pp, pair, label="", color=crel)
    #  end

    pp = plot!(pp,
               psol[1], psol[2],
               label=["A" "R" "I"],
               color = [:blue :orange :red],
               xlabel="Time (days)",
               ylabel="Fraction of Cells",
               yticks = 0:0.5:1.0,
               xticks = 0:5:15,
               #  xlims = (-Inf, 15),
               #  ylims=(-Inf, 1.0),
               linewidth=4,
               legend = :outerright,
               order = 0,
               size=(500, s))
    return pp
end

function get_pulsed_vs_continuous(s=200)
    pulse_pairs = [[k_recruit, (0.0, 1.0)],
                   [k_release, (1.0, 2.0)],
                   [k_recruit, (2.0, 3.0)],
                   [k_release, (3.0, 4.0)],
                   [k_recruit, (4.0, 5.0)],
                   [k_release, (5.0, 6.0)],
                   [k_recruit, (6.0, 7.0)],
                   [k_release, (7.0, 8.0)],
                   [k_recruit, (8.0, 9.0)],
                   [k_release, (9.0, 10.0)],
                   [k_release, (10.0, 15.0)]
                   ]
    gpp = [[0.0, 1.0], [2.0, 3.0], [4.0, 5.0], [6.0, 7.0], [8.0, 9.0]]

    cont_pairs = [[k_release, (0.0, 5.0)],
                  [k_recruit, (5.0, 10.0)],
                  [k_release, (10.0, 15.0)]]
    gpc = [[5.0, 10.0]]

    psol = run_recruitment(pulse_pairs)
    csol = run_recruitment(cont_pairs)

    pc = plot(
              psol[1], psol[2][:, 3],
              label="Pulsed",
              color=:purple,
              xlabel="Time (days)",
              ylabel="I(t)",
              ylims=(0.0, 1.0),
              yticks = 0:0.5:1.0,
              linewidth=3,
              legend = :topright,
              size=(500, 300))
    pc = plot!(pc,
               csol[1], csol[2][:, 3],
               label="Continuous",
               color=:deepskyblue,
               linewidth=3)

    for pair in gpp
        pc = vspan!(pc, pair, label="", color=:violet, alpha=0.2)
    end
    for pair in gpc
        pc = vspan!(pc, pair, label="", color=:cyan, alpha=0.2)
    end
    return pc
end

function plot_pulsed_comparison(s=200)
    l = @layout [a ; b ; c]

    pp = get_pulsed_silencing_curve(s)
    pc = get_pulsed_vs_continuous(s)

    savefig(pp, "fig_3b.pdf")
    savefig(pc, "fig_3c.pdf")
    # p = plot(pp, pc, ps, layout=l)
    # f2 = "fig_2c.pdf"
    # savefig(p, f2)
end

plot_pulsed_comparison()

################################################################################
### FIGURE 3D
################################################################################

function three_state_model(dc, c, k, t)
    dc[1] = k['a'] * c[2] - k['s'] * c[1]
    dc[2] = k['s'] * c[1] - k['a'] * c[2] - k['i'] * c[2] + k['p'] * c[3]
    dc[3] = k['i'] * c[2] - k['p'] * c[3]

    @assert abs(dc[1]+dc[2]+dc[3])<=epsilon
end

function isoutofdomain(c, k, t)
    for i in 1:3:length(c)
        if abs(1-(c[i] + c[i+1] + c[i+2])) >= epsilon
            println(c)
            println(k)
            println(t)
            return false
        end
    end
    return true
end

function run_recruitment(k_time_pairs, c0=Nothing)
    t = fill(Float64, 0)
    u = fill(Float64[], 0)

    if (c0 == Nothing)
        c0 = [1.0; 0.0; 0.0]
    end

    for pair in k_time_pairs
        k = pair[1]
        times = pair[2]
        prob = ODEProblem(three_state_model, c0, times, k)
        sol = solve(prob, #AutoTsit5(Rosenbrock23());
                    # isoutofdomain=isoutofdomain,
                    saveat=0.01,
                    d_discontinuities=[times[1], times[2]])
        t = [t; sol.t]
        u = [u; sol.u]
        c0 = sol.u[end]
    end
    as = [x[1] for x in u]
    rs = [x[2] for x in u]
    is = [x[3] for x in u]
    data = [as rs is]

    return [t, data]
end

k_strong_sil_recruit = Dict(
    'a' => 0.31,
    's' => 11,
    'i' => 0.3,
    'p' => 0)

k_strong_sil_release = Dict(
    'a' => 0.31,
    's' => 0.1,
    'i' => 0,
    'p' => 0)

k_weaka_recruit = Dict(
    'a' => 0.69,
    's' => 0.1,
    'i' => 0,
    'p' => 0.1)

k_nomema_recruit = Dict('a' => 1.0,
                        's' => 0.1,
                        'i' => 0,
                        'p' => 0)

k_nomema_release = Dict('a' => 0,
                        's' => 0.1,
                        'i' => 0,
                        'p' => 0)

function plot_act(kp, kr, fname) # ka ks
    k_rec = Dict('a' => kp, 's' => kr, 'i' => 0, 'p' => 0)
    k_rel = Dict('a' => 0,  's' => kr, 'i' => 0, 'p' => 0)
    pair1 = [k_rec, (0, 5.0)]
    pair2 = [k_rel, (5.0, 15.0)]
    c0 = [0.0; 1.0; 0.0]
    sol = run_recruitment([pair1, pair2], c0)
    seq = plot()
    seq = vspan!(seq, [0, 5], label="", color=parse(Colorant, "#F5EBCE"))
    seq = vspan!(seq, [5, 15], label="", color=parse(Colorant, "#DEE6F0"))
    seq = plot!(seq,
                sol[1], sol[2][:, 2],
                color = :grey,
                label = "",
                xlabel="Days",
                xlims=[0, 15],
                xticks = 0:5:15,
                ylabel="Fraction Silent",
                yticks = 0:0.5:1.0,
                ylims = [0, 1.0],
                borderpad = 0.0,
                legendfontsize=13,
                linewidth=2)

    #  seq.o[:legend](borderaxespad=-1.0)
    seq = plot!(seq,
                size=(600, 300),
                legend=true)
    savefig(seq, fname)
    return seq
end

function plot_sil(ka, ks, ki, fname)
    k_rec = Dict('a' => ka, 's' => ks, 'i' => ki, 'p' => 0)
    k_rel = Dict('a' => ka, 's' => 0,  'i' => 0,  'p' => 0)
    pair1 = [k_rec, (0, 5.0)]
    pair2 = [k_rel, (5.0, 15.0)]
    c0 = [1.0; 0.0; 0.0]
    sol = run_recruitment([pair1, pair2], c0)
    seq = plot()
    seq = vspan!(seq, [0, 5], label="", color=parse(Colorant, "#F5EBCE"))
    seq = vspan!(seq, [5, 15], label="", color=parse(Colorant, "#DEE6F0"))
    seq = plot!(seq,
                sol[1], [1.0-a for a in sol[2][:, 1]],
                color = :grey,
                label = "",
                xlabel="Days",
                xlims=[0, 15],
                xticks = 0:5:15,
                ylabel="Fraction Silent",
                yticks = 0:0.5:1.0,
                ylims=[0, 1.0],
                borderpad = 0.0,
                legendfontsize=13,
                linewidth=2)

    #  seq.o[:legend](borderaxespad=-1.0)
    seq = plot!(seq,
                size=(600, 300),
                legend=true)
    savefig(seq, fname)
    return seq
end

plot_act(1.0, 1.0, "fig_3d_1.pdf")
plot_sil(1.1, 15.0, 0.5, "fig_3d_2.pdf")

################################################################################
### FIGURE 3E
################################################################################

function unpack_values(ari)
    return [ari['a'], ari['r'], ari['i']]
end

function unpack_parms(k)
    return [k['a'], k['s'], k['i'], k['p']]
end

function hill(k, a, exp=5, kd=0.5)
    if exp==Nothing
        exp = 5
    end
    if kd==Nothing
        kd = 0.5
    end
    return k * a #(a^exp/(a^exp + kd^exp))
end

function get_contrib(src_vals, dst_vals, kvals, he=Nothing, hk=Nothing, is_input=false, input_fn_val=Nothing)
    sa, sr, si     = src_vals
    da, dr, di     = dst_vals
    ka, ks, ki, kp = unpack_parms(kvals)

    if (sa == 0)
        return [0, 0, 0]
    end

    if (is_input)&&(input_fn_val != Nothing)
        sa = input_fn_val
    end

    dda =  dr*ka*sa - da*ks*sa
    ddr = -dr*ka*sa + da*ks*sa - dr*ki*sa + di*kp*sa
    ddi =                     dr*ki*sa - di*kp*sa

    # dda =  dr*ka - da*hill(ks, sa, he, hk)
    # ddr = -dr*ka + da*hill(ks, sa, he, hk) - dr*hill(ki, sa, he, hk) + di*hill(kp, sa, he, hk)
    # ddi =                                    dr*hill(ki, sa, he, hk) - di*hill(kp, sa, he, hk)

    @assert dda + ddr + ddi <= epsilon

    return [dda, ddr, ddi]
end

function get_bkgd_contrib(src_vals, kvals)
    sa, sr, si     = src_vals
    ka, ks, ki, kp = unpack_parms(kvals)

    dda =  sr*ka - sa*ks
    ddr = -sr*ka + sa*ks - sr*ki + si*kp
    ddi =                  sr*ki - si*kp

    @assert dda + ddr + ddi <= epsilon

    return [dda, ddr, ddi]
end

function get_k_bkgd_func(a=0.0, s=0.0, i=0.0, p=0.0)
    return function()
        k_bkgd = Dict('a' => a,
                      's' => s,
                      'i' => i,
                      'p' => p)
        return k_bkgd
    end
end

function get_kval_func(k_sil_delta, k_act_delta)
    f = function get_k_vals(src, dst, weight, t, input_id, t_spans, kvspec)
        k_bkgd = Dict('a' => 0.0,
                      's' => 0.0,
                      'i' => 0.0,
                      'p' => 0.0)

        function add_silencing(k)
            k['s'] += k_sil_delta['s']
            k['i'] += k_sil_delta['i']
            return k
        end

        function add_activation(k)
            k['a'] += k_act_delta['a']
            k['p'] += k_act_delta['p']
            return k
        end

        function get_k(weight)
            if (weight==1)
                return add_activation(k_bkgd)
            elseif (weight==-1)
                return add_silencing(k_bkgd)
            end
        end

        # if src is the input node, and t is not within the timezones, k values are 0
        # else return based on wghts
        if (src == input_id)
            input_on = false
            for span in t_spans
                t_start = span[1]
                t_end = span[2]
                if (t_start <= t) && (t <= t_end)
                    input_on = true
                end
            end
            if (! input_on)
                return  k_bkgd
            end
        end

        if kvspec != Nothing
            kvs = get(kvspec, src, Nothing)
            if kvs != Nothing
                return kvs
            end
        end
        return get_k(weight)
    end
    return f
end

function digraph_model(du, u, p, t)
    digraph    = p["digraph"]
    input_id   = p["input_id"]
    t_spans    = p["t_spans"]
    hill_exp   = p["hill_exp"]
    hill_kd    = p["hill_kd"]
    get_k_vals = p["get_k_vals"]
    get_k_bkgd = p["get_k_bkgd"]
    kvspec = get(p, "k_spec", Nothing)
    input_fn = get(p, "input_fn", Nothing)

    if (input_fn != Nothing)
        input_fn_val = input_fn(t)
    else
        input_fn_val = Nothing
    end

    if get_k_vals == Nothing
        get_k_vals = get_kval_func(10.0, 1.0)
    end
    if get_k_bkgd == Nothing
        get_k_bkgd = get_k_bkgd_func()
    end

    vs = collect(vertices(digraph))
    es = collect(edges(digraph))
    # vertex ==> a, r, i
    vdict = Dict(v => Dict("da" => 0.0,
                           "dr" => 0.0,
                           "di" => 0.0) for v in vs)

    # compute deltas
    for v in vs
        if (v != input_id)
            s_vals = u[3 * (v - 1) + 1 : (3 * (v - 1)) + 3] # compute background activation
            d_vals = s_vals
            k_vals = get_k_bkgd()
            deltas = get_bkgd_contrib(s_vals, k_vals)
            vdict[v]["da"] += deltas[1]
            vdict[v]["dr"] += deltas[2]
            vdict[v]["di"] += deltas[3]
            for e in [x for x in es if x.dst==v] # compute recruitment-driven activation/silencing
                src = e.src
                dst = e.dst

                @assert dst==v

                is_input = (src==input_id)

                weight = e.weight
                src_vals = u[(3 * (src - 1)) + 1 : (3 * (src - 1)) + 3]
                dst_vals = u[(3 * (dst - 1)) + 1 : (3 * (dst - 1)) + 3]
                k_vals = get_k_vals(src, dst, weight, t, input_id, t_spans, kvspec)
                deltas = get_contrib(src_vals, dst_vals, k_vals, hill_exp, hill_kd, is_input, input_fn_val)

                # if (deltas[1] > 0)&&(dst==length(vs))
                #     println(src)
                #     println(dst)
                #     println(weight)
                #     println(src_vals)
                #     println(dst_vals)
                #     println(k_vals)
                #     println(deltas)
                #     println()
                # end
                #
                vdict[v]["da"] += deltas[1]
                vdict[v]["dr"] += deltas[2]
                vdict[v]["di"] += deltas[3]
            end
        end
    end

    # return deltas inplace
    for idx in 1:length(vs)
        v = vs[idx]
        du[3 * (v-1) + 1] = vdict[v]["da"]
        du[3 * (v-1) + 2] = vdict[v]["dr"]
        du[3 * (v-1) + 3] = vdict[v]["di"]
    end
end

function run_digraph(parms, tlim, tstep=Nothing)
    digraph = parms["digraph"]
    num_vertices = length(collect(vertices(digraph)))

    if (parms["u0"] == Nothing)
        u0 = reduce(append!, [[1.0, 0.0, 0.0] for x in 1:num_vertices])
    else
        u0 = parms["u0"]
    end

    @assert length(u0) == 3*num_vertices

    tspan = (0.0, tlim)

    prob = ODEProblem(digraph_model, u0, tspan, parms)

    if (tstep==Nothing)
        sol = solve(prob, d_discontinuities=collect(Iterators.flatten(parms["t_spans"])))#, saveat=0.001)
    else
        sol = solve(prob, d_discontinuities=collect(Iterators.flatten(parms["t_spans"])), saveat=tstep)
    end

    t = sol.t
    u = sol.u

    avals = [[x[(3*(i-1))+1] for x in u] for i in 1:num_vertices]
    data = hcat(avals...)

    return [t, data]
end

function get_digraph_data(iff_parms, t_span)
    sol = run_digraph(iff_parms, t_span)
    labels = [string("CR ", i) for i in 1:length(collect(vertices(iff_parms["digraph"])))]
    labels = hcat(labels...)
    pdata = sol[2]
    return [sol[1], pdata, labels]
end

function plot_digraph(iff_parms, t_span, plot_input=true, plot_outputs_only=false, plot_last_only=false)
    sol = run_digraph(iff_parms, t_span)
    labels = [string("CR ", i) for i in 1:length(collect(vertices(iff_parms["digraph"])))]
    labels = hcat(labels...)
    if plot_input
        pdata = sol[2]
    elseif plot_outputs_only
        pdata = sol[2][:, end-1:end]
    elseif plot_last_only
        pdata = sol[2][:, end]
    else
        pdata = sol[2][:, 2:end]
    end
    p1 = plot(sol[1], pdata,
              label=labels,
              palette=:seaborn_colorblind,
              xlabel="Days",
              ylabel="Fraction Active",
              ylims=(-0.1, 1.1),
              yticks=0:0.5:1.0,
              xticks=0:t_span/2.0:t_span,
              linewidth=2,
              legend=:outerright,
              size=(700, 300))

    for span in iff_parms["t_spans"]
        p1 = vspan!(p1, span, label="", color=:gray, alpha=0.2)
    end
    return p1
end

function construct_parms(digraph, iid, spans,
                          u0=Nothing, hill_exp=Nothing, hill_kd=Nothing,
                          get_k_bkgd=Nothing, get_k_vals=Nothing)
    parms = Dict("digraph"  => digraph,
                  "input_id" => iid,
                  "t_spans"  => spans,
                  "u0" => u0,
                  "hill_exp" => hill_exp,
                  "hill_kd" => hill_kd,
                  "get_k_vals" => get_k_vals,
                  "get_k_bkgd" => get_k_bkgd)
    return parms
end

colors = get(ColorSchemes.viridis, reverse([0.3, 0.6, 0.9]))

sources = [1, 1, 2]
destinations = [2, 3, 3]
wghts = [1, 1, -1]
digraph = SimpleWeightedDiGraph(sources, destinations, wghts)
iid = 1
spans = [[20, 40], [60, 80]]
u0 = [1.0, 0.0, 0.0,
      0.0, 1.0, 0.0,
      0.0, 1.0, 0.0]
get_k_bkgd = get_k_bkgd_func(0, 1.0)

function ifn(t)
    if t <= 20
        return 0
    elseif t <= 40
        return 1
    elseif t<= 60
        return 0
    elseif t<= 80
        return 1
    else
        return 0
    end
end

iff_parms = construct_parms(digraph, iid, spans, u0, Nothing, Nothing, get_k_bkgd, Nothing)
push!(iff_parms, "k_spec" => [Dict('a' => 1.0, 's' => 0.0,  'i' => 0, 'p' => 0),
                               Dict('a' => 0.1, 's' => 15.0, 'i' => 0.5, 'p' => 0),
                               Dict('a' => 2.0, 's' => 0.0,  'i' => 0, 'p' => 0)])
# push!(iff_parms, "input_fn" => ifn)

colors  = ["turquoise4" "chocolate2" "rebeccapurple"]
colors = ColorSchemes.Dark2_3
t, data, ls = get_digraph_data(iff_parms, 100.0)
t_span = 100

p1 = plot(ifn, 0, 100,
          label="",
          xlabel="",
          ylabel="Input",
          guidefontsize = 18,
          ylims=(-0.1, 1.1),
          yticks=0:0.5:1.0,
          xticks=0:t_span/2.0:t_span,
          color=colors[1],
          linewidth=4)

p2 = plot(t, data[:, end-1], label = "",
          xlabel="",
          ylabel="Neg. Arm",
          guidefontsize = 18,
          ylims=(-0.1, 1.1),
          yticks=0:0.5:1.0,
          xticks=0:t_span/2.0:t_span,
          linewidth=4,
          legend=:topright,
          color=colors[2],
          size=(600, 400))

p3 = plot(t, data[:, end], label = "",
          xlabel="Days",
          ylabel="Output",
          guidefontsize = 18,
          ylims=(-0.02, 0.2),
          yticks=0:0.1:0.2,
          xticks=0:t_span/2.0:t_span,
          linewidth=4,
          legend=:topright,
          color=colors[3],
          size=(600, 500))

lay = @layout [a; b; c]
p = plot(p1, p2, p3, layout=lay, link=:all)
savefig(p, "fig_3e.pdf")

################################################################################
### FIGURE 3 + S2 Databases
################################################################################

function get_adapt_response(ka, ks, ki)
    sources = [1, 1, 2]
    destinations = [2, 3, 3]
    wghts = [1, 1, -1]
    digraph = SimpleWeightedDiGraph(sources, destinations, wghts)
    iid = 1
    spans = [[20, 40], [60, 80]]
    u0 = [1.0, 0.0, 0.0,
          0.0, 1.0, 0.0,
          0.0, 1.0, 0.0]
    get_k_bkgd = get_k_bkgd_func(0, 1.0)
    iff_parms = construct_parms(digraph, iid, spans, u0, Nothing, Nothing, get_k_bkgd, Nothing)
    push!(iff_parms, "k_spec" => [Dict('a' => ka, 's' => 0.0,  'i' => 0, 'p' => 0),
                                  Dict('a' => 0.1, 's' => ks, 'i' => ki, 'p' => 0),
                                  Dict('a' => 0.0, 's' => 0.0,  'i' => 0, 'p' => 0)])

    t, data, ls = get_digraph_data(iff_parms, 100.0)
    output = data[:, end]

    peak = maximum(output)
    pkind = findfirst(o->o==peak, output)

    sind = findfirst(time->time>=30, t)
    set = output[sind]

    @assert t[sind] >= t[pkind]
    @assert set >= 0
    @assert peak > 0

    k = -log(set/peak)/(t[sind] - t[pkind])

    pk2ind = findfirst(time->time>=50, t)
    pk2 = maximum(output[pk2ind:end])

    return [peak, log(2)/k, pk2]
end

function get_response_dbase_3f(small = false)
    karange = 0.01:0.1:20.01
    ksrange = 0.01:0.1:20.01
    kirange = [0.5]

    if small
        karange = 1:2
        ksrange = 1:10
        kirange = 0.25:0.25:0.5
    end

    dblen = length(karange) * length(ksrange) * length(kirange)
    println("Should take approx. " * string(dblen * 0.08) * " seconds")

    kas = zeros(dblen)
    kss = zeros(dblen)
    kis = zeros(dblen)
    pks = zeros(dblen)
    pk2s = zeros(dblen)
    ads = zeros(dblen)
    global idx = 1

    for ks in ksrange
        println("Working on ks ", ks)
        for ka in karange
            for ki in kirange
                pk, ad, pk2 = get_adapt_response(ka, ks, ki)
                kas[idx] = ka
                kss[idx] = ks
                kis[idx] = ki
                pks[idx] = pk
                pk2s[idx] = pk2
                ads[idx] = ad
                global idx = idx + 1
            end
        end
    end

    df = DataFrame(Dict("ka" => kas,
                        "ks" => kss,
                        "ki" => kis,
                        "ad" => ads,
                        "pk" => pks,
                        "pk2" => pk2s))
    if (! small)
        CSV.write("../loop_simulations/3f_response_peaks.csv", df)
    end
end

function get_response_dbase(small = false)
    #  karange = sort([collect(0.01:0.5:20.01); collect(0.01:0.05:0.5); 0.5])
    #  ksrange = sort([collect(0.01:0.5:20.01); collect(0.01:0.05:0.5); 0.5])
    #  kirange = sort([collect(0.01:0.2:2.01); collect(0.01:0.02:0.21); 0.5])
    #
    #  # spanning in log space
    #  karange = sort([exp10.(range(log10(0.01), log10(20.01), length=20)); 0.5])
    #  ksrange = sort([exp10.(range(log10(0.01), log10(20.01), length=20)); 0.5])
    #
    karange = 0.01:1.0:20.01
    ksrange = 0.01:1.0:20.01

    karange = sort([collect(karange); 0.05; 0.10; 0.25; 0.5; 1.0])
    ksrange = sort([collect(ksrange); 0.05; 0.10; 0.25; 0.5; 1.0])
    kirange = sort([exp10.(range(log10(0.01), log10(2.01), length=length(karange))); 0.5])

    if small
        karange = 1:2
        ksrange = 1:10
        kirange = 0.25:0.25:0.5
    end

    dblen = length(karange) * length(ksrange) * length(kirange)
    println("Should take approx. " * string(dblen * 0.08) * " seconds")

    kas = zeros(dblen)
    kss = zeros(dblen)
    kis = zeros(dblen)
    pks = zeros(dblen)
    pk2s = zeros(dblen)
    ads = zeros(dblen)
    global idx = 1

    for ks in ksrange
        println("Working on ks ", ks)
        for ka in karange
            for ki in kirange
                pk, ad, pk2 = get_adapt_response(ka, ks, ki)
                kas[idx] = ka
                kss[idx] = ks
                kis[idx] = ki
                pks[idx] = pk
                pk2s[idx] = pk2
                ads[idx] = ad
                global idx = idx + 1
            end
        end
    end

    df = DataFrame(Dict("ka" => kas,
                        "ks" => kss,
                        "ki" => kis,
                        "ad" => ads,
                        "pk" => pks,
                        "pk2" => pk2s))
    if (! small)
        CSV.write("../loop_simulations/20_response_peaks.csv", df)
    end
end

t_span = 3

function get_cascade(n_mediators, ks, ka, input_level)
    input = 1
    mediators = collect(2:1+n_mediators)
    control = n_mediators+2
    output = n_mediators+3

    sources = zeros(0)
    destinations = zeros(0)
    wghts = zeros(0)

    for m in mediators
        push!(sources, input)
        push!(destinations, m)
        push!(wghts, 1.0)

        push!(sources, m)
        push!(destinations, output)
        push!(wghts, -1.0)
    end

    if (length(mediators)==0)
        push!(sources, input)
        push!(destinations, output)
        push!(wghts, -1.0)
    end

    push!(sources, input)
    push!(destinations, control)
    push!(wghts, -1.0)

    sources = [Int(s) for s in sources]
    destinations = [Int(d) for d in destinations]

    digraph = SimpleWeightedDiGraph(sources, destinations, wghts)

    iid = 1
    spans = [[0, t_span]]

    u0 = zeros(0)
    append!(u0, [input_level, 1.0-input_level, 0.0])
    for m in mediators
        append!(u0, [0.0, 1.0, 0.0])
    end
    append!(u0, [1.0, 0.0, 0.0])
    append!(u0, [1.0, 0.0, 0.0])

    get_k_bkgd = get_k_bkgd_func(0.0)
    get_k_vals = get_kval_func(Dict('s' => ks, 'i' => 0), Dict('a' => ka, 'p' => 0))
    iff_parms = construct_parms(digraph, iid, spans, u0, Nothing, Nothing, get_k_bkgd, get_k_vals)
    sol = run_digraph(iff_parms, t_span, 0.001)
    return sol
end

function plot_solutions(solutions, labels, fname="none")
    @assert length(solutions) == length(labels)

    cscheme = get(ColorSchemes.viridis, reverse(0.1:0.8/(length(labels)):0.9))

    @. hfn(x,p) = p[1] + p[2] * ((x+0im)^p[3])/((x+0im)^p[3] + (p[4]+0im)^p[3])
    hcoeffs = zeros(length(labels))
    kd50s = zeros(length(labels))

    p = plot()
    for i in 1:length(labels)
        sol = solutions[i]
        lab = labels[i]
        xvals = sol[1]
        yvals = [y for y in sol[2][:, end]]

        fit = curve_fit(hfn, xvals, yvals, [1.0, -1.0, 1.0, 1.0])
        if fit.converged
            hcoeffs[i] = fit.param[3]
            kd50s[i] = fit.param[4]
        end

        if lab == "0 Repressor"
            p = plot!(p,
                      xvals, yvals,
                      label="Direct Repression", color=:grey,
                      linewidth=2, linestyle=:dash)
        else
            if lab != "1 Repressor"
                lab = lab * "s"
            end
            p = plot!(p,
                      xvals, yvals,
                      label=lab, color=cscheme[i],
                      linewidth=2, linestyle=:solid)
        end
    end

    p = plot!(p,
              #  color_palette=ColorSchemes.viridis,
              xlabel="Days of Recruitment",
              ylabel="Fraction Active",
              ylims=(-0.1, 1.1),
              yticks=0:0.5:1.0,
              xlims=(0, t_span),
              xticks=0:t_span/2.0:t_span,
              legend=:topright,
              #  legendtitle="# Repressors",
              size=(450, 300))

    if fname != "none"
        savefig(p, fname)
    end

    return [hcoeffs, kd50s]
end

function get_coeffs(n_mediators, ks, ka, input_level)
    sol = get_cascade(n_mediators, ks, ka, input_level)

    @. hfn(x,p) = 1 - ((x)^p[1])/((x)^p[1] + (p[2])^p[1])

    xvals = sol[1]
    yvals = [y for y in sol[2][:, end]]
    fit = curve_fit(
                    hfn, xvals, yvals, [1.0, 1.0],
                    lower=[0.0, 0.0], autodiff=:forwarddiff)

    if fit.converged
        return fit.param
    else
        return [NaN, NaN]
    end
end

n_reps = [0, 2, 5, 20]
parms = [[n, 2, 1.0, 1.0] for n in n_reps]
cascades = [get_cascade(p...) for p in parms]
plot_solutions(cascades, [string(n)*" Repressor" for n in n_reps], "fig_4j.pdf");

function get_hill_dbase(small = false)
    if ! small
        nmrange = [1, 2, 5, 10, 20]
        ksrange = 0:0.25:20
        karange = 0:0.1:2
        karange = 0:0.25:20
    else
        nmrange = 19:20
        ksrange = 19:20
        karange = 19:20
    end

    dflen = length(nmrange) * length(ksrange) * length(karange)
    println("Should take about " * string(dflen * 0.2) * " seconds")
    nms = zeros(dflen)
    kss = zeros(dflen)
    kas = zeros(dflen)
    ins = zeros(dflen)
    hcs = zeros(dflen)
    kds = zeros(dflen)
    idx = 1

    t1 = time()
    for n_m in reverse(nmrange)
        for ks in reverse(ksrange)
            for ka in reverse(karange)
                @time h, k = get_coeffs(n_m, ks, ka, 1.0)
                nms[idx] = n_m
                kss[idx] = ks
                kas[idx] = ka
                ins[idx] = 1.0
                hcs[idx] = h
                kds[idx] = k
                idx = idx + 1
            end
        end
        GC.gc()
    end
    t2 = time()

    #  println(t2 - t1)

    @assert idx - 1 == dflen

    cascade_df = DataFrame(Dict("nmediators" => nms,
                                "ks" => kss,
                                "ka" => kas,
                                "input" => ins,
                                "hillcoeff" => hcs,
                                "kd50" => kds))
    CSV.write("../loop_simulations/hill_coeffs.csv", cascade_df)
end

#  hill_dbase(true)
#  @time hill_dbase(false)

println("To generate full .csv files run:")
println("    get_response_dbase(false)")
println("    get_hill_dbase(false)")


text:
@testset "Interpolate" begin
    @testset "Linear1D" begin
        β = π
        tgrid = SimpleGrid.Uniform{Float64}([0.0, β], 33)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(t) = t
        data = zeros(tgrid.size)

        for (ti, t) in enumerate(tgrid.grid)
            data[ti] = f(t)
        end

        for ti = 1:tgrid.size - 1
            t = tgrid[ti] + 1.e-6
            fbar = Interp.interp1D(data, tgrid, t)
            @test abs(f(tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test f(tgrid[ti]) < fbar
            @test f(tgrid[ti + 1]) > fbar
        end
        for ti = 2:tgrid.size
            t = tgrid[ti] - 1.e-6
            fbar = Interp.interp1D(data, tgrid, t)
            @test abs(f(tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test f(tgrid[ti]) > fbar
            @test f(tgrid[ti - 1]) < fbar
        end

        t = tgrid[1] + eps(Float64)*1e3
        fbar = Interp.interp1D(data, tgrid, t)
        @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt

        t = tgrid[tgrid.size] - eps(Float64)*1e3
        fbar = Interp.interp1D(data, tgrid, t)
        @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt

        tlist = rand(10) * β
        # println(tlist)

        for (ti, t) in enumerate(tlist)
            fbar = Interp.interp1D(data, tgrid, t)
            # println("$k, $t, $fbar, ", f(k, t))
            @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
        end

        data2 = zeros((2,tgrid.size))

        for i in 1:2
            for (ti, t) in enumerate(tgrid.grid)
                data2[i,ti] = f(t)
            end
        end
        tlist = rand(10) * β

        fbars=Interp.interp1DGrid(data2, tgrid, tlist, axis=2)
        for (ti, t) in enumerate(tlist)
            fbar = Interp.interp1D(data2, tgrid, t, axis=2)
            @test abs(f(t) - fbar[1]) < 3.e-6 # linear interpolation, so error is δK+δt
            @test abs(f(t) - fbar[2]) < 3.e-6 # linear interpolation, so error is δK+δt
            @test abs(f(t) - fbars[1, ti]) < 3.e-6 # linear interpolation, so error is δK+δt
            @test abs(f(t) - fbars[2, ti]) < 3.e-6 # linear interpolation, so error is δK+δt
        end

    end

    @testset "LinearND" begin
        β = π
        tgrid = SimpleGrid.Uniform{Float64}([0.0, β], 22)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(t1,t2) = 1.0+t1+t2+t1*t2
        data = zeros((tgrid.size,tgrid.size))

        for (ti, t) in enumerate(tgrid.grid)
            for (di, d) in enumerate(tgrid.grid)
                data[ti, di] = f(t,d)
            end
        end

        for ti = 1:tgrid.size - 1
            for di = 1:tgrid.size - 1
                t = tgrid[ti] + 1.e-6
                d = tgrid[di] + 1.e-6
                fbar = Interp.linearND(data, [tgrid,tgrid], [t,d])
                @test abs(f(tgrid[ti],tgrid[di]) - fbar) < 9.e-6 # linear interpolation, so error is dδt+tδd+δt+δd
                @test f(tgrid[ti],tgrid[di]) < fbar
                @test f(tgrid[ti + 1],tgrid[di+1]) > fbar
            end
        end
        for ti = 2:tgrid.size
            for di = 2:tgrid.size
                t = tgrid[ti] - 1.e-6
                d = tgrid[di] - 1.e-6
                fbar = Interp.linearND(data, [tgrid,tgrid], [t,d])
                @test abs(f(tgrid[ti],tgrid[di]) - fbar) < 9.e-6 # linear interpolation, so error is dδt+tδd+δt+δd
                @test f(tgrid[ti],tgrid[di]) > fbar
                @test f(tgrid[ti - 1],tgrid[di - 1]) < fbar
            end
        end

        tlist = rand(10) * β
        dlist = rand(10) * β
        # println(tlist)

        for (ti, t) in enumerate(tlist)
            for (di, d) in enumerate(tlist)
                fbar = Interp.linearND(data, [tgrid,tgrid], [t,d])
                @test abs(f(t,d) - fbar) < 9.e-6 # linear interpolation, so error is δK+δt
            end
        end
    end

    @testset "interpND" begin
        β = π
        tgrid = SimpleGrid.Uniform{Float64}([0.0, β], 22)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(t1,t2) = 1.0+t1+t2+t1*t2
        data = zeros((tgrid.size,tgrid.size))

        for (ti, t) in enumerate(tgrid.grid)
            for (di, d) in enumerate(tgrid.grid)
                data[ti, di] = f(t,d)
            end
        end

        for ti = 1:tgrid.size - 1
            for di = 1:tgrid.size - 1
                t = tgrid[ti] + 1.e-6
                d = tgrid[di] + 1.e-6
                fbar = Interp.interpND(data, [tgrid,tgrid], [t,d])
                @test abs(f(tgrid[ti],tgrid[di]) - fbar) < 9.e-6 # linear interpolation, so error is dδt+tδd+δt+δd
                @test f(tgrid[ti],tgrid[di]) < fbar
                @test f(tgrid[ti + 1],tgrid[di+1]) > fbar
            end
        end
        for ti = 2:tgrid.size
            for di = 2:tgrid.size
                t = tgrid[ti] - 1.e-6
                d = tgrid[di] - 1.e-6
                fbar = Interp.interpND(data, [tgrid,tgrid], [t,d])
                @test abs(f(tgrid[ti],tgrid[di]) - fbar) < 9.e-6 # linear interpolation, so error is dδt+tδd+δt+δd
                @test f(tgrid[ti],tgrid[di]) > fbar
                @test f(tgrid[ti - 1],tgrid[di - 1]) < fbar
            end
        end

        tlist = rand(10) * β
        dlist = rand(10) * β
        # println(tlist)

        for (ti, t) in enumerate(tlist)
            for (di, d) in enumerate(tlist)
                fbar = Interp.interpND(data, [tgrid,tgrid], [t,d])
                @test abs(f(t,d) - fbar) < 9.e-6 # linear interpolation, so error is δK+δt
            end
        end
    end

    @testset "BaryCheb" begin
        β = π
        tgrid = SimpleGrid.BaryCheb{Float64}([0.0, β], 16)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(t) = t
        data = zeros(tgrid.size)

        for (ti, t) in enumerate(tgrid.grid)
            data[ti] = f(t)
        end

        for ti = 1:tgrid.size - 1
            t = tgrid[ti] + 1.e-6
            fbar = Interp.interp1D(data, tgrid, t)
            @test abs(f(tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test f(tgrid[ti]) < fbar
            @test f(tgrid[ti + 1]) > fbar
        end
        for ti = 2:tgrid.size
            t = tgrid[ti] - 1.e-6
            fbar = Interp.interp1D(data, tgrid, t)
            @test abs(f(tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test f(tgrid[ti]) > fbar
            @test f(tgrid[ti - 1]) < fbar
        end

        t = tgrid[1] + eps(Float64)*1e3
        fbar = Interp.interp1D(data, tgrid, t)
        @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt

        t = tgrid[tgrid.size] - eps(Float64)*1e3
        fbar = Interp.interp1D(data, tgrid, t)
        @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt

        tlist = rand(10) * β
        # println(tlist)

        for (ti, t) in enumerate(tlist)
            fbar = Interp.interp1D(data, tgrid, t)
            # println("$k, $t, $fbar, ", f(k, t))
            @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
        end
    end

    @testset "DensedLog" begin
        β = 4
        tgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, β], [0.0, 0.5β, β], 4, 0.001, 4)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(t) = t
        data = zeros(tgrid.size)

        for (ti, t) in enumerate(tgrid.grid)
            data[ti] = f(t)
        end

        for ti = 1:tgrid.size - 1
            t = tgrid[ti] + 1.e-6
            fbar = Interp.interp1D(data, tgrid, t)
            @test abs(f(tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test f(tgrid[ti]) < fbar
            @test f(tgrid[ti + 1]) > fbar
        end
        for ti = 2:tgrid.size
            t = tgrid[ti] - 1.e-6
            fbar = Interp.interp1D(data, tgrid, t)
            @test abs(f(tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test f(tgrid[ti]) > fbar
            @test f(tgrid[ti - 1]) < fbar
        end

        t = tgrid[1] + eps(Float64)*1e3
        fbar = Interp.interp1D(data, tgrid, t)
        @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt

        t = tgrid[tgrid.size] - eps(Float64)*1e3
        fbar = Interp.interp1D(data, tgrid, t)
        @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt

        tlist = rand(10) * β
        tlist = sort(tlist)
        println(tlist)
        ff = Interp.interp1DGrid(data, tgrid, tlist)
        ff_c = similar(ff)

        for (ti, t) in enumerate(tlist)
            fbar = Interp.interp1D(data, tgrid, t)
            ff_c[ti] = fbar
            # println("$k, $t, $fbar, ", f(k, t))
            @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test abs(f(t) - ff[ti]) < 3.e-6 # linear interpolation, so error is δK+δt
        end
        # println(tgrid.grid)
        # println("ff_c:",ff_c)
        # println("ff:",ff)
    end

    @testset "Integrate" begin
        β = 1.0
        tgrid = SimpleGrid.GaussLegendre{Float64}([0.0, β], 4)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(t) = t
        data = zeros(tgrid.size)
        for (ti, t) in enumerate(tgrid.grid)
            data[ti] = f(t)
        end
        println(tgrid.grid)
        println(data)
        println(tgrid.weight)
        println(sum(data.*tgrid.weight))
        int_result = Interp.integrate1D(data, tgrid)
        @test abs(int_result - 0.5) < 3.e-6

        β = 1.0
        tgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, β], [0.0, 0.5β, β], 2, 0.001, 3)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        data = zeros(2, tgrid.size)
        for (ti, t) in enumerate(tgrid.grid)
            data[1,ti] = f(t)
            data[2,ti] = f(t)
        end
        println(tgrid.grid)
        println(data)
        int_result = Interp.integrate1D(data, tgrid; axis=2)
        @test abs(int_result[1] - 0.5) < 3.e-6
        @test abs(int_result[2] - 0.5) < 3.e-6
    end

    @testset "Find Neighbor and Interp Sliced" begin
        β = 4
        tgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, β], [0.0, 0.5β, β], 4, 0.001, 4)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(t) = t
        data = zeros(tgrid.size)

        for (ti, t) in enumerate(tgrid.grid)
            data[ti] = f(t)
        end

        for ti = 1:tgrid.size - 1
            t = tgrid[ti] + 1.e-6
            neighbor = Interp.findneighbor(Interp.LinearInterp(),tgrid, t)
            @test neighbor.index[1] == floor(tgrid, t)
            @test neighbor.index[2] == floor(tgrid, t) + 1

            neighbor = Interp.findneighbor(tgrid, t)
            data_slice = Interp.dataslice(data,neighbor.index)#data[neighbor.index[1]:neighbor.index[2]]
            fbar = Interp.interpsliced(neighbor, data_slice)
            @test abs(f(tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test f(tgrid[ti]) < fbar
            @test f(tgrid[ti + 1]) > fbar
        end
        for ti = 2:tgrid.size
            t = tgrid[ti] - 1.e-6
            fbar = Interp.interp1D(data, tgrid, t)
            neighbor = Interp.findneighbor(Interp.LinearInterp(),tgrid, t)
            @test neighbor.index[1] == floor(tgrid, t)
            @test neighbor.index[2] == floor(tgrid, t) + 1

            neighbor = Interp.findneighbor(tgrid, t)
            data_slice = Interp.dataslice(data,neighbor.index)#data[neighbor.index[1]:neighbor.index[2]]
            fbar = Interp.interpsliced(neighbor, data_slice)
            @test abs(f(tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test f(tgrid[ti]) > fbar
            @test f(tgrid[ti - 1]) < fbar
        end

        t = tgrid[1] + eps(Float64)*1e3
        fbar = Interp.interp1D(data, tgrid, t)
        @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt

        t = tgrid[tgrid.size] - eps(Float64)*1e3
        fbar = Interp.interp1D(data, tgrid, t)
        @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt

    end

end

ext: jl


text:
using IDSGraphs
using Graphs

using Test

function test_graph()
    g = IDSGraph()
    add_edge!(g, 'a', 'b')
    add_edge!(g, 'a', 'e')
    add_edge!(g, 'b', 'c')
    add_edge!(g, 'b', 'd')
    add_edge!(g, 'e', 'f')
    add_edge!(g, 'f', 'd')
    add_edge!(g, 'c', 'i')
    add_edge!(g, 'c', 'j')
    add_edge!(g, 'd', 'g')
    add_edge!(g, 'd', 'h')

    return g
end

@testset "unit tests" begin
    @testset "linear graph creation" begin
        g = IDSGraph()
        @test nv(g) == 0
        @test ne(g) == 0

        for (n, c) in enumerate('a':'z')
            add_vertex!(g, c)
            @test nv(g) == n
            @test ne(g) == 0
        end

        add_edge!(g, 'a', 'z')
        @test ne(g) == 1
        @test has_edge(g, 'a', 'z')
    end

    @testset "verify test graph" begin
        g = test_graph()
        @test nv(g) == 10
        @test ne(g) == 10

        for c in 'a':'j'
            @test has_vertex(g, c)
        end

        @test length(symdiff(unique(vertices(g)), 'a':'j')) == 0

        # test inneighbors
        @test length(inneighbors(g, 'a')) == 0
        @test length(inneighbors(g, 'b')) == 1
        @test length(inneighbors(g, 'c')) == 1
        @test length(inneighbors(g, 'e')) == 1
        @test length(inneighbors(g, 'f')) == 1
        @test length(inneighbors(g, 'g')) == 1
        @test length(inneighbors(g, 'h')) == 1
        @test length(inneighbors(g, 'i')) == 1
        @test length(inneighbors(g, 'j')) == 1
        @test length(inneighbors(g, 'd')) == 2

        # test outneighbors
        @test length(outneighbors(g, 'a')) == 2
        @test length(outneighbors(g, 'b')) == 2
        @test length(outneighbors(g, 'c')) == 2
        @test length(outneighbors(g, 'e')) == 1
        @test length(outneighbors(g, 'f')) == 1
        @test length(outneighbors(g, 'g')) == 0
        @test length(outneighbors(g, 'h')) == 0
        @test length(outneighbors(g, 'i')) == 0
        @test length(outneighbors(g, 'j')) == 0
        @test length(outneighbors(g, 'd')) == 2
    end

    @testset "graph properties/invariants" begin
        g = test_graph()

        # zero (identity) graph
        zero_graph = zero(g)
        @test nv(zero_graph) == 0
        @test ne(zero_graph) == 0
        @test length(collect(vertices(zero_graph))) == 0
        @test length(collect(edges(zero_graph))) == 0

        # types
        @test eltype(typeof(g)) == Char
        @test edgetype(g) == IDSGraphs.CharEdge

        # directed graph type
        @test is_directed(typeof(g))
    end

    @testset "conditional subgraphs" begin
        subg = IDSGraphs.subgraph(v -> v in 'b':'g', test_graph())

        @test nv(subg) == 6
        @test ne(subg) == 5

        @test length(symdiff(unique(vertices(subg)), 'b':'g')) == 0

        @test has_edge(subg, 'b', 'c')
        @test has_edge(subg, 'b', 'd')
        @test has_edge(subg, 'e', 'f')
        @test has_edge(subg, 'f', 'd')
        @test has_edge(subg, 'd', 'g')
    end
end


text:
"Holds the tableau of an Partitioned Additive Runge-Kutta method for Hamiltonian systems."
TableauHPARK = TableauVPARK

"Parameters for right-hand side function of Hamiltonian partitioned additive Runge-Kutta methods."
mutable struct ParametersHPARK{DT,TT,D,S,R,VT,FT,UT,GT,ϕT} <: Parameters{DT,TT}
    f_v::VT
    f_f::FT
    f_u::UT
    f_g::GT
    f_ϕ::ϕT

    Δt::TT

    t_q::CoefficientsARK{TT}
    t_p::CoefficientsARK{TT}
    t_q̃::CoefficientsPRK{TT}
    t_p̃::CoefficientsPRK{TT}
    t_λ::CoefficientsMRK{TT}

    t::TT
    q::Vector{DT}
    p::Vector{DT}
    λ::Vector{DT}

    function ParametersHPARK{DT,TT,D,S,R,VT,FT,UT,GT,ϕT}(f_v, f_f, f_u, f_g, f_ϕ, Δt, t_q, t_p, t_q̃, t_p̃, t_λ) where {DT,TT,D,S,R,VT,FT,UT,GT,ϕT}
        # create solution vectors
        q = zeros(DT,D)
        p = zeros(DT,D)
        λ = zeros(DT,D)

        new(f_v, f_f, f_u, f_g, f_ϕ, Δt,
            t_q, t_p, t_q̃, t_p̃, t_λ,
            zero(TT), q, p, λ)
    end
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheVPARK{ST,TT,D,S,R},
                                        params::ParametersHPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local tqᵢ::TT
    local tpᵢ::TT
    local tλᵢ::TT

    for i in 1:S
        for k in 1:D
            # copy x to Y, Z
            cache.Yi[i][k] = x[2*(D*(i-1)+k-1)+1]
            cache.Zi[i][k] = x[2*(D*(i-1)+k-1)+2]

            # compute Q and P
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
            cache.Pi[i][k] = params.p[k] + params.Δt * cache.Zi[i][k]
        end

        # compute f(X)
        tqᵢ = params.t + params.Δt * params.t_q.c[i]
        tpᵢ = params.t + params.Δt * params.t_p.c[i]
        params.f_v(tqᵢ, cache.Qi[i], cache.Pi[i], cache.Vi[i])
        params.f_f(tpᵢ, cache.Qi[i], cache.Pi[i], cache.Fi[i])
    end

    for i in 1:R
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+2]
            cache.Λp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+3]

            # compute Q and V
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        tλᵢ = params.t + params.Δt * params.t_λ.c[i]
        params.f_u(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Up[i])
        params.f_g(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Gp[i])
        params.f_ϕ(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Φp[i])
    end
end

"Compute stages of variational partitioned additive Runge-Kutta methods."
@generated function function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersHPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    cache = IntegratorCacheVPARK{ST,TT,D,S,R}()

    quote
        compute_stages!(y, $cache, params)

        # compute b = - [(Y-AV-AU), (Z-AF-AG)]
        for i in 1:S
            for k in 1:D
                b[2*(D*(i-1)+k-1)+1] = - $cache.Yi[i][k]
                b[2*(D*(i-1)+k-1)+2] = - $cache.Zi[i][k]
                for j in 1:S
                    b[2*(D*(i-1)+k-1)+1] += params.t_q.a[i,j] * $cache.Vi[j][k]
                    b[2*(D*(i-1)+k-1)+2] += params.t_p.a[i,j] * $cache.Fi[j][k]
                end
                for j in 1:R
                    b[2*(D*(i-1)+k-1)+1] += params.t_q.α[i,j] * $cache.Up[j][k]
                    b[2*(D*(i-1)+k-1)+2] += params.t_p.α[i,j] * $cache.Gp[j][k]
                end
            end
        end

        # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
        for i in 1:R
            for k in 1:D
                b[2*D*S+3*(D*(i-1)+k-1)+1] = - $cache.Yp[i][k]
                b[2*D*S+3*(D*(i-1)+k-1)+2] = - $cache.Zp[i][k]
                b[2*D*S+3*(D*(i-1)+k-1)+3] = - $cache.Φp[i][k]
                for j in 1:S
                    b[2*D*S+3*(D*(i-1)+k-1)+1] += params.t_q̃.a[i,j] * $cache.Vi[j][k]
                    b[2*D*S+3*(D*(i-1)+k-1)+2] += params.t_p̃.a[i,j] * $cache.Fi[j][k]
                end
                for j in 1:R
                    b[2*D*S+3*(D*(i-1)+k-1)+1] += params.t_q̃.α[i,j] * $cache.Up[j][k]
                    b[2*D*S+3*(D*(i-1)+k-1)+2] += params.t_p̃.α[i,j] * $cache.Gp[j][k]
                end
            end
        end

        # compute b = - [Λ₁-λ]
        if params.t_λ.c[1] == 0
            for k in 1:D
                b[2*D*S+3*(k-1)+3] = - $cache.Λp[1][k] + params.λ[k]
            end
        end
    end
end

@doc raw"""
Partitioned Additive Runge-Kutta integrator for Hamiltonian systems subject
to Dirac constraints.

This integrator solves the following system of equations for the internal stages,
```math
\begin{align}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= 0 , & i &= 1, ..., r ,
\end{align}
```
with definitions
```math
\begin{align}
V_{n,i} &= \hphantom{-} \frac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\
F_{n,i} &=           -  \frac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\
U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
G_{n,i} &=           -  \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r ,
\end{align}
```
and update rule
```math
\begin{align}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} .
\end{align}
```
"""
struct IntegratorHPARK{DT, TT, ET <: PDAE{DT,TT},
                               PT <: ParametersHPARK{DT,TT},
                               ST <: NonlinearSolver{DT},
                               IT <: InitialGuessPODE{DT,TT}} <: AbstractIntegratorSPARK{DT, TT}
    equation::ET
    tableau::TableauHPARK{TT}

    params::PT
    solver::ST
    iguess::IT
end

function IntegratorHPARK(equation::PDAE{DT,TT,FT,PT,UT,GT,ϕT,VT},
                         tableau::TableauHPARK{TT}, Δt::TT) where {DT,TT,FT,PT,UT,GT,ϕT,VT}
    D = equation.d
    S = tableau.s
    R = tableau.r

    N = 2*D*S + 3*D*R

    # create params
    params = ParametersHPARK{DT,TT,D,S,R,FT,PT,UT,GT,ϕT}(
                                equation.v, equation.f, equation.u, equation.g, equation.ϕ, Δt,
                                tableau.q, tableau.p, tableau.q̃, tableau.p̃, tableau.λ)

    # create solver
    solver = create_nonlinear_solver(DT, N, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create integrator
    IntegratorHPARK{DT, TT, typeof(equation), typeof(params), typeof(solver), typeof(iguess)}(
                                        equation, tableau, params, solver, iguess)
end

equation(int::IntegratorHPARK) = int.equation
timestep(int::IntegratorHPARK) = int.params.Δt
tableau(int::IntegratorHPARK) = int.tableau
nstages(int::IntegratorHPARK) = int.tableau.s
pstages(int::IntegratorHPARK) = int.tableau.r

function create_integrator_cache(int::IntegratorHPARK{DT,TT}) where {DT,TT}
    IntegratorCacheVPARK{DT, TT, ndims(int), nstages(int), pstages(int)}()
end

function update_params!(params::ParametersHPARK, cache::IntegratorCacheVPARK)
    # set time for nonlinear solver and copy previous solution
    params.t  = cache.t
    params.q .= cache.q
    params.p .= cache.p
    params.λ .= cache.λ
end

function initialize!(int::IntegratorHPARK, cache::IntegratorCacheVPARK)
    cache.t̅ = cache.t - timestep(int)

    equation(int).v(cache.t, cache.q, cache.p, cache.v)
    equation(int).f(cache.t, cache.q, cache.p, cache.f)

    initialize!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f,
                            cache.t̅, cache.q̅, cache.p̅, cache.v̅, cache.f̅)
end

function initial_guess!(int::IntegratorHPARK, cache::IntegratorCacheVPARK)
    for i in 1:nstages(int)
        evaluate!(int.iguess, cache.q, cache.p, cache.v, cache.f,
                              cache.q̅, cache.p̅, cache.v̅, cache.f̅,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q.c[i], tableau(int).p.c[i])

        for k in 1:ndims(int)
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - cache.q[k])/timestep(int)
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - cache.p[k])/timestep(int)
        end
    end

    for i in 1:pstages(int)
        evaluate!(int.iguess, cache.q, cache.p, cache.v, cache.f,
                              cache.q̅, cache.p̅, cache.v̅, cache.f̅,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in 1:ndims(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - cache.q[k])/timestep(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - cache.p[k])/timestep(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+3] = 0
        end
    end

    if int.params.t_λ.c[1] == 0
        for k in 1:ndims(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(k-1)+3] = cache.λ[k]
        end
    end
end

"Integrate DAE with variational partitioned additive Runge-Kutta integrator."
function integrate_step!(int::IntegratorHPARK{DT,TT}, cache::IntegratorCacheVPARK{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, cache)

    # compute initial guess
    initial_guess!(int, cache)

    # reset cache
    reset!(cache, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params, cache.n)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params, cache.n)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache, int.params)

    # compute final update
    update_solution!(cache.q, cache.Vi, int.params.t_q.b, timestep(int))
    update_solution!(cache.p, cache.Fi, int.params.t_p.b, timestep(int))

    # compute projection
    update_solution!(cache.q, cache.Up, int.params.t_q.β, timestep(int))
    update_solution!(cache.p, cache.Gp, int.params.t_p.β, timestep(int))
    # TODO # update_multiplier!(cache.λ, cache.Λp, int.params.t_λ.b)

    # copy solution to initial guess
    update!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f)

    # take care of periodic solutions
    cut_periodic_solution!(cache, equation(int).periodicity)
end


text:
#
# QWTwPlot module
# qwt - based 2D plotting
#

#__precompile__(true)
module QWTWPlot

using Printf
using Libdl

import Qt_jll
import qwtw_jll
import marble_jll

function __init__()
	# this is not OK  qwtwStart(Int64(0)) # start in normal mode
end

# DLLs and function handles below:
qwtwLibHandle = 0
qwtwFigureH = 0
qwtwClipGroupH = 0
qwtwRemoveLineH = 0
qwtwDebugMode = false
qwtwChangeLineH = 0
#qwtwFigure3DH = 0 # not supported because of license restrictions..    maybe later...
qwtwMapViewH = 0
qwtwsetimpstatusH = 0
qwtwCLearH = 0
qwtwPlotH = 0
qwtwPlot2H = 0
#qwtwPlot3DH = 0 # not supported because of license restrictions..    maybe later...
qwtwXlabelH = 0
qwtwYlabelH = 0
qwywTitleH = 0
qwtwVersionH = 0
qwtwMWShowH = 0
qwtEnableBroadcastH = 0

qwtMglH = 0
qwtMglLine = 0
qwtMglMesh = 0

qwtStartH = 0
qwtStartDebugH = 0
qwtStopH = 0
started = false

old_path = ""
old_qtPath = ""
oldLdLibPath = ""

"""
function saveEnv()
	saves some of the env vars
"""
function saveEnv()
	global old_path, old_qtPath, oldLdLibPath
	try old_path = ENV["PATH"]; catch end
	try old_qtPath = ENV["QT_PLUGIN_PATH"]; catch end
	try oldLdLibPath = ENV["LD_LIBRARY_PATH"]; catch end
	return
end

"""
function restoreEnv()
	restores env vars back
"""
function restoreEnv()
	global old_path, old_qtPath, oldLdLibPath
	ENV["PATH"] = old_path
	ENV["QT_PLUGIN_PATH"] = old_qtPath
	ENV["LD_LIBRARY_PATH"] = oldLdLibPath
	return
end

"""
function printEnv()
	simply print some of the env variables
"""
function printEnv()
	try
		@printf "\n\tPATH = %s \n\n" String(ENV["PATH"])
		@printf "\n\tQT_PLUGIN_PATH = %s \n\n" String(ENV["QT_PLUGIN_PATH"])
		try
			ldp = String(ENV["LD_LIBRARY_PATH"])
			@printf "\n\tLD_LIBRARY_PATH = %s \n\n" ldp
		catch
			@printf "no LD_LIBRARY_PATH\n\n"
		end
	catch
		@printf "printEnv: not everything was printed \n"
	end
end

"""
function addEnvItem(item, var)
	add something to the ENV from the beginning
	item -  what to add
	var - where to add
"""
function addEnvItem(item, var::String; debug = false)
	dv = ":"
	@static if Sys.iswindows()
		dv = ";"
	end
	useDV = true
	try # check: do we have it?
		test = ENV[var]  # yes
	catch ex # no such variable
		useDV = false
	end

	item2Add = "" # what we will actually add
	if typeof(item) == String # simple case
		item2Add = item
		if debug
			@printf "adding [%s] to [%s] \n" item2Add  var
		end
	elseif typeof(item) == Base.RefValue{String} # this can happen also
		item2Add = item[]
		if debug
			@printf "adding [%s] to [%s] \n" item2Add  var
		end
	else  # strange case
		@printf "WARNING: trying to add following item to %s : " var
		print(item)
		try
			item2Add = String(item)
		catch ex
			item2Add = item
		end
	end

	try
		if useDV
			ENV[var] = item2Add * dv * ENV[var]
		else
			ENV[var] = item2Add
			if debug
				@printf " new ENV %s created \n " var
			end
		end
	catch ex
		@printf "ERROR while adding item to [%s]: " var
		print(ex)
		print(item)
	end
end

"""
	qstart(;debug = false, qwtw_test = false, libraryName = "libqwtw")::Int32
starts the qwtw "C" library.

Without it, nothing will work. Call it before any other functions.\\
if debug, will try to enable debug print out. A lot of info. Just for debugging.\\
qwtw_test	-> if true, will try to not modify env variables. \\
"""
function qstart(;debug = false, qwtw_test = false, libraryName = "libqwtw")::Int32
	global qwtwDebugMode
	qwtwDebugMode = debug
	qwtw_libName = "nolib"
	if debug
		@printf "startint qwtw; current path: %s\n\n" ENV["PATH"]
	end

	# this could be still useful:
	#if debug qwtw_libName *= "d"; end

	global qwtwLibHandle, qwtwFigureH, qwtwMapViewH,  qwtwsetimpstatusH, qwtwCLearH, qwtwPlotH
	global qwtwPlot2H, qwtwXlabelH, qwtwYlabelH, qwywTitleH, qwtwVersionH, qwtwMWShowH, qwtwRemoveLineH
	global qwtEnableBroadcastH, qwtDisableBroadcastH, qwtwChangeLineH, qwtwClipGroupH
	#global qwtwPlot3DH, qwtwFigure3DH
	global qwtStartH, qwtStopH, started
	global old_path, old_qtPath, oldLdLibPath
	global qwtMglH, qwtMglLine, qwtMglMesh

	if started
		@printf "qwtw already started\n"
		return 0
	end

	@static if Sys.iswindows() #  this part will handle OS differences
		qwtw_libName = libraryName * ".dll"
	else # hopefully this may be Linux
		if debug
			@printf "\t non-Windows detected\n"
		end
		qwtw_libName = libraryName * ".so"
	end

	saveEnv()

	marbleDataPath = joinpath(marble_jll.artifact_dir, "data")
	marblePluginPath = joinpath(marble_jll.artifact_dir, "plugins")

	if qwtw_test
		@printf "qwtw debug; loading %s .. \n" qwtw_libName
		@printf "\nPATH = %s\n\n" ENV["PATH"]
	else
		@static if Sys.iswindows()
			ENV["QT_PLUGIN_PATH"]=Qt_jll.artifact_dir * "\\plugins"
			addEnvItem(qwtw_jll.PATH, "PATH", debug = debug)
			addEnvItem(qwtw_jll.LIBPATH, "PATH", debug = debug)

			#ENV["PATH"]= boost_jll.artifact_dir * "\\bin;" *  ENV["PATH"]
			#ENV["PATH"]= qwt_jll.artifact_dir * "\\bin;" *  ENV["PATH"]
			#ENV["PATH"]= Qt_jll.artifact_dir * "\\bin;" *  ENV["PATH"]

			#ENV["PATH"]= CompilerSupportLibraries_jll.artifact_dir * "\\bin;" *  ENV["PATH"]
			#new_env["PATH"]= FreeType2_jll.artifact_dir * "\\bin;" *  ENV["PATH"]
		else
			ENV["QT_PLUGIN_PATH"]=string(Qt_jll.artifact_dir) * "/plugins"
			addEnvItem(qwtw_jll.PATH, "PATH", debug = debug)
			addEnvItem(qwtw_jll.LIBPATH, "LD_LIBRARY_PATH", debug = debug)

			#ENV["PATH"]= boost_jll.artifact_dir * "/bin;" *  ENV["PATH"]

			#ENV["LD_LIBRARY_PATH"] = string(Qt_jll.LIBPATH) * ":" * ENV["LD_LIBRARY_PATH"] # do not sure why this is not already there
			#ENV["LD_LIBRARY_PATH"] = string(CompilerSupportLibraries_jll.LIBPATH) * ":" * ENV["LD_LIBRARY_PATH"] # do not sure why this is not already there
			#ENV["LD_LIBRARY_PATH"] = string(FreeType2_jll.LIBPATH) * ":" * ENV["LD_LIBRARY_PATH"] # do not sure why this is not already there
		end
	end

	if debug
		@printf "\nloading %s \n" qwtw_libName
		@printf "corrected ENV:\n"
		@static if Sys.iswindows()

		else
			@printf "LD_LIBRARY_PATH: %s \n\n" String(ENV["LD_LIBRARY_PATH"])
		end
		@printf "PATH: %s \n\n" String(ENV["PATH"])
		@printf "\n"
	end

	try
		qwtwLibHandle = Libdl.dlopen(qwtw_libName)
	catch ex
		if debug
			printEnv();
		end
		restoreEnv()
		@printf "Sorry, dlopen for %s failed; something is wrong\n" qwtw_libName
		throw(ex)
	end

	if debug
		@printf "\nlibrary %s opened from %s \n" qwtw_libName  Libdl.dlpath(qwtwLibHandle)
	end
	qwtwFigureH = Libdl.dlsym(qwtwLibHandle, "qwtfigure")
	#qwtwFigure3DH = Libdl.dlsym(qwtwLibHandle, "qwtfigure3d")

	try
		qwtwMapViewH = Libdl.dlsym(qwtwLibHandle, "qwtmap")
	catch
		qwtwMapViewH = 0
		@printf "WARNING: topview functions disabled (looks like no [marble] support)\n"
	end

	qwtwsetimpstatusH = Libdl.dlsym(qwtwLibHandle, "qwtsetimpstatus")
	qwtwCLearH = Libdl.dlsym(qwtwLibHandle, "qwtclear")
	qwtwPlotH = Libdl.dlsym(qwtwLibHandle, "qwtplot")
	#qwtwPlot3DH = Libdl.dlsym(qwtwLibHandle, "qwtplot3d")
	qwtwPlot2H = Libdl.dlsym(qwtwLibHandle, "qwtplot2")
	qwtwXlabelH = Libdl.dlsym(qwtwLibHandle, "qwtxlabel")
	qwtwYlabelH = Libdl.dlsym(qwtwLibHandle, "qwtylabel")
	qwywTitleH = Libdl.dlsym(qwtwLibHandle, "qwttitle")
	qwtwVersionH = Libdl.dlsym(qwtwLibHandle, "qwtversion")
	qwtwMWShowH = Libdl.dlsym(qwtwLibHandle, "qwtshowmw")
	qwtStartH = Libdl.dlsym(qwtwLibHandle, "qtstart")
	qwtStartDebugH = Libdl.dlsym(qwtwLibHandle, "qtstart_debug")
	qwtStopH = Libdl.dlsym(qwtwLibHandle, "qwtclose")
	qwtwRemoveLineH = Libdl.dlsym(qwtwLibHandle, "qwtremove")
	qwtwChangeLineH = Libdl.dlsym(qwtwLibHandle, "qwtchange")

	try
		qwtwClipGroupH = Libdl.dlsym(qwtwLibHandle, "qwtclipgroup")
	catch
		@printf "WARNING: clip groups disabled (wrong qwtw version?) \n"
	end

	try
		qwtEnableBroadcastH = Libdl.dlsym(qwtwLibHandle, "qwtEnableCoordBroadcast")
		qwtDisableBroadcastH = Libdl.dlsym(qwtwLibHandle, "qwtDisableCoordBroadcast")
	catch
		@printf "WARNING: UDP broacast disabled \n"
	end

	try
		qwtMglH = Libdl.dlsym(qwtwLibHandle, "qwtmgl")
		qwtMglLine = Libdl.dlsym(qwtwLibHandle, "qwtmgl_line")
		qwtMglMesh = Libdl.dlsym(qwtwLibHandle, "qwtmgl_mesh")
	catch
		@printf "WARNING: 3D features disabled \n"
	end

	# hangs tight!!!!
	if debug
		@printf "starting qstart in debug mode \n"
		@printf "\tmarbleDataPath %s \n" marbleDataPath
		@printf "\tmarblePluginPath %s \n" marblePluginPath
	end
	test = 1
	try
		if debug
			test = ccall(qwtStartDebugH, Int32, (Ptr{UInt8}, Ptr{UInt8}, Int32),
				marbleDataPath, marblePluginPath, 10);
		else
			test = ccall(qwtStartH, Int32, (Ptr{UInt8}, Ptr{UInt8}),
				marbleDataPath, marblePluginPath); # very important to call this in the very beginning
		end
	catch ex
		if debug
			printEnv();
		end
		restoreEnv()
		@printf "Sorry, cannot start qwtw. Something is wrong\n"
		throw(ex)
	end
	#end
	#@printf "qtstart = %d \n" test

	#version = qversion();
	#println(version);
	if test == 0
		started = true
		if debug
			version = qversion();
			println(version);
		end
	else
		@printf "\nERROR: library was not started (return code %d) \n" test
		if debug
			printEnv();
		end
	end

	restoreEnv()

	return test
end

"""
	qstop()
close everything and detach from qwtw library.

It maybe useful for debugging. What it does actually, it sends a command to the QT process to exit.\\
Have to call `qstart()` before, though.
"""
function qstop()
	global qwtwLibHandle, qwtStopH
	global started

	if qwtwLibHandle != 0
		if (qwtStopH != 0)
			ccall(qwtStopH,  Cvoid, ())
			#@printf "qwtw stopped\n"
		else
			@printf "error: was not started correctly\n"
		end
		Libdl.dlclose(qwtwLibHandle)
	else
		@printf "not started (was qstart() called?)\n"
		return
	end
	qwtwLibHandle = 0
	started = false
	return
end

"""
	qversion() :: String

useful for debugging.\\
return version info (as a string);
"""
function qversion() :: String
	global qwtwVersionH, qwtwLibHandle
	if qwtwLibHandle == 0
		return "not started yet"
	end
	v = zeros(UInt8, 1024)
	#@printf "qwtwc version: "
	bs = ccall(qwtwVersionH, Int32, (Ptr{UInt8}, Int32), v, 1024);
	#return bytestring(pointer(v))
	cmd = unsafe_string(pointer(v), bs);
	#@printf " bs = %d .... " bs
	return cmd
end

# looks like not used anymore
function traceit( msg )
     # global g_bTraceOn
      if ( true )
         bt = backtrace() ;
         s = sprint(io->Base.show_backtrace(io, bt))
         println( "debug: $s: $msg" )
      end
end

"""
	qfigure(n)

create a new plot window, OR make plot (with this ID `n`) 'active'.

looks like `n` have to be an integer number (this plot ID, or zero if you do not care).\\
after this function, this plot is the "active plot". \\
If `n == 0` then another new plot will be created.
xAxisType and yAxisType could be :aLinear (default) or :aLog (logarithmic)

Returns ID of the created plot
"""
function qfigure(n::Integer = 0; xAxisType=:aLinear, yAxisType=:aLinear)::Int32
	global qwtwFigureH, started
	if !started
		@printf "not started (was qstart() called?)\n"
		return 0
	end
	flags::UInt32 = 0
	if xAxisType == :aLog
		flags |= 1
	end
	if yAxisType == :aLog
		flags |= 2
	end
	test = ccall(qwtwFigureH, Int32, (Int32, UInt32), n, flags);
	return test
end;

"""
	qclipgrp(gr)

set current 'clip group'.

gr is the ID of this new group.
all (new) subsequent plots/figires will belong to this group. if press 'clip' button,
then if will work only for plots from the same group.

useful in case some of the plots have another time range.
"""
function qclipgrp(gr)
	global qwtwClipGroupH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end
	if qwtwClipGroupH == 0
		@printf "clip groups are not supported, sorry\n"
		return
	end

	ccall(qwtwClipGroupH, Cvoid, (Int32,), Int32(gr));
end
export qclipgrp

"""
	qremove(id::Int32)

remove a line from a plot.
"""
function qremove(id::Int32)
	global qwtwRemoveLineH, started
	if !started
		@printf "not started (was qstart() called?)\n"
		return
	end
	ccall(qwtwRemoveLineH, Cvoid, (Int32,), id);
	return
end
export qremove

"""
	qfmap(n)

create a new  plot window to draw on a map (with specific window ID)\\
'n' is Int32
"""
function qfmap(n)::Int32
	global qwtwMapViewH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return 0
	end
	if qwtwMapViewH == 0
		@printf "map vew not supported, sorry\n"
		return 0
	end

	test = ccall(qwtwMapViewH, Int32, (Int32,), n);
	return test;
end;

"""
	qfmap()

create a new  plot window to draw on a map.
"""
function qfmap()::Int32
	return qfmap(0)
end;

"""
	qmgl(n = 0)

create a new  plot window to draw 3D lines / surfaces.

currently works only for Linux (?)
"""
function qmgl(n = 0)::Int32
	global qwtMglH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return 0
	end
	if qwtMglH == 0
		@printf "3D not supported, sorry\n"
		return 0
	end

	test = ccall(qwtMglH, Int32, (Int32,), n);
	return test
end;

export qmgl

"""
	qmgline(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, style::String = "-sb"; name = "")

draw a 3D line.

x, y, and z are the vectors with point coordinates.\\
style - how to draw a line. Explanation is here: http://mathgl.sourceforge.net/doc_en/Line-styles.html\\
name - a legend for this line
"""
function qmgline(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, style::String = "-sb"; name = "")
	global qwtMglLine, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end
	if qwtMglLine == 0
		@printf "3D not supported, sorry\n"
		return
	end
	if length(x) < 1
		@printf "QWTWPlot::qmgline empty X value\n"
		return
	end
	n = length(x)
	@assert length(x) == length(y)
	@assert length(z) == length(x)

	try
	ccall(qwtMglLine, Cvoid, (Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{UInt8}, Ptr{UInt8}),
		n, x, y, z, name, style);

		sleep(0.025)
	catch
		@printf "ERROR in QWTWPlot::qmgline\n"
	end
	return
end
export qmgline

"""
	qmglmesh(data::Array{Float64, 2}, xMin = -10.0, xMax = 10.0, yMin = -10.0, yMax = 10.0,  style::String = ""; name = "", type = 0)

draw a 3D mesh/surface.   Currently works only for Linux (is this true?). \\

data - vertical (z) coordinates of the surface. should be Array{Float64, 2}\\
xMin, xMax, yMin, yMax is the definition of the range where to draw a surface. All the
coordinates from 'data' are evenly distributed inside this range.\\
style how to draw the surface. Some hints could be taken from here: http://mathgl.sourceforge.net/doc_en/Color-scheme.html \\
type  -   0 or 1;  mesh/grid  or surface\\
name - not used yet (maybe will add later)

"""
function qmglmesh(data::Array{Float64, 2}, xMin = -10.0, xMax = 10.0, yMin = -10.0, yMax = 10.0,  style::String = ""; name = "", type = 0)
	global qwtMglMesh, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end
	if qwtMglMesh == 0
		@printf "3D not supported, sorry\n"
		return
	end
	if length(data) < 1
		@printf "QWTWPlot::qmgline empty X value\n"
		return
	end
	d = size(data)
	if (d[1] < 1) || (d[2] < 1)
		@printf "QWTWPlot::qmglmesh empty data\n"
		return
	end
	xSize = d[1]
	ySize = d[2]

	try
	ccall(qwtMglMesh, Cvoid, (Int32, Int32,
		Float64, Float64, Float64, Float64,
		Ptr{Float64},
		Ptr{UInt8}, Ptr{UInt8}, Int32),
		xSize, ySize,
		xMin, xMax, yMin, yMax, data,  name, style, type);

		sleep(0.025)
	catch
		@printf "ERROR in QWTWPlot::qwtMglMesh\n"
	end
	return
end

export qmglmesh

# create a new  window to draw a 3D points (QT engine)
#function qf3d(n)
#	global qwtwFigure3DH
#	ccall(qwtwFigure3DH, Cvoid, (Int32,), n);
#end;

"""
	qimportant(i)

set up an importance status for next lines.

looks like `0` means 'not important', `1` means "important.\\
'not important' will not participate in 'clipping':\\
	'not important' lines may be not completely inside the window.
"""
function qimportant(i)
	global qwtwsetimpstatusH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtwsetimpstatusH, Cvoid, (Int32,), i);
	return
end

"""
	qclear()

close all the plots
"""
function qclear()
	global qwtwCLearH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtwCLearH, Cvoid, ());
	return
end

"""
	qsmw()

open/show "main control window".
"""
function qsmw()
	global qwtwMWShowH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtwMWShowH, Cvoid, ());
	return
end

"""
qplot(x::Vector{Float64}, y::Vector{Float64}, name::String, style::String, lineWidth, symSize):: Int32

plot normal lines.

#Parameters:

x and y:   the points.\\

name:      	name for this line.\\

style: 	 	how to draw a line.\\

lineWidth:	is a line width.\\

symSize:	size of the symbols, if they are used in 'style' spec.\\

what does 'style' parameter means? It's a string which has 1 or 2 or 3 symbols. \\

Look at two places for the detail explanation:\\
* example code  https://github.com/ig-or/QWTWPlot.jl/blob/master/src/qwexample.jl
* https://github.com/ig-or/QWTWPlot.jl/blob/master/docs/line-styles.md

this function returns ID of the line created. hopefully you can use it in some other functions\\
ID supposed to be >=0 is all is OK, or <0 in case of error
"""
function qplot(x::Vector{Float64}, y::Vector{Float64}, name::String, style::String="-b", lineWidth=1, symSize=1) :: Int32
	global qwtwPlotH, qwtwLibHandle, started, qwtwDebugMode
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return -88
	end

	if length(x) != length(y)
		@printf "qplot: x[%d], y[%d]\n" length(x) length(y)
		traceit("error")
		return -89
	end
	#@assert (length(x) == length(y))

	n = length(x)
	if n < 1
		@printf "QWTWPlot qplot: bad input array length (%d) \n" n
		return 578
	end
	ww::Int32 = lineWidth;
	s::Int32 = symSize
	test::Int32 = -50

	if qwtwDebugMode
		@printf "qplot: (%s) (%s) n=%d \n " name style n
	end
	try
		test = ccall(qwtwPlotH, Int32, (Ptr{Float64}, Ptr{Float64}, Int32, Ptr{UInt8}, Ptr{UInt8}, Int32, Int32),
			x, y, n, name, style, ww, s);
		sleep(0.025)
	catch
		@printf "qplot: error #2  n = %d;  name = %s style = %s\n" n name style
		traceit("error #2")
		return -42
	end
	return test
end

"""
qplot1(x::Vector{Float64}, y::Vector{Float64}, name::String, style::String, symSize):: Int32

plot lines with line width == 1.

```

x and y:   the points
name:      	name for this line
style: 	 	how to draw a line
symSize:	size of the symbols

```

what does 'style' parameter means? It's a string which has 1 or 2 or 3 symbols. \\
Look at two places for the detail explanation:

* example code  https://github.com/ig-or/QWTWPlot.jl/blob/master/src/qwexample.jl
* https://github.com/ig-or/QWTWPlot.jl/blob/master/docs/line-styles.md

this function returns ID of the line created. hopefully you can use it in some other functions\\
ID supposed to be >=0 is all is OK, or <0 in case of error

"""
function qplot1(x::Vector{Float64}, y::Vector{Float64}, name::String, style::String, symSize) :: Int32
	global qwtwPlotH, qwtwLibHandle, started, qwtwDebugMode
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return -88
	end

	n = length(x)
	n1 = length(y)
	#@assert n > 0
	#@assert n1 > 0
	#@assert (n == n1)
	if (n < 1) || (n1 < 1) || (n != n1)
		@printf "QWTWPlot qplot: bad input array length (%d) \n" n
		error("qplot1: wrong array length")
		return -89
	end
	ww::Int32 = symSize;
	test::Int32 = -50
	if qwtwDebugMode
		@printf "qplot1: (%s) (%s) n=%d \n " name style n
	end

	try
		test = ccall(qwtwPlotH, Int32, (Ptr{Float64}, Ptr{Float64}, Int32, Ptr{UInt8}, Ptr{UInt8}, Int32, Int32),
			x, y, n, name, style, 1, ww);
	catch ex
		@printf "qplot1 ERROR\n"
		throw(ex)
		return -42
	end
	sleep(0.025)
	return test
end

"""
qchange(id::Int32, x::Vector{Float64}, y::Vector{Float64}, t::Vector{Float64} = []) :: Int32

change existing line.
"""
function qchange(id::Int32, x::Vector{Float64}, y::Vector{Float64}, t::Vector{Float64} = Vector{Float64}([])) :: Int32
	global qwtwChangeLineH, started
	if !started
		@printf "not started (was qstart() called?)\n"
		return
	end
	n = length(x)
	n1 = length(y)
	n2 = length(t)
	if (n < 1) || (n1 < 1) || (n != n1)
		@printf "QWTWPlot qchange: bad input array length (%d) \n" n
		error("qchange: wrong array length")
		return -89
	end

	zz = Ptr{Float64}(C_NULL)
	tt = Ptr{Float64}(C_NULL)
	if n2 > 0
		tt = t
		if n2 != n
			error("qchange: wrong _t_ array length")
			return 84
		end
	end
	test = -50
	try
		test = ccall(qwtwChangeLineH, Int32, (Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32),
			id, x, y, zz, tt, n);
	catch ex
		@printf "qchange ERROR\n"
		throw(ex)
		return -42
	end
	return test
end
export qchange

# draw symbols in 3D space
# currently style and 'w' are not used
#function qplot3d(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64},
#	 		name::String, style::String, w, time::Vector{Float64})
#	global qwtwPlotH
#	@assert (length(x) == length(y))
#	n = length(x)
#	ww::Int32 = w;
#	ccall(qwtwPlot3DH, Cvoid, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
#			Int32, Ptr{UInt8}, Ptr{UInt8}, Int32, Int32,  Ptr{Float64}),
#		x, y, z, n, name, style, 1, ww, time);
#	sleep(0.025)
#
#end;
#

"""
	qEnableCoordBroadcast(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, time::Vector{Float64})

Enable UDP client.\\

In case you'd like to connect to some external software from this library,
it has UDP client and server inside. If/when you will move a marker, it will
broadcast 'marker information' via UDP. \\
Other way should also work: if some other application will broadcast UDP with marker info,
all the markers in this library supposed to be updated.

parameters: x, y,  and z are the vectors describing the line in 3D space;
and 'time' is an additional time information, for every point of the line.

It will send out marker info to UDP port 49561.\\
Incoming (to this library) UDP  port number is 49562, and IP address is "127.0.0.1".

Incoming message format:  'CRDS"<x><y><z>"FFFF', where
<x>, <y>, and <z> are  Float64 point coordinates, 8 bytes each,
so this message contains 4 + 3*8 + 4 bytes\\
outgoing message format is approximately the same: "EEEE"<x><y><z>"FFFF"

If UDP client is enabled, server will also start. For all this to work, some firewall rules have to be added probably,
line 'allow 49561 and 49562 ports for the local host'.
"""
function qEnableCoordBroadcast(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64},
	 		time::Vector{Float64})
	global qwtEnableBroadcastH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	@assert (length(x) == length(y))
	@assert (length(x) == length(z))
	@assert (length(x) == length(time))

	n = length(time)

	if n < 1
		@printf "QWTWPlot qEnableCoordBroadcast: bad input array length (%d) \n" n
	end

	ccall(qwtEnableBroadcastH, Cvoid, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
			 Ptr{Float64}, Int32),
		x, y, z , time, n);
	sleep(0.025)
	return
end;

"""
	qDisableCoordBroadcast()
disable all this UDP -related stuff.\\

This function will stop UDP server and client.
"""
function qDisableCoordBroadcast()
	global qwtDisableBroadcastH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtDisableBroadcastH, Cvoid, ());
	sleep(0.025);
	return
end;

"""
	qplot2(x::Array{Float64}, y::Array{Float64}, time::Array{Float64}, name::String, style::String, lineWidth=1, symSize=1):: Int32

plot with additional parameter (time?) info.

## Parameters:

```

'x', 'y': the line
'name' name of this line
'style'  the line style
'lineWidth' line width
'symSize' symbol size;
'time' time info, for every point

```

## Example

```

julia-repl
julia> time=collect(0.0:0.01:10.0)
x = sin.(time)
y = cos.(time)
qfigure()
qplot(time, x + y, "function 1", "-b", 3)
qfigure()
qplot2(x, y, time, "function 2", "-m", 3)

```
Now use marker on both plots and see that it moves on both plots.

this function returns ID of the line created. hopefully you can use it in some other functions\\
ID supposed to be >=0 is all is OK, or <0 in case of error

"""
function qplot2(x::Array{Float64}, y::Array{Float64}, time::Array{Float64}, name::String, style::String, lineWidth=1, symSize=1):: Int32
	global qwtwPlot2H, qwtwLibHandle, started, qwtwDebugMode
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return -88
	end
	if ((length(x) != length(y)) || (length(y) != length(time)))
		@printf "QWTWPlot qplot2 error: different array length \n"
		return -75
	end

	n = length(x)
	if n < 1
		@printf "QWTWPlot qplot2 error: bad array length \n"
		return -76
	end
	ww::Int32 = lineWidth;
	s::Int32 = symSize;
	test::Int32 = -50
	if qwtwDebugMode
		@printf "qplot2: (%s) (%s) n=%d \n " name style n
	end
	try
		test = ccall(qwtwPlot2H, Int32, (Ptr{Float64}, Ptr{Float64}, Int32, Ptr{UInt8}, Ptr{UInt8}, Int32, Int32, Ptr{Float64}),
			x, y, n, name, style, ww, s, time);
	catch ex
		@printf "qplot1 ERROR\n"
		throw(ex)
		return -42
	end

	sleep(0.025)
	return test
end

"""
	qxlabel(s::String)
put a label on the horizontal axis on the bottom.
"""
function qxlabel(s::String)
	global qwtwXlabelH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtwXlabelH, Cvoid, (Ptr{UInt8},), s);
	return
end;

"""
	qylabel(s::String)
put a label on the left vertical axis
"""
function qylabel(s::String)
	global qwtwYlabelH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtwYlabelH, Cvoid, (Ptr{UInt8},), s);
	return
end;

"""
	qtitle(s::String)
put a title on the current plot.
"""
function qtitle(s::String)
	global qwywTitleH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwywTitleH, Cvoid, (Ptr{UInt8},), s);
	return
end

"""
	qStarted()
does it started or not yet.
"""
function qStarted()
	global started
	return started
end

export qfigure, qfmap, qplot, qplot1, qplot2, qxlabel,  qylabel, qtitle
export qimportant, qclear, qstart, qstop, qversion, qsmw
export traceit
export  qEnableCoordBroadcast, qDisableCoordBroadcast
export qStarted
#export qplot3d, qf3d,

end # module


text:
abstract type Basis end

Base.size(basis::Basis)::Int64 = basis.o.size
statistics(basis::Basis)::Statistics = basis.o.statistics

const SVEResultType = Tuple{PiecewiseLegendrePoly,Vector{Float64},PiecewiseLegendrePoly}

"""
IRBasis
"""

#==
struct IRBasis <: Basis
    u::PiecewiseLegendrePoly
    uhat::PiecewiseLegendreFT
    v::PiecewiseLegendrePoly
    s::Vector{Float64}
    statistics::Statistics
    size::Int64
end

IRBasis(o::PyObject) = IRBasis(
        o, o.u, o.uhat, o.v, o.s,
        o.statistics == "F" ? fermion : boson,
        o.size
    )

const SVEResultType = Tuple{Matrix{Float64}, Vector{Float64}, Matrix{Float64}}

function get_kernel(statistics::Statistics, lambda::Float64, kernel::Union{Kernel,Nothing})
    if kernel === nothing
        kernel = LogisticKernel(lambda)
    else
        if hasproperty(kernel, :lambda)
            if kernel.lambda != lambda
                throw(RuntimeError("lambda mismatch!"))
            end
        end
    end
    kernel
end

function IRBasis(
    statistics::Statistics, lambda::Float64, eps=None;
    kernel::Union{Kernel,Nothing}=nothing,
    sve_result::Union{SVEResultType,Nothing}=nothing)

    kernel = get_kernel(statistics, lambda, kernel)
    if sve_result === nothing
        u, s, v = pysve.compute(self.kernel.o, eps)
    else
        u, s, v = sve_result
    end

    _even_odd = statistics == fermion ? (:odd) : (:even)

    IRBasis(
        u, hat(uhat, _even_odd, n_asymp=conv_radius(kernel)),
        v, s, statistics
    )
end
==#

"""
FiniteTempBasis
"""
struct FiniteTempBasis <: Basis
    beta::Float64
    o::PyObject
    u::PiecewiseLegendrePoly
    uhat::PiecewiseLegendreFT
    v::PiecewiseLegendrePoly
    s::Vector{Float64}
    statistics::Statistics
    size::Int64
end

function FiniteTempBasis(o::PyObject)
    return FiniteTempBasis(o.beta, o, o.u, o.uhat, o.v, o.s,
                           o.statistics == "F" ? fermion : boson,
                           o.size)
end

"""
Create a FiniteTempBasis object by decomposing a given kernel
"""
function FiniteTempBasis(kernel::KernelBase, statistics::Statistics, beta::Real,
                         eps::Union{Float64,Nothing}=nothing;
                         sve_result::Union{SVEResultType,Nothing}=nothing)
    if sve_result !== nothing
        sve_result = (sve_result[1].o, sve_result[2], sve_result[3].o)
    end
    o = sparse_ir.FiniteTempBasis(statistics == fermion ? "F" : "B", Float64(beta),
                                  Float64(kernel.o.lambda_ / beta); eps=eps,
                                  kernel=kernel.o, sve_result=sve_result)
    return FiniteTempBasis(o)
end

function FiniteTempBasis(statistics::Statistics, beta::Real, wmax::Real,
                         eps::Union{Float64,Nothing}=nothing
                         ;
                         kernel::Union{KernelBase,Nothing}=nothing,
                         sve_result::Union{SVEResultType,Nothing}=nothing)
    kernel_py = kernel !== nothing ? kernel.o : nothing
    stat = statistics == fermion ? "F" : "B"
    if eps !== nothing
        eps > 1.6e-8 || error("xprec is not supported yet.")
    end
    if sve_result !== nothing
        sve_result = (sve_result[1].o, sve_result[2], sve_result[3].o)
    end
    o = sparse_ir.FiniteTempBasis(stat, Float64(beta), Float64(wmax); eps=eps,
                                  kernel=kernel_py, sve_result=sve_result)
    return FiniteTempBasis(o)
end

"""
Create a FiniteTempBasis object using LogisticKernel
"""
function FiniteTempBasis(beta::Real,
                         wmax::Real,
                         statistics::Statistics,
                         eps::Union{Float64,Nothing})
    return FiniteTempBasis(LogisticKernel(Float64(beta * wmax)), statistics, Float64(beta),
                           eps)
end

"""
Construct FiniteTempBasis objects for fermion and bosons using
the same LogisticKernel instance.
"""
function finite_temp_bases(beta::Real, wmax::Real, eps::Real,
                           sve_result::Union{SVEResultType,Nothing}=nothing)::Tuple{FiniteTempBasis,
                                                                                FiniteTempBasis}
    if sve_result !== nothing
        sve_result = (sve_result[1].o, sve_result[2], sve_result[3].o)
    end
    basis_f, basis_b = sparse_ir.basis.finite_temp_bases(Float64(beta),
                                                         Float64(wmax),
                                                         Float64(eps),
                                                         sve_result)
    return (FiniteTempBasis(basis_f), FiniteTempBasis(basis_b))
end

function default_omega_sampling_points(basis::FiniteTempBasis)::Vector{Float64}
    return basis.o.default_omega_sampling_points()
end

function default_tau_sampling_points(basis::FiniteTempBasis)::Vector{Float64}
    return basis.o.default_tau_sampling_points()
end

function default_matsubara_sampling_points(basis::FiniteTempBasis)::Vector{Int64}
    return basis.o.default_matsubara_sampling_points()
end

function Base.getproperty(basis::FiniteTempBasis, d::Symbol)
    if d === :sve_result
        u_py, s::Vector{Float64}, v_py = getfield(basis, :o).sve_result
        u = PiecewiseLegendrePoly(u_py)
        v = PiecewiseLegendrePoly(v_py)
        return u, s, v
    else
        return getfield(basis, d)
    end
end

function Base.propertynames(basis::FiniteTempBasis, private::Bool=false)
    return (:sve_result,
            fieldnames(FiniteTempBasis)...)
end


text:
# TODO:

# terms.jl:
# How to handle adjoints? Most of the time, the terms will be Hermitian, but some of the
# times, they are not. In some cases it works well (e.g. `MomentumTwoBodyTerm`), but adding
# the option for all terms would be annoying at best.
# Perhaps the "functions" you pass to the terms could know whether permuting their arguments
# changes anything - this way the terms would know what structure they have.

# sum.jl:
# Idea: summing should merge like terms (e.g. summing a harmonic potential with a delta
# potential should be able to merge the `FullOneBodyTerm`s.

# abstract.jl:
# Idea: perhaps the type hirearchy should be reworked a bit. These terms don't really need
# to know the starting address. Hamiltonians could still be required to provide them.
# A potential (or more general Hamiltonian addon) could inherit the starting address from
# the Hamiltonian it's being added to.

# abstract.jl:
# Idea: add IsDiagonal to LOStructure. Then they can be skipped when generating
# offdiagonals. Another option would be to add IsSymmetric, which would provide an adjoint
# implementation.

module NewRimuOperators

using StaticArrays
using LinearAlgebra

using Rimu
# Import some Rimu internals:
using Rimu.Hamiltonians: momentum_transfer_excitation
using Rimu.Hamiltonians: momentum_transfer_diagonal
using Rimu.Hamiltonians: momentum_external_potential_excitation
using Rimu.Hamiltonians: momentum_external_potential_diagonal
using Rimu.Hamiltonians: momentum_space_harmonic_potential
using Rimu.Hamiltonians: AbstractOffdiagonals
using Rimu.BitStringAddresses: update_component
# These will be extended
import Rimu: num_offdiagonals, get_offdiagonal, diagonal_element, offdiagonals, starting_address, LOStructure

export ConstFunction, InteractionMatrix
export column
export Hubbard, Transcorrelated
export HarmonicPotential, DeltaPotential, TranscorrelatedDeltaPotential

include("utilities.jl")
include("abstract.jl")
include("column.jl")
include("sum.jl")
include("terms.jl")
include("extensions.jl")

include("hubbard.jl")
include("transcorrelated.jl")

include("potentials.jl")

end # module


text:
# Needed for sampling discrete probability distributions efficiently
using StatsBase

# Uniform sample from the d-simplex; optional 2nd argument makes 1st element 1-eps.
randp(d::Integer)::Vector{Float64} = diff([0; sort(rand(d-1)); 1])
randp(d::Integer,ε::Float64)::Vector{Float64} = [(1-ε); ε.*randp(d-1)]

# Define a random Pauli channel on n qubits with support size at most k and
# weight 1-𝛿 on the identity Pauli. Returns a tuple (p,P).
function RandomPauliChannel(n::Integer,k::Integer,𝛿::Float64)::Tuple{Vector{Float64}, Matrix{Int8}}
    P = unique(vcat(zeros(Int8,1,n),rand(Int8.(0:3),k-1,n)),dims=1)
    (randp(size(P)[1],𝛿), P)
end

# Sample the channel given by (p,P) m times and return the associated Paulis
function SampleChannel(p::Vector{Float64},P::Matrix{Int8},m::Integer)::Matrix{Int8}
    (k,n) = size(P)
    w = Weights(p)
    s = sample(1:k, w, m) # sampling
    [ P[a,b] for a in s, b in 1:n ]
end

# This is the star product from the paper, aka the symplectic product.
# If A, B = {0,1,2,3} are single-qubit Paulis, then star(A,B) == 0 iff A and B commute.
star(A::Matrix{Int8},B::Matrix{Int8})::Matrix{Int8} = ((A .>> 1) .& B) .⊻ ((B .>> 1) .& A)

# Given the channel (p,P) and a list of states A, record the measurement outcomes R
function Measure(A::Matrix{Int8},p::Vector{Float64},P::Matrix{Int8})::Matrix{Int8}
    m = size(A)[1]
    C = SampleChannel(p,P,m)
    star(A,C)
end

# Given a list of measurement outcomes R, and a list A of "B adjusted" inputs,
# estimate the probability of all-0 using the estimator with lower variance.
# Here w[k] = (-1/2)^(k-1) / m is the weighing factor for the estimator.
# This estimator assumes B is not the all-0 vector already.
function EstProb(B::Vector{Int8},A::Matrix{Int8},R::Matrix{Int8},w::Vector{Float64},t::Integer)::Float64
    if all(B .== 0)
        return w' * (counts(sum(R, dims = 2),0:t))
    end
    AB = star(A, B' .+ zeros(Int8,size(A)) )
    w' * (counts(sum(AB .⊻  R, dims = 2),0:t) - counts(sum(AB, dims = 2),0:t) )
end

# Estimate the whole channel using branch and prune from probe states A
# and measurement record R, with 𝛿 the pruning threshold.
# This implementation has a bad runtime of O(m*n^2/𝛿). We could save a factor of n
# by recycling the marginal estimates from round j and using them in round j+1.
function EstChan(A::Matrix{Int8},R::Matrix{Int8},𝛿::Float64)::Tuple{Vector{Float64}, Matrix{Int8}}
    (m,n) = size(A)
    p = Float64[]
    P = Vector{Int8}[[]]
    w = (-1/2).^(0:n) / m
    for j = 1:n
        q = Float64[]
        Q = Array{Int8,1}[]
        AA = A[:,1:j]
        RR = R[:,1:j]
        ww = w[1:(j+1)]
        for s = 1:length(P)
            prefx = P[s]
            for c = Int8.(0:3)
                B = vcat(prefx, c)
                est = EstProb(B, AA, RR, ww, j)
                if est > 𝛿
                    q = vcat(q, est)
                    push!(Q,B)
                end
            end
        end
        p = q
        P = Q
    end
    (p,vcat(P'...))
end

# Estimate 1 std error bars using the jackknife estimator for variance
function EstChanErrorBars(A::Matrix{Int8},Q::Matrix{Int8},R::Matrix{Int8})::Vector{Float64}
    k = size(Q)[1]
    (m,n) = size(A)
    s = zeros(k)
    for j = 1:k
        AQ = star(A, Q[j,:]' .+ zeros(Int8, (m,n)) )
        s[j] = std(dropdims((-1/2).^sum(AQ .⊻  R,dims=2) - (-1/2).^sum(AQ,dims=2),dims=2))
    end
    s/sqrt(m)
end

# represent a sparse Pauli channel as a dictionary
sparseP(p::Vector{Float64},P::Matrix{Int8})::Dict{Vector{Int8},Float64} =
    Dict{Vector{Int8},Float64}( [P[j,:] for j in 1:size(P)[1] ] .=> p )

# compute the total variation distance between two Pauli dictionaries
tvd(pP::Dict{Vector{Int8},Float64},qQ::Dict{Vector{Int8},Float64})::Float64 =
    sum(abs,values(merge(-,pP,qQ)))/2

# Brute force compute the estimator on *all* the Paulis. Not the same as EstChan(A,R,0.0).
# It only keeps the final results with positive estimates. Don't use with large n!
function EstAll(A::Matrix{Int8},R::Matrix{Int8})::Tuple{Vector{Float64}, Matrix{Int8}}
    (m,n) = size(A)
    @assert n < 6
    w = (-1/2).^(0:n) / m
    p0 = w' * counts(sum(R, dims = 2),0:n)
    s = pushfirst!([EstProb(reverse(digits(Int8, j,base=4,pad=n)),A,R,w,n) for j = 1:4^n-1], p0)
    S = vcat(reverse.(digits.(Int8, findall(x -> x > 0, s) .- 1, base=4,pad=n))'...)
    s = s[s .> 0]
    (s,S)
end

max_stars_repo_name: sflammia/PauliPopRec

max_issues_repo_name: sflammia/PauliPopRec

max_forks_repo_name: sflammia/PauliPopRec


text:
function DFS!(G)
	v = 1:size(G.AdjMatrix, 1);
	for u = v
		G.V[u].color = 'W';
		G.V[u].pred = -1;
	end
	for u = v
		if G.V[u].color == 'W';
			DFS_Visit!(G, u);
		end
	end
end
function DFS_Visit!(G, u)
	G.V[u].color = 'G';
        for v in findall(G.AdjMatrix[u,:] .==1)
		if G.V[v].color == 'W'
			G.V[v].pred = u;
			DFS_Visit!(G, v);
		end
	end
	G.V[u].color = 'B';
end


text:
module OpenTelemetryAPI

include("common.jl")
include("propagator.jl")
include("trace/trace.jl")
include("metric/metric.jl")
include("log.jl")

end # module