text:
 # ESA PVIR specific stuff

include("itmix-setup.jl")

sd = "esa-14nov"
data, data_summarized = FileIO.load("data/ITMIX2/results_Werder/$sd/data.jld", "data", "data_summarized");

if !isdefined(:gls)
    gls = open("data/ITMIX2/results_Werder/$sd/gls.jls", "r") do io
        deserialize(io)
    end
end
runit = false
include("itmix2-evaluate.jl")

using Mustache
"""
Turn a DataFrame into a latex table.  Copied from Mustache.jl README.
"""
function df_to_table(df; label="label", caption="caption", fmt=repeat("c", length(df)),
                     digits=4, header=replace(join(names(df), " & ")*"\\\\ \\hline",
                                              "_", "\\_")
                     )
    # row = join(["{{$x}}" for x in map(string, names(df))], " & ")

    # tpl="""
    # \\begin{table}
    #   \\centering
    #   \\begin{tabular}{$fmt}
    #   $header
    # {{#:DF}}    $row\\\\
    # {{/:DF}}  \\end{tabular}
    #   \\caption{$caption}
    #   \\label{tab:$label}
    # \\end{table}
    # """
    # render(tpl, DF=df)

    df_to_table_header(df, fmt=fmt, header=header) *
    df_to_table_body(df, digits) *
    df_to_table_footer(df, label=label, caption=caption)
end

function df_to_table_footer(df; label="label", caption="caption")
    tpl="""
      \\end{tabular}
      \\caption{$caption}
      \\label{tab:$label}
    \\end{table}
    """
    render(tpl, DF=df)
end

function df_to_table_header(df; fmt=repeat("c", length(df)), header="")
    tpl="""
    \\begin{table}
      \\centering
      \\begin{tabular}{$fmt}
        $header
    """
end

function df_to_table_body(df, digits)
    df = signif_df(df, digits)

    row = join(["{{$x}}" for x in map(string, names(df))], " & ")

    tpl="""
    {{#:DF}}    $row\\\\
    {{/:DF}}"""
    render(tpl, DF=df)
end

function signif_df(df, digits=3)
    out = DataFrame()
    colwise(x -> x isa AbstractFloat ? signif(x, digits) : x, df)

    for (n,c) in zip(names(df), DataFrames.columns(df))
        if eltype(c) <: AbstractFloat
            out[n] = signif.(c, (digits,))
        else
            out[n] = c
        end
    end
    out
end

# Compare tables, for exp0, 4, 10, 15 of U, S1, A
# exp4: longitudinal
# exp 10 50% retained
# exp 15 20% retained
exps = [0,4,10,15]
gls_PVIR = [:Unteraar, :Synthetic1, :Austfonna, :Tasman]

"""
Make a table like
```
Glacier & RMSE  & no IV & max err. & no IV  & improv. & ΔRMSE & eff_sample_size \\
```
Note that the RMSE is from the control radar lines only.
"""
function make_table0(data, sig=4)
    error("Remove Rhat")
    first = true
    last = false
    tab = 1
    for e in exps
        !first && println("    Exp $e\\\\")
        for gl in gls_PVIR
            tab = by(data[(data[:expnr].==e) .& (data[:gl].==gl),:], :gl,
                     d -> DataFrame(rmse = signif.(d[d[:rm_iv],:rmse], sig),
                                    rms_rmiv = signif.(d[.!d[:rm_iv],:rmse], sig),
                                    max = signif.(d[d[:rm_iv],:max], sig),
                                    max_rmiv = signif.(d[.!d[:rm_iv],:max], sig),
                                    iv_helps = d[d[:rm_iv],:rmse] .< d[.!d[:rm_iv],:rmse],
                                    iv_rmse_change = signif.(d[d[:rm_iv],:rmse] .- d[.!d[:rm_iv],:rmse], sig),
                                    Rhat = signif.(median(d[:Rhat]), 2)
                                    ))
            if first
                print(df_to_table_header(tab,
                                         fmt="l"*repeat("r", length(tab)-1),
                                         header="Glacier & RMSE & w/o IV & max err. & w/o IV & improv. & \$\\Delta\$ RMSE & \$\\hat{R}\$\\\\ \\hline"
                                         ))
                println("    Exp $e\\\\")
                first = false
            end
            print(df_to_table_body(tab))
        end
    end
    print(df_to_table_footer(tab,
                             label="tab:sum-$e",
                             caption="Summary of results for ITMIX2 experiments $exps for selected glaciers with IV-data.",
                             ))
end

function make_tables1(data, sig=4)
    error("remove Rhat")
    for e in exps
        tab = by(data[(data[:expnr].==e),:], :gl,
                 d -> DataFrame(rmse = signif.(d[d[:rm_iv],:rmse], sig),
                                rms_rmiv = signif.(d[.!d[:rm_iv],:rmse], sig),
                                max = signif.(d[d[:rm_iv],:max], sig),
                                max_rmiv = signif.(d[.!d[:rm_iv],:max], sig),
                                iv_helps = d[d[:rm_iv],:rmse] .< d[.!d[:rm_iv],:rmse],
                                iv_rmse_change = signif.(d[d[:rm_iv],:rmse] .- d[.!d[:rm_iv],:rmse], sig),
                                Rhat = signif.(median(d[:Rhat]), 2)
                                ))
        print(df_to_table(tab,
                          label="tab:sum-$e",
                          fmt="l"*repeat("r", length(tab)-1),
                          caption="Summary of results for ITMIX2 experiment $e for glaciers with IV-data.",
                          header="Glacier & RMSE & w/o IV & max err. & w/o IV & improv. & \$\\Delta\$ RMSE & \$\\hat{R}\$\\\\ \\hline")
              )
    end
end

# esa-14nov
#  expnr = 0
# 11×7 DataFrames.DataFrame
# │ Row │ gl             │ rmse  │ rms_rmiv │ max    │ max_rmiv │ iv_helps │ iv_rmse_change │
# ├─────┼────────────────┼───────┼──────────┼────────┼──────────┼──────────┼────────────────┤
# │ 1   │ Synthetic1     │ 39.1  │ 50.3     │ 145.0  │ 172.0    │ true     │ -11.2          │
# │ 2   │ Unteraar       │ 86.0  │ 145.0    │ 255.0  │ 393.0    │ true     │ -59.0          │
# │ 3   │ Austfonna      │ 226.0 │ 136.0    │ 1050.0 │ 686.0    │ false    │ 90.0           │
# │ 4   │ SouthGlacier   │ 61.5  │ 38.4     │ 164.0  │ 93.7     │ false    │ 23.1           │
# │ 5   │ Synthetic2     │ 131.0 │ 135.0    │ 308.0  │ 318.0    │ true     │ -4.0           │
# │ 6   │ Brewster       │ 56.3  │ 53.6     │ 99.0   │ 92.5     │ false    │ 2.7            │
# │ 7   │ Devon          │ 337.0 │ 523.0    │ 1190.0 │ 1930.0   │ true     │ -186.0         │
# │ 8   │ NorthGlacier   │ 148.0 │ 166.0    │ 375.0  │ 400.0    │ true     │ -18.0          │
# │ 9   │ Synthetic3     │ 64.5  │ 47.1     │ 231.0  │ 160.0    │ false    │ 17.4           │
# │ 10  │ Hellstugubreen │ 78.1  │ 79.1     │ 146.0  │ 146.0    │ true     │ -1.0           │
# │ 11  │ Tasman         │ 295.0 │ 223.0    │ 450.0  │ 411.0    │ false    │ 72.0           │

#  expnr = 4
# 11×7 DataFrames.DataFrame
# │ Row │ gl             │ rmse  │ rms_rmiv │ max    │ max_rmiv │ iv_helps │ iv_rmse_change │
# ├─────┼────────────────┼───────┼──────────┼────────┼──────────┼──────────┼────────────────┤
# │ 1   │ Synthetic1     │ 41.0  │ 40.2     │ 126.0  │ 122.0    │ false    │ 0.8            │
# │ 2   │ Unteraar       │ 82.0  │ 83.3     │ 246.0  │ 248.0    │ true     │ -1.3           │
# │ 3   │ Austfonna      │ 151.0 │ 131.0    │ 804.0  │ 611.0    │ false    │ 20.0           │
# │ 4   │ SouthGlacier   │ 21.4  │ 23.3     │ 66.7   │ 75.1     │ true     │ -1.9           │
# │ 5   │ Synthetic2     │ 25.3  │ 23.8     │ 89.3   │ 80.9     │ false    │ 1.5            │
# │ 6   │ Brewster       │ 26.0  │ 26.3     │ 65.7   │ 66.0     │ true     │ -0.3           │
# │ 7   │ Devon          │ 200.0 │ 259.0    │ 1000.0 │ 1670.0   │ true     │ -59.0          │
# │ 8   │ NorthGlacier   │ 35.5  │ 38.5     │ 98.2   │ 124.0    │ true     │ -3.0           │
# │ 9   │ Synthetic3     │ 21.6  │ 21.6     │ 63.0   │ 66.7     │ false    │ 0.0            │
# │ 10  │ Hellstugubreen │ 28.8  │ 29.1     │ 71.8   │ 72.0     │ true     │ -0.3           │
# │ 11  │ Tasman         │ 180.0 │ 251.0    │ 201.0  │ 398.0    │ true     │ -71.0          │

#  expnr = 10
# 11×7 DataFrames.DataFrame
# │ Row │ gl             │ rmse  │ rms_rmiv │ max   │ max_rmiv │ iv_helps │ iv_rmse_change │
# ├─────┼────────────────┼───────┼──────────┼───────┼──────────┼──────────┼────────────────┤
# │ 1   │ Synthetic1     │ 28.2  │ 27.8     │ 54.7  │ 53.8     │ false    │ 0.4            │
# │ 2   │ Unteraar       │ 62.1  │ 62.1     │ 164.0 │ 165.0    │ false    │ 0.0            │
# │ 3   │ Austfonna      │ 72.8  │ 73.2     │ 288.0 │ 291.0    │ true     │ -0.4           │
# │ 4   │ SouthGlacier   │ 20.8  │ 21.2     │ 65.4  │ 66.8     │ true     │ -0.4           │
# │ 5   │ Synthetic2     │ 25.9  │ 26.0     │ 61.2  │ 61.5     │ true     │ -0.1           │
# │ 6   │ Brewster       │ 27.5  │ 27.3     │ 50.8  │ 49.9     │ false    │ 0.2            │
# │ 7   │ Devon          │ 137.0 │ 138.0    │ 499.0 │ 487.0    │ true     │ -1.0           │
# │ 8   │ NorthGlacier   │ 35.1  │ 34.9     │ 88.1  │ 88.7     │ false    │ 0.2            │
# │ 9   │ Synthetic3     │ 22.6  │ 22.7     │ 70.3  │ 70.8     │ true     │ -0.1           │
# │ 10  │ Hellstugubreen │ 28.6  │ 28.7     │ 96.8  │ 96.7     │ true     │ -0.1           │
# │ 11  │ Tasman         │ 206.0 │ 230.0    │ 280.0 │ 255.0    │ true     │ -24.0          │

#  expnr = 15
# 11×7 DataFrames.DataFrame
# │ Row │ gl             │ rmse  │ rms_rmiv │ max   │ max_rmiv │ iv_helps │ iv_rmse_change │
# ├─────┼────────────────┼───────┼──────────┼───────┼──────────┼──────────┼────────────────┤
# │ 1   │ Synthetic1     │ 31.7  │ 31.7     │ 120.0 │ 119.0    │ false    │ 0.0            │
# │ 2   │ Unteraar       │ 52.7  │ 52.5     │ 184.0 │ 184.0    │ false    │ 0.2            │
# │ 3   │ Austfonna      │ 66.6  │ 66.8     │ 288.0 │ 291.0    │ true     │ -0.2           │
# │ 4   │ SouthGlacier   │ 19.9  │ 20.4     │ 67.6  │ 66.7     │ true     │ -0.5           │
# │ 5   │ Synthetic2     │ 23.5  │ 23.3     │ 52.7  │ 52.5     │ false    │ 0.2            │
# │ 6   │ Brewster       │ 26.0  │ 23.7     │ 54.7  │ 50.0     │ false    │ 2.3            │
# │ 7   │ Devon          │ 152.0 │ 160.0    │ 479.0 │ 451.0    │ true     │ -8.0           │
# │ 8   │ NorthGlacier   │ 37.3  │ 37.9     │ 99.7  │ 103.0    │ true     │ -0.6           │
# │ 9   │ Synthetic3     │ 22.7  │ 22.9     │ 79.2  │ 78.1     │ true     │ -0.2           │
# │ 10  │ Hellstugubreen │ 40.2  │ 41.0     │ 73.8  │ 74.6     │ true     │ -0.8           │
# │ 11  │ Tasman         │ 224.0 │ 220.0    │ 350.0 │ 313.0    │ false    │ 4.0            │

using Plots

function plot_PVIR_2x8(name, expnr, subdir=sd;
                       reuse=true,
                       xyscale=1e3,
                       kwargs...)
    gid = BM.ITMIXGlacier(name,2)
    gl, gb, pp, pm, pn, pl = gls[gid];
    ps = Array{Any}(2,4)

    # with IV
    for i=1:2
        xf = i==2 ? x->round(Int, x) : x->""
        xl = i==2 ? "χ (km)" : ""

        # https://discourse.julialang.org/t/avoid-repeating-tick-labels-with-subplots-in-plots-jl/351/4
        data_out, hs2d, hs2ds, iv2d, iv2ds, r_lines = load_one(gid, expnr, [true,false][i]);

        ps[i,1] = BM.plot2d_h(gl, hs2d.v, xyscale=xyscale, colorbar=false,
                              xformatter=xf, xlabel=xl)
        ps[i,2] = BM.plot2d_h(gl, hs2ds.v, xyscale=xyscale, colorbar=false,
                              xformatter=xf, yformatter=x->"", xlabel=xl, ylabel="")
        ps[i,3] = heatmap(gl.dem.x/xyscale, gl.dem.y/xyscale, iv2d.v', size=(600,600), #, zlims=(0,50))
                          aspect_ratio=:equal, colorbar=false,
                          title= i==1 ? "iv" :"",
                          xformatter=xf, xlabel=xl,
                          kwargs...);
        ps[i,4] = heatmap(gl.dem.x/xyscale, gl.dem.y/xyscale, iv2ds.v', size=(600,600), #, zlims=(0,50))
                          aspect_ratio=:equal, colorbar=false,
                          title= i==1 ? "std(iv)" :"",
                          xformatter=xf, yformatter=x->"", xlabel=xl, ylabel="",
                          kwargs...);
    end
    c1 = scatter([0,0], [0,1], zcolor=[clims...], clims=clims,# aspect_ratio=:equal,
                 xlims=(1,1.1), axis=false, label="", colorbar_title=tit,
                 seriescolor=seriescolor
                 )#, markeralpha=0, markerstrokewidth=0)

    l = @layout [grid(2,2) a{0.01w} grid(2,2) a{0.01w}]
    plot(permutedims(ps, [2,1])[1:4]..., c1, permutedims(ps, [2,1])[5:8]..., c2, layout=l, reuse=reuse; size=(900,600), kwargs...)
