text:
# ------------------------------------------------------------------------------------------
# # An Intro to DifferentialEquations.jl
#
# ## Basic Introduction Via Ordinary Differential Equations
#
# This notebook will get you started with DifferentialEquations.jl by introducing you to the
# functionality for solving ordinary differential equations (ODEs). The corresponding
# documentation page is the [ODE
# tutorial](http://docs.juliadiffeq.org/latest/tutorials/ode_example.html). While some of
# the syntax may be different for other types of equations, the same general principles hold
# in each case. Our goal is to give a gentle and thorough introduction that highlights these
# principles in a way that will help you generalize what you have learned.
#
# ### Background
#
# If you are new to the study of differential equations, it can be helpful to do a quick
# background read on [the definition of ordinary differential
# equations](https://en.wikipedia.org/wiki/Ordinary_differential_equation). We define an
# ordinary differential equation as an equation which describes the way that a variable $u$
# changes, that is
#
# $$ u' = f(u,p,t) $$
#
# where $p$ are the parameters of the model, $t$ is the time variable, and $f$ is the
# nonlinear model of how $u$ changes. The initial value problem also includes the
# information about the starting value:
#
# $$ u(t_0) = u_0 $$
#
# Together, if you know the starting value and you know how the value will change with time,
# then you know what the value will be at any time point in the future. This is the
# intuitive definition of a differential equation.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### First Model: Exponential Growth
#
# Our first model will be the canonical exponential growth model. This model says that the
# rate of change is proportional to the current value, and is this:
#
# $$ u' = au $$
#
# where we have a starting value $u(0)=u_0$. Let's say we put 1 dollar into Bitcoin which is
# increasing at a rate of $98\%$ per year. Then calling now $t=0$ and measuring time in
# years, our model is:
#
# $$ u' = 0.98u $$
#
# and $u(0) = 1.0$. We encode this into Julia by noticing that, in this setup, we match the
# general form when
#
# $$ f(u,p,t) = 0.98u $$
#
# with $ u_0 = 1.0 $. If we want to solve this model on a time span from `t=0.0` to `t=1.0`,
# then we define an `ODEProblem` by specifying this function `f`, this initial condition
# `u0`, and this time span as follows:
# ------------------------------------------------------------------------------------------

using DifferentialEquations
f(u,p,t) = 0.98u
u0 = 1.0
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)

# ------------------------------------------------------------------------------------------
# To solve our `ODEProblem` we use the command `solve`.
# ------------------------------------------------------------------------------------------

sol = solve(prob)

# ------------------------------------------------------------------------------------------
# and that's it: we have succesfully solved our first ODE!
#
# #### Analyzing the Solution
#
# Of course, the solution type is not interesting in and of itself. We want to understand
# the solution! The documentation page which explains in detail the functions for analyzing
# the solution is the [Solution
# Handling](http://docs.juliadiffeq.org/latest/basics/solution.html) page. Here we will
# describe some of the basics. You can plot the solution using the plot recipe provided by
# [Plots.jl](http://docs.juliaplots.org/latest/):
# ------------------------------------------------------------------------------------------

using Plots; gr()
plot(sol)

# ------------------------------------------------------------------------------------------
# From the picture we see that the solution is an exponential curve, which matches our
# intuition. As a plot recipe, we can annotate the result using any of the [Plots.jl
# attributes](http://docs.juliaplots.org/latest/attributes/). For example:
# ------------------------------------------------------------------------------------------

plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false

# ------------------------------------------------------------------------------------------
# Using the mutating `plot!` command we can add other pieces to our plot. For this ODE we
# know that the true solution is $u(t) = u_0 exp(at)$, so let's add some of the true
# solution to our plot:
# ------------------------------------------------------------------------------------------

plot!(sol.t, t->1.0*exp(0.98t),lw=3,ls=:dash,label="True Solution!")

# ------------------------------------------------------------------------------------------
# In the previous command I demonstrated `sol.t`, which grabs the array of time points that
# the solution was saved at:
# ------------------------------------------------------------------------------------------

sol.t

# ------------------------------------------------------------------------------------------
# We can get the array of solution values using `sol.u`:
# ------------------------------------------------------------------------------------------

sol.u

# ------------------------------------------------------------------------------------------
# `sol.u[i]` is the value of the solution at time `sol.t[i]`. We can compute arrays of
# functions of the solution values using standard comprehensions, like:
# ------------------------------------------------------------------------------------------

[t+u for (u,t) in tuples(sol)]

# ------------------------------------------------------------------------------------------
# However, one interesting feature is that, by default, the solution is a continuous
# function. If we check the print out again:
# ------------------------------------------------------------------------------------------

sol

# ------------------------------------------------------------------------------------------
# you see that it says that the solution has a 4th order interpolation, meaning that it is a
# continuous function of 4th order accuracy. We can call the solution as a function of time
# `sol(t)`. For example, to get the value at `t=0.45`, we can use the command:
# ------------------------------------------------------------------------------------------

sol(0.45)

# ------------------------------------------------------------------------------------------
# #### Controlling the Solver
#
# DifferentialEquations.jl has a common set of solver controls among its algorithms which
# can be found [at the Common Solver
# Options](http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html) page. We will
# detail some of the most widely used options.
#
# The most useful options are the tolerances `abstol` and `reltol`. These tell the internal
# adaptive time stepping engine how precise of a solution you want. Generally, `reltol` is
# the relative accuracy while `abstol` is the accuracy when `u` is near zero. These
# tolerances are local tolerances and thus are not global guarantees. However, a good rule
# of thumb is that the total solution accuracy is 1-2 digits less than the relative
# tolerances. Thus for the defaults `abstol=1e-6` and `reltol=1e-3`, you can expect a global
# accuracy of about 1-2 digits. If we want to get around 6 digits of accuracy, we can use
# the commands:
# ------------------------------------------------------------------------------------------

sol = solve(prob,abstol=1e-8,reltol=1e-8)

# ------------------------------------------------------------------------------------------
# Now we can see no visible difference against the true solution:
# ------------------------------------------------------------------------------------------

plot(sol)
plot!(sol.t, t->1.0*exp(0.98t),lw=3,ls=:dash,label="True Solution!")

# ------------------------------------------------------------------------------------------
# Notice that by decreasing the tolerance, the number of steps the solver had to take was
# `9` instead of the previous `5`. There is a trade off between accuracy and speed, and it
# is up to you to determine what is the right balance for your problem.
#
# Another common option is to use `saveat` to make the solver save at specific time points.
# For example, if we want the solution at an even grid of `t=0.1k` for integers `k`, we
# would use the command:
# ------------------------------------------------------------------------------------------

sol = solve(prob,saveat=0.1)

# ------------------------------------------------------------------------------------------
# Notice that when `saveat` is used the continuous output variables are no longer saved and
# thus `sol(t)`, the interpolation, is only first order. We can save at an uneven grid of
# points by passing a collection of values to `saveat`. For example:
# ------------------------------------------------------------------------------------------

sol = solve(prob,saveat=[0.2,0.7,0.9])

# ------------------------------------------------------------------------------------------
# By default it always saves the first and last values, but we can turn this off as well:
# ------------------------------------------------------------------------------------------

sol = solve(prob,saveat=[0.2,0.7,0.9],save_start = false, save_end = false)

# ------------------------------------------------------------------------------------------
# If we need to reduce the amount of saving, we can also turn off the continuous output
# directly via `dense=false`:
# ------------------------------------------------------------------------------------------

sol = solve(prob,dense=false)

# ------------------------------------------------------------------------------------------
# and to turn off all intermediate saving we can use `save_everystep=false`:
# ------------------------------------------------------------------------------------------

sol = solve(prob,save_everystep=false)

# ------------------------------------------------------------------------------------------
# More advanced saving behaviors, such as saving functionals of the solution, are handled
# via the `SavingCallback` in the [Callback Library](http://docs.juliadiffeq.org/latest/feat
# ures/callback_library.html#SavingCallback-1) which will be addressed later in the
# tutorial.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Choosing Solver Algorithms
#
# There is no best algorithm for numerically solving a differential equation. When you call
# `solve(prob)`, DifferentialEquations.jl makes a guess at a good algorithm for your
# problem, given the properties that you ask for (the tolerances, the saving information,
# etc.). However, in many cases you may want more direct control. A later notebook will help
# introduce the various *algorithms* in DifferentialEquations.jl, but for now let's
# introduce the *syntax*.
#
# The most crucial determining factor in choosing a numerical method is the stiffness of the
# model. Stiffness is roughly characterized by a Jacobian `f` with large eigenvalues. That's
# quite mathematical, and we can think of it more intuitively: if you have big numbers in
# `f` (like parameters of order `1e5`), then it's probably stiff. Or, as the creator of the
# MATLAB ODE Suite, Lawrence Shampine, likes to define it, if the standard algorithms are
# slow, then it's stiff. We will go into more depth about diagnosing stiffness in a later
# tutorial, but for now note that if you believe your model may be stiff, you can hint this
# to the algorithm chooser via `alg_hints = [:stiff]`.
# ------------------------------------------------------------------------------------------

sol = solve(prob,alg_hints=[:stiff])

# ------------------------------------------------------------------------------------------
# Stiff algorithms have to solve implicit equations and linear systems at each step so they
# should only be used when required.
#
# If we want to choose an algorithm directly, you can pass the algorithm type after the
# problem as `solve(prob,alg)`. For example, let's solve this problem using the `Tsit5()`
# algorithm, and just for show let's change the relative tolerance to `1e-6` at the same
# time:
# ------------------------------------------------------------------------------------------

sol = solve(prob,Tsit5(),reltol=1e-6)

# ------------------------------------------------------------------------------------------
# ### Systems of ODEs: The Lorenz Equation
#
# Now let's move to a system of ODEs. The [Lorenz
# equation](https://en.wikipedia.org/wiki/Lorenz_system) is the famous "butterfly attractor"
# that spawned chaos theory. It is defined by the system of ODEs:
#
# $$ \frac{dx}{dt} = \sigma (y - x) $$
# $$ \frac{dy}{dt} = x (\rho - z) -y $$
# $$ \frac{dz}{dt} = xy - \beta z $$
#
# To define a system of differential equations in DifferentialEquations.jl, we define our
# `f` as a vector function with a vector initial condition. Thus, for the vector `u =
# [x,y,z]'`, we have the derivative function:
# ------------------------------------------------------------------------------------------

function lorenz!(du,u,p,t)
    σ,ρ,β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
end

# ------------------------------------------------------------------------------------------
# Notice here we used the in-place format which writes the output to the preallocated vector
# `du`. For systems of equations the in-place format is faster. We use the initial condition
# $u_0 = [1.0,0.0,0.0]$ as follows:
# ------------------------------------------------------------------------------------------

u0 = [1.0,0.0,0.0]

# ------------------------------------------------------------------------------------------
# Lastly, for this model we made use of the parameters `p`. We need to set this value in the
# `ODEProblem` as well. For our model we want to solve using the parameters $\sigma = 10$,
# $\rho = 28$, and $\beta = 8/3$, and thus we build the parameter collection:
# ------------------------------------------------------------------------------------------

p = (10,28,8/3) # we could also make this an array, or any other type!

# ------------------------------------------------------------------------------------------
# Now we generate the `ODEProblem` type. In this case, since we have parameters, we add the
# parameter values to the end of the constructor call. Let's solve this on a time span of
# `t=0` to `t=100`:
# ------------------------------------------------------------------------------------------

tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)

# ------------------------------------------------------------------------------------------
# Now, just as before, we solve the problem:
# ------------------------------------------------------------------------------------------

