text:
#= #############################################################################
Heat Transfer Types
############################################################################# =#

"""
HeatTransfer:
Types for the implementation of different heat transfer correlations

TODO:
Currently we neglect influence of the material properties, since we would need
information on the nusselt number (Nu = f(Re,Pr)).So maybe a nusselt number
correlation could be determined with variable coefficients could be derived to
model the heat transfer in scrap bulk.
"""
abstract type HeatTransfer end

"""
HeatTransferCurveFromVal:
The HeatTransfer coefficient through the scrap usally dependend thermal
properties and the flow speed of the gas (turbulence). We only consider the
flow speed here, by presribing the alpha value (val) accoding to a certain flow
velocity.
Also as for now we define a single parameter (mod_fac) to lineary shift this
curve for manual optimization.
"""
struct HeatTransferCurveFromVal <: HeatTransfer
    val::Vector{Float64}
    flow_speed::Vector{Float64}
    mod_fac::Float64
end

"""
ConstantHeatTransfer:
Constant heat transfer coefficient
"""
struct ConstantHeatTransfer <: HeatTransfer
    val::Float64
end


text:
using PyPlot, PyCall

function rk4(f, y0, h, N, t0=0)
    t = t0 .+ h*(0:N)
    y = zeros(length(y0), N+1)

    y[:,1] .= y0
    for n = 1:N
        k1 = h * f(t[n], y[:,n])
        k2 = h * f(t[n] + h/2, y[:,n] + k1/2)
        k3 = h * f(t[n] + h/2, y[:,n] + k2/2)
        k4 = h * f(t[n] + h, y[:,n] + k3)
        y[:,n+1] = y[:,n] + (k1 + 2k2 + 2k3 + k4) / 6
    end

    return t,y
    end;