end

# ## hack to place a colorbar willy-nilly
# ## https://github.com/tbreloff/Plots.jl/issues/412
# h1 = contour(rand(3,3), cbar=false, clims=(0,3), fill=true)
# scatter!(1:3,1:3, zcolor=1:3, cbar=false, clims=(0,3))
# h3 = contour(rand(3,3), cbar=false, clims=(0,3), fill=true)
# h2 = scatter([0,0], [0,1], zcolor=[0,3], clims=(0,3),
#              xlims=(1,1.1), axis=false, grid=false, label="", colorbar_title="cbar")
# l = @layout [grid(2,1) a{0.01w}]
# plot(h1,h3,h2,layout=l)

#     h2 = scatter([0,0], [0,1], zcolor=[clims...], clims=clims,# aspect_ratio=:equal,
#                  xlims=(1,1.1), axis=false, label="", colorbar_title=tit,
#                  seriescolor=seriescolor
#                  )#, markeralpha=0, markerstrokewidth=0)
#     l = @layout [grid(1,1) a{0.01w}]
#     plot(h1,h2,layout=l, size=fsize, palette=:viridis; kwargs...)

function plot_PVIR_2x2(name, expnr, subdir=sd;
                       plotiv=false,
                       reuse=true,
                       xyscale=1e3,
                       seriescolor=[:inferno, :pu_or][1],
                       fsize=(1800,900),
                       kwargs...)
    gid = BM.ITMIXGlacier(name,2)
    gl, gb, pp, pm, pn, pl = gls[gid];
    ps = Array{Any}(2,3)

    clims_h = (1, 1.0)
    clims_sh = (0, 1e-5)
    clims_err = (-1e-5, 1e-5)
    # just to the plot limits:
    for rm_iv = [true,false]
        data_out, hs2d, hs2ds, iv2d, iv2ds, r_lines_fitting = load_one(gid, expnr, rm_iv);
        d, ds = plotiv ? (iv2d, iv2ds) : (hs2d, hs2ds)
        rd = plotiv ? gl.iv.vals : gl.h.vals
        err, err_ = if plotiv
            err = BM.error2d(d.v, gl, :iv, gl.ivmask)
            (err, VAWTools.extremanan(err))
        else
            err = BM.error_h(d.v, gl)
            (err, VAWTools.extremanan(err.v))
        end
        h = VAWTools.extremanan(d.v)
        sh = VAWTools.extremanan(ds.v)
        clims_h = (min(clims_h[1], h[1]), max(clims_h[2], h[2]))
        clims_sh = (min(clims_sh[1], sh[1]), max(clims_sh[2], sh[2]))
        clims_err = (min(clims_err[1], err_[1]), max(clims_err[2], err_[2]))
    end

    labels = ["Exp. $expnr with IV", "Exp. $expnr w/o IV"]
    letters = [["A", "B", "C"], ["D", "E", "F"]]
    # plot with IV and without:
    for i=1:2 # rows, corresponding to rm_iv false/true
        xf = i==2 ? x->round(Int, x) : x->""
        xl = i==2 ? "χ (km)" : ""
        # https://discourse.julialang.org/t/avoid-repeating-tick-labels-with-subplots-in-plots-jl/351/4
        data_out, hs2d, hs2ds, iv2d, iv2ds, r_lines_fitting = load_one(gid, expnr, [true,false][i]);
        d, ds = plotiv ? (iv2d, iv2ds) : (hs2d, hs2ds)
        rd = plotiv ? gl.iv.vals : gl.h.vals
        err, rmse, linf = if plotiv
            err = BM.error2d(d.v, gl, :iv, gl.ivmask)
            vv = err[.!isnan.(err)]
            rmse = round(Int, sqrt(sum(vv.^2)/length(vv)) )
            (err, rmse, round(Int, maximum(abs.(vv))) )
        else
            err = BM.error_h(d.v, gl)
            vv = err.v[.!isnan.(err.v)]
            # RMSE, max-error on control lines
            rmse, maxerr = BM.compare_to_radar(gl, d, ds, r_lines_fitting, thin= (name==:Tasman ? 1 : 10)) # TODO hack...
            # maxerr = (maximum(abs.(vv)),)

            (err, round(Int, rmse[1]), round(Int, maxerr[1]))
        end

        ## Plot value
        ps[i,1] = BM.plot2d_h(gl, d.v, xyscale=xyscale, colorbar=false,
                              plotiv=plotiv,
                              xformatter=xf, xlabel=xl, clims=clims_h,
                              seriescolor=seriescolor,
                              plotline=plotiv ? false : I2.get_rlines(gid, expnr),
                              label=labels[i],
                              outlines=true
                              )
        # add a legend
        scatter!([d.x[1]/xyscale], [d.y[1]/xyscale], label="_", legendtitle=labels[i], markeralpha=0, markerstrokewidth=0,
                 xformatter=xf, xlabel=xl)

        ## Plot std
        ps[i,2] = BM.plot2d_h(gl, ds.v, xyscale=xyscale, colorbar=false,
                              plotiv=plotiv,
                              xformatter=xf, yformatter=x->"", xlabel=xl, ylabel="",
                              clims=clims_sh,
                              plotline=plotiv ? I2.get_rlines(gid, expnr) : false,
                              seriescolor=seriescolor,
                              #plotline=true,
                              plottraj=false,
                              outlines=true)

        ## Plot error
        lt = if rd isa Traj
            # plot error on the radar tracks
            ps[i,3] = Plots.plot();
            BM.plot_traj!(err; xyscale=xyscale, plotline=I2.get_rlines(gid, expnr),
                          colorbar=false, aspect_ratio=:equal,
                          yformatter=x->"", ylabel="",
                          xformatter=xf, xlabel=xl,
                          clims=clims_err)
           "RMSE = $(rmse)m\nMax err = $(linf)m"
        else
            # this does not work for point measurements, such as SouthGlacier
            err[.!gl.ivmask] = NaN
            ps[i,3] = heatmap(gl.dem.x/xyscale, gl.dem.y/xyscale, err',
                              colorbar=false, aspect_ratio=:equal,
                              clims=clims_err,
                              kwargs...);
            "RMSE = $(rmse)m/a\nMax err = $(linf)m/a"
        end
        # make a label
        scatter!([d.x[1]/xyscale], [d.y[1]/xyscale], label="_", legendtitle=lt,
                 markeralpha=0, markerstrokewidth=0,
                 xformatter=xf, xlabel=xl, yformatter=x->"")
        BM.plot_outlines!(gl, xyscale=xyscale, ylabel="", xlabel=xl)
    end
    xls,yls = xlims(), ylims()
    for p in ps
        xlims!(p, xls)
        ylims!(p, yls)
    end
    # ## hack to place a colorbar willy-nilly
    # ## https://github.com/tbreloff/Plots.jl/issues/412

    c1 = scatter([0,0], [0,1], zcolor=[clims_h...], clims=clims_h,# aspect_ratio=:equal,
                 xlims=(1,1.1), axis=false, label="", colorbar_title=plotiv ? "iv (m/a)" : "h (m)",
                 seriescolor=seriescolor, grid=false
                 )#, markeralpha=0, markerstrokewidth=0)
    c2 = scatter([0,0], [0,1], zcolor=[clims_sh...], clims=clims_sh,# aspect_ratio=:equal,
                 xlims=(1,1.1), axis=false, label="", colorbar_title=plotiv ? "std(iv) (m/a)" : "std(h) (m)",
                 seriescolor=seriescolor, grid=false
                 )#, markeralpha=0, markerstrokewidth=0)
    c3 = scatter([0,0], [0,1], zcolor=[clims_err...], clims=clims_err,# aspect_ratio=:equal,
                 xlims=(1,1.1), axis=false, label="", colorbar_title=plotiv ? "error(iv) (m/a)" : "error(h) (m)",
                 seriescolor=seriescolor, grid=false
                 )#, markeralpha=0, markerstrokewidth=0)

    l = @layout [grid(2,1) a{0.01w} grid(2,1) b{0.01w} grid(2,1) b{0.01w}]
    # l = @layout [grid(2,1) [a{0.01w}
    #                         b{0.01w}
    #                         c{0.01w}] grid(2,1) d{0.01w}]
    #l = @layout [grid(2,1) grid(3,1,widths=[0.01,0.01,0.01]) grid(2,1) d{0.01w}]
    plot(ps[1:2]..., c1, ps[3:4]..., c2, ps[5:6]..., c3, layout=l, reuse=reuse; size=fsize, kwargs...)
end

# dir_ = "/home/mauro/projects/ESA-ice-volume/BITEModel/scripts/data/ITMIX2/results_Werder/esa-14nov/plots/esa/"
# error()
# plotiv = false
# for name in gls_PVIR
#     @show name
#     for e in exps
#         ee = e<10 ? "0$e" : "$e"
#         @show ee
#         fln = plotiv ? "iv" : "h"
#         plot_PVIR_2x2(name, e, plotiv=plotiv)
#         savefig(dir_*"$(fln)_$name-$ee.png")
# #        error()
#     end
# end

############
# get Matthias' results
function get_matthias()
    data_huss = make_data()

    # ITMIX2
    dir = "/home/mauro/projects/ESA-ice-volume/reporting/n6_Product_Validation_Intercomparison/itmix2-results-matthias"
    for name in map(BM.getid, keys(gls))
        gid = BM.ITMIXGlacier(name,2)
        for e in 1:16
            push!(data_huss, read_others_results(gid, e, dir, "Huss"))
        end
    end

    # ITMIX1
    dir = "/home/mauro/projects/ESA-ice-volume/reporting/n6_Product_Validation_Intercomparison/itmix1-results-matthias"
    for name in map(BM.getid, keys(gls))
        name in [:Austfonna,:Devon] && continue
        gid = BM.ITMIXGlacier(name,2)
        for e in 0
            push!(data_huss, read_others_results(gid, e, dir, "Huss"))
            # if name==:Austfonna
            #     push!(data_huss, read_others_results(gid, e, dir, "Huss", thin=6, start=(190, 523), stop=(2927, ))
            # else
        #     push!(data_huss, read_others_results(gid, e, dir, "Huss"))
            # end
        end
    end
    sort!(data_huss)
end

data_huss = get_matthias()
## Compare
glids0 = [BM.getid(k) for k in keys(gls) if !(BM.getid(k) in [:Austfonna, :Devon])]
glids = [BM.getid(k) for k in keys(gls)]

data_huss0 = data_huss[in.(data_huss[:gl], (glids0,)), :];
huss = summarize_data(data_huss0, expnrs=0:16, person=true);
huss0 = summarize_data(data_huss0, expnrs=0, person=true);
huss1_16 = summarize_data(data_huss, expnrs=1:16, person=true);
huss4 = summarize_data(data_huss, expnrs=4, person=true);

data0 = data[in.(data[:gl], (glids0,)), :];
me = summarize_data(data0[data0[:rm_iv].==false,:], expnrs=0:16, person=true);
me0 = summarize_data(data0[data0[:rm_iv].==false,:], expnrs=0, person=true);
me1_16 = summarize_data(data[data[:rm_iv].==false,:], expnrs=1:16, person=true);
me4 = summarize_data(data[data[:rm_iv].==false,:], expnrs=4, person=true);

me_rmiv = summarize_data(data0[data0[:rm_iv].==true,:], expnrs=0:16, person=true);
me0_rmiv = summarize_data(data0[data0[:rm_iv].==true,:], expnrs=0, person=true);
me1_16_rmiv = summarize_data(data[data[:rm_iv].==true,:], expnrs=1:16, person=true);

function compare_two(summarized_huss, summarized_me,
                     median_mean_h = (sort(summarized_huss)[:median_mean_h] .+ sort(summarized_me)[:median_mean_h])/2)
    huss, me = sort(summarized_huss), sort(summarized_me)
    fls = [(:rel_mean_rmse, :mean_rmse)]
    out = DataFrame()
    out[:gl] = huss[:gl]
    out[:median_mean_h_huss] = huss[:median_mean_h]
    out[:median_mean_h_me] = me[:median_mean_h]
    out[:median_mean_h] = median_mean_h

    out[:rel_mean_rmse_huss] = huss[:mean_rmse]./median_mean_h
    out[:rel_mean_rmse_me] = me[:mean_rmse]./median_mean_h

    for (flo,fl) in fls
        rmse_diff = (huss[fl] .- me[fl])./median_mean_h
        out[flo] = rmse_diff
        out[:better] = significatly_better.(rmse_diff, 0.05)
    end
    out
end

function significatly_better(diff, thresh)
    if diff>thresh
        return "true"
    elseif diff<-thresh
        return "false"
    else
        return "-"
    end
end

# tables of Intercomparison
df_to_table(compare_two(huss1_16, me1_16)[[:gl, :median_mean_h_huss, :median_mean_h_me, :rel_mean_rmse_huss, :rel_mean_rmse_me, :better]]
                                           , label="comp1-16", caption="Comparing ITMIX2 results between Huss and Werder",
            fmt="lrrrrr", digits=3,
            header="Glacier & \$\\bar{h}_{Huss}\$ & \$\\bar{h}_{BITE}\$ & \$\\overline{\\rm RMSE}_{H\\&F}\$ & \$\\overline{\\rm RMSE}_{BITE}\$ & better\\\\ \\hline") |> print

df_to_table(compare_two(huss4, me4)[[:gl, :median_mean_h_huss, :median_mean_h_me, :rel_mean_rmse_huss, :rel_mean_rmse_me, :better]]
                                           , label="comp4", caption="Comparing ITMIX2 results between Huss and Werder for Exp04",
            fmt="lrrrrr", digits=3,
            header="Glacier & \$\\bar{h}_{Huss}\$ & \$\\bar{h}_{BITE}\$ & \$\\overline{\\rm RMSE}_{H\\&F}\$ & \$\\overline{\\rm RMSE}_{BITE}\$ & better\\\\ \\hline") |> print

df_to_table(compare_two(huss0, me0,
                        (sort(huss0)[:median_mean_h] .+ sort(me0)[:median_mean_h])/2
                        )[[:gl, :median_mean_h_huss, :median_mean_h_me, :rel_mean_rmse_huss, :rel_mean_rmse_me, :better]]
                                           , label="comp4", caption="Comparing ITMIX results between Huss and Werder (i.e. Exp00)",
            fmt="lrrrrr", digits=3,
            header="Glacier & \$\\bar{h}_{Huss}\$ & \$\\bar{h}_{BITE}\$ & \$\\overline{\\rm RMSE}_{H\\&F}\$ & \$\\overline{\\rm RMSE}_{BITE}\$ & better\\\\ \\hline") |> print

df_to_table(compare_two(huss0, me0_rmiv,
                        (sort(huss0)[:median_mean_h] .+ sort(me0)[:median_mean_h])/2
                        )[[:gl, :median_mean_h_huss, :median_mean_h_me, :rel_mean_rmse_huss, :rel_mean_rmse_me, :better]]
                                           , label="comp4", caption="Comparing ITMIX results between Huss and Werder without using IV (i.e. Exp00)",
            fmt="lrrrrr", digits=3,
            header="Glacier & \$\\bar{h}_{Huss}\$ & \$\\bar{h}_{BITE}\$ & \$\\overline{\\rm RMSE}_{H\\&F}\$ & \$\\overline{\\rm RMSE}_{BITE}\$ & better\\\\ \\hline") |> print

function get_mean_thicknesses(data)
    d = summarize_data(data)[[:gl, :median_mean_h]]
    Dict(k=>v for (k,v) in zip(d[:gl], d[:median_mean_h]))
end

function plot_whiskers(data_huss, data_me, data_mean_h=vcat(data_huss, data_me); kw...)
    # join
    all = vcat(data_huss, data_me);
    all[:gl] = CategoricalArray(string.(all[:gl]));

    # get median mean thicknesses
    mts = get_mean_thicknesses(data_mean_h)

    # explode (manual stack)
    # need [person, gl, expnr, err_h]
    #    T = CategoricalArrays.CategoricalArray{String,1,UInt32,String,CategoricalArrays.CategoricalString{UInt32},Union{}}
    T = CategoricalArrays.CategoricalString{UInt32}
    newer = DataFrame([T, T, Int64, Float64], [:person, :gl, :rm_iv, :err_h], 0)

    for g in groupby(all, [:person, :gl, :rm_iv])
        person, gl, rm_iv = g[:person][1], g[:gl][1], g[:rm_iv][1]
        person = person==:Huss ? "H&F" : "BITE"
        gl = gl=="Hellstugubreen" ? "Hellstugub." : gl
        mt = mts[Symbol(g[:gl][1])]
        err_h = filter(x->!isnan(x), vcat([gg.v for gg in g[:err_h]]...))
        newer = vcat(newer, DataFrame(person=person, gl=gl, rm_iv=rm_iv, err_h=err_h))
    end

    StatPlots.@df newer boxplot(:gl, :err_h, group=:person; kw...)
end

dir_ = "/home/mauro/projects/ESA-ice-volume/BITEModel/scripts/data/ITMIX2/results_Werder/esa-14nov/plots/esa/"

plot_whiskers(data_huss0[data_huss0[:expnr].==0,:],
              data0[(data0[:expnr].==0) .&  (data0[:rm_iv].==false),:], vcat(data_huss, data);
              ylims=(-250,250), ylabel="rel. error h (m)",
              size=(900,600), yticks=-500:50:500);
savefig(dir_*"intercomp-huss-exp00.png")

plot_whiskers(data_huss0[data_huss0[:expnr].==0,:],
              data0[(data0[:expnr].==0) .&  (data0[:rm_iv].==true),:], vcat(data_huss, data);
              ylims=(-250,250), ylabel="rel. error h (m)",
              size=(900,600), yticks=-500:50:500);
savefig(dir_*"intercomp-huss-exp00-rmiv.png")

plot_whiskers(data_huss[data_huss[:expnr].==4,:],
              data[(data[:expnr].==4) .&  (data[:rm_iv].==false),:], vcat(data_huss, data);
              ylims=(-250,250), ylabel="rel. error h (m)",
              size=(900,600), yticks=-500:50:500);
savefig(dir_*"intercomp-huss-exp04.png")

dh = data_huss[ in.(data_huss[:expnr], (1:16,)), :];
dm = data[ in.(data[:expnr], (1:16,))  .&  (data[:rm_iv].==false), :];
plot_whiskers(dh, dm,
              vcat(data_huss, data);
              ylims=(-250,250), ylabel="rel. error h (m)",
              size=(900,600), yticks=-500:50:500);
savefig(dir_*"intercomp-huss-exp01-16.png")


text:
export rule

@rule GCV(:x, Marginalisation) (m_y::UniNormalOrExpLinQuad, q_z::Any, q_κ::Any, q_ω::Any) = begin

    y_mean, y_var = mean_var(m_y)
    z_mean, z_var = mean_var(q_z)
    κ_mean, κ_var = mean_var(q_κ)
    ω_mean, ω_var = mean_var(q_ω)

    ksi = κ_mean ^ 2 * z_var + z_mean ^ 2 * κ_var + z_var * κ_var
    A = exp(-ω_mean + ω_var / 2)
    B = exp(-κ_mean * z_mean + ksi / 2)

    return NormalMeanVariance(y_mean, y_var + inv(A * B))
end

@rule GCV(:x, Marginalisation) (q_y::Any, q_z::Any, q_κ::Any, q_ω::Any) = begin

    y_mean        = mean(q_y)
    z_mean, z_var = mean_var(q_z)
    κ_mean, κ_var = mean_var(q_κ)
    ω_mean, ω_var = mean_var(q_ω)

    ksi = κ_mean ^ 2 * z_var + z_mean ^ 2 * κ_var + z_var * κ_var
    A = exp(-ω_mean + ω_var / 2)
    B = exp(-κ_mean * z_mean + ksi / 2)

    return NormalMeanVariance(y_mean, inv(A * B))
end

ext: jl


text:
module GGThemr

using AbstractPlotting, Colors

export ggthemr, show_ggthemr, ggthemr_colorthemes

function show_ggthemr(theme::Symbol)
    AbstractPlotting.set_theme!(ggthemr(theme))
    scene = scatter(randn(20), randn(20), markersize = 0.2)
    scatter!(scene, randn(20), randn(20), markersize = 0.2)
    scatter!(scene, randn(20), randn(20), markersize = 0.2)
    scatter!(scene, randn(20), randn(20), markersize = 0.2)
    scatter!(scene, randn(20), randn(20), markersize = 0.2)
end

ggthemr_colorthemes() = collect(keys(ColorTheme))

const ggthemr_style = Theme(
    linewidth = 2,
    font = "NotoSans",
    axis2d = Theme(
        frame = Theme(
            linewidth = 1.5,
            frames = ((true, false), (true, false))
        ),
        grid = Theme(
            linewidth = (2, 2),
            linestyle = (:dash, :dash),
        ),
        ticks = Theme(
            textsize = (4.5, 4.5),
        ),
    ),
    scatter = Theme(
        markersize = 6,
        strokewidth = 2,
    )
)

function ggthemr(theme::Symbol)
    ct = ColorTheme[theme]
    palettes = merge(AbstractPlotting.default_palettes, Attributes(color = ct[:swatch][2:end]))
    merge(ggthemr_style, Theme(
        palette = palettes,
        color = ct[:swatch][2], # maybe it should be ct[:swatch][1], this is the color to use in absence of grouping
        backgroundcolor = parse(Color, ct[:background]),
        #colorgradient = parse(Color, ct[:gradient]),
        axis = Theme(
            frame = Theme(
                linecolor = parse(Color, ct[:line][2])
            ),
            grid = Theme(
                linecolor = (parse(Color, ct[:gridline]), parse(Color, ct[:gridline]))
            ),
            ticks = Theme(
                textcolor = (parse(Color, ct[:text][2]), parse(Color, ct[:text][2]))
            ),
            names = Theme(
                textcolor = (parse(Color, ct[:text][2]), parse(Color, ct[:text][2]))
            ),
        ),
        scatter = Theme(
            strokecolor = parse(Color, ct[:background]),
        ),
    ))
end

const ColorTheme = Dict{Symbol, Dict{Symbol, Any}}()

include("palettes.jl")

end #module


text:
# Computes the s_2 function.
include("ellpic_bulirsch.jl")

function s2(r::T,b::T,Kofk::T,Eofk::T) where {T <: Real}
# For now, just compute linear component:
Lambda1 = zero(typeof(b))
if b > 1.0+r ||  r == 0.0
  # No occultation:
  Lambda1 = 0  # Case 1
elseif b < r-1.0
  # Full occultation:
  Lambda1 = 0  # Case 11
else
  if b == 0
    Lambda1 = -2/3*(1.0-r^2)^(3/2) # Case 10
  elseif b==r
    if r == 0.5
      Lambda1 = 1/3-4/(9pi) # Case 6
    elseif r < 0.5
#      Kofk = ellk_bulirsch(4r^2); Eofk = ellec_bulirsch(4r^2)
#      Lambda1 = 1/3+2/(9pi)*(4*(2r^2-1)*ellec_bulirsch(4r^2)+(1-4r^2)*ellk_bulirsch(4r^2)) # Case 5
      Lambda1 = 1/3+2/(9pi)*(4*(2r^2-1)*Eofk+(1-4r^2)*Kofk) # Case 5
    else
#      Kofk = ellk_bulirsch(inv(4r^2)); Eofk = ellec_bulirsch(inv(4r^2))
#      Lambda1 = 1/3+16r/(9pi)*(2r^2-1)*ellec_bulirsch(inv(4r^2))-(1-4r^2)*(3-8r^2)/(9pi*r)*ellk_bulirsch(inv(4r^2)) # Case 7
      Lambda1 = 1/3+16r/(9pi)*(2r^2-1)*Eofk-(1-4r^2)*(3-8r^2)/(9pi*r)*Kofk # Case 7
    end
  else
    k2 = 1.0+(1.0-(b+r)^2)/(4*b*r)
    if (b+r) > 1.0 # k^2 < 1, Case 2, Case 8
#      Kofk = ellk_bulirsch(k2)
#      Eofk = ellec_bulirsch(k2)
      Piofnk = ellpic_bulirsch(-k2*(b+r)^2,k2)
      xi = 2*b*r*(4-7r^2-b^2)
      Lambda1 = (((r+b)^2-1)/(r+b)*(-2r*(2*(r+b)^2+(r+b)*(r-b)-3)*Kofk+3*(b-r)*Piofnk)
            -2*xi*Eofk)/(9*pi*sqrt(b*r))
    elseif (b+r) < 1.0  # k^2 > 1, Case 3, Case 9
      k2inv = inv(k2)
#      Kofk = ellk_bulirsch(k2inv)
#      Eofk = ellec_bulirsch(k2inv)
      if abs(b-r) != 1.0
        Piofnk = ellpic_bulirsch(-inv(k2*(b+r)^2),k2inv)/sqrt(1-(b-r)^2)
      else
        Piofnk = 0.0
      end
      Lambda1 = 2*((1-(r+b)^2)*(sqrt(1-(b-r)^2)*Kofk+3*(b-r)/(b+r)*Piofnk)
            -sqrt(1-(b-r)^2)*(4-7r^2-b^2)*Eofk)/(9*pi)
    else
      # b+r = 1 or k^2=1, Case 4 (extending r up to 1)
      Lambda1 = 2/(3pi)*acos(1.-2*r)-4/(9pi)*(3+2r-8r^2)*sqrt(r*b)-2/3*convert(typeof(b),r>.5)
    end
  end
end
s2 = 1.0-1.5*Lambda1-convert(typeof(b),r>b)
return s2*2pi/3
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
"""
    couple_models(aero::Peters, stru::TypicalSection, flap::SimpleFlap)

Create an aerostructural model using the unsteady aerodynamic model defined by
Peters et al., a two-degree of freedom typical section model, and a linear
steady-state control surface model.  This model introduces the freestream
velocity ``U_\\infty``, air density ``\\rho_\\infty``, and flap deflection
``\\delta`` as additional parameters.
"""
couple_models(aero::Peters, stru::TypicalSection, flap::SimpleFlap) = (aero, stru, flap)

# --- Traits --- #

number_of_additional_parameters(::Type{<:Peters}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = 3

coupling_inplaceness(::Type{<:Peters}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = OutOfPlace()

coupling_rate_jacobian_type(::Type{<:Peters}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Linear()
coupling_state_jacobian_type(::Type{<:Peters}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Nonlinear()
coupling_parameter_jacobian_type(::Type{<:Peters}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Nonlinear()
coupling_time_gradient_type(::Type{<:Peters}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Zeros()

# --- Methods --- #

function get_coupling_inputs(aero::Peters{N,TF,SV,SA}, stru::TypicalSection,
    flap::SimpleFlap, dx, x, p, t) where {N,TF,SV,SA}

    # extract rate variables
    dλ = dx[SVector{N}(1:N)]
    dh, dθ, dhdot, dθdot = dx[SVector{4}(N+1:N+4)]

    # extract state variables
    λ = x[SVector{N}(1:N)]
    h, θ, hdot, θdot = x[SVector{4}(N+1:N+4)]

    # extract parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p

    # extract model constants
    bbar = aero.b

    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    udot, vdot, ωdot = section_accelerations(dhdot, dθdot)

    # aerodynamic loads
    Na, Aa, Ma = peters_loads(a, b, ρ, a0, α0, cd0, cm0, bbar, u, v, ω, vdot, ωdot, λ)

    # loads due to flap deflections
    Nf, Af, Mf = simpleflap_loads(b, u, ρ, cnδ, caδ, cmδ, δ)

    # loads per unit span
    L = Na + Nf
    M = Ma + Mf

    # return inputs
    return SVector(u, ω, vdot, ωdot, L, M)
end

# --- Performance Overloads --- #

# TODO: State rate, state, and parameter jacobians

# --- Convenience Methods --- #

function set_additional_parameters!(padd, aero::Peters, stru::TypicalSection,
    flap::SimpleFlap; U, rho, delta)

    padd[1] = U
    padd[2] = rho
    padd[3] = delta

    return padd
end

function separate_additional_parameters(aero::Peters, stru::TypicalSection,
    flap::SimpleFlap, padd)

    return (U = padd[1], rho = padd[2], delta = padd[3])
end

# --- Plotting --- #

@recipe function f(aero::Peters{N,TF,SV,SA}, stru::TypicalSection, flap::SimpleFlap,
    dx, x, y, p, t) where {N,TF,SV,SA}

    framestyle --> :origin
    grid --> false
    xlims --> (-1.0, 1.0)
    ylims --> (-1.5, 1.5)
    label --> @sprintf("t = %6.3f", t)

    # extract state variables
    λ = x[SVector{N}(1:N)]
    h, θ, hdot, θdot = x[SVector{4}(N+1:N+4)]

    # extract parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p

    xplot = [-(0.5 + a*b)*cos(θ),    (0.5 - a*b)*cos(θ)]
    yplot = [ (0.5 + a*b)*sin(θ)-h, -(0.5 - a*b)*sin(θ)-h]

    return xplot, yplot
end


text:
"""
   `ecc_anom_init_guess_danby(mean_anomaly, eccentricity)`

Returns initial guess for the eccentric anomaly for use by itterative solvers of Kepler's equation for bound orbits.

Based on "The Solution of Kepler's Equations - Part Three"
Danby, J. M. A. (1987) Journal: Celestial Mechanics, Volume 40, Issue 3-4, pp. 303-312 (1987CeMec..40..303D)
"""
function ecc_anom_init_guess_danby(M::Real, ecc::Real)
	@assert -2π<= M <= 2π
	@assert 0 <= ecc <= 1.0
    if  M < zero(M)
		M += 2π
	end
    E = (M<π) ? M + 0.85*ecc : M - 0.85*ecc
end;

"""
   `update_ecc_anom_laguerre(eccentric_anomaly_guess, mean_anomaly, eccentricity)`

Update the current guess for solution to Kepler's equation

Based on "An Improved Algorithm due to Laguerre for the Solution of Kepler's Equation"
   Conway, B. A.  (1986) Celestial Mechanics, Volume 39, Issue 2, pp.199-211 (1986CeMec..39..199C)
"""
function update_ecc_anom_laguerre(E::Real, M::Real, ecc::Real)
  (es, ec) = ecc .* sincos(E)
  F = (E-es)-M
  Fp = one(M)-ec
  Fpp = es
  n = 5
  root = sqrt(abs((n-1)*((n-1)*Fp*Fp-n*F*Fpp)))
  denom = Fp>zero(E) ? Fp+root : Fp-root
  return E-n*F/denom
end;

"""
   `calc_ecc_anom( mean_anomaly, eccentricity )`
   `calc_ecc_anom( param::Vector )``
	Estimates eccentric anomaly for given 'mean_anomaly' and 'eccentricity'.
If passed a parameter vector, param[1] = mean_anomaly and param[2] = eccentricity.
Optional parameter `tol` specifies tolerance (default 1e-8)
"""
function calc_ecc_anom end
function calc_ecc_anom(mean_anom::Real, ecc::Real; tol::Real = 1.0e-8)
	@assert 0 <= ecc <= 1.0
	@assert 1e-16 <= tol < 1
  	M = rem2pi(mean_anom,RoundNearest)
    E = ecc_anom_init_guess_danby(M,ecc)
	local E_old
    max_its_laguerre = 200
    for i in 1:max_its_laguerre
       E_old = E
       E = update_ecc_anom_laguerre(E_old, M, ecc)
       if abs(E-E_old) < tol break end
    end
    return E
end

function calc_ecc_anom(param::Vector; tol::Real = 1.0e-8)
	@assert length(param) == 2
	calc_ecc_anom(param[1], param[2], tol=tol)
end


text:
include("input.jl")
include("dist_input.jl")

nusers{T<:Inputs}(inp::T) = inp.nusers
nitems{T<:Inputs}(inp::T) = inp.nitems

users_and_ratings{T<:Inputs}(inp::T, i::Int64) = _sprowsvals(get(inp.R), i)
all_user_ratings{T<:Inputs}(inp::T, i::Int64) = _spvals(get(inp.R), i)
all_users_rated{T<:Inputs}(inp::T, i::Int64) = _sprows(get(inp.R), i)

items_and_ratings{T<:Inputs}(inp::T, u::Int64) = _sprowsvals(get(inp.RT), u)
all_item_ratings{T<:Inputs}(inp::T, u::Int64) = _spvals(get(inp.RT), u)
all_items_rated{T<:Inputs}(inp::T, u::Int64) = _sprows(get(inp.RT), u)

function _sprowsvals{R<:InputRatings}(sp::R, col::Int64)
    s = SparseVector(sp[:,col])
    nonzeroinds(s), nonzeros(s)
end

function _sprows{R<:InputRatings}(sp::R, col::Int64)
    s = SparseVector(sp[:,col])
    nonzeroinds(s)
end

function _spvals{R<:InputRatings}(sp::R, col::Int64)
    s = SparseVector(sp[:,col])
    nonzeros(s)
end

function _sprowsvals(sp::ParallelSparseMatMul.SharedSparseMatrixCSC{Float64,Int64}, col::Int64)
    rowstart = sp.colptr[col]
    rowend = sp.colptr[col+1] - 1
    sp.rowval[rowstart:rowend], sp.nzval[rowstart:rowend]
end

function _sprows(sp::ParallelSparseMatMul.SharedSparseMatrixCSC{Float64,Int64}, col::Int64)
    rowstart = sp.colptr[col]
    rowend = sp.colptr[col+1] - 1
    sp.rowval[rowstart:rowend]
end

function _spvals(sp::ParallelSparseMatMul.SharedSparseMatrixCSC{Float64,Int64}, col::Int64)
    rowstart = sp.colptr[col]
    rowend = sp.colptr[col+1] - 1
    sp.nzval[rowstart:rowend]
end


text:
using FractionalTimeDG
import FractionalTimeDG
using SpecialFunctions: erfcx
import GaussQuadrature
using OffsetArrays

T = Float64
two = parse(T, "2")
α = 1 / two
λ = 1 / two
tmax = two

Ehalf(x) = erfcx(-x)

f(t) = cos(π*t)
u0 = one(T)

function u(t::T, λ::T, u0::T, f::Function, M::Integer,
           store::Store{T}) where T <: AbstractFloat
    y, wy = GaussQuadrature.legendre(T, M)
    y .= ( y .+ 1 ) / 2
    wy .= wy / 2
    s = zero(T)
    for m = 1:M
        s += wy[m] * Ehalf(-λ*sqrt(t)*y[m]) * f(t*(1-y[m]^2)) * y[m]
    end
    return u0 * Ehalf(-λ*sqrt(t)) + 2t * s
end


text:
T=6
tsname="joLinearFunction - Difference"
@testset "$tsname" begin
for t=1:T # start test loop

if t<=2
    (m,n)=(5,5)
else
    (m,n)=(rand(3:7),rand(3:7))
end

if t%2==1
    tname="loop $t for real($m,$n)"
    a=rand(m,n)
    b=rand(m,n)
    vn=rand(n)
    mvn=rand(n,2)
    vm=rand(m)
    mvm=rand(m,2)
    A=joLinearFunctionAll(m,n,v->a*v,v->transpose(a)*v,v->adjoint(a)*v,v->conj(a)*v,v->a\v,v->transpose(a)\v,v->adjoint(a)\v,v->conj(a)\v,eltype(vn),eltype(vm);fMVok=true,iMVok=true)
    B=joLinearFunctionAll(m,n,v->b*v,v->transpose(b)*v,v->adjoint(b)*v,v->conj(b)*v,v->b\v,v->transpose(b)\v,v->adjoint(b)\v,v->conj(b)\v,eltype(vn),eltype(vm))
else
    tname="loop $t for complex($m,$n)"
    a=rand(ComplexF64,m,n)
    b=rand(ComplexF64,m,n)
    vn=rand(ComplexF64,n)
    mvn=rand(ComplexF64,n,2)
    vm=rand(ComplexF64,m)
    mvm=rand(ComplexF64,m,2)
    A=joLinearFunctionAll(m,n,v->a*v,v->transpose(a)*v,v->adjoint(a)*v,v->conj(a)*v,v->a\v,v->transpose(a)\v,v->adjoint(a)\v,v->conj(a)\v,eltype(vn),eltype(vm))
    B=joLinearFunctionAll(m,n,v->b*v,v->transpose(b)*v,v->adjoint(b)*v,v->conj(b)*v,v->b\v,v->transpose(b)\v,v->adjoint(b)\v,v->conj(b)\v,eltype(vn),eltype(vm);fMVok=true,iMVok=true)
end
c=a-b
C=A-B

verbose && println("$tsname $tname")
    @testset "$tname A*B" begin
        @test size(C)==size(c)
        @test length(C)==length(c)
        @test norm(jo_full(C)-c)<joTol
        @test norm(elements(C)-c)<joTol
        @test norm(C)-norm(c)<joTol
        for i=1:2
            @test norm(C,i)-norm(c,i)<joTol
        end
        @test norm(C,Inf)-norm(c,Inf)<joTol
        @test norm(C)-norm(c)<joTol
        for i=1:2
            @test norm(C,i)-norm(c,i)<joTol
        end
        @test norm(C,Inf)-norm(c,Inf)<joTol
        @test norm(elements(adjoint(C))-adjoint(c))<joTol
        @test norm(elements(transpose(C))-transpose(c))<joTol
        @test norm(elements(+C)-(+c))<joTol
        @test norm(elements(-C)-(-c))<joTol
        @test norm(C*vn-c*vn)<joTol
        @test norm(C*mvn-c*mvn)<joTol
        @test norm(adjoint(C)*vm-adjoint(c)*vm)<joTol
        @test norm(adjoint(C)*mvm-adjoint(c)*mvm)<joTol
        @test norm(transpose(C)*vm-transpose(c)*vm)<joTol
        @test norm(transpose(C)*mvm-transpose(c)*mvm)<joTol
        @test norm(adjoint(transpose(C))*vn-adjoint(transpose(c))*vn)<joTol
        @test norm(adjoint(transpose(C))*mvn-adjoint(transpose(c))*mvn)<joTol
        @test norm(transpose(adjoint(C))*vn-transpose(adjoint(c))*vn)<joTol
        @test norm(transpose(adjoint(C))*mvn-transpose(adjoint(c))*mvn)<joTol
        @test norm(elements(transpose(C)*C)-(transpose(c)*c))<joTol
        @test norm(elements(C*transpose(C))-(c*transpose(c)))<joTol
        @test norm(elements(adjoint(C)*C)-(adjoint(c)*c))<joTol
        @test norm(elements(C*adjoint(C))-(c*adjoint(c)))<joTol
        if issquare(C)
            @test norm(C\vn-c\vn)<joTol
            @test norm(C\mvn-c\mvn)<joTol
            @test norm(adjoint(C)\vm-adjoint(c)\vm)<joTol
            @test norm(adjoint(C)\mvm-adjoint(c)\mvm)<joTol
            @test norm(transpose(C)\vm-transpose(c)\vm)<joTol
            @test norm(transpose(C)\mvm-transpose(c)\mvm)<joTol
        end
    end

end # end test loop
end


text:
function plothistogram(res, model = 0; removelowfrequencies = true, plotsbackend = nothing)
  if plotsbackend != nothing
      plotsbackend()
  end
  if removelowfrequencies == true
      fmin = res.ABCsetup.Models[1].constants[9]
  else
      fmin = 0.0
  end

  model = model + 1

  DFres = res.Posterior[model].MeanHistogram
  DFres[:VAF] = collect(0.01:0.01:1.0)
  DF = DataFrame(VAF = res.VAF)
  DFres = DFres[1:75, :]
  DFreshist = deepcopy(DFres)
  DFres = DFres[DFres[:VAF].>=fmin, :]

  if model - 1 == 0
      postcolor = RGBA(0.0, 76/255, 153/255) # plot neutral in blue and selection in red as in paper
  else
      postcolor = RGBA(0.75, 0.3, 0.3)
  end

  bar(DFreshist[:VAF], DFreshist[:truecounts], linecolor = RGBA(0.431, 0.431, 0.431, 0.9),
  fillcolor = RGBA(0.431, 0.431, 0.431, 0.9))
  plot!(DFres[:VAF], DFres[:mean], color = postcolor, w = 2)
  plot!(DFres[:VAF], DFres[:upperq95], fillrange = DFres[:lowerq95],
               fillalpha = 0.5,
               fillcolor = postcolor, linecolor = false,
               markerstrokecolor=:white, titlefont = font(8, "Calibri"), ytickfont = font(6, "Calibri"), xtickfont = font(6, "Calibri"), legend = false, grid = false,
               yaxis = ("Counts"), xaxis = ("VAF"))
  if model == 2
      xint = median(res.Posterior[model].Parameters[:frequency], Weights(res.Posterior[model].Parameters[:weight])) *
      median(res.Posterior[model].Parameters[:cellularity], Weights(res.Posterior[model].Parameters[:weight]))
      xint = xint/2
      vline!([xint], line=(2,:dash,0.6,:black))
  end

  if model == 3
      xint1 = median(res.Posterior[model].Parameters[:frequency1], Weights(res.Posterior[model].Parameters[:weight])) *
      median(res.Posterior[model].Parameters[:cellularity], Weights(res.Posterior[model].Parameters[:weight]))
      xint1 = xint1 / 2
      xint2 = median(res.Posterior[model].Parameters[:frequency2], Weights(res.Posterior[model].Parameters[:weight])) *
      median(res.Posterior[model].Parameters[:cellularity], Weights(res.Posterior[model].Parameters[:weight]))
      xint12 = xint2 / 2
      vline!([xint1, xint2], line=(2,:dash,0.6,[:black, :black]))
  end

end

function plotmodelposterior(res; plotsbackend = nothing)
    if plotsbackend != nothing
        plotsbackend()
    end
    DF = DataFrame(Model = map(x -> "$x", res.ModelProb[:Model]), Probability = res.ModelProb[:Probability])

    Plots.bar(DF[:Model], DF[:Probability],
    title="Model Probabilities", yaxis = ("Probability"), xaxis = ("Number of subclones"),
    linecolor = :white, fillcolor = RGBA(0.5, 0.5, 0.5, 0.8),
    markerstrokecolor=:white, titlefont = font(10, "Calibri"), ytickfont = font(8, "Calibri"), xtickfont = font(8, "Calibri"), legend = false, grid = false)
end

function plotparameterposterior(res, model = 1; plotsbackend = nothing)
    if plotsbackend != nothing
        plotsbackend()
    end
    if model == 0
        model = model + 1
        Plots.histogram(Matrix(res.Posterior[model].Parameters)[:, 1:3], nbins = 20, layout = 3, weights = Array(res.Posterior[model].Parameters[:weight]),
        title=["Mutation rate" "# Clonal mutations" "Cellularity"],
        linecolor = :white, fillcolor = RGBA(0.75, 0.3, 0.3),
        markerstrokecolor=:white, titlefont = font(10, "Calibri"), ytickfont = font(6, "Calibri"), xtickfont = font(6, "Calibri"), legend = false)
    elseif model == 1
        model = model + 1
        Plots.histogram(Matrix(res.Posterior[model].Parameters)[:, 1:7],
        nbins = 20, layout = 7,weights = Array(res.Posterior[model].Parameters[:weight]),
        title=["Mutation rate" "# Clonal mutations" "s" "t" "Cellularity" "Subclone frequency" "# Mutations in subclone"],
        linecolor = :white, fillcolor = RGBA(0.75, 0.3, 0.3),
        markerstrokecolor=:white, titlefont = font(10, "Calibri"),
        xtickfont = font(6, "Calibri"), ytickfont = font(6, "Calibri"), legend = false)
    elseif model == 2
        model = model + 1
        Plots.histogram(Matrix(res.Posterior[model].Parameters)[:, 1:11],
        nbins = 20, layout = 11,weights = Array(res.Posterior[model].Parameters[:weight]),
        title=["Mutation rate" "# Clonal mutations" "s1" "t1" "s2" "t2" "Cellularity" "Subclone 1 frequency" "Subclone 2 frequency" "# Mutations in \nsubclone 1" "# Mutations in subclone 2"],
        linecolor = :white, titlefont = font(8, "Calibri"),xtickfont = font(6, "Calibri"), ytickfont = font(6, "Calibri"), fillcolor = RGBA(0.75, 0.3, 0.3), legend = false)
    end
end

"""
    saveallplots(res; <keyword arguments>)

Create and save all plots. For each model, the VAF histogram with model results overlayed is plotted as is the posterior distributions for parameters and models. Function takes a fitABCmodels output type.
...
## Arguments
- `resultsdirectory = "output"`: Directory to save the plots.
- `outputformat = ".pdf"`: Format to save the plots as. Default is pdf, other common formats are available such png etc.
- `plotsbackend = nothing`: Backend to use for plotting in Plots.jl
...
"""
function saveallplots(res; resultsdirectory = "output", outputformat = ".pdf", plotsbackend = nothing, savepopulations = false, popnum = 1)
    if plotsbackend != nothing
        plotsbackend()
    end
  sname = res.SampleName
  if savepopulations
      dir = joinpath(resultsdirectory, res.SampleName, "populations", "population_$popnum")
      makeplotsdirectories(dir)
  else
      dir = joinpath(resultsdirectory, res.SampleName, "finalpopulation")
      makedirectory(resultsdirectory)
      makeplotsdirectories(dir)
  end
  p = plotmodelposterior(res, plotsbackend = plotsbackend)
  savefig(joinpath(dir, "plots", "$(sname)-modelposterior$(outputformat)"))

  model = 0
  for post in res.Posterior
    if post.Probability > 0.0
      p = plothistogram(res, model, plotsbackend = plotsbackend)
      savefig(joinpath(dir, "plots", "$(sname)-histogram-$(model)clone$(outputformat)"))
      p = plotparameterposterior(res, model, plotsbackend = plotsbackend);
      savefig(joinpath(dir, "plots", "$(sname)-posterior-$(model)clone$(outputformat)"))
    end
    model = model + 1
  end

end


text:
#
# Handles NBody problem states.
#

"""
Stores the state of each body in the NBody problem.
"""
struct Body{F<:AbstractFloat} <: AbstractBody

    r::SVector{3, Unitful.Length{F}}
    v::SVector{3, Unitful.Velocity{F}}
    m::Unitful.Mass{F}
    name::String

    function Body(r::VR, v::VV, m::M, name::String="") where {
            T1  <: Real,
            T2  <: Real,
            T3  <: Real,
            R   <: Unitful.Length{T1},
            V   <: Unitful.Velocity{T2},
            VR  <: AbstractVector{R},
            VV  <: AbstractVector{V},
            M   <: Unitful.Mass{T3}
        }
        if length(r) ≢ length(v) ≢ 3
            throw(ArgumentError("The `Body` constructor requires 3 element vectors for position `r` and velocity `v`"))
        else
            T = promote_type(T1, T2, T3)
            if !(T <: AbstractFloat)
                @warn "Non-float parameters provided. Defaulting to Float64."
                T = Float64
            end
            return new{T}(SVector{3}(T.(r)), SVector{3}(T.(v)), T(m), name)
        end
    end

    function Body(r::VR, v::VV, m::M, name::String="") where {
            R   <: Real,
            V   <: Real,
            VR  <: AbstractVector{R},
            VV  <: AbstractVector{V},
            M   <: Real
        }
        @warn "No units provided! Assuming km, km/s, and kg."
        return Body(r * u"km", v * u"km/s", m * u"kg", name)
    end
end

Base.convert(::Type{T}, b::Body) where {T<:AbstractFloat} = Body(T.(b.r), T.(b.v), T(b.m), b.name)
Base.promote(::Type{Body{A}}, ::Type{Body{B}}) where {A<:AbstractFloat, B<:AbstractFloat} = Body{promote_type(A,B)}
Core.Float16(o::Body) = convert(Float16, o)
Core.Float32(o::Body) = convert(Float32, o)
Core.Float64(o::Body) = convert(Float64, o)
Base.MPFR.BigFloat(o::Body) = convert(BigFloat, o)

"""
Describes a system of `n` `NBodyStates`'s.
"""
struct NBodySystem{N, T<:AbstractFloat} <: AbstractOrbitalSystem

    bodies::SVector{N, Body{T}}

    function NBodySystem(b::B...) where B<:Body
        N = length(b)
        bodies = promote(b...)
        T = promote_type([typeof(bᵢ).parameters[1] for bᵢ ∈ b]...)
        return new{length(b),T}(SVector{N, Body{T}}(bodies...))
    end
    NBodySystem(b::VB) where {B<:Body, VB<:AbstractVector{B}} = NBodySystem(b...)

end

"""
The `length` of an `NBodySystem` is the number of bodies in the system.
"""
Base.@pure Base.length(sys::NBodySystem{N,T}) where N where T = N

"""
The n-th `index` of an `NBodySystem` is the n-th body in the system.
"""
Base.getindex(sys::NBodySystem, i) = sys.bodies[i]
Base.convert(::Type{T}, sys::NBodySystem) where {T<:AbstractFloat} = NBodySystem(convert.(T, sys.bodies)...)
Base.promote(::Type{NBodySystem{N, A}}, ::Type{NBodySystem{N, B}}) where {A<:AbstractFloat, B<:AbstractFloat,N} = NBodySystem{promote_type(A,B)}
Core.Float16(o::NBodySystem) = convert(Float16, o)
Core.Float32(o::NBodySystem) = convert(Float32, o)
Core.Float64(o::NBodySystem) = convert(Float64, o)
Base.MPFR.BigFloat(o::NBodySystem) = convert(BigFloat, o)


text:
# Counting Inversions
#
# This implementation assumes all the values in the collection are distinct.

function count_inversions_naive(a::AbstractVector)::Int
    cnt = 0
    n = length(a)
    for i in 1:(n - 1)
        cnt += sum(a[i] .> @view a[(i + 1):n])
    end
    cnt
end

"Merge sorts `a` and returns its number of inversions."
function count_inversions_mergesort!(a::AbstractVector)::Int
    aux = zeros(eltype(a), length(a))
    ms_inv_r!(a, aux)
end

function ms_inv_r!(a::AbstractVector, aux::AbstractVector)::Int
    n = length(a)
    if n <= 1
        return 0
    end

    m = div(n, 2)
    cnt = @views ms_inv_r!(a[1:m], aux[1:m]) +
                 ms_inv_r!(a[(m + 1):n], aux[(m + 1):n])

    i = 1
    j = m + 1
    for k in 1:n
        if i > m
            aux[k] = a[j]
            j += 1
        elseif j > n
            aux[k] = a[i]
            i += 1
        elseif a[i] < a[j]
            aux[k] = a[i]
            i += 1
        else
            aux[k] = a[j]
            j += 1
            cnt += m - i + 1
        end
    end
    a .= aux

    cnt
end

# Plargarized from https://stackoverflow.com/a/23201616
function count_inversions_fenwick(a::AbstractVector)::Int
    n = length(a)
    cnt = 0
    count_tree = zeros(Int, n)
    ranks = Dict(v => i for (i, v) in enumerate(sort(a)))
    for x in Iterators.reverse(a)
        r = ranks[x]
        i = r - 1
        while i > 0
            cnt += count_tree[i]
            i -= i & (-i)
        end

        i = r
        while i <= n
            count_tree[i] += 1
            i += i & (-i)
        end
    end
    cnt
end

#===== Unit Tests =====#

using Test
using Random: seed!, shuffle!

@testset "Count inversions" begin
    n = 1000
    seed!()
    a = collect(1:n)
    shuffle!(a)

    n_inversions_naive = count_inversions_naive(a)
    n_inversions_fenwick = count_inversions_fenwick(a)
    n_inversions_mergesort = count_inversions_mergesort!(a)
    @test n_inversions_naive == n_inversions_mergesort == n_inversions_fenwick
    @test issorted(a)
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra
using Catlab.Theories
using Catlab.Syntax
using Catlab.GAT
using Catlab.Present
using Catlab.Programs
using Catlab.Graphics
import Catlab.Graphics: Graphviz
import Catlab.Theories: dom, codom, id, compose, otimes, braid, munit, mcopy, delete, pair, proj1, proj2, ⊗, ⋅, σ

#export CircuitDom, Circuit, show_diagram, apply, iNOT, iAND, iOR, Impl, Circuits

# Helper function to show graphviz diagrams
show_diagram(d::WiringDiagram) = to_graphviz(
    add_junctions(d),
    orientation=LeftToRight,
    labels=false, label_attr=:xlabel,
    node_attrs=Graphviz.Attributes(
        :fontname => "Courier",
    ),
    edge_attrs=Graphviz.Attributes(
        :fontname => "Courier",
    )
)

# The skeleton of n-dimensional bool spaces which are the input/output types of logic circuits
struct CircuitDom
    n::Int
end

# Implement a logic circuit as a function from Vector{Bool} to Vector{Bool}
struct Circuit
    dom::CircuitDom
    codom::CircuitDom
    impl::Function # Vector{Bool} of size dom -> Vector{Bool} of size codom
end

# Convenience constructor
Circuit(nins::Int, nouts::Int, impl) = Circuit(CircuitDom(nins), CircuitDom(nouts), impl)

# Apply a circuit to the specified input vector
apply(circuit::Circuit, input::Vector{Bool})::Vector{Bool} = circuit.impl(input)

# Basic logic gate implementations
iNOT = Circuit(1, 1, x->map(b->!b, x))

iAND = Circuit(2, 1, x->[x[1] && x[2]])

iOR = Circuit(2, 1, x->[x[1] || x[2]])

#iδ = Circuit(1, 2, x->[x[1], x[1]])

# Make circuits compositional by implementing them as an SMC instance
@instance CartesianCategory{CircuitDom, Circuit} begin
    id(A::CircuitDom) = Circuit(A,A, x->x)
    dom(f::Circuit) = f.dom
    codom(f::Circuit) = f.codom

    # Circuit composition is just function composition
    compose(f::Circuit, g::Circuit) = begin
        @assert(f.codom == g.dom)
        return Circuit(f.dom, g.codom, x->g.impl(f.impl(x)))
    end

    otimes(A::CircuitDom, B::CircuitDom) = CircuitDom(A.n + B.n)
    # Monoidal product of circuits runs both circuits (TODO: run them in parallel to improve performance) and concatenates the results
    otimes(f::Circuit, g::Circuit) = begin
        impl = x -> vcat(f.impl(x[1:f.dom.n]), g.impl(x[f.dom.n + 1:end]))
        return Circuit(f.dom.n + g.dom.n, f.codom.n + g.codom.n, impl)
    end

    # A symmetric braiding just swaps the A and B parts of the input vector
    braid(A::CircuitDom, B::CircuitDom) = begin
        impl = x -> vcat(x[A.n+1:end], x[1:A.n])
        n = A.n + B.n
        Circuit(n, n, impl)
    end

    # Monoidal unit is a 0-dimensional bool space (i.e. a point or empty list)
    munit(::Type{CircuitDom}) = CircuitDom(0)

    # Stuff for Cartesian Category
    mcopy(A::CircuitDom) = Circuit(A.n, A.n+A.n, x->vcat(x,x))
    delete(A::CircuitDom) = Circuit(A.n, 0, x->Bool[])

    pair(f::Circuit, g::Circuit) = begin
        @assert(f.dom == g.dom)
        return Circuit(f.dom.n, f.codom.n+g.codom.n, x->vcat(f.impl(x), g.impl(x)))
    end
    proj1(A::CircuitDom, B::CircuitDom) = Circuit(A.n+B.n, A.n, x->x[1:A.n])
    proj2(A::CircuitDom, B::CircuitDom) = Circuit(A.n+B.n, B.n, x->x[A.n+1:end])
end

# Syntactic representation of logic circuits

@present Circuits(FreeCartesianCategory) begin
    B::Ob
    NOT::Hom(B,B)
    AND::Hom(B⊗B,B)
    OR::Hom(B⊗B,B)
    #XOR::Hom(B⊗B,B)
    #FA::Hom(B⊗B⊗B,B⊗B)
end

# Functor from circuit diagrams to circuit implementations
gens = Dict(
    Circuits[:B]=>CircuitDom(1),
    Circuits[:NOT]=>iNOT,
    Circuits[:AND]=>iAND,
    Circuits[:OR]=>iOR
)
Impl(expr) = functor((CircuitDom, Circuit), expr, generators=gens)


text:
module Xerces_jll

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 0
end

if VERSION < v"1.3.0-rc4"
    # We lie a bit in the registry that JLL packages are usable on Julia 1.0-1.2.
    # This is to allow packages that might want to support Julia 1.0 to get the
    # benefits of a JLL package on 1.3 (requiring them to declare a dependence on
    # this JLL package in their Project.toml) but engage in heroic hacks to do
    # something other than actually use a JLL package on 1.0-1.2.  By allowing
    # this package to be installed (but not loaded) on 1.0-1.2, we enable users
    # to avoid splitting their package versions into pre-1.3 and post-1.3 branches
    # if they are willing to engage in the kinds of hoop-jumping they might need
    # to in order to install binaries in a JLL-compatible way on 1.0-1.2. One
    # example of this hoop-jumping being to express a dependency on this JLL
    # package, then import it within a `VERSION >= v"1.3"` conditional, and use
    # the deprecated `build.jl` mechanism to download the binaries through e.g.
    # `BinaryProvider.jl`.  This should work well for the simplest packages, and
    # require greater and greater heroics for more and more complex packages.
    error("Unable to import Xerces_jll on Julia versions older than 1.3!")
end

using Pkg, Pkg.BinaryPlatforms, Pkg.Artifacts, Libdl
import Base: UUID

# We put these inter-JLL-package API values here so that they are always defined, even if there
# is no underlying wrapper held within this JLL package.
const PATH_list = String[]
const LIBPATH_list = String[]

# Load Artifacts.toml file
artifacts_toml = joinpath(@__DIR__, "..", "Artifacts.toml")

# Extract all platforms
artifacts = Pkg.Artifacts.load_artifacts_toml(artifacts_toml; pkg_uuid=UUID("637d83c4-b86a-5d90-b82d-5cf0573a8cfc"))
platforms = [Pkg.Artifacts.unpack_platform(e, "Xerces", artifacts_toml) for e in artifacts["Xerces"]]

# Filter platforms based on what wrappers we've generated on-disk
filter!(p -> isfile(joinpath(@__DIR__, "wrappers", replace(triplet(p), "arm-" => "armv7l-") * ".jl")), platforms)

# From the available options, choose the best platform
best_platform = select_platform(Dict(p => triplet(p) for p in platforms))

# Silently fail if there's no binaries for this platform
if best_platform === nothing
    @debug("Unable to load Xerces; unsupported platform $(triplet(platform_key_abi()))")
else
    # Load the appropriate wrapper.  Note that on older Julia versions, we still
    # say "arm-linux-gnueabihf" instead of the more correct "armv7l-linux-gnueabihf",
    # so we manually correct for that here:
    best_platform = replace(best_platform, "arm-" => "armv7l-")
    include(joinpath(@__DIR__, "wrappers", "$(best_platform).jl"))
end

end  # module Xerces_jll


text:
# Set up the environment to run this example. Make sure you're within the folder that this
# file lives in.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Plots, Random, Stheno, Flux, Zygote
rng = MersenneTwister(5)
T = Float64

# define the step function
# step_func(x) = 0.0 if x<=0, 1.0 if x>0
function step_func(x)
    ϵ=T(0.01)*randn(rng, T)
    return x>zero(T) ? one(T) + ϵ : zero(T) + ϵ
end

# prepare data
## training data, drawn from Normal(0, 1)
train_X = randn(rng, T, 100)
## Flux requires data to be at least two dimension
Xtrain = reshape(train_X, 1, :)
## evaluate to get the function value
train_y = step_func.(train_X)

## test data drawn from Uniform(-5, 5)
test_X = collect(-5.0:0.02:5.0)
Xtest = reshape(test_X, 1, :)
test_y = step_func.(test_X)

# Build an MLP that transform the input data
# for Flux usage, please refer to: https://fluxml.ai/Flux.jl/stable/
mlp = Chain(Dense(1, 6, relu), Dense(6, 2, relu)) |> (T==Float32 ? f32 : f64)
# extract MLP parameters
θ_mlp = params(mlp)

# Build a GP model with Stheno
# here we consider using anisotropic EQ kernel, this kernel contains to hyperparameters: length scale l & scaling factor γ
# these hyperparameters are positive, we add on this restriction by resetting them in log-scale
# for Stheno usage, please refer to: https://github.com/willtebbutt/Stheno.jl/tree/wct/example-revamp/examples/getting_started
logl = randn(rng, T, 2)
logγ = T[0.0]

function build_gp(logl, logγ)
    ard_eq_kernel = exp(T(2.0) * logγ[1]) * stretch(EQ(), exp.(-logl))
    return GP(T(0.0), ard_eq_kernel, GPC())
end

# Since we always assume our data to be noisy, we model this noise by λ, also in log-scale
logλ = T[0.01]

# Collect MLP and GP parameters
ps = Params([logl, logγ, logλ, θ_mlp...])

# When training, we always specifies a loss function to optimize, for GP we use negative log-likelihood
# !!! NOTE:
# Stheno and Flux have slightly different conventions regarding their data. In particular
# - Flux assumes that you'll provide a `Matrix` of data, in which each column is an
#   observation.
# - Stheno assumes that you'll provide an `AbstractVector`, where each element is an
#   observation. To handle multi-dimensional inputs we have the `ColVecs` object, that
#   literally just wraps a `Matrix` and tells Stheno to pretend is a vector of vectors. This
#   is helpful to remove some ambiguities that arise if you don't do this.
function NLL(X, y)
    Z = mlp(X)
    gp = build_gp(logl, logγ)
    gp_Z = gp(ColVecs(Z), exp(T(2.0)*logλ[1]))
    return -logpdf(gp_Z, y)
end

# Training the overall model with Flux optimizers
using Flux.Optimise: update!

train_data = (Xtrain, train_y)
opt = ADAGrad()
nlls = []
for i in 1:1500
    nll = NLL(train_data...)
    push!(nlls, nll)
    gs = gradient(()->NLL(train_data...), ps)
    for p in ps
	update!(opt, p, gs[p])
    end
end

loss_plot = plot(xlabel="Epoches", ylabel="Negative log-likelihood", legend=false)
plot!(loss_plot, nlls)
png(loss_plot, "loss.png")

# Visualize the performance of our model
function predict(X, Xtrain, ytrain)
    Z = mlp(X); Ztrain = mlp(Xtrain)
    gp = build_gp(logl, logγ)
    noisy_prior = gp(ColVecs(Ztrain), exp(T(2.0)*logλ[1]))
    posterior = gp | Obs(noisy_prior, ytrain)
    return posterior(ColVecs(Z))
end

posterior = predict(Xtest, Xtrain, train_y)
post_dist = marginals(posterior)
pred_y = mean.(post_dist)
var_y = std.(post_dist)

predict_plot = plot(legend=true, xlabel="x", ylabel="y", ylim=(-0.5, 1.5))
plot!(predict_plot, test_X, pred_y, ribbons=3*var_y, st=:line, fillalpha=0.3, color=:blue, label="NNGP")
plot!(predict_plot, train_X, train_y, st=:scatter, color=:red, lw=5, label="Training set")
png(predict_plot, "predict.png")


text:
using Distributed

@everywhere using Pkg
@everywhere Pkg.activate(".")
@everywhere using LightGraphs

sample_graphs(n, p, k) = [erdos_renyi(n, p, seed=i) for i in 1:k]

function sample_graphs_t(n, p, k)
    gs = @sync @distributed for i in 1:k
        erdos_renyi(n, p, seed=i)
    end
    gs
end

println("processes: ", nworkers())
@time sample_graphs(10^5, 10^-4, 32)
@time sample_graphs(10^5, 10^-4, 32)
@time sample_graphs_t(10^5, 10^-4, 32)
@time sample_graphs_t(10^5, 10^-4, 32)


text:
function calibrate(x, poly)
    return sum(x.^i .* poly[i] for i in 1:length(poly))
end

function calibrate(x, shift, gain, gain²)
    return shift + x*gain + x^2*gain²
end

function calibrate(x::T, shift, gain) where T <: Integer
    return shift + (x+rand()*2-1)*gain
end

function calibrate(x::T, shift, gain) where T <: AbstractFloat
    return shift + (x+rand()*2-1)*gain
end

function ucalibrate(x, shift, gain)
    return shift + x*gain
end

function ucalibrate(x, shift, gain, gain²)
    return shift + x*gain + x^2 * gain²
end

function ucalibrate(x, shift, gain, gain², gain³)
    return shift + x*gain + x^2 * gain² + x^3*gain³
end

function gainshift!(m::Array{T, 1}, shift, gain, gain²) where T <: Integer
    tmp = zeros(T, size(m))
    first = firstindex(m)
    last = lastindex(m)
    @inbounds @fastmath for i in eachindex(m)
        j = unsafe_trunc.(T, calibrate.(i .+ 2 .*rand(m[i]) .- 1, shift, gain, gain²))
        for _j in j[first .≤ j .≤ last]
            tmp[_j] += 1
        end
    end
    m .= tmp
end

function gainshift!(m::Array{T, 1}, shift, gain) where T <: Integer
    tmp = zeros(T, size(m))
    first = firstindex(m)
    last = lastindex(m)
    @inbounds @fastmath for i in eachindex(m)
        j = unsafe_trunc.(T, ucalibrate.(i .+ 2 .*rand(m[i]) .- 1, shift, gain))
        for _j in j[first .≤ j .≤ last]
            tmp[_j] += 1
        end
    end
    m .= tmp
end

function gainshift(m::Array{T, 1}, shift, gain, gain²) where T <: Integer
    tmp = zeros(T, size(m))
    first = firstindex(m)
    last = lastindex(m)
    @inbounds @fastmath for i in eachindex(m)
        j = unsafe_trunc.(T, ucalibrate.(i .+ 2 .*rand(m[i]) .- 1, shift, gain, gain²))
        for _j in j[first .≤ j .≤ last]
            tmp[_j] += 1
        end
    end
    tmp
end

function gainshift(m::Array{T, 1}, shift, gain) where T <: Integer
    tmp = zeros(T, size(m))
    first = firstindex(m)
    last = lastindex(m)
    @inbounds @fastmath for i in eachindex(m)
        j = unsafe_trunc.(T, ucalibrate.(i .+ 2 .*rand(m[i]) .- 1, shift, gain))
        for _j in j[first .≤ j .≤ last]
            tmp[_j] += 1
        end
    end
    tmp
end

function gainshift(m::Array{T, 1}, shift, gain, gain², gain³) where T <: Integer
    tmp = zeros(T, size(m))
    first = firstindex(m)
    last = lastindex(m)
    @inbounds @fastmath for i in eachindex(m)
        j = unsafe_trunc.(T, ucalibrate.(i .+ 2 .*rand(m[i]) .- 1, shift, gain, gain², gain³))
        for _j in j[first .≤ j .≤ last]
            tmp[_j] += 1
        end
    end
    tmp
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
using Distances, NearestNeighbors
mns = TransD_GP.init(opt, log10λ)
##
mnsold = deepcopy(mns)
log10λold = deepcopy(log10λ)
TransD_GP.birth!(log10λ, optlog10λ, mns, opt.δ, doall=true)
TransD_GP.sync_model!(mns, opt)
λ² = log10λ.fstar
idxs = TransD_GP.gettrainidx(opt.kdtree, mnsold.xtrain, mnsold.n)
ftest, = GP.GPfit(K, mnsold.ftrain[:,1:mnsold.n], mnsold.xtrain[:,1:mnsold.n],
        opt.xall, λ², λ²[:,idxs], opt.δ, p=2, demean=false, nogetvars=true,
        my=mean(opt.fbounds, dims=2))
##
m = deepcopy(log10λold)
m.fstar = 0.5log10.(m.fstar)
fig,ax = plt.subplots(1,2, sharex=true, sharey=true)
vmin, vmax = minimum(m.fstar[1,:]), maximum(m.fstar[1,:])
im1 = ax[1].imshow(reshape(m.fstar[1,:],length(y), length(x)), extent=[x[1],x[end],y[end],y[1]],
    vmin=vmin, vmax=vmax)
colorbar(im1, ax=ax[1])
ax[1].scatter(m.xtrain[1,1:m.n], m.xtrain[2,1:m.n],marker="+",c="r")
ax[1].scatter(m.xtrain[1,1:m.n], m.xtrain[2,1:m.n],c=m.ftrain[1,1:m.n], alpha=0.8,
vmin=vmin, vmax=vmax)
vmin, vmax = minimum(m.fstar[2,:]), maximum(m.fstar[2,:])
im2 = ax[2].imshow(reshape(m.fstar[2,:],length(y), length(x)), extent=[x[1],x[end],y[end],y[1]],
    vmin=vmin, vmax=vmax)
colorbar(im2, ax=ax[2])
ax[2].scatter(m.xtrain[1,1:m.n], m.xtrain[2,1:m.n],marker="+",c="r")
ax[2].scatter(m.xtrain[1,1:m.n], m.xtrain[2,1:m.n],c=m.ftrain[2,1:m.n], alpha=0.8,
    vmin=vmin, vmax=vmax)
##
fig,ax = plt.subplots(1,2, sharex=true, sharey=true)
vmin, vmax = minimum(mnsold.fstar), maximum(mnsold.fstar)
im1 = ax[1].imshow(reshape(mnsold.fstar,length(y), length(x)), extent=[x[1],x[end],y[end],y[1]],
    vmin=vmin, vmax=vmax)
colorbar(im1, ax=ax[1])
ax[1].scatter(mnsold.xtrain[1,1:mnsold.n], mnsold.xtrain[2,1:mnsold.n],marker="+",c="r")
ax[1].scatter(mnsold.xtrain[1,1:mnsold.n], mnsold.xtrain[2,1:mnsold.n],c=mnsold.ftrain[1:mnsold.n],
        alpha=0.8, vmin=vmin, vmax=vmax)
vmin, vmax = minimum(mns.fstar), maximum(mns.fstar)
im2 = ax[2].imshow(reshape(mns.fstar,length(y), length(x)), extent=[x[1],x[end],y[end],y[1]],
    vmin=vmin, vmax=vmax)
colorbar(im2, ax=ax[2])
ax[2].scatter(mns.xtrain[1,1:mns.n], mns.xtrain[2,1:mns.n],marker="+",c="r")
ax[2].scatter(mns.xtrain[1,1:mns.n], mns.xtrain[2,1:mns.n],c=mns.ftrain[1:mns.n], alpha=0.8,
    vmin=vmin, vmax=vmax)


text:
module TestIteration
    using Base.Test, DataTables, Compat

    dv = NullableArray(Nullable{Int}[1, 2, Nullable()])
    dm = NullableArray([1 2; 3 4])
    dt = NullableArray(zeros(2, 2, 2))

    dt = DataTable(A = NullableArray(1:2), B = NullableArray(2:3))

    for row in eachrow(dt)
        @test isa(row, DataTableRow)
        @test isequal(row[:B]-row[:A], Nullable(1))

        # issue #683 (https://github.com/JuliaStats/DataFrames.jl/pull/683)
        @test typeof(collect(row)) == @compat Array{Tuple{Symbol, Any}, 1}
    end

    for col in eachcol(dt)
        @test isa(col, @compat Tuple{Symbol, NullableVector})
    end

    @test isequal(map(x -> minimum(convert(Array, x)), eachrow(dt)), Any[1,2])
    @test isequal(map(minimum, eachcol(dt)), DataTable(A = Nullable{Int}[1], B = Nullable{Int}[2]))

    row = DataTableRow(dt, 1)

    row[:A] = 100
    @test isequal(dt[1, :A], Nullable(100))

    row[1] = 101
    @test isequal(dt[1, :A], Nullable(101))

    dt = DataTable(A = NullableArray(1:4), B = NullableArray(["M", "F", "F", "M"]))

    s1 = view(dt, 1:3)
    s1[2,:A] = 4
    @test isequal(dt[2, :A], Nullable(4))
    @test isequal(view(s1, 1:2), view(dt, 1:2))

    s2 = view(dt, 1:2:3)
    s2[2, :B] = "M"
    @test isequal(dt[3, :B], Nullable("M"))
    @test isequal(view(s2, 1:1:2), view(dt, [1,3]))

    # @test_fail for x in dt; end # Raises an error
end


text:
using ShardSort
using Test

@testset "ShardSort.jl" begin
    # Write your tests here.
end

max_stars_repo_name: pthariensflame/shard-sort

max_issues_repo_name: pthariensflame/shard-sort

max_forks_repo_name: pthariensflame/shard-sort


text:
# Synaptically connected STG model network built as in Prinz et al Nature Neurosci 2004

using OrdinaryDiffEq, Plots

include(joinpath(@__DIR__, "prinz_kinetics.jl"));
include(joinpath(@__DIR__, "prinz_synapses.jl"));

# Figure 3e Prinz (qualitatively close to expected output)
# AB/PD 2
@named ABPD = Soma([NaV(100mS/cm^2),
                    CaT(2.5mS/cm^2),
                    CaS(6mS/cm^2),
                    KA(50mS/cm^2),
                    KCa(5mS/cm^2),
                    Kdr(100mS/cm^2),
                    H(.01mS/cm^2),
                    leak(0mS/cm^2)],
                    gradients, area = area, V0 = -50mV, aux = [calcium_conversion]);

ABPD.sys # display

# LP 4
@named LP = Soma([NaV( 100mS/cm^2),
                  CaT(   0mS/cm^2),
                  CaS(   4mS/cm^2),
                  KA(   20mS/cm^2),
                  KCa(   0mS/cm^2),
                  Kdr(  25mS/cm^2),
                  H(   .05mS/cm^2),
                  leak(.03mS/cm^2)],
                  gradients, area = area, V0 = -50mV, aux = [calcium_conversion]);

LP.sys # display

# PY 1
@named PY = Soma([NaV( 100mS/cm^2),
                  CaT( 2.5mS/cm^2),
                  CaS(   2mS/cm^2),
                  KA(   50mS/cm^2),
                  KCa(   0mS/cm^2),
                  Kdr( 125mS/cm^2),
                  H(   .05mS/cm^2),
                  leak(.01mS/cm^2)],
                  gradients, area = area, V0 = -50mV, aux = [calcium_conversion]);

PY.sys # display

topology = [ABPD => (LP, Glut(30nS)),
            ABPD => (LP, Chol(30nS)),
            ABPD => (PY, Glut(10nS)),
            ABPD => (PY, Chol(3nS)),
            LP   => (ABPD, Glut(30nS)),
            LP   => (PY, Glut(1nS)),
            PY   => (LP, Glut(30nS))];

network = Network([ABPD, LP, PY], topology)

t = 10000
sim = Simulation(network, time = t*ms)
solution = solve(sim, Rosenbrock23())

# Plot at 5kHz sampling
fig = plot(solution; plotdensity=Int(t*5), size=(1200,800), vars = [ABPD.sys.Vₘ, LP.sys.Vₘ, PY.sys.Vₘ])
fig

# Uncomment and eval `png(...)` to save as PNG
# png(fig, "figure3e_simulated")


text:
module PeriodicSystems

#using SLICOTMath
using SLICOT_jll
using DescriptorSystems
using LinearAlgebra
using FFTW
using MatrixEquations
using MatrixPencils
using Polynomials
using Random
using Interpolations
using Symbolics
#using StaticArrays
#using DifferentialEquations
using OrdinaryDiffEq
using IRKGaussLegendre
using Primes
using ApproxFun

include("SLICOTtools.jl")
using .SLICOTtools: mb03vd!, mb03vy!, mb03bd!, mb03wd!, mb03vw!

import LinearAlgebra: BlasInt, BlasFloat, BlasReal, BlasComplex, copy_oftype, transpose, adjoint, opnorm, normalize, rdiv!
import Base: +, -, *, /, \, (==), (!=), ^, isapprox, iszero, convert, promote_op, size, length, ndims,
             hcat, vcat, hvcat, inv, show, lastindex, require_one_based_indexing, print, show, one, zero, eltype
import MatrixPencils: isregular, rmeval
import MatrixEquations: sylvd2!, luslv!
import DescriptorSystems: isstable
import Polynomials: AbstractRationalFunction, AbstractPolynomial, poles, isconstant, variable, degree, pqs

export PeriodicStateSpace, pschur, pschur1, pschur2, phess, phess1, psreduc_reg, psreduc_fast, check_psim, mshift, pseig, tvmeval, hreval, tvstm
export ts2hr, ts2pfm, ts2ffm, pfm2hr, hr2psm, psm2hr, pm2pa, ffm2hr, pmaverage
export monodromy, psceig
export PeriodicArray, PeriodicMatrix
export PeriodicTimeSeriesMatrix, HarmonicArray, FourierFunctionMatrix, PeriodicFunctionMatrix,  PeriodicSymbolicMatrix
export isperiodic, isconstant, iscontinuous, islti, set_period
export mb03vd!, mb03vy!, mb03bd!, mb03wd!
export ps
export psaverage, psc2d, psmrc2d, psteval
export ps2fls, hr2bt, hr2btupd, phasemat, ps2frls, DiagDerOp, ps2ls
export pspole, pszero, isstable
export prlyap, pflyap, pslyapd, pdlyaps!, dpsylv2, pslyapdkr
export pmshift

abstract type AbstractDynamicalSystem end
abstract type AbstractLTISystem <: AbstractDynamicalSystem end
abstract type AbstractLTVSystem <: AbstractDynamicalSystem end
abstract type AbstractDescriptorStateSpace <: AbstractLTISystem end
abstract type AbstractPeriodicStateSpace <: AbstractLTVSystem end

abstract type AbstractPeriodicArray{Domain,T} end

include("types/PeriodicMatrices.jl")
include("types/PeriodicStateSpace.jl")
include("ps.jl")
include("conversions.jl")
include("pslifting.jl")
include("psanalysis.jl")
include("pslyap.jl")
include("pmops.jl")
include("psutils.jl")
end


text:
module ReadSPSS

using Compat
using Compat.Mmap, Missings, CategoricalArrays

export SPSSDataFrame, read_sav

include("headerrecord.jl")
include("variablerecord.jl")
include("valuelabels.jl")
include("documentrecord.jl")
include("extensionrecord.jl")
include("bytecodedata.jl")
include("readers.jl")

end # module


text:
module CreepLaw

# This implements viscous creep laws and routines to compute with them
#
# Note that various simple creep laws are defined in this file;
# more complex ones (such as DislocationCreep) are in separate files
# in the same directory
#
# In case you want to add new creep laws, have a look at how the ones
# here are implemented. Please add tests as well!

using Base: Float64
using Parameters, LaTeXStrings, Unitful
using ..Units
using GeoParams: AbstractMaterialParam
using BibTeX
import Base.show

abstract type AbstractCreepLaw <: AbstractMaterialParam end

export  ComputeCreepLaw_EpsII, ComputeCreepLaw_TauII,       # calculation routines
        CreepLawVariables,                                     # holds additional parameters required for calculations
        LinearViscous,
        PowerlawViscous

"""
 Struct that holds parameters required for creep law calculations (such as P,T, grain size etc.)

You can assign the struct as
```julia-repl
 p = CreepLawVariables(P=0.0, T=0.0, d=1.0)
```
where you can also pass vectors or arrays as values.

Note that, if needed, this can be extended, w/out interfering with existing calculation
"""
@with_kw struct CreepLawVariables
    P::GeoUnit  =   100.0MPa   # pressure
    T::GeoUnit  =   500.0C     # temperature
    d::GeoUnit  =   1.0cm      # grainsize
    f::GeoUnit  =   0.0MPa     # water-fugacity
end
include("DislocationCreep.jl")
# Linear viscous rheology ------------------------------------------------
"""
    LinearViscous(η=1e20Pa*s)

Defines a linear viscous creeplaw

The (isotopic) linear viscous rheology is given by
```math
    \\tau_{ij} = 2 \\eta \\dot{\\varepsilon}_{ij}
```
or
```math
    \\dot{\\varepsilon}_{ij}  = {\\tau_{ij}  \\over 2 \\eta }
```

where ``\\eta_0`` is the reference viscosity [Pa*s] at reference strain rate ``\\dot{\\varepsilon}_0``[1/s], and ``n`` the power law exponent [].
"""
@with_kw_noshow mutable struct LinearViscous <: AbstractCreepLaw
    equation::LaTeXString   =   L"\tau_{ij} = 2 \eta  \dot{\varepsilon}_{ij}"
    η::GeoUnit            =   1e20Pa*s                # viscosity
end

# Calculation routines for linear viscous rheologies
function ComputeCreepLaw_EpsII(TauII, s::LinearViscous, p::CreepLawVariables)
    @unpack η   = s

    return EpsII = TauII/(2.0*η);
end

function ComputeCreepLaw_TauII(EpsII, s::LinearViscous, p::CreepLawVariables)
    @unpack η   = s

    return TauII = 2.0*η*EpsII;
end

# Print info
function show(io::IO, g::LinearViscous)
    print(io, "Linear viscosity: η=$(g.η.val)")
end
#-------------------------------------------------------------------------

# Powerlaw viscous rheology ----------------------------------------------
"""
    PowerlawViscous(η0=1e18Pa*s, n=2.0NoUnits, ε0=1e-15/s)

Defines a power law viscous creeplaw as:

```math
    \\tau_{ij}^n  = 2 \\eta_0 \\left( \\dot{\\varepsilon}_{ij} \\over \\dot{\\varepsilon}_0 \\right)
```
where ``\\eta`` is the effective viscosity [Pa*s].
"""
@with_kw_noshow mutable struct PowerlawViscous <: AbstractCreepLaw
    η0::GeoUnit     =  1e18Pa*s            # reference viscosity
    n::GeoUnit      =  2.0*NoUnits         # powerlaw exponent
    ε0::GeoUnit     =  1e-15*1/s;          # reference strainrate
end

# Print info
function show(io::IO, g::PowerlawViscous)
    print(io, "Powerlaw viscosity: η0=$(g.η0.val), n=$(g.n.val), ε0=$(g.ε0.val) ")
end
#-------------------------------------------------------------------------

# Help info for the calculation routines
"""
    ComputeCreepLaw_EpsII(TauII, s:<AbstractCreepLaw, p::CreepLawVariables)

Returns the strainrate invariant ``\\dot{\\varepsilon}_{II}`` for a given deviatoric stress
invariant ``\\tau_{II}`` for any of the viscous creep laws implemented.
``p`` is a struct that can hold various parameters (such as ``P``, ``T``) that the creep law
may need for the calculations

```math
    \\dot{\\varepsilon}_{II}   = f( \\tau_{II} )
```

"""
ComputeCreepLaw_EpsII

"""
    ComputeCreepLaw_TauII(EpsII, s:<AbstractCreepLaw, p::CreepLawVariables)

Returns the deviatoric stress invariant ``\\tau_{II}`` for a given strain rate
invariant ``\\dot{\\varepsilon}_{II}`` for any of the viscous creep laws implemented.
``p`` is a struct that can hold various parameters (such as ``P``, ``T``) that the creep law
may need for the calculations

```math
    \\tau_{II}  = f( \\dot{\\varepsilon}_{II} )
```

"""
ComputeCreepLaw_TauII

end


text:
module TestFGenerators
using Test

@testset "$file" for file in sort([
    file for file in readdir(@__DIR__) if match(r"^test_.*\.jl$", file) !== nothing
])

    if file == "test_doctest.jl"
        if lowercase(get(ENV, "JULIA_PKGEVAL", "false")) == "true"
            @info "Skipping doctests on PkgEval."
            continue
        elseif VERSION >= v"1.5-"
            @info "Skipping doctests on Julia $VERSION."
            continue
        end
    end

    include(file)
end
end  # module


text:
using MesocosmMetabolism, CSV, Turing, DataFrames, MCMCChains, StatsPlots, Interpolations, ReverseDiff

data = CSV.read("/Users/elizabethmohr/Documents/MSU/RProjects/stream-analogue-mesocosms/data/DO.csv", DataFrame)
Q = 1.794 #L/hr
V = 27.0 # L
dayStart = 4

# Note: these priors are taken from the default values in streamMetabolizer,
# converted to units of hours and converted to a per volume basis
# by dividing by mean depth in SAMS of 0.15 m.
GPP_daily_mu = 1.72 # g/(L*hr) - only applies to 12 hours when lights are on
GPP_daily_sigma = 3.33
ER_daily_mu = -2.0 # g/(L*hr)
ER_daily_sigma = 2.0
k600_meanlog = -0.69 # hr^-1
k600_sdlog = 0.53
K02_conv = 1.0
err_obs_iid_sigma_scale = 0.03
err_proc_iid_sigma_scale = 5.0

setadbackend(:reversediff)

results = Dict()
for i in unique(data[:,:loggerID])[1:4]
    subset = filter(r -> r[:loggerID] == i, data)
    times = subset.modelTime
    PAR = Interpolations.ConstantInterpolation((times,), subset.PAR)
    OSat = Interpolations.interpolate((times,), subset."DO.sat", Interpolations.Gridded(Interpolations.Linear()))
    days = [convert(Int64, floor(times[i]/24.0 - (dayStart/24.0)) + 1) for i in 1:length(times)]
    day = Interpolations.ConstantInterpolation((times,), days)
    u0 = [subset.DO_mgL[1]]
    nDays = convert(Int64,ceil((maximum(times) - minimum(times))/24))
    problem = initialize_process_model(
        times = times,
        PAR = PAR,
        OSat = OSat,
        day = day,
        GPP_daily = fill(GPP_daily_mu, nDays),
        ER_daily = fill(ER_daily_mu, nDays),
        kGas = exp(k600_meanlog),
        Q = Q,
        V = V,
        dayStart = dayStart,
        u0 = u0)
    model = statisticalModel(
        nDays,
        GPP_daily_mu,
        GPP_daily_sigma,
        ER_daily_mu,
        ER_daily_sigma,
        k600_meanlog,
        k600_sdlog,
        K02_conv,
        err_obs_iid_sigma_scale,
        err_proc_iid_sigma_scale)

    sampleMe = model(DOdata = subset.DO_mgL, ODE = problem,times = times)
    chains = sample(sampleMe, NUTS(), MCMCThreads(), 1000, 3)
    results[i] = chains
end

gelmandiag(chains)

df = DataFrame(chains)

# GPP trace plot
@df df plot(
    :iteration,
    cols(columnindex(df,"GPP_daily[1]")),
    group = (:chain),
    lw = 3,
    xlabel = "Iteration",
    ylabel = "GPP",
    legendtitle = "Chain",
    xguidefontsize=15,
    yguidefontsize = 15)

# GPP posterior and prior
@df df density(
    cols(columnindex(df,"GPP_daily[1]")),
    group = (:chain),
    lw = 3,
    xlabel = "GPP",
    ylabel = "Probability Density",
    legendtitle = "Chain",
    xguidefontsize=15,
    yguidefontsize = 15)

plot!(Normal(GPP_daily_mu, GPP_daily_sigma),lw = 3)

# ER trace plot
@df df plot(
    :iteration,
    cols(columnindex(df,"ER_daily[1]")),
    group = (:chain),
    lw = 3,
    xlabel = "Iteration",
    ylabel = "ER",
    legendtitle = "Chain",
    xguidefontsize=15,
    yguidefontsize = 15)

# ER posterior and prior
@df df density(
    cols(columnindex(df,"ER_daily[1]")),
    group = (:chain),
    lw = 3,
    xlabel = "ER",
    ylabel = "Probability Density",
    legendtitle = "Chain",
    xguidefontsize= 15,
    yguidefontsize = 15)

plot!(Normal(ER_daily_mu, ER_daily_sigma),lw = 3)

# k600 trace plot
@df df plot(
    :iteration,
    cols(columnindex(df,"k600")),
    group = (:chain),
    lw = 3,
    xlabel = "Iteration",
    ylabel = "k600",
    legendtitle = "Chain",
    xguidefontsize=15,
    yguidefontsize = 15)

# k600 posterior and prior
@df df density(
    cols(columnindex(df,"k600")),
    group = (:chain),
    lw = 3,
    xlabel = "k600",
    ylabel = "Probability Density",
    legendtitle = "Chain",
    xguidefontsize= 15,
    yguidefontsize = 15)

plot!(LogNormal(k600_meanlog, k600_sdlog),lw = 3)

# Plot data with modeled DO from multiple samples
using DifferentialEquations
logger = 415148
MesocosmData = filter(r -> r[:loggerID] == logger, data)
df = DataFrame(results[logger])
times = MesocosmData.modelTime
PAR = Interpolations.ConstantInterpolation((times,), MesocosmData.PAR)
OSat = Interpolations.interpolate((times,), MesocosmData."DO.sat", Interpolations.Gridded(Interpolations.Linear()))
days = [convert(Int64, floor(times[i]/24.0 - (dayStart/24.0)) + 1) for i in 1:length(times)]
day = Interpolations.ConstantInterpolation((times,), days)
u0 = [MesocosmData.DO_mgL[1]]
nDays = convert(Int64,ceil((maximum(times) - minimum(times))/24))

DOPlot = StatsPlots.scatter(
    MesocosmData.modelTime,
    MesocosmData.DO_mgL,
    ylabel = "DO (mg/L)",
    xguidefontsize= 15,
    yguidefontsize = 15);

n = 100
rows = rand(1:nrow(df), n)

for k in 1:n
    GPP = Array(df[k, r"^GPP_daily"])
    ER = Array(df[k, r"^ER_daily"])
    k600 = df[k, "k600"]
    problem = initialize_process_model(
        times = MesocosmData.modelTime,
        PAR = PAR,
        OSat = OSat,
        day = day,
        GPP_daily = GPP,
        ER_daily = ER,
        kGas = k600,
        Q = Q,
        V = V,
        dayStart = dayStart,
        u0 = u0)
    sol = solve(problem, alg = Euler(), dt = 0.01)
    plot!(sol, alpha=0.3, color = "#BBBBBB", legend = :false)
end
xlabel!("Hour of day")

display(DOPlot)

# Plot prior distributions
using Plots, Distributions, StatsPlots
plot(Normal(GPP_daily_mu, GPP_daily_sigma))
plot(Normal(ER_daily_mu, ER_daily_sigma))
plot(LogNormal(k600_meanlog, k600_sdlog))
plot(truncated(Cauchy(0, err_obs_iid_sigma_scale), 0, Inf))
plot(truncated(Cauchy(0, err_proc_iid_sigma_scale), 0, Inf))

max_stars_repo_name: elizabethjmohr/mesocosm-metabolism

max_issues_repo_name: elizabethjmohr/mesocosm-metabolism

max_forks_repo_name: elizabethjmohr/mesocosm-metabolism


text:
"""
Check that dfs can find AbstractDelayed objects nested inside non-AbstractDelayed objects
"""
module TestDFS

using Test
using TreeParzen
import TreeParzen: Delayed, Graph

a = Dict(
    :c => 11,
    :d => HP.Normal(:test, 1.0, 2.0),
)
order = Graph.dfs(a)
@test all(
    isa.(
        order,
        [
            Delayed.Normal,
            Delayed.AbstractParam,
        ]
    )
)
end
true


text:
using InlineStrings, BenchmarkTools, Random

const types = [String1, String3, String7, String15, String31, String63, String127, String255]

ascii = [InlineString(randstring(i)) for i in map(x -> sizeof(x) - 1, types)]
str_ascii = [String(x) for x in ascii]
utf8 = [InlineString(string("π", x[3:end])) for x in ascii]
str_utf8 = [String(x) for x in utf8]

function bench(ascii, )
    for (x, y) in zip(ascii, str_ascii)
        @show typeof(x), @btime reverse($x)
        @show typeof(y), @btime reverse($y)
    end

    for (x, y) in zip(utf8, str_utf8)
        @show typeof(x), @btime reverse($x)
        @show typeof(y), @btime reverse($y)
    end
end
bench()


text:
function test_mf()
    println("-- Testing MF-based recommender")

    m = [NaN 3 NaN 1 2 1 NaN 4
         1 2 NaN NaN 3 2 NaN 3
         NaN 2 3 3 NaN 5 NaN 1]
    data = DataAccessor(m)

    recommender = MF(data, 2)
    build!(recommender, learning_rate=15e-4, max_iter=100)

    # top-4 recommantion list should be same as CF/SVD-based recommender
    rec = recommend(recommender, 1, 4, [i for i in 1:8])
    @test Set([first(r) for r in rec]) == Set([2, 5, 6, 8])
end

test_mf()


text:
include("../src/global/ImportFolders.jl")

import Config
import Server
import Connection
import Event
import Helper

function onBeforeResponse(sock::Base.TCPSocket, cmd::Connection.Command, ans::Connection.Answer)
  println("before-response\n", cmd)
  ans.data = rand(1:300)
  #for i = 1:1000 Server.request(sock, 1, Helper.Point(rand(1:200), rand(1:200))) end
end

function onAfterResponse(ans::Connection.Answer)
  println("after-serponse\n", ans)
end

con = Server.create(ip"127.0.0.1", Config.val(:conServerPort))
Event.on(con.observer, Connection.EVENT_BEFORE_RESPONSE, onBeforeResponse)
Event.on(con.observer, Connection.EVENT_AFTER_REQUEST, onAfterResponse)
Server.run(con)

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
import Base: cat, hcat, vcat
import AutoGrad: back
using AutoGrad: Arg, Value, forw, uncat
using Base: dims2cat, cat_shape, __cat
using CUDA: CuArray

# Concatenation:

# Need to extend cat definitions from AutoGrad/src/cat.jl:
const NAVK = Union{Number,AbstractArray,Value,KnetArray}
cat(X::NAVK...; dims) = forw(cat,X...;dims=dims)
back(::typeof(cat),::Type{Arg{N}},y1::NAVK,y::NAVK,x::NAVK...; dims) where {N}=uncat(y1,N,dims,x...)

# Benchmarks in μs for hcat and vcat: a=rand(1000,1000) v=rand(1000), t=v'
#		cpu	gpu	g->c->g	vkernel
# hcat(a,a)	2350	225	16160
# hcat(a,v)	1230	115	6490
# hcat(v,a)	1220	120	6490
# hcat(v,v)	3.53	12.53	48.49
# vcat(a,a)	2630	10980	16590	665
# vcat(a,t)	1350	10860	6550	338
# vcat(t,a)	1360	10850	6570	338
# vcat(v,v)	2.13	12.33	45.40	13.58

# setindex! methods called by hcat/vcat:
# hcat(v,v): I = (Colon(),1:1) I = (Colon(),2:2)
# vcat(v,v): uses single index
# hcat(m,m): I = (Colon(),1:5) I = (Colon(),6:10)
# vcat(m,m): I = (1:3,Colon()) I = (4:6,Colon())

# based on _typed_hcat{T}(::Type{T}, A::AbstractVecOrMat...) in base/abstractarray.jl:1316 julia-1.2.0
function hcat(A::KnetVecOrMat{T}...) where {T}
    nargs = length(A)
    nrows = size(A[1], 1)
    ncols = 0
    for j = 1:nargs
        Aj = A[j]
        if size(Aj, 1) != nrows
            throw(ArgumentError("number of rows of each array must match (got $(map(x->size(x,1), A)))"))
        end
        nd = ndims(Aj)
        ncols += (nd==2 ? size(Aj,2) : 1)
    end
    B = similar(A[1], nrows, ncols)
    pos = 1
    for k = 1:nargs
        Ak = A[k]
        n = length(Ak)
        copyto!(B, pos, Ak, 1, n)
        pos += n
    end
    return B
end

function vcat(A::KnetVector{T}...) where {T}
    nargs = length(A)
    nrows = 0
    for a in A
        nrows += length(a)
    end
    B = similar(A[1], nrows)
    pos = 1
    for k = 1:nargs
        Ak = A[k]
        n = length(Ak)
        copyto!(B, pos, Ak, 1, n)
        pos += n
    end
    return B
end

# based on _typed_vcat{T}(::Type{T}, A::AbstractVecOrTuple{AbstractVecOrMat}) in base/abstractarray.jl:1353 julia-1.2.0
function vcat(A::KnetVecOrMat{T}...) where {T}
    nargs = length(A)
    nrows = sum(a->size(a, 1), A)::Int
    ncols = size(A[1], 2)
    for j = 2:nargs
        if size(A[j], 2) != ncols
            throw(ArgumentError("number of columns of each array must match (got $(map(x->size(x,2), A)))"))
        end
    end
    B = similar(A[1], T, nrows, ncols)
    pos = 1
    for k = 1:nargs
        Ak = A[k]
        p1 = pos+size(Ak,1)-1
        B[pos:p1, :] = Ak
        pos = p1+1
    end
    return B
end

function vcat(A::KnetArray{T}...) where {T}
    nargs = length(A)
    size2 = size(A[1])[2:end]
    for j = 2:nargs
        if size(A[j])[2:end] != size2
            throw(ArgumentError("size(a)[2:end] of each array must match (got $(map(size, A)))"))
        end
    end
    A2 = (reshape(a, size(a,1), :) for a in A)
    B = vcat(A2...)
    reshape(B, size(B,1), size2...)
end

function cat_old(a1::KnetVecOrMat{T}, a::KnetVecOrMat{T}...; dims) where {T}
    if     dims==1 || dims==Val(1); vcat(a1, a...)
    elseif dims==2 || dims==Val(2); hcat(a1, a...)
    else error("cat(a...;dims=$dims) not implemented.")
    end
end

# Avoid using Base for unimplemented cat methods:

# using AutoGrad: NA # Union{Number,AbstractArray}
# const NAK = Union{Number,AbstractArray,KnetArray}
# cat(a::NA, as::NA...; dims)=Base._cat(dims, a, as...)
# cat(a::NAK, as::NAK...; dims)=throw(MethodError(cat, (a, as...)))
# hcat(a::NA, as::NA...)=cat(a,as...; dims=2)
# hcat(a::NAK, as::NAK...)=throw(MethodError(hcat, (a, as...)))
# vcat(a::NA, as::NA...)=cat(a,as...; dims=1)
# vcat(a::NAK, as::NAK...)=throw(MethodError(vcat, (a, as...)))

# # Ambiguity fix for abstractarray.jl:1066-1072
# using Base: hvcat_fill, promote_typeof
# vcat(X::Number, Xs::Number...) = hvcat_fill(Array{promote_typeof(X, Xs...)}(undef,1+length(Xs)), (X, Xs...))
# hcat(X::Number, Xs::Number...) = hvcat_fill(Array{promote_typeof(X, Xs...)}(undef,1,1+length(Xs)), (X, Xs...))

# vcat(X::KnetArray...)=cat(X...; dims=Val(1)) # karray.jl version is 30%-80% faster
hcat(X::KnetArray...)=cat(X...; dims=Val(2))   # This should only kick in for dims > 2, karray.jl 1/2D versions are 100% faster

# Based on _cat_t, abstractarray.jl:1439, julia 1.2.0
function cat(X::KnetArray{T}...; dims) where {T}
    catdims = dims2cat(dims)
    # catdims == (true,) && return vcat_old(X...) # 10-30% faster
    shape = cat_shape(catdims, (), map(size, X)...)
    # length(shape) <= 2 && catdims == (false,true) && return hcat_old(X...) # 50% faster
    A = similar(X[1], T, shape) # cat_similar(X[1], T, shape)
    if T <: Number && count(!iszero, catdims) > 1
        fill!(A, zero(T))
    end
    __cat(CuArray(A), shape, catdims, map(CuArray,X)...)
    return A
end


text:
names = readdir(".")

for name=names
    if startswith(name, "test_") && endswith(name, ".jl")
        println("Run $name")
        run(`julia $name`)
    end
end


text:
# Recreate figure 3.3 from the original paper
using Plots
using DRAGO

t = 0:0.05:100
s = -sin.(t) .- cos.(0.5t)
(y, x, w, S) = simulate_data(s, percentage_gaps = 0.45, max_gap_ratio = 0.025)
(fRes, xRes, residual) = drago(y, 0.6, 500, S, ΔFmin = 10e-2)

p1 = plot(s, lab = "Simulated")
p2 = scatter(get_present(S), S * y, lab = "Noise + Gaps (y)")
p3 = scatter(get_present(S), S * y, lab = "", color = :lightgrey, alpha = 0.2)
plot!(fRes, color = :black, linewidth = 2, lab = "Smoothed signal (f)")
p4 = plot(x, color = :red, lab = "True outlier")
plot!(S' * xRes, color = :black, lab = "Estimated outlier (x)")
p5 = plot(residual, lab = "Residual y - (x + f)", color = :black)
plot(p1, p2, p3, p4, p5, ylims = (-4, 4), layout = (5, 1), size = (800, 800))

savefig("./test.svg")