sol = solve(prob)

# ------------------------------------------------------------------------------------------
# The same solution handling features apply to this case. Thus `sol.t` stores the time
# points and `sol.u` is an array storing the solution at the corresponding time points.
#
# However, there are a few extra features which are good to know when dealing with systems
# of equations. First of all, `sol` also acts like an array. `sol[i]` returns the solution
# at the `i`th time point.
# ------------------------------------------------------------------------------------------

sol.t[10],sol[10]

# ------------------------------------------------------------------------------------------
# Additionally, the solution acts like a matrix where `sol[j,i]` is the value of the `j`th
# variable at time `i`:
# ------------------------------------------------------------------------------------------

sol[2,10]

# ------------------------------------------------------------------------------------------
# We can get a real matrix by performing a conversion:
# ------------------------------------------------------------------------------------------

A = convert(Array,sol)

# ------------------------------------------------------------------------------------------
# This is the same as sol, i.e. `sol[i,j] = A[i,j]`, but now it's a true matrix. Plotting
# will by default show the time series for each variable:
# ------------------------------------------------------------------------------------------

plot(sol)

# ------------------------------------------------------------------------------------------
# If we instead want to plot values against each other, we can use the `vars` command. Let's
# plot variable `1` against variable `2` against variable `3`:
# ------------------------------------------------------------------------------------------

plot(sol,vars=(1,2,3))

# ------------------------------------------------------------------------------------------
# This is the classic Lorenz attractor plot, where the `x` axis is `u[1]`, the `y` axis is
# `u[2]`, and the `z` axis is `u[3]`. Note that the plot recipe by default uses the
# interpolation, but we can turn this off:
# ------------------------------------------------------------------------------------------

plot(sol,vars=(1,2,3),denseplot=false)

# ------------------------------------------------------------------------------------------
# Yikes! This shows how calculating the continuous solution has saved a lot of computational
# effort by computing only a sparse solution and filling in the values! Note that in vars,
# `0=time`, and thus we can plot the time series of a single component like:
# ------------------------------------------------------------------------------------------

plot(sol,vars=(0,2))

# ------------------------------------------------------------------------------------------
# ### A DSL for Parameterized Functions
#
# In many cases you may be defining a lot of functions with parameters. There exists the
# domain-specific language (DSL) defined by the `@ode_def` macro for helping with this
# common problem. For example, we can define the Lotka-Volterra equation:
#
# $$ \frac{dx}{dt} = ax - bxy $$
# $$ \frac{dy}{dt} = -cy + dxy $$
#
# as follows:
# ------------------------------------------------------------------------------------------

function lotka_volterra!(du,u,p,t)
  du[1] = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = -p[3]*u[2] + p[4]*u[1]*u[2]
end

# ------------------------------------------------------------------------------------------
# However, that can be hard to follow since there's a lot of "programming" getting in the
# way. Instead, you can use the `@ode_def` macro:
# ------------------------------------------------------------------------------------------

lv! = @ode_def LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d

# ------------------------------------------------------------------------------------------
# We can then use the result just like an ODE function from before:
# ------------------------------------------------------------------------------------------

u0 = [1.0,1.0]
p = (1.5,1.0,3.0,1.0)
tspan = (0.0,10.0)
prob = ODEProblem(lv!,u0,tspan,p)
sol = solve(prob)
plot(sol)

# ------------------------------------------------------------------------------------------
# Not only is the DSL convenient syntax, but it does some magic behind the scenes. For
# example, further parts of the tutorial will describe how solvers for stiff differential
# equations have to make use of the Jacobian in calculations. Here, the DSL uses symbolic
# differentiation to automatically derive that function:
# ------------------------------------------------------------------------------------------

lv!.Jex

# ------------------------------------------------------------------------------------------
# The DSL can derive many other functions; this ability is used to speed up the solvers. An
# extension to DifferentialEquations.jl, [Latexify.jl](https://korsbo.github.io/Latexify.jl/
# latest/tutorials/parameterizedfunctions.html), allows you to extract these pieces as LaTeX
# expressions.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Internal Types
#
# The last basic user-interface feature to explore is the choice of types.
# DifferentialEquations.jl respects your input types to determine the internal types that
# are used. Thus since in the previous cases, when we used `Float64` values for the initial
# condition, this meant that the internal values would be solved using `Float64`. We made
# sure that time was specified via `Float64` values, meaning that time steps would utilize
# 64-bit floats as well. But, by simply changing these types we can change what is used
# internally.
#
# As a quick example, let's say we want to solve an ODE defined by a matrix. To do this, we
# can simply use a matrix as input.
# ------------------------------------------------------------------------------------------

A  = [1. 0  0 -5
      4 -2  4 -3
     -4  0  0  1
      5 -2  2  3]
u0 = rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)

# ------------------------------------------------------------------------------------------
# There is no real difference from what we did before, but now in this case `u0` is a `4x2`
# matrix. Because of that, the solution at each time point is matrix:
# ------------------------------------------------------------------------------------------

sol[3]

# ------------------------------------------------------------------------------------------
# In DifferentialEquations.jl, you can use any type that defines `+`, `-`, `*`, `/`, and has
# an appropriate `norm`. For example, if we want arbitrary precision floating point numbers,
# we can change the input to be a matrix of `BigFloat`:
# ------------------------------------------------------------------------------------------

big_u0 = big.(u0)

# ------------------------------------------------------------------------------------------
# and we can solve the `ODEProblem` with arbitrary precision numbers by using that initial
# condition:
# ------------------------------------------------------------------------------------------

prob = ODEProblem(f,big_u0,tspan)
sol = solve(prob)

sol[1,3]

# ------------------------------------------------------------------------------------------
# To really make use of this, we would want to change `abstol` and `reltol` to be small!
# Notice that the type for "time" is different than the type for the dependent variables,
# and this can be used to optimize the algorithm via keeping multiple precisions. We can
# convert time to be arbitrary precision as well by defining our time span with `BigFloat`
# variables:
#
# ------------------------------------------------------------------------------------------

prob = ODEProblem(f,big_u0,big.(tspan))
sol = solve(prob)

# ------------------------------------------------------------------------------------------
# Let's end by showing a more complicated use of types. For small arrays, it's usually
# faster to do operations on static arrays via the package
# [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl). The syntax is similar
# to that of normal arrays, but for these special arrays we utilize the `@SMatrix` macro to
# indicate we want to create a static array.
# ------------------------------------------------------------------------------------------

using StaticArrays
A  = @SMatrix [ 1.0  0.0 0.0 -5.0
                4.0 -2.0 4.0 -3.0
               -4.0  0.0 0.0  1.0
                5.0 -2.0 2.0  3.0]
u0 = @SMatrix rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)

sol[3]

# ------------------------------------------------------------------------------------------
# ## Conclusion
#
# These are the basic controls in DifferentialEquations.jl. All equations are defined via a
# problem type, and the `solve` command is used with an algorithm choice (or the default) to
# get a solution. Every solution acts the same, like an array `sol[i]` with `sol.t[i]`, and
# also like a continuous function `sol(t)` with a nice plot command `plot(sol)`. The Common
# Solver Options can be used to control the solver for any equation type. Lastly, the types
# used in the numerical solving are determined by the input types, and this can be used to
# solve with arbitrary precision and add additional optimizations (this can be used to solve
# via GPUs for example!). While this was shown on ODEs, these techniques generalize to other
# types of equations as well.
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
@testset "AbstractValEdge" begin
    @test ValEdge <: AbstractValEdge
    @test ValDiEdge <: AbstractValEdge
end

@testset "Edge Constructors" for
    V       in TEST_VERTEX_TYPES_SMALL,
    E_VALS  in TEST_EDGEVAL_TYPES_SMALL,
    u       in V[1, 2, typemax(V)],
    v       in V[1, 2, typemax(V)]

    values = rand_sample(E_VALS)

    # We enforce that the constructor always puts the smaller vertex as src
    @testset "ValEdge($V($u), $V($v), $values)" begin

        src_should, dst_should = minmax(u, v)
        e = ValEdge(u, v, values)
        e_rev = ValEdge(v, u, values)

        @test e.src == e_rev.src == src_should
        @test e.dst == e_rev.dst == dst_should
        @test e.src <= e.dst
        @test e.vals == values
        @test e isa ValEdge{V, E_VALS}
    end

    @testset "ValDiEdge($V($u), $V($v), $values)" begin

        e = ValDiEdge(u, v, values)

        @test e.src == u
        @test e.dst == v
        @test e.vals == values
        @test e isa ValDiEdge{V, E_VALS}
    end
end

@testset "edge access functions" for
    V       in TEST_VERTEX_TYPES_SMALL,
    E_VALS  in TEST_EDGEVAL_TYPES_SMALL,
    u       in V[1, 2],
    v       in V[1, 2]

    values = rand_sample(E_VALS)

    e = ValEdge(u, v, values)
    e_rev = ValEdge(v, u, values)
    @testset "e == ValEdge($V($u), $V($v), $values)" begin

        @testset "src, dst" begin
            @test issetequal((src(e), dst(e)), (u, v))
            @test issetequal((src(e_rev), dst(e_rev)), (u, v))
            @test src(e) isa V
            @test dst(e) isa V
        end

        @testset "vals" begin
            @test vals(e) == values
            @test vals(e) isa E_VALS
        end

        @testset "val(e, key=$key)" for
            key in (keys(values) ∪ Base.OneTo(length(values))) # number & symbols

            @test val(e, key=key) == values[key]
            key isa Integer && @test val(e, key=key) isa E_VALS.types[key]
        end
    end

    e = ValDiEdge(u, v, values)
    @testset "e == ValDiEdge($V($u), $V($v), $values)" begin

        @testset "src, dst" begin
            @test src(e) == u
            @test dst(e) == v
            @test src(e) isa V
            @test dst(e) isa V
        end

        @testset "vals" begin
            @test vals(e) == values
            @test vals(e) isa E_VALS
        end

        @testset "val(e, key=$key)" for
            key in (keys(values) ∪ Base.OneTo(length(values))) # number & symbols

            @test val(e, key=key) == values[key]
            @test val(e, key=key) == vals(e)[key]
            key isa Integer && @test val(e, key=key) isa E_VALS.types[key]
        end
    end

end

@testset "edge access functions val(e) [no key]" begin

    @testset "exactly one value, E_VALS == $E_VALS" for
        E_VALS in TEST_EDGEVAL_TYPES_SINGE_VALUE_SMALL

        values = rand_sample(E_VALS)
        e_undir = ValEdge(1, 2, values)
        e_dir = ValEdge(1, 2, values)

        @test val(e_undir) == val(e_undir, key=1)
        @test val(e_undir) == vals(e_undir)[1]
        @test val(e_undir) isa E_VALS.types[1]

        @test val(e_dir) == val(e_dir, key=1)
        @test val(e_dir) == vals(e_dir)[1]
        @test val(e_dir) isa E_VALS.types[1]

    end

    @testset "less ore more than one value, E_VALS == $E_VALS" for
        E_VALS in TEST_EDGEVAL_TYPES_NON_SINGE_VALUE_SMALL

        values = rand_sample(E_VALS)
        e_undir = ValEdge(1, 2, values)
        e_dir = ValEdge(1, 2, values)

        @test_throws Exception val(e_undir)
        @test_throws Exception val(e_dir)
    end
end

@testset "reverse edge" begin

    @testset "edge == $e" for
        V       in TEST_VERTEX_TYPES_SMALL,
        E_VALS  in TEST_EDGEVAL_TYPES_SMALL,
        u       in V[1, 2],
        v       in V[1, 2],
        e       in ( ValEdge(u, v, rand_sample(E_VALS)), )

        e_rev = reverse(e)
        @test src(e) == src(e_rev)
        @test dst(e) == dst(e_rev)
        @test vals(e) == vals(e_rev)
        @test typeof(e) == typeof(e_rev)
    end

    @testset "edge == $e" for
        V       in TEST_VERTEX_TYPES_SMALL,
        E_VALS  in TEST_EDGEVAL_TYPES_SMALL,
        u       in V[1, 2],
        v       in V[1, 2],
        e       in ( ValDiEdge(u, v, rand_sample(E_VALS)), )

        e_rev = reverse(e)
        @test src(e) == dst(e_rev)
        @test dst(e) == src(e_rev)
        @test vals(e) == vals(e_rev)
        @test typeof(e) == typeof(e_rev)
    end

end

ext: jl


text:
struct StandardAggregation
end

function (::StandardAggregation)(S::SparseMatrixCSC{T,R}) where {T,R}

    n = size(S, 1)
    x = zeros(R, n)
    y = zeros(R, n)

    next_aggregate = 1

    for i = 1:n
        if x[i] != 0
            continue
        end

        has_agg_neighbors = false
        has_neighbors = false

        for j in nzrange(S, i)
            row = S.rowval[j]
            if row != i
                has_neighbors = true
                if x[row] != 0
                    has_agg_neighbors = true
                    break
                end
            end
        end

        if !has_neighbors
            x[i] = -n
        elseif !has_agg_neighbors
            x[i] = next_aggregate
            y[next_aggregate] = i

            for j in nzrange(S, i)
                row = S.rowval[j]
                x[row] = next_aggregate
            end

            next_aggregate += 1
        end
    end

    # Pass 2
    for i = 1:n
        if x[i] != 0
            continue
        end

        for j in nzrange(S, i)
            row = S.rowval[j]
            x_row = x[row]
            if x_row > 0
                x[i] = -x_row
                break
            end
        end
    end

    next_aggregate -= 1

    # Pass 3
    for i = 1:n
        xi = x[i]
        if xi != 0
            if xi > 0
                x[i] = xi - 1
            elseif xi == -n
                x[i] = -1
            else
                x[i] = -xi - 1
            end
            continue
        end

        x[i] = next_aggregate
        y[next_aggregate + 1] = i

        for j in nzrange(S, i)
            row = S.rowval[j]

            if x[row] == 0
                x[row] = next_aggregate
            end
        end

        next_aggregate += 1
    end

    y = y[1:next_aggregate]
    M,N = (n, next_aggregate)

    # Some nodes not aggregated
    if minimum(x) == -1
        mask = x .!= -1
        I = collect(R, 1:n)[mask]
        J = x[mask] .+ R(1)
        #J = x[mask] + 1
        V = ones(eltype(S), length(J))
        AggOp = sparse(J,I,V,N,M)
    else
        Tp = collect(R, 1:n+1)
        x .= x .+ R(1)
        #x .= x .+ 1
        Tx = ones(eltype(S), length(x))
        AggOp = SparseMatrixCSC(N, M, Tp, x, Tx)
    end

    AggOp
end


text:
module OpenCV
using OpenCV_jll
include(joinpath(OpenCV_jll.artifact_dir, "OpenCV", "src", "OpenCV.jl"))
end


text:
# GPS GLONASS Gelileo BeiDou etc

include("SatID.jl")
include("ObsID.jl")
include("Constant.jl")


text:
using NCEI
using Test: @test, @test_throws, @testset

const cdo_token = "ijrAfEVJJAefHdruUsTkjRyWLqSZpBhY" # Dedicated CDO token for testing

@testset "Testing NCEI" begin
    @testset "Endpoint: Data Categories" begin
    all_datacategories = cdo_datacategories(cdo_token)
    @test size(all_datacategories) == (42, 2)
    one_data_category = cdo_datacategories(cdo_token, "ANNAGR")
    @test size(one_data_category) == (1, 2)
    datacategory_location = cdo_datacategories(cdo_token, locations = "CITY:US390029")
    @test size(datacategory_location) == (38, 2)
    @test_throws ArgumentError("N7 is not a valid data category. For a complete list of valid data categories run `cdo_datacategories(CDO_token::AbstractString)`.") cdo_datacategories(cdo_token, "N7")
    end
    @testset "Endpoint: Datasets" begin
        all_datasets = cdo_datasets(cdo_token)
        @test size(all_datasets) == (11, 6)
        one_dataset = cdo_datasets(cdo_token, "GSOY")
        @test size(one_dataset) == (1, 5)
        one_datatype = cdo_datasets(cdo_token, datatypes = "TOBS")
        @test size(one_datatype) == (1, 6)
        station = cdo_datasets(cdo_token, stations = "GHCND:USC00010008")
        @test size(station) == (6, 6)
        @test_throws ArgumentError("N7 is not a valid dataset. For a complete list of valid datasets run `cdo_datasets(CDO_token::AbstractString)`.") cdo_datasets(cdo_token, "N7")
    end
    @testset "Endpoint: Data Types" begin
        all_datatypes = cdo_datatypes(cdo_token)
        @test size(all_datatypes) == (1527, 5)
        one_datatype = cdo_datatypes(cdo_token, "ACMH")
        @test size(one_datatype) == (1, 4)
        datacategory = cdo_datatypes(cdo_token, datacategories = "TEMP")
        @test size(datacategory) == (59, 5)
        datacategories = cdo_datatypes(cdo_token,
                                       stations = ["COOP:310090", "COOP:310184", "COOP:310212"])
        @test size(datacategories) == (21, 5)
        @test_throws ArgumentError("N7 is not a valid data type. For a complete list of valid data types run `cdo_datatypes(CDO_token::AbstractString)`.") cdo_datatypes(cdo_token, "N7")
    end
    @testset "Endpoint: Location Categories" begin
        all_locationcategories = cdo_locationcategories(cdo_token)
        @test size(all_locationcategories) == (12, 2)
        one_locationcategory = cdo_locationcategories(cdo_token, "CLIM_REG")
        @test size(one_locationcategory) == (1, 2)
        @test_throws ArgumentError("N7 is not a valid location category. For a complete list of valid location categories run `cdo_locationcategories(CDO_token::AbstractString)`.") cdo_locationcategories(cdo_token, "N7")
    end
    @testset "Endpoint: Locations" begin
        all_locations = cdo_locations(cdo_token,
                                      datasets = "GHCND",
                                      locationcategories = "ST")
        @test size(all_locations) == (51, 5)
        one_location = cdo_locations(cdo_token, "FIPS:37")
        @test size(one_location) == (1, 5)
        @test_throws ArgumentError("N7 is not a valid location. For a complete list of valid locations run `cdo_locations(CDO_token::AbstractString)`.") cdo_locations(cdo_token, "N7")
    end
    @testset "Endpoint: Stations" begin
        station = cdo_stations(cdo_token, "COOP:010008")
        @test size(station) == (1, 9)
        location = cdo_stations(cdo_token,
                                locations = "FIPS:37",
                                startdate = Date(2000, 1, 1),
                                enddate = Date(2000, 1, 1))
        @test size(location) == (417, 9)
        @test_throws ArgumentError("N7 is not a valid weather station. For a complete list of valid stations run `cdo_stations(CDO_token::AbstractString)`.") cdo_stations(cdo_token, "N7")
    end
    @testset "Endpoint: Data" begin
        GHCND = cdo_data(cdo_token,
                         "GHCND",
                         Date(2010, 5, 1),
                         Date(2010, 5, 1),
                         locations = "ZIP:28801")
        @test size(GHCND) == (8, 5)
        station = cdo_data(cdo_token,
                           "PRECIP_15",
                           Date(2010, 5, 1),
                           Date(2010, 5, 31),
                           stations = "COOP:010008")
        @test size(station) == (63, 5)
        GSOM = cdo_data(cdo_token,
                        "GSOM",
                        Date(2010, 5, 1),
                        Date(2010, 5, 31),
                        stations = "GHCND:USC00010008")
        @test size(GSOM) == (10, 5)
        Yearly = cdo_data(cdo_token,
                          "GHCND",
                          Date(2000, 1, 1),
                          Date(2002, 12, 31),
                          datatypes = "TAVG",
                          locations = "ZIP:91711")
        @test size(Yearly) == (1093, 5)
        Decade = cdo_data(cdo_token,
                          "GSOY",
                          Date(1990, 1, 1),
                          Date(2001, 12, 31),
                          datatypes = "TAVG",
                          locations = "ZIP:91711")
        @test size(Decade) == (3, 5)
    end
    @testset "Are these valid arguments?" begin
        @test_throws ArgumentError cdo_datacategories(cdo_token, locations = "California")
        @test_throws ArgumentError cdo_datasets(cdo_token, stations = "USC00010008")
        @test_throws ArgumentError cdo_datatypes(cdo_token, datacategories = "temp")
    end
end


text:
@testset "Shuffle" begin

    msa_types = (
        Matrix{Residue},
        NamedResidueMatrix{Array{Residue,2}},
        MultipleSequenceAlignment,
        AnnotatedMultipleSequenceAlignment
        )

    N = length(msa_types)

    pf09645_sto = joinpath(pwd(), "data", "PF09645_full.stockholm")

    msas = [ read(pf09645_sto, Stockholm, T) for T in msa_types ]
    gaps = [ msa .== GAP for msa in msas ]
    lcol = [ mean(msa .== Residue('L'), dims=1) for msa in msas ]
    lseq = [ mean(msa .== Residue('L'), dims=2) for msa in msas ]

    Random.seed!(42)

    @testset "General" begin

        msa = msas[1]

        for dim in [1,2]

            aln = copy(msa)
            shuffle!(aln, dim)
            @test aln != msa
            @test (aln .== GAP) == gaps[1] # default: fixed gaps

            aln = copy(msa)
            shuffle!(aln, dim, true)
            @test aln != msa
            @test (aln .== GAP) == gaps[1]

            aln = copy(msa)
            shuffle!(aln, dim, false)
            @test aln != msa
            @test (aln .== GAP) != gaps[1]
        end

        @test_throws AssertionError shuffle(msa, 0, true)
        @test_throws AssertionError shuffle(msa, 3, true)
    end

    @testset "Columns" begin

        for i in 1:N
            # Fixed gaps
            msa = msas[i]
            aln = shuffle(msa, 2, true)
            @test aln != getresidues(msa)
            @test (aln .== GAP) == gaps[i]
            @test lcol[i] == mean(aln .== Residue('L'), dims=1)
            @test lseq[i] != mean(aln .== Residue('L'), dims=2)
            # Change gap positions
            aln = shuffle(msa, 2, false)
            @test aln != getresidues(msa)
            @test (aln .== GAP) != gaps[i]
            @test lcol[i] == mean(aln .== Residue('L'), dims=1)
        end
    end

    @testset "Sequences" begin

        for  i in 1:N
            # Fixed gaps
            msa = msas[i]
            aln = shuffle(msa, 1, true)
            @test aln != getresidues(msa)
            @test (aln .== GAP) == gaps[i]
            @test lcol[i] != mean(aln .== Residue('L'), dims=1)
            @test lseq[i] == mean(aln .== Residue('L'), dims=2)
            # Change gap positions
            aln = shuffle(msa, 1, false)
            @test aln != getresidues(msa)
            @test (aln .== GAP) != gaps[i]
            @test lseq[i] == mean(aln .== Residue('L'), dims=2)
        end
    end
end


text:
using ForneyLab
using LinearAlgebra
import ForneyLab: SoftFactor, @ensureVariables, generateId, addNode!, associate!,
                  averageEnergy, differentialEntropy, Interface, Variable, slug, ProbabilityDistribution,
                  unsafeLogMean, unsafeMean, unsafeCov, unsafePrecision, unsafeMeanCov, cholinv
import SpecialFunctions: digamma
export GeneralisedFilterX, GFX

"""
Description:

    A Gaussian node for transitions between states in generalised coordinates:

    𝒩(y | A(θ)x + B(η)u, V(τ)),

    where for instance for an order-2 differential

        A(x) = (S + sθ')x,

    with S = | 1 Δt | and s = |  0  | ,
             | 0  1 |         | -Δt |

        B(η) = b η ,

    with b = |  0 | ,
             | Δt |
    and
        V(τ) = τ^-1 *[Δt^3/3   Δt^2/2;
                      Δt^2/2       Δt];

    Interfaces:
        1. y (output vector)
        2. θ (state coefficients)
        3. x (generalised coordinates)
        4. η (control coefficients)
        5. u (exogenous input)
        6. τ (precision)

    Construction:
        GeneralisedFilterX(y, θ, x, η, u, τ, Δt=1.0, id=:some_id)

"""

mutable struct GeneralisedFilterX <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    # Sampling time
    Δt::Float64

    function GeneralisedFilterX(y, θ, x, η, u, τ; Δt::Float64=1.0, id=generateId(GeneralisedFilterX))
        @ensureVariables(y, x, θ, η, u, τ)
        self = new(id, Array{Interface}(undef, 6), Dict{Symbol,Interface}(), Δt)
        addNode!(currentGraph(), self)
        self.i[:y] = self.interfaces[1] = associate!(Interface(self), y)
        self.i[:θ] = self.interfaces[2] = associate!(Interface(self), θ)
        self.i[:x] = self.interfaces[3] = associate!(Interface(self), x)
        self.i[:η] = self.interfaces[4] = associate!(Interface(self), η)
        self.i[:u] = self.interfaces[5] = associate!(Interface(self), u)
        self.i[:τ] = self.interfaces[6] = associate!(Interface(self), τ)
        return self
    end
end

slug(::Type{GeneralisedFilterX}) = "GFX"

function averageEnergy(::Type{GeneralisedFilterX},
                       marg_y::ProbabilityDistribution{Multivariate},
                       marg_x::ProbabilityDistribution{Multivariate},
                       marg_θ::ProbabilityDistribution{Multivariate},
                       marg_η::ProbabilityDistribution{Univariate},
                       marg_u::ProbabilityDistribution{Univariate},
                       marg_τ::ProbabilityDistribution{Univariate})

    #TODO
    error("not implemented yet")
end

max_stars_repo_name: wmkouw/gfx-node

max_issues_repo_name: wmkouw/GFX

max_forks_repo_name: wmkouw/GFX


text:
# Packages needed for parallel processing need to be loaded "twice":s
using Bio.Seq: DNASequence, DNAKmer, canonical, FASTASeqRecord, FASTAReader, each, neighbors
@everywhere using Bio.Seq: DNASequence, DNAKmer, canonical, FASTASeqRecord, FASTAReader, each, neighbors
using DataStructures: DefaultDict
@everywhere using DataStructures: DefaultDict
using FastaIO
@everywhere using FastaIO

@everywhere function twin{k}(::Type{DNAKmer{k}}, km)
    DNAKmer{k}(Bio.Seq.reverse_complement(DNASequence("$km")))
end

function contig_to_string(c)
  return "$(c[1])" * join(["$(x[end])" for x in c[2:end]])
end

@everywhere function get_contig{k}(::Type{DNAKmer{k}}, kmers, km)
  contig_fw = get_contig_forward(DNAKmer{k}, kmers, km)
  contig_bw = get_contig_forward(DNAKmer{k}, kmers, twin(DNAKmer{k}, km))

  if km in neighbors(contig_fw[end])
    return contig_fw
  else
    return [[twin(DNAKmer{k}, x) for x in reverse(contig_bw)[1:end-1]] ; contig_fw]
  end
end

@everywhere function get_contig_forward{k}(::Type{DNAKmer{k}}, kmers, km)
    c_fw = DNAKmer{k}[km]
    while true
        # check if forward exists, and only 1:
        fw_neighbors = [kmer for kmer in neighbors(c_fw[end]) if canonical(kmer) in kmers]
        if length(fw_neighbors) != 1
          break
        end
        candidate = fw_neighbors[1]
        if candidate == km || candidate == twin(DNAKmer{k}, km) || candidate == twin(DNAKmer{k}, c_fw[end])
          break
        end
        bw_neighbors = [kmer for kmer in neighbors(twin(DNAKmer{k}, candidate)) if canonical(kmer) in kmers]
        if length(bw_neighbors) != 1
          break
        end
        push!(c_fw, candidate)
    end
    return c_fw
end

@everywhere function all_contigs{k}(::Type{DNAKmer{k}}, kmers)
  done = Set{DNAKmer{k}}()
  contigs = []
  for kmer in kmers
    if kmer in done
      continue
    end
    contig = get_contig(DNAKmer{k}, kmers, kmer)
      for c_kmer in contig
      push!(done, canonical(c_kmer))
    end
    push!(contigs, contig)
  end
  return contigs
end

### classifies each kmer with his colors (present alleles)
@everywhere function find_all_kmer_colors{k}(::Type{DNAKmer{k}}, fastafile)
  record = FASTASeqRecord{DNASequence}()
  kmer_class = DefaultDict{DNAKmer{k}, Vector{Int16}}(() -> Int16[])
  allele_ids = Int16[]
  allele_idx::Int16 = 1
  open(FASTAReader{DNASequence}, fastafile) do reader
    while !eof(reader)
      read!(reader, record)
      seen = Set{DNAKmer{k}}()
      for (pos, kmer) in each(DNAKmer{k}, record.seq)
        can_kmer = canonical(kmer)
        if !in(can_kmer, seen)
          push!(kmer_class[can_kmer], allele_idx)
          push!(seen, can_kmer)
        end
      end
      # find the separator; will assume that if I see a "_", that's it, otherwise try "-";
      separator = in('_', record.name) ? "_" : "-"
      # update idx; the counter idx is incremental (1,2, ...) because we need the array sorted.
      # But this is not always sin the allele ordering, so we have to save the original id to restore it later;
      allele_id = parse(Int16,split(record.name, separator)[end])
      push!(allele_ids, allele_id)
      allele_idx += 1
    end
  end
  return kmer_class, allele_ids
end

@everywhere function build_db_graph{k}(::Type{DNAKmer{k}}, fastafile)
  # get kmers and colors (alleles)
  kmer_class, allele_ids = find_all_kmer_colors(DNAKmer{k}, fastafile)

  contig_list = all_contigs(DNAKmer{k}, keys(kmer_class))
  # select 1 kmer per contig, and set weights:
  kmer_weights = Dict{UInt64, Int}() # Kmer is converted to UInt64, because DNAKmer apparently does not support write(), needed for the parallel pmap() call;
  filtered_kmer_class = Dict{UInt64, Vector{Int16}}()
  for contig in contig_list
      kmer = canonical(contig[1])
      kmer_uint = convert(UInt64,kmer)
      kmer_weights[kmer_uint] = length(contig)
      filtered_kmer_class[kmer_uint] = kmer_class[kmer]
  end
  return (filtered_kmer_class, allele_ids, kmer_weights)
end

max_stars_repo_name: matnguyen/MentaLiST

max_issues_repo_name: matnguyen/MentaLiST

max_forks_repo_name: matnguyen/MentaLiST


text:
#=
Code to support input and output of data and configuration.
Input data can be loaded from NetCDF files.
Output data can be written to NetCDF or CSV files.
For configuration files we use TOML.
=#

"""Turn "a.aa.aaa" into (:a, :aa, :aaa)"""
symbols(s) = Tuple(Symbol(x) for x in split(s, '.'))

"""Turn symbols"a.aa.aaa" into (:a, :aa, :aaa)"""
macro symbols_str(s)
    Tuple(Symbol(x) for x in split(s, '.'))
end

"Get a nested field using a tuple of Symbols"
param(obj, fields) = foldl(getproperty, fields; init = obj)
param(obj, fields::AbstractString) = param(obj, symbols(fields))
function param(obj, fields, default)
    try
        return param(obj, fields)
    catch
        return default
    end
end

"""
    Config(path::AbstractString)
    Config(dict::AbstractDict)
    Config(dict::Dict{String,Any}, path::Union{String,Nothing})

Struct that contains the parsed TOML configuration, as well as a reference to the TOML path,
if it exists. It behaves largely like a distionary, but it overloads `getproperty` and
`setproperty` to support syntax like `config.model.reinit = false`.
"""
struct Config
    dict::Dict{String,Any}  # nested key value mapping of all settings
    path::Union{String,Nothing}  # path to the TOML file, or nothing
end

Config(path::AbstractString) = Config(TOML.parsefile(path), path)
Config(dict::AbstractDict) = Config(dict, nothing)

# allows using getproperty, e.g. config.input.time instead of config["input"]["time"]
function Base.getproperty(config::Config, f::Symbol)
    dict = Dict(config)
    path = pathof(config)
    a = dict[String(f)]
    # if it is a Dict, wrap the result in Config to keep the getproperty behavior
    return a isa AbstractDict ? Config(a, path) : a
end

function Base.setproperty!(config::Config, f::Symbol, x)
    dict = Dict(config)
    return dict[String(f)] = x
end

function get_alias(config::Config, key, alias, default)
    dict = Dict(config)
    path = pathof(config)
    alias_or_default = get(dict, alias, default)
    a = get(dict, key, alias_or_default)
    # if it is a Dict, wrap the result in Config
    return a isa AbstractDict ? Config(a, path) : a
end

# also used in autocomplete
Base.propertynames(config::Config) = collect(keys(Dict(config)))
Base.haskey(config::Config, key) = haskey(Dict(config), key)
Base.keys(config::Config) = keys(Dict(config))
Base.values(config::Config) = values(Dict(config))
Base.pairs(config::Config) = pairs(Dict(config))
Base.get(config::Config, key, default) = get(Dict(config), key, default)
Base.getindex(config::Config, key) = getindex(Dict(config), key)
Base.setindex!(config::Config, val, key) = setindex!(Dict(config), val, key)
Base.pop!(config::Config, key) = pop!(Dict(config), key)
Base.pop!(config::Config, key, default) = pop!(Dict(config), key, default)
Base.Dict(config::Config) = getfield(config, :dict)
Base.pathof(config::Config) = getfield(config, :path)
Base.dirname(config::Config) = dirname(pathof(config))
Base.iterate(config::Config) = iterate(Dict(config))
Base.iterate(config::Config, state) = iterate(Dict(config), state)

"Extract NetCDF variable name `ncname` from `var` (type `String` or `Config`). If `var` has
type `Config`, either `scale` and `offset` are expected (with `ncname`) or a `value` (uniform
value), these are stored as part of `NamedTuple` `modifier`"
function ncvar_name_modifier(var; verbose = true)
    ncname = nothing
    modifier = (scale = 1.0, offset = 0.0, value = nothing)
    if isa(var, Config)
        if haskey(var, "netcdf") &&
           haskey(var.netcdf, "variable") &&
           haskey(var.netcdf.variable, "name")
            ncname = var.netcdf.variable.name
            scale = param(var, "scale", 1.0)
            offset = param(var, "offset", 0.0)
            modifier = (scale = scale, offset = offset, value = nothing)
            if verbose
                @info "NetCDF parameter $ncname is modified with scale $scale and offset $offset"
            end
        elseif haskey(var, "value")
            modifier = (scale = 1.0, offset = 0.0, value = param(var, "value"))
        else
            error("Unrecognized modifier $(Dict(var))")
        end
    elseif isa(var, String)
        ncname = var
    else
        error("Unknown type")
    end
    return ncname, modifier
end

"Extract a NetCDF variable at a given time"
function get_at(
    ds::CFDataset,
    varname::AbstractString,
    times::AbstractVector{<:TimeType},
    t::TimeType,
)
    # this behaves like a forward fill interpolation
    i = findfirst(>=(t), times)
    t < first(times) && throw(DomainError("time $t before dataset begin $(first(times))"))
    i === nothing && throw(DomainError("time $t after dataset end $(last(times))"))
    return get_at(ds, varname, i)
end

function get_at(ds::CFDataset, varname::AbstractString, i)
    return read_standardized(ds, varname, (x = :, y = :, time = i))
end

function get_param_res(model)
    Dict(
        symbols"vertical.precipitation" => model.lateral.river.reservoir.precipitation,
        symbols"vertical.potential_evaporation" =>
            model.lateral.river.reservoir.evaporation,
    )
end

function get_param_lake(model)
    Dict(
        symbols"vertical.precipitation" => model.lateral.river.lake.precipitation,
        symbols"vertical.potential_evaporation" => model.lateral.river.lake.evaporation,
    )
end

mover_params = (symbols"vertical.precipitation", symbols"vertical.potential_evaporation")

function load_fixed_forcing(model)
    @unpack reader, network, config = model
    @unpack forcing_parameters = reader

    do_reservoirs = get(config.model, "reservoirs", false)::Bool
    do_lakes = get(config.model, "lakes", false)::Bool

    mover_params =
        (symbols"vertical.precipitation", symbols"vertical.potential_evaporation")
    reverse_indices = network.land.reverse_indices
    if do_reservoirs
        sel_reservoirs = network.reservoir.indices_coverage
        param_res = get_param_res(model)
    end
    if do_lakes
        sel_lakes = network.lake.indices_coverage
        param_lake = get_param_lake(model)
    end

    for (par, ncvar) in forcing_parameters
        if ncvar.name === nothing
            val = ncvar.value * ncvar.scale + ncvar.offset
            param_vector = param(model, par)
            param_vector .= val
            # set fixed precipitation and evaporation over the lakes and reservoirs and put
            # these into the lakes and reservoirs structs and set the precipitation and
            # evaporation to 0 in the vertical model
            if par in mover_params
                if do_reservoirs
                    for (i, sel_reservoir) in enumerate(sel_reservoirs)
                        param_vector[reverse_indices[sel_reservoir]] .= 0
                        param_res[par][i] = val
                    end
                end
                if do_lakes
                    for (i, sel_lake) in enumerate(sel_lakes)
                        param_vector[reverse_indices[sel_lake]] .= 0
                        param_lake[par][i] = val
                    end
                end
            end
        end
    end
end

"Get dynamic NetCDF input for the given time"
function update_forcing!(model)
    @unpack vertical, clock, reader, network, config = model
    @unpack dataset, forcing_parameters = reader
    nctimes = dataset["time"][:]

    do_reservoirs = get(config.model, "reservoirs", false)::Bool
    do_lakes = get(config.model, "lakes", false)::Bool

    if do_reservoirs
        sel_reservoirs = network.reservoir.indices_coverage
        param_res = get_param_res(model)
    end
    if do_lakes
        sel_lakes = network.lake.indices_coverage
        param_lake = get_param_lake(model)
    end

    # load from NetCDF into the model according to the mapping
    for (par, ncvar) in forcing_parameters
        # no need to update fixed values
        ncvar.name === nothing && continue

        time = convert(eltype(nctimes), clock.time)
        data = get_at(dataset, ncvar.name, nctimes, time)

        if ncvar.scale != 1.0 || ncvar.offset != 0.0
            data .= data .* ncvar.scale .+ ncvar.offset
        end

        # calculate the mean precipitation and evaporation over the lakes and reservoirs
        # and put these into the lakes and reservoirs structs
        # and set the precipitation and evaporation to 0 in the vertical model
        if par in mover_params
            if do_reservoirs
                for (i, sel_reservoir) in enumerate(sel_reservoirs)
                    avg = mean(data[sel_reservoir])
                    data[sel_reservoir] .= 0
                    param_res[par][i] = avg
                end
            end
            if do_lakes
                for (i, sel_lake) in enumerate(sel_lakes)
                    avg = mean(data[sel_lake])
                    data[sel_lake] .= 0
                    param_lake[par][i] = avg
                end
            end
        end

        param_vector = param(model, par)
        sel = active_indices(network, par)
        data_sel = data[sel]
        if any(ismissing, data_sel)
            print(par)
            msg = "Forcing data has missing values on active model cells for $(ncvar.name)"
            throw(ArgumentError(msg))
        end
        param_vector .= data_sel
    end

    return model
end

"""
    monthday_passed(curr, avail)

Given two monthday tuples such as (12, 31) and (12, 15), return true if the first argument
falls after or on the same day as the second argument, assuming the same year. The tuples
generally come from `Dates.monthday`.

# Examples
```julia-repl
julia> monthday_passed((12, 31), (12, 15))
true
```
"""
monthday_passed(curr, avail) = (curr[1] >= avail[1]) && (curr[2] >= avail[2])

"Get dynamic and cyclic NetCDF input"
function load_dynamic_input!(model)
    update_forcing!(model)
    if haskey(model.config.input, "cyclic")
        update_cyclic!(model)
    end
end

"Get cyclic NetCDF input for the given time"
function update_cyclic!(model)
    @unpack vertical, clock, reader, network, config = model
    @unpack cyclic_dataset, cyclic_times, cyclic_parameters = reader
    #sel = network.land.indices

    # pick up the data that is valid for the past 24 hours
    month_day = monthday(clock.time - Day(1))

    is_first_timestep = clock.iteration == 1
    if is_first_timestep || (month_day in cyclic_times)
        # time for an update of the cyclic forcing
        i = findlast(t -> monthday_passed(month_day, t), cyclic_times)
        isnothing(i) && error("Could not find applicable cyclic timestep for $month_day")

        # load from NetCDF into the model according to the mapping
        for (par, ncvar) in cyclic_parameters
            data = get_at(cyclic_dataset, ncvar.name, i)
            param_vector = param(model, par)
            sel = active_indices(network, par)
            param_vector .= data[sel]
            if ncvar.scale != 1.0 || ncvar.offset != 0.0
                param_vector .= param_vector .* ncvar.scale .+ ncvar.offset
            end
        end
    end
end

"""
    nc_handles::Dict{String, NCDataset{Nothing}}

For each NetCDF file that will be opened for writing, store an entry in this Dict from the
absolute path of the file to the NCDataset. This allows us to close the NCDataset if we try
to create them twice in the same session, and thus providing a workaround for this issue:
https://github.com/Alexander-Barth/NCDatasets.jl/issues/106

Note that using this will prevent automatic garbage collection and thus closure of the
NCDataset.
"""
const nc_handles = Dict{String,NCDataset{Nothing}}()

"Safely create a NetCDF file, even if it has already been opened for creation"
function create_tracked_netcdf(path)
    abs_path = abspath(path)
    # close existing NCDataset if it exists
    if haskey(nc_handles, abs_path)
        # fine if it was already closed
        close(nc_handles[abs_path])
    end
    # create directory if needed
    mkpath(dirname(path))
    ds = NCDataset(path, "c")
    nc_handles[abs_path] = ds
    return ds
end

"prepare an output dataset for scalar data"
function setup_scalar_netcdf(
    output_path,
    ncvars,
    calendar,
    time_units,
    float_type = Float32,
)
    ds = create_tracked_netcdf(output_path)
    defDim(ds, "time", Inf)  # unlimited
    defVar(
        ds,
        "time",
        Float64,
        ("time",),
        attrib = ["units" => time_units, "calendar" => calendar],
    )
    for nc in ncvars
        # Delft-FEWS requires the attribute :cf_role = "timeseries_id" when a NetCDF file
        # contains more than one location list
        defVar(
            ds,
            nc.location_dim,
            nc.locations,
            (nc.location_dim,),
            attrib = ["cf_role" => "timeseries_id"],
        )
        defVar(
            ds,
            nc.var,
            float_type,
            (nc.location_dim, "time"),
            attrib = ["_FillValue" => float_type(NaN)],
        )
    end
    return ds
end

"prepare an output dataset for grid data"
function setup_grid_netcdf(
    output_path,
    ncx,
    ncy,
    parameters,
    calendar,
    time_units,
    maxlayers,
    sizeinmetres;
    float_type = Float32,
)

    ds = create_tracked_netcdf(output_path)
    defDim(ds, "time", Inf)  # unlimited
    if sizeinmetres
        defVar(
            ds,
            "x",
            ncx,
            ("x",),
            attrib = [
                "_FillValue" => NaN,
                "long_name" => "x coordinate of projection",
                "standard_name" => "projection_x_coordinate",
                "axis" => "X",
                "units" => "m",
            ],
        )
        defVar(
            ds,
            "y",
            ncy,
            ("y",),
            attrib = [
                "_FillValue" => NaN,
                "long_name" => "y coordinate of projection",
                "standard_name" => "projection_y_coordinate",
                "axis" => "Y",
                "units" => "m",
            ],
        )

    else
        defVar(
            ds,
            "lon",
            ncx,
            ("lon",),
            attrib = [
                "_FillValue" => NaN,
                "long_name" => "longitude",
                "standard_name" => "longitude",
                "axis" => "X",
                "units" => "degrees_east",
            ],
        )
        defVar(
            ds,
            "lat",
            ncy,
            ("lat",),
            attrib = [
                "_FillValue" => NaN,
                "long_name" => "latitude",
                "standard_name" => "latitude",
                "axis" => "Y",
                "units" => "degrees_north",
            ],
        )
    end
    if isnothing(maxlayers) == false
        defVar(ds, "layer", collect(1:maxlayers), ("layer",))
    end
    defVar(
        ds,
        "time",
        Float64,
        ("time",),
        attrib = ["units" => time_units, "calendar" => calendar],
    )
    if sizeinmetres
        for (key, val) in pairs(parameters)
            if eltype(val.vector) <: AbstractFloat
                # all floats are saved as Float32
                defVar(
                    ds,
                    key,
                    float_type,
                    ("x", "y", "time"),
                    attrib = ["_FillValue" => float_type(NaN)],
                )
            elseif eltype(val.vector) <: SVector
                # SVectors are used to store layers
                defVar(
                    ds,
                    key,
                    float_type,
                    ("x", "y", "layer", "time"),
                    attrib = ["_FillValue" => float_type(NaN)],
                )
            else
                error("Unsupported output type: ", typeof(val.vector))
            end
        end
    else
        for (key, val) in pairs(parameters)
            if eltype(val.vector) <: AbstractFloat
                # all floats are saved as Float32
                defVar(
                    ds,
                    key,
                    float_type,
                    ("lon", "lat", "time"),
                    attrib = ["_FillValue" => float_type(NaN)],
                )
            elseif eltype(val.vector) <: SVector
                # SVectors are used to store layers
                defVar(
                    ds,
                    key,
                    float_type,
                    ("lon", "lat", "layer", "time"),
                    attrib = ["_FillValue" => float_type(NaN)],
                )
            else
                error("Unsupported output type: ", typeof(val.vector))
            end
        end
    end
    return ds
end

"Add a new time to the unlimited time dimension, and return the index"
function add_time(ds, time)
    i = length(ds["time"]) + 1
    ds["time"][i] = time
    return i
end

struct NCReader
    dataset::CFDataset
    cyclic_dataset::Union{NCDataset,Nothing}
    cyclic_times::Vector{Tuple{Int,Int}}
    forcing_parameters::Dict{Tuple{Symbol,Vararg{Symbol}},NamedTuple}
    cyclic_parameters::Dict{Tuple{Symbol,Vararg{Symbol}},NamedTuple}
end

struct Writer
    dataset::Union{NCDataset,Nothing}           # dataset (NetCDF) for grid data
    parameters::Dict{String,Any}                # mapping of NetCDF variable names to model parameters (arrays)
    nc_path::Union{String,Nothing}              # path NetCDF file (grid data)
    csv_path::Union{String,Nothing}             # path of CSV file
    csv_cols::Vector                            # model parameter (arrays) and associated reducer function for CSV output
    csv_io::IO                                  # file handle to CSV file
    state_dataset::Union{NCDataset,Nothing}     # dataset with model states (NetCDF)
    state_parameters::Dict{String,Any}          # mapping of NetCDF variable names to model states (arrays)
    state_nc_path::Union{String,Nothing}        # path NetCDF file with states
    dataset_scalar::Union{NCDataset,Nothing}    # dataset(NetCDF) for scalar data
    nc_scalar::Vector                           # model parameter (arrays) and associated reducer function for NetCDF scalar output
    ncvars_dims::Vector                         # NetCDF variable, location dimension and location name for scalar data
    nc_scalar_path::Union{String,Nothing}       # path NetCDF file (scalar data)
end

function prepare_reader(config)
    tomldir = dirname(config)
    path_forcing = config.input.path_forcing
    cyclic_path = joinpath(tomldir, config.input.path_static)

    # absolute paths are not supported, see Glob.jl#2
    if isabspath(path_forcing)
        parts = splitpath(path_forcing)
        # use the root/drive as the dir, to support * in directory names as well
        glob_dir = parts[1]
        glob_path = joinpath(parts[2:end])
    else
        glob_dir = tomldir
        glob_path = path_forcing
    end

    dynamic_paths = glob(glob_path, glob_dir)  # expand "data/forcing-year-*.nc"
    dataset = NCDataset(dynamic_paths, aggdim = "time", deferopen = false)

    # check for cyclic parameters
    do_cyclic = haskey(config.input, "cyclic")

    # TODO:include leaf_area_index climatology in update() vertical SBM model
    # we currently assume the same dimension ordering as the forcing
    if do_cyclic == true
        cyclic_dataset = NCDataset(cyclic_path)
        cyclic_nc_times = collect(cyclic_dataset["time"])
        cyclic_times = timecycles(cyclic_nc_times)
    else
        cyclic_dataset = nothing
        cyclic_times = Tuple{Int,Int}[]
    end

    # create map from internal location to NetCDF variable name for forcing parameters
    forcing_parameters = Dict{Tuple{Symbol,Vararg{Symbol}},NamedTuple}()
    for par in config.input.forcing
        fields = symbols(par)
        ncname, mod = ncvar_name_modifier(param(config.input, fields))
        forcing_parameters[fields] =
            (name = ncname, scale = mod.scale, offset = mod.offset, value = mod.value)
    end

    # create map from internal location to NetCDF variable name for cyclic parameters
    if do_cyclic == true
        cyclic_parameters = Dict{Tuple{Symbol,Vararg{Symbol}},NamedTuple}()
        for par in config.input.cyclic
            fields = symbols(par)
            ncname, mod = ncvar_name_modifier(param(config.input, fields))
            cyclic_parameters[fields] =
                (name = ncname, scale = mod.scale, offset = mod.offset)
        end
    else
        cyclic_parameters = Dict{Tuple{Symbol,Vararg{Symbol}},NamedTuple}()
    end

    # check if there is overlap
    if do_cyclic == true
        forcing_keys = keys(forcing_parameters)
        cyclic_keys = keys(cyclic_parameters)
        for forcing_key in forcing_keys
            if forcing_key in cyclic_keys
                error("parameter specified in both forcing and cyclic")
            end
        end
    end

    return NCReader(
        dataset,
        cyclic_dataset,
        cyclic_times,
        forcing_parameters,
        cyclic_parameters,
    )
end

"Get a Vector of all unique location ids from a 2D map"
function locations_map(ds, mapname, config)
    map_2d = ncread(
        ds,
        param(config.input, mapname);
        type = Union{Int,Missing},
        allow_missing = true,
    )
    ids = unique(skipmissing(map_2d))
    return ids
end

"Get a Vector{Tuple} with NetCDF variable, dimension and location names for scalar data"
function nc_variables_dims(nc_variables, dataset, config)
    ncvars_dims = []
    for nc_var in nc_variables
        var = nc_var["name"]::String
        if haskey(nc_var, "map")
            mapname = nc_var["map"]
            ids = string.(locations_map(dataset, mapname, config))
            location_dim = string(var, '_', nc_var["map"])
            push!(ncvars_dims, (var = var, location_dim = location_dim, locations = ids))
        else
            push!(
                ncvars_dims,
                (
                    var = var,
                    location_dim = nc_var["location"],
                    locations = [nc_var["location"]],
                ),
            )
        end
    end
    return ncvars_dims
end

"Get a Vector{String} of all columns names for the CSV header, exept the first, time"
function csv_header(cols, dataset, config)
    header = [col["header"] for col in cols]
    header = String[]
    for col in cols
        h = col["header"]::String
        if haskey(col, "map")
            mapname = col["map"]
            ids = locations_map(dataset, mapname, config)
            hvec = [string(h, '_', id) for id in ids]
            append!(header, hvec)
        else
            push!(header, h)
        end
    end
    return header
end

"""
Flatten a nested dictionary, keeping track of the full address of the keys.
Useful for converting TOML of the format:

    [a.b]
    field_of_b = "name_in_output"

    [a.d.c]
    field_of_c = "other_name_in_output"

to a non-nested format:

Dict(
    symbols"a.b.field_of_b" => "name_in_output,
    symbols"a.d.c.field_of_c" => "other_name_in_output,
)
"""
function flat!(d, path, el::Dict)
    for (k, v) in pairs(el)
        flat!(d, string(path, '.', k), v)
    end
    return d
end

function flat!(d, path, el)
    k = symbols(path)
    d[k] = el
    return d
end

"""
    ncnames(dict)

Create a flat mapping from internal parameter locations to NetCDF variable names.

Ignores top level values in the Dict. This function is used to convert a TOML such as:

```toml
[output]
path = "path/to/file.nc"

[output.vertical]
canopystorage = "my_canopystorage"

[output.lateral.river]
q = "my_q"
```

To a dictionary of the flattened parameter locations and NetCDF names. The top level
values are ignored since the output path is not a NetCDF name.

```julia
Dict(
    (:vertical, :canopystorage) => "my_canopystorage,
    (:lateral, :river, :q) => "my_q,
)
```
"""
function ncnames(dict)
    ncnames_dict = Dict{Tuple{Symbol,Vararg{Symbol}},String}()
    for (k, v) in dict
        if v isa Dict  # ignore top level values (e.g. output.path)
            flat!(ncnames_dict, k, v)
        end
    end
    return ncnames_dict
end

"""
    out_map(ncnames_dict, modelmap)

Create a Dict that maps parameter NetCDF names to arrays in the Model.
"""
function out_map(ncnames_dict, modelmap)
    output_map = Dict{String,Any}()
    for (par, ncname) in ncnames_dict
        A = param(modelmap, par)
        output_map[ncname] = (par = par, vector = A)
    end
    return output_map
end

function get_reducer_func(col, rev_inds, args...)
    parameter = col["parameter"]
    if occursin("reservoir", parameter)
        reducer_func = reducer(col, rev_inds.reservoir, args...)
    elseif occursin("lake", parameter)
        reducer_func = reducer(col, rev_inds.lake, args...)
    elseif occursin("river", parameter)
        reducer_func = reducer(col, rev_inds.river, args...)
    elseif occursin("drain", parameter)
        reducer_func = reducer(col, rev_inds.drain, args...)
    else
        reducer_func = reducer(col, rev_inds.land, args...)
    end
end

function prepare_writer(
    config,
    reader,
    modelmap,
    rev_inds,
    x_nc,
    y_nc,
    nc_static;
    maxlayers = nothing,
)
    tomldir = dirname(config)
    sizeinmetres = get(config.model, "sizeinmetres", false)::Bool

    calendar = get(config, "calendar", "standard")::String
    time_units = get(config, "time_units", CFTime.DEFAULT_TIME_UNITS)

    # create an output NetCDF that will hold all timesteps of selected parameters for grid
    # data but only if config.output.path has been set
    if haskey(config, "output") && haskey(config.output, "path")
        nc_path = joinpath(tomldir, config.output.path)
        # create a flat mapping from internal parameter locations to NetCDF variable names
        output_ncnames = ncnames(config.output)
        # fill the output_map by mapping parameter NetCDF names to arrays
        output_map = out_map(output_ncnames, modelmap)
        ds = setup_grid_netcdf(
            nc_path,
            x_nc,
            y_nc,
            output_map,
            calendar,
            time_units,
            maxlayers,
            sizeinmetres,
        )
    else
        nc_path = nothing
        output_map = Dict{String,Any}()
        ds = nothing
    end

    # create a separate state output NetCDF that will hold the last timestep of all states
    # but only if config.state.path_output has been set
    if haskey(config, "state") && haskey(config.state, "path_output")
        state_ncnames = ncnames(config.state)
        state_map = out_map(state_ncnames, modelmap)
        nc_state_path = joinpath(tomldir, config.state.path_output)
        ds_outstate = setup_grid_netcdf(
            nc_state_path,
            x_nc,
            y_nc,
            state_map,
            calendar,
            time_units,
            maxlayers,
            sizeinmetres;
            float_type = Float,
        )
    else
        ds_outstate = nothing
        state_map = Dict{String,Any}()
        nc_state_path = nothing
    end

    # create an output NetCDF that will hold all timesteps of selected parameters for scalar
    # data, but only if config.netcdf.variable has been set.
    if haskey(config, "netcdf") && haskey(config.netcdf, "variable")
        nc_scalar_path = joinpath(tomldir, config.netcdf.path)
        # get NetCDF info for scalar data (variable name, locationset (dim) and
        # location ids)
        ncvars_dims = nc_variables_dims(config.netcdf.variable, nc_static, config)
        ds_scalar = setup_scalar_netcdf(nc_scalar_path, ncvars_dims, calendar, time_units)
        # create a vector of (parameter, reducer) named tuples which will be used to
        # retrieve and reduce data during a model run
        nc_scalar = []
        for var in config.netcdf.variable
            parameter = var["parameter"]
            reducer_func = get_reducer_func(var, rev_inds, x_nc, y_nc, config, nc_static)
            push!(nc_scalar, (parameter = parameter, reducer = reducer_func))
        end
    else
        ds_scalar = nothing
        nc_scalar = []
        ncvars_dims = []
        nc_scalar_path = nothing
    end

    if haskey(config, "csv") && haskey(config.csv, "column")
        # open CSV file and write header
        csv_path = joinpath(tomldir, config.csv.path)
        # create directory if needed
        mkpath(dirname(csv_path))
        csv_io = open(csv_path, "w")
        print(csv_io, "time,")
        header = csv_header(config.csv.column, nc_static, config)
        println(csv_io, join(header, ','))
        flush(csv_io)

        # create a vector of (parameter, reducer) named tuples which will be used to
        # retrieve and reduce the CSV data during a model run
        csv_cols = []
        for col in config.csv.column
            parameter = col["parameter"]
            reducer_func = get_reducer_func(col, rev_inds, x_nc, y_nc, config, nc_static)
            push!(csv_cols, (parameter = parameter, reducer = reducer_func))
        end
    else
        # no CSV file is checked by isnothing(csv_path)
        csv_path = nothing
        csv_cols = []
        csv_io = IOBuffer()
    end

    return Writer(
        ds,
        output_map,
        nc_path,
        csv_path,
        csv_cols,
        csv_io,
        ds_outstate,
        state_map,
        nc_state_path,
        ds_scalar,
        nc_scalar,
        ncvars_dims,
        nc_scalar_path,
    )
end

"Write a new timestep with scalar data to a NetCDF file"
function write_netcdf_timestep(model, dataset)
    @unpack writer, clock = model

    time_index = add_time(dataset, clock.time)
    for (nt, nc) in zip(writer.nc_scalar, writer.ncvars_dims)
        A = param(model, nt.parameter)
        # could be a value, or a vector in case of map
        v = nt.reducer(A)
        dataset[nc.var][:, time_index] .= v
    end
    return model
end

"Write a new timestep with grid data to a NetCDF file"
function write_netcdf_timestep(model, dataset, parameters)
    @unpack vertical, clock, reader, network = model

    time_index = add_time(dataset, clock.time)

    buffer = zeros(Union{Float,Missing}, size(model.network.land.reverse_indices))
    for (key, val) in parameters
        @unpack par, vector = val
        sel = active_indices(network, par)
        # write the active cells vector to the 2d buffer matrix
        elemtype = eltype(vector)
        if elemtype <: AbstractFloat
            # ensure no other information is written
            fill!(buffer, missing)
            buffer[sel] .= vector
            dataset[key][:, :, time_index] = buffer
        elseif elemtype <: SVector
            nlayer = length(first(vector))
            for i = 1:nlayer
                # ensure no other information is written
                fill!(buffer, missing)
                buffer[sel] .= getindex.(vector, i)
                dataset[key][:, :, i, time_index] = buffer
            end
        else
            error("Unsupported output type: ", elemtype)
        end
    end

    return model
end

# don't do anything for no dataset, used if no output NetCDF is needed
write_netcdf_timestep(model, dataset::Nothing, parameters) = model
write_netcdf_timestep(model, dataset::Nothing) = model

"Write model output"
function write_output(model)
    @unpack vertical, clock, reader, network, writer = model
    @unpack dataset, dataset_scalar, parameters = writer

    write_csv_row(model)
    write_netcdf_timestep(model, dataset, parameters)
    write_netcdf_timestep(model, dataset_scalar)

    return model
end

"""
    timecycles(times)

Given a vector of times, return a tuple of (month, day) for each time entry, to use as a
cyclic time series that repeats every year. By using `monthday` rather than `dayofyear`,
leap year offsets are avoided.

It can generate such a series from eiher TimeTypes given that the year is constant, or
it will interpret integers as either months or days of year if possible.
"""
function timecycles(times)
    if eltype(times) <: TimeType
        # all timestamps are from the same year
        year1 = year(first(times))
        if !all(==(year1), year.(times))
            error("unsupported cyclic timeseries")
        end
        # returns a (month, day) tuple for each date
        return monthday.(times)
    elseif eltype(times) <: Integer
        if length(times) == 12
            months = Date(2000, 1, 1):Month(1):Date(2000, 12, 31)
            return monthday.(months)
        elseif length(times) == 365
            days = Date(2001, 1, 1):Day(1):Date(2001, 12, 31)
            return monthday.(days)
        elseif length(times) == 366
            days = Date(2000, 1, 1):Day(1):Date(2000, 12, 31)
            return monthday.(days)
        else
            error("unsupported cyclic timeseries")
        end
    else
        error("unsupported cyclic timeseries")
    end
end

"Close input and output datasets that are opened on model initialization"
function close_files(model; delete_output::Bool = false)
    @unpack reader, writer, config = model

    close(reader.dataset)
    if haskey(config.input, "cyclic")
        close(reader.cyclic_dataset)
    end
    writer.dataset === nothing || close(writer.dataset)
    writer.dataset_scalar === nothing || close(writer.dataset_scalar)
    close(writer.csv_io)  # can be an IOBuffer
    writer.state_dataset === nothing || close(writer.state_dataset)

    if delete_output
        if writer.nc_path !== nothing
            isfile(writer.nc_path) && rm(writer.nc_path)
        end
        if writer.csv_path !== nothing
            isfile(writer.csv_path) && rm(writer.csv_path)
        end
        if writer.state_nc_path !== nothing
            isfile(writer.state_nc_path) && rm(writer.state_nc_path)
        end
        if writer.nc_scalar_path !== nothing
            isfile(writer.nc_scalar_path) && rm(writer.nc_scalar_path)
        end
    end
    return nothing
end

"Mapping from reducer strings in the TOML to functions"
function reducerfunction(reducer::AbstractString)
    functionmap = Dict{String,Function}(
        "maximum" => maximum,
        "minimum" => minimum,
        "mean" => mean,
        "median" => median,
        "first" => first,
        "last" => last,
        "only" => only,
    )
    f = get(functionmap, reducer, nothing)
    isnothing(f) && error("unknown reducer")
    return f
end

"Get a reducer function based on output settings for scalar data defined in a dictionary"
function reducer(col, rev_inds, x_nc, y_nc, config, dataset)
    if haskey(col, "map")
        # assumes the parameter in "map" has a 2D input map, with
        # integers indicating the points or zones that are to be aggregated
        mapname = col["map"]
        # if no reducer is given, pick "only", this is the only safe reducer,
        # and makes sense in the case of a gauge map
        reducer_name = get(col, "reducer", "only")
        f = reducerfunction(reducer_name)
        map_2d = ncread(
            dataset,
            param(config.input, mapname);
            type = Union{Int,Missing},
            allow_missing = true,
        )
        ids = unique(skipmissing(map_2d))
        # from id to list of internal indices
        inds = Dict{Int,Vector{Int}}(id => Vector{Int}() for id in ids)
        for i in eachindex(map_2d)
            v = map_2d[i]
            ismissing(v) && continue
            v::Int
            vector = inds[v]
            ind = rev_inds[i]
            if iszero(ind)
                param = col["parameter"]
                error("""inactive cell found in requested CSV output
                    map \"$mapname\" value $v for parameter $param""")
            end
            push!(vector, ind)
        end
        return A -> (f(A[inds[id]]) for id in ids)
    elseif haskey(col, "reducer")
        # reduce over all active cells
        # needs to be behind the map if statement, because it also can use a reducer
        return reducerfunction(col["reducer"])
    elseif haskey(col, "index")
        index = col["index"]
        if index isa Int
            # linear index into the internal vector of active cells
            # this one mostly makes sense for debugging, or for vectors of only a few elements
            return x -> getindex(x, index)
        elseif index isa Dict
            # index into the 2D input/output arrays
            # the first always corresponds to the x dimension, then the y dimension
            # this is 1-based
            i = index["x"]::Int
            j = index["y"]::Int
            ind = rev_inds[i, j]
            iszero(ind) && error("inactive loc specified for output")
            return A -> getindex(A, ind)
        else
            error("unknown index used")
        end
    elseif haskey(col, "coordinate")
        x = col["coordinate"]["x"]::Float64
        y = col["coordinate"]["y"]::Float64
        # find the closest cell center index
        _, iy = findmin(abs.(y_nc .- y))
        _, ix = findmin(abs.(x_nc .- x))
        I = CartesianIndex(ix, iy)
        i = rev_inds[I]
        iszero(i) && error("inactive coordinate specified for output")
        return A -> getindex(A, i)
    else
        error("unknown reducer")
    end
end

function write_csv_row(model)
    @unpack writer, clock = model
    isnothing(writer.csv_path) && return nothing
    io = writer.csv_io
    print(io, string(clock.time))
    for nt in writer.csv_cols
        A = param(model, nt.parameter)
        # could be a value, or a vector in case of map
        v = nt.reducer(A)
        # numbers are also iterable
        for el in v
            print(io, ',', el)
        end
    end
    println(io)
end

"From a time and a calendar, create the right CFTime DateTimeX type"
function cftime(time, calendar)
    timetype = CFTime.timetype(calendar)
    # invalid Gregorian dates like 2020-02-30 cannot be made into a DateTime
    # even though they may exist in the target calendar
    # though this constructor does offer support for string formatted dates
    dt = DateTime(time)
    cal_time = reinterpret(timetype, dt)
    return cal_time
end

function reset_clock!(clock::Clock, config)
    new_clock = Clock(config)
    # we want this method to be mutating
    clock.time = new_clock.time
    clock.iteration = new_clock.iteration
    clock.Δt = new_clock.Δt
    return clock
end

function advance!(clock)
    clock.iteration += 1
    clock.time += clock.Δt
    return clock
end

function rewind!(clock)
    clock.iteration -= 1
    clock.time -= clock.Δt
    return clock
end

"Read a rating curve from CSV into a NamedTuple of vectors"
function read_sh_csv(path)
    data, header = readdlm(path, ',', Float, header = true)
    names = vec(uppercase.(header))
    idx_h = findfirst(==("H"), names)
    idx_s = findfirst(==("S"), names)

    if isnothing(idx_h) || isnothing(idx_s)
        error("$path needs to provide H and S columns, got $names")
    end

    return (H = data[:, idx_h], S = data[:, idx_s])
end

"Read a specific storage curve from CSV into a NamedTuple of vectors"
function read_hq_csv(path)
    data = readdlm(path, ',', Float, skipstart = 1)
    # Q is a matrix with 365 columns, one for each day in the year
    return (H = data[:, 1], Q = data[:, 2:end])
end

# these represent the type of the rating curve and specific storage data
const SH = NamedTuple{(:H, :S),Tuple{Vector{Float},Vector{Float}}}
const HQ = NamedTuple{(:H, :Q),Tuple{Vector{Float},Matrix{Float}}}

is_increasing(v) = last(v) > first(v)

"""
    nc_dim_name(dims::Vector{Symbol}, name::Symbol)
    nc_dim_name(ds::CFDataset, name::Symbol)

Given a NetCDF dataset or list of dimensions, and an internal dimension name, return the
corresponding NetCDF dimension name. Certain common alternatives are supported, e.g. :lon or
:longitude instead of :x.
"""
function nc_dim_name(dims::Vector{Symbol}, name::Symbol)
    # direct naming
    name in dims && return name
    # list of common alternative names
    if name == :x
        for candidate in (:lon, :longitude)
            candidate in dims && return candidate
        end
        error("No x dimension candidate found")
    elseif name == :y
        for candidate in (:lat, :latitude)
            candidate in dims && return candidate
        end
        error("No y dimension candidate found")
    else
        error("Unknown dimension $name")
    end
end

nc_dim_name(ds::CFDataset, name::Symbol) = nc_dim_name(Symbol.(keys(ds.dim)), name)

"""
    nc_dim(ds::CFDataset, name::Symbol)

Return the dimension coordinate, based on the internal name (:x, :y, :layer, :time),
which will map to the correct NetCDF name using `nc_dim_name`.
"""
nc_dim(ds::CFDataset, name) = ds[nc_dim_name(ds, name)]

"""
    internal_dim_name(name::Symbol)

Given a NetCDF dimension name string, return the corresponding internal dimension name.
"""
function internal_dim_name(name::Symbol)
    if name in (:x, :lon, :longitude)
        return :x
    elseif name in (:y, :lat, :latitude)
        return :y
    elseif name in (:time, :layer)
        return name
    else
        error("Unknown dimension $name")
    end
end

"""
    read_dims(A::NCDatasets.CFVariable, dim_sel)

Return the data of a NetCDF data variable as an Array. Only dimensions in `dim_sel`, a
NamedTuple like (x=:, y=:, time=1). Other dimensions that may be present need to be size 1,
otherwise an error is thrown.

`dim_sel` keys should be the internal dimension names; :x, :y, :time, :layer.
"""
function read_dims(A::NCDatasets.CFVariable, dim_sel::NamedTuple)
    dimsizes = dimsize(A)
    indexer = []
    data_dim_order = Symbol[]
    # need to iterate in this order, to get the indices in order
    for (dim_name, dim_size) in pairs(dimsizes)
        if dim_size == 1
            push!(indexer, 1)
        else
            dim = internal_dim_name(dim_name)
            # each dim should either be in dim_sel or be size 1
            if dim in keys(dim_sel)
                idx = dim_sel[dim]
                push!(indexer, idx)
                # Would be nice to generalize this to anything that drops the dimension
                # when used to index. Not sure how we could support `begin` or `end` here.
                if !(idx isa Int)
                    push!(data_dim_order, dim)
                end
            else
                throw(ArgumentError("""NetCDF dimension $dim_name has length $dim_size.
                    Only extra dimensions of length 1 are supported."""))
            end
        end
    end
    data = A[indexer...]
    order = Tuple(data_dim_order)
    @assert ndims(data) == length(order) "dimensions not properly recorded"
    return data, order
end

"""
    dim_directions(ds, dim_names)

Given a NCDataset and a list of internal dimension names, return a NamedTuple, which maps
the dimension name to `true` if it is increasing, or `false` otherwise.

For the layer dimension, we allow coordinate arrays to be missing, in which case we
consider it increasing, going from the top layer (1) to deeper layers. This is to keep
accepting data that we have accepted before.
"""
function dim_directions(ds::CFDataset, dim_names)
    pairs = Pair{Symbol,Bool}[]
    for d in dim_names
        if d == :layer && !(haskey(ds, "layer"))
            inc = true
        else
            inc = is_increasing(nc_dim(ds, d))
        end
        push!(pairs, d => inc)
    end
    return NamedTuple(pairs)
end

"""
    permute_data(data, dim_names)

Given an Array of data, and a list of its dimension names, return a permuted version
of the data such that the dimension order will follow (:x, :y, :layer). No permutation
is done if this is not needed.
"""
function permute_data(data, dim_names)
    @assert ndims(data) == length(dim_names)
    desired_order = (:x, :y, :layer)
    @assert dim_names ⊆ desired_order
    if length(dim_names) == 2
        if first(dim_names) == :x
            return data, dim_names
        else
            return permutedims(data), reverse(dim_names)
        end
    elseif length(dim_names) == 3
        permutation = [findfirst(==(d), dim_names) for d in desired_order]
        if permutation == (1, 2, 3)
            return data, dim_names
        else
            return permutedims(data, permutation), dim_names[permutation]
        end
    else
        error("Unsupported number of dimensions")
    end
end

"""
    reverse_data!(data, dims_increasing)

Reverse the data as needed, such that it is increasing along all dimensions.
"""
function reverse_data!(data, dims_increasing)
    # for the reverse call it is important that the dims_increasing tuple is ordered in the
    # desired internal ordering, just like the data is after permutation
    if length(dims_increasing) == 2
        dims_increasing_ordered = (x = dims_increasing.x, y = dims_increasing.y)
    elseif length(dims_increasing) == 3
        dims_increasing_ordered =
            (x = dims_increasing.x, y = dims_increasing.y, layer = dims_increasing.layer)
    else
        error("Unsupported number of dimensions")
    end
    # the non increasing dimensions should be reversed, to make them all increasing
    dims = Tuple(findall(.!values(dims_increasing_ordered)))
    return reverse!(data; dims)
end

"""
    read_standardized(ds::CFDataset, varname::AbstractString, dim_names)

Read the dimensions listed in dim_names from a variable with name `varname` from a NetCDF
dataset `ds`. `dim_sel` should be a NamedTuple like (x=:, y=:, time=1), which will return
a 2 dimensional array with x and y axes, representing the first index in the time dimension.
"""
function read_standardized(ds::CFDataset, varname::AbstractString, dim_sel::NamedTuple)
    data, data_dim_order = read_dims(ds[varname], dim_sel)
    data, new_dim_order = permute_data(data, data_dim_order)
    dims_increasing = dim_directions(ds, new_dim_order)
    data = reverse_data!(data, dims_increasing)
    return data
end

"""
    read_x_axis(ds::CFDataset)

Return the x coordinate Vector{Float64}, whether it is called x, lon or longitude.
Also sorts the vector to be increasing, to match `read_standardized`.
"""
function read_x_axis(ds::CFDataset)
    candidates = ("x", "lon", "longitude")
    for candidate in candidates
        if haskey(ds, candidate)
            return sort!(Float64.(ds[candidate][:]))
        end
    end
    error("no x axis found in $(path(ds))")
end

"""
    read_y_axis(ds::CFDataset)

Return the y coordinate Vector{Float64}, whether it is called y, lat or latitude.
Also sorts the vector to be increasing, to match `read_standardized`.
"""
function read_y_axis(ds::CFDataset)
    candidates = ("y", "lat", "latitude")
    for candidate in candidates
        if haskey(ds, candidate)
            return sort!(Float64.(ds[candidate][:]))
        end
    end
    error("no y axis found in $(path(ds))")
end


text:
using StatsBase
using Base.Test

##### weighted var & std

x = rand(10)
wv = weights(rand(10))
m = mean(x, wv)

@test_approx_eq var(x, wv) sum(abs2(x .- m), wv) ./ sum(wv)
@test_approx_eq var(x, wv; mean=0) sum(abs2(x), wv) ./ sum(wv)
@test_approx_eq var(x, wv; mean=1.0) sum(abs2(x .- 1.0), wv) ./ sum(wv)

@test_approx_eq std(x, wv) sqrt(var(x, wv))
@test_approx_eq std(x, wv; mean=0) sqrt(var(x, wv; mean=0))
@test_approx_eq std(x, wv; mean=1.0) sqrt(var(x, wv; mean=1.0))

(m, v) = mean_and_var(x)
@test m == mean(x)
@test v == var(x)

(m, s) = mean_and_std(x)
@test m == mean(x)
@test s == std(x)

(m, v) = mean_and_var(x, wv)
@test m == mean(x, wv)
@test v == var(x, wv)

(m, s) = mean_and_std(x, wv)
@test m == mean(x, wv)
@test s == std(x, wv)

x = rand(5, 6)
w1 = rand(5)
w2 = rand(6)
wv1 = weights(w1)
wv2 = weights(w2)
m1 = mean(x, wv1, 1)
m2 = mean(x, wv2, 2)

@test_approx_eq var(x, wv1, 1; mean=0) sum(abs2(x) .* w1, 1) ./ sum(wv1)
@test_approx_eq var(x, wv2, 2; mean=0) sum(abs2(x) .* w2', 2) ./ sum(wv2)

@test_approx_eq var(x, wv1, 1; mean=m1) sum(abs2(x .- m1) .* w1, 1) ./ sum(wv1)
@test_approx_eq var(x, wv2, 2; mean=m2) sum(abs2(x .- m2) .* w2', 2) ./ sum(wv2)

@test_approx_eq var(x, wv1, 1) sum(abs2(x .- m1) .* w1, 1) ./ sum(wv1)
@test_approx_eq var(x, wv2, 2) sum(abs2(x .- m2) .* w2', 2) ./ sum(wv2)

@test_approx_eq std(x, wv1, 1) sqrt(var(x, wv1, 1))
@test_approx_eq std(x, wv2, 2) sqrt(var(x, wv2, 2))
@test_approx_eq std(x, wv1, 1; mean=0) sqrt(var(x, wv1, 1; mean=0))
@test_approx_eq std(x, wv2, 2; mean=0) sqrt(var(x, wv2, 2; mean=0))
@test_approx_eq std(x, wv1, 1; mean=m1) sqrt(var(x, wv1, 1; mean=m1))
@test_approx_eq std(x, wv2, 2; mean=m2) sqrt(var(x, wv2, 2; mean=m2))

for d in 1:2
    (m, v) = mean_and_var(x, d)
    @test m == mean(x, d)
    @test v == var(x, d)

    (m, s) = mean_and_std(x, d)
    @test m == mean(x, d)
    @test s == std(x, d)
end

(m, v) = mean_and_var(x, wv1, 1)
@test m == mean(x, wv1, 1)
@test v == var(x, wv1, 1)

(m, v) = mean_and_var(x, wv2, 2)
@test m == mean(x, wv2, 2)
@test v == var(x, wv2, 2)

(m, s) = mean_and_std(x, wv1, 1)
@test m == mean(x, wv1, 1)
@test s == std(x, wv1, 1)

(m, s) = mean_and_std(x, wv2, 2)
@test m == mean(x, wv2, 2)
@test s == std(x, wv2, 2)

##### skewness & kurtosis

wv = weights(ones(5) * 2.0)

@test_approx_eq skewness(1:5) 0.0
@test_approx_eq skewness([1, 2, 3, 4, 5]) 0.0
@test_approx_eq skewness([1, 2, 2, 2, 5]) 1.1731251294063556
@test_approx_eq skewness([1, 4, 4, 4, 5]) -1.1731251294063556

@test_approx_eq skewness([1, 2, 2, 2, 5], wv) 1.1731251294063556

@test_approx_eq kurtosis(1:5) -1.3
@test_approx_eq kurtosis([1, 2, 3, 4, 5]) -1.3
@test_approx_eq kurtosis([1, 2, 3, 3, 2]) -1.1530612244897953

@test_approx_eq kurtosis([1, 2, 3, 4, 5], wv) -1.3

##### general moments

x = collect(2.0:8.0)
@test_approx_eq moment(x, 2) sum((x .- 5).^2) / length(x)
@test_approx_eq moment(x, 3) sum((x .- 5).^3) / length(x)
@test_approx_eq moment(x, 4) sum((x .- 5).^4) / length(x)
@test_approx_eq moment(x, 5) sum((x .- 5).^5) / length(x)

@test_approx_eq moment(x, 2, 4.0) sum((x .- 4).^2) / length(x)
@test_approx_eq moment(x, 3, 4.0) sum((x .- 4).^3) / length(x)
@test_approx_eq moment(x, 4, 4.0) sum((x .- 4).^4) / length(x)
@test_approx_eq moment(x, 5, 4.0) sum((x .- 4).^5) / length(x)

w = weights([1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0])
x2 = collect(2.0:6.0)
@test_approx_eq moment(x, 2, w) sum((x2 .- 4).^2) / 5
@test_approx_eq moment(x, 3, w) sum((x2 .- 4).^3) / 5
@test_approx_eq moment(x, 4, w) sum((x2 .- 4).^4) / 5
@test_approx_eq moment(x, 5, w) sum((x2 .- 4).^5) / 5