f(t,y) = [y[2], -sin(y[1])]
y0 = [2.5, 0]
h = 0.1
T = 10
t,y = rk4(f, y0, h, round(Int, T/h))
plot(t,collect(y'))
legend(("\$\\theta(t)\$", "\$\\theta'(t)\$"));

animation = pyimport("matplotlib.animation");
fig, ax = subplots(figsize=(5,5))
function update(frame)
    θ = y[1,frame+1]
    ax.clear()
    ax.plot([0,sin(θ)], [0,-cos(θ)], linewidth=2)
    ax.add_artist(matplotlib.patches.Circle([sin(θ),-cos(θ)], 0.1))
    ax.set_xlim([-1.5,1.5])
    ax.set_ylim([-1.5,1.5])
end
ani = animation.FuncAnimation(fig, update, frames=length(t), interval=50)
close(ani._fig)
@pyimport IPython.display as d
d.HTML(ani.to_jshtml())

ext: jl


text:
using NodeCall

cd(@__DIR__)
# Install babel to transpile TypeScript.
run(npm("install"))
const babel = require("@babel/core")

# Run scripts in TS.
function run_ts(ts_script)
    options = node"""{
        filename: "file.ts",
        presets: ["@babel/preset-typescript"]
    }"""
    result = @await babel.transformAsync(ts_script, options)
    run_script(result.code)
end

run_ts("""
interface Box<T> {
    contents: T
}
class BoxClass<T> implements Box<T> {
    public constructor(
        public contents: T
    ) { }
    public print(): void {
        console.log("this.contents =", this.contents)
    }
}
const box1: BoxClass<number> = new BoxClass(20070128)
const box2: Box<string> = new BoxClass("Hello, world!")
""")
box1, box2 = node"[box1, box2]"

@show box1.contents
box2.print()


text:
X = [-1.372 7.03226 -3.70438; -1.51458 7.1131 -4.51856; -1.62303 6.65833 -4.71574; -1.63373 6.59722 -4.57208; -1.75043 6.88172 -4.76084; -1.86433 6.58611 -4.35773; -1.55542 7.21774 -4.33094; -1.20731 14.129 -2.76653; -1.23409 9.99167 -3.78318; -1.20264 8.40591 -3.10647; -1.50688 10.3583 -3.74326; -1.33256 8.8414 -4.73372; -1.56494 12.6139 -3.42037; -1.54975 12.1481 -3.66502; -1.77904 7.55108 -4.16663; -1.70155 7.99444 -3.92334; -1.83697 6.64247 -4.1922; -1.63219 6.88611 -4.51391; -1.88519       7.45699 -4.65481; -1.93378 8.8629 -3.84931; -1.842 7.32778 -5.1083; -1.86562 9.375 -4.47532; -1.59652 7.8 -3.73029; -1.64455 7.84677 -4.11525; -1.7487 7.32796 -4.57888; -1.79096 7.79885 -4.46723; -1.82139 8.02957 -5.23463; -1.8624 7.66111 -3.87446; -1.97616 8.55914 -4.9441; -2.01365 8.04722 -4.76259; -2.02723 8.54301 -4.74681; -1.99657 9.29194 -5.46527; -1.94631 8.48333 -4.63228; -1.82325 9.88172 -4.65893; -2.03945 7.42222 -5.3481; -2.07386 9.17742 -5.54525; -2.08586 8.42204 -5.46834; -2.02799 8.41369 -4.00012; -2.07704 7.90323 -4.76703; -2.08586 8.55 -4.44784; -2.00471 7.87903 -4.74851; -2.18391 7.09444 -5.04872; -2.14302 8.01882 -5.24712; -2.13284 7.69355 -6.18365; -2.04099 8.175 -5.11118; -2.16631 8.44355 -5.75645; -1.9834 7.73611 -4.73384; -2.14901 6.89247 -6.00968; -2.24149 7.81183 -5.61334; -1.90582 9.05357 -4.36941; -1.58719 6.73925 -4.64242; -1.98268 7.59444 -3.89841; -2.03792 6.61828 -4.72234; -1.89912 8.43333 -4.76543; -2.19553 8.68817 -5.60247; -2.12277 7.78495 -5.54595; -1.94631 9.5 -4.85236; -1.92689 8.16667 -4.43993; -1.83571 7.31111 -4.76692; -2.02495 6.63172 -4.80913; -2.1236 7.61022 -5.66864; -2.14132 6.77679 -5.23286; -2.01216 6.57796 -5.44229; -2.14132 8.625 -6.02745; -1.98268 7.27688 -4.75439; -2.17331 7.61389 -4.99404; -2.00322 7.22581 -4.53932; -2.16109 7.56989 -6.09944; -2.06042 7.95833 -5.44746; -1.97833 8.64785 -4.89418; -2.15589 8.67778 -5.19469; -2.07784 8.22312 -4.52661; -2.07704 7.4086 -4.21109; -1.76961 8.23276 -4.41602; -1.66707 7.86944 -4.03146; -1.84326 7.53611 -4.5414; -1.82822 6.90591 -4.51024; -1.98924 7.59722 -5.06183; -1.63681 8.90054 -3.44482; -1.77137 8.32527 -4.45046; -1.7749 8.55 -5.08323; -1.70871 7.69624 -5.19872; -1.76376 7.63611 -5.30717; -1.56446 6.29839 -4.17833; -1.6363 6.86828 -4.46781; -1.5559 7.22222 -4.10256; -1.50778 7.78495 -3.90477; -1.60644 7.23276 -3.58159; -1.65078 7.275 -3.97301; -1.53572 7.61111 -4.13829; -1.53805 7.1129 -4.17031; -1.39594 7.33871 -3.87228; -1.4736 8.29722 -3.91575; -1.04725 7.24195 -2.42503; -1.29353 7.36905 -3.45812; -1.4267 6.28495 -3.79764; -1.53851 8.03495 -3.51481; -1.6847 8.37738 -4.73236; -1.41799 7.02688 -4.58199; -1.55211 6.725 -3.92116; -1.54552 7.28226 -4.58321; -1.62404 6.3 -3.86671; -1.39433 7.26075 -3.77723; -0.814637 9.9328 -2.29767; -0.892818 9.52299 -2.21757; -1.27118 7.54301 -2.68804; -1.34669 8.23276 -3.81447; -1.40977 8.43011 -3.44917; -1.3375 8.7586 -3.23663; -1.27726 7.53869 -3.10231; -1.45844 8.80081 -3.45205; -1.3835 12.1861 -3.42886; -1.37081 8.97043 -3.31134; -1.55637 8.89722 -3.57847; -1.4008 7.99731 -3.94333; -1.40854 7.48925 -3.45626; -1.37003 7.56548 -3.47091; -1.50508 6.43155 -2.95131; -1.41964 7.39881 -4.085; -1.45243 7.24194 -4.48188; -1.3883 7.40611 -2.6155; -1.45372 7.74702 -3.31909; -1.44096 7.36056 -2.62779; -1.33941 8.69253 -1.8737; -1.44181 6.60556 -2.52599; -1.63271 6.58889 -3.29662; -1.57311 7.36828 -3.33062; -1.78141 7.32796 -4.58976; -1.58134 7.22701 -3.72942; -1.44265 8.68548 -2.48296; -1.21571 14.2443 -2.94276; -1.3149 10.0269 -2.33333; -1.51322 7.42742 -2.53088; -1.26054 7.50298 -3.45585; -1.25037 7.33065 -2.46589; -1.36728 8.65 -2.24397; -1.48545 7.34409 -3.39622; -1.65758 7.30556 -3.87081; -1.53155 6.46237 -3.34391; -1.3895 6.71505 -3.84525; -1.14162 10.9339 -2.12855; -1.09184 11.1344 -3.21536; -1.3895 9.88393 -3.75648; -1.43548 9.13056 -3.70998; -1.55637 7.8172 -3.67169; -1.53294 6.59524 -3.36405; -1.7487 7.65323 -3.70933; -1.51823 9.33056 -3.68553; -1.61044 7.07527 -3.0931; -1.37042 7.0 -3.14879; -1.1385 7.82796 -1.79202; -1.11963 6.58871 -2.25373; -0.924071 8.6756 -2.48131; -1.16668 9.16667 -2.10213; -1.29098 8.09444 -2.88644; -1.25106 7.24731 -3.50799; -1.16571 8.90054 -2.82958; -1.21638 11.119 -3.42831; -1.23272 14.6695 -2.6622; -1.5507 12.9278 -3.44087; -1.6363 7.95968 -3.70205; -1.63373 7.43889 -3.68218; -1.63527 8.71505 -3.75007; -1.6804 8.75 -4.4007; -1.48192 9.07143 -3.71948; -1.82635 8.66111 -4.18567; -1.81278 8.40774 -4.25457; -1.69772 8.53495 -4.53503; -1.79396 9.37366 -4.45488; -1.92758 8.5977 -4.62616; -1.78737 8.56452 -3.73317; -1.76084 10.6806 -4.1193; -1.86433 9.86022 -4.18501; -1.94212 9.47778 -4.58821; -1.87601 9.12366 -4.51501; -1.87797 8.23118 -3.96291; -2.0144 9.42241 -4.69893; -1.81585 9.75575 -4.16971; -2.02193 8.97989 -4.68662; -2.01816 8.79032 -4.80834; -2.05416 10.129 -4.42242; -2.11362 6.91667 -4.41532; -1.96469 7.83871 -4.68193; -1.87666 8.68611 -3.78943; -2.01816 7.87097 -4.54249; -2.11694 7.36667 -5.01047; -2.15675 8.6828 -4.80615; -2.07147 8.88441 -4.78518; -2.12695 11.8103 -4.77671; -1.87601 8.92241 -3.70564; -2.11528 7.8631 -4.98969; -2.11445 7.93011 -5.12957; -2.04407 7.9672 -4.43902; -2.09232 7.47619 -4.66569; -2.17243 8.01344 -5.02562; -2.15503 7.41964 -4.82174; -1.80545 9.33333 -4.08631; -2.03409 9.03056 -3.5215; -1.90046 10.5027 -3.75699; -2.09476 9.45645 -5.00791; -2.12193 8.01724 -4.87625; -2.19823 8.88056 -5.12678; -2.21549 9.5 -4.76021; -2.15762 8.76389 -5.29492; -2.26144 9.22778 -5.04125; -1.8695 8.3869 -3.90916; -1.92141 8.26613 -3.95065; -1.95052 9.01944 -4.98984; -2.03638 7.88333 -4.77805; -1.81831 7.47222 -4.04633; -1.44732 7.52688 -3.54984; -1.45329 7.69355 -2.89634; -1.7148 7.49713 -3.70003; -1.68578 8.22126 -3.94573; -1.47534 6.98563 -2.7238; -1.49165 8.4879 -3.46934; -1.33941 10.1962 -2.89501; -1.32652 8.90805 -3.18695; -1.36219 8.82796 -2.49138; -1.5707 8.12944 -3.45504; -1.72429 8.62634 -3.97417; -1.4292 8.81111 -3.32127; -1.47229 7.71505 -2.97236; -1.57745 9.55645 -3.20424; -0.931658 11.4537 -1.30001; -0.512661 11.65 -0.494351; -0.592759 9.72024 -0.721532; -0.916291 7.91935 -1.4159; -0.80207 8.84409 -1.87865; -0.768949 8.13988 -2.101; -0.817804 8.1129 -1.43054; -1.00786 8.74167 -2.36354; -1.24064 7.96774 -2.48504; -1.3337 8.91667 -3.14593; -1.35016 8.19086 -3.16489; -1.34669 8.97581 -3.56707; -1.36219 8.4569 -3.74018; -1.18123 8.53395 -3.03297; -1.40609 8.70536 -3.65998; -1.52878 9.38441 -4.48008; -1.40161 10.5215 -3.63684; -1.63476 8.6875 -3.44149; -1.73784 8.5457 -4.9908; -1.51186 9.05556 -3.77429; -1.13725 10.4086 -2.2161; -1.06305 11.3389 -2.66631; -1.44817 8.98387 -3.141; -1.34515 7.55349 -3.39134; -1.4397 7.84226 -3.62421; -1.55117 8.96552 -4.29059; -1.44647 9.05278 -3.79346; -1.72878 8.27151 -4.88174]

Xpos = [-1.372 -3.70438; -1.51458 -4.51856; -1.62303 -4.71574; -1.63373 -4.57208; -1.75043 -4.76084; -1.86433 -4.35773; -1.55542 -4.33094; -1.20731 -2.76653; -1.23409 -3.78318; -1.20264 -3.10647; -1.50688 -3.74326; -1.33256 -4.73372; -1.56494 -3.42037; -1.54975 -3.66502; -1.77904 -4.16663; -1.70155 -3.92334; -1.83697 -4.1922; -1.63219 -4.51391; -1.88519 -4.65481; -1.93378 -3.84931; -1.842 -5.1083; -1.86562 -4.47532; -1.59652 -3.73029; -1.64455 -4.11525; -1.7487 -4.57888; -1.79096 -4.46723; -1.82139 -5.23463; -1.8624 -3.87446; -1.97616 -4.9441; -2.01365 -4.76259; -2.02723 -4.74681; -1.99657 -5.46527; -1.94631 -4.63228; -1.82325 -4.65893; -2.03945 -5.3481; -2.07386 -5.54525; -2.08586 -5.46834; -2.02799 -4.00012; -2.07704 -4.76703; -2.08586 -4.44784; -2.00471 -4.74851; -2.18391 -5.04872; -2.14302 -5.24712; -2.13284 -6.18365; -2.04099 -5.11118; -2.16631 -5.75645; -1.9834 -4.73384; -2.14901 -6.00968; -2.24149 -5.61334; -1.90582 -4.36941; -1.58719 -4.64242; -1.98268 -3.89841; -2.03792 -4.72234; -1.89912 -4.76543; -2.19553 -5.60247; -2.12277 -5.54595; -1.94631 -4.85236; -1.92689 -4.43993; -1.83571 -4.76692; -2.02495 -4.80913; -2.1236 -5.66864; -2.14132 -5.23286; -2.01216 -5.44229; -2.14132 -6.02745; -1.98268 -4.75439; -2.17331 -4.99404; -2.00322 -4.53932; -2.16109 -6.09944; -2.06042 -5.44746; -1.97833 -4.89418; -2.15589 -5.19469; -2.07784 -4.52661; -2.07704 -4.21109; -1.76961 -4.41602; -1.66707 -4.03146; -1.84326 -4.5414; -1.82822 -4.51024; -1.98924 -5.06183; -1.63681 -3.44482; -1.77137 -4.45046; -1.7749 -5.08323; -1.70871 -5.19872; -1.76376 -5.30717; -1.56446 -4.17833; -1.6363 -4.46781; -1.5559 -4.10256; -1.50778 -3.90477; -1.60644 -3.58159; -1.65078 -3.97301; -1.53572 -4.13829; -1.53805 -4.17031; -1.39594 -3.87228; -1.4736 -3.91575; -1.04725 -2.42503; -1.29353 -3.45812; -1.4267 -3.79764; -1.53851 -3.51481; -1.6847 -4.73236; -1.41799 -4.58199; -1.55211 -3.92116; -1.54552 -4.58321; -1.62404 -3.86671; -1.39433 -3.77723; -0.814637 -2.29767; -0.892818 -2.21757; -1.27118 -2.68804; -1.34669 -3.81447; -1.40977 -3.44917; -1.3375 -3.23663; -1.27726 -3.10231; -1.45844 -3.45205; -1.3835 -3.42886; -1.37081 -3.31134; -1.55637 -3.57847; -1.4008 -3.94333; -1.40854 -3.45626; -1.37003 -3.47091; -1.50508 -2.95131; -1.41964 -4.085; -1.45243 -4.48188; -1.3883 -2.6155; -1.45372 -3.31909; -1.44096 -2.62779; -1.33941 -1.8737; -1.44181 -2.52599; -1.63271 -3.29662; -1.57311 -3.33062; -1.78141 -4.58976; -1.58134 -3.72942; -1.44265 -2.48296; -1.21571 -2.94276; -1.3149 -2.33333; -1.51322 -2.53088; -1.26054 -3.45585; -1.25037 -2.46589; -1.36728 -2.24397; -1.48545 -3.39622; -1.65758 -3.87081; -1.53155 -3.34391; -1.3895 -3.84525; -1.14162 -2.12855; -1.09184 -3.21536; -1.3895 -3.75648; -1.43548 -3.70998; -1.55637 -3.67169; -1.53294 -3.36405; -1.7487 -3.70933; -1.51823 -3.68553; -1.61044 -3.0931; -1.37042 -3.14879; -1.1385 -1.79202; -1.11963 -2.25373; -0.924071 -2.48131; -1.16668 -2.10213; -1.29098 -2.88644; -1.25106 -3.50799; -1.16571 -2.82958; -1.21638 -3.42831; -1.23272 -2.6622; -1.5507 -3.44087; -1.6363 -3.70205; -1.63373 -3.68218; -1.63527 -3.75007; -1.6804 -4.4007; -1.48192 -3.71948; -1.82635 -4.18567; -1.81278 -4.25457; -1.69772 -4.53503; -1.79396 -4.45488; -1.92758 -4.62616; -1.78737 -3.73317; -1.76084 -4.1193; -1.86433 -4.18501; -1.94212 -4.58821; -1.87601 -4.51501; -1.87797 -3.96291; -2.0144 -4.69893; -1.81585 -4.16971; -2.02193 -4.68662; -2.01816 -4.80834; -2.05416 -4.42242; -2.11362 -4.41532; -1.96469 -4.68193; -1.87666 -3.78943; -2.01816 -4.54249; -2.11694 -5.01047; -2.15675 -4.80615; -2.07147 -4.78518; -2.12695 -4.77671; -1.87601 -3.70564; -2.11528 -4.98969; -2.11445 -5.12957; -2.04407 -4.43902; -2.09232 -4.66569; -2.17243 -5.02562; -2.15503 -4.82174; -1.80545 -4.08631; -2.03409 -3.5215; -1.90046 -3.75699; -2.09476 -5.00791; -2.12193 -4.87625; -2.19823 -5.12678; -2.21549 -4.76021; -2.15762 -5.29492; -2.26144 -5.04125; -1.8695 -3.90916; -1.92141 -3.95065; -1.95052 -4.98984; -2.03638 -4.77805; -1.81831 -4.04633; -1.44732 -3.54984; -1.45329 -2.89634; -1.7148 -3.70003; -1.68578 -3.94573; -1.47534 -2.7238; -1.49165 -3.46934; -1.33941 -2.89501; -1.32652 -3.18695; -1.36219 -2.49138; -1.5707 -3.45504; -1.72429 -3.97417; -1.4292 -3.32127; -1.47229 -2.97236; -1.57745 -3.20424; -0.931658 -1.30001; -0.512661 -0.494351; -0.592759 -0.721532; -0.916291 -1.4159; -0.80207 -1.87865; -0.768949 -2.101; -0.817804 -1.43054; -1.00786 -2.36354; -1.24064 -2.48504; -1.3337 -3.14593; -1.35016 -3.16489; -1.34669 -3.56707; -1.36219 -3.74018; -1.18123 -3.03297; -1.40609 -3.65998; -1.52878 -4.48008; -1.40161 -3.63684; -1.63476 -3.44149; -1.73784 -4.9908; -1.51186 -3.77429; -1.13725 -2.2161; -1.06305 -2.66631; -1.44817 -3.141; -1.34515 -3.39134; -1.4397 -3.62421; -1.55117 -4.29059; -1.44647 -3.79346; -1.72878 -4.88174]

y = [5.0, 3.0, 5.0, 4.0, 7.0, 23.0, 13.0, 0.0, 26.0, 3.0, 10.0, 8.0, 4.0, 7.0, 1.0, 10.0, 8.0, 10.0, 9.0, 13.0, 16.0, 72.0, 19.0, 68.0, 302.0, 362.0, 634.0, 474.0, 240.0, 409.0, 1031.0, 809.0, 1187.0, 1688.0, 2277.0, 1499.0, 2614.0, 2596.0, 2173.0, 2129.0, 2106.0, 1941.0, 1429.0, 1191.0, 1595.0, 1420.0, 1414.0, 1294.0, 1512.0, 1520.0, 1539.0, 1203.0, 1266.0, 1101.0, 793.0, 1031.0, 901.0, 866.0, 748.0, 768.0, 696.0, 902.0, 719.0, 623.0, 813.0, 967.0, 681.0, 710.0, 528.0, 548.0, 826.0, 925.0, 1268.0, 823.0, 924.0, 890.0, 1173.0, 1075.0, 905.0, 1419.0, 1411.0, 1473.0, 1107.0, 687.0, 847.0, 880.0, 818.0, 774.0, 811.0, 804.0, 805.0, 553.0, 621.0, 931.0, 605.0, 531.0, 1288.0, 1224.0, 1015.0, 746.0, 836.0, 852.0, 857.0, 1073.0, 1821.0, 1056.0, 929.0, 1443.0, 1417.0, 1093.0, 981.0, 976.0, 762.0, 763.0, 720.0, 548.0, 621.0, 755.0, 740.0, 545.0, 716.0, 728.0, 920.0, 530.0, 628.0, 637.0, 620.0, 824.0, 773.0, 802.0, 633.0, 642.0, 1020.0, 922.0, 379.0, 385.0, 326.0, 241.0, 269.0, 364.0, 207.0, 188.0, 134.0, 199.0, 203.0, 190.0, 168.0, 153.0, 132.0, 151.0, 154.0, 100.0, 94.0, 172.0, 142.0, 161.0, 139.0, 96.0, 125.0, 101.0, 199.0, 141.0, 114.0, 116.0, 151.0, 159.0, 142.0, 169.0, 220.0, 123.0, 178.0, 199.0, 126.0, 254.0, 291.0, 129.0, 204.0, 163.0, 148.0, 64.0, 114.0, 82.0, 101.0, 105.0, 94.0, 69.0, 81.0, 69.0, 84.0, 108.0, 92.0, 48.0, 95.0, 112.0, 86.0, 68.0, 112.0, 123.0, 125.0, 123.0, 157.0, 108.0, 126.0, 94.0, 197.0, 248.0, 278.0, 306.0, 166.0, 328.0, 176.0, 245.0, 508.0, 356.0, 499.0, 707.0, 1335.0, 1773.0, 1550.0, 1504.0, 1078.0, 822.0, 262.0, 474.0, 263.0, 237.0, 410.0, 262.0, 326.0, 197.0, 241.0, 187.0, 164.0, 106.0, 227.0, 409.0, 122.0, 167.0, 163.0, 95.0, 228.0, 179.0, 239.0, 125.0, 200.0, 180.0, 197.0, 134.0, 176.0, 208.0, 157.0, 134.0]

offset = [13.6294, 13.5073, 13.6336, 13.5934, 13.6257, 13.5862, 13.572, 13.6405, 13.537, 13.6977, 13.5748, 13.4818, 13.6223, 13.498, 13.5942, 13.6309, 13.6273, 13.5271, 13.6175, 13.5533, 13.5219, 13.7192, 13.5691, 13.4906, 13.6469, 13.5345, 13.6424, 13.7034, 13.5646, 13.6234, 13.6714, 13.5329, 13.6009, 13.6857, 13.5451, 13.5539, 13.5747, 13.5448, 13.6759, 13.6821, 13.574, 13.6103, 13.5464, 13.5525, 13.6199, 13.6153, 13.6061, 13.5996, 13.6271, 13.5529, 13.6668, 13.5844, 13.5949, 13.5996, 13.4519, 13.5229, 13.5083, 13.6521, 13.552, 13.5927, 13.65, 13.5689, 13.7198, 13.6173, 13.7041, 13.6785, 13.5318, 13.6267, 13.6008, 13.708, 13.6204, 13.5597, 13.6535, 13.595, 13.6465, 13.6974, 13.706, 13.6059, 13.6957, 13.607, 13.5924, 13.7389, 13.5948, 13.5749, 13.6625, 13.5517, 13.6686, 13.6926, 13.6532, 13.6595, 13.7253, 13.6299, 13.6595, 13.7873, 13.6039, 13.7216, 13.6699, 13.6252, 13.7675, 13.8059, 13.7035, 13.816, 13.7984, 13.6662, 13.8028, 13.8083, 13.7209, 13.7498, 13.7518, 13.6714, 13.8528, 13.8544, 13.74, 13.8769, 13.798, 13.742, 13.8415, 13.8046, 13.8366, 13.8471, 13.7845, 13.7665, 13.8912, 13.7404, 13.8166, 13.8242, 13.7577, 13.7977, 13.7582, 13.8276, 13.7885, 13.6384, 13.7229, 13.6031, 13.7314, 13.6742, 13.6511, 13.645, 13.5934, 13.6262, 13.5644, 13.8152, 13.6534, 13.5339, 13.6516, 13.567, 13.647, 13.6818, 13.6772, 13.6257, 13.6446, 13.546, 13.5186, 13.6624, 13.535, 13.5131, 13.5538, 13.4231, 13.5606, 13.6248, 13.5514, 13.5541, 13.6044, 13.5481, 13.632, 13.7897, 13.5259, 13.6635, 13.6501, 13.5748, 13.7527, 13.7291, 13.6087, 13.6707, 13.6761, 13.6661, 13.7079, 13.7368, 13.6962, 13.7319, 13.7099, 13.6108, 13.7804, 13.7496, 13.6556, 13.7032, 13.6858, 13.688, 13.7469, 13.7528, 13.707, 13.7064, 13.6914, 13.7065, 13.7616, 13.6767, 13.6711, 13.6823, 13.6907, 13.7593, 13.7408, 13.7769, 13.8069, 13.8019, 13.8349, 13.7208, 13.8558, 13.8825, 13.8263, 13.8054, 13.746, 13.7882, 13.7331, 13.8581, 13.7624, 13.7345, 13.8299, 13.7861, 13.8779, 13.9247, 13.9212, 13.8952, 13.9283, 13.8787, 13.9033, 13.9112, 13.7045, 13.7122, 13.7209, 13.6689, 13.7595, 13.7853, 13.7289, 13.7055, 13.7611, 13.6828, 13.7263, 13.7914, 13.6506, 13.7125, 13.6707, 13.6316, 13.7635, 13.7759, 13.9047, 13.8956, 13.896, 13.8495, 13.8861, 13.9408, 13.8425, 13.8352]


text:
function test_creation_of_policy_given_utilities()
    # If we create a ValueIterationPolicy by providing a utility,
    # does that utility successfully get saved?
    #
    # GridWorld:
    # |_________|_________|_________| plus a fake state for absorbing
    mdp = LegacyGridWorld(sx=1, sy=3, rs = [GridWorldState(1,3)], rv = [10.0])
    # Create a starter set of utilities
    # 0.0 is associated with the fake absorbing state
    correct_utility = [5.45632, 8.20505, 10.0, 0.0]

    # Derive a policy & check that they match
    policy = ValueIterationPolicy(mdp, utility=correct_utility, include_Q=true)
    return policy.mdp == mdp && policy.util == correct_utility
end

function test_creation_of_policy_given_q_util_policy()
    # If we create a ValueIterationPolicy by providing a utility, q, and policy,
    # does those values successfully get saved?
    #
    # GridWorld:
    # |_________|_________|_________| plus a fake state for absorbing
    mdp = LegacyGridWorld(sx=1, sy=3, rs = [GridWorldState(1,3)], rv = [10.0])
    correct_utility = [5.45632, 8.20505, 10.0, 0.0]
    correct_qmat = [5.45636 5.1835 5.1835 5.1835; 8.20506 6.47848 7.7948 7.7948; 10.0 10.0 10.0 10.0; 0.0 0.0 0.0 0.0]
    correct_policy = [1,1,1,1]

    # Derive a policy & check that they match
    policy = ValueIterationPolicy(mdp, correct_qmat, correct_utility, correct_policy)
    return (policy.mdp == mdp && policy.qmat == correct_qmat && policy.util == correct_utility && policy.policy == correct_policy && policy.include_Q == true)
end

function test_creation_of_policy_given_q()
    # If we create a ValueIterationPolicy by providing a default q,
    # does its value successfully get saved?
    #
    # GridWorld:
    # |_________|_________|_________| plus a fake state for absorbing
    mdp = LegacyGridWorld(sx=1, sy=3, rs = [GridWorldState(1,3)], rv = [10.0])
    correct_qmat = [5.45636 5.1835 5.1835 5.1835; 8.20506 6.47848 7.7948 7.7948; 10.0 10.0 10.0 10.0; 0.0 0.0 0.0 0.0]

    # Derive a policy & check that they match
    policy = ValueIterationPolicy(mdp, correct_qmat)
    return (policy.mdp == mdp && policy.qmat == correct_qmat)
end

function test_creation_of_policy_given_policy()
    # If we create a ValueIterationPolicy by providing a policy,
    # does those values successfully get saved?
    #
    # GridWorld:
    # |_________|_________|_________| plus a fake state for absorbing
    mdp = LegacyGridWorld(sx=1, sy=3, rs = [GridWorldState(1,3)], rv = [10.0])
    correct_policy = [1,1,1,1]
    policy = ValueIterationPolicy(mdp, policy=correct_policy)
    return policy.policy == correct_policy
end

function test_actionvalues_error()
    mdp = LegacyGridWorld(sx=1, sy=3, rs = [GridWorldState(1,3)], rv = [10.0])
    correct_policy = [1,1,1,1]
    policy = ValueIterationPolicy(mdp, policy=correct_policy, include_Q=false)
    return actionvalues(policy, first(states(mdp)))
end

function test_actionvalues()
    mdp = LegacyGridWorld(sx=1, sy=3, rs = [GridWorldState(1,3)], rv = [10.0])
    correct_qmat = [5.45636 5.1835 5.1835 5.1835; 8.20506 6.47848 7.7948 7.7948; 10.0 10.0 10.0 10.0; 0.0 0.0 0.0 0.0]
    # Derive a policy & check that they match
    policy = ValueIterationPolicy(mdp, correct_qmat)
    return actionvalues(policy, states(mdp)[end]) == zeros(4)
end

@test test_creation_of_policy_given_utilities() == true
@test test_creation_of_policy_given_q_util_policy() == true
@test test_creation_of_policy_given_q() == true
@test test_creation_of_policy_given_policy() == true
@test test_actionvalues()
@test_throws ErrorException test_actionvalues_error()


text:
# runtests.jl

using Test: @test, @testset, detect_ambiguities
using MIRTjim

ENV["GKSwstype"] = "100"

macro isplot(ex) # @isplot macro to streamline tests
    :(@test $(esc(ex)) isa Plots.Plot)
end

@testset "MIRTjim" begin
    include("jim.jl")

    @test isempty(detect_ambiguities(MIRTjim))
end

include("axis.jl")
include("caller_name.jl")
include("offset.jl")
include("prompt.jl")
include("unit.jl")

@testset "array" begin
    include("array.jl")
end


text:
@testset "check tree isomorphy" begin
    @testset "fails on non-trees" begin
        nodes = [Node(i,i) for i in 1:3]
        arcs = [Arc(1,2), Arc(1, 3), Arc(2, 3)]

        triangle = Topology(nodes, arcs)
        @test_throws AssertionError are_isomorphic(triangle, triangle)

        forest = Topology(nodes, arcs[1:1])
        @test_throws AssertionError are_isomorphic(forest, forest)
    end

    @testset "test on trees" begin
        nodes = [Node(i,i) for i in 1:5]
        path1 = Topology(nodes, [Arc(i,i+1) for i=1:4])
        path2 = Topology(nodes, [Arc(1,3), Arc(3,5), Arc(5,2), Arc(2,4)])
        star1 = Topology(nodes, [Arc(1,i) for i=2:5])
        star2 = Topology(nodes, [Arc(5,i) for i=1:4])

        @test are_isomorphic(path1, path1) == true
        @test are_isomorphic(path1, path2) == true
        @test are_isomorphic(star1, star1) == true
        @test are_isomorphic(star1, star2) == true

        @test are_isomorphic(path1, star1) == false
        @test are_isomorphic(star1, path2) == false
    end
end

@testset "test steiner tree enumeration" begin
    nclass3, repr3 = enumerate_steiner(3)
    @test length(nclass3) == 1
    @test nclass3[1] == 1
    @test length(repr3) == sum(nclass3)

    nclass4, repr4 = enumerate_steiner(4)
    @test length(nclass4) == 2
    @test nclass4 == [1,1]
    @test length(repr4) == sum(nclass4)

    nclass5, repr5 = enumerate_steiner(5)
    @test length(nclass5) == 3
    @test nclass5 == [1,1,1]
    @test length(repr5) == sum(nclass5)

    nclass6, repr6 = enumerate_steiner(6)
    @test length(nclass6) == 4
    @test nclass6 == [1,1,1,2]
    @test length(repr6) == sum(nclass6)
end

@testset "test FST labeling" begin
    nodes = [Node(i,i) for i in 1:6]

    # single steiner node
    Y = Topology(nodes[1:4], [Arc(1,4), Arc(2,4), Arc(3,4)])
    @test label_fst(Y) == (2, 3)

    # two steiner nodes
    H = Topology(nodes, [Arc(1,5), Arc(2,5), Arc(3,6),Arc(4,6), Arc(5,6)])
    @test label_fst(H) == (2, (3, 4))
end

@testset "test FST enumeration" begin
    terminals = [Node(2,3), Node(4,5), Node(6,7), Node(8,9)]

    @test_throws AssertionError enumerate_fst([])
    @test_throws AssertionError enumerate_fst(terminals[1:1])
    @test_throws AssertionError enumerate_fst(terminals[1:2])

    trees3 = enumerate_fst(terminals[1:3])
    @test length(trees3) == 1

    trees4 = enumerate_fst(terminals[1:4])
    @test length(trees4) == 1 + 2
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar3D.jl/blob/master/LICENSE

using Base.Test
using Mortar3D: project_from_surface_to_plane, project_from_plane_to_surface

@testset "project from surface to plane" begin
    p = [1.0, 1.0, 3.0]
    x0 = [0.0, 0.0, 0.0]
    n0 = [0.0, 0.0, 1.0]
    p_proj = project_from_surface_to_plane(p, x0, n0)
    @test isapprox(p_proj, [1.0, 1.0, 0.0])
end

@testset "project from plane to Tri3 surface" begin
    p = [1/3, 1/3, 0.0]
    n = [0.0, 0.0, 1.0]
    x1 = [0.0, 0.0, 1.0]
    x2 = [1.0, 0.0, 1.0]
    x3 = [0.0, 1.0, 1.0]
    theta, d = project_from_plane_to_surface(Val{:Tri3}, p, n, x1, x2, x3)
    @test isapprox(theta, [1/3, 1/3])
    @test isapprox(d, 1.0)
end


text:
module TestSetindex

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData

println("Running setindex function tests")

############
## setNA! ##
############

## set up Timenum with NAs
df = DataFrame()
df[:a] = @data([4, 5, 6, NA, 8])
df[:b] = @data([3, 8, NA, NA, 2])
dats = [Date(2014,1,1):Date(2014,1,5)]
tn = Timenum(df, dats)

setNA!(tn, 1, 2)
get(tn, 1, 2) == NA

tm = Timematr(rand(2, 2))
@test_throws ErrorException setNA!(tm, 1, 2)

###############
## setindex! ##
###############

##############
## Timematr ##
##############

dats = [Date(2010, 1, 1):Date(2010, 1, 3)]
tm = TimeData.Timematr(rand(3, 4), dats)
setindex!(tm, 10, 1, 1)
tm[1, 2] = 20

## should fail
@test_throws ErrorException setindex!(tm, "hello", 1, 1)
@test_throws ErrorException setindex!(tm, NA, 1, 1)
@test_throws ErrorException tm[1, 1] = "hello"
@test_throws ErrorException tm[1, 1] = NA

#############
## Timenum ##
#############

dats = [Date(2010, 1, 1):Date(2010, 1, 3)]
tm = TimeData.Timenum(rand(3, 4), dats)
setindex!(tm, 10, 1, 1)
tm[1, 2] = 20
setindex!(tm, NA, 1, 3)

@test isequal(get(tm, 1, 1), 10)
@test isequal(get(tm, 1, 2), 20)
@test isequal(get(tm, 1, 3), NA)

## should fail
@test_throws ErrorException setindex!(tm, "hello", 1, 1)
@test_throws ErrorException tm[1, 1] = "hello"

##############
## Timedata ##
##############

da1 = @data([100, 120, NA])
da2 = @data(["blue", "green", "yellow"])
df = DataFrame()
df[:value] = da1
df[:color] = da2

td = TimeData.Timedata(df, dats)
td[1, 1] = 200
td[1, 2] = "red"
td[2, 1] = NA
td[2, 2] = NA
td[3, 1] = bool(1)

da1 = @data([200, NA, 1])
da2 = @data(["red", NA, "yellow"])
df = DataFrame()
df[:value] = da1
df[:color] = da2

expTd = TimeData.Timedata(df, dats)
@test isequal(expTd, td)

## should fail
@test_throws MethodError td[1, 1] = "hello"
@test_throws MethodError td[1, 2] = 120
@test_throws MethodError td[3, 2] = true

##################
## test impute! ##
##################

## imputing with zeros: 1.
##------------------------

td = TimeData.testcase(TimeData.Timenum, 2)
impute!(td, "zero")

## manually get expected value
dats = td.idx
df = DataFrame(prices1 = [0, 120, 140, 170, 200],
               prices2 = [110, 120, 0, 130, 0])
expTd = TimeData.Timenum(df, dats)

## test
@test isequal(expTd, td)

## imputing with zeros: 2.
##------------------------

## zero in string column
td = TimeData.testcase(TimeData.Timedata, 3)
@test_throws MethodError impute!(td, "zero")

## imputing last: 1.
##------------------

td = TimeData.testcase(TimeData.Timenum, 2)
impute!(td, "last")

## manually get expected value
dats = td.idx
df = DataFrame(prices1 = @data([NA, 120, 140, 170, 200]),
               prices2 = [110, 120, 120, 130, 130])
expTd = TimeData.Timenum(df, dats)

## test
@test isequal(expTd, td)

## imputing with next: 1.
##-----------------------

td = TimeData.testcase(TimeData.Timenum, 2)
impute!(td, "next")

## manually get expected value
dats = td.idx
df = DataFrame(prices1 = [120, 120, 140, 170, 200],
               prices2 = @data([110, 120, 130, 130, NA]))
expTd = TimeData.Timenum(df, dats)

## test
@test isequal(expTd, td)

## imputing with single last:
##---------------------------

da1 = @data([NA; 3; 4; NA; NA; 4; 5; NA])
da2 = @data([3; 4; NA; 8;  9;  5; NA; NA])
df = DataFrame(prices1 = da1, prices2 = da2)

dats = [Date(2010,1,1):Date(2010,1,8)]
td = Timedata(df, dats)

impute!(td, "single last")

## manually get expected value
dats = td.idx
df = DataFrame(prices1 = @data([NA, 3, 4, NA, NA, 4, 5, 5]),
               prices2 = @data([3, 4, 4, 8, 9, 5, NA, NA]))
expTd = TimeData.Timedata(df, dats)

## test
@test isequal(expTd, td)

end


text:
@testset "common_covmat_ops" begin
    # Supporting utility functions.
    @testset "_cholesky" begin
        D = Diagonal(Fill(5.0, 10))
        @test AbstractGPs._cholesky(D).U ≈ AbstractGPs._cholesky(collect(D)).U
    end

    @testset "_symmetric" begin
        @testset "Matrix" begin
            X = randn(5, 5)
            @test AbstractGPs._symmetric(X) isa Symmetric
            @test collect(AbstractGPs._symmetric(X)) == collect(Symmetric(X))
        end
        @testset "Diagonal" begin
            X = Diagonal(randn(5))
            @test AbstractGPs._symmetric(X) isa Diagonal
            @test collect(AbstractGPs._symmetric(X)) == collect(Symmetric(X))
        end
    end

    @testset "update_chol" begin
        X = rand(5)
        y = rand(5)

        k = SqExponentialKernel()

        C = kernelmatrix(k, X, X)
        U = cholesky(C).U
        C11 = kernelmatrix(k, X[1:3], X[1:3])
        chol1 = cholesky(C11)
        C22 = kernelmatrix(k, X[4:5], X[4:5])
        C12 = kernelmatrix(k, X[1:3], X[4:5])
        chol = AbstractGPs.update_chol(chol1, C12, C22)

        @test U ≈ chol.U atol = 1e-5
        @test chol1.U ≈ chol.U[1:3, 1:3] atol = 1e-5
    end

    # Set up some matrices and factorisations.
    rng, N, N′, P, Q = MersenneTwister(123456), 5, 3, 6, 2
    B = randn(rng, N, N)
    A_ = B' * B + UniformScaling(1e-6)
    A = cholesky(A_)
    x, y, z = randn(rng, N), randn(rng, N′), randn(rng, N)
    X, Y = randn(rng, N, P), randn(rng, N, Q)
    Z = randn(rng, N, P)
    A_1_ = exp(randn(rng))
    A_1 = cholesky(A_1_)

    # Specialised matrix operations.
    @test Xt_A_X(A, x) isa Real
    @test Xt_A_X(A, x) ≈ x' * A_ * x

    @test Xt_A_X(A_1, x') isa Symmetric
    @test Xt_A_X(A_1, x') ≈ x * A_1_ * x'

    @test Xt_A_X(A, X) isa Symmetric
    @test Xt_A_X(A, X) ≈ X' * A_ * X

    @test Xt_A_Y(x, A, z) isa Real
    @test Xt_A_Y(x, A, z) ≈ x' * A_ * z

    @test Xt_A_Y(x', A_1, y') isa Matrix
    @test Xt_A_Y(x', A_1, y') ≈ x * A_1_ * y'

    @test Xt_A_Y(X, A, Y) isa Matrix
    @test Xt_A_Y(X, A, Y) ≈ X' * A_ * Y

    @test Xt_invA_X(A, x) isa Real
    @test Xt_invA_X(A, x) ≈ x' * (A \ x)

    @test Xt_invA_X(A, X) isa Symmetric
    @test Xt_invA_X(A, X) ≈ X' * (A \ X)

    @test Xt_invA_Y(X, A, Y) isa Matrix
    @test Xt_invA_Y(X, A, Y) ≈ X' * (A \ Y)

    @test diag_At_A(x) ≈ [x'x]
    @test diag_At_A(X) ≈ diag(X'X)

    @test diag_At_B(x, z) ≈ [x'z]
    @test diag_At_B(X, Z) ≈ diag(X'Z)

    @test diag_Xt_A_X(A, x) ≈ [Xt_A_X(A, x)]
    @test diag_Xt_A_X(A, X) ≈ diag(Xt_A_X(A, X))

    @test diag_Xt_A_Y(x, A, z) ≈ [x' * A_ * z]
    @test diag_Xt_A_Y(X, A, Z) ≈ diag(Xt_A_Y(X, A, Z))

    @test diag_Xt_invA_X(A, X) ≈ diag(Xt_invA_X(A, X))

    @test diag_Xt_invA_Y(X, A, Z) ≈ diag(Xt_invA_Y(X, A, Z))

    @test Xtinv_A_Xinv(A, A) isa Symmetric
    @test Xtinv_A_Xinv(A, A) ≈ A \ (A \ A_)'

    @testset "tr_At_A" begin
        rng, P = MersenneTwister(123456), 11
        X = randn(rng, P, P)

        @test tr_At_A(X) ≈ tr(X'X)
        adjoint_test(tr_At_A, randn(rng), X)
    end
end


text:
#=
Assorted utility functions used throughout the code.
Dependencies: none
=#

using LinearAlgebra, Gurobi, Distributions

"""
	cost_least_squares(X, c, intercept=false)

Solve a least-squares problem to fit a d x p matrix B such that c approx == B*X

# Arguments
- `X`: p x n feature matrix
- `c`: d x n matrix of cost vectors
"""
function cost_least_squares(X, c; intercept=false)
	(p, n) = size(X)
	(d, n2) = size(c)

	if n != n2
		error("dimensions are mismatched")
	end

	if intercept
		X = [X; ones(1,n)]
	end

	Xt = X'
	ct = c'

	Bt = Xt\ct

	return Bt'
end

"""
	ridge(X, c, reg_param)

Solve a ridge regression problem min_B (1/n)0.5*||c - B*X||^2 + 0.5*reg_param*||B||^2
Here ||.|| is the Frobenius norm

# Arguments
- `X`: p x n feature matrix
- `c`: d x n matrix of cost vectors
"""
function ridge(X, c, reg_param)
	(p, n) = size(X)
	(d, n2) = size(c)

	if n != n2
		error("dimensions are mismatched")
	end

	Xt = X'
	ct = c'

	Bt = (X*Xt + n*reg_param*eye(p))\(X*ct)

	return Bt'
end

"""
	oracle_dataset(c, oracle)

Apply the optimization oracle to each column of the d x n matrix c, return
an n vector of z_star values and a d x n matrix of optimal solutions
"""
function oracle_dataset(c, oracle)
	(d, n) = size(c)
	z_star_data = zeros(n)
    w_star_data = zeros(d, n)
    for i = 1:n
        (z_i, w_i) = oracle(c[:,i])
        z_star_data[i] = z_i
        w_star_data[:,i] = w_i
    end
	return (z_star_data, w_star_data)
end

"""
	spo_loss(B_new, X, c, oracle; z_star=[])

Compute the SPO loss of `B_new` on the training/holdout/test set (`X`, `c`),
given the optimization oracle, `oracle`.

# Arguments
- `B_new`: d x p model matrix B
- `X`: p x n feature matrix
- `c`: d x n matrix of cost vectors
- `oracle`: the optimization oracle function
- `z_star`: optional vector of the z^\ast(c_i) values
"""
function spo_loss(B_new, X, c, oracle; z_star=[])

	if z_star == []
		(z_star, w_star) = oracle_dataset(c, oracle)
	end

	n = length(z_star)

	spo_sum = 0
	for i = 1:n
		c_hat = B_new*X[:,i]
		(z_oracle, w_oracle) = oracle(c_hat)
		spo_loss_cur = dot(c[:,i], w_oracle) - z_star[i]
		spo_sum = spo_sum + spo_loss_cur
	end
	spo_loss_avg = spo_sum/n
	return spo_loss_avg
end

"""
	compare_models_percent(B_modA, B_modB, X, c, oracle; eps=0.000001)

Compute the percentage of times on the training/test set (X, c) such that modA produces a
solution wtih smaller true cost than modB.

# Arguments
- `B_modA`: d x p model matrix B of modA
- `B_modB`: d x p model matrix B of modB
- `X`: p x n feature matrix
- `c`: d x n matrix of (true) cost vectors
- `oracle`: the optimization oracle function
- `eps`: tolerance paramter. modA beats modB if cost(modA) <= cost(modB) + eps
"""
function compare_models_percent(B_modA, B_modB, X, c, oracle; eps=0.000001)
	(d, n) = size(c)

	count = 0
	for i = 1:n
		c_hatA = B_modA*X[:,i]
		c_hatB = B_modB*X[:,i]

		(z_oracleA, w_oracleA) = oracle(c_hatA)
		(z_oracleB, w_oracleB) = oracle(c_hatB)

		diff = dot(c[:,i], w_oracleA) - dot(c[:,i], w_oracleB)
		if diff <= eps
			count = count + 1
		end
	end
	return count/n
end

"""
	spo_plus_loss(B_new, X, c, oracle; z_star=[], w_star=[])

Compute the SPO plus loss of `B_new` on the training/holdout/test set (`X`, `c`),
given the optimization oracle, `oracle`.

# Arguments
- `B_new`: d x p model matrix B
- `X`: p x n feature matrix
- `c`: d x n matrix of cost vectors
- `oracle`: the optimization oracle function
- `z_star`: optional vector of the z^\ast(c_i) values
- `w_star`: optional d x n matrix of w^\ast(c_i)
"""
function spo_plus_loss(B_new, X, c, oracle; z_star=[], w_star=[])

	if z_star == [] || w_star == []
		(z_star, w_star) = oracle_dataset(c, oracle)
	end

	n = length(z_star)

	spo_plus_sum = 0
	for i = 1:n
		c_hat = B_new*X[:,i]
		spoplus_cost_vec = 2*c_hat - c[:,i]
		(z_oracle, w_oracle) = oracle(spoplus_cost_vec)
		spo_plus_cost = -z_oracle + 2*dot(c_hat, w_star[:,i]) - z_star[i]
		spo_plus_sum = spo_plus_sum + spo_plus_cost
	end
	spo_plus_avg = spo_plus_sum/n
	return spo_plus_avg
end

"""
	least_squares_loss(B_new, X, c)

Compute the least squares loss of `B_new` on the training/holdout/test set (`X`, `c`).

# Arguments
- `B_new`: d x p model matrix B
- `X`: p x n feature matrix
- `c`: d x n matrix of cost vectors
"""
function least_squares_loss(B_new, X, c)
	(p, n) = size(X)
	residuals = B_new*X - c
	error = (1/n)*(1/2)*(norm(residuals)^2)

	return error
end

"""
Absoulte (i.e. ell_1) loss
"""
function absolute_loss(B_new, X, c)
	(p, n) = size(X)
	residuals = B_new*X - c
	error = (1/n)*norm(residuals, 1)

	return error
end

"""
Huber loss with parameter delta
"""
function huber_loss(B_new, X, c, delta)
	(p, n) = size(X)
	(d, n_2) = size(c)
	if n != n_2
		error("Dimension mismatch in Huber loss calucation.")
	end

	residuals = B_new*X - c
	total_error = 0
	for i = 1:n
		for j = 1:d
			a_res = abs(residuals[j, i])
			if a_res <= delta
				cur_error = 0.5*a_res^2
			else
				cur_error = delta*(a_res - 0.5*delta)
			end
			total_error = total_error + cur_error
		end
	end

	error = (1/n)*total_error
	return error
end

"""
	generate_poly_kernel_data(B_true, n, degree; inner_constant=1, outer_constant = 1, kernel_damp_normalize=true,
		kernel_damp_factor=1, noise=true, noise_half_width=0, normalize_c=true)

Generate (X, c) from the polynomial kernel model X_{ji} ~ N(0, 1) and
c_i(j) = ( (alpha_j * B_true[j,:] * X[:,i]  + inner_constant)^degree + outer_constant ) * epsilon_{ij} where
alpha_j is a damping term and epsilon_{ij} is a noise term.

# Arguments
- `kernel_damp_normalize`: if true, then set
alpha_j = kernel_damp_factor/norm(B_true[j,:]). This results in
(alpha_j * B_true[j,:] * X[:,i]  + inner_constant) being normally distributed with
mean inner_constant and standard deviation kernel_damp_factor.
- `noise`:  if true, generate epsilon_{ij} ~  Uniform[1 - noise_half_width, 1 + noise_half_width]
- `normalize_c`:  if true, normalize c at the end of everything
"""
function generate_poly_kernel_data(B_true, n, degree; inner_constant=1, outer_constant = 1, kernel_damp_normalize=true,
	kernel_damp_factor=1, noise=true, noise_half_width=0, normalize_c=true, normalize_small_threshold = 0.0001)

    (d, p) = size(B_true)
    X_observed = randn(p, n)
    dot_prods = B_true*X_observed

	# first generate c_observed without noise
    c_observed = zeros(d, n)
    for j = 1:d
		if kernel_damp_normalize
			cur_kernel_damp_factor = kernel_damp_factor/norm(B_true[j,:])
		else
			cur_kernel_damp_factor = kernel_damp_factor
		end
        for i = 1:n
			c_observed[j, i] = (cur_kernel_damp_factor*dot_prods[j, i] + inner_constant)^degree + outer_constant
			if noise
				epsilon = (1 - noise_half_width) + 2*noise_half_width*rand()
				c_observed[j, i] = c_observed[j, i]*epsilon
			end
        end
    end

	if normalize_c
		for i = 1:n
			c_observed[:, i] = c_observed[:, i]/norm(c_observed[:, i])
			for j = 1:d
				if abs(c_observed[j, i]) < normalize_small_threshold
					c_observed[j, i] = 0
				end
			end
		end
	end

    return (X_observed, c_observed)
end

function generate_poly_kernel_data_simple(B_true, n, polykernel_degree, noise_half_width)
	(d, p) = size(B_true)
	alpha_factor = 1/sqrt(p)

	if noise_half_width == 0
		noise_on = false
	else
		noise_on = true
	end

	return generate_poly_kernel_data(B_true, n, polykernel_degree;
        kernel_damp_factor = alpha_factor,
		noise = noise_on, noise_half_width = noise_half_width,
        kernel_damp_normalize = false,
        normalize_c = false,
		inner_constant = 3)
end

"""
Computes sigmoid function of z, z can be a scalar or vector or matrix
"""
function sigmoid(z)
	return 1.0 ./ (1.0 .+ exp.(-z))
end

function generate_sigmoid_data(p_features, d_feasibleregion, n_sigmoid, n_samples, bernoulli_prob)
	# True matrices parameterizing network
	B_true_1 = rand(Bernoulli(bernoulli_prob), n_sigmoid, p_features)

	# This step ensures that the input to the sigmoid functions is a normal RV with SD of 5
	col_norms = sqrt.(vec(sum(B_true_1.^2, dims = 2)))
	for j = 1:n_sigmoid
		if col_norms[j] < 0.01
			col_norms[j] = 5
		end
	end
	B_true_1 = Diagonal(5 ./ col_norms)*B_true_1

	# Second layer matrix
	B_true_2 = rand(Bernoulli(bernoulli_prob), d_feasibleregion, n_sigmoid)

	# Generate data
	X_observed = randn(p_features, n_samples)
	c_observed = B_true_2*sigmoid(B_true_1*X_observed)

	return (X_observed, c_observed)
end

function generate_sigmoid_returns_data(p_features, d_feasibleregion, n_sigmoid, n_samples, bernoulli_prob,
	L_matrix, sigma_noise, return_intercept)

	(X, r) = generate_sigmoid_data(p_features, d_feasibleregion, n_sigmoid, n_samples, bernoulli_prob)

	(d, f) = size(L_matrix)
	for i = 1:n_samples
		r[:,i] = return_intercept*ones(d) + r[:,i] + L_matrix*randn(f) + sigma_noise*randn(d)
	end

	c = -r
	return (X, c)
end

"""
	generate_rbf_kernel_data(B_true, n, sigma)

Generate (X, c) from the RBF kernel model X_{ji} ~ N(0, 1) and
c[j,i] = exp(-||B[j,:] -  X[:,i]||^2 / 2sigma^2)
"""
function generate_rbf_kernel_data(B_true, n, sigma)
	(d, p) = size(B_true)
    X_observed = randn(p, n)

	c_observed = zeros(d, n)
	for j = 1:d
		for i = 1:n
			dist = norm(B_true[j,:] - X_observed[:,i])
			c_observed[j,i] = exp(-dist^2/(2*sigma^2))
		end
	end

	return (X_observed, c_observed)
end

"""
	generate_rbf_kernel_returns_data(B_true, n, sigma_kernel, L_matrix,
		sigma_noise, return_intercept)

Returns are generated according to:
1. (X, r_kernel) = generate_rbf_kernel_data(B_true, n, sigma_kernel)
2. r = return_intercept + r_kernel + L_matrix*f + sigma_noise*eps where f and eps are
vectors of i.i.d. normals. This ensures that, conditional on the features x, the returns
have the desired covariance structure Sigma = L_matrix*L_matrix' + sigma_noise^2*eye(d),
which is independent of x.

The method returns (X, c) where c = -r
"""
function generate_rbf_kernel_returns_data(B_true, n, sigma_kernel, L_matrix,
	sigma_noise, return_intercept)

	(X, r) = generate_rbf_kernel_data(B_true, n, sigma_kernel)

	(d, f) = size(L_matrix)
	for i = 1:n
		r[:,i] = return_intercept*ones(d) + r[:,i] + L_matrix*randn(f) + sigma_noise*randn(d)
	end

	c = -r
	return (X, c)
end

function generate_poly_kernel_returns_data(B_true, n, polykernel_degree,
	L_matrix, sigma_noise, return_intercept)

	(d, p) = size(B_true)
	alpha_factor = 0.05/sqrt(p)

	(X, r) = generate_poly_kernel_data(B_true, n, polykernel_degree;
        kernel_damp_factor = alpha_factor,
		noise = false,
        kernel_damp_normalize = false,
        normalize_c = false,
		inner_constant = 0.10^(1/polykernel_degree),
		outer_constant = 0)

	(d, f) = size(L_matrix)
	for i = 1:n
		r[:,i] = return_intercept*ones(d) + r[:,i] + L_matrix*randn(f) + sigma_noise*randn(d)
	end

	c = -r
	return (X, c)
end

function setup_gurobi_env(; quiet_mode = true, method_type = :barrier, use_time_limit = true, time_limit = 60.0)
	env = Gurobi.Env()

	if quiet_mode
		setparams!(env; OutputFlag = 0)
	end

	if method_type == :barrier
		setparams!(env; Method = 2)
	elseif method_type == :method3
		setparams!(env; Method = 3)
	elseif method_type != :default
		error("Enter a valid method type for Gurobi.")
	end

	if use_time_limit
		setparams!(env; TimeLimit = time_limit)
	end

	return env
end

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize


text:
## REGISTERING LABELS OF OBJECTS DURING ASSIGNMENT

const HANDLE_GIVEN_ID = Dict{UInt64,Symbol}()
SHOW_COLOR = true
"""
    color_on()

Enable color and bold output at the REPL, for enhanced display of MLJ objects.

"""
color_on() = (global SHOW_COLOR=true;)
"""
    color_off()

Suppress color and bold output at the REPL for displaying MLJ objects.

"""
color_off() = (global SHOW_COLOR=false;)

macro colon(p)
    Expr(:quote, p)
end

"""
    @constant x = value

Equivalent to `const x = value` but registers the binding thus:

    MLJBase.HANDLE_GIVEN_ID[objectid(value)] = :x

Registered objects get displayed using the variable name to which it
was bound in calls to `show(x)`, etc.

WARNING: As with any `const` declaration, binding `x` to new value of
the same type is not prevented and the registration will not be updated.

"""
macro constant(ex)
    ex.head == :(=) || throw(error("Expression must be an assignment."))
    handle = ex.args[1]
    value = ex.args[2]
    quote
        const $(esc(handle)) = $(esc(value))
        id = objectid($(esc(handle)))
        HANDLE_GIVEN_ID[id] = @colon $handle
        $(esc(handle))
    end
end

"""to display abbreviated versions of integers"""
function abbreviated(n)
    as_string = string(n)
    return as_string[1]*"…"*as_string[end-1:end]
end

"""return abbreviated object id (as string)  or it's registered handle (as string) if this exists"""
function handle(X)
    id = objectid(X)
    if id in keys(HANDLE_GIVEN_ID)
        return string(HANDLE_GIVEN_ID[id])
    else
        return abbreviated(id)
    end
end

## SHOW METHOD FOR NAMED TUPLES

# long version of showing a named tuple:
Base.show(stream::IO, ::MIME"text/plain", t::NamedTuple) = pretty_nt(stream, t)
pretty_nt(t) = pretty_nt(stdout, t)
pretty_nt(stream::IO, t) = pretty_nt(stream, t, 0)
pretty_nt(stream, t, n) = print(stream, t)
function pretty_nt(stream, t::NamedTuple, n)
    print(stream, "(")
    first_item = true
    for k in keys(t)
        value =  getproperty(t, k)
        if !first_item
            print(stream, crind(n + 1))
        else
            first_item = false
        end
        print(stream, "$k = ")
        pretty_nt(stream, value, n + length("$k = ") + 1)
        print(stream, ",")
    end
    print(stream, ")")
end

## OTHER EXPOSED SHOW METHODS

# string consisting of carriage return followed by indentation of length n:
crind(n) = "\n"*repeat(' ', n)

# trait to tag those objects to be displayed as constructed:
show_as_constructed(::Any) = false
show_as_constructed(::Type{<:Model}) = true

# short version of showing a `MLJType` object:
function Base.show(stream::IO, object::MLJType)
    id = objectid(object)
    description = string(typeof(object).name.name)
    parameters = typeof(object).parameters
    p_string = ""
    if length(parameters) == 1
        p = parameters[1]
        p isa DataType && (p_string = string(p.name.name))
    end
    description *= "{"*p_string*"}"
    str = "$description @ $(handle(object))"
    if !isempty(fieldnames(typeof(object)))
        printstyled(IOContext(stream, :color=> SHOW_COLOR), str, bold=false, color=:blue)
    else
        print(stream, str)
    end
end

# longer versions of showing objects
function Base.show(stream::IO, T::MIME"text/plain", object::MLJType)
    show(stream, T, object, Val(show_as_constructed(typeof(object))))
end

# fallback:
function Base.show(stream::IO, ::MIME"text/plain", object, ::Val{false})
    show(stream, MIME("text/plain"), object)
end

# fallback for MLJType:
function Base.show(stream::IO, ::MIME"text/plain", object::MLJType, ::Val{false})
    _recursive_show(stream, object, 1, DEFAULT_SHOW_DEPTH)
end

function Base.show(stream::IO, ::MIME"text/plain", object, ::Val{true})
    pretty(stream, object)
end
pretty(stream::IO, object) = pretty(stream, object, 0, DEFAULT_SHOW_DEPTH + 1, 0)
pretty(stream, object, current_depth, depth, n) = show(stream, object)
function pretty(stream, object::M, current_depth, depth, n) where M<:MLJType
    if current_depth == depth
        show(stream, object)
    else
        prefix = string(coretype(typeof(object)))
        print(stream, prefix, "(")
        first_item = true
        for k in fieldnames(M)
            value =  getproperty(object, k)
            if !first_item
                print(stream, crind(n + length(prefix) + 1))
            else
                first_item = false
            end
            print(stream, "$k = ")
            pretty(stream, value, current_depth + 1, depth, n + length(prefix) + 1 + length("$k = "))
            print(stream, ",")
        end
        print(stream, ")")
        if current_depth == 0
            description = " @ $(handle(object))"
            printstyled(IOContext(stream, :color=> SHOW_COLOR), description, bold=false, color=:blue)
        end
    end
end

# version showing a `MLJType` object to arbitrary depth:
Base.show(stream::IO, object::M, depth::Int) where M<:MLJType =
    show(stream, object, depth, Val(show_as_constructed(M)))
Base.show(stream::IO, object::MLJType, depth::Int, ::Val{false}) =
    _recursive_show(stream, object, 1, depth)
Base.show(stream::IO, object::MLJType, depth::Int, ::Val{true}) =
    pretty(stream, object, 0, 100, 0)

# for convenience:
Base.show(object::MLJType, depth::Int) = show(stdout, object, depth)

"""
    @more

Entered at the REPL, equivalent to `show(ans, 100)`. Use to get a
recursive description of all fields of the last REPL value.

"""
macro more()
    esc(quote
        show(Main.ans, 100)
    end)
end

## METHODS TO SUPRESS THE DISPLAY OF LARGE NON-BASETYPE OBJECTS

istoobig(::Any) = true
istoobig(::DataType) = false
istoobig(::UnionAll) = false
istoobig(::Union) = false
istoobig(::Number) = false
istoobig(::Char) = false
istoobig(::Function) = false
istoobig(s::Symbol) = false
istoobig(str::AbstractString) = length(str) > 50

## THE `_show` METHOD

# Note: The `_show` method controls how field values are displayed in
# the table generated by `_recursive_show`. See top of file.

# _show fallback:
function _show(stream::IO, object)
    if !istoobig(object)
        show(stream, MIME("text/plain"), object)
        println(stream)
    else
        println(stream, "(omitted ", typeof(object), ")")
    end
end

_show(stream::IO, object::MLJType) = println(stream, object)

# _show for other types:

istoobig(t::Tuple{Vararg{T}}) where T<:Union{Number,Symbol,Char,MLJType} =
    length(t) > 5
function _show(stream::IO, t::Tuple)
    if !istoobig(t)
        show(stream, MIME("text/plain"), t)
        println(stream)
    else
        println(stream, "(omitted $(typeof(t)) of length $(length(t)))")
    end
end

istoobig(A::AbstractArray{T}) where T<:Union{Number,Symbol,Char,MLJType} =
    maximum(size(A)) > 5
function _show(stream::IO, A::AbstractArray)
    if !istoobig(A)
        show(stream, MIME("text/plain"), A)
        println(stream)
    else
        println(stream, "(omitted $(typeof(A)) of size $(size(A)))")
    end
end

istoobig(d::Dict{T,Any}) where T <: Union{Number,Symbol,Char,MLJType} =
    length(keys(d)) > 5
function _show(stream::IO, d::Dict{T, Any}) where T <: Union{Number,Symbol}
    if isempty(d)
        println(stream, "empty $(typeof(d))")
    elseif !istoobig(d)
        println(stream, "omitted $(typeof(d)) with keys: ")
        show(stream, MIME("text/plain"), collect(keys(d)))
        println(stream)
    else
        println(stream, "(omitted $(typeof(d)))")
    end
end

function _show(stream::IO, v::Array{T, 1}) where T
    if !istoobig(v)
        show(stream, MIME("text/plain"), v)
        println(stream)
    else
        println(stream, "(omitted Vector{$T} of length $(length(v)))")
    end
end

_show(stream::IO, T::DataType) = println(stream, T)

_show(stream::IO, ::Nothing) = println(stream, "nothing")

## THE RECURSIVE SHOW METHOD

"""
    _recursive_show(stream, object, current_depth, depth)

Generate a table of the field values of the `MLJType` object,
dislaying each value by calling the method `_show` on it. The
behaviour of `_show(stream, f)` is as follows:

1. If `f` is itself a `MLJType` object, then its short form is shown
and `_recursive_show` generates as separate table for each of its
field values (and so on, up to a depth of argument `depth`).

2. Otherwise `f` is displayed as "(omitted T)" where `T = typeof(f)`,
unless `istoobig(f)` is false (the `istoobig` fall-back for arbitrary
types being `true`). In the latter case, the long (ie,
MIME"plain/text") form of `f` is shown. To override this behaviour,
overload the `_show` method for the type in question.

"""
function _recursive_show(stream::IO, object::MLJType, current_depth, depth)
    if depth == 0 || isempty(fieldnames(typeof(object)))
        println(stream, object)
    elseif current_depth <= depth
        fields = fieldnames(typeof(object))
        print(stream, "#"^current_depth, " ")
        show(stream, object)
        println(stream, ": ")
#        println(stream)
        if isempty(fields)
            println(stream)
            return
        end
        for fld in fields
            fld_string = string(fld)*" "^(max(0,COLUMN_WIDTH - length(string(fld))))*"=>   "
            print(stream, fld_string)
            if isdefined(object, fld)
                _show(stream, getfield(object, fld))
                #           println(stream)
            else
                println(stream, "(undefined)")
                #           println(stream)
            end
        end
        println(stream)
        for fld in fields
            if isdefined(object, fld)
                subobject = getfield(object, fld)
                if isa(subobject, MLJType) && !isempty(fieldnames(typeof(subobject)))
                    _recursive_show(stream, getfield(object, fld),
                                    current_depth + 1, depth)
                end
            end
        end
    end
end


text:
using Documenter
import TypeStability

if length(ARGS) > 0
    tag = ARGS[1]
else
    tag = nothing
end

makedocs(modules = [TypeStability])

deploydocs(deps = Deps.pip("mkdocs", "python-markdown-math"),
           repo = "github.com/Collegeville/TypeStability.jl",
           julia  = "0.6.4") #Currently uses the old package manager


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

using SLiDE

READ_DIR = joinpath("data", "readfiles")

function read_to_check(x::String)
    y = read_file(x)
    d = Dict(k => sort(read_file(joinpath(SLIDE_DIR, y["Path"]..., ensurearray(v["name"])...)))
        for (k,v) in y["Input"])
    return Dict(Symbol(k) => edit_with(v, Rename.(propertynames(v), Symbol.(y["Input"][k]["col"])))
        for (k,v) in d)
end

y = read_file(joinpath("data", "readfiles", "list_sets.yml"));
set = Dict((length(ensurearray(k)) == 1 ? Symbol(k) : Tuple(Symbol.(k))) =>
    sort(read_file(joinpath(y["Path"]..., ensurearray(v)...)))[:,1] for (k,v) in y["Input"])

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using JuMP
using Gurobi

type OPModel
    m
    status
end

function create_model(op::OrienteeringProblem; output=0)
    m = Model(solver=Gurobi.GurobiSolver(OutputFlag=output))
    N = length(op)

    without_start = [1:op.start-1; op.start+1:N]
    without_stop = [1:op.stop-1; op.stop+1:N]
    without_both = intersect(without_start, without_stop)

    @defVar(m, x[1:N,1:N], Bin)
    @defVar(m, 2 <= u[without_start] <= N)

    @setObjective(m, Max, sum{ sum{op.r[i]*x[i,j], j=1:N}, i=1:N })

    @addConstraint(m, sum{x[op.start,j], j=without_start} == 1)
    @addConstraint(m, sum{x[i,op.stop], i=without_stop} == 1)

    @addConstraint(m, connectivity[k=without_both], sum{x[i,k], i=1:N} == sum{x[k,j], j=1:N})

    @addConstraint(m, once[k=1:N], sum{x[k,j], j=1:N} <= 1)
    @addConstraint(m, sum{ sum{op.distances[i,j]*x[i,j], j=1:N}, i=1:N } <= op.distance_limit)
    @addConstraint(m, nosubtour[i=without_start,j=without_start], u[i]-u[j]+1 <= (N-1)*(1-x[i,j]))

    if op.start != op.stop
        @addConstraint(m, sum{x[op.stop,i],i=1:N} == 0)
    end

    return OPModel(m, :NotSolved)
end

function solve!(opm::OPModel)
    opm.status = solve(opm.m)
end

type LeiferModel
    m
end

# returns a model with the phase one constraints
function create_leifer_model(op::OrienteeringProblem)
    m = Model(solver=Gurobi.GurobiSolver(OutputFlag=output))
    N = length(op)

    V = [1:N]
    without_start = [1:op.start-1; op.start+1:N]
    without_stop = [1:op.stop-1; op.stop+1:N]
    Vprime = intersect(without_start, without_stop)

    # relaxed
    @defVar(m, 0 <= x[1:N,1:N] <= 1) # (7')
    @defVar(m, 0 <= y[1:N] <= 1) # (8')

    @setObjective(m, Max, sum{ op.r[j]*y[j] j=1:N }) # (1)

    ## Initial constraints

    @addConstraint(m, sum{x[op.start,j], j=without_start} == 1) # (2)
    @addConstraint(m, sum{x[i,op.stop], i=without_stop} == 1) # (3)

    @addConstraint(m, conservation[j=Vprime], sum{x[i,j],i=1:N} + sum{x[j,k],k=1:N} == 2*y[j]) # (4)

    @addConstraint(m, sum{ sum{op.distances[i,j]*x[i,j], j=1:N}, i=1:N } <= op.distance_limit) # (5)

    ## Phase one constraints

    # (skip constraints from Theorm 1 for now)

    # Theorem 2

    # create suboptimal path

    # Theorem 3

    return LeiferModel(m)
end

function nearest_neighbor_heur(op::OrienteeringProblem)
    p = [op.start, op.stop]
    ti = 1
    ui = 2

    c = op.distances
    s = Set{Int}(Vprime)
    len = distance(op,p)
    ctu = c[p[ti],p[ui]]
    margin = op.distance_limit - len + ctu
    filter!(k->c[p[ti],k]+c[k,p[ui]]<=margin, s)

    while !isempty(S)
        (imaxval, imax) = findmax([op.r[i]/c[p[ti],i] for i in s])
        (jmaxval, jmax) = findmax([op.r[j]/c[j,p[tu]] for j in s])
        if imaxval > jmaxval
            insert!(p, ti+1, imax)
            ti += 1
        else
            insert!(p, ui, jmax)
        end

        len = distance(op,p)
        ctu = c[p[ti],p[ui]]
        margin = op.distance_limit - len + ctu
        filter!(k->c[p[ti],k]+c[k,p[ui]]<=margin, s)
    end

    return p
end


text:
using Base: @propagate_inbounds

using OffsetArrays
using JLD2

using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Fields

using Oceananigans.Grids: topology, total_size, interior_parent_indices, parent_index_range
using Oceananigans.Fields: show_location, interior_view_indices, data_summary, reduced_location

import Oceananigans.Fields: Field, set!, interior, indices
import Oceananigans.Architectures: architecture

struct FieldTimeSeries{LX, LY, LZ, K, I, D, G, T, B, χ} <: AbstractField{LX, LY, LZ, G, T, 4}
                   data :: D
                   grid :: G
    boundary_conditions :: B
                indices :: I
                  times :: χ

    function FieldTimeSeries{LX, LY, LZ, K}(data::D, grid::G, bcs::B,
                                            times::χ, indices::I) where {LX, LY, LZ, K, D, G, B, χ, I}
        T = eltype(data)
        return new{LX, LY, LZ, K, I, D, G, T, B, χ}(data, grid, bcs, indices, times)
    end
end

architecture(fts::FieldTimeSeries) = architecture(fts.grid)

#####
##### Constructors
#####

"""
    FieldTimeSeries{LX, LY, LZ}(grid, times, [FT=eltype(grid);]
                                indices = (:, :, :),
                                boundary_conditions = nothing)

Return a `FieldTimeSeries` at location `(LX, LY, LZ)`, on `grid`, at `times`.
"""
function FieldTimeSeries{LX, LY, LZ}(grid, times, FT=eltype(grid);
                                     indices = (:, :, :),
                                     boundary_conditions = nothing) where {LX, LY, LZ}

    Nt = length(times)
    arch = architecture(grid)
    loc = (LX, LY, LZ)
    space_size = total_size(loc, grid, indices)
    underlying_data = zeros(FT, arch, space_size..., Nt)
    data = offset_data(underlying_data, grid, loc, indices)

    return FieldTimeSeries{LX, LY, LZ, InMemory}(data, grid, boundary_conditions, times, indices)
end

"""
    FieldTimeSeries(path, name;
                    backend = InMemory(),
                    grid = nothing,
                    iterations = nothing,
                    times = nothing)

Returns a `FieldTimeSeries` for the field `name` describing a field's time history from a JLD2 file
located at `path`.

Keyword arguments
=================

- `backend`: `InMemory()` to load data into a 4D array or `OnDisk()` to lazily load data from disk
             when indexing into `FieldTimeSeries`.

- `grid`: A grid to associated with data, in the case that the native grid was not serialized
          properly.

- `iterations`: Iterations to load. Defaults to all iterations found in the file.

- `times`: Save times to load, as determined through an approximate floating point
           comparison to recorded save times. Defaults to times associated with `iterations`.
           Takes precedence over `iterations` if `times` is specified.
"""
FieldTimeSeries(path, name; backend=InMemory(), kwargs...) =
    FieldTimeSeries(path, name, backend; kwargs...)

#####
##### InMemory time serieses
#####

const InMemoryFieldTimeSeries{LX, LY, LZ} = FieldTimeSeries{LX, LY, LZ, InMemory}

struct UnspecifiedBoundaryConditions end

function FieldTimeSeries(path, name, backend::InMemory;
                         architecture = nothing,
                         grid = nothing,
                         location = nothing,
                         boundary_conditions = UnspecifiedBoundaryConditions(),
                         iterations = nothing,
                         times = nothing)

    file = jldopen(path)

    # Defaults
    isnothing(iterations)   && (iterations = parse.(Int, keys(file["timeseries/t"])))
    isnothing(times)        && (times      = [file["timeseries/t/$i"] for i in iterations])
    isnothing(location)     && (location   = file["timeseries/$name/serialized/location"])

    if boundary_conditions isa UnspecifiedBoundaryConditions
        boundary_conditions = file["timeseries/$name/serialized/boundary_conditions"]
    end

    indices = try
        file["timeseries/$name/serialized/indices"]
    catch
        (:, :, :)
    end

    isnothing(grid) && (grid = file["serialized/grid"])

    # Default to CPU if neither architecture nor grid is specified
    architecture = isnothing(architecture) ? (isnothing(grid) ? CPU() : Architectures.architecture(grid)) : architecture

    # This should be removed in a month or two (4/5/2022).
    grid = try
        on_architecture(architecture, grid)
    catch err # Likely, the grid has CuArrays in it...
        if grid isa RectilinearGrid # we can try...
            Nx = file["grid/Nx"]
            Ny = file["grid/Ny"]
            Nz = file["grid/Nz"]
            Hx = file["grid/Hx"]
            Hy = file["grid/Hy"]
            Hz = file["grid/Hz"]
            xᶠᵃᵃ = file["grid/xᶠᵃᵃ"]
            yᵃᶠᵃ = file["grid/yᵃᶠᵃ"]
            zᵃᵃᶠ = file["grid/zᵃᵃᶠ"]
            x = file["grid/Δxᶠᵃᵃ"] isa Number ? (xᶠᵃᵃ[1], xᶠᵃᵃ[Nx+1]) : xᶠᵃᵃ
            y = file["grid/Δyᵃᶠᵃ"] isa Number ? (yᵃᶠᵃ[1], yᵃᶠᵃ[Ny+1]) : yᵃᶠᵃ
            z = file["grid/Δzᵃᵃᶠ"] isa Number ? (zᵃᵃᶠ[1], zᵃᵃᶠ[Nz+1]) : zᵃᵃᶠ
            topo = topology(grid)

            # Reduce for Flat dimensions
            domain = NamedTuple((:x, :y, :z)[i] => (x, y, z)[i] for i=1:3 if topo[i] !== Flat)
            size = Tuple((Nx, Ny, Nz)[i] for i=1:3 if topo[i] !== Flat)
            halo = Tuple((Hx, Hy, Hz)[i] for i=1:3 if topo[i] !== Flat)

            RectilinearGrid(architecture; size, halo, topology=topo, domain...)
        else
            throw(err)
        end
    end

    close(file)

    LX, LY, LZ = location
    time_series = FieldTimeSeries{LX, LY, LZ}(grid, times; indices, boundary_conditions)

    set!(time_series, path, name)

    return time_series
end

Base.parent(fts::FieldTimeSeries) = parent(fts.data)

function Base.getindex(fts::InMemoryFieldTimeSeries, n::Int)
    underlying_data = view(parent(fts), :, :, :, n)
    data = offset_data(underlying_data, fts.grid, location(fts), fts.indices)
    boundary_conditions = fts.boundary_conditions
    indices = fts.indices
    return Field(location(fts), fts.grid; data, boundary_conditions, indices)
end

#####
##### set!
#####

"""
    Field(location, path, name, iter;
          grid = nothing,
          architecture = nothing,
          indices = (:, :, :),
          boundary_conditions = nothing)

Load a field called `name` saved in a JLD2 file at `path` at `iter`ation.
Unless specified, the `grid` is loaded from `path`.
"""
function Field(location, path::String, name::String, iter;
               grid = nothing,
               architecture = nothing,
               indices = (:, :, :),
               boundary_conditions = nothing)

    file = jldopen(path)

    # Default to CPU if neither architecture nor grid is specified
    architecture = isnothing(architecture) ?
        (isnothing(grid) ? CPU() : Architectures.architecture(grid)) :
        architecture

    grid = isnothing(grid) ?
        on_architecture(architecture, file["serialized/grid"]) : grid

    raw_data = arch_array(architecture, file["timeseries/$name/$iter"])

    close(file)

    data = offset_data(raw_data, grid, location, indices)

    return Field(location, grid; boundary_conditions, indices, data)
end

function set!(time_series::InMemoryFieldTimeSeries, path::String, name::String)

    file = jldopen(path)
    file_iterations = parse.(Int, keys(file["timeseries/t"]))
    file_times = [file["timeseries/t/$i"] for i in file_iterations]
    close(file)

    for (n, time) in enumerate(time_series.times)
        file_index = findfirst(t -> t ≈ time, file_times)
        file_iter = file_iterations[file_index]

        field_n = Field(location(time_series), path, name, file_iter,
                        indices = time_series.indices,
                        boundary_conditions = time_series.boundary_conditions,
                        grid = time_series.grid)

        set!(time_series[n], field_n)
    end

    return nothing
end

function set!(fts::FieldTimeSeries, fields_vector::AbstractVector{<:AbstractField})
    raw_data = parent(fts)
    file = jldopen(path)

    for (n, field) in enumerate(fields_vector)
        raw_data[:, :, :, n] .= parent(field)
    end

    close(file)

    return nothing
end

function interior(fts::FieldTimeSeries)
    loc = location(fts)
    topo = topology(fts.grid)
    sz = size(fts.grid)
    halo_sz = halo_size(fts.grid)

    i_interior = interior_parent_indices.(loc, topo, sz, halo_sz)

    indices = fts.indices
    i_view = interior_view_indices.(indices, i_interior)

    return view(parent(fts), i_view..., :)
end

interior(fts::FieldTimeSeries, I...) = view(interior(fts), I...)

indices(fts::FieldTimeSeries) = fts.indices

#####
##### OnDisk time serieses
#####

struct OnDiskData
    path :: String
    name :: String
end

function FieldTimeSeries(path, name, backend::OnDisk; architecture=nothing, grid=nothing)
    file = jldopen(path)

    if isnothing(grid)
        grid = on_architecture(architecture, file["serialized/grid"])
    end

    iterations = parse.(Int, keys(file["timeseries/t"]))
    times = [file["timeseries/t/$i"] for i in iterations]

    data = OnDiskData(path, name)
    LX, LY, LZ = file["timeseries/$name/serialized/location"]
    bcs = file["timeseries/$name/serialized/boundary_conditions"]
    indices = file["timeseries/$name/serialized/indices"]

    close(file)

    return FieldTimeSeries{LX, LY, LZ, OnDisk}(data, grid, bcs, times, indices)
end

#####
##### Methods
#####

# Include the time dimension.
@inline Base.size(fts::FieldTimeSeries) = (size(location(fts), fts.grid, fts.indices)..., length(fts.times))

@propagate_inbounds Base.getindex(f::FieldTimeSeries{LX, LY, LZ, InMemory}, i, j, k, n) where {LX, LY, LZ} = f.data[i, j, k, n]

function Base.getindex(fts::FieldTimeSeries{LX, LY, LZ, OnDisk}, n::Int) where {LX, LY, LZ}
    # Load data
    arch = architecture(fts)
    file = jldopen(fts.data.path)
    iter = keys(file["timeseries/t"])[n]
    raw_data = arch_array(architecture(fts), file["timeseries/$(fts.data.name)/$iter"])
    close(file)

    # Wrap Field
    loc = (LX, LY, LZ)
    field_data = offset_data(raw_data, fts.grid, loc, fts.indices)

    return Field(loc, fts.grid; indices=fts.indices, boundary_conditions=fts.boundary_conditions, data=field_data)
end

Base.setindex!(fts::FieldTimeSeries, val, inds...) = Base.setindex!(fts.data, val, inds...)
Base.parent(fts::FieldTimeSeries{LX, LY, LZ, OnDisk}) where {LX, LY, LZ} = nothing

#####
##### Basic support for reductions
#####
##### TODO: support for reductions across _time_ (ie when 4 ∈ dims)
#####

const FTS = FieldTimeSeries

for reduction in (:sum, :maximum, :minimum, :all, :any, :prod)
    reduction! = Symbol(reduction, '!')

    @eval begin

        # Allocating
        function Base.$(reduction)(f::Function, fts::FTS; dims=:, kw...)
            if dims isa Colon
                return Base.$(reduction)($(reduction)(f, fts[n]; kw...) for n in 1:length(fts.times))
            else
                T = filltype(Base.$(reduction!), fts)
                loc = LX, LY, LZ = reduced_location(location(fts); dims)
                times = fts.times
                rts = FieldTimeSeries{LX, LY, LZ}(grid, times, T; indices=fts.indices)
                return Base.$(reduction!)(f, rts, fts; kw...)
            end
        end

        Base.$(reduction)(fts::FTS; kw...) = Base.$(reduction)(identity, fts; kw...)

        function Base.$(reduction!)(f::Function,rts::FTS, fts::FTS; dims=:, kw...)
            dims isa Tuple && 4 ∈ dims && error("Reduction across the time dimension (dim=4) is not yet supported!")
            times = rts.times
            for n = 1:length(times)
                Base.$(reduction!)(f, rts[i], fts[i]; dims, kw...)
            end
            return rts
        end

        Base.$(reduction!)(rts::FTS, fts::FTS; kw...) = Base.$(reduction!)(identity, rts, fts; kw...)
    end
end

#####
##### Show methods
#####

backend_str(::InMemory) = "InMemory"
backend_str(::OnDisk) = "OnDisk"

#####
##### show
#####

function Base.summary(fts::FieldTimeSeries{LX, LY, LZ, K}) where {LX, LY, LZ, K}
    arch = architecture(fts)
    A = typeof(arch)
    return string("$(join(size(fts), "×")) FieldTimeSeries{$(backend_str(K()))} located at ", show_location(fts), " on ", A)
end

function Base.show(io::IO, fts::FieldTimeSeries)
    prefix = string(summary(fts), '\n',
                   "├── grid: ", summary(fts.grid), '\n',
                   "├── indices: ", fts.indices, '\n')

    suffix = string("└── data: ", summary(fts.data), '\n',
                    "    └── ", data_summary(fts))

    return print(io, prefix, suffix)
end


text:
using Org
using Test

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


text:
using SequentialMonteCarlo
using RNGPool
using SMCExamples.FiniteFeynmanKac
using Test

setRNGs(0)

d = 3
n = 10
ffk = FiniteFeynmanKac.randomFiniteFK(d, n)
ffkout = FiniteFeynmanKac.calculateEtasZs(ffk)
model = FiniteFeynmanKac.makeSMCModel(ffk)
smcio = SMCIO{model.particle, model.pScratch}(2^13, n, Threads.nthreads(),
  false, 0.5)
smc!(model, smcio)

@test smcio.logZhats ≈ ffkout.logZhats atol=0.2

lM = FiniteFeynmanKac.makelM(ffk)
p1 = smcio.internal.zetaAncs[1]
p2 = smcio.zetas[1]
lM(n, p1, p2, smcio.internal.particleScratch)

## just test that the commands run; actual testing is part of the tests in
## SequentialMonteCarlo.jl

FiniteFeynmanKac.fullDensity(ffk, FiniteFeynmanKac.Int642Path(1, d, n))
FiniteFeynmanKac.allEtas(ffkout, p -> p.x, false)
FiniteFeynmanKac.allEtas(ffkout, p -> p.x, true)
FiniteFeynmanKac.allGammas(ffkout, p -> p.x, false)
FiniteFeynmanKac.allGammas(ffkout, p -> p.x, true)
FiniteFeynmanKac.vpns(ffk, ffkout, p -> p.x, true, true, smcio.resample)
FiniteFeynmanKac.avar(ffk, ffkout, p -> p.x, true, true, smcio.resample)
FiniteFeynmanKac.allavarhat1s(ffk, ffkout, smcio.resample)

@test FiniteFeynmanKac.Path2Int64(FiniteFeynmanKac.Int642Path(1, d, n), d) == 1


text:
function modalparami(A,B,C,dt,dofu)
	#
	# Compute natural frequencies, damping ratios, mode shapes
	# and modal masses from A, C, D and dt
	#
	# dofu: degrees of freedon where inputs are applied
	# 	example : consider we measure at points (a,b,c,d,e) and
	# 	the inputs are located at points b,e => dofu = [2,5]
	#
	# javier.cara@upm.es, 2016-02
	#

	# size
	no,ns = size(C)

	# invariants
	# ---------------------------------
	# eigenvalues and eigenvectors
	D,V = eigen(A)
	# L
	Lambda = log.(D)/dt

	# invariante V^{-1}*B
	VinvB = V\B

	# invariante C*V
	CV = C*V

	# modal parameters
	# ---------------------------------
	# allocating
	wm1 = zeros(ns)
	zm1 = zeros(ns)
	vm1 = zeros(no,ns)
	mm1 = zeros(ns)

	s = 1
	ss = 1
	while s < ns
		if abs(imag(D[s])) > 1e-10 # complex number
			wm1[ss] = abs(Lambda[s])
			zm1[ss] = -real(Lambda[s])/wm1[ss]

			vm0 = CV[:,s]/Lambda[s]^2
			# maximum element = 1
			maxval = vm0[1]
			for j in 2:no
				if abs(vm0[j]) > abs(maxval)
					maxval = vm0[j]
				end
			end
			vm1[:,ss] = real(vm0/maxval) # in theory, it should be a real vector

			c = (D[s] - 1.0)/( Lambda[s]*(Lambda[s] - conj(Lambda[s])) ) * (vm1[dofu,ss:ss]')/(VinvB[s:s,:]*maxval) # c is a vector with 1 component
			mm1[ss] = abs(real( c[1] ) )

			# next value
			s = s+2
			ss = ss+1

		 else
		    # real eigenvalue
		    s = s+1
		 end
	end

	# deleting zero values
	wm = wm1[1:ss-1]
	zm = zm1[1:ss-1]
	vm = vm1[:,1:ss-1]
	mm = mm1[1:ss-1]

	# sorting frequencies
	pos = sortperm(wm)
	wm = wm[pos]
	zm = zm[pos]
	vm = vm[:,pos]
	mm = mm[pos]

	return wm,zm,vm,mm

end


text:
# Tabix Index
# ===========

type Tabix
    format::Int32
    columns::NTuple{3,Int}
    meta::Char
    skip::Int
    names::Vector{String}
    indexes::Vector{Tuple{BinIndex,LinearIndex,Nullable{PseudoBin}}}
    n_no_coor::Nullable{UInt64}
end

function Base.show(io::IO, index::Tabix)
    println(io, summary(index), ":")
    println(io, "  format: ", format2str(index.format))
    println(io, "  columns: ", index.columns)
    println(io, "  meta char: '", index.meta, "'")
    println(io, "  skip lines: ", index.skip)
      print(io, "  names: ", index.names)
end

"""
    Tabix(filename::AbstractString)
    Tabix(input::IO)
    read(filename::AbstractString, ::Type{Tabix})
    read(input::IO, ::Type{Tabix})

Load a Tabix index from a file.
"""
function Tabix(filename::AbstractString)
    return read(filename, Tabix)
end

function Tabix(input::IO)
    return read(input, Tabix)
end

# The BED rule is half-closed-half-open and 0-based like Python.
function is_bed_rule(format)
    return format & 0x10000 != 0
end

function format2str(format)
    if format == 1
        return "SAM"
    elseif format == 2
        return "VCF"
    else
        if is_bed_rule(format)
            return "generic (BED rule)"
        else
            return "generic"
        end
    end
end


text:
# Julia wrapper for header: cvode.h
# Automatically generated using Clang.jl

function CVodeCreate(lmm::Cint)
    ccall((:CVodeCreate, libsundials_cvode), CVODEMemPtr, (Cint,), lmm)
end

function CVodeCreate(lmm)
    CVodeCreate(convert(Cint, lmm))
end

function CVodeInit(cvode_mem, f::CVRhsFn, t0::realtype, y0::N_Vector)
    ccall((:CVodeInit, libsundials_cvode), Cint, (CVODEMemPtr, CVRhsFn, realtype, N_Vector), cvode_mem, f, t0, y0)
end

function CVodeInit(cvode_mem, f, t0, y0)
    __y0 = convert(NVector, y0)
    CVodeInit(cvode_mem, CVRhsFn_wrapper(f), t0, convert(N_Vector, __y0))
end

function CVodeReInit(cvode_mem, t0::realtype, y0::N_Vector)
    ccall((:CVodeReInit, libsundials_cvode), Cint, (CVODEMemPtr, realtype, N_Vector), cvode_mem, t0, y0)
end

function CVodeReInit(cvode_mem, t0, y0)
    __y0 = convert(NVector, y0)
    CVodeReInit(cvode_mem, t0, convert(N_Vector, __y0))
end

function CVodeSStolerances(cvode_mem, reltol::realtype, abstol::realtype)
    ccall((:CVodeSStolerances, libsundials_cvode), Cint, (CVODEMemPtr, realtype, realtype), cvode_mem, reltol, abstol)
end

function CVodeSVtolerances(cvode_mem, reltol::realtype, abstol::N_Vector)
    ccall((:CVodeSVtolerances, libsundials_cvode), Cint, (CVODEMemPtr, realtype, N_Vector), cvode_mem, reltol, abstol)
end

function CVodeSVtolerances(cvode_mem, reltol, abstol)
    __abstol = convert(NVector, abstol)
    CVodeSVtolerances(cvode_mem, reltol, convert(N_Vector, __abstol))
end

function CVodeWFtolerances(cvode_mem, efun::CVEwtFn)
    ccall((:CVodeWFtolerances, libsundials_cvode), Cint, (CVODEMemPtr, CVEwtFn), cvode_mem, efun)
end

function CVodeSetErrHandlerFn(cvode_mem, ehfun::CVErrHandlerFn, eh_data)
    ccall((:CVodeSetErrHandlerFn, libsundials_cvode), Cint, (CVODEMemPtr, CVErrHandlerFn, Ptr{Cvoid}), cvode_mem, ehfun, eh_data)
end

function CVodeSetErrFile(cvode_mem, errfp)
    ccall((:CVodeSetErrFile, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{FILE}), cvode_mem, errfp)
end

function CVodeSetUserData(cvode_mem, user_data)
    ccall((:CVodeSetUserData, libsundials_cvode), Cint, (CVODEMemPtr, Any), cvode_mem, user_data)
end

function CVodeSetMaxOrd(cvode_mem, maxord::Cint)
    ccall((:CVodeSetMaxOrd, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, maxord)
end

function CVodeSetMaxOrd(cvode_mem, maxord)
    CVodeSetMaxOrd(cvode_mem, convert(Cint, maxord))
end

function CVodeSetMaxNumSteps(cvode_mem, mxsteps::Clong)
    ccall((:CVodeSetMaxNumSteps, libsundials_cvode), Cint, (CVODEMemPtr, Clong), cvode_mem, mxsteps)
end

function CVodeSetMaxNumSteps(cvode_mem, mxsteps)
    CVodeSetMaxNumSteps(cvode_mem, convert(Clong, mxsteps))
end

function CVodeSetMaxHnilWarns(cvode_mem, mxhnil::Clong)
    ccall((:CVodeSetMaxHnilWarns, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, mxhnil)
end

function CVodeSetMaxHnilWarns(cvode_mem, mxhnil)
    CVodeSetMaxHnilWarns(cvode_mem, convert(Clong, mxhnil))
end

function CVodeSetStabLimDet(cvode_mem, stldet::Cint)
    ccall((:CVodeSetStabLimDet, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, stldet)
end

function CVodeSetStabLimDet(cvode_mem, stldet)
    CVodeSetStabLimDet(cvode_mem, convert(Cint, stldet))
end

function CVodeSetInitStep(cvode_mem, hin::realtype)
    ccall((:CVodeSetInitStep, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, hin)
end

function CVodeSetMinStep(cvode_mem, hmin::realtype)
    ccall((:CVodeSetMinStep, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, hmin)
end

function CVodeSetMaxStep(cvode_mem, hmax::realtype)
    ccall((:CVodeSetMaxStep, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, hmax)
end

function CVodeSetStopTime(cvode_mem, tstop::realtype)
    ccall((:CVodeSetStopTime, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, tstop)
end

function CVodeSetMaxErrTestFails(cvode_mem, maxnef::Cint)
    ccall((:CVodeSetMaxErrTestFails, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, maxnef)
end

function CVodeSetMaxErrTestFails(cvode_mem, maxnef)
    CVodeSetMaxErrTestFails(cvode_mem, convert(Cint, maxnef))
end

function CVodeSetMaxNonlinIters(cvode_mem, maxcor::Cint)
    ccall((:CVodeSetMaxNonlinIters, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, maxcor)
end

function CVodeSetMaxNonlinIters(cvode_mem, maxcor)
    CVodeSetMaxNonlinIters(cvode_mem, convert(Cint, maxcor))
end

function CVodeSetMaxConvFails(cvode_mem, maxncf::Cint)
    ccall((:CVodeSetMaxConvFails, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, maxncf)
end

function CVodeSetMaxConvFails(cvode_mem, maxncf)
    CVodeSetMaxConvFails(cvode_mem, convert(Cint, maxncf))
end

function CVodeSetNonlinConvCoef(cvode_mem, nlscoef::realtype)
    ccall((:CVodeSetNonlinConvCoef, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, nlscoef)
end

function CVodeSetConstraints(cvode_mem, constraints::N_Vector)
    ccall((:CVodeSetConstraints, libsundials_cvode), Cint, (CVODEMemPtr, N_Vector), cvode_mem, constraints)
end

function CVodeSetConstraints(cvode_mem, constraints)
    __constraints = convert(NVector, constraints)
    CVodeSetConstraints(cvode_mem, convert(N_Vector, __constraints))
end

function CVodeSetNonlinearSolver(cvode_mem, NLS::SUNNonlinearSolver)
    ccall((:CVodeSetNonlinearSolver, libsundials_cvode), Cint, (CVODEMemPtr, SUNNonlinearSolver), cvode_mem, NLS)
end

function CVodeRootInit(cvode_mem, nrtfn::Cint, g::CVRootFn)
    ccall((:CVodeRootInit, libsundials_cvode), Cint, (CVODEMemPtr, Cint, CVRootFn), cvode_mem, nrtfn, g)
end

function CVodeRootInit(cvode_mem, nrtfn, g)
    CVodeRootInit(cvode_mem, convert(Cint, nrtfn), CVRootFn_wrapper(g))
end

function CVodeSetRootDirection(cvode_mem, rootdir)
    ccall((:CVodeSetRootDirection, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Cint}), cvode_mem, rootdir)
end

function CVodeSetNoInactiveRootWarn(cvode_mem)
    ccall((:CVodeSetNoInactiveRootWarn, libsundials_cvode), Cint, (CVODEMemPtr,), cvode_mem)
end

function CVode(cvode_mem, tout::realtype, yout::N_Vector, tret, itask::Cint)
    ccall((:CVode, libsundials_cvode), Cint, (CVODEMemPtr, realtype, N_Vector, Ptr{realtype}, Cint), cvode_mem, tout, yout, tret, itask)
end

function CVode(cvode_mem, tout, yout::N_Vector, tret, itask)
    __yout = convert(NVector, yout)
    CVode(cvode_mem, tout, convert(N_Vector, __yout), tret, convert(Cint, itask))
end

function CVodeGetDky(cvode_mem, t::realtype, k::Cint, dky::N_Vector)
    ccall((:CVodeGetDky, libsundials_cvode), Cint, (CVODEMemPtr, realtype, Cint, N_Vector), cvode_mem, t, k, dky)
end

function CVodeGetDky(cvode_mem, t, k, dky)
    __dky = convert(NVector, dky)
    CVodeGetDky(cvode_mem, t, convert(Cint, k), convert(N_Vector, __dky))
end

function CVodeGetWorkSpace(cvode_mem, lenrw, leniw)
    ccall((:CVodeGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrw, leniw)
end

function CVodeGetNumSteps(cvode_mem, nsteps)
    ccall((:CVodeGetNumSteps, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nsteps)
end

function CVodeGetNumRhsEvals(cvode_mem, nfevals)
    ccall((:CVodeGetNumRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevals)
end

function CVodeGetNumLinSolvSetups(cvode_mem, nlinsetups)
    ccall((:CVodeGetNumLinSolvSetups, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nlinsetups)
end

function CVodeGetNumErrTestFails(cvode_mem, netfails)
    ccall((:CVodeGetNumErrTestFails, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, netfails)
end

function CVodeGetLastOrder(cvode_mem, qlast)
    ccall((:CVodeGetLastOrder, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Cint}), cvode_mem, qlast)
end

function CVodeGetCurrentOrder(cvode_mem, qcur)
    ccall((:CVodeGetCurrentOrder, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Cint}), cvode_mem, qcur)
end

function CVodeGetCurrentGamma(cvode_mem, gamma)
    ccall((:CVodeGetCurrentGamma, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, gamma)
end

function CVodeGetNumStabLimOrderReds(cvode_mem, nslred)
    ccall((:CVodeGetNumStabLimOrderReds, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nslred)
end

function CVodeGetActualInitStep(cvode_mem, hinused)
    ccall((:CVodeGetActualInitStep, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, hinused)
end

function CVodeGetLastStep(cvode_mem, hlast)
    ccall((:CVodeGetLastStep, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, hlast)
end

function CVodeGetCurrentStep(cvode_mem, hcur)
    ccall((:CVodeGetCurrentStep, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, hcur)
end

function CVodeGetCurrentState(cvode_mem, y)
    ccall((:CVodeGetCurrentState, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{N_Vector}), cvode_mem, y)
end

function CVodeGetCurrentTime(cvode_mem, tcur)
    ccall((:CVodeGetCurrentTime, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, tcur)
end

function CVodeGetTolScaleFactor(cvode_mem, tolsfac)
    ccall((:CVodeGetTolScaleFactor, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, tolsfac)
end

function CVodeGetErrWeights(cvode_mem, eweight::N_Vector)
    ccall((:CVodeGetErrWeights, libsundials_cvode), Cint, (CVODEMemPtr, N_Vector), cvode_mem, eweight)
end

function CVodeGetErrWeights(cvode_mem, eweight::N_Vector)
    __eweight = convert(NVector, eweight)
    CVodeGetErrWeights(cvode_mem, convert(N_Vector, __eweight))
end

function CVodeGetEstLocalErrors(cvode_mem, ele::N_Vector)
    ccall((:CVodeGetEstLocalErrors, libsundials_cvode), Cint, (CVODEMemPtr, N_Vector), cvode_mem, ele)
end

function CVodeGetEstLocalErrors(cvode_mem, ele)
    __ele = convert(NVector, ele)
    CVodeGetEstLocalErrors(cvode_mem, convert(N_Vector, __ele))
end

function CVodeGetNumGEvals(cvode_mem, ngevals)
    ccall((:CVodeGetNumGEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, ngevals)
end

function CVodeGetRootInfo(cvode_mem, rootsfound)
    ccall((:CVodeGetRootInfo, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Cint}), cvode_mem, rootsfound)
end

function CVodeGetIntegratorStats(cvode_mem, nsteps, nfevals, nlinsetups, netfails, qlast, qcur, hinused, hlast, hcur, tcur)
    ccall((:CVodeGetIntegratorStats, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Cint}, Ptr{Cint}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}), cvode_mem, nsteps, nfevals, nlinsetups, netfails, qlast, qcur, hinused, hlast, hcur, tcur)
end

function CVodeGetNumNonlinSolvIters(cvode_mem, nniters)
    ccall((:CVodeGetNumNonlinSolvIters, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nniters)
end

function CVodeGetNumNonlinSolvConvFails(cvode_mem, nncfails)
    ccall((:CVodeGetNumNonlinSolvConvFails, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nncfails)
end

function CVodeGetNonlinSolvStats(cvode_mem, nniters, nncfails)
    ccall((:CVodeGetNonlinSolvStats, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, nniters, nncfails)
end

function CVodeGetReturnFlagName(flag::Clong)
    ccall((:CVodeGetReturnFlagName, libsundials_cvode), Cstring, (Clong,), flag)
end

function CVodeGetReturnFlagName(flag)
    CVodeGetReturnFlagName(convert(Clong, flag))
end

function CVodeFree(cvode_mem)
    ccall((:CVodeFree, libsundials_cvode), Cvoid, (Ref{CVODEMemPtr},), cvode_mem)
end
# Julia wrapper for header: cvode_bandpre.h
# Automatically generated using Clang.jl

function CVBandPrecInit(cvode_mem, N::sunindextype, mu::sunindextype, ml::sunindextype)
    ccall((:CVBandPrecInit, libsundials_cvode), Cint, (CVODEMemPtr, sunindextype, sunindextype, sunindextype), cvode_mem, N, mu, ml)
end

function CVBandPrecGetWorkSpace(cvode_mem, lenrwLS, leniwLS)
    ccall((:CVBandPrecGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwLS, leniwLS)
end

function CVBandPrecGetNumRhsEvals(cvode_mem, nfevalsBP)
    ccall((:CVBandPrecGetNumRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevalsBP)
end
# Julia wrapper for header: cvode_bbdpre.h
# Automatically generated using Clang.jl

function CVBBDPrecInit(cvode_mem, Nlocal::sunindextype, mudq::sunindextype, mldq::sunindextype, mukeep::sunindextype, mlkeep::sunindextype, dqrely::realtype, gloc::CVLocalFn, cfn::CVCommFn)
    ccall((:CVBBDPrecInit, libsundials_cvode), Cint, (CVODEMemPtr, sunindextype, sunindextype, sunindextype, sunindextype, sunindextype, realtype, CVLocalFn, CVCommFn), cvode_mem, Nlocal, mudq, mldq, mukeep, mlkeep, dqrely, gloc, cfn)
end

function CVBBDPrecReInit(cvode_mem, mudq::sunindextype, mldq::sunindextype, dqrely::realtype)
    ccall((:CVBBDPrecReInit, libsundials_cvode), Cint, (CVODEMemPtr, sunindextype, sunindextype, realtype), cvode_mem, mudq, mldq, dqrely)
end

function CVBBDPrecGetWorkSpace(cvode_mem, lenrwBBDP, leniwBBDP)
    ccall((:CVBBDPrecGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwBBDP, leniwBBDP)
end

function CVBBDPrecGetNumGfnEvals(cvode_mem, ngevalsBBDP)
    ccall((:CVBBDPrecGetNumGfnEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, ngevalsBBDP)
end
# Julia wrapper for header: cvode_diag.h
# Automatically generated using Clang.jl

function CVDiag(cvode_mem)
    ccall((:CVDiag, libsundials_cvode), Cint, (CVODEMemPtr,), cvode_mem)
end

function CVDiagGetWorkSpace(cvode_mem, lenrwLS, leniwLS)
    ccall((:CVDiagGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwLS, leniwLS)
end

function CVDiagGetNumRhsEvals(cvode_mem, nfevalsLS)
    ccall((:CVDiagGetNumRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevalsLS)
end

function CVDiagGetLastFlag(cvode_mem, flag)
    ccall((:CVDiagGetLastFlag, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, flag)
end

function CVDiagGetReturnFlagName(flag::Clong)
    ccall((:CVDiagGetReturnFlagName, libsundials_cvode), Cstring, (Clong,), flag)
end

function CVDiagGetReturnFlagName(flag)
    CVDiagGetReturnFlagName(convert(Clong, flag))
end
# Julia wrapper for header: cvode_direct.h
# Automatically generated using Clang.jl

function CVDlsSetLinearSolver(cvode_mem, LS::SUNLinearSolver, A::SUNMatrix)
    ccall((:CVDlsSetLinearSolver, libsundials_cvode), Cint, (CVODEMemPtr, SUNLinearSolver, SUNMatrix), cvode_mem, LS, A)
end

function CVDlsSetJacFn(cvode_mem, jac::CVDlsJacFn)
    ccall((:CVDlsSetJacFn, libsundials_cvode), Cint, (CVODEMemPtr, CVDlsJacFn), cvode_mem, jac)
end

function CVDlsGetWorkSpace(cvode_mem, lenrwLS, leniwLS)
    ccall((:CVDlsGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwLS, leniwLS)
end

function CVDlsGetNumJacEvals(cvode_mem, njevals)
    ccall((:CVDlsGetNumJacEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, njevals)
end

function CVDlsGetNumRhsEvals(cvode_mem, nfevalsLS)
    ccall((:CVDlsGetNumRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevalsLS)
end

function CVDlsGetLastFlag(cvode_mem, flag)
    ccall((:CVDlsGetLastFlag, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, flag)
end

function CVDlsGetReturnFlagName(flag::Clong)
    ccall((:CVDlsGetReturnFlagName, libsundials_cvode), Cstring, (Clong,), flag)
end

function CVDlsGetReturnFlagName(flag)
    CVDlsGetReturnFlagName(convert(Clong, flag))
end
# Julia wrapper for header: cvode_ls.h
# Automatically generated using Clang.jl

function CVodeSetLinearSolver(cvode_mem, LS::SUNLinearSolver, A::SUNMatrix)
    ccall((:CVodeSetLinearSolver, libsundials_cvode), Cint, (CVODEMemPtr, SUNLinearSolver, SUNMatrix), cvode_mem, LS, A)
end

function CVodeSetJacFn(cvode_mem, jac::CVLsJacFn)
    ccall((:CVodeSetJacFn, libsundials_cvode), Cint, (CVODEMemPtr, CVLsJacFn), cvode_mem, jac)
end

function CVodeSetMaxStepsBetweenJac(cvode_mem, msbj::Clong)
    ccall((:CVodeSetMaxStepsBetweenJac, libsundials_cvode), Cint, (CVODEMemPtr, Clong), cvode_mem, msbj)
end

function CVodeSetMaxStepsBetweenJac(cvode_mem, msbj)
    CVodeSetMaxStepsBetweenJac(cvode_mem, convert(Clong, msbj))
end

function CVodeSetEpsLin(cvode_mem, eplifac::realtype)
    ccall((:CVodeSetEpsLin, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, eplifac)
end

function CVodeSetPreconditioner(cvode_mem, pset::CVLsPrecSetupFn, psolve::CVLsPrecSolveFn)
    ccall((:CVodeSetPreconditioner, libsundials_cvode), Cint, (CVODEMemPtr, CVLsPrecSetupFn, CVLsPrecSolveFn), cvode_mem, pset, psolve)
end

function CVodeSetJacTimes(cvode_mem, jtsetup::CVLsJacTimesSetupFn, jtimes::CVLsJacTimesVecFn)
    ccall((:CVodeSetJacTimes, libsundials_cvode), Cint, (CVODEMemPtr, CVLsJacTimesSetupFn, CVLsJacTimesVecFn), cvode_mem, jtsetup, jtimes)
end

function CVodeSetLinSysFn(cvode_mem, linsys::CVLsLinSysFn)
    ccall((:CVodeSetLinSysFn, libsundials_cvode), Cint, (CVODEMemPtr, CVLsLinSysFn), cvode_mem, linsys)
end

function CVodeGetLinWorkSpace(cvode_mem, lenrwLS, leniwLS)
    ccall((:CVodeGetLinWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwLS, leniwLS)
end

function CVodeGetNumJacEvals(cvode_mem, njevals)
    ccall((:CVodeGetNumJacEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, njevals)
end

function CVodeGetNumPrecEvals(cvode_mem, npevals)
    ccall((:CVodeGetNumPrecEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, npevals)
end

function CVodeGetNumPrecSolves(cvode_mem, npsolves)
    ccall((:CVodeGetNumPrecSolves, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, npsolves)
end

function CVodeGetNumLinIters(cvode_mem, nliters)
    ccall((:CVodeGetNumLinIters, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nliters)
end

function CVodeGetNumLinConvFails(cvode_mem, nlcfails)
    ccall((:CVodeGetNumLinConvFails, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nlcfails)
end

function CVodeGetNumJTSetupEvals(cvode_mem, njtsetups)
    ccall((:CVodeGetNumJTSetupEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, njtsetups)
end

function CVodeGetNumJtimesEvals(cvode_mem, njvevals)
    ccall((:CVodeGetNumJtimesEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, njvevals)
end

function CVodeGetNumLinRhsEvals(cvode_mem, nfevalsLS)
    ccall((:CVodeGetNumLinRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevalsLS)
end

function CVodeGetLastLinFlag(cvode_mem, flag)
    ccall((:CVodeGetLastLinFlag, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, flag)
end

function CVodeGetLinReturnFlagName(flag::Clong)
    ccall((:CVodeGetLinReturnFlagName, libsundials_cvode), Cstring, (Clong,), flag)
end

function CVodeGetLinReturnFlagName(flag)
    CVodeGetLinReturnFlagName(convert(Cint, flag))
end
# Julia wrapper for header: cvode_spils.h
# Automatically generated using Clang.jl

function CVSpilsSetLinearSolver(cvode_mem, LS::SUNLinearSolver)
    ccall((:CVSpilsSetLinearSolver, libsundials_cvode), Cint, (CVODEMemPtr, SUNLinearSolver), cvode_mem, LS)
end

function CVSpilsSetEpsLin(cvode_mem, eplifac::realtype)
    ccall((:CVSpilsSetEpsLin, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, eplifac)
end

function CVSpilsSetPreconditioner(cvode_mem, pset::CVSpilsPrecSetupFn, psolve::CVSpilsPrecSolveFn)
    ccall((:CVSpilsSetPreconditioner, libsundials_cvode), Cint, (CVODEMemPtr, CVSpilsPrecSetupFn, CVSpilsPrecSolveFn), cvode_mem, pset, psolve)
end

function CVSpilsSetJacTimes(cvode_mem, jtsetup::CVSpilsJacTimesSetupFn, jtimes::CVSpilsJacTimesVecFn)
    ccall((:CVSpilsSetJacTimes, libsundials_cvode), Cint, (CVODEMemPtr, CVSpilsJacTimesSetupFn, CVSpilsJacTimesVecFn), cvode_mem, jtsetup, jtimes)
end

function CVSpilsGetWorkSpace(cvode_mem, lenrwLS, leniwLS)
    ccall((:CVSpilsGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwLS, leniwLS)
end

function CVSpilsGetNumPrecEvals(cvode_mem, npevals)
    ccall((:CVSpilsGetNumPrecEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, npevals)
end

function CVSpilsGetNumPrecSolves(cvode_mem, npsolves)
    ccall((:CVSpilsGetNumPrecSolves, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, npsolves)
end

function CVSpilsGetNumLinIters(cvode_mem, nliters)
    ccall((:CVSpilsGetNumLinIters, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nliters)
end

function CVSpilsGetNumConvFails(cvode_mem, nlcfails)
    ccall((:CVSpilsGetNumConvFails, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nlcfails)
end

function CVSpilsGetNumJTSetupEvals(cvode_mem, njtsetups)
    ccall((:CVSpilsGetNumJTSetupEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, njtsetups)
end

function CVSpilsGetNumJtimesEvals(cvode_mem, njvevals)
    ccall((:CVSpilsGetNumJtimesEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, njvevals)
end

function CVSpilsGetNumRhsEvals(cvode_mem, nfevalsLS)
    ccall((:CVSpilsGetNumRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevalsLS)
end

function CVSpilsGetLastFlag(cvode_mem, flag)
    ccall((:CVSpilsGetLastFlag, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, flag)
end

function CVSpilsGetReturnFlagName(flag::Clong)
    ccall((:CVSpilsGetReturnFlagName, libsundials_cvode), Cstring, (Clong,), flag)
end

function CVSpilsGetReturnFlagName(flag)
    CVSpilsGetReturnFlagName(convert(Clong, flag))
end


text:
# Next, you need to find the largest basins so you know what areas are most
# important to avoid.
#
# A basin is all locations that eventually flow downward to a single low point.
# Therefore, every low point has a basin, although some basins are very small.
# Locations of height 9 do not count as being in any basin, and all other
# locations will always be part of exactly one basin.
#
# The size of a basin is the number of locations within the basin, including
# the low point. The example above has four basins.
#
# The top-left basin, size 3:
#
# 2199943210
# 3987894921
# 9856789892
# 8767896789
# 9899965678
#
# The top-right basin, size 9:
#
# 2199943210
# 3987894921
# 9856789892
# 8767896789
# 9899965678
#
# The middle basin, size 14:
#
# 2199943210
# 3987894921
# 9856789892
# 8767896789
# 9899965678
#
# The bottom-right basin, size 9:
#
# 2199943210
# 3987894921
# 9856789892
# 8767896789
# 9899965678
#
# Find the three largest basins and multiply their sizes together. In the above
# example, this is 9 * 14 * 9 = 1134.
#
# What do you get if you multiply together the sizes of the three largest
# basins?

# read input
const width = 100
const height = 100
field = Array{UInt16}(undef, height, width)
foreach(eachline("input.txt"), Iterators.countfrom()) do l,y
  field[y,:] = parse.(UInt16, split(l, ""))
end

# find the three largest basins
basins = Int[]
for x in 1:width, y in 1:height
  field[y, x] == 9 && continue

  if y > 1 && field[y-1, x] > 9
    if x > 1 && field[y, x-1] > 9 && field[y-1, x] != field[y, x-1]
      # merge two basins
      (consumed_basin, basinid) = (field[y-1, x], field[y, x-1])
      field[y, x] = basinid

      # set all elements from the consumed basin to the combined basin id -
      # because our loop is running in column-major order, we can iterate over
      # the linear indices of field to find our consumed basin
      foreach(1:LinearIndices(field)[y, x]) do i
        field[i] == consumed_basin && (field[i] = basinid)
      end

      # increase basinid's count and set consumed_basin to 0
      basins[basinid - 9] += basins[consumed_basin - 9]
      basins[consumed_basin - 9] = 0
    else
      field[y, x] = field[y-1, x]
    end
  elseif x > 1 && field[y, x-1] > 9
    field[y, x] = field[y, x-1]
  else
    # new basin - basin id is position in basins vector plus 9
    push!(basins, 0)
    field[y, x] = length(basins) + 9
  end

  # increase basin's size
  basins[field[y, x] - 9] += 1
end
sort!(basins; rev=true)

println("Result: ", prod(basins[1:3